/********************************
*
*                    FILE   TEFF.C
*
*  This file contains the functions necessary to calculate the
*  local T_eff at each point on the surface of the objects in
*  the binary system.
*
************************************/

#include "header.h"

void Tstars( char *object )
/***************
*
*  This function calculates the local effective temperature at each
*  point on the surfaces of the two stars from the effective
*  temperature at the pole and the gravity darkening.
*
****************/
{
   long int itheta, iphi;
   long int maxthetaindex1, maxphiindex1, maxthetaindex2, maxphiindex2;
   double gravratio;

   if( strcmp( object, "star1") == 0 ) {
      maxthetaindex1 = star1.thetapoints - 1;
      maxphiindex1 = star1.phipoints - 1;
      for( itheta=0; itheta <= maxthetaindex1; itheta++) {
         for( iphi=0; iphi <= maxphiindex1; iphi++) {
            gravratio = gravity1[itheta][iphi] / star1.polegrav;
            localT1[itheta][iphi] = star1.poletemp *
                                        pow( gravratio, star1.gravbeta);
         }
      }
   }

   if( strcmp( object, "star2") == 0 ) {
      maxthetaindex2 = star2.thetapoints - 1;
      maxphiindex2 = star2.phipoints - 1;
      for( itheta=0; itheta <= maxthetaindex2; itheta++) {
         for( iphi=0; iphi <= maxphiindex2; iphi++) {
            gravratio = gravity2[itheta][iphi] / star2.polegrav;
            localT2[itheta][iphi] = star2.poletemp *
                                        pow( gravratio, star2.gravbeta);
         }
      }
   }
   return;
}


void Tstar1spot( double theta, double phi, double radius, 
                 double SpotToverStarT )
/******************************************
*
*  This function puts a spot on star 1.  Note that if two
*  spots overlap, the temperature reduction in the overlap
*  region is the product of the individual reductions.
*
********************************************/
{
   long itheta, iphi, maxthetaindex1, maxphiindex1;
   double lspot, mspot, nspot, cosRadius, cosAngle;

   lspot = sin(theta) * cos(phi);
   mspot = sin(theta) * sin(phi);
   nspot = cos(theta);
   cosRadius = cos(radius);

   maxthetaindex1 = star1.thetapoints - 1;
   maxphiindex1 = star1.phipoints - 1;

   for( itheta=0; itheta <= maxthetaindex1; itheta++) {
      for( iphi=0; iphi <= maxphiindex1; iphi++) {
         cosAngle =   lspot * l1[itheta][iphi]
                    + mspot * m1[itheta][iphi]
                    + nspot * n1[itheta][iphi];
         if( cosAngle >= cosRadius ) {
            localT1[itheta][iphi] = SpotToverStarT 
                                       * localT1[itheta][iphi];
         }
      }
   }

   return;
}


void Tstar2spot( double theta, double phi, double radius, 
                 double SpotToverStarT )
/******************************************
*
*  This function puts a spot on star 2.  Note that if two
*  spots overlap, the temperature reduction in the overlap
*  region is the product of the individual reductions.
*
********************************************/
{
   long itheta, iphi, maxthetaindex2, maxphiindex2;
   double lspot, mspot, nspot, cosRadius, cosAngle;

   lspot = sin(theta) * cos(phi);
   mspot = sin(theta) * sin(phi);
   nspot = cos(theta);
   cosRadius = cos(radius);

   maxthetaindex2 = star2.thetapoints - 1;
   maxphiindex2 = star2.phipoints - 1;

   for( itheta=0; itheta <= maxthetaindex2; itheta++) {
      for( iphi=0; iphi <= maxphiindex2; iphi++) {
         cosAngle =   lspot * l2[itheta][iphi]
                    + mspot * m2[itheta][iphi]
                    + nspot * n2[itheta][iphi];
         if( cosAngle >= cosRadius ) {
            localT2[itheta][iphi] = SpotToverStarT 
                                       * localT2[itheta][iphi];
         }
      }
   }

   return;
}


