/********************************
*
*                    FILE     GEOMETRY.C
*
*  All the functions concerned with the basic geometry of the stars
*  and disk are in this file.
*
************************************/

#include "header.h"


void StarGeom ( struct starparams *star, 
           double       r [SIZEDIMEN1] [SIZEDIMEN2],
           double       l [SIZEDIMEN1] [SIZEDIMEN2],
           double       m [SIZEDIMEN1] [SIZEDIMEN2],
           double       n [SIZEDIMEN1] [SIZEDIMEN2],
           double gravity [SIZEDIMEN1] [SIZEDIMEN2],
           double normalx [SIZEDIMEN1] [SIZEDIMEN2],
           double normaly [SIZEDIMEN1] [SIZEDIMEN2],
           double normalz [SIZEDIMEN1] [SIZEDIMEN2],
           double      dS [SIZEDIMEN1] [SIZEDIMEN2]
         )
/*******************************************
*
*  This function is the main function that sets up the geometry
*  of the stars.
*
*******************************************/
{
   double x, y, z, rdum, ldum, mdum, ndum, rguess, pi, rsquared,
           mu2, theta, deltatheta, phi, deltaphi, cosbeta, sintheta;
   long int itheta, iphi, maxthetaindex, maxphiindex;

   maxthetaindex = star->thetapoints - 1;
   maxphiindex = star->phipoints - 1;
   pi = 3.14159265358979323;
   
   /*********************************
   *
   *  Calculate a few special numbers
   *
   **********************************/
   mu2 = star->massFrac;
   (*star).rL1 = findL1( mu2 );
      ldum =  1.0;
      ndum =  0.0;
   (*star).psiLobe = psi( star->rL1, ldum, ndum, mu2);
   (*star).frontradius = star->fracRadius * (star->rL1 - 1.0e-10);
      ldum =  1.0;
      ndum =  0.0;
   (*star).psiSurface = psi( star->frontradius, ldum, ndum, mu2);
      rguess = 0.9 * star->frontradius;
      ldum =  0.0;
      ndum =  0.0;
   (*star).sideradius = 
               findr( star->psiSurface, star->rL1, rguess, ldum, ndum, mu2);
      ldum = -1.0;
      ndum =  0.0;
   (*star).backradius = 
               findr( star->psiSurface, star->rL1, rguess, ldum, ndum, mu2);
      ldum =  0.0;
      ndum =  1.0;
   (*star).poleradius = 
               findr( star->psiSurface, star->rL1, rguess, ldum, ndum, mu2);

   /****************
   *
   *  Calculate the direction cosines at the grid points.
   *  Do not pile up points at theta = 0 or pi! 
   *  Do not double count the endpoints of the range!
   *
   *****************/

   deltatheta = pi / (maxthetaindex + 1.0);
   deltaphi   = 2.0 * pi / (maxphiindex +1.0);   
   for( itheta=0; itheta <= maxthetaindex; itheta++) {
      theta = deltatheta/2.0 + itheta * deltatheta;
      for( iphi=0; iphi <= maxphiindex; iphi++) {
      	 phi = iphi * deltaphi;
         l[itheta][iphi] = sin(theta) * cos(phi);
         m[itheta][iphi] = sin(theta) * sin(phi);
         n[itheta][iphi] = cos(theta);
      }
   }

   /***************
   *
   *  Calculate the value of r at the grid points.
   *
   ****************/

   ldum = 0.0;
   ndum = 0.0;
   rguess = 0.5 * star->rL1;
   rguess = findr( star->psiSurface, star->rL1, rguess, ldum, ndum, mu2);
   for( itheta=0; itheta <= maxthetaindex; itheta++) {
      for( iphi=0; iphi <= maxphiindex; iphi++) {
         ldum = l[itheta][iphi];
         ndum = n[itheta][iphi];
         if( (ldum > 0.999999) && (star->fracRadius > 0.99999) )
            r[itheta][iphi] = star->rL1;
         else
            r[itheta][iphi] =
              findr( star->psiSurface, star->rL1, rguess, ldum, ndum, mu2);
      }
   }

   /*************
   *
   *   Calculate the gravity and the surface normal at each of the
   *   grid points.  The gravity at the pole is set equal to the
   *   gravity on the z axis at the surface of the star.
   *   This is not exactly the gravity at the pole, but the
   *   difference is irrelevant.
   *
   *************/

   ldum = 0.0;
   mdum = 0.0;
   ndum = 1.0;
   rdum = findr( star->psiSurface, star->rL1, r[0][0], ldum, ndum, mu2);
   x = -dpsidx( rdum, ldum, mdum, ndum, mu2);
   y = -dpsidy( rdum, ldum, mdum, ndum, mu2);
   z = -dpsidz( rdum, ldum, mdum, ndum, mu2);
   (*star).polegrav = sqrt( x*x + y*y + z*z );
   for( itheta=0; itheta <= maxthetaindex; itheta++) {
      for( iphi=0; iphi <= maxphiindex; iphi++) {
         ldum = l[itheta][iphi];
         mdum = m[itheta][iphi];
         ndum = n[itheta][iphi];
         rdum = r[itheta][iphi];
         if( (ldum > 0.9999) && (star->fracRadius > 0.9995) ) {
            gravity[itheta][iphi] = 1.0e-8;
            normalx[itheta][iphi] = 1.0;
            normaly[itheta][iphi] = 0.0;
            normalz[itheta][iphi] = 0.0;
         }
         else {
            x = -dpsidx( rdum, ldum, mdum, ndum, mu2);
            y = -dpsidy( rdum, ldum, mdum, ndum, mu2);
            z = -dpsidz( rdum, ldum, mdum, ndum, mu2);
            gravity[itheta][iphi] = sqrt( x*x + y*y + z*z );
            normalx[itheta][iphi] = x / gravity[itheta][iphi];
            normaly[itheta][iphi] = y / gravity[itheta][iphi];
            normalz[itheta][iphi] = z / gravity[itheta][iphi];
         }
      }
   }

   /****************
   *
   *  And finally, calculate the scaler element of surface area dS2
   *
   *****************/

   for( itheta=0; itheta <= maxthetaindex; itheta++) {
      for( iphi=0; iphi <= maxphiindex; iphi++) {
         cosbeta =   normalx[itheta][iphi] * l[itheta][iphi]
                   + normaly[itheta][iphi] * m[itheta][iphi]
                   + normalz[itheta][iphi] * n[itheta][iphi];
         sintheta = sqrt( l[itheta][iphi] * l[itheta][iphi]
                        + m[itheta][iphi] * m[itheta][iphi] );
         rsquared = r[itheta][iphi] * r[itheta][iphi];
         dS[itheta][iphi] =
                     rsquared * sintheta * deltatheta * deltaphi / cosbeta;
      }
   }
   return;
}

