/********************************
*
*                 FILE   EDGES.C
*
*   The functions in this file have to do with edges,
*      --  finding edges of objects
*      --  deciding whether points or lines are inside or outside
*               an edge.
*
*********************************/

#include "header.h"

long int staredge( long ntheta, long nphi, 
                   double  x[SIZEDIMEN1][SIZEDIMEN2], 
                   double  y[SIZEDIMEN1][SIZEDIMEN2],
                   double mu[SIZEDIMEN1][SIZEDIMEN2],
                   double *edgex, double *edgey)
/***************************
*
*  This function returns a set of points defining the edge of 
*  a three dimensional star that has been projected onto the
*  two dimensional plane of the sky.
*  Input data:
*     x[][], y[][]     = arrays holding the projected (ordered) points
*                           on the surface of the star.
*     mu[][]           = the angle between the normal to the surface of
*                           the star and the direction to the earth.
*     ntheta,nphi      = the sizes of the arrays.
*  Output data:
*     edgex[], edgey[] = the edge of the projected star.  The points
*                           are ordered, so that drawing lines between
*                           successive points outlines the star.
*     negdepoints      = the number of points in the arrays defining the
*                           edge curve.
*
************************************/
{
   double twopi, xlow, xhigh, ylow, yhigh, centerx, centery, edgemu, xa, ya,
               deltar, edger[LINEPOINTS], edgetheta[LINEPOINTS];
   long int maxthetaindex, maxphiindex, i, j, jminus,
               nedgepoints, extrapoints;

   edgemu = -1.0e-8;
   twopi = 2.0 * 3.1415926535897932385;
   maxthetaindex = ntheta -1;
   maxphiindex = nphi -1;
   nedgepoints = 0;

   /*************************
   *
   *  First find a preliminary (but good) guess at which points define
   *  the projected edge of the star by looking for places where mu changes
   *  sign.  Note that the curves in the phi direction are assumed to be
   *  closed.
   *
   *************************/
   for( i = 0; i <= maxthetaindex; i++) {
      for( j = 0; j <= maxphiindex; j++) {
      	 jminus = j - 1;
      	 if( j == 0 )
      	    jminus = maxphiindex;
         if( (mu[i][j] >= edgemu) && (mu[i][jminus] < edgemu) ) {
            edgex[nedgepoints] = x[i][j];
            edgey[nedgepoints] = y[i][j];
            nedgepoints = nedgepoints + 1;
         }
         if( (mu[i][j] < edgemu)  && (mu[i][jminus] >= edgemu) ) {
            edgex[nedgepoints] = x[i][jminus];
            edgey[nedgepoints] = y[i][jminus];
            nedgepoints = nedgepoints + 1;
         }
      }
   }
   for( j = 0; j <= maxphiindex; j++) {
      for( i = 1; i <= maxthetaindex; i++) {
         if( (mu[i][j] >= edgemu) && (mu[i-1][j] <  edgemu) ) {
            edgex[nedgepoints] = x[i][j];
            edgey[nedgepoints] = y[i][j];
            nedgepoints = nedgepoints + 1;
         }
         else if( (mu[i][j] <  edgemu) && (mu[i-1][j] >= edgemu) ) {
            edgex[nedgepoints] = x[i-1][j];
            edgey[nedgepoints] = y[i-1][j];
            nedgepoints = nedgepoints + 1;
         }
      }
   }

   /*******************************
   *
   *  Now sort the points so that plotting them in order draws 
   *  the edge.  The edge is assumed to be single valued as seen
   *  from its midpoint.
   *
   ********************************/
   xlow  = edgex[0];
   xhigh = edgex[0];
   ylow  = edgey[0];
   yhigh = edgey[0];
   for( i = 0; i < nedgepoints; i++) {
      if( edgex[i] < xlow  ) xlow  = edgex[i];
      if( edgex[i] > xhigh ) xhigh = edgex[i];
      if( edgey[i] < ylow  ) ylow  = edgey[i];
      if( edgey[i] > yhigh ) yhigh = edgey[i];
   }
   centerx = (xlow + xhigh) / 2.0;
   centery = (ylow + yhigh) / 2.0;

   for( i = 0; i < nedgepoints; i++) {
      xa = edgex[i] - centerx;
      ya = edgey[i] - centery;
      edger[i] = sqrt( xa*xa + ya*ya );
      edgetheta[i] = 0.0;
      if( edger[i] != 0.0 )
          edgetheta[i] = acos( xa / edger[i] );
      if( ya < 0.0 )
          edgetheta[i]= twopi - edgetheta[i];
   }
   sort( nedgepoints, edgetheta, edger);

   /********************************
   *
   *  Now check for those few strange extra points that can
   *  lie outside the edge found so far.
   *  (yes, reader, this can happen).
   *
   *********************************/
   extrapoints = 0;
   for( i = 0; i <=maxthetaindex; i++) {
      for( j = 0; j <= maxthetaindex; j++) {
      	 if( mu[i][j] >= edgemu ) {
            xa = x[i][j] - centerx;
            ya = y[i][j] - centery;
               deltar = edgedistance( nedgepoints, edger, edgetheta, xa, ya);
            if( deltar > 0.0 ) {
               edgex[nedgepoints + extrapoints] = x[i][j];
               edgey[nedgepoints + extrapoints] = y[i][j];
               extrapoints = extrapoints + 1;
            }
         }
      }
   }
   nedgepoints = nedgepoints + extrapoints;

   /********************
   *
   *  Resort with the new points, if any.
   *
   ********************/
   for( i = 0; i < nedgepoints; i++) {
      xa = edgex[i] - centerx;
      ya = edgey[i] - centery;
      edger[i] = sqrt( xa*xa + ya*ya );
      edgetheta[i] = 0.0;
      if( edger[i] != 0.0 )
          edgetheta[i] = acos( xa / edger[i] );
      if( ya < 0.0 )
          edgetheta[i]= twopi - edgetheta[i];
   }
   sort( nedgepoints, edgetheta, edger);
   
   /*********************************
   *
   *  Finally, get rid of all the duplicate points.
   *
   **********************************/
   for( i = 1; i < nedgepoints; i++) {
      if ( (edger[i] == edger[i-1]) && (edgetheta[i] == edgetheta[i-1]) ) {
         for( j = i; j < nedgepoints; j++) {
            edger[j] = edger[j+1];
            edgetheta[j] = edgetheta[j+1];
         }
         nedgepoints = nedgepoints - 1;
      }
   }
   for( i = 0; i < nedgepoints; i++) {
      edgex[i] = edger[i] * cos( edgetheta[i] ) + centerx;
      edgey[i] = edger[i] * sin( edgetheta[i] ) + centery;
   }
   return( nedgepoints );
}


