/* 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"
#ifdef _WITH_PGPLOT
#ifndef _WITH_GNUPLOT
#include "cpgplot.h"
#endif
#endif
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Compute fractional visibility
@param (int) Comp The stellar component
@return (void)
@heading Star Setup
****************************************************************************/
void LightFractionalPot(int Comp)
{
long NElem; /* # of surface elements */
long i, j; /* loop variables */
static long LastEclipsed[MAXELEMENTS]; /* array of indices to eclipsed */
static long NextEclipsed[MAXELEMENTS]; /* array of indices to eclipsed */
long Last; /* element eclipsed before */
long Next; /* element eclipsed after */
double lzero, mzero, nzero; /* surface normal vector */
int Comp2; /* the other star */
double LOS_PotLast; /* line-of-sight potential min */
double LOS_PotNext; /* line-of-sight potential min */
double LOS_PotInt; /* line-of-sight potential min */
double x_z, y_z, z_z; /* surface vector */
double t_1, t_2; /* brackets for minfind */
double tmin = 0; /* minimum along LOS */
double Phase; /* orbital phase */
int Acomp; /* sign */
double Qi, Li, QLi; /* smalles circle around star */
double t_l; /* temporary storage */
double SqrXp; /* misc computational */
double temp; /* temporary storage */
double fratio = 1.0; /* the async rotation */
j = 0; t_1 = 0.03; t_2 = 1.0;
/* secondary is mirrored - we turn by PI and mirror y->(-y) */
/* to have secondary at (0,0,0) with correct orientation */
/* we are interested in the potential surface of the OTHER star */
if (Comp == Primary) {
Phase = Orbit.Phase;
Comp2 = 1; Acomp = 1;
SqrXp = SQR(Binary[Secondary].Xp);
/* SqrPol = SQR(Binary[Secondary].Radius); */
if (Flags.asynchron2 == ON) fratio = Binary[Secondary].Fratio;
} else {
Phase = Orbit.Phase + PI;
Comp2 = 0; Acomp = (-1);
SqrXp = SQR(Binary[Primary].Xp);
/* SqrPol = SQR(Binary[Primary].Radius); */
if (Flags.asynchron1 == ON) fratio = Binary[Primary].Fratio;
}
/* ------- vector towards observer (LOS = Line Of Sight) --------- */
lzero = sin(Orbit.Inclination) * cos(Phase);
mzero = sin(Orbit.Inclination) * sin(Phase);
nzero = cos(Orbit.Inclination);
/* ----------- search for eclipsed/non-eclipsed pairs ----------- */
NElem = Binary[Comp].NumElem;
for (i = 0; i < NElem; ++i) {
LastEclipsed[i] = 0;
NextEclipsed[i] = 0;
}
j = 0;
for (i = Binary[Comp].N_PhiStep[0]; i < NElem; ++i) {
if (fabs((double)Surface[Comp][i].EclipseFlag -
(double)Surface[Comp][i].next_ptr->EclipseFlag - 1.0)
<= FLT_EPSILON) {
NextEclipsed[j] = Surface[Comp][i].next_ptr->MyNum;
if (NextEclipsed[j] > (Binary[Comp].NumElem - 1))
NextEclipsed[j] = Binary[Comp].NumElem - 1;
LastEclipsed[j] = i;
++j;
}
if ( fabs((double)Surface[Comp][i].EclipseFlag -
(double)Surface[Comp][i].last_ptr->EclipseFlag - 1.0)
<= FLT_EPSILON) {
NextEclipsed[j] = Surface[Comp][i].last_ptr->MyNum;
if (NextEclipsed[j] > (Binary[Comp].NumElem - 1))
NextEclipsed[j] = Binary[Comp].NumElem - 1;
LastEclipsed[j] = i;
++j;
}
if (fabs((double)Surface[Comp][i].EclipseFlag -
(double) Surface[Comp][i].phi_ptr->EclipseFlag - 1.0)
<= FLT_EPSILON) {
NextEclipsed[j] = Surface[Comp][i].phi_ptr->MyNum;
if (NextEclipsed[j] > (Binary[Comp].NumElem - 1))
NextEclipsed[j] = Binary[Comp].NumElem - 1;
LastEclipsed[j] = i;
++j;
}
if (fabs((double)Surface[Comp][i].EclipseFlag -
(double)Surface[Comp][i].phi1_ptr->EclipseFlag - 1.0)
<= FLT_EPSILON) {
NextEclipsed[j] = Surface[Comp][i].phi1_ptr->MyNum;
if (NextEclipsed[j] > (Binary[Comp].NumElem - 1))
NextEclipsed[j] = Binary[Comp].NumElem - 1;
LastEclipsed[j] = i;
++j;
}
}
/* ---------- compute fractional visibility ------------------- */
for (i = 0; i < j; ++i) {
Last = LastEclipsed[i]; Next = NextEclipsed[i];
if (fabs(Surface[Comp][Last].LOS_Pot) >= FLT_EPSILON) {
LOS_PotLast = Surface[Comp][Last].LOS_Pot;
} else {
x_z = Surface[Comp][Last].lambda;
y_z = Acomp * (Surface[Comp][Last].mu);
z_z = Surface[Comp][Last].nu;
Qi = 2.0*(x_z*lzero+y_z*mzero+z_z*nzero-lzero);
Li = 1.0 + (x_z*x_z)+(y_z*y_z)+(z_z*z_z)-2*x_z-SqrXp;
QLi = (Qi*Qi)-4.0*Li;
if (QLi >= FLT_EPSILON) { /* intersection */
t_l = sqrt(QLi); /* just to hold temporary result */
t_1 = (-Qi + t_l)/2.0;
t_2 = (-Qi - t_l)/2.0;
} else {
t_1 = 0.2;
t_2 = 1.0;
}
/* changed Binary[Comp].Mq to Binary[Comp2].Mq Jul 6 */
(void) MinFinder(&Binary[Comp2].Mq, &fratio,
&t_1, &t_2,
&lzero, &mzero, &nzero,
&x_z, &y_z, &z_z,
&tmin, &LOS_PotLast);
LOS_PotLast = -LOS_PotLast;
Surface[Comp][Last].LOS_Pot = (float) LOS_PotLast;
}
if (fabs(Surface[Comp][Next].LOS_Pot) >= FLT_EPSILON) {
LOS_PotNext = Surface[Comp][Next].LOS_Pot;
} else {
x_z = Surface[Comp][Next].lambda;
y_z = Acomp * (Surface[Comp][Next].mu);
z_z = Surface[Comp][Next].nu;
Qi = 2.0*(x_z*lzero+y_z*mzero+z_z*nzero-lzero);
Li = 1.0 + (x_z*x_z)+(y_z*y_z)+(z_z*z_z)-2*x_z-SqrXp;
QLi = Qi*Qi-4.0*Li;
if (QLi >= FLT_EPSILON) { /* intersection */
t_l = sqrt(QLi); /* just to hold temporary result */
t_1 = (-Qi + t_l)/2.0;
t_2 = (-Qi - t_l)/2.0;
} else {
t_1 = 0.2;
t_2 = 1.0;
}
/* changed Binary[Comp].Mq to Binary[Comp2].Mq Jul 6 */
(void) MinFinder(&Binary[Comp2].Mq, &fratio,
&t_1, &t_2,
&lzero, &mzero, &nzero,
&x_z, &y_z, &z_z,
&tmin, &LOS_PotNext);
LOS_PotNext = -LOS_PotNext;
Surface[Comp][Next].LOS_Pot = (float) LOS_PotNext;
}
LOS_PotInt = fabs((LOS_PotNext - LOS_PotLast)/2.);
/* EclipseFlag == 1 means IsEclipsed */
temp = fabs(LOS_PotNext - Binary[Comp2].RochePot);
if ((temp+FLT_EPSILON) <= fabs(LOS_PotInt)) {
if (Surface[Comp][Next].EclipseFlag == 1)
Surface[Comp][Next].visibility =
(float) (0.0 + 0.5 * (1.0 - temp/ LOS_PotInt ));
if (Surface[Comp][Next].EclipseFlag == 0)
Surface[Comp][Next].visibility =
(float) (1.0 - 0.5 * (1.0 - temp/ LOS_PotInt ));
} else {
Surface[Comp][Next].visibility =
(float) fabs(Surface[Comp][Next].EclipseFlag-1.0);
}
temp = fabs(LOS_PotLast - Binary[Comp2].RochePot);
if ((temp+FLT_EPSILON) <= fabs(LOS_PotInt)) {
if (Surface[Comp][Last].EclipseFlag == 1)
Surface[Comp][Last].visibility =
(float) (0.0 + 0.5 * (1.0 - temp/ LOS_PotInt ));
if (Surface[Comp][Last].EclipseFlag == 0)
Surface[Comp][Last].visibility =
(float) (1.0 - 0.5 * (1.0 - temp/ LOS_PotInt ));
} else {
Surface[Comp][Last].visibility =
(float) fabs(Surface[Comp][Last].EclipseFlag-1.0);
}
}
}
syntax highlighted by Code2HTML, v. 0.9.1