/*******************
*                   
*  The following functions calculate various quantities in the
*  restricted three-body system.
*
*  Two equivalent coordinate systems are used in this set of
*  functions, a Cartesian coordinate system and a spherical
*  coordinate system.  The reason for using both is that some
*  calculations are easier in the spherical coordinate system,
*  some in the Cartesian.
*  The Cartesian coordinate system:
*    -- Position of a point is ( x, y, z)
*    -- Fixed in the rotating frame so that the stars are motionless
*           in the coordinate system.
*    -- Has its origin at the center of mass of a star, which is called
        star 2 throughout.
*            x : is along the line joining the centers of mass
*            y : in the orbital plane, perpendicular to x, in
*                  the direction opposite the orbital motion of star 2
*            z : perpendicular to orbital plane, making a right
*                  handed system.
*
*  The spherical coordinate system:
*    -- Position of a point is ( r, theta, phi)
*    -- Origin is the same as the origin of the Cartesian system.
*    -- Aligned with the Cartesian coordinates system such that
*           1)  theta = 0  is along the z axis
*           2)  theta = pi/2, phi = 0 is along the x axis
*
*  The conversion of the position of a point from the spherical to
*  the Cartesian system is given by
*       x = |r|l,  y = |r|m,   z = |r|n
*  Where the direction cosines are
*                 l = sin(theta) * cos(phi)
*                 m = sin(theta) * sin(phi)
*                 n = cos(theta)
*
*  The reverse conversion is given by
*      r**2 = (x**2 + y**2 + z**2)
*        tan(phi) = y / x
*        tan(theta)  = sqrt(x**2 + y**2) / z
*
*******************/        