long int diskedges( long edgetype, long maxrindex, long maxphiindex, 
                    double  x [SIZEDIMEN1] [SIZEDIMEN2],
                    double  y [SIZEDIMEN1] [SIZEDIMEN2],
                    double mu [SIZEDIMEN1] [SIZEDIMEN2],
                    double *edgex, double *edgey)
/******************************
*
*   This function returns a set of points defining the edge
*   of the disk as it is projected onto the plane of the sky.
*   The edge is the inner or outer edge depending on the switch isouter.
*   Input data:
*      edgetype         =  1, find the projected outer edge of the disk
*                          2, find the projected top edge of the disk
*                          3, find the projected inner edge of the disk
*      x[][], y[][]     =  arrays holding the projected (ordered) points
*                             on the surface of the disk.
*      mu[][]           =  the angle between the normal to the surface of
*                             the original disk and the direction to the
*                             earth.
*      maxrindex,
*         maxphiindex   =  the maximum sizes of the indices on the arrays.
*
*   Output data:
*      edgex[], edgey[] =  the edge of the projected disk.  The points
*                             are ordered so that drawing lines between
*                             successive points outlines the object.
*      nedgepoints      =  the number of points in the edge arrays.
*
*******************************/
{
   double twopi, edgemu, xa, ya, maxy, xlow, xhigh, ylow, yhigh,
        centerx, centery, edger[LINEPOINTS], edgetheta[LINEPOINTS];
   long int nedgepoints, iphi, last, i, j, k;

   edgemu = -1.0e-7;
   twopi = 2.0 * 3.14159265358979323846;

   /*****************************
   *
   *   if edgetype == 1,  find the outer edge of the disk.
   *
   *****************************/
   if( edgetype == 1 ) {
      nedgepoints = 0;
      last = 0;
      for( iphi = 0; iphi <= maxphiindex; iphi++) {
      	 k = iphi;
      	 if( k > maxphiindex )
      	    k = 0;
         if( mu[maxrindex][k] < edgemu ) {
            edgex[nedgepoints] = x[maxrindex - 3][k];
            edgey[nedgepoints] = y[maxrindex - 3][k];
            nedgepoints = nedgepoints + 1;
            if( iphi > 0 ) {
               if( last < 0 ) {
                  edgex[nedgepoints] = x[maxrindex - 3][iphi-1];
                  edgey[nedgepoints] = y[maxrindex - 3][iphi-1];
                  nedgepoints = nedgepoints + 1;
               }
            }
            last = 1;
         }
         if( mu[maxrindex][k] >= edgemu ) {
            edgex[nedgepoints] = x[maxrindex][k];
            edgey[nedgepoints] = y[maxrindex][k];
            nedgepoints = nedgepoints + 1;             
            if( iphi > 0 ) {
               if( last > 0 ) {
                  edgex[nedgepoints] = x[maxrindex - 3][k];
                  edgey[nedgepoints] = y[maxrindex - 3][k];
                  nedgepoints = nedgepoints + 1;
               }
            }
            last = -1;
         }
      }
   }

   /**************************
   *
   *  If edgetype == 2,  find the edge defined by that part of the
   *  projected top edge of the disk that is closest to the earth.
   *  A few dummy points are added to make it a closed curve.
   *
   **************************/
   else if( edgetype == 2 ) {
      nedgepoints = 0;
      centerx = 0.0;
      maxy = 0.0;
      for( iphi = 0; iphi <= maxphiindex; iphi++) {
         if( mu[maxrindex-2][iphi] >= 0.0 ) {
            edgex[nedgepoints] = x[maxrindex-3][iphi];
            edgey[nedgepoints] = y[maxrindex-3][iphi];
            centerx = centerx + edgex[nedgepoints];
            if( edgey[nedgepoints] > maxy ) maxy = edgey[nedgepoints];
            nedgepoints = nedgepoints + 1;
         }
      }
      centerx = centerx / nedgepoints;
      edgex[nedgepoints]   = centerx + 10.0;
      edgey[nedgepoints]   = maxy;
      edgex[nedgepoints+1] = centerx;
      edgey[nedgepoints+1] = maxy + 5.0;
      edgex[nedgepoints+2] = centerx - 10.0;
      edgey[nedgepoints+2] = maxy;
      nedgepoints = nedgepoints + 3;
   }

   /**************************
   *
   *  If  edgetype == 3,  find the edge defined by the projected
   *  edge of that half of inner edge of the disk closest to the
   *  the earth, and a few dummy points inserted to make it a
   *  closed curve.
   *
   ***************************************/
   else if( edgetype == 3 ) {
      nedgepoints = 0;
      centerx = 0.0;
      maxy = 0.0;
      for( iphi = 0; iphi <= maxphiindex; iphi++) {
         if( mu[maxrindex][iphi] >= 0.0 ) {
            edgex[nedgepoints] = x[0][iphi];
            edgey[nedgepoints] = y[0][iphi];
            centerx = centerx + edgex[nedgepoints];
            if( edgey[nedgepoints] > maxy ) maxy = edgey[nedgepoints];
            nedgepoints = nedgepoints + 1;
         }
      }
      centerx = centerx / nedgepoints;
      edgex[nedgepoints]   = centerx + 10.0;
      edgey[nedgepoints]   = maxy;
      edgex[nedgepoints+1] = centerx;
      edgey[nedgepoints+1] = maxy + 1.0;
      edgex[nedgepoints+2] = centerx - 10.0;
      edgey[nedgepoints+2] = maxy;
      nedgepoints = nedgepoints + 3;
   }
   else {
      quit("Unrecoginized edgetype in function diskedges");
   }

   /*******************************
   *
   *  Now sort the points so that plotting them in order draws 
   *  the edge.  The edge is assumed to be single valued as seen
   *  from its midpoint.
   *
   ********************************/
   xlow  = edgex[0];
   xhigh = edgex[0];
   ylow  = edgey[0];
   yhigh = edgey[0];
   for( i = 1; i < nedgepoints; i++) {
      if( edgex[i] < xlow  ) xlow  = edgex[i];
      if( edgex[i] > xhigh ) xhigh = edgex[i];
      if( edgey[i] < ylow  ) ylow  = edgey[i];
      if( edgey[i] > yhigh ) yhigh = edgey[i];
   }
   centerx = (xlow + xhigh) / 2.0;
   centery = (ylow + yhigh) / 2.0;
   
   for( i = 0; i < nedgepoints; i++) {
      xa = edgex[i] - centerx;
      ya = edgey[i] - centery;
      edger[i] = sqrt( xa*xa + ya*ya );
      edgetheta[i] = 0.0;
      if( edger[i] != 0.0 )
          edgetheta[i] = acos( xa / edger[i] );
      if( ya < 0.0 )
          edgetheta[i]= twopi - edgetheta[i];
   }
   sort( nedgepoints, edgetheta, edger);

   /*********************************
   *
   *  Finally, get rid of the duplicate point.
   *
   **********************************/
   for( i = 1; i < nedgepoints; i++) {
      if ( (edger[i] == edger[i-1]) && (edgetheta[i] == edgetheta[i-1]) ) {
         for( j = i; j < nedgepoints; j++) {
            edger[j] = edger[j+1];
            edgetheta[j] = edgetheta[j+1];
         }
         nedgepoints = nedgepoints - 1;
      }
   }
   for( i = 0; i < nedgepoints; i++) {
      edgex[i] = edger[i] * cos( edgetheta[i] ) + centerx;
      edgey[i] = edger[i] * sin( edgetheta[i] ) + centery;
   }
   return( nedgepoints );
}


