/*******************************************
*
*              FILE   LINES.C
*
*  This file contains the functions unique to the creation of a line
*  image of the binary system. 
*
*  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;
*               the names of the output files have the orbital phase
*               incorporated into their names, but in this case the
*               phase is between 0.0 and 1.0
*
**********************************************/

#include "header.h"

long LineNumber;

void lineImage( double orbphase )
/*************************
*
*  Make a picture of the binary at one orbital phase by drawing
*  the grid describing its geometry.
*
**************************/
{
   FILE *out;
   char plotfile[40], s[10], color[20];
   double centerofmass, edgemu, inclination;
   double x, y, z, lsky[4], msky[4], nsky[4];
   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          linex[LINEPOINTS],          liney[LINEPOINTS],
                 visible[LINEPOINTS], originx, coordsign;
   long dummyint;
   long int isclosed, keepside, nline;
   long int i, itheta, iphi, ir, nedge1, nedge2,
               nouteredgedisk, ntopedgedisk, ninneredgedisk;
   long int maxthetaindex1, maxphiindex1, maxthetaindex2, maxphiindex2,
              maxrindex, maxphiindexdisk;

   /********************************
   *
   *  A few preliminaries.  Note that the output file has the
   *  orbital phase encoded in its name.
   *
   *********************************/
   centerofmass = 1.0 / ( 1.0 + starsys.q);
   x = orbphase;
   if( x < 0.0 ) 
      x = x + 1.0;
   dummyint = 10000 + 1000 * x;
   sprintf( s, "%-ld", dummyint);
   strcpy( plotfile, "image");
   strcat( plotfile, s);
   strcat( plotfile, ".mac");
   if( ( out = fopen( plotfile, "w")) == NULL )
      quit("Cannot open that file.");

   inclination = starsys.inclination;
   coordsys( lsky, msky, nsky, inclination, orbphase);
   edgemu = -1.0e-7;
   LineNumber = 0;

   /********************************
   *
   *  Project all the objects onto the plane of the sky and put the
   *  projected images into the coordinate system centered on 
   *  the center of mass.
   *  Then find the edges of all the active objects.
   *
   *********************************/
   nedge1 = 0;
   nedge2 = 0;
   nouteredgedisk = 0;
   ntopedgedisk = 0; 
   ninneredgedisk = 0;
   if( strcmp( mode.disk, "ON") == 0 ) {
      maxrindex = disktop.rpoints - 1;
      maxphiindexdisk = disktop.phipoints - 1;
      coordsign = -1.0;
      originx = centerofmass;
      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];
         }
      }
      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 = centerofmass;
      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];
         }
      }
      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 = centerofmass - 1.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];
         }
      }
      nedge2 = staredge( star2.thetapoints, star2.phipoints,
                              x2, y2, mu2, edge2x, edge2y);
   }

   /****************************
   *
   *  Now plot the image of the system.
   *
   ******************************/
   /*************************************
   *
   *  First plot star 2.
   *
   **************************************/
   if( strcmp( mode.star2, "ON") == 0 ) {
      strcpy( color, "red");
      isclosed = 1;
      for( i = 0; i < nedge2; i++) {
         linex[i] = edge2x[i];
         liney[i] = edge2y[i];
         visible[i] = 1.0;
      }
      nline = nedge2;
      if( fabs( orbphase ) > 0.25 ) {
         if( strcmp( mode.star1, "ON") == 0) {
            keepside = 1;
            nline = cutoff( keepside, nedge1, edge1x, edge1y,
                           isclosed, nline, linex, liney, visible);
         }
         if( strcmp( mode.disk, "ON") == 0) {
            keepside = 1;
            nline = cutoff( keepside,
                        nouteredgedisk, outeredgediskx, outeredgedisky,
                        isclosed, nline, linex, liney, visible);
         }
      }
      plotline( out, color, isclosed, nline, linex, liney, visible);
      for( itheta = 0; itheta <= maxthetaindex2; itheta++) {
         for( iphi = 0; iphi <= maxphiindex2; iphi++) {
            linex[iphi] = x2[itheta][iphi];
            liney[iphi] = y2[itheta][iphi];
            visible[iphi] = 1.0;
            if( mu2[itheta][iphi] < edgemu )
               visible[iphi] = 0.0;
         }            
         nline = maxphiindex2 + 1;
         isclosed = 1;
         if( fabs( orbphase ) > 0.25 ) {
            if( strcmp( mode.star1, "ON") == 0) {
               keepside = 1;
               nline = cutoff( keepside, nedge1, edge1x, edge1y,
                              isclosed, nline, linex, liney, visible);
            }
            if( strcmp( mode.disk, "ON") == 0) {
               keepside = 1;
               nline = cutoff( keepside,
                           nouteredgedisk, outeredgediskx, outeredgedisky,
                           isclosed, nline, linex, liney, visible);
            }
         }
         plotline( out, color, isclosed, nline, linex, liney, visible);
      }
      for( iphi = 0; iphi <= maxphiindex2; iphi++) {
         for( itheta = 0; itheta <= maxthetaindex2; itheta++) {
            linex[itheta] = x2[itheta][iphi];
            liney[itheta] = y2[itheta][iphi];
            visible[itheta] = 1.0;
            if( mu2[itheta][iphi] < edgemu )
               visible[itheta] = 0.0;
         }
         nline = maxthetaindex2 + 1;
         isclosed = 0;
         if( fabs( orbphase ) > 0.25 ) {
            if( strcmp( mode.star1, "ON") == 0) {
               keepside = 1;
               nline = cutoff( keepside, nedge1, edge1x, edge1y,
                              isclosed, nline, linex, liney, visible);
            }
            if( strcmp( mode.disk, "ON") == 0) {
               keepside = 1;
               nline = cutoff( keepside,
                           nouteredgedisk, outeredgediskx, outeredgedisky,
                           isclosed, nline, linex, liney, visible);
            }
         }
         plotline( out, color, isclosed, nline, linex, liney, visible);
      }
   }
      
   /**************************
   *
   *  Then plot the disk.
   *
   ****************************/
   /********************
   *
   *  Plot the edge of the disk and the surface of the disk
   *  as if they are separate entities.  First do the edge:
   *
   **********************/
   if( strcmp( mode.disk, "ON") == 0 ) {
      strcpy( color, "cyan");
      for( ir = (maxrindex-2); ir <= maxrindex; ir++) {
         for( iphi = 0; iphi <= maxphiindexdisk; iphi++) {
            linex[iphi] = xdisk[ir][iphi];
            liney[iphi] = ydisk[ir][iphi];
            visible[iphi] = 1.0;
            if( mudisk[ir][iphi] < edgemu )
               visible[iphi] = 0.0;            
         }
         nline = maxphiindexdisk + 1;
         isclosed = 1;
         if( fabs( orbphase ) <= 0.25 ) {
            if( strcmp( mode.star2, "ON") == 0) {
               keepside = 1;
               nline = cutoff( keepside, nedge2, edge2x, edge2y,
                              isclosed, nline, linex, liney, visible);
            }
         }
         plotline( out, color, isclosed, nline, linex, liney, visible);
      }
      for( iphi = 0; iphi <= maxphiindexdisk; iphi++) {
         for( ir = (maxrindex-2); ir <= maxrindex; ir++) {
            linex[ir - maxrindex + 2] = xdisk[ir][iphi];
            liney[ir - maxrindex + 2] = ydisk[ir][iphi];
            visible[ir - maxrindex + 2] = 1.0;
            if( mudisk[ir][iphi] < edgemu )
               visible[ir - maxrindex + 2] = 0.0;
         }
         nline = 3;
         isclosed = 0;
         if( fabs( orbphase ) <= 0.25 ) {
            if( strcmp( mode.star2, "ON") == 0) {
               keepside = 1;
               nline = cutoff( keepside, nedge2, edge2x, edge2y,
                           isclosed, nline, linex, liney, visible);
            }
         }
         plotline( out, color, isclosed, nline, linex, liney, visible);
      }

      /***************************
      *
      *  Now plot the top of the disk.  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.
      *
      ****************************/
      for( ir = 0; ir <= (maxrindex-3); ir++) {
         for( iphi = 0; iphi <= maxphiindexdisk; iphi++) {
            linex[iphi] = xdisk[ir][iphi];
            liney[iphi] = ydisk[ir][iphi];
            visible[iphi] = 1.0;
            if( mudisk[ir][iphi] < edgemu )
               visible[iphi] = 0.0;
         }            
         nline = maxphiindexdisk + 1;
         isclosed = 1;
         if( strcmp( mode.star1, "ON") == 0) {
            keepside = 1;
            nline = cutoff( keepside, nedge1, edge1x, edge1y,
                           isclosed, nline, linex, liney, visible);
         }
         keepside = 0;
         nline = cutoff( keepside, ntopedgedisk, topedgediskx, topedgedisky,
                           isclosed, nline, linex, liney, visible);
         if( fabs( orbphase ) <= 0.25 ) {
            if( strcmp( mode.star2, "ON") == 0) {
               keepside = 1;
               nline = cutoff( keepside, nedge2, edge2x, edge2y,
                              isclosed, nline, linex, liney, visible);
            }
         }
         plotline( out, color, isclosed, nline, linex, liney, visible);
      }
      for( iphi = 0; iphi <= maxphiindexdisk; iphi++) {
         for( ir = 0; ir <= (maxrindex-3); ir++) {
            linex[ir] = xdisk[ir][iphi];
            liney[ir] = ydisk[ir][iphi];
            visible[ir] = 1.0;
            if( mudisk[ir][iphi] < edgemu )
               visible[ir] = 0.0;
         }
         nline = maxrindex - 2;
         isclosed = 0;
         if( strcmp( mode.star1, "ON") == 0) {
            keepside = 1;
            nline = cutoff( keepside, nedge1, edge1x, edge1y,
                         isclosed, nline, linex, liney, visible);
         }
         keepside = 0;
         nline = cutoff( keepside, ntopedgedisk, topedgediskx, topedgedisky,
                           isclosed, nline, linex, liney, visible);
         if( fabs( orbphase ) <= 0.25 ) {
            if( strcmp( mode.star2, "ON") == 0) {
               keepside = 1;
               nline = cutoff( keepside, nedge2, edge2x, edge2y,
                           isclosed, nline, linex, liney, visible);
            }
         }
         plotline( out, color, isclosed, nline, linex, liney, visible);
      }
   }

   /********************************
   *
   *  Finally, plot star 1.
   *
   *********************************/
   if( strcmp( mode.star1, "ON" ) == 0 ) {
      strcpy( color, "black");
      isclosed = 1;
      for( i = 0; i < nedge1; i++) {
         linex[i] = edge1x[i];
         liney[i] = edge1y[i];
         visible[i] = 1.0;
      }
      nline = nedge1;
      if( fabs( orbphase ) <= 0.25 ) {
         if( strcmp( mode.star2, "ON") == 0) {
            keepside = 1;
            nline = cutoff( keepside, nedge2, edge2x, edge2y,
                              isclosed, nline, linex, liney, visible);
         }
      }
      plotline( out, color, isclosed, nline, linex, liney, visible);
      for( itheta = 0; itheta <= maxthetaindex1; itheta++) {
         for( iphi = 0; iphi <= maxphiindex1; iphi++) {
            linex[iphi] = x1[itheta][iphi];
            liney[iphi] = y1[itheta][iphi];
            visible[iphi] = 1.0;
            if( mu1[itheta][iphi] < edgemu )
               visible[iphi] = 0.0;
         }            
         nline = maxphiindex1 + 1;
         isclosed = 1;
         if( strcmp( mode.disk, "ON") == 0) {
            keepside = 0;
            nline = cutoff( keepside, nedge1, edge1x, edge1y,
                           isclosed, nline, linex, liney, visible);
         }
         if( fabs( orbphase ) <= 0.25 ) {
            if( strcmp( mode.star2, "ON") == 0) {
               keepside = 1;
               nline = cutoff( keepside, nedge2, edge2x, edge2y,
                              isclosed, nline, linex, liney, visible);
            }
         }
         plotline( out, color, isclosed, nline, linex, liney, visible);
      }
      for( iphi = 0; iphi <= maxphiindex1; iphi++) {
         for( itheta = 0; itheta <= maxthetaindex1; itheta++) {
            linex[itheta] = x1[itheta][iphi];
            liney[itheta] = y1[itheta][iphi];
            visible[itheta] = 1.0;
            if( mu1[itheta][iphi] < edgemu )
               visible[itheta] = 0.0;
         }
         nline = maxthetaindex1 + 1;
         isclosed = 0;
         if( strcmp( mode.disk, "ON") == 0) {
            keepside = 0;
            nline = cutoff( keepside, nedge1, edge1x, edge1y,
                              isclosed, nline, linex, liney, visible);
         }
         if( fabs( orbphase ) <= 0.25 ) {
            if( strcmp( mode.star2, "ON") == 0) {
               keepside = 1;
               nline = cutoff( keepside, nedge2, edge2x, edge2y,
                              isclosed, nline, linex, liney, visible);
            }
         }
         plotline( out, color, isclosed, nline, linex, liney, visible);
      }
   }

   nline = 0;
   strcpy( color, "black");
   plotline( out, color, isclosed, nline, linex, liney, visible);
   fclose( out );
   return;
}