double psi( double r, double l, double n, double mu2)
/************
*
*   This function calculates the potential of the zero velocity
*   surface at a point in space, and is the heart of the geometry
*   calculation.  The input data are:
*      r = distance from the center of star 2 to the point.
*      l = the direction cosine along the x direction
*      n = the direction cosine along the z direction
*      mu2 = the mass fraction = m2 / (m1+m2)
*
*************/
{
   double a, b, c, d, e, sum;
   a = mu2 / r;
   b = (1.0 - mu2) / sqrt( r*r - 2.0*r*l + 1.0 );
   c = 0.5 * (1.0 - n*n) * r*r;
   d = 0.5 * (1.0 - mu2) * (1.0 - mu2);
   e = -(1.0 - mu2) * r * l;
   sum = a + b + c + d + e;
   return( sum );
}

double findL1( double mu2 )
/******************
*
*  This function returns the distance of the inner Lagrangian
*  point from the center of star 2.  It uses a Newton's iteration
*  on d (psi) / dr with l=1 and n=0.  The input data are:
*     mu2 = the mass fraction = m2 / (m1+m2)
*
******************/
{
   double rzero, l, n, epsilonr, deriv, deltar, r;
   int i;

   rzero = mu2;  /* First guess at r. */
   l = 1.0;
   n = 0.0;
   /*  epsilonr = 1.0e-6 gives r to an accuracy of 1.0e-10 or better
          and will take 0 to 10 iterations over mu2 = ??????? */
   epsilonr = 1.0e-6;
   
   for( i=0; ; i++) {
      deriv = (dpsidr( (rzero+epsilonr), l, n, mu2) - dpsidr(rzero, l, n, mu2))
                         / epsilonr;
      deltar = -dpsidr(rzero, l, n, mu2) / deriv;
      r = rzero + deltar;
      if( fabs(deltar) < epsilonr ) break;
      if( i > 100 )
         quit("Too many iterations in function findL1!");
      rzero = r;
   }
   return( r );
}

double dpsidr( double r, double l, double n, double mu2)
/*****************
*
*  This function calculates d (psi) / dr.  The input data are:
*      r = distance from the center of star 2 to the point.
*      l = the direction cosine along the x direction
*      n = the direction cosine along the z direction
*      mu2 = the mass fraction = m2 / (m1+m2)
*
*******************/
{
   double a, b, c, d, x, sum;

   a = -mu2 / (r*r);
   x = fabs( r*r - 2.0*r*l + 1.0 );
   b = (1.0 - mu2) * (l - r) / ( sqrt(x*x*x) );
   c = (1.0 - n*n) * r;
   d = -(1.0 - mu2) * l;
   sum = a + b + c + d;
   return( sum );
}

double findr( double psiTarget, double rmax, double rguess, 
                double l, double n, double mu2)