double edgedistance( long nedgepoints, double *edger, double *edgetheta, 
                     double x, double y)
/***************************
*
*  This function calculates the distance from point (x,y) to the 
*  the nearest edge of a closed curve described by the arrays
*  (edger[], edgetheta[]).  The points in the curve must be arranged
*  in order of increasing theta.
*
*  The origin of the (r,theta) coordinate system and the (x,y)
*  coordinate system must be the same and must be somewhere near
*  the center of the closed curve.
*
*  The distance is positive if the point is outside the
*  curve and negative if it is inside.
*
*  Great care was expended in optimizing this function because it
*  threatens to be the majority consumer of cpu cycles in the light
*  curve synthesis program.
*
*****************************/
{
   int kmax, kmin, kmid;
   double twopi, r, theta, costheta, sintheta;
   double lowerr, upperr, lowertheta, uppertheta;
   double xlower, ylower, xupper, yupper, deltax, deltay;
   double tprime, xcross, ycross, rcross, deltar;

   /********************************
   *
   *  First find the (r,theta) coordinates of the point to
   *  be examined.
   *
   ***********************************/
   twopi= 6.2831853072;
   r = sqrt( x*x + y*y );
   if( r == 0.0 ) {
      theta = 0.0;
      costheta = 1.0;
      sintheta = 0.0;
   }
   else {
      costheta = x / r;
      sintheta = y / r;
      theta = acos( costheta );
      if( y < 0.0 ) {
         theta = twopi - theta;
      }
   }

   /******************************
   *
   *  The following algorithm for deciding where theta crosses
   *  the array thetadum[] is more efficient than checking
   *  successive pairs of points in the array.
   *
   ***********************************/
   kmax = nedgepoints;
   kmin = -1;
   kmid = (kmax + kmin) / 2;
   while( (kmax - kmin) > 1 ) {
      if( theta > edgetheta[kmid] ) {
         kmin = kmid;
      }
      else {
         kmax = kmid;
      }
      kmid = (kmax + kmin) / 2;
   }
   if( kmin == -1 ) {
      lowerr = edger[nedgepoints-1];
      lowertheta = edgetheta[nedgepoints-1] - twopi;
   }
   else {
      lowerr = edger[kmin];
      lowertheta = edgetheta[kmin];
   }
   if( kmax == nedgepoints ) {
      upperr = edger[0];
      uppertheta = edgetheta[0] + twopi;
   }
   else {
      upperr = edger[kmax];
      uppertheta = edgetheta[kmax];
   }
   if( uppertheta > lowertheta ) {
      /********************
      *
      *  NOTE:  the following interpolation is put into
      *  a Cartesian coordinate system because the edges are
      *  straight lines in (x,y), not (r,theta).  It is parametrized
      *  because the line can have an infinite slope.  The equation
      *  for tprime is not obvious, but it is correct.
      *
      *********************/
      xlower = lowerr * cos( lowertheta );
      ylower = lowerr * sin( lowertheta );
      xupper = upperr * cos( uppertheta );
      yupper = upperr * sin( uppertheta );
      deltax = xupper - xlower;
      deltay = yupper - ylower;
      tprime = ( ylower*costheta - xlower*sintheta )
                         / ( deltax*sintheta - deltay*costheta );
      xcross = tprime*deltax + xlower;
      ycross = tprime*deltay + ylower;
      rcross = sqrt( xcross*xcross + ycross*ycross );
   }
   else {
      rcross = 0.5 * ( upperr + lowerr );
   }
   deltar = r - rcross;
   return( deltar );
}


