/* 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 <math.h>
#include <stdio.h>
#include <string.h>
#include <float.h>
#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 > 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
syntax highlighted by Code2HTML, v. 0.9.1