/***************
*
*  This function finds the value of r that gives psi = psiTarget
*  using the Newton's iteration method.
*  THIS FUNCTION ONLY WORKS FOR psiTarget CORRESPONDING TO RADII
*  INSIDE OR ON THE ROCHE LOBE, NOT OUTSIDE THE ROCHE LOBE.
*  Input data:
*    psiTarget = the target value of psi
*    rmax      = maximum allowed value of r.  This is put in to
*                 avoid getting into the wrong Roche Lobe.
*    rguess    = the first guess at the value of r
*    l         = direction cosine in the x direction
*    n         = direction cosine in the z direction
*    mu2       = mass fraction = m2 / (m1+m2)
*
*****************/
{
   double deltar, r, epsilonr;
   int i;

   /*************
   *
   *    epsilonr = 1.0e-6 gives r to an accuracy better than 1.0e-10
   *    in less than 5 iterations.
   *
   **************/
   epsilonr = 1.0e-7;

   for( i=0; ; i++) {
      deltar = (psiTarget - psi(rguess, l, n, mu2)) / dpsidr(rguess, l, n, mu2);
      /*  Force r to stay in bounds */
      if( (rguess + deltar) > rmax ) deltar = 0.99 * rmax - rguess;
      if( (rguess + deltar) <= 0.0 ) deltar = -0.9 * rguess;
      r = rguess + deltar;
      if( fabs(deltar) < epsilonr ) break;
      if( i > 10 )
         quit("Too many iterations in function findr!");
      rguess = r;
   }   
   return( r );
}

/*******************
*
*  The following three functions, dpsidx, dpsidy, and dpsidz
*  are fairly straight forward.
*
********************/
double dpsidx( double r, double l, double m, double n, double mu2)
{
   double x, a, b, c, d, deriv;

   x = sqrt( r*r - 2.0*r*l + 1.0 );
   a = -(1.0 - mu2) * (r*l - 1.0) / (x*x*x);
   b = -mu2 * l / (r*r);
   c = r*l;
   d = -(1.0 - mu2);   
   deriv = a + b + c + d;
   return( deriv );
}
double dpsidy( double r, double l, double m, double n, double mu2)
{
   double x, a, b, c, deriv;
   
   x = sqrt( r*r - 2.0*r*l + 1.0 );
   a = -(1.0 - mu2) * r * m / (x*x*x);
   b = -mu2 * m / (r*r);
   c = r * m;
   deriv = a + b + c;
   return( deriv );
}
double dpsidz( double r, double l, double m, double n, double mu2)
{
   double x, a, b, deriv;

   x = sqrt( r*r - 2.0*r*l + 1.0 );
   a = -(1.0 - mu2) * r * n / (x*x*x);
   b = -mu2 * n / (r*r);
   deriv = a + b;
   return( deriv );
}


