/* NIGHTFALL Light Curve Synthesis Program */ /* Copyright (C) 1998 Rainer Wichmann */ /* */ /* This program is free software; you can redistribute it */ /* and/or modify */ /* it under the terms of the GNU General Public License as */ /* published by */ /* the Free Software Foundation; either version 2 of the License, or */ /* (at your option) any later version. */ /* */ /* This program is distributed in the hope that it will be useful, */ /* but WITHOUT ANY WARRANTY; without even the implied warranty of */ /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */ /* GNU General Public License for more details. */ /* */ /* You should have received a copy of the GNU General Public License */ /* along with this program; if not, write to the Free Software */ /* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */ #include #include #include #include #include "Light.h" #ifdef _WITH_PGPLOT #ifndef _WITH_GNUPLOT #include "cpgplot.h" #endif #endif /**************************************************************************** @package nightfall @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de) @version 1.0 @short Compute fractional visibility @param (int) Comp The stellar component @return (void) @heading Star Setup ****************************************************************************/ void LightFractionalPot(int Comp) { long NElem; /* # of surface elements */ long i, j; /* loop variables */ static long LastEclipsed[MAXELEMENTS]; /* array of indices to eclipsed */ static long NextEclipsed[MAXELEMENTS]; /* array of indices to eclipsed */ long Last; /* element eclipsed before */ long Next; /* element eclipsed after */ double lzero, mzero, nzero; /* surface normal vector */ int Comp2; /* the other star */ double LOS_PotLast; /* line-of-sight potential min */ double LOS_PotNext; /* line-of-sight potential min */ double LOS_PotInt; /* line-of-sight potential min */ double x_z, y_z, z_z; /* surface vector */ double t_1, t_2; /* brackets for minfind */ double tmin = 0; /* minimum along LOS */ double Phase; /* orbital phase */ int Acomp; /* sign */ double Qi, Li, QLi; /* smalles circle around star */ double t_l; /* temporary storage */ double SqrXp; /* misc computational */ double temp; /* temporary storage */ double fratio = 1.0; /* the async rotation */ j = 0; t_1 = 0.03; t_2 = 1.0; /* secondary is mirrored - we turn by PI and mirror y->(-y) */ /* to have secondary at (0,0,0) with correct orientation */ /* we are interested in the potential surface of the OTHER star */ if (Comp == Primary) { Phase = Orbit.Phase; Comp2 = 1; Acomp = 1; SqrXp = SQR(Binary[Secondary].Xp); /* SqrPol = SQR(Binary[Secondary].Radius); */ if (Flags.asynchron2 == ON) fratio = Binary[Secondary].Fratio; } else { Phase = Orbit.Phase + PI; Comp2 = 0; Acomp = (-1); SqrXp = SQR(Binary[Primary].Xp); /* SqrPol = SQR(Binary[Primary].Radius); */ if (Flags.asynchron1 == ON) fratio = Binary[Primary].Fratio; } /* ------- vector towards observer (LOS = Line Of Sight) --------- */ lzero = sin(Orbit.Inclination) * cos(Phase); mzero = sin(Orbit.Inclination) * sin(Phase); nzero = cos(Orbit.Inclination); /* ----------- search for eclipsed/non-eclipsed pairs ----------- */ NElem = Binary[Comp].NumElem; for (i = 0; i < NElem; ++i) { LastEclipsed[i] = 0; NextEclipsed[i] = 0; } j = 0; for (i = Binary[Comp].N_PhiStep[0]; i < NElem; ++i) { if (fabs((double)Surface[Comp][i].EclipseFlag - (double)Surface[Comp][i].next_ptr->EclipseFlag - 1.0) <= FLT_EPSILON) { NextEclipsed[j] = Surface[Comp][i].next_ptr->MyNum; if (NextEclipsed[j] > (Binary[Comp].NumElem - 1)) NextEclipsed[j] = Binary[Comp].NumElem - 1; LastEclipsed[j] = i; ++j; } if ( fabs((double)Surface[Comp][i].EclipseFlag - (double)Surface[Comp][i].last_ptr->EclipseFlag - 1.0) <= FLT_EPSILON) { NextEclipsed[j] = Surface[Comp][i].last_ptr->MyNum; if (NextEclipsed[j] > (Binary[Comp].NumElem - 1)) NextEclipsed[j] = Binary[Comp].NumElem - 1; LastEclipsed[j] = i; ++j; } if (fabs((double)Surface[Comp][i].EclipseFlag - (double) Surface[Comp][i].phi_ptr->EclipseFlag - 1.0) <= FLT_EPSILON) { NextEclipsed[j] = Surface[Comp][i].phi_ptr->MyNum; if (NextEclipsed[j] > (Binary[Comp].NumElem - 1)) NextEclipsed[j] = Binary[Comp].NumElem - 1; LastEclipsed[j] = i; ++j; } if (fabs((double)Surface[Comp][i].EclipseFlag - (double)Surface[Comp][i].phi1_ptr->EclipseFlag - 1.0) <= FLT_EPSILON) { NextEclipsed[j] = Surface[Comp][i].phi1_ptr->MyNum; if (NextEclipsed[j] > (Binary[Comp].NumElem - 1)) NextEclipsed[j] = Binary[Comp].NumElem - 1; LastEclipsed[j] = i; ++j; } } /* ---------- compute fractional visibility ------------------- */ for (i = 0; i < j; ++i) { Last = LastEclipsed[i]; Next = NextEclipsed[i]; if (fabs(Surface[Comp][Last].LOS_Pot) >= FLT_EPSILON) { LOS_PotLast = Surface[Comp][Last].LOS_Pot; } else { x_z = Surface[Comp][Last].lambda; y_z = Acomp * (Surface[Comp][Last].mu); z_z = Surface[Comp][Last].nu; Qi = 2.0*(x_z*lzero+y_z*mzero+z_z*nzero-lzero); Li = 1.0 + (x_z*x_z)+(y_z*y_z)+(z_z*z_z)-2*x_z-SqrXp; QLi = (Qi*Qi)-4.0*Li; if (QLi >= FLT_EPSILON) { /* intersection */ t_l = sqrt(QLi); /* just to hold temporary result */ t_1 = (-Qi + t_l)/2.0; t_2 = (-Qi - t_l)/2.0; } else { t_1 = 0.2; t_2 = 1.0; } /* changed Binary[Comp].Mq to Binary[Comp2].Mq Jul 6 */ (void) MinFinder(&Binary[Comp2].Mq, &fratio, &t_1, &t_2, &lzero, &mzero, &nzero, &x_z, &y_z, &z_z, &tmin, &LOS_PotLast); LOS_PotLast = -LOS_PotLast; Surface[Comp][Last].LOS_Pot = (float) LOS_PotLast; } if (fabs(Surface[Comp][Next].LOS_Pot) >= FLT_EPSILON) { LOS_PotNext = Surface[Comp][Next].LOS_Pot; } else { x_z = Surface[Comp][Next].lambda; y_z = Acomp * (Surface[Comp][Next].mu); z_z = Surface[Comp][Next].nu; Qi = 2.0*(x_z*lzero+y_z*mzero+z_z*nzero-lzero); Li = 1.0 + (x_z*x_z)+(y_z*y_z)+(z_z*z_z)-2*x_z-SqrXp; QLi = Qi*Qi-4.0*Li; if (QLi >= FLT_EPSILON) { /* intersection */ t_l = sqrt(QLi); /* just to hold temporary result */ t_1 = (-Qi + t_l)/2.0; t_2 = (-Qi - t_l)/2.0; } else { t_1 = 0.2; t_2 = 1.0; } /* changed Binary[Comp].Mq to Binary[Comp2].Mq Jul 6 */ (void) MinFinder(&Binary[Comp2].Mq, &fratio, &t_1, &t_2, &lzero, &mzero, &nzero, &x_z, &y_z, &z_z, &tmin, &LOS_PotNext); LOS_PotNext = -LOS_PotNext; Surface[Comp][Next].LOS_Pot = (float) LOS_PotNext; } LOS_PotInt = fabs((LOS_PotNext - LOS_PotLast)/2.); /* EclipseFlag == 1 means IsEclipsed */ temp = fabs(LOS_PotNext - Binary[Comp2].RochePot); if ((temp+FLT_EPSILON) <= fabs(LOS_PotInt)) { if (Surface[Comp][Next].EclipseFlag == 1) Surface[Comp][Next].visibility = (float) (0.0 + 0.5 * (1.0 - temp/ LOS_PotInt )); if (Surface[Comp][Next].EclipseFlag == 0) Surface[Comp][Next].visibility = (float) (1.0 - 0.5 * (1.0 - temp/ LOS_PotInt )); } else { Surface[Comp][Next].visibility = (float) fabs(Surface[Comp][Next].EclipseFlag-1.0); } temp = fabs(LOS_PotLast - Binary[Comp2].RochePot); if ((temp+FLT_EPSILON) <= fabs(LOS_PotInt)) { if (Surface[Comp][Last].EclipseFlag == 1) Surface[Comp][Last].visibility = (float) (0.0 + 0.5 * (1.0 - temp/ LOS_PotInt )); if (Surface[Comp][Last].EclipseFlag == 0) Surface[Comp][Last].visibility = (float) (1.0 - 0.5 * (1.0 - temp/ LOS_PotInt )); } else { Surface[Comp][Last].visibility = (float) fabs(Surface[Comp][Last].EclipseFlag-1.0); } } }