void sort( int npoints, double *x, double *y)
/*****************
*
*  Sort the points pairs (x,y) in the arrays x and y by 
*  increasing x.
*
******************/
{
   long int i, j, jsmall;
   double small, dummy;

   for( i = 0; i < (npoints -1); i++) {
      small = x[i];
      jsmall = i;
      for( j = i+1; j < npoints; j++) {
         if ( x[j] < small ) {
            small = x[j];
            jsmall = j;
         }
      }
      x[jsmall] = x[i];
      x[i] = small;
      dummy = y[jsmall];
      y[jsmall] = y[i];
      y[i] = dummy;
   }
   return;
}


long int cutoff( long keepside, long nedge, double *edgex, double *edgey,
                 long isclosed, long npoints, 
                 double *x, double *y, double *visible)
/****************************
*
*  This function calculates what part of a line is visible
*  and what part has been cut off by a closed curve.
*  Input data:
*     keepside         =  0, points on the inside of the edge are visible
*                      =  1, points on the outside of the edge are visible
*     edgex[], edgey[] =  array of points describing the edge of the
*                            occulting curve.  The occulting curve is
*                            assumed to be closed.
*     nedge            =  number of points in the occulting curve.
*     isclosed         =  1, if the line is closed.
*                         0, if the line is open.
*     x[], y[]         =  array of points describing the line.
*     npoints          =  the number of points in the line.
*     visible[]        =  1.0, if the point is visible
*                      =  0.0, if the point is invisible.
*
*  Output:
*      1.  The function inserts extra points in the line at the positiions
*          where the line crosses the occulting edge. 
*      2.  The array visible[] is set to 0.0 for points made invisible by
*          the occulting edge and to 1.0 elsewhere, including 1.0 at the
*          occulting edge itself.
*      3.  The function returns the revised (increased) number of points in
*          the line.
*      
****************************/
{
   double centerx, centery, xa, ya, edger[LINEPOINTS], edgetheta[LINEPOINTS];
   double twopi, foundedge, deltar[LINEPOINTS];
   double xlow, xhigh, ylow, yhigh, xlimb[LINEPOINTS], ylimb[LINEPOINTS];
   double xplus, xminus, yplus, yminus, newdeltar;
   long int i, j, k, limit,
              limbpoints, place[LINEPOINTS];

   if( nedge == 0 ) return( npoints );
              
   /*****************************
   *
   *  First put the edge points into a polar coordinate system
   *  with its origin near the center of the figure described
   *  by the edge.  Then sort the points in order of increasing
   *  theta.
   *
   **********************************/
   twopi = 2.0 * 3.1415926535897932385;
   xlow  = edgex[0];
   xhigh = edgex[0];
   ylow  = edgey[0];
   yhigh = edgey[0];
   for( i = 1; i < nedge; i++) {
      if( edgex[i] < xlow  ) xlow  = edgex[i];
      if( edgex[i] > xhigh ) xhigh = edgex[i];
      if( edgey[i] < ylow  ) ylow  = edgey[i];
      if( edgey[i] > yhigh ) yhigh = edgey[i];
   }
   centerx = (xlow + xhigh) / 2.0;
   centery = (ylow + yhigh) / 2.0;
   for( i = 0; i < nedge; i++) {
      xa = edgex[i] - centerx;
      ya = edgey[i] - centery;
      edger[i] = sqrt( xa*xa + ya*ya );
      edgetheta[i] = 0.0;
      if( edger[i] != 0.0 )
         edgetheta[i] = acos( xa /edger[i] );
      if( ya < 0.0 )
         edgetheta[i] = twopi - edgetheta[i];
   }
   sort( nedge, edgetheta, edger );

   /*********************************
   *
   *  Calculate the distance of each point in the line
   *  from the edge of the occulting curve.
   *  Negative is inside the curve, positive outside.
   *
   **************************************/
   for( i = 0; i < npoints; i++ ) {
      xa = x[i] - centerx;
      ya = y[i] - centery;
      deltar[i] = edgedistance( nedge, edger, edgetheta, xa, ya);
      if( (keepside == 1) && (deltar[i] < 0.0) )
         visible[i] = 0.0;
      if( (keepside == 0) && (deltar[i] >= 0.0) )
         visible[i] = 0.0;
   }
   if( isclosed == 1 ) {
      limit = npoints;
      x[npoints] = x[0];
      y[npoints] = y[0];
      deltar[npoints] = deltar[0];
      visible[npoints] = visible[0];
   }
   else 
      limit = npoints - 1;

   /*********************************
   *
   *  The following code for interpolating the position of the
   *  limb was chosen because it is failsafe, even near sharp
   *  corners, and is guaranteed to converge, albeit slowly.
   *  An edge is defined as a place where
   *     1)  the points go from visibile to invisible.
   *     2)  the distance of points from the edge go from
   *          positive to negative.
   *  The edge crossing is found by repetitively cutting the
   *  line segment with an edge crossing in half, and looking at
   *  the half that crosses the edge.
   *
   **********************************/
   limbpoints = 0;
   for( i = 1; i <= limit; i++) {
      if( fabs( visible[i] - visible[i-1]) > 0.001 ) {
         foundedge = 0;
         if( (deltar[i] < 0.0) && (deltar[i-1] > 0.0) ) {
            xplus  = x[i-1] - centerx;
            yplus  = y[i-1] - centery;
            xminus = x[i]   - centerx;
            yminus = y[i]   - centery;
            foundedge = 1;
         }
         if( (deltar[i] > 0.0) && (deltar[i-1] < 0.0) ) {
            xplus  = x[i]   - centerx;
            yplus  = y[i]   - centery;
            xminus = x[i-1] - centerx;
            yminus = y[i-1] - centery;
            foundedge = 1;
         }
         if( foundedge == 1 ) {
            limbpoints = limbpoints + 1;
            for( j = 0; j <= 5; j++) {
               xa = 0.5 * ( xplus + xminus );
               ya = 0.5 * ( yplus + yminus );
               newdeltar = edgedistance( nedge, edger, edgetheta, xa, ya);
               if( newdeltar >= 0.0 ) {
                  xplus = xa;
                  yplus = ya;
               }
               else {
                  xminus = xa;
                  yminus = ya;
               }
            }
            xlimb[limbpoints] = 0.5 * (xplus + xminus) + centerx;
            ylimb[limbpoints] = 0.5 * (yplus + yminus) + centery;
            place[limbpoints] = i - 1;
         }
      }
   }

   /******************************
   *
   *  Insert crossing points into the line at the appropriate
   *  places by changing the entries in array visible[].
   *
   ********************************/
   for( i = limbpoints; i > 0; i--) {
      j = place[i];
      for( k = (npoints-1); k > j; k--) {
         x[k+1] = x[k];
         y[k+1] = y[k];
         visible[k+1] = visible[k];
      }
      npoints = npoints + 1;
      x[j+1] = xlimb[i];
      y[j+1] = ylimb[i];
      visible[j+1] = 1.0;
   }
   return( npoints );
}