void DiskGeom( void )
/*******************************
*
*  This is the main function that sets up the geometry of the disk.
*
*********************************/
{
   double twopi, alpha, hangle, deltarho, deltaphi, phi, rho, deltaz,
      L1distance, maxstar1rad, massfraction;
   long int i, j, maxrindex, maxphiindex, ir, iphi;
   
   twopi = 2.0 * 3.14159265358979323846;

   /*******************************
   *
   *  The user inputs the inner and outer radii of the disk in units
   *  of RL1.  Convert the radii to units of the orbital separation
   *  and then calculate a few other useful numbers.
   *
   ********************************/
   massfraction = starsys.q / ( 1.0 + starsys.q );
   L1distance = findL1( massfraction );
   disktop.innerR   = disktop.innerR   * L1distance;
   disktop.outerR   = disktop.outerR   * L1distance;
   disktop.spotr[0] = disktop.spotr[0] * L1distance;
   disktop.spotr[1] = disktop.spotr[1] * L1distance;
   diskrim.outerR   = diskrim.outerR   * L1distance;
   maxstar1rad = 0.0;
   if( strcmp( mode.star1, "ON") == 0) {
      for( i = 0; i < star1.thetapoints; i++) {
         for( j = 0; j < star1.phipoints; j++) {
            if( r1[i][j] > maxstar1rad ) maxstar1rad = r1[i][j];
         }
      }
      if( disktop.innerR < maxstar1rad ) {
         printf("   RDISKIN is less than the radius of star 1.\n");
         printf("   It has been increased to %8.5lf\n", maxstar1rad);
         disktop.innerR = maxstar1rad;
      }
   }

   /********************************
   *
   *  Calculate the values of r and the direction cosines
   *  at each of the grid points.  Note that 3 points are put on
   *  the edge of the disk and the rest on the top of the disk.
   *  There are no grid points on the bottom.  Also note that the
   *  highest of the edge grid points is at the same place as the
   *  top grid point with the largest r.  They are distinguished
   *  by different surface normals.
   *
   *  Note that rho is the radius in the plane of the orbit, while
   *  rdisk[][] is the distance from the center of the disk to the
   *  point on the surface of the disk.
   *
   *********************************/
   maxrindex = disktop.rpoints - 1;
   maxphiindex = disktop.phipoints - 1;
   alpha = twopi * (disktop.flareangle / 360.0);
   deltarho = (disktop.outerR - disktop.innerR) / (disktop.rpoints - 4);
   deltaphi = twopi / (maxphiindex + 1.0);
   deltaz = disktop.outerR * tan( alpha );

   for( iphi = 0; iphi <= maxphiindex; iphi++) {
      phi = iphi * deltaphi;
      for( ir = 0; ir <= maxrindex; ir++) {
         if( ir < (maxrindex - 2) ) {
            rho = ir * deltarho + disktop.innerR;
            hangle = alpha;
         }
         else {
            rho = disktop.outerR;
            hangle = alpha * ( maxrindex - ir - 1 );
         }
         rdisk[ir][iphi] = rho / cos( hangle );
         ldisk[ir][iphi] = cos( hangle ) * cos( phi );
         mdisk[ir][iphi] = cos( hangle ) * sin( phi );
         ndisk[ir][iphi] = sin( hangle );
      }
   }
   /***********************
   *
   *  Now calculate the normals to the surface elements on the disk.
   *
   ***********************/
   for( iphi = 0; iphi <= maxphiindex; iphi++) {
      phi = iphi * deltaphi;
      for( ir = 0; ir <= maxrindex; ir++) {
         if( ir < (maxrindex - 2) ) {
            normdiskx[ir][iphi] = -sin( alpha ) * cos( phi );
            normdisky[ir][iphi] = -sin( alpha ) * sin( phi );
            normdiskz[ir][iphi] =  cos( alpha );
         }
         else {
            normdiskx[ir][iphi] = cos( phi );
            normdisky[ir][iphi] = sin( phi );
            normdiskz[ir][iphi] = 0.0;
         }
      }
   }   
   /***************************
   *
   *  Now calculate the area elements at each grid point.  Note
   *  that the second order corrections at the edges of the grid
   *  are important.
   *
   ****************************/
   for( iphi = 0; iphi <= maxphiindex; iphi++) {
      phi = iphi * deltaphi;
      for( ir = 0; ir <= maxrindex; ir++) {
         if( ir == 0 ) {
            rho = disktop.innerR;
            dSdisk[ir][iphi] =   0.5 * rho * deltarho * deltaphi
                               + 0.125 * deltarho * deltarho * deltaphi;
            dSdisk[ir][iphi] = dSdisk[ir][iphi] / cos( alpha );
         }
         else if( (ir > 0) && (ir < (maxrindex - 3)) ) {
            rho = ir * deltarho + disktop.innerR;
            dSdisk[ir][iphi] = rho * deltarho * deltaphi;
            dSdisk[ir][iphi] = dSdisk[ir][iphi] / cos( alpha );
         }
         else if( ir == (maxrindex - 3) ) {
            rho = disktop.outerR;
            dSdisk[ir][iphi] =   0.5 * rho * deltarho * deltaphi
                               - 0.125 * deltarho * deltarho * deltaphi;
            dSdisk[ir][iphi] = dSdisk[ir][iphi] / cos( alpha );
         }
         else if( ir == (maxrindex - 2) ) {
            dSdisk[ir][iphi] = 0.5 * disktop.outerR * deltaphi * deltaz;
         }
         else if( ir == (maxrindex - 1) ) {
            dSdisk[ir][iphi] =       disktop.outerR * deltaphi * deltaz;
         }
         else if( ir == (maxrindex - 0) ) {
            dSdisk[ir][iphi] = 0.5 * disktop.outerR * deltaphi * deltaz;
         }
         else {
            quit("Error in DiskGeom.");
         }
      }
   }
   return;
}