void Tdisk( void )
/******************************
*
*  This function sets up the distribution of T across the
*  accretion disk and its rim.
*
*  NOTE:  The values of r in the array rdisk[][] are the
*  distances from the center of the disk to a point on the
*  surface of the disk.  The radius of the disk as measured
*  in the plane of the disk is  rdisk[][] * cos( flareangle );
*
***************************************/
{
   double angle;
   double r, deltar, r1spot, r2spot, phi, deltaphi, phi1, phi2;
   long int ir, minrindex, maxrindex, iphi, maxphiindex;

   angle = ( disktop.flareangle / 360.0 ) * (2.0 * 3.14159265358979);
   maxphiindex = disktop.phipoints - 1;

   /***************************
   *
   *  First put the temperature distribution on the top of the disk
   *
   *****************************/
   if( strcmp( disktop.Tsource, "POLYNOMIAL") == 0 ) {
      maxrindex = disktop.rpoints - 4;
      for( ir = 0; ir <= maxrindex; ir++) {
      	 r = rdisk[ir][0] * cos( angle );
      	 localTdisk[ir][0] =   disktop.Tcoef[0]
                             + disktop.Tcoef[1] * r
                             + disktop.Tcoef[2] * r * r
                             + disktop.Tcoef[3] * r * r * r;
         if( localTdisk[ir][0] < 0.0 ) {
            printf("   The polynomial coefficients for the disk\n");
            printf("   temperature make the temperature go negative.");
            quit("");
         }
         for( iphi = 0; iphi <= maxphiindex; iphi++) {
            localTdisk[ir][iphi] = localTdisk[ir][0];
         }
      }
   }
   else if( strcmp( disktop.Tsource, "STEADYSTATE") == 0 ) {
      quit("TEMPDISK = STEADYSTATE not yet implimented.");
   }
   else if( strcmp( disktop.Tsource, "FILE") == 0 ) {
      readDiskT( disktop.Tfile );
   }
   else
      quit("Unrecognized source for TEMPDISK.");

   /***********************
   *
   *  Put a spot on top of the disk if desired.  The spot is inserted 
   *      if its temperature is greater than or equal to zero.
   *  The two radii are the radii in the plane of the orbit, not along
   *     the surface of the disk.
   *  The spot is fit onto the grid points by rounding rin, rout,
   *     phi1 and phi2 to values at grid lines.
   *
   ************************/
   if( disktop.spotT >= 0.0 ) {
      deltaphi = 360.0 / disktop.phipoints;
      phi1 = disktop.spotphi[0] - 0.5 * deltaphi;
      phi2 = disktop.spotphi[1] + 0.5 * deltaphi;
      deltar = ( rdisk[1][0] - rdisk[0][0] ) / cos( angle );
      r1spot = disktop.spotr[0] / cos( angle ) - 0.5*deltar;
      r2spot = disktop.spotr[1] / cos( angle ) + 0.5*deltar;
      maxrindex = disktop.rpoints - 4;
      for( ir = 0; ir <= maxrindex; ir++) {
         if( (rdisk[ir][0] >= r1spot) && (rdisk[ir][0] <= r2spot) ) {
            for( iphi = 0; iphi <= maxphiindex; iphi++) {
               phi = iphi * deltaphi;
               if( (phi2 >= phi1) && (phi >= phi1) && (phi <= phi2) )
                  localTdisk[ir][iphi] = disktop.spotT;
               if( (phi2 <= phi1) && ( (phi >= phi1) || (phi <= phi2) ) )
                  localTdisk[ir][iphi] = disktop.spotT;
            }
         }
      }
   }

   /*************************
   *
   *  Give a temperature to the rim of the disk.
   *
   **************************/
   minrindex = disktop.rpoints - 3;
   maxrindex = disktop.rpoints - 1;
   for( ir = minrindex; ir <= maxrindex; ir++) {
      for( iphi = 0; iphi <= maxphiindex; iphi++) {
         localTdisk[ir][iphi] = diskrim.T;
      }
   }

   /***********************
   *
   *  Now put a spot on the rim of the disk if desired.  The spot is 
   *  inserted if its temperature is greater than or equal to zero.
   *  The spot is fit onto the grid points by rounding phi1 and phi2 to 
   *     values at grid lines.
   *
   ************************/
   if( diskrim.spotT >= 0.0 ) {
      minrindex = disktop.rpoints - 3;
      maxrindex = disktop.rpoints - 1;
      deltaphi = 360.0 / disktop.phipoints;
      phi1 = diskrim.spotphi[0] - 0.5 * deltaphi;
      phi2 = diskrim.spotphi[1] + 0.5 * deltaphi;
      for( ir = minrindex; ir <= maxrindex; ir++) {
         for( iphi = 0; iphi <= maxphiindex; iphi++) {
            phi = iphi * deltaphi;
            if( (phi2 >= phi1) && (phi >= phi1) && (phi <= phi2) )
               localTdisk[ir][iphi] = diskrim.spotT;
            if( (phi2 <= phi1) && ( (phi >= phi1) || (phi <= phi2) ) )
               localTdisk[ir][iphi] = diskrim.spotT;
         }
      }
   }

   return;
}