void plotline( FILE *out, char *color, long isclosed, long npoints, 
                 double *x, double *y, double *visible)
/*********************************
*
*  This function creates an output file to be used as part of the
*  command file in a plotting program.  It plots a single line
*  traced out by the set of points (x,y).  Note that a trick has to
*  be used for supermongo because supermongo only accepts macros with
*  a rather short length.
*  Input data:
*     out       =  pointer pointing at the output file.
*     x, y      =  arrays containing the (x,y) positions of the points.
*     npoints   >  0, the number of data points in the (x,y) arrays.
*               =  0, close the file with the plot commands.
*     isclosed  =  1, if this is a closed curve.
*               =  0, if this is an open curve.
*     visible   =  array containing saying whether the points are to
*                     be plotted.  If visible[i] >= 0.5, point i is
*                     plotted; otherwise it is not. 
*
************************************/
{
   char movePenUp[20], movePenDown[20];
   long int i, isfirstpoint;

   if(      strcmp(mode.plotprogram,  "gle")   == 0) {
      strcpy( movePenUp,   "amove");
      strcpy( movePenDown, "aline");
   }
   else if( strcmp(mode.plotprogram, "mongo")  == 0) {
      strcpy( movePenUp,   "relocate");
      strcpy( movePenDown, "draw    ");
   }
   else if( strcmp(mode.plotprogram, "supermongo") == 0) {
      strcpy( movePenUp,   "relocate");
      strcpy( movePenDown, "draw    ");
   }
   else
      quit("Unrecognized plotting program in plotline().");

   /*********************************
   *
   *   Close the supemongo file by adding the global macro.
   *
   ************************************/
   if( npoints == 0) {
      if( strcmp(mode.plotprogram, "supermongo") == 0) {
        fprintf( out, "Plotall\n");
	for( i = 1; i <= LineNumber; i++ ) {
           fprintf( out, "   Line%-ld\n", i);
        }
      }
      fprintf( out, "   ctype black\n");
      return;
   }
   if( LineNumber == 0 ) {
      fprintf( out, "plotfig\n");
      fprintf( out, "   erase\n");
      fprintf( out, "   lweight 2\n");
      fprintf( out, "   limits -1.3 1.3 -1.3 1.3\n");
      fprintf( out, "   box 3 3 3 3\n");
      fprintf( out, "   Plotall\n");
   }
   isfirstpoint = 1;
   for( i = 0; i < npoints; i++) {
      if( visible[i] >= 0.5 ) {
         if( isfirstpoint == 1 ) {
            LineNumber += 1;
            fprintf( out, "Line%-ld\n", LineNumber);
            fprintf( out, "   ctype %s\n", color);     
            fprintf( out, "%11s %10.5lf  %10.5lf\n", movePenUp, x[i], y[i]);
            isfirstpoint = 0;
         }
         else {
            fprintf( out, "%11s %10.5lf  %10.5lf\n", movePenDown, x[i], y[i]);
         }
      }
      else {
         isfirstpoint = 1;
      }
   }

   if( isclosed == 1 )
      if( visible[0] >= 0.5 )
         if( isfirstpoint == 0 )
            fprintf( out, "%11s %10.5lf  %10.5lf\n", movePenDown, x[0], y[0]);

   return;
}