/* 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" /**************************************************************************** @package nightfall @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de) @version 1.0 @short Attach the throat of star a to star b @tip This is to verify eclipse by 'own' star @param (void) @return (void) @heading Light Curve ****************************************************************************/ void LightCopyThroat() { long j; /* loop counter */ long Nplus = 0; /* surface elements + throat */ long NP = 0; /* surface elements primary */ long NS = 0; /* surface elements secondary */ /* -------- attach Secondary throat to Primary ---------------------- */ Nplus = Binary[Secondary].N_PhiStep[StepsPerPi-1] + Binary[Secondary].N_PhiStep[StepsPerPi-2]; NP = Binary[Primary].NumElem; NS = Binary[Secondary].NumElem; Binary[Primary].NumPlus = NP + Nplus; if (Binary[Primary].NumPlus > (long) MaximumElements) nf_error(_(errmsg[10])); for (j = 0; j < Nplus; ++j) { Surface[Primary][NP+j] = Surface[Secondary][NS-Nplus+j]; /* --------------- INVERT X --------------------------------------- */ Surface[Primary][NP+j].lambda = (float) (1.0 - Surface[Primary][NP+j].lambda); Surface[Primary][NP+j].l = (-Surface[Primary][NP+j].l); } /* ---------------- attach Primary throat to Secondary ---------------- */ Nplus = Binary[Primary].N_PhiStep[StepsPerPi-1] + Binary[Primary].N_PhiStep[StepsPerPi-2]; Binary[Secondary].NumPlus = NS + Nplus; if (Binary[Primary].NumPlus > (long) MaximumElements) nf_error(_(errmsg[10])); for (j = 0; j < Nplus; ++j) { Surface[Secondary][NS+j] = Surface[Primary][NP-Nplus+j]; /* --------------------- INVERT X --------------------------------- */ Surface[Secondary][NS+j].lambda = (float) (1.0 - Surface[Secondary][NS+j].lambda); Surface[Secondary][NS+j].l = (-Surface[Secondary][NS+j].l); } return; } /**************************************************************************** @package nightfall @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de) @version 1.0 @short Copy the throat back @param (void) @return (void) @heading Light Curve ****************************************************************************/ void LightCopyThroatBack() { long j; /* loop counter */ long Nplus = 0; /* surface elements + throat */ long NP = 0; /* surface elements primary */ long NS = 0; /* surface elements secondary */ NP = Binary[Primary].NumElem; NS = Binary[Secondary].NumElem; /* ----------- copy back Primary throat to Primary ------------------- */ Nplus = Binary[Secondary].NumPlus - NS; for (j = 0; j < Nplus; ++j) { if(Surface[Secondary][NS+j].EclipseFlag == ON){ Surface[Primary][NP-Nplus+j].EclipseFlag = Surface[Secondary][NS+j].EclipseFlag; Surface[Primary][NP-Nplus+j].visibility = Surface[Secondary][NS+j].visibility; } } /* ------------- copy back Secondary throat to Secondary ------------- */ Nplus = Binary[Primary].NumPlus - NP; for (j = 0; j < Nplus; ++j) { if (Surface[Primary][NP+j].EclipseFlag == ON){ Surface[Secondary][NS-Nplus+j].EclipseFlag = Surface[Primary][NP+j].EclipseFlag; Surface[Secondary][NS-Nplus+j].visibility = Surface[Primary][NP+j].visibility; } } return; } /******************************************************************** @package nightfall @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de) @version 1.0 @short Calculate various quantities for simple geometric tests @param (void) @return (int) Exit status @heading Light Curve ********************************************************************/ int LightSetupTests() { int testerr = 0; /* Exit status RootFind */ double RLag1; /* Lagrange One */ /* double Xl, h1, h2; */ /* Xl lagrange one */ double Xp1, Xp2; /* max. distance surface-center of star */ /* double PhiLagrange; */ /* tangent angle of roche lobes at L1 */ double CosPhiLagrange; /* tangent angle at Lagrange One (cos) */ double CosPhiOne = 0.0; /* perhaps more stringent if Roche lobe */ /* not filled */ /* ---------- simple geometric tests for eclipse possible ------ */ /* osculating cone angle at L1 */ /* Xl = Binary[Primary].RLag1; h1 = 1.0/(Xl*Xl*Xl); h2 = 1.0/((1.0-Xl)*(1.0-Xl)*(1.0-Xl)); PhiLagrange = sqrt((2*h1 + 2*Mq*h2 + (1.0 - Mq))/(h1 + Mq*h2 - (1.0 - Mq))); PhiLagrange = atan(PhiLagrange); CosPhiLagrange = cos(PhiLagrange); */ Xp1 = 1.0; Xp2 = 1.0; /* calculate Xp == distance from center to surface towards L1 */ /* lambda, nu global; potential on x-axis */ lambda = 1.0; nu = 0.0; RochePot = Binary[Primary].RochePot; Mq = Binary[Primary].Mq; if (Flags.asynchron1 == ON) F = Binary[Primary].Fratio; else F = 1.0; RLag1 = Binary[Primary].RLag1; if (fabs(Binary[Primary].RocheFill - 1.0) <= FLT_EPSILON) { Xp1 = Binary[Primary].RXCrit; } else { if (Flags.fill == OFF) { Xp1 = RootFind(RocheSurface, 0.00001, Binary[Primary].RXCrit, 1.0e-8, "LightSetupTests", &testerr); if (testerr == 1) return(8); } else { Xp1 = Binary[Primary].LimRad; } } Binary[Primary].Xp = Xp1; RochePot = Binary[Secondary].RochePot; Mq = Binary[Secondary].Mq; if (Flags.asynchron2 == ON) F = Binary[Secondary].Fratio; else F = 1.0; RLag1 = Binary[Secondary].RLag1; if (fabs(Binary[Secondary].RocheFill - 1.0) <= FLT_EPSILON) { Xp2 = Binary[Secondary].RXCrit; } else { if (Flags.fill == OFF) { Xp2 = RootFind(RocheSurface, 0.00001,Binary[Secondary].RXCrit, 1.0e-8, "LightSetupTests", &testerr ); if (testerr == 8) return(8); } else { Xp2 = Binary[Secondary].LimRad; } } Binary[Secondary].Xp = Xp2; if (Flags.debug[GEO] == ON) { printf("\nXp[P]: %6.2f Xp[S]: %6.2f\n", Binary[Primary].Xp, Binary[Secondary].Xp); } /* ------ sanity check if F < 1.0 ------------------------------ */ if ( ( (Binary[Secondary].Fratio - 1.0) <= (-FLT_EPSILON) && Flags.asynchron1 == ON ) || ( (Binary[Primary].Fratio - 1.0) <= (-FLT_EPSILON) && Flags.asynchron2 == ON ) ) { if ((Binary[Secondary].Xp + Binary[Primary].Xp) >= (1.0+DBL_EPSILON) ) { WARNING(_("Fill Factor too large in Async Rotating Star --> Stars intersect")); return(1); } } if (Flags.fill == OFF) { CosPhiOne = sqrt( MAX(0.0, (1.0 - SQR(Xp1 + Xp2))) ); } else { CosPhiOne = 0.1; } /* osculating cone angle at L1 */ /* Chanan et al. (1972), ApJ 208, 512 */ if ( MIN(Binary[Primary].Mq, Binary[Secondary].Mq) > 0.4) CosPhiLagrange = cos (DTOR*57.88); else if ( MIN(Binary[Primary].Mq, Binary[Secondary].Mq) > 0.2) CosPhiLagrange = cos (DTOR*59.00); else if ( MIN(Binary[Primary].Mq, Binary[Secondary].Mq) > 0.1) CosPhiLagrange = cos (DTOR*60.56); else if ( MIN(Binary[Primary].Mq, Binary[Secondary].Mq) > 0.01) CosPhiLagrange = cos (DTOR*67.17); else CosPhiLagrange = 0.0; /* exact contact system */ if ( (fabs(Binary[Primary].RocheFill - 1.0) <= FLT_EPSILON) && (fabs(Binary[Secondary].RocheFill - 1.0) <= FLT_EPSILON) ) { CosPhiLagrange = CosPhiLagrange; } else { CosPhiLagrange = CosPhiOne; } Binary[Primary].CosPhiLagrange = CosPhiLagrange; Binary[Secondary].CosPhiLagrange = CosPhiLagrange; #ifdef HAVE_DISK /* * CosPhiLagrange for Disk; same formula as for CosPhiOne, but * now with the outer disk radius. */ Xp2 = Binary[Disk].Rout * Binary[Secondary].RLag1; Binary[Disk].CosPhiLagrange = sqrt( MAX(0.0, (1.0 - SQR(Xp1 + Xp2))) ); #endif return(0); } #ifdef HAVE_DISK void SecondaryEclipsedByDisk (double CosPhase2, double SinPhase2); void PrimaryEclipsedByDisk (double CosPhase, double SinPhase); void DiskEclipsedByDisk (double CosPhase2, double SinPhase2); /******************************************************************** @package nightfall @author Patrick Risse (risse@astro.uni-tuebingen.de) @version 1.0 @short Eclipse Verification for a CB-System with a Disk @param (int) j The phase index @return (void) @heading Light Curve ********************************************************************/ void LightCurveDisk(int j) { long NElem, NElemP, NElemS, NElemD; /* # of surface elements */ long i; /* loop counter */ int eclipseS=0; int eclipseP=0; int eclipseD=0; int count=0; int invis=0; int intersect = 0; int testsec = 0; int tested = 0; int EclipseImpossible = 0; int eclipsed = 0; /* test variable */ const int No = 0; /* test value */ const int Invisible = 10; /* test value */ int MinErr; /* exit status MinFinder */ double Phase, Phase2; /* orbital phase */ double OmegaDotPrimary; /* surface rotation async */ double OmegaDotSecondary; /* surface rotation async */ double lzero, mzero, nzero; /* LOS vector */ double lzero2, mzero2, nzero2; /* LOS vector */ double x_z, y_z, z_z; double x_tl, y_tl, z_tl; double Qi, Li, QLi; /* for sphere test */ double t_l, t_1, t_2, t_dist; /* intersection w/ sphere */ double tmin = 0; /* minimum along LOS */ double PotMin; /* minimum along LOS */ double SqrXp1, SqrXp2; /* square of (center - lagrange one) */ double SqrPol1,SqrPol2; /* square of polar radius */ double Sini, Cosi; /* cos, sin of orbital inclination */ double SinPhase, CosPhase; /* cos, sin of phase (primary) */ double SinPhase2, CosPhase2; /* cos, sin of phase (secondary) */ double CosPhiLagrange; /* tangent angle at Lagrange One (cos) */ double Xp1, Xp2; /* max. distance surface-center of star */ double Pot_One, Pot_Two; /* surface potentials */ SurfaceElement *SurfPtrD; /* pointer to surface Disk */ BinaryComponent *BinPtr; /* pointer to binary */ double fratio1 = 1.0; /* async rotation */ double fratio2 = 1.0; /* async rotation */ /* COMMENTS IN ENGLISH PLEASE !!!!!!!!! MK */ /*double vec_e1[3];*/ /* Ebenenvektor 1 */ /*double vec_e2[3];*/ /* Ebenenvektor 2 */ /*double vec_p1[3];*/ /* Eckpunkt der Ebene */ /*double vec_en[3];*/ /* Normalenvektor der Ebene */ /*double vec_g1[3];*/ /* Geradenvektor */ /*double vec_p2[3];*/ /* Eckpunkt der Gerade */ /*double vec_dif_p2p1[3];*/ /* Differenze der beiten Ortsvektoren */ /*double det1,det2;*/ /* Determinanten */ /* double vec_intersection[3];*/ /* Schnittpunkt zwischen Gerade und Ebene */ /* double Distance[3]; */ /* Distanz zwischen Schnittpunkt und Flaechenmittelpunkt */ /* double r; */ /* Faktor */ /* double rad_angle,angle; */ /* Angle between two vectors */ /* Recalculate the position of the Disk */ //DivideDisk(Disk,j); if (Flags.asynchron1 == ON) fratio1 = Binary[Primary].Fratio; if (Flags.asynchron2 == ON) fratio2 = Binary[Secondary].Fratio; OmegaDotPrimary = fratio1 * 2.0 * PI / Orbit.TruePeriod; OmegaDotSecondary = fratio2 * 2.0 * PI / Orbit.TruePeriod; Sini = sin(Orbit.Inclination); if (fabs(Sini) <= DBL_EPSILON) Sini = 0.0; Cosi = cos(Orbit.Inclination); if (fabs(Cosi) <= DBL_EPSILON) Cosi = 0.0; NElem = Binary[Primary].NumElem; Mq = Binary[Primary].Mq; /* Mq is global variable */ Pot_One = Binary[Primary].RochePot; Pot_Two = Binary[Secondary].RochePot; Xp1 = Binary[Primary].Xp; Xp2 = Binary[Secondary].Xp; CosPhiLagrange = Binary[Disk].CosPhiLagrange; SqrXp1 = SQR(Xp1); SqrXp2 = SQR(Xp2); SqrPol1 = SQR(Binary[Primary].Radius); /* square of polar radius */ SqrPol2 = SQR(Binary[Secondary].Radius); /* ------------------ eclipse test ------------------ */ /* to allow for generalization, we dont insist on */ /* same number of surface elements on both components */ NElem = IMAX(Binary[Primary].NumPlus, Binary[Secondary].NumPlus); NElemP = IMAX(Binary[Primary].NumElem, Binary[Primary].NumPlus); NElemS = IMAX(Binary[Secondary].NumElem, Binary[Secondary].NumPlus); NElemD = Binary[Disk].NumElem; /* We start at 90 deg, Phase = -0.25, both stars fully visible */ if (Flags.elliptic == OFF) { Phase = (PI/2 + (2*PI/PhaseSteps) * j ); FluxOut[j].Phase = (float) Phase; Orbit.Phase = Phase; Orbit.Nu = Phase; } else { Phase = Orbit.Phase; /* Orbit.OmegaZero is mean anomaly at secondary eclipse */ FluxOut[j].Phase = (float) (Orbit.M + PI + Orbit.OmegaZero); } SinPhase = sin(Phase); if (fabs(SinPhase) <= DBL_EPSILON) SinPhase = 0.0; CosPhase = cos(Phase); if (fabs(CosPhase) <= DBL_EPSILON) CosPhase = 0.0; /* vector towards observer (LOS = Line Of Sight) */ lzero = Sini * CosPhase; mzero = Sini * SinPhase; nzero = Cosi; /* Secondary turn by PI and mirror y->(-y) */ /* this is equivalent to mirror by x->(-x) */ Phase2 = (PI + Phase); SinPhase2 = sin(Phase2); if (fabs(SinPhase2) <= DBL_EPSILON) SinPhase2 = 0.0; CosPhase2 = cos(Phase2); if (fabs(CosPhase2) <= DBL_EPSILON) CosPhase2 = 0.0; /* vector towards observer (LOS = Line Of Sight) */ lzero2 = Sini * CosPhase2; mzero2 = Sini * SinPhase2; nzero2 = Cosi; SurfPtrD = Surface[Disk]; BinPtr = &Binary[Disk]; /* >>>>>>>>>>>>>>>>> the disk <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */ /****** NEW CODE by MPR 29.04.2002 **********************/ /* Disk surrounds Secondary. Non-tilted display in Gnuplot ok. */ /* loop over surface elements */ for (i = 0; i < NElemD; ++i) { #if 0 if (fabs(SurfPtrD->lambda) <= DBL_EPSILON) SurfPtrD->lambda = 0.0; if (fabs(SurfPtrD->mu) <= DBL_EPSILON) SurfPtrD->mu = 0.0; if (fabs(SurfPtrD->nu) <= DBL_EPSILON) SurfPtrD->nu = 0.0; #endif /* gamma = angle surface normal - LOS */ SurfPtrD->CosGamma = ( lzero2 * SurfPtrD->l + mzero2 * (-SurfPtrD->m) + nzero2 * SurfPtrD->n ); /* if (i == 0) { printf ("lzero2: %5.3f mzero2: %5.3f nzero2: %5.3f CosGamma: %5.3f\n", lzero2, mzero2, nzero2, SurfPtrD->CosGamma); } */ /*=========================================================== */ /* Default is visible */ eclipsed = No; /* Start checking. */ /* 1. Does the surface normal point away from observer ? * (Tested, works.) */ if (SurfPtrD->CosGamma < DBL_EPSILON ) { eclipsed = Invisible; } if (eclipsed == No ) { /* 2. Is element eclipsed by secondary ? * (Tested, works for non-inclined disk.) * * Test for an intersection with surface of Secondary. * Both Secondary and disk are at (0,0,0). Phase is shifted * by PI, and y->(-y). */ x_z = SurfPtrD->lambda; y_z = (-SurfPtrD->mu); z_z = SurfPtrD->nu; Qi = 2.0 * ( x_z*lzero2 + y_z*mzero2 + z_z*nzero2 ); Li = x_z*x_z + y_z*y_z + z_z*z_z - (Xp2)*(Xp2); QLi = Qi*Qi - 4.0*Li; t_l = sqrt(QLi); t_1 = (-Qi+t_l)/2; t_2 = (-Qi-t_l)/2; /* If true there is an intersection with the * sphere of the secondary star. Further tests * are required */ if (QLi <= 0.0) { /* no intersection */; } else if (t_1 > 0.0 && t_2 > 0.0) { /* check whether distance is greater than polar radius */ /* -- maximum inscribed sphere */ t_l = (t_1 + t_2)/2.0; /* t_l is Midpoint of ray intersecting Xp2-Sphere */ x_tl = x_z + t_l * lzero2; y_tl = y_z + t_l * mzero2; z_tl = z_z + t_l * nzero2; t_dist = SQR(x_tl) + y_tl * y_tl + z_tl * z_tl; if (t_dist <= (SqrPol2-DBL_EPSILON)) { /* definitely eclipsed */ eclipsed = EclipsedBySecondary; } else { /* Find Potential Minimum between t_1 and t_2 * final resort * changed Binary[Pri].Mq to Binary[Sec].Mq Jul 6 */ MinErr = MinFinder(&Binary[Secondary].Mq, &fratio2, &t_1, &t_2, &lzero, &mzero, &nzero, &x_z, &y_z, &z_z, &tmin, &PotMin); if (MinErr == 1 && Flags.fill == OFF) WARNING(_("LightCurve: No Minimum Found")); PotMin = -PotMin; SurfPtrD->LOS_Pot = (float) PotMin; /* added condition tmin >= 0 */ if (PotMin >= Pot_Two && tmin >= 0) { eclipsed = EclipsedBySecondary; } } } /* QLi > 0.0 */ } if (CosPhase2 > (CosPhiLagrange-DBL_EPSILON)) { ++testsec; if (eclipsed == No) { /* 3. Is the element eclipsed by the primary? * (Tested, works for non-inclined disk.) */ x_z = SurfPtrD->lambda; y_z = (-SurfPtrD->mu); z_z = SurfPtrD->nu; Qi = 2.0*( x_z*lzero2 + y_z*mzero2 + z_z*nzero2 - lzero2); Li = 1.0 + x_z*x_z + y_z*y_z + z_z*z_z - 2*x_z - SqrXp1; QLi = Qi*Qi - 4.0*Li; t_l = sqrt(QLi); /* just to hold temporary result */ t_1 = (-Qi + t_l)/2.0; t_2 = (-Qi - t_l)/2.0; if (QLi <= 0.0) { /* no intersection */; } else if (t_1 > 0.0 && t_2 > 0.0) /* intersection */ { ++intersect; /* check whether distance is greater than polar radius */ /* -- maximum inscribed sphere */ t_l = sqrt(QLi); /* just to hold temporary result */ t_1 = (-Qi + t_l)/2.0; t_2 = (-Qi - t_l)/2.0; t_l = (t_1 + t_2)/2.0; /* t_l is Midpoint of ray intersecting Xp1-Sphere */ x_tl = x_z + t_l * lzero2; y_tl = y_z + t_l * mzero2; z_tl = z_z + t_l * nzero2; t_dist = SQR(x_tl-1.0) + y_tl * y_tl + z_tl * z_tl; /* we do not need the sqrt, compare the squares */ if (t_dist <= (SqrPol1-DBL_EPSILON)) { /* definitely eclipsed */ eclipsed = EclipsedByPrimary ; } else { ++tested; /* Find Potential Minimum between t_1 and t_2 */ /* final resort */ /* changed Binary[Pri].Mq to Binary[Sec].Mq Jul 6 */ MinErr = MinFinder(&Binary[Primary].Mq, &fratio1, &t_1, &t_2, &lzero2, &mzero2, &nzero2, &x_z, &y_z, &z_z, &tmin, &PotMin); if (MinErr == 1) WARNING(_("LightCurve: No Minimum Found")); PotMin = -PotMin; SurfPtrD->LOS_Pot = (float) PotMin; /* added condition tmin >= 0 */ if (PotMin >= Pot_One && tmin >= 0) { eclipsed = EclipsedByPrimary; } } /* t_dist > (SqrPol1-DBL_EPSILON */ } /* QLi > 0.0) */ } /* (eclipsed < 1 )*/ } else { ++EclipseImpossible; } if (eclipsed == EclipsedByPrimary) eclipseP = eclipseP+1; if (eclipsed == EclipsedByDisk) eclipseD = eclipseD+1; if (eclipsed == EclipsedBySecondary) eclipseS = eclipseS+1; if (eclipsed == Invisible) invis = invis+1; if (eclipsed == No) { count=count+1; } SurfPtrD->EclipseFlag = eclipsed; if (SurfPtrD->EclipseFlag > 0 ) SurfPtrD->visibility = 0; else SurfPtrD->visibility = 1; if (eclipsed == 10) SurfPtrD->CosGamma = -1; ++SurfPtrD; } /* end loop over surface elements */ if (Flags.debug[GEO]) { printf("=========================================\n"); printf("Step %02d\n", j); printf("Total %d visible %d invisible %d eclipsed_tot %d\n", count+invis+eclipseP+eclipseS+eclipseD, count,invis,eclipseP+eclipseS+eclipseD); printf("eclipsed: by Prim %d, by Sec %d, by Disk %d\n", eclipseP,eclipseS,eclipseD); printf("testsec : %d\n", testsec); printf("intersect: %d\n", intersect); printf("tested : %d\n", tested); printf("imposs. : %d\n", EclipseImpossible); } count=0; eclipseP=0; eclipseS=0; invis=0; intersect=0; SecondaryEclipsedByDisk (CosPhase2, SinPhase2); DiskEclipsedByDisk (CosPhase2, SinPhase2); if (-CosPhase2 > (CosPhiLagrange-DBL_EPSILON)) PrimaryEclipsedByDisk (CosPhase, SinPhase); return; } /* end LightCurveDisk*/ #endif /******************************************************************** @package nightfall @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de) @version 1.0 @short Eclipse Verification @param (int) j The phase index @return (void) @heading Light Curve ********************************************************************/ void LightCurve(int j) { long NElem, NElemP, NElemS; /* # of surface elements */ long i; /* loop counter */ int eclipsed = 1; /* test variable */ const int Yes = 1; /* test value */ const int No = 0; /* test value */ const int Invisible = 10; /* test value */ int MinErr; /* exit status MinFinder */ double Phase, Phase2; /* orbital phase */ double OmegaDotPrimary; /* surface rotation async */ double OmegaDotSecondary; /* surface rotation async */ double lzero, mzero, nzero; /* LOS vector */ double lzero2, mzero2, nzero2; /* LOS vector */ double x_z, y_z, z_z; double x_tl, y_tl, z_tl; double Qi, Li, QLi; /* for sphere test */ double t_l, t_1, t_2, t_dist; /* intersection w/ sphere */ double tmin = 0; /* minimum along LOS */ double PotMin; /* minimum along LOS */ double SqrXp1, SqrXp2; /* square of (center - lagrange one) */ double SqrPol1,SqrPol2; /* square of polar radius */ double Sini, Cosi; /* cos, sin of orbital inclination */ double SinPhase, CosPhase; /* cos, sin of phase (primary) */ double SinPhase2, CosPhase2; /* cos, sin of phase (secondary) */ double Cosgamma; /* angle between LOS and surface normal */ double CosPhiLagrange; /* tangent angle at Lagrange One (cos) */ double Xp1, Xp2; /* max. distance surface-center of star */ double Pot_One, Pot_Two; /* surface potentials */ SurfaceElement *SurfPtrP; /* pointer to surface Primary */ SurfaceElement *SurfPtrS; /* pointer to surface Secondary */ double fratio1 = 1.0; /* async rotation */ double fratio2 = 1.0; /* async rotation */ if (Flags.asynchron1 == ON) fratio1 = Binary[Primary].Fratio; if (Flags.asynchron2 == ON) fratio2 = Binary[Secondary].Fratio; OmegaDotPrimary = fratio1 * 2.0 * PI / Orbit.TruePeriod; OmegaDotSecondary = fratio2 * 2.0 * PI / Orbit.TruePeriod; Sini = sin(Orbit.Inclination); Cosi = cos(Orbit.Inclination); NElem = Binary[Primary].NumElem; Mq = Binary[Primary].Mq; /* Mq is global variable */ Pot_One = Binary[Primary].RochePot; Pot_Two = Binary[Secondary].RochePot; Xp1 = Binary[Primary].Xp; Xp2 = Binary[Secondary].Xp; /* CosPhiLagrange is same for Primary and Secondary */ CosPhiLagrange = Binary[Primary].CosPhiLagrange; SqrXp1 = SQR(Xp1); SqrXp2 = SQR(Xp2); SqrPol1 = SQR(Binary[Primary].Radius); /* square of polar radius */ SqrPol2 = SQR(Binary[Secondary].Radius); /* ------------------ eclipse test ------------------ */ /* to allow for generalization, we dont insist on */ /* same number of surface elements on both components */ NElem = IMAX(Binary[Primary].NumPlus,Binary[Secondary].NumPlus); NElemP = IMAX(Binary[Primary].NumElem, Binary[Primary].NumPlus); NElemS = IMAX(Binary[Secondary].NumElem, Binary[Secondary].NumPlus); /* NElem = IMAX(Binary[Primary].NumElem,Binary[Secondary].NumElem); NElemP = Binary[Primary].NumElem; NElemS = Binary[Secondary].NumElem; */ /* We start at 90 deg, Phase = -0.25, both stars fully visible */ if (Flags.elliptic == OFF) { Phase = (PI/2 + (2*PI/PhaseSteps) * j ); FluxOut[j].Phase = (float) Phase; Orbit.Phase = Phase; Orbit.Nu = Phase; } else { Phase = Orbit.Phase; /* Orbit.OmegaZero is mean anomaly at secondary eclipse */ FluxOut[j].Phase = (float) (Orbit.M + PI + Orbit.OmegaZero); } SinPhase = sin(Phase); CosPhase = cos(Phase); /* vector towards observer (LOS = Line Of Sight) */ lzero = Sini * CosPhase; mzero = Sini * SinPhase; nzero = Cosi; /* Secondary turn by PI and mirror y->(-y) */ /* this is equivalent to mirror by x->(-x) */ Phase2 = (PI + Phase); SinPhase2 = sin(Phase2); CosPhase2 = cos(Phase2); /* vector towards observer (LOS = Line Of Sight) */ lzero2 = Sini * CosPhase2; mzero2 = Sini * SinPhase2; nzero2 = Cosi; /* test whether eclipse of Primary definitely impossible */ /* this flag is (Flags.InEclipse1) required for fractional visibility */ if ( CosPhase <= (CosPhiLagrange-DBL_EPSILON) ) { Flags.InEclipse1 = OFF; } else { Flags.InEclipse1 = ON; } /* test whether eclipse of Secondary definitely impossible */ /* this flag (Flags.InEclipse2) is required for fractional visibility */ if ( CosPhase2 <= (CosPhiLagrange-DBL_EPSILON) ) { Flags.InEclipse2 = OFF; } else { Flags.InEclipse2 = ON; } SurfPtrP = Surface[Primary]; SurfPtrS = Surface[Secondary]; for (i = 0; i < NElem; ++i) { /* loop over surface elements */ /* START WITH PRIMARY */ if (i < NElemP) { SurfPtrP->Velocity = - OmegaDotPrimary * Orbit.TrueDistance * Orbit.Dist * (-SurfPtrP->mu * lzero + SurfPtrP->lambda * mzero); SurfPtrP->LOS_Pot = 0.0; /* angle surface normal - LOS */ Cosgamma = lzero * SurfPtrP->l + mzero * SurfPtrP->m + nzero * SurfPtrP->n; /* IF (Cosgamma > 0) THEN Visible Side, Test for Eclipse */ if (Cosgamma <= 0) eclipsed = Invisible; /* invisible */ if (Cosgamma >= DBL_EPSILON){ /* Visible Side */ eclipsed = Yes; /* test whether eclipse definitely impossible */ if (CosPhase <= (CosPhiLagrange-DBL_EPSILON)) eclipsed = No; if (Flags.fill == ON) { if (i > Binary[Primary].NumElem) eclipsed = Yes; } /* IF (eclipsed == Yes) MORE TESTING */ if (eclipsed == Yes) { /* intersect with circle of rad Xp2 centered on second */ /* -- minimum sphere containing the star */ x_z = SurfPtrP->lambda; y_z = SurfPtrP->mu; z_z = SurfPtrP->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 - SqrXp2; QLi = Qi*Qi - 4.0*Li; if (QLi <= 0.0) { /* No intersection */ eclipsed = No; } else { /* Intersection, next test */ /* check whether distance is greater than polar radius */ /* -- maximum inscribed sphere */ t_l = sqrt(QLi); /* just to hold temporary result */ t_1 = (-Qi + t_l)/2.0; t_2 = (-Qi - t_l)/2.0; t_l = (t_1 + t_2)/2.0; /* t_l is Midpoint of ray intersecting Xp2-Sphere */ x_tl = x_z + t_l * lzero; y_tl = y_z + t_l * mzero; z_tl = z_z + t_l * nzero; t_dist = SQR(x_tl-1.0) + y_tl * y_tl + z_tl * z_tl; if (Flags.fill == ON) { if (i > Binary[Primary].NumElem) { /* rw 06.02.2007 */ /* eclipsed = Yes; */ /* superfluous */ /* incorrect, will test against wrong star * t_1 = MAX(t_1, t_2); * t_2 = t_1 + 1.0; */ /* next test would yield incorrect result, thus jump * to the MinFinder() routine */ t_dist = SqrPol2+1.; } } if (t_dist <= (SqrPol2-DBL_EPSILON)) { eclipsed = Yes; /* definitely eclipsed */ } else { /* Find Potential Minimum between t_1 and t_2 */ /* final resort */ /* changed Binary[Pri].Mq to Binary[Sec].Mq Jul 6 */ MinErr = MinFinder(&Binary[Secondary].Mq, &fratio2, &t_1, &t_2, &lzero, &mzero, &nzero, &x_z, &y_z, &z_z, &tmin, &PotMin); if (MinErr == 1 && Flags.fill == OFF) WARNING(_("LightCurve: No Minimum Found")); PotMin = -PotMin; SurfPtrP->LOS_Pot = (float) PotMin; /* added condition tmin >= 0 */ if (PotMin >= Pot_Two && tmin >= 0) { eclipsed = Yes; } else { eclipsed = No; } } } } SurfPtrP->EclipseFlag = eclipsed; SurfPtrP->CosGamma = (float) Cosgamma; if (eclipsed > 0) SurfPtrP->visibility = 0; else SurfPtrP->visibility = 1; } else { /* invisible side */ SurfPtrP->EclipseFlag = eclipsed; SurfPtrP->CosGamma = -1; if (eclipsed > 0) SurfPtrP->visibility = 0; else SurfPtrP->visibility = 1; } ++SurfPtrP; } /* END PRIMARY */ /* HERE SECONDARY */ if (i < NElemS) { SurfPtrS->Velocity = - OmegaDotSecondary * Orbit.TrueDistance * Orbit.Dist * (SurfPtrS->mu * lzero2 + SurfPtrS->lambda * mzero2); SurfPtrS->LOS_Pot = 0.0; /* angle surface normal - LOS */ /* dont forget to mirror in y->(-y) */ Cosgamma = lzero2 * SurfPtrS->l + mzero2 * (-SurfPtrS->m) + nzero2 * SurfPtrS->n; /* IF (Cosgamma > 0) THEN Visible Side, Test for Eclipse */ if (Cosgamma <= 0) eclipsed = Invisible; /* invisible */ if (Cosgamma >= DBL_EPSILON){ /* Visible Side */ eclipsed = Yes; /* test whether eclipse definitely impossible */ if (CosPhase2 <= (CosPhiLagrange-DBL_EPSILON)) eclipsed = No; if (Flags.fill == ON) { if (i > Binary[Secondary].NumElem) eclipsed = Yes; } /* IF (eclipsed == Yes) MORE TESTING */ if (eclipsed == Yes) { /* intersect with circle of rad Xp1 centered on prim */ /* -- minimum sphere containing the star */ x_z = SurfPtrS->lambda; y_z = (-SurfPtrS->mu); z_z = SurfPtrS->nu; Qi = 2.0*( x_z*lzero2 + y_z*mzero2 + z_z*nzero2 - lzero2); Li = 1.0 + x_z*x_z + y_z * y_z + z_z*z_z - 2*x_z - SqrXp1; QLi = Qi*Qi - 4.0*Li; if (QLi <= 0.0) { /* No intersection */ eclipsed = No; } else { /* Intersection, next test */ /* check whether distance is greater than polar radius */ /* -- maximum inscribed sphere */ t_l = sqrt(QLi); /* just to hold temporary result */ t_1 = (-Qi + t_l)/2.0; t_2 = (-Qi - t_l)/2.0; t_l = (t_1 + t_2)/2.0; /* t_l is Midpoint of ray intersecting Xp2-Sphere */ x_tl = x_z + t_l * lzero2; y_tl = y_z + t_l * mzero2; z_tl = z_z + t_l * nzero2; t_dist = SQR(x_tl-1.0) + y_tl*y_tl + z_tl*z_tl; if (Flags.fill == ON) { if (i > Binary[Secondary].NumElem) { /* rw 06.02.2007 */ /* eclipsed = Yes; */ /* superfluous */ /* incorrect, will test against wrong star * t_1 = MAX(t_1, t_2); * t_2 = t_1 + 1.0; */ /* next test would yield incorrect result, thus jump * to the MinFinder() routine */ t_dist = SqrPol1+1.; } } /* we do not need the sqrt, compare the squares */ if (t_dist <= (SqrPol1-DBL_EPSILON)) { eclipsed = Yes; /* definitely eclipsed */ } else { /* Find Potential Minimum between t_1 and t_2 */ /* final resort */ /* changed Binary[Sec].Mq to Binary[Pri].Mq Jul 6 */ MinErr = MinFinder(&Binary[Primary].Mq, &fratio1, &t_1, &t_2, &lzero2, &mzero2, &nzero2, &x_z, &y_z, &z_z, &tmin, &PotMin); if (MinErr == 1) WARNING(_("LightCurve: No Minimum Found")); PotMin = -PotMin; SurfPtrS->LOS_Pot = (float) PotMin; /* added condition tmin >= 0 */ if (PotMin >= (Pot_One-DBL_EPSILON) && tmin >= 0) { eclipsed = Yes; } else { eclipsed = No; } } } } SurfPtrS->EclipseFlag = eclipsed; SurfPtrS->CosGamma = (float) Cosgamma; if (eclipsed > 0) SurfPtrS->visibility = 0; else SurfPtrS->visibility = 1; } else { SurfPtrS->EclipseFlag = eclipsed; SurfPtrS->CosGamma = -1; if (eclipsed > 0) SurfPtrS->visibility = 0; else SurfPtrS->visibility = 1; } ++SurfPtrS; } /* END SECONDARY */ } /* end loop over surface elements */ return; } /* end LightCurve */ /******************************************************************** @package nightfall @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de) @version 1.0 @short Flux computation @param (int) Comp The stellar component @param (int) Phase The phase index @return (void) @heading Light Curve ********************************************************************/ void LightFlux(int Comp, int Phase) { register long i; /* loop counter */ double Cosgamma; /* LOS angle */ double Common; /* common part */ double Inverse; /* if limb brightening */ double Weight; /* weight factor */ double SumOfVelocity; /* normalization */ double SumOfWeight; /* normalization */ double OneMinCosgamma; /* 1.0 - Cosgamma */ SurfaceElement *SurfPtr; /* pointer to surface */ long nelem; /* loop limit */ PhotoOut *OFlux; /* pointer to outflux */ /* ------------------ flux computing --------------- */ /* We start at 90 deg, Phase = 0.25, both stars fully visible */ /* loops over passbands are unrolled, all if's in loop eliminated */ SumOfVelocity = 0.0; SumOfWeight = 0.0; nelem = Binary[Comp].NumElem; OFlux = &FluxOut[Phase]; SurfPtr = Surface[Comp]; if (Flags.limb == 0 || Flags.limb == 3) { for (i = 0; i < nelem; ++i) { /* loop over surface elements */ /* add up Flux of visible elements */ /* IF NOT eclipsed THEN */ if (SurfPtr->CosGamma >= 0.0) { Cosgamma = SurfPtr->CosGamma; Common = SurfPtr->visibility * Cosgamma; if (Flags.limb == 0) Inverse = 1.0; else Inverse = (SurfPtr->AreaType == FRONTSIDE) ? -3.0 : 1.0; OneMinCosgamma = (1.0 - Cosgamma); Weight = SurfPtr->f_[0] * Common * (1.0 - Inverse * Limb[Comp][0][0]*OneMinCosgamma); OFlux->Mag[0] = OFlux->Mag[0] + Weight; Weight = SurfPtr->f_[1] * Common * (1.0 - Inverse * Limb[Comp][1][0]*OneMinCosgamma); OFlux->Mag[1] = OFlux->Mag[1] + Weight; Weight = SurfPtr->f_[2] * Common * (1.0 - Inverse * Limb[Comp][2][0]*OneMinCosgamma); OFlux->Mag[2] = OFlux->Mag[2] + Weight; SumOfVelocity = SumOfVelocity + Weight * SurfPtr->Velocity; SumOfWeight = SumOfWeight + Weight; Weight = SurfPtr->f_[3] * Common * (1.0 - Inverse * Limb[Comp][3][0]*OneMinCosgamma); OFlux->Mag[3] = OFlux->Mag[3] + Weight; Weight = SurfPtr->f_[4] * Common * (1.0 - Inverse * Limb[Comp][4][0]*OneMinCosgamma); OFlux->Mag[4] = OFlux->Mag[4] + Weight; Weight = SurfPtr->f_[5] * Common * (1.0 - Inverse * Limb[Comp][5][0]*OneMinCosgamma); OFlux->Mag[5] = OFlux->Mag[5] + Weight; Weight = SurfPtr->f_[6] * Common * (1.0 - Inverse * Limb[Comp][6][0]*OneMinCosgamma); OFlux->Mag[6] = OFlux->Mag[6] + Weight; Weight = SurfPtr->f_[7] * Common * (1.0 - Inverse * Limb[Comp][7][0]*OneMinCosgamma); OFlux->Mag[7] = OFlux->Mag[7] + Weight; Weight = SurfPtr->f_[8] * Common * (1.0 - Inverse * Limb[Comp][8][0]*OneMinCosgamma); OFlux->Mag[8] = OFlux->Mag[8] + Weight; Weight = SurfPtr->f_[9] * Common * (1.0 - Inverse * Limb[Comp][9][0]*OneMinCosgamma); OFlux->Mag[9] = OFlux->Mag[9] + Weight; Weight = SurfPtr->f_[10] * Common * (1.0 - Inverse * Limb[Comp][10][0]*OneMinCosgamma); OFlux->Mag[10] = OFlux->Mag[10] + Weight; Weight = SurfPtr->f_[11] * Common * (1.0 - Inverse * Limb[Comp][10][0]*OneMinCosgamma); OFlux->Mag[11] = OFlux->Mag[11] + Weight; } ++SurfPtr; } } else if ( Flags.limb == 1) { for (i = 0; i < nelem; ++i) { /* loop over surface elements */ /* add up Flux of visible elements */ /* IF NOT eclipsed THEN */ if (SurfPtr->CosGamma >= 0.0) { Cosgamma = SurfPtr->CosGamma; Common = SurfPtr->visibility * Cosgamma; OneMinCosgamma = (1.0 - Cosgamma); Weight = SurfPtr->f_[0] * Common *(1.0 - Limb[Comp][0][0]*OneMinCosgamma - Limb[Comp][0][1]*(1.0- Cosgamma*Cosgamma)); OFlux->Mag[0] = OFlux->Mag[0] + Weight; Weight = SurfPtr->f_[1] * Common *(1.0 - Limb[Comp][1][0]*OneMinCosgamma - Limb[Comp][1][1]*(1.0- Cosgamma*Cosgamma)); OFlux->Mag[1] = OFlux->Mag[1] + Weight; Weight = SurfPtr->f_[2] * Common *(1.0 - Limb[Comp][2][0]*OneMinCosgamma - Limb[Comp][2][1]*(1.0- Cosgamma*Cosgamma)); OFlux->Mag[2] = OFlux->Mag[2] + Weight; SumOfVelocity = SumOfVelocity + Weight * SurfPtr->Velocity; SumOfWeight = SumOfWeight + Weight; Weight = SurfPtr->f_[3] * Common *(1.0 - Limb[Comp][3][0]*OneMinCosgamma - Limb[Comp][3][1]*(1.0- Cosgamma*Cosgamma)); OFlux->Mag[3] = OFlux->Mag[3] + Weight; Weight = SurfPtr->f_[4] * Common *(1.0 - Limb[Comp][4][0]*OneMinCosgamma - Limb[Comp][4][1]*(1.0- Cosgamma*Cosgamma)); OFlux->Mag[4] = OFlux->Mag[4] + Weight; Weight = SurfPtr->f_[5] * Common *(1.0 - Limb[Comp][5][0]*OneMinCosgamma - Limb[Comp][5][1]*(1.0- Cosgamma*Cosgamma)); OFlux->Mag[5] = OFlux->Mag[5] + Weight; Weight = SurfPtr->f_[6] * Common *(1.0 - Limb[Comp][6][0]*OneMinCosgamma - Limb[Comp][6][1]*(1.0- Cosgamma*Cosgamma)); OFlux->Mag[6] = OFlux->Mag[6] + Weight; Weight = SurfPtr->f_[7] * Common *(1.0 - Limb[Comp][7][0]*OneMinCosgamma - Limb[Comp][7][1]*(1.0- Cosgamma*Cosgamma)); OFlux->Mag[7] = OFlux->Mag[7] + Weight; Weight = SurfPtr->f_[8] * Common *(1.0 - Limb[Comp][8][0]*OneMinCosgamma - Limb[Comp][8][1]*(1.0- Cosgamma*Cosgamma)); OFlux->Mag[8] = OFlux->Mag[8] + Weight; Weight = SurfPtr->f_[9] * Common *(1.0 - Limb[Comp][9][0]*OneMinCosgamma - Limb[Comp][9][1]*(1.0- Cosgamma*Cosgamma)); OFlux->Mag[9] = OFlux->Mag[9] + Weight; Weight = SurfPtr->f_[10] * Common *(1.0 - Limb[Comp][10][0]*OneMinCosgamma - Limb[Comp][10][1]*(1.0- Cosgamma*Cosgamma)); OFlux->Mag[10] = OFlux->Mag[10] + Weight; Weight = SurfPtr->f_[11] * Common *(1.0 - Limb[Comp][11][0]*OneMinCosgamma - Limb[Comp][11][1]*(1.0- Cosgamma*Cosgamma)); OFlux->Mag[11] = OFlux->Mag[11] + Weight; } ++SurfPtr; } } else if ( Flags.limb == 2) { for (i = 0; i < nelem; ++i) { /* loop over surface elements */ /* add up Flux of visible elements */ /* IF NOT eclipsed THEN */ if (SurfPtr->CosGamma >= 0.0) { Cosgamma = SurfPtr->CosGamma; Common = SurfPtr->visibility * Cosgamma; OneMinCosgamma = (1.0 - Cosgamma); Weight = SurfPtr->f_[0] * Common *(1.0 - Limb[Comp][0][0]*OneMinCosgamma - Limb[Comp][0][1]*(1.0- sqrt(Cosgamma))); OFlux->Mag[0] = OFlux->Mag[0] + Weight; Weight = SurfPtr->f_[1] * Common *(1.0 - Limb[Comp][1][0]*OneMinCosgamma - Limb[Comp][1][1]*(1.0- sqrt(Cosgamma))); OFlux->Mag[1] = OFlux->Mag[1] + Weight; Weight = SurfPtr->f_[2] * Common *(1.0 - Limb[Comp][2][0]*OneMinCosgamma - Limb[Comp][2][1]*(1.0- sqrt(Cosgamma))); OFlux->Mag[2] = OFlux->Mag[2] + Weight; SumOfVelocity = SumOfVelocity + Weight * SurfPtr->Velocity; SumOfWeight = SumOfWeight + Weight; Weight = SurfPtr->f_[3] * Common *(1.0 - Limb[Comp][3][0]*OneMinCosgamma - Limb[Comp][3][1]*(1.0- sqrt(Cosgamma))); OFlux->Mag[3] = OFlux->Mag[3] + Weight; Weight = SurfPtr->f_[4] * Common *(1.0 - Limb[Comp][4][0]*OneMinCosgamma - Limb[Comp][4][1]*(1.0- sqrt(Cosgamma))); OFlux->Mag[4] = OFlux->Mag[4] + Weight; Weight = SurfPtr->f_[5] * Common *(1.0 - Limb[Comp][5][0]*OneMinCosgamma - Limb[Comp][5][1]*(1.0- sqrt(Cosgamma))); OFlux->Mag[5] = OFlux->Mag[5] + Weight; Weight = SurfPtr->f_[6] * Common *(1.0 - Limb[Comp][6][0]*OneMinCosgamma - Limb[Comp][6][1]*(1.0- sqrt(Cosgamma))); OFlux->Mag[6] = OFlux->Mag[6] + Weight; Weight = SurfPtr->f_[7] * Common *(1.0 - Limb[Comp][7][0]*OneMinCosgamma - Limb[Comp][7][1]*(1.0- sqrt(Cosgamma))); OFlux->Mag[7] = OFlux->Mag[7] + Weight; Weight = SurfPtr->f_[8] * Common *(1.0 - Limb[Comp][8][0]*OneMinCosgamma - Limb[Comp][8][1]*(1.0- sqrt(Cosgamma))); OFlux->Mag[8] = OFlux->Mag[8] + Weight; Weight = SurfPtr->f_[9] * Common *(1.0 - Limb[Comp][9][0]*OneMinCosgamma - Limb[Comp][9][1]*(1.0- sqrt(Cosgamma))); OFlux->Mag[9] = OFlux->Mag[9] + Weight; Weight = SurfPtr->f_[10] * Common *(1.0 - Limb[Comp][10][0]*OneMinCosgamma - Limb[Comp][10][1]*(1.0- sqrt(Cosgamma))); OFlux->Mag[10] = OFlux->Mag[10] + Weight; Weight = SurfPtr->f_[11] * Common *(1.0 - Limb[Comp][11][0]*OneMinCosgamma - Limb[Comp][11][1]*(1.0- sqrt(Cosgamma))); OFlux->Mag[11] = OFlux->Mag[11] + Weight; } ++SurfPtr; } } /* ------------ normalize velocity ---------------------- */ if (SumOfWeight >= DBL_EPSILON) { FluxOut[Phase].D_RV[Comp] = (float) (SumOfVelocity / SumOfWeight); } else { FluxOut[Phase].D_RV[Comp] = 0.0; } FluxOut[Phase].RV[Comp] = FluxOut[Phase].D_RV[Comp] + FluxOut[Phase].RV[Comp]; return; } /**************************************************************************** @package nightfall @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de) @version 1.0 @short Find minimum of the potential along some line-of-sight vector. @long Uses Brents algorithm, which is a combination of golden section search and parabolic interpolation. Adopted from procedure 'localmin' in: Richard Brent, Algorithms for minimization without derivatives, Prentice-Hall, Inc. (1973) @tip Input is a parametric form of LOS vector + brackets. @param (double*) q Mass ratio @param (double*) ff Async rotation @param (double*) t1 Minimum bracket @param (double*) t3 Minimum bracket @param (double*) l0 x-component LOS vector @param (double*) m0 y-component LOS vector @param (double*) n0 z-component LOS vector @param (double*) x0 x-component of point on LOS vector @param (double*) y0 y-component of point on LOS vector @param (double*) z0 z-component of point on LOS vector @param (double*) tmin The position of the minimum @param (double*) Cmin The LOS potential minimum @return (int) The error status @heading Light Curve ****************************************************************************/ int MinFinder(double *q, double *ff, double *t1, double *t3, double *l0, double *m0, double *n0, double *x0, double *y0, double *z0, double *tmin, /*@out@*/ double *Cmin) { register int i; /* loop variable */ double t2; /* middle value */ double tmp; /* temp var */ double step, stepstore; /* stepsizes */ double Fe = 0.0, Ff = 0.0; /* function values */ double Fg = 0.0, Fh = 0.0; /* function values */ double e = 0.0, f, g, h; /* positions (parametric) */ double delta1, delta2; /* test convergence */ double Start, End, Middle; /* brackets */ double nextstep = 0.0; /* stepsize */ double x, y, z; /* position in cartesian coordinates */ double P1 = 0.0, P2, P3; /* fit parabola */ double precise = 0.01; /* convergence criterion */ if ( (*t1) <= ((*t3)-DBL_EPSILON) ) { Start = (*t1); End = (*t3); } else { Start = (*t3); End = (*t1); } step = 0.0; /* we have a parametric description of a vector X(h) */ /* we want to find min(func) = min(X(h)) as function of h */ t2 = 0.62*(*t1) + 0.38*(*t3); /* initial value of X(h) */ x = 1.0 - ((*x0) + (*l0)*(t2)); y = -(*y0) - (*m0)*(t2); z = (*z0) + (*n0)*(t2); /* NEW jul 6,98 -- change of coordinate system, calculate */ /* in system of eclipsing star */ /* (i.e. x = 1.0 - x; y = -y; z = z;) */ /* t2 is the initial guess for the minimum between t1 and t3 */ h = g = f = t2; /* initial value of F */ tmp = 1.0-x; tmp = sqrt(tmp*tmp + y*y + z*z ); Ff = Fg = Fh = -( 1.0/sqrt(x * x + y * y + z * z) + (*q) * (1.0/tmp - x) + ( x * x + y * y ) * (*ff)*(*ff) * (0.5*( (*q) + 1.0))); /* we do a maximum of MAXITER iterations - convergence usually */ /* after less than ten iterations */ /*@i@*/ for (i = MAXITER; i; --i) { Middle = 0.5 * (Start + End); delta1 = precise * fabs(h) + DBL_EPSILON; delta2 = 2.0 * delta1; /* test for convergence */ /* RETURN if converged */ if (fabs(h-Middle) <= (delta2 - 0.5 * (End-Start) ) ) { *tmin = h; *Cmin = Fh; return(0); } if (fabs(step) >= delta1) { /* fit parabola */ P3 = (h-g)*(Fh-Ff); P2 = (h-f)*(Fh-Fg); P1 = (h-f)*P2 - (h-g)*P3; P2 = 2.0 * (P2 - P3); if (P2 >= DBL_EPSILON) { P1 = -P1; } else { P2 = -P2; } stepstore = step; step = nextstep; if ( (fabs(P1) >= fabs(0.5 * P2 * stepstore)) || (P1 <= P2*(Start-h)) || (P1 >= P2*(End-h))) { /* a golden-section step */ if (h >= Middle) { step = (Start-h); } else { step = (End-h); } nextstep = 0.381966 * step; } else { /* do a parabolic interpolation step */ nextstep = P1/P2; e = h + nextstep; /* f must not be evaluated too close to Start or End */ if ((e-Start) <= delta2 || (End-e) <= delta2) { if ((Middle-h) <= (-DBL_EPSILON)) { nextstep = nextstep-delta1; } else { nextstep = nextstep+delta1; } } } } else { /* a golden-section step */ if (h >= Middle) { step = (Start-h); } else { step = (End-h); } nextstep = 0.381966 * step; } if (fabs(nextstep) >= delta1) { e = h + nextstep; } else { /* f must not be evaluated too close */ if ( nextstep <= (-DBL_EPSILON) ) { e = h - delta1; } else { e = h + delta1; } } /* e ist the current best guess for the minimum */ /* calculate X(e) */ x = 1.0 - ((*x0) + (*l0) * e); y = - (*y0) - (*m0) * e; z = (*z0) + (*n0) * e; /* NEW jul 6 -- change of coordinate system, calculate */ /* in system of eclipsing star */ /* i.e x = 1.0 - x; y = -y; z = z; */ /* calculate f(X(e)) */ tmp = 1.0-x; tmp = sqrt(tmp*tmp + y*y + z*z ); Fe = -( 1.0/sqrt(x * x + y * y + z * z) + (*q) * (1.0/tmp - x) + ( x * x + y * y )* (*ff)*(*ff) * (0.5*( (*q) + 1.0))); /* update */ if (Fe <= Fh) { if ( e >= h ) Start = h; else End = h; f = g; g = h; h = e; Ff = Fg; Fg = Fh; Fh = Fe; } else { if ( e <= (h-DBL_EPSILON) ) Start = e; else End = e; if ( (Fe <= Fg) || (fabs (g - x) <= FLT_EPSILON) ) { f = g; g = e; Ff = Fg; Fg = Fe; } else if ( (Fe <= Ff) || (fabs (f- g) <= FLT_EPSILON) || (fabs (f - h) <= FLT_EPSILON) ) { f = e; Ff = Fe; } } /* end of main loop */ } /* only reach here if no convergence */ *tmin = h; *Cmin = Fh; return(1); } #ifdef HAVE_DISK /******************************************************************** @package nightfall @author Patrick Risse (risse@astro.uni-tuebingen.de) @version 1.0 @short Eclipse checking @param (int) Comp The stellar component @return (int) eclipsed @heading Light Curve ********************************************************************/ void eclipsedbydisk(double lzero2, double mzero2, double nzero2, long Diskelements, int Comp, SurfaceElement *SurfPtrD, int *eclipsed) { long Element_number = 0; /* Variable for counting */ double vec_e1[3]; /* Ebenenvektor 1 */ double vec_e2[3]; /* Ebenenvektor 2 */ double vec_p1[3]; /* Eckpunkt der Ebene */ /*double vec_en[3];*/ /* Normalenvektor der Ebene */ double vec_g1[3]; /* Geradenvektor*/ double vec_p2[3]; /* Eckpunkt der Gerade */ double vec_dif_p1p2[3]; /* Differenze der beiten Ortsvektoren */ double det1,det2; /* Determinanten */ double vec_intersection[3];/* Schnittpunkt zwischen Gerade und Ebene */ double Distance[3]; /* Distanz zwischen Schnittpunkt und Flaechenmittelpunkt */ double r; /* Faktor */ double angle; /* Angle between two vectors */ double delta = 0; double AreaAngle1,AreaAngle2; /* */ const int Invisible = 10; SurfaceElement *Store; BinaryComponent *BinPtr; /* pointer to binary */ BinPtr = &Binary[Disk]; Store=SurfPtrD; SurfPtrD=Surface[Disk]; //printf("ELEMENT: %d\n",Store->MyNum); while (Element_number < Diskelements && *eclipsed < 1) { if ((Element_number != Store->MyNum) && (SurfPtrD->EclipseFlag != Invisible)) { /* Gleichungssystem basteln */ vec_p1[0]=SurfPtrD->lambda; vec_p1[1]=SurfPtrD->mu; vec_p1[2]=SurfPtrD->nu; vec_e1[0]=SurfPtrD->VecE1[0]; vec_e1[1]=SurfPtrD->VecE1[1]; vec_e1[2]=SurfPtrD->VecE1[2]; vec_e2[0]=SurfPtrD->VecE2[0]; vec_e2[1]=SurfPtrD->VecE2[1]; vec_e2[2]=SurfPtrD->VecE2[2]; /* Vektor fuer Gerade (Sichtstrahl) */ vec_g1[0]=lzero2; vec_g1[1]=mzero2; vec_g1[2]=nzero2; vec_p2[0]=Store->lambda; vec_p2[1]=Store->mu; vec_p2[2]=Store->nu; /* Wir haben nun folgendes Gleichungssystem */ /* X = vec_p1 * u * vec_e1 * v * vec_e2 = vec_p2 * r * vec_g1 */ /* Determinante des Gleichungssystems berechnen */ det1=vec_e1[0]*vec_e2[1]*vec_g1[2] +vec_e1[1]*vec_e2[2]*vec_g1[0] +vec_e1[2]*vec_e2[0]*vec_g1[1] -vec_e1[2]*vec_e2[1]*vec_g1[0] -vec_e1[0]*vec_e2[2]*vec_g1[1] -vec_e1[1]*vec_e2[0]*vec_g1[2]; if (fabs(det1) != 0.0) { //printf("det1 %20.20f Element: %d \n",det1,SurfPtrD->MyNum); /*Wenn die Determinante ungleich Null ist muss der Abstand des Schnittpunktes Element bestimmt werden */ vec_dif_p1p2[0]=vec_p1[0]-vec_p2[0]; vec_dif_p1p2[1]=vec_p1[1]-vec_p2[1]; vec_dif_p1p2[2]=vec_p1[2]-vec_p2[2]; det2=vec_e1[0]*vec_e2[1]*vec_dif_p1p2[2] +vec_e1[1]*vec_e2[2]*vec_dif_p1p2[0] +vec_e1[2]*vec_e2[0]*vec_dif_p1p2[1] -vec_e1[2]*vec_e2[1]*vec_dif_p1p2[0] -vec_e1[0]*vec_e2[2]*vec_dif_p1p2[1] -vec_e1[1]*vec_e2[0]*vec_dif_p1p2[2]; /* in Geraden Gleichung einsetzen und ausrechen ist der Abstand Groesser als das Flaechenelement so wird das Untersuchte Element nicht bedeckt */ r=det2/det1; //printf("verglichen mit Element %d r: %f\n",SurfPtrD->MyNum,r); if (r > 0.0 +DBL_EPSILON) { /* Damit liegt der Schnittpunkt bei: */ vec_intersection[0]=vec_p2[0]+r*vec_g1[0]; vec_intersection[1]=vec_p2[1]+r*vec_g1[1]; vec_intersection[2]=vec_p2[2]+r*vec_g1[2]; Distance[0]=vec_intersection[0];//-vec_p1[0]; Distance[1]=vec_intersection[1];//-vec_p1[1]; Distance[2]=vec_intersection[2];//-vec_p1[2]; /* Dann rechnen wir noch etwas aus */ if (Distance[1] >= 0.0 && Distance[0] >= 0.0) angle = ((atan(Distance[1]/Distance[0]))/(2.*PI))*360.; else if (Distance[1] >= 0.0 && Distance[0] < 0.0) angle = (((atan(Distance[1]/Distance[0]))/(2.*PI))*360.)+180.; else if (Distance[1] < 0.0 && Distance[0] <= 0.0) angle = (((atan(Distance[1]/Distance[0]))/(2.*PI))*360.)+180; else /* if (Distance[1] < 0.0 && Distance[0] >= 0.0) */ angle = (((atan(Distance[1]/Distance[0]))/(2.*PI))*360.)+360; AreaAngle1=SurfPtrD->AreaAngle + delta; AreaAngle2=SurfPtrD->AreaAngle - delta; if (AreaAngle2 < 0.0) AreaAngle2=AreaAngle2+360.0; //printf("angle %f\n",angle); delta=360.0/(2*BinPtr->Segments); if (SurfPtrD->AreaType == TOP_SEGMENT || SurfPtrD->AreaType == BOTTOM_SEGMENT) { if (sqrt(Distance[0]*Distance[0]+Distance[1]*Distance[1]) > SurfPtrD->RadiusIn && sqrt(Distance[0]*Distance[0]+Distance[1]*Distance[1]) <= SurfPtrD->RadiusOut ) { /* Vektor vom Ursprung zum Schnittpunkt bauen */ /* Da Ursprung (0/0/0) vektor = vec_intersection */ if (angle <= (SurfPtrD->AreaAngle + delta) && angle > (SurfPtrD->AreaAngle - delta)) { //printf("Vorher %d\n",SurfPtrD->MyNum); //printf("angle %f Areaangle %f \n",angle,SurfPtrD->AreaAngle); //printf("Oben/Unten: Element: %d\n",SurfPtrD->MyNum); *eclipsed = EclipsedByDisk; } } } if (SurfPtrD->AreaType == OUT_RECTANGLE || SurfPtrD->AreaType == IN_RECTANGLE ) { if (Store->MyNum == -1 ) // debug test: { printf("Seite: Element: %ld\n",SurfPtrD->MyNum); printf("Distance:[%f,%f,%f]\n",Distance[0],Distance[1],Distance[2]); printf("Determinanten: %f, %f; r: %f\n",det1, det2,r); printf("angle %f Areaangle %f \n",angle,SurfPtrD->AreaAngle); printf("n[%f, %f, %f]\n", SurfPtrD -> l, SurfPtrD->m,SurfPtrD->n); printf("vec_e1[%f, %f, %f]\n", vec_e1[0],vec_e1[1],vec_e1[2]); printf("vec_e2[%f, %f, %f]\n", vec_e2[0],vec_e2[1],vec_e2[2]); printf("vec_p1[%f, %f, %f]\n", vec_p1[0],vec_p1[1],vec_p1[2]); printf("vec_p2[%f, %f, %f]\n", vec_p2[0],vec_p2[1],vec_p2[2]); printf("vec_g1[%f, %f, %f]\n", vec_g1[0],vec_g1[1],vec_g1[2]); /* printf("delta:[%f,%f,%f]\n"); */ printf("------------------------------\n"); } /* Mittelpunkt von element vorhanden lambda,mu,nu */ if (Distance[2] < SurfPtrD->nu+SurfPtrD->AreaHeight/2 && Distance[2] > SurfPtrD->nu-SurfPtrD->AreaHeight/2) { if ((AreaAngle1>AreaAngle2 && angle < AreaAngle1 && angle > AreaAngle2)|| (AreaAngle1 AreaAngle2 || angle < AreaAngle1))) { *eclipsed = EclipsedByDisk; } else { } } } } /* if r >=0.0 */ } /* If det1 !=0 */ } ++SurfPtrD; Element_number++; } if (*eclipsed == 0) //printf("Num: %d \n",Store->MyNum); SurfPtrD=Store; } void diskHeight (double * Height_in, double * Height_out) { double diskExp = 1.0; /* Exponent for r dependency of scale height */ double diskProp = 1.0; /* Proportionality factor for disk height */ double diskAdd = 1.0; /* Additive constant for disk height */ double Rout = Binary[Disk].RLag1 * Binary[Disk].Rout; double Rin = Binary[Disk].RLag1 * Binary[Disk].Rin; if (Flags.tdisk == 0) { diskAdd = Binary[Disk].Thick; diskProp = Binary[Disk].HR; diskExp = 1.0; } else { if (Flags.tdisk == 1) /* isothermal */ diskExp = 1.5; else diskExp = (9.0/8.0); /* reprocessing */ diskAdd = 0.0; diskProp = (Binary[Disk].Thick / pow(Rin, diskExp)); } /* Height of the Disk on the outer side */ *Height_out = (diskAdd + diskProp * pow(Rout, diskExp)) / 2.0; /* Height of the Disk on the inner side */ *Height_in = (diskAdd + diskProp * pow(Rin, diskExp)) / 2.0; } /******************************************************************** @package nightfall @author Rainer Wichmann @version 1.55 @short Test eclipse of secondary by surrounding disk @param (double) CosPhase2 Cosine of phase @param (double) SinPhase2 Sine of phase @return (void) @heading Light Curve ********************************************************************/ void SecondaryEclipsedByDisk (double CosPhase2, double SinPhase2) { /* Radii of outer and inner cylinder between which * the disk is situated. */ double Rout = Binary[Disk].RLag1 * Binary[Disk].Rout; double Rin = Binary[Disk].RLag1 * Binary[Disk].Rin; /* Height of the Disk on the outer side */ double HeightIn; double HeightOut; double x_z; double y_z; double z_z; double A; double B; double A2B; SurfaceElement *SurfPtrS = Surface[Secondary]; unsigned long NElemS = Binary[Secondary].NumElem; unsigned long countS = 0; double Tani = tan(Orbit.Inclination); double Zin, Zout; double t_1; double t_2; /* View from above */ if (Tani < FLT_EPSILON) return; diskHeight (&HeightIn, &HeightOut); while (countS < NElemS) { /* Invisible anyway */ if (SurfPtrS->visibility == 0) goto next; /* Above disk => visible */ if (SurfPtrS->nu > HeightOut) goto next; x_z = SurfPtrS->lambda; y_z = (-SurfPtrS->mu); z_z = SurfPtrS->nu; /* Compute intersection with inner cylinder */ A = x_z * CosPhase2 + y_z * SinPhase2; B = (x_z * x_z) + (y_z * y_z) - (Rin * Rin); A2B = (A*A) - B; if (A2B <= 0) goto next; /* z-coordinate of intersection */ t_1 = ((-A) - sqrt(A2B)) / Tani; t_2 = ((-A) + sqrt(A2B)) / Tani; Zin = MAX(t_1,t_2) + z_z; /* Compute intersection with outer cylinder */ B = (x_z * x_z) + (y_z * y_z) - (Rout * Rout); A2B = (A*A) - B; if (A2B <= 0) goto next; /* z-coordinate of intersection */ t_1 = ((-A) - sqrt(A2B)) / Tani; t_2 = ((-A) + sqrt(A2B)) / Tani; Zout = MAX(t_1,t_2) + z_z; if ( (Zout < (-HeightOut)) || (Zin > HeightIn && Zout > HeightOut) ) goto next; SurfPtrS->EclipseFlag = 1; SurfPtrS->visibility = 0; next: ++SurfPtrS; ++countS; } return; } /******************************************************************** @package nightfall @author Rainer Wichmann @version 1.55 @short Test eclipse of disk by itself @param (double) CosPhase2 Cosine of phase @param (double) SinPhase2 Sine of phase @return (void) @heading Light Curve ********************************************************************/ void DiskEclipsedByDisk (double CosPhase2, double SinPhase2) { /* Radii of outer and inner cylinder between which * the disk is situated. */ double Rout = Binary[Disk].RLag1 * Binary[Disk].Rout; double Rin = Binary[Disk].RLag1 * Binary[Disk].Rin; /* Height of the Disk on the outer side */ double HeightIn; double HeightOut; double x_z; double y_z; double z_z; double A; double B; double A2B; SurfaceElement *SurfPtr = Surface[Disk]; unsigned long NElem = Binary[Disk].NumElem; unsigned long count = 0; double Tani = tan(Orbit.Inclination); double Zin, Zout; double t_1; double t_2; /* View from above */ if (Tani < FLT_EPSILON) return; diskHeight (&HeightIn, &HeightOut); while (count < NElem) { /* Invisible anyway */ if (SurfPtr->visibility == 0) goto next; /* Above disk => visible */ if (SurfPtr->nu > HeightOut) goto next; /* Outer side => never eclipsed by disk */ if (SurfPtr->AreaType == OUT_RECTANGLE) goto next; x_z = SurfPtr->lambda; y_z = (-SurfPtr->mu); z_z = SurfPtr->nu; A = x_z * CosPhase2 + y_z * SinPhase2; if (SurfPtr->AreaType == IN_RECTANGLE) { /* Compute intersection with inner cylinder */ B = (x_z * x_z) + (y_z * y_z) - (Rin * Rin); A2B = (A*A) - B; if (A2B <= 0) goto outer; /* z-coordinate of intersection */ t_1 = ((-A) - sqrt(A2B)) / Tani; t_2 = ((-A) + sqrt(A2B)) / Tani; Zin = MAX(t_1,t_2) + z_z; if (Zin < HeightIn && Zin > (-HeightIn)) goto set_eclipsed; } outer: /* Compute intersection with outer cylinder */ B = (x_z * x_z) + (y_z * y_z) - (Rout * Rout); A2B = (A*A) - B; if (A2B <= 0) goto next; /* z-coordinate of intersection */ t_1 = ((-A) - sqrt(A2B)) / Tani; t_2 = ((-A) + sqrt(A2B)) / Tani; Zout = MAX(t_1,t_2) + z_z; if (Zout < (-HeightOut) || Zout > (HeightOut)) goto next; set_eclipsed: SurfPtr->EclipseFlag = 1; SurfPtr->visibility = 0; next: ++SurfPtr; ++count; } return; } /******************************************************************** @package nightfall @author Rainer Wichmann @version 1.55 @short Test eclipse of primary by disk surrounding secondary @param (double) CosPhase2 Cosine of phase @param (double) SinPhase2 Sine of phase @return (void) @heading Light Curve ********************************************************************/ void PrimaryEclipsedByDisk (double CosPhase, double SinPhase) { /* Radii of outer and inner cylinder between which * the disk is situated. */ double Rout = Binary[Disk].RLag1 * Binary[Disk].Rout; double Rin = Binary[Disk].RLag1 * Binary[Disk].Rin; /* Height of the Disk on the outer side */ double HeightIn; double HeightOut; double x_z; double y_z; double z_z; double A; double B; double A2B; SurfaceElement *SurfPtr = Surface[Primary]; unsigned long NElem = Binary[Primary].NumElem; unsigned long count = 0; double Tani = tan(Orbit.Inclination); double Z; double t_1; double t_2; /* View from above */ if (Tani < FLT_EPSILON) return; diskHeight (&HeightIn, &HeightOut); while (count < NElem) { /* Invisible anyway */ if (SurfPtr->visibility == 0) goto next; /* Above disk => visible */ if (SurfPtr->nu > HeightOut) goto next; x_z = SurfPtr->lambda; y_z = SurfPtr->mu; z_z = SurfPtr->nu; /* Compute intersection with outer cylinder */ A = x_z * CosPhase + y_z * SinPhase - CosPhase; B = 1.0 + (x_z * x_z) + (y_z * y_z) - (2.0 * x_z) - (Rout * Rout); A2B = (A*A) - B; if (A2B <= 0) goto next; /* z-coordinate of intersection */ t_1 = ((-A) - sqrt(A2B)) / Tani; t_2 = ((-A) + sqrt(A2B)) / Tani; Z = MIN(t_1,t_2) + z_z; if (Z > HeightOut) goto next; Z = MAX(t_1,t_2) + z_z; if (Z < (-HeightOut)) goto next; /* Compute intersection with inner cylinder */ B = 1.0 + (x_z * x_z) + (y_z * y_z) - (2.0 * x_z) - (Rin * Rin); A2B = (A*A) - B; if (A2B > 0) { /* z-coordinate of intersection */ t_1 = ((-A) - sqrt(A2B)) / Tani; t_2 = ((-A) + sqrt(A2B)) / Tani; Z = MIN(t_1,t_2) + z_z; if (Z < (-HeightIn)) { Z = MAX(t_1,t_2) + z_z; if (Z > HeightIn) goto next; } } SurfPtr->EclipseFlag = 1; SurfPtr->visibility = 0; next: ++SurfPtr; ++count; } return; } #endif