void irradiate( void )
/****************************************************
*
*  This function calculates the so-called reflection effect
*  between the two stars.  The effective temperature of star 1, 
*  T1_eff, due to the reflection effect is calculated from
*
*     (T1_eff)^4  = (1/sigma) * F1_out
*
*  where,
*
*     F1_out = F1_intrinsic + albedo*F_irr
*
*  and
*
*     d( F_irr ) = sigma*(T2_eff)^4 * limbdarkening 
*                             * mu1 * mu2 * d( S2 ) / distance^2
*
*  Note that the limb darkening is for the total flux, not for
*  monochromatic flux.  I have chosen to use
*
*     limbdarkening = (1/pi) * (1 - u2 + u2*mu2)
*
*  where the factor of 1/pi is a normalization constant chosen so
*
*     integral( limbdarkening * mu * d(omega)) = 1
*
*  and the limb darkening coefficient is taken to be appropriate for
*  the sun:
*
*     u = ~0.73   (see Allen, AQ, 3rd edition, p.170)
*
*  One might also consider using a grey atmosphere limb darkening.
*
***************************************************/               
{
   double localB1 [SIZEDIMEN1] [SIZEDIMEN2],
          localB2 [SIZEDIMEN1] [SIZEDIMEN2],
            Firr1 [SIZEDIMEN1] [SIZEDIMEN2],
            Firr2 [SIZEDIMEN1] [SIZEDIMEN2],
          x1prime [SIZEDIMEN1] [SIZEDIMEN2],
          y1prime [SIZEDIMEN1] [SIZEDIMEN2],
          z1prime [SIZEDIMEN1] [SIZEDIMEN2],
               x2 [SIZEDIMEN1] [SIZEDIMEN2],
               y2 [SIZEDIMEN1] [SIZEDIMEN2],
               z2 [SIZEDIMEN1] [SIZEDIMEN2];
   double sigma, u1, u2, oneoverpi, dx, dy, dz, oneoverd, oneoverdsqr;
   double normalx1prime, normaly1prime, normalz1prime, project1, project2;
   double mu1, mu2, t4, sumflux;
   long int itheta1, iphi1, itheta2, iphi2;
   long int maxthetaindex1, maxphiindex1, maxthetaindex2, maxphiindex2;

   if( (star1.albedo == 0.0) && (star2.albedo == 0) )
                 return;

   sigma = 5.66956e-5;   /*  in cgs units   */
   u1 = 0.73;
   u2 = 0.73;
   oneoverpi = 1.0 / 3.14159265358979;
   
   maxthetaindex1 = star1.thetapoints - 1;
   maxphiindex1 = star1.phipoints - 1;
   maxthetaindex2 = star2.thetapoints - 1;
   maxphiindex2 = star2.phipoints - 1;

   for( itheta1 = 0; itheta1 <= maxthetaindex1; itheta1++) {
      for( iphi1 = 0; iphi1 <= maxphiindex1; iphi1++) {
         localB1[itheta1][iphi1] = sigma * pow( localT1[itheta1][iphi1], 4.0);
         Firr1[itheta1][iphi1] = 0.0;
         x1prime[itheta1][iphi1] = 1.0 -r1[itheta1][iphi1] * l1[itheta1][iphi1];
         y1prime[itheta1][iphi1] =     -r1[itheta1][iphi1] * m1[itheta1][iphi1];
         z1prime[itheta1][iphi1] =      r1[itheta1][iphi1] * n1[itheta1][iphi1];
      }
   }
   for( itheta2 = 0; itheta2 <= maxthetaindex2; itheta2++) {
      for( iphi2 = 0; iphi2 <= maxphiindex2; iphi2++) {
         localB2[itheta2][iphi2] = sigma * pow( localT2[itheta2][iphi2], 4.0);
         Firr2[itheta2][iphi2] = 0.0;
         x2[itheta2][iphi2] = r2[itheta2][iphi2] * l2[itheta2][iphi2];
         y2[itheta2][iphi2] = r2[itheta2][iphi2] * m2[itheta2][iphi2];
         z2[itheta2][iphi2] = r2[itheta2][iphi2] * n2[itheta2][iphi2];
      }
   }

   /*********************
   *
   *   The following 4-deep set of loops is where the irradiation
   *   on each surface element is added.
   *
   ***********************/   
   for( itheta1 = 0; itheta1 <= maxthetaindex1; itheta1++) {
      for( iphi1 = 0; iphi1 <= maxphiindex1; iphi1++) {
         for( itheta2 = 0; itheta2 <= maxthetaindex2; itheta2++) {
            for( iphi2 = 0; iphi2 <= maxphiindex2; iphi2++) {
               dx = x1prime[itheta1][iphi1] - x2[itheta2][iphi2];
               dy = y1prime[itheta1][iphi1] - y2[itheta2][iphi2];
               dz = z1prime[itheta1][iphi1] - z2[itheta2][iphi2];
               normalx1prime = -normalx1[itheta1][iphi1];
               normaly1prime = -normaly1[itheta1][iphi1];
               normalz1prime =  normalz1[itheta1][iphi1];
               project2 =   dx * normalx2[itheta2][iphi2]
                          + dy * normaly2[itheta2][iphi2]
                          + dz * normalz2[itheta2][iphi2];
               project1 =   dx * normalx1prime
                          + dy * normaly1prime
                          + dz * normalz1prime;
               if( (project2 > 0.0) && (project1 < 0.0) ) {
                  oneoverdsqr = 1.0 / ( dx*dx + dy*dy + dz*dz);
                  oneoverd = sqrt( oneoverdsqr );
                  mu1 = -project1 * oneoverd;
                  mu2 =  project2 * oneoverd;
                  Firr1[itheta1][iphi1] =  Firr1[itheta1][iphi1]
                                             + localB2[itheta2][iphi2]
                                             * oneoverpi*(1.0-u2+u2*mu2) 
                                             * mu1 * mu2 * oneoverdsqr
                                             * dS2[itheta2][iphi2];
                  Firr2[itheta2][iphi2] = Firr2[itheta2][iphi2]
                                             + localB1[itheta1][iphi1]
                                             * oneoverpi*(1.0-u1+u1*mu1) 
                                             * mu1 * mu2 * oneoverdsqr
                                             * dS1[itheta1][iphi1];
               }
            }
         }
      }
   }

   for( itheta1 = 0; itheta1 <= maxthetaindex1; itheta1++) {
      for( iphi1 = 0; iphi1 <= maxphiindex1; iphi1++) {
         sumflux = localB1[itheta1][iphi1]
                     + star1.albedo * Firr1[itheta1][iphi1];
         t4 = sumflux / sigma;
         localT1[itheta1][iphi1] = pow( t4, 0.25);
      }
   }
   for( itheta2 = 0; itheta2 <= maxthetaindex2; itheta2++) {
      for( iphi2 = 0; iphi2 <= maxphiindex2; iphi2++) {
         sumflux = localB2[itheta2][iphi2]
                  + star2.albedo * Firr2[itheta2][iphi2];
         t4 = sumflux / sigma;
         localT2[itheta2][iphi2] = pow( t4, 0.25);
      }
   }

   return;
}