/*******************************************
*
*              FILE   POINTS.C
*
*  This file contains all the functions unique to the calculation
*  of the image of the binary at a specific orbital phase made from
*  a plot of the visible grid points.
*
*  NOTE:  Orbital phase 0.0 is defined to be inferior conjunction
*               of star2 (= star2 closest to the earth).
*         Also, internally, the program forces -0.5 < orbphase <= 0.5.
*
**********************************************/

#include "header.h"

void pointImage( double orbphase )
/*************************
*
*  This function creates the file "points.dat" with the (x,y) positions
*  of all the visible grid points in the system.  It should be identical
*  to the function totflux.c except for the final three double for()
*  loops, which write the positions of the points to a file instead of
*  adding their flux to the total flux.
*
**************************/
{
   FILE *out;
   char plotfile[40];
   double          x1 [SIZEDIMEN1] [SIZEDIMEN2],
                   y1 [SIZEDIMEN1] [SIZEDIMEN2],
                  mu1 [SIZEDIMEN1] [SIZEDIMEN2],
                   x2 [SIZEDIMEN1] [SIZEDIMEN2],
                   y2 [SIZEDIMEN1] [SIZEDIMEN2],
                  mu2 [SIZEDIMEN1] [SIZEDIMEN2],
                xdisk [SIZEDIMEN1] [SIZEDIMEN2],
                ydisk [SIZEDIMEN1] [SIZEDIMEN2],
               mudisk [SIZEDIMEN1] [SIZEDIMEN2];
   double         edge1x[LINEPOINTS],         edge1y[LINEPOINTS],
                  edge2x[LINEPOINTS],         edge2y[LINEPOINTS],
          outeredgediskx[LINEPOINTS], outeredgedisky[LINEPOINTS],
            topedgediskx[LINEPOINTS],   topedgedisky[LINEPOINTS],
          inneredgediskx[LINEPOINTS], inneredgedisky[LINEPOINTS];
   double inclination, x, y, z, lsky[4], msky[4], nsky[4],
               originx, coordsign, transmission;
   long int keepside, itheta, iphi, ir, nedge1, nedge2,
               max, nouteredgedisk, ntopedgedisk, ninneredgedisk;
   long int maxthetaindex1, maxphiindex1, maxthetaindex2, maxphiindex2,
              maxrindex, maxphiindexdisk;
   double    visible1 [SIZEDIMEN1] [SIZEDIMEN2],
             visible2 [SIZEDIMEN1] [SIZEDIMEN2],
          visibledisk [SIZEDIMEN1] [SIZEDIMEN2];

   /********************************
   *
   *  Project all the objects onto the plane of the sky and put the
   *  projected images into the coordinate system centered on star 2.
   *  Then find the edges of all the active objects.  Set visible[][]
   *  equal to zero for those points beyond the limbs of the objects.
   *
   *********************************/   
   if( (orbphase < -0.5) || (orbphase >0.5) )
      quit("Orbital phase must lie between -0.5 and 0.5 in totFlux.");
   inclination = starsys.inclination;
   coordsys( lsky, msky, nsky, inclination, orbphase);

   nedge1 = 0;
   nedge2 = 0;
   ntopedgedisk = 0;
   ninneredgedisk = 0;
   nouteredgedisk = 0;
   if( strcmp( mode.disk, "ON") == 0 ) {
      maxrindex = disktop.rpoints - 1;
      maxphiindexdisk = disktop.phipoints - 1;
      coordsign = -1.0;
      originx = 1.0;
      for( ir = 0; ir <= maxrindex; ir++) {
         for( iphi = 0; iphi <= maxphiindexdisk; iphi++) {
            x = originx + coordsign * rdisk[ir][iphi] * ldisk[ir][iphi];
            y =           coordsign * rdisk[ir][iphi] * mdisk[ir][iphi];
            z =                       rdisk[ir][iphi] * ndisk[ir][iphi];
            project( &xdisk[ir][iphi], &ydisk[ir][iphi],
                                           x, y, z, lsky, msky);
            mudisk[ir][iphi] =     coordsign * normdiskx[ir][iphi] * nsky[1]
                                 + coordsign * normdisky[ir][iphi] * nsky[2]
                                 +             normdiskz[ir][iphi] * nsky[3];
            visibledisk[ir][iphi] = 1.0;
            if( mudisk[ir][iphi] < 0.0 )
                   visibledisk[ir][iphi] = 0.0;
         }
      }
      nouteredgedisk = diskedges( 1, maxrindex, maxphiindexdisk,
                                           xdisk, ydisk, mudisk,
                                           outeredgediskx, outeredgedisky);
      if( (90.0 - starsys.inclination) < disktop.flareangle )
         ntopedgedisk   = diskedges( 2, maxrindex, maxphiindexdisk,
                                              xdisk, ydisk, mudisk,
                                              topedgediskx, topedgedisky);
      if( starsys.inclination > 0.1 )
         ninneredgedisk = diskedges( 3, maxrindex, maxphiindexdisk,
                                              xdisk, ydisk, mudisk,
                                              inneredgediskx, inneredgedisky);
   }

   if( strcmp( mode.star1, "ON") == 0 ) {
      maxthetaindex1 = star1.thetapoints - 1;
      maxphiindex1 = star1.phipoints - 1;
      coordsign = -1.0;
      originx = 1.0;
      for( itheta = 0; itheta <= maxthetaindex1; itheta++) {
         for( iphi = 0; iphi <= maxphiindex1; iphi++) {
            x = originx + coordsign * r1[itheta][iphi] * l1[itheta][iphi];
            y =           coordsign * r1[itheta][iphi] * m1[itheta][iphi];
            z =                       r1[itheta][iphi] * n1[itheta][iphi];
            project( &x1[itheta][iphi], &y1[itheta][iphi],
                         x, y, z, lsky, msky);
            mu1[itheta][iphi] =   coordsign * normalx1[itheta][iphi] * nsky[1]
                                + coordsign * normaly1[itheta][iphi] * nsky[2]
                                +             normalz1[itheta][iphi] * nsky[3];
            visible1[itheta][iphi] = 1.0;
            if( mu1[itheta][iphi] < 0.0 )
                   visible1[itheta][iphi] = 0.0;
         }
      }
      nedge1 = staredge( star1.thetapoints, star1.phipoints,
                                    x1, y1, mu1, edge1x, edge1y);
      cutoffstar1( ninneredgedisk, inneredgediskx, inneredgedisky,
                      nedge1, edge1x, edge1y);
      cutoffstar1( ntopedgedisk, topedgediskx, topedgedisky,
                      nedge1, edge1x, edge1y);
   }

   if( strcmp( mode.star2, "ON") == 0 ) {
      maxthetaindex2 = star2.thetapoints - 1;
      maxphiindex2 = star2.phipoints - 1;
      coordsign = 1.0;
      originx = 0.0;
      for( itheta = 0; itheta <= maxthetaindex2; itheta++) {
         for( iphi = 0; iphi <= maxphiindex2; iphi++) {
            x = originx + coordsign * r2[itheta][iphi] * l2[itheta][iphi];
            y =           coordsign * r2[itheta][iphi] * m2[itheta][iphi];
            z =                       r2[itheta][iphi] * n2[itheta][iphi];
            project( &x2[itheta][iphi], &y2[itheta][iphi],
                         x, y, z, lsky, msky);
            mu2[itheta][iphi] =   coordsign * normalx2[itheta][iphi] * nsky[1]
                                + coordsign * normaly2[itheta][iphi] * nsky[2]
                                +             normalz2[itheta][iphi] * nsky[3];
            visible2[itheta][iphi] = 1.0;
            if( mu2[itheta][iphi] < 0.0 )
                    visible2[itheta][iphi] = 0.0;
         }
      }
      nedge2 = staredge( star2.thetapoints, star2.phipoints,
                              x2, y2, mu2, edge2x, edge2y);
   }

   /****************************
   *
   *  Now find which points are occulted by other objects in 
   *  the system.
   *
   ******************************/
   /*************************************
   *
   *  First occult star 2 by star 1 and the disk.  Note that this
   *  version does not handle a disk with a transparent center correctly.
   *
   **************************************/
   if( strcmp( mode.star2, "ON") == 0 ) {
      if( fabs( orbphase ) > 0.25 ) {
      	 if( strcmp( mode.star1, "ON") == 0 ) {
      	    keepside = 1;
            transmission = 0.0;
            occult( keepside, transmission,
                    maxthetaindex2, maxphiindex2, x2, y2, visible2,
                    nedge1, edge1x, edge1y);
         }
      	 if( strcmp( mode.disk, "ON") == 0 ) {
      	    keepside = 1;
      	    transmission = disktop.transmission;
            occult( keepside, transmission,
                    maxthetaindex2, maxphiindex2, x2, y2, visible2,
                    nouteredgedisk, outeredgediskx, outeredgedisky);
         }
      }
   }
      
   /**************************
   *
   *  Then the disk.
   *
   ****************************/
   /********************
   *
   *  Treat the edge of the disk and the surface of the disk
   *  as if they are separate entities.
   *  The difference between the top of the disk and the edge is
   *  that the top can be occulted by the top edge of the disk
   *  and by star 1.
   *
   **********************/
   if( strcmp( mode.disk, "ON") == 0 ) {
      /*************************************
      *
      *  First do the edge of the disk.  The exchange and then
      *  re-exchange of elements is done to avoid writing yet
      *  another special function to handle the edge.
      *
      **************************************/
      if( strcmp( mode.star2, "ON") == 0 ) {
         if( fabs( orbphase ) <= 0.25 ) {
            max = 3;
            for( ir = 0; ir < max; ir++) {
               for( iphi = 0; iphi <= maxphiindexdisk; iphi++) {
                  z = xdisk[ir][iphi];
                  xdisk[ir][iphi] = xdisk[maxrindex - ir][iphi];
                  xdisk[maxrindex - ir][iphi] = z;
                  z = ydisk[ir][iphi];
                  ydisk[ir][iphi] = ydisk[maxrindex - ir][iphi];
                  ydisk[maxrindex - ir][iphi] = z;
                  z = visibledisk[ir][iphi];
                  visibledisk[ir][iphi] = visibledisk[maxrindex - ir][iphi];
                  visibledisk[maxrindex - ir][iphi] = z;
               }
            }
            keepside = 1;
            transmission = 0.0;
            occult( keepside, transmission,
                    max, maxphiindexdisk, xdisk, ydisk, visibledisk,
                    nedge2, edge2x, edge2y);
            for( ir = 0; ir < max; ir++) {
               for( iphi = 0; iphi <= maxphiindexdisk; iphi++) {
                  z = xdisk[ir][iphi];
                  xdisk[ir][iphi] = xdisk[maxrindex - ir][iphi];
                  xdisk[maxrindex - ir][iphi] = z;
                  z = ydisk[ir][iphi];
                  ydisk[ir][iphi] = ydisk[maxrindex - ir][iphi];
                  ydisk[maxrindex - ir][iphi] = z;
                  z = visibledisk[ir][iphi];
                  visibledisk[ir][iphi] = visibledisk[maxrindex - ir][iphi];
                  visibledisk[maxrindex - ir][iphi] = z;
               }
            }
         }
      }      
     
      /********************************
      *
      *  Now do the top of the disk.
      *
      ****************************/
      max = maxrindex - 3;
      if( strcmp( mode.star1, "ON") == 0 ) {
         keepside = 1;
         transmission = 0.0;
         occult( keepside, transmission,
                 max, maxphiindexdisk, xdisk, ydisk, visibledisk,
                 nedge1, edge1x, edge1y);
      }
      keepside = 0;
      transmission = disktop.transmission;
      occult( keepside, transmission,
              max, maxphiindexdisk, xdisk, ydisk, visibledisk,
              ntopedgedisk, topedgediskx, topedgedisky);
      if( strcmp( mode.star2, "ON") == 0 ) {
      	 if( fabs( orbphase ) <= 0.25 ) {
            keepside = 1;
            transmission = 0.0;
            occult( keepside, transmission,
                    max, maxphiindexdisk, xdisk, ydisk, visibledisk,
                    nedge2, edge2x, edge2y);
         }
      }      
   }

   /********************************
   *
   *  Finally, do star 1.
   *
   *********************************/
   if( strcmp( mode.star1, "ON" ) == 0 ) {
      keepside = 0;
      transmission = 0.0;
      occult( keepside, transmission,
              maxthetaindex1, maxphiindex1, x1, y1, visible1,
              nedge1, edge1x, edge1y);      
      if( strcmp( mode.star2, "ON") == 0 ) {
         if( fabs( orbphase ) <= 0.25 ) {
            keepside = 1;
            transmission = 0.0;
            occult( keepside, transmission,
                    maxthetaindex1, maxphiindex1, x1, y1, visible1,
                    nedge2, edge2x, edge2y);      
         }
      }
   }

   /*********************************************
   *
   *   Plot all points that have visibility > 0.0.
   *
   ***********************************************/
   strcpy( plotfile, "points.dat");
   if( ( out = fopen( plotfile, "w")) == NULL )
      quit("Cannot open file points.dat for writing the point plot.");
      
   if( strcmp( mode.star2, "ON") == 0 ) {
      for( itheta = 0; itheta <= maxthetaindex2; itheta++) {
         for( iphi = 0; iphi <= maxphiindex2; iphi++) {
            if( visible2[itheta][iphi] > 0.0 ) {
               fprintf( out," %3ld %3ld  %8.5lf %8.5lf\n",
                        itheta, iphi, x2[itheta][iphi], y2[itheta][iphi]);
            }
         }
      }
   }
   if( strcmp( mode.star1, "ON") == 0 ) {
      for( itheta = 0; itheta <= maxthetaindex1; itheta++) {
         for( iphi = 0; iphi <= maxphiindex1; iphi++) {
            if( visible1[itheta][iphi] > 0.0 ) {
               fprintf( out," %3ld %3ld  %8.5lf %8.5lf\n",
                        itheta, iphi, x1[itheta][iphi], y1[itheta][iphi]);
            }
         }
      }
   }
   if( strcmp( mode.disk, "ON") == 0 ) {
      for( ir = 0; ir <= maxrindex; ir++) {
         for( iphi = 0; iphi <= maxphiindexdisk; iphi++) {
            if( visibledisk[ir][iphi] > 0.0 ) {
               fprintf( out," %3ld %3ld  %8.5lf %8.5lf\n",
                    ir, iphi, xdisk[ir][iphi], ydisk[ir][iphi]);
            }
         }
      }
   }
   fclose( out );
   return;
}