long int cutoffstar1( long nedge, double *edgex, double *edgey,
                      long npoints, double *x, double *y)
/****************************
*
*  This function modifies the edge of star 1 by cutting it off at either
*  the inner or the outer edge of the accretion disk.  The cut off points
*  are replaced by new points on the cut off edge of the star.
*  Input data:
*     edgex[], edgey[] =  array of points describing the inner or outer edge
*                            of the disk.  The curve is assumed to be closed.
*     nedge            =  number of points in the arrays describing the
*                            inner edge of the disk.
*     x[], y[]         =  array of points describing the edge of star 1.
*     npoints          =  the number of points in the arrays describing
*                            edge of star 1.
*
*  Output:
*      1.  The revised edge of star 1.
*      
****************************/
{
   double centerx, centery, xa, ya, twopi, r, theta, xlow, xhigh, deltar,
                   ylow, yhigh, epsilonr;
   double edger[LINEPOINTS], edgetheta[LINEPOINTS];
   long int i;

   if( nedge == 0 ) return( npoints );
   epsilonr = 1.0e-7;

   /*****************************
   *
   *  First put the edge points into a polar coordinate system
   *  with its origin near the center of the figure described
   *  by the edge.  Then sort the points in order of increasing
   *  theta.
   *
   **********************************/
   twopi = 2.0 * 3.1415926535897932385;
   xlow  = edgex[0];
   xhigh = edgex[0];
   ylow  = edgey[0];
   yhigh = edgey[0];
   for( i = 1; i < nedge; i++) {
      if( edgex[i] < xlow  ) xlow  = edgex[i];
      if( edgex[i] > xhigh ) xhigh = edgex[i];
      if( edgey[i] < ylow  ) ylow  = edgey[i];
      if( edgey[i] > yhigh ) yhigh = edgey[i];
   }
   centerx = (xlow + xhigh) / 2.0;
   centery = (ylow + yhigh) / 2.0;

   for( i = 0; i < nedge; i++) {
      xa = edgex[i] - centerx;
      ya = edgey[i] - centery;
      edger[i] = sqrt( xa*xa + ya*ya );
      edgetheta[i] = 0.0;
      if( edger[i] != 0.0 )
         edgetheta[i] = acos( xa /edger[i] );
      if( ya < 0.0 )
         edgetheta[i] = twopi - edgetheta[i];
   }
   sort( nedge, edgetheta, edger );
 
   /*****************************
   *
   *  Find the distance from each point on the edge of star1
   *  to the edge of the surrounding occulting object.  If the
   *  distance is positive, replace the point on the edge of star1
   *  inwards at constant theta to the edge of the occulting object.
   *
   ***********************************/
   for( i = 0; i <= npoints; i++ ) {
      xa = x[i] - centerx;
      ya = y[i] - centery;
      deltar = edgedistance( nedge, edger, edgetheta, xa, ya);
      if( deltar > epsilonr ) {
         r = sqrt( xa*xa + ya*ya);
         theta = 0.0;
         if( r != 0.0 )
            theta = acos( xa / r );
         if( ya < 0.0 )
            theta = twopi - theta;
         x[i] = (r - deltar) * cos( theta ) + centerx;
         y[i] = (r - deltar) * sin( theta ) + centery;
      }
   }

   /**************************
   *
   *  Resort the points along the edge of star1 because the
   *  origin has changed.
   *
   ***************************/
   xlow  = x[0];
   xhigh = x[0];
   ylow  = y[0];
   yhigh = y[0];
   for( i = 1; i < npoints; i++) {
      if( x[i] < xlow  ) xlow  = x[i];
      if( x[i] > xhigh ) xhigh = x[i];
      if( y[i] < ylow  ) ylow  = y[i];
      if( y[i] > yhigh ) yhigh = y[i];
   }
   centerx = (xlow + xhigh) / 2.0;
   centery = (ylow + yhigh) / 2.0;

   for( i = 0; i < npoints; i++) {
      xa = x[i] - centerx;
      ya = y[i] - centery;
      edger[i] = sqrt( xa*xa + ya*ya );
      edgetheta[i] = 0.0;
      if( edger[i] != 0.0 )
         edgetheta[i] = acos( xa /edger[i] );
      if( ya < 0.0 )
         edgetheta[i] = twopi - edgetheta[i];
   }
   sort( npoints, edgetheta, edger );

   for( i = 0; i < npoints; i++) {
         x[i] = edger[i] * cos( edgetheta[i] ) + centerx;
         y[i] = edger[i] * sin( edgetheta[i] ) + centery;
   }
   
   return( 0 );
}