/* 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. */
/* ascending node is at +pi/2 counted clockwise from LOS */
/* at secondary eclipse */
/* periastron length is counted clockwise from ascending node */
/* eg. periastron is at secondary eclipse for 270 deg length */
#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 Find the phase of periastron
@long Given Omega ( == length of periastron) and Eccentricity
(both in global Orbit structure), will determine the
phase of periastron passage (Orbit.OmegaZero), the
mean anomaly of the quadrature (Orbit.MQuad), and the
minimum distance Orbit.MinDist. Orbit.MQuad
Orbit.MQuad is the starting point for lightcurve
@param (void)
@return (void)
@heading Eccentric Orbit
****************************************************************************/
void FindPeriastron()
{
/* find the phase of the periastron */
/* we have t = 0, M = -t0, nu = 3pi/2-Omega */
/* at periastron */
double temp; /* temporary storage */
OrbitType *OrbPtr; /* orbit pointer */
OrbPtr = &Orbit;
OrbPtr->Nu = (3.0 * PI / 2.0) - (DTOR * OrbPtr->Omega);
/* excentric anomaly */
temp = sqrt( (1.0 + OrbPtr->Excentricity) / (1.0 - OrbPtr->Excentricity));
temp = tan(OrbPtr->Nu/2.0) / temp;
OrbPtr->E = 2.0 * atan(temp);
if (OrbPtr->E <= (-FLT_EPSILON)) OrbPtr->E = OrbPtr->E + 2.0 * PI;
/* passage of periastron OmegaZero = -(mean anomaly) */
OrbPtr->OmegaZero = -(OrbPtr->E - OrbPtr->Excentricity * sin(OrbPtr->E));
/* minimum distance */
OrbPtr->MinDist = (1.0 - SQR(OrbPtr->Excentricity))/( 1.0 + OrbPtr->Excentricity);
/* starting point */
OrbPtr->Nu = PI - (DTOR * OrbPtr->Omega);
temp = sqrt( (1.0 + OrbPtr->Excentricity) / (1.0 - OrbPtr->Excentricity));
temp = tan(OrbPtr->Nu/2.0) / temp;
OrbPtr->E = 2.0 * atan(temp);
/* mean anomaly of quadrature */
OrbPtr->MQuad = OrbPtr->E - OrbPtr->Excentricity * sin(OrbPtr->E);
if (Flags.debug[VERBOSE] == ON) {
printf("\n");
printf(_("Eccentric orbit: \n"));
printf(_(" Mean Anomaly at Secondary Eclipse %7.4f,\n"),OrbPtr->OmegaZero);
printf(_(" Minimum Distance %6.3f,\n"),OrbPtr->MinDist);
printf(_(" Mean Anomaly at Quadrature %7.4f\n"),OrbPtr->MQuad);
}
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Solve the Kepler equation
@long given j = index of phase, will determine phase and distance,
and update Orbit.Phase and Orbit.Dist accordingly
@param (int) j The index of the phase
@return (void)
@heading Eccentric Orbit
****************************************************************************/
void SolveKepler(int j)
{
double Old, New; /* for iterative solution */
OrbitType *OrbPtr; /* orbit pointer */
OrbPtr = &Orbit;
/* we have phase zero at secondary eclipse, and start at quadrature */
/* rotation symmetry about LOS, thus length of ascending node redundant */
/* OmegaQuad = MQuad + OmegaZero, M = OmegaQuad + d(Omega) - OmegaZero */
/* => M = MQuad + d(Omega) */
OrbPtr->M = ( OrbPtr->MQuad + (2.0*PI/PhaseSteps)*j );
/* solve Kepler */
/* start value */
New = OrbPtr->M + OrbPtr->Excentricity * sin(OrbPtr->M);
do {
Old = New;
New = OrbPtr->M + OrbPtr->Excentricity * sin(Old);
} while (fabs(Old-New) >= FLT_EPSILON);
OrbPtr->E = New;
OrbPtr->Nu = 2.0 * atan( tan(OrbPtr->E/2.0)
* ( sqrt ( (1.0 + OrbPtr->Excentricity)
/ (1.0 - OrbPtr->Excentricity))));
/* set actual phase angle - PI/2 for quadrature w/ primary left side */
OrbPtr->Phase = OrbPtr->Nu + (DTOR * OrbPtr->Omega) - (PI/2.0);
if (OrbPtr->Phase <= (-FLT_EPSILON)) OrbPtr->Phase = OrbPtr->Phase + 2.0*PI;
OrbPtr->Dist = (1.0 - SQR(OrbPtr->Excentricity) )
/ ( 1.0 + OrbPtr->Excentricity * cos (OrbPtr->E));
/* normalize distance to minimum, check for roundoff error */
OrbPtr->Dist = OrbPtr->Dist / OrbPtr->MinDist;
if ((OrbPtr->Dist - 1.0) <= (-FLT_EPSILON)) OrbPtr->Dist = 1.0;
return;
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Update radius and potential
@param (int) Component The stellar component
@return (int) The exit Status
@heading Eccentric Orbit
****************************************************************************/
int UpdateGeometry(int Comp)
{
double radphase; /* mean radius */
double radpol; /* polar radius */
double rochepot; /* Roch potential */
int testerr = 0; /* exit status */
OrbitType *OrbPtr = &Orbit; /* orbit pointer */
/* scale the volume */
Binary[Comp].ScaledVolume =
Binary[Comp].Volume / (OrbPtr->Dist * OrbPtr->Dist * OrbPtr->Dist);
/* next seach the mean radius as function of scaled volume */
/* set global variables */
PhaseScaledVolume = Binary[Comp].ScaledVolume;
Mq = Binary[Comp].Mq;
/* mean radius */
radphase = RootFind(VolOfR, FLT_EPSILON,
Binary[Comp].Radius+0.04, FLT_MIN,
"UpdateGeometry", &testerr);
if (testerr == 1) return(8);
/* then compute polar radius as function of mean radius */
radpol = PolRad(Binary[Comp].Mq, radphase);
Binary[Comp].Radius = radpol;
/* adjust fill factor, do sanity check */
Binary[Comp].RocheFill = Binary[Comp].Radius / Binary[Comp].RCrit;
if ((Binary[Comp].RocheFill - 1.0) >= FLT_EPSILON)
Binary[Comp].RocheFill = 1.0;
/* finally compute the roche potential */
rochepot = 1.0/radpol
+ Binary[Comp].Mq/sqrt(1 + SQR(radpol));
Binary[Comp].RochePot = rochepot;
return(0);
}
syntax highlighted by Code2HTML, v. 0.9.1