/* 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 <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <ctype.h>
#include <float.h>
#include "Light.h"
/* contains miscellaneous small routines */
/******************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Get data path
@param (void)
@return (char) * The path
@heading GUI
*******************************************************************/
char * data_data_fls()
{
static char fls_path[1024];
int i;
for (i=0; i<1024; ++i) fls_path[i] = '\0';
if (getenv("NIGHTFALL_DATA_DIR") != NULL) {
strncpy(fls_path, getenv("NIGHTFALL_DATA_DIR"), 1023);
return fls_path;
}
if (getenv("NIGHTFALL_DATAROOT") != NULL) {
strncpy(fls_path, getenv("NIGHTFALL_DATAROOT"), (1023 - 5) );
strcat(fls_path, "/data");
return fls_path;
}
strncpy(fls_path, DATAFLS, (1023 - 5) );
return fls_path;
}
/******************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Get cfg path
@param (void)
@return (char) * The path
@heading GUI
*******************************************************************/
char * data_cfg_fls()
{
static char fls_path[1024];
int i;
for (i=0; i<1024; ++i) fls_path[i] = '\0';
if (getenv("NIGHTFALL_CFG_DIR") != NULL) {
strncpy(fls_path, getenv("NIGHTFALL_CFG_DIR"), 1023);
return fls_path;
}
if (getenv("NIGHTFALL_DATAROOT") != NULL) {
strncpy(fls_path, getenv("NIGHTFALL_DATAROOT"), (1023 - 5) );
strcat(fls_path, "/cfg");
return fls_path;
}
strncpy(fls_path, CFGFLS, (1023 - 5) );
return fls_path;
}
/******************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Get doc path
@param (void)
@return (char) * The path
@heading GUI
*******************************************************************/
char * data_doc_fls()
{
static char fls_path[1024];
int i;
for (i=0; i<1024; ++i) fls_path[i] = '\0';
if (getenv("NIGHTFALL_DOC_DIR") != NULL) {
strncpy(fls_path, getenv("NIGHTFALL_DOC_DIR"), 1023);
return fls_path;
}
if (getenv("NIGHTFALL_DATAROOT") != NULL) {
strncpy(fls_path, getenv("NIGHTFALL_DATAROOT"), (1023 - 5) );
strcat(fls_path, "/doc");
return fls_path;
}
strncpy(fls_path, DOCFLS, (1023 - 5) );
return fls_path;
}
/******************************************************************
@package nightfall
@author Markus Kuster (kuster@astro.uni-tuebingen.de)
@version 1.0
@short Get pixmap path
@param (void)
@return (char) * The path
@heading GUI
*******************************************************************/
char * data_pix_fls()
{
static char fls_path[1024];
int i;
for (i=0; i<1024; ++i) fls_path[i] = '\0';
if (getenv("NIGHTFALL_PIX_DIR") != NULL) {
strncpy(fls_path, getenv("NIGHTFALL_PIX_DIR"), 1023);
return fls_path;
}
if (getenv("NIGHTFALL_DATAROOT") != NULL) {
strncpy(fls_path, getenv("NIGHTFALL_DATAROOT"), (1023 - 5) );
strcat(fls_path, "/pixmaps");
return fls_path;
}
strncpy(fls_path, PIXFLS, (1023 - 5) );
return fls_path;
}
static int gd_size = 32;
static float gd_tt[] = {
2012.6, 2322.9, 2646.3, 2824.5, 3054.2,
3238.7, 3302.6, 3434.3, 3548.1, 3665.6,
3938, 4230.6, 4634.8, 5077.5, 5562.5,
5709.5, 6093.9, 6254.9, 6504.3, 6676,
7033.3, 7125.6, 7313.8, 7506.9, 7605.4,
7705.2, 7806.3, 7908.7, 8012.4, 8117.5,
8835.2, 10467, };
static float gd_beta[] = {
0.22016, 0.22326, 0.22016, 0.2031, 0.1938,
0.1845, 0.17054, 0.18295, 0.1845, 0.1938,
0.23256, 0.34419, 0.4155, 0.43411, 0.4031,
0.3876, 0.34574, 0.32248, 0.29922, 0.27287,
0.24806, 0.2031, 0.15814, 0.17364, 0.31318,
0.58605, 0.72868, 0.8155, 0.93643, 0.94574,
0.9876, 1, };
/******************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Compute gravitational darkening exponent
@param (float) Temperature
@return (float) Beta, the gravitational darkening exponent
@heading Light Curve
*******************************************************************/
float getGravDarkExp (float temp)
{
int i;
float val = 1.0;
gd_beta[gd_size-1] = 1.0;
if (temp <= gd_tt[0])
return gd_beta[0];
if (temp >= gd_tt[gd_size-1])
return 1.0;
for (i = 1; i < gd_size; ++i)
{
if (fabs(gd_tt[i] - temp) < FLT_EPSILON)
return gd_beta[i];
if (gd_tt[i] > temp)
{
val = gd_beta[i-1] * ((gd_tt[i]-temp)/(gd_tt[i] - gd_tt[i-1]))
+ gd_beta[i] * ((temp-gd_tt[i-1])/(gd_tt[i] - gd_tt[i-1]));
return val;
}
}
return val;
}
static float Albedo[NUM_COMP] = { 0.0 };
#define ALBEDO_OFFSET 1.0
void SetAlbedo (int Comp, const char * str)
{
float ff = atof(str);
if (Comp >=0 && Comp < NUM_COMP && 0.0 <= ff && ff <= 1.0)
Albedo[Comp] = (ff + ALBEDO_OFFSET);
else
WARNERR(_("Bad albedo value: "), str);
return;
}
void ComputeGravDark()
{
#ifdef HAVE_DISK
int ncomp = 3;
#else
int ncomp = 2;
#endif
int i;
double limit = 7700.0;
if (getenv("NIGHTFALL_RADIATIVE") != NULL)
limit = atof(getenv("NIGHTFALL_RADIATIVE"));
if (limit <= 4000. || limit >= 12000.) limit = 7700.0;
/* >>>>>>>>>>> determine albedo and GravDarkCoeff <<<<<<<< */
for (i = 0; i < ncomp; ++i) {
if (Binary[i].Temperature <= limit) {
/* -------------- convective ---------------------- */
if (Albedo[i] >= ALBEDO_OFFSET)
Binary[i].Albedo = Albedo[i] - ALBEDO_OFFSET;
else
Binary[i].Albedo = 0.5;
Binary[i].GravDarkCoeff = getGravDarkExp(Binary[i].Temperature);
Binary[i].GravDarkCoeff = Binary[i].GravDarkCoeff / 4.0;
} else {
/* ---------------- radiative ------------------------ */
if (Albedo[i] >= ALBEDO_OFFSET)
Binary[i].Albedo = Albedo[i] - ALBEDO_OFFSET;
else
Binary[i].Albedo = 1.0;
Binary[i].GravDarkCoeff = 0.25;
}
if (Flags.debug[VERBOSE] == ON)
printf(_(" Comp %2d GravDarkCoeff %8.4f Albedo %8.4f\n"),
i+1, Binary[i].GravDarkCoeff, Binary[i].Albedo);
}
return;
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Set Xwindow size for PGPLOT
@tip Allocates 21 * sizeof(char) without freeing (called only once)
@param (void)
@return (void)
@heading Plotting
****************************************************************************/
void SetPgplotWindow()
{
int width = 551; /* width and height */
char *windowsize; /* the window size */
float frac; /* fractional display width */
char * frac_size; /* fractional display width */
static char hasbuf[32];
char gnu_geometry[256];
/* this function is only called at startup, we never free frac_size */
/* unless program exits (not really memory leak) */
/*@i@*/ frac_size = malloc(21);
if (frac_size == NULL) nf_error(_(errmsg[0]));
/* check whether environment variable is defined */
windowsize = getenv("PGPLOT_XW_WIDTH");
/* if not, set fractional display size */
if (windowsize == NULL) {
if (getenv("GNUPLOT_GEOMETRY") == NULL)
strncpy(gnu_geometry, GNU_GEOMETRY, 255);
else
strncpy(gnu_geometry, getenv("GNUPLOT_GEOMETRY"), 255);
gnu_geometry[255] = '\0';
(void) sscanf(gnu_geometry, "%dx", &width);
frac = (float)width/867.0; /* 867 pix is builtin default for PGPLOT */
frac = frac / 1.2;
sprintf(frac_size,"PGPLOT_XW_WIDTH=%4.2f", frac);
#ifdef HAVE_PUTENV
/* is putenv() POSIX ? */
/*@i@*/ putenv(frac_size);
#else
free(frac_size);
#endif
}
if (NULL == getenv("PGPLOT_BUFFER")) {
sprintf(hasbuf,"PGPLOT_BUFFER=yes");
putenv(hasbuf);
}
/*@i@*/ return;
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Print error message and exit
@param (char) error[] The error text
@return (void)
@heading Error handling (non-interactive)
****************************************************************************/
void nf_error(char *error_msg)
{
int i = 0, status = 0;
char command[256+64];
char command_one[256] = " \n";
while (i < 255 && *error_msg != '\0') {
command_one[i] = *error_msg;
++error_msg;
++i;
}
if (myrank == 0)
{
status =
system("which xmessage > /dev/null");
if (status == 0) {
sprintf(command,
_("xmessage \"ERROR: %s\" "), command_one);
status =
system(command);
}
}
fprintf(stderr, _("**ERROR**: %s\n"), command_one);
fprintf(stderr, _(".. exit to system\n"));
MPI_Abort(MPI_COMM_WORLD, EXIT_FAILURE);
exit (EXIT_FAILURE);
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Convert string to lowercase
@tip strlwr is not ANSI
@param (char) *instring The input string
@return (void)
@heading String Handling
****************************************************************************/
void nf_strlwr(char * instring)
{
if (instring == NULL) return;
while ( (*instring) != '\0') {
if (isalpha(*instring) != 0) (*instring) = tolower(*instring);
++instring;
}
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Allocate a float 2D matrix
@tip Code snippet from comp.lang.c FAQ, modified
@param (long) nrh number of rows
@param (long) nch number of collumns
@return (float **) pointer to allocated memory
@heading Memory Handling
****************************************************************************/
float **matrix(long nrows, long ncols)
{
long i;
float **m;
/* allocate pointers to rows */
/*@ignore@*/
/* cast argument of malloc b/o NEC SX-5 compiler warning */
m = malloc((size_t) (nrows) * sizeof(float*));
if (m == NULL) return (NULL);
/* allocate rows and set pointers to them */
for (i = 0; i < nrows; ++i) {
m[i] = malloc((size_t)(ncols)*sizeof(float));
if (!m[i]) return (NULL);
}
/* return pointer to array of pointers to rows */
return m;
/*@end@*/
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Set normalization phase for output flux (from data file)
@param (int) band the bandpass
@param (int) phase the phase at which to normalize
@return (void)
@heading Light Curve
****************************************************************************/
static float NormalPhase[NUM_MAG] = {
-1.0, -1.0, -1.0, -1.0,
-1.0, -1.0, -1.0, -1.0,
-1.0, -1.0, -1.0, -1.0
};
double P90[NUM_MAG] = {
0.75, 0.75, 0.75, 0.75,
0.75, 0.75, 0.75, 0.75,
0.75, 0.75, 0.75, 0.75
};
void SetNormalPhase (int band, float phase)
{
/* internal consistency check
*/
if (band >= NUM_MAG || band < 0)
{
if (Flags.debug[VERBOSE] == ON || Flags.debug[WARN] == ON)
fprintf(stderr,
_("**Warning**: NormalPhase(): Bandpass (%d) out of range (0, %d)\n"),
band, NUM_MAG-1);
return;
}
/* map phase to 0..1
*/
if (phase > 1.0)
phase = phase - ((int) phase);
else if (phase < -1.0)
phase = phase - ((int) phase);
phase = (phase < 0.0) ? (phase + 1.0) : phase;
NormalPhase[band] = phase;
return;
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Reset normalization phase for output flux
@param (void)
@return (void)
@heading Light Curve
****************************************************************************/
void ResetNormalPhase ()
{
register int i;
for (i = 0; i < NUM_MAG; ++i)
NormalPhase[i] = -1.0;
return;
}
void ResetNormalPhaseOne (int band)
{
if (band >= 0 && band < NUM_MAG)
NormalPhase[band] = -1.0;
return;
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Set normalization points for output flux
@param (int) band the bandpass
@param (int) phase the phase at which to normalize
@return (void)
@heading Light Curve
****************************************************************************/
static int NormalPoint[NUM_MAG] = { 0 };
void NormalPoints ()
{
register int j, i;
float mindist, dist;
float phase;
for (i = 0; i < NUM_MAG; ++i)
{
/* NormalPhase is the user-supplied (in data file) phase
* at which normalization should occur. Ideally this should be
* at quadrature.
*/
if (NormalPhase[i] < 0.0)
{
NormalPoint[i] = 0;
phase = (FluxOut[0].Phase / (PI+PI)) - 0.5;
/* map phase to 0..1
*/
if (phase > 1.0)
phase = phase - ((int) phase);
else if (phase < -1.0)
phase = phase - ((int) phase);
phase = (phase < 0.0) ? (phase + 1.0) : phase;
P90[i] = phase;
}
else
{
mindist = 1.0;
for (j = 0; j < PhaseSteps; ++j)
{
phase = (FluxOut[j].Phase / (PI+PI)) - 0.5;
/* map phase to 0..1
*/
if (phase > 1.0)
phase = phase - ((int) phase);
else if (phase < -1.0)
phase = phase - ((int) phase);
phase = (phase < 0.0) ? (phase + 1.0) : phase;
/* search for better minimum
*/
if (mindist > (dist = fabs(phase - NormalPhase[i])))
{
mindist = dist;
NormalPoint[i] = j;
P90[i] = phase;
}
}
}
}
return;
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.37
@short Luminosity for components
@param (int) component the bandpass
@param (int) phaseindex the current phase index
@return (void)
@heading Light Curve
****************************************************************************/
void LightLuminosity(int component, int phaseindex)
{
static int check = -1;
int band;
/* check whether all normalization points are -1 */
if (check == -1)
{
check = 0;
for (band = 0; band < NUM_MAG; ++band)
{
if ((NormalPhase[band] != -1))
{
check = 1;
break;
}
}
}
if ((check == 0) && (phaseindex == 0))
{
/* Normalization point is phase index 0 for all bands
*/
for (band = 0; band < NUM_MAG; ++band)
{
L90[component][band] = FluxOut[0].Mag[band];
if (component == Secondary)
L90[component][band] -= L90[Primary][band];
#ifdef HAVE_DISK
if (component == Disk)
L90[component][band] -= (L90[Primary][band]+L90[Secondary][band]);
#endif
}
}
if (check == 1)
{
float mindist = 1.0, dist;
float phase;
int j, jnorm = -1;
for (band = 0; band < NUM_MAG; ++band)
{
/* find best normalization phase jnorm
*/
if (NormalPhase[band] != -1.0)
{
for (j = 0; j <= phaseindex; ++j)
{
phase = (FluxOut[j].Phase / (PI+PI)) - 0.5;
/* map phase to 0..1
*/
if (phase > 1.0)
phase = phase - ((int) phase);
else if (phase < -1.0)
phase = phase - ((int) phase);
phase = (phase < 0.0) ? (phase + 1.0) : phase;
/* search for better minimum
*/
if (mindist > (dist = fabs(phase - NormalPhase[band])))
{
mindist = dist;
jnorm = j;
}
}
}
else
{
jnorm = 0;
}
/* copy flux to L90 if we are at best phase
*/
if (jnorm == phaseindex)
{
L90[component][band] = FluxOut[jnorm].Mag[band];
if (component == Secondary)
L90[component][band] -=
L90[Primary][band];
#ifdef HAVE_DISK
if (component == Disk)
L90[component][band] -=
(L90[Primary][band]+L90[Secondary][band]);
#endif
}
}
}
return;
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Normalize output flux
@param (void)
@return (void)
@heading Light Curve
****************************************************************************/
void NormalizeFlux()
{
int j, band;
double ThirdAdd[NUM_MAG];
/* ---------------- initialize -------------------------------------- */
NormalPoints ();
for (band = 0; band < NUM_MAG; ++band)
{
F90[band] = FluxOut[NormalPoint[band]].Mag[band];
if ( ( 1.0 - Orbit.Third[band]) >= FLT_EPSILON )
{
ThirdAdd[band] = F90[band]
* ( Orbit.Third[band] / ( 1.0 - Orbit.Third[band]) );
F90[band] = F90[band] + ThirdAdd[band];
}
}
/* ---------------- Normalize Output ------------------------------------ */
for (j = 0; j < PhaseSteps; ++j)
{
for (band = 0; band < NUM_MAG; ++band)
{
FluxOut[j].Mag[band] = -2.5
* log10((ThirdAdd[band] + FluxOut[j].Mag[band])/F90[band]);
}
}
return;
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Shift phase to correct interval
@param (void)
@return (void)
@heading Light Curve
****************************************************************************/
void PhaseShift()
{
int i, test;
test = ON;
do {
if ( (FluxOut[0].Phase - PI) >= FLT_EPSILON ) {
for (i = 0; i < PhaseSteps; ++i) {
FluxOut[i].Phase=FluxOut[i].Phase-2.0*PI;
}
test = ON; /* there was something to do */
} else {
if ( (FluxOut[0].Phase - (-PI)) <= (-FLT_EPSILON) ) {
for (i = 0; i < PhaseSteps; ++i) {
FluxOut[i].Phase=FluxOut[i].Phase+2.0*PI;
}
test = ON; /* there was something to do */
} else { test = OFF; }
}
} while (test == ON);
return;
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Read a single line from an open file
@param (char*) s The buffer string
@param (int) lim Size of buffer string
@param (FILE*) fpe Input filehandle
@return (int) Length of line read
@heading Light Curve
****************************************************************************/
int LireLigne(char *s, int lim, FILE *fpe)
{
if (fgets(s, lim, fpe) == NULL){
return 0;
}
return (int) strlen(s);
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Reset the data file list
@param (void)
@return (void)
@heading Data Fitting
****************************************************************************/
void ClearList()
{
FileType *item_ptr;
FileType *current;
item_ptr = FileList;
current = FileList;
while (current != NULL) {
item_ptr = current->nextFile;
free(current);
current = item_ptr;
}
FileList = NULL;
/* Reset all normalization points
*/
ResetNormalPhase ();
return;
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Add an item to the data file list
@param (char*) s The item
@return (void)
@heading Data Fitting
****************************************************************************/
void AddToList(const char * item, int Format)
{
FileType *new_ptr;
/*@ignore@*/
new_ptr = malloc(sizeof(FileType));
if (new_ptr != NULL) {
strncpy(new_ptr->DataFile, item, MAX_CFG_INLINE);
new_ptr->DataFormat = Format;
new_ptr->nextFile = FileList;
FileList = new_ptr;
}
return;
/*@end@*/
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Potential in Y direction at LagrangeOne
@tip Uses Global Variable Mq (Mass Ratio)
@param (double) y Distance in y-direction
@return (double)
@heading Light Curve
****************************************************************************/
double RocheYPerpendL1(double y)
{
/* Uses Global Variable Mq (Mass Ratio) */
/* Uses Global Variable RochePot */
/* Potential in Y direction at LagrangeOne */
/* used to define the throat */
double Eval; /* return value */
double r2, r1; /* radii from Prim./Second. */
/* pointer */
BinaryComponent *Bptr = &Binary[Primary];
r1 = sqrt(SQR(Bptr->RLag1) + (y*y));
r2 = sqrt(SQR(1.0-Bptr->RLag1) + (y*y));
Eval = -RochePot
+ (1.0/r1)
+ (Mq * (1.0/r2 - Bptr->RLag1))
+ ((Mq+1.0)/2.0) * (SQR(Bptr->RLag1) + (y*y));
return(Eval);
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Potential in Z direction at LagrangeOne
@tip Uses Global Variable Mq (Mass Ratio)
@param (double) z Distance in z-direction
@return (double)
@heading Light Curve
****************************************************************************/
double RocheZPerpendL1(double z)
{
/* Uses Global Variable Mq (Mass Ratio) */
/* Uses Global Variable RochePot */
/* Potential in Z direction at LagrangeOne */
/* used to define the throat */
double Eval; /* return value */
double r2, r1; /* radii from Prim./Second. */
/* pointer */
BinaryComponent *Bptr = &Binary[Primary];
r1 = sqrt( (Bptr->RLag1*Bptr->RLag1) + (z*z));
r2 = sqrt( (1.0-Bptr->RLag1)*(1.0-Bptr->RLag1) + (z*z));
Eval = -RochePot
+ (1.0/r1)
+ (Mq * (1.0/r2 - Bptr->RLag1))
+ ((Mq+1.0)/2.0) * (Bptr->RLag1 * Bptr->RLag1 ) ;
return(Eval);
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Potential in Z direction at Secondary
@tip Uses Global Variable Mq (Mass Ratio)
@param (double) z Distance in z-direction
@return (double)
@heading Light Curve
****************************************************************************/
double RochePerpendSecond(double z)
{
/* Uses Global Variable Mq (Mass Ratio) */
/* Uses Global Variable RochePot */
double Eval; /* return value */
double r2, r1; /* radii from Prim./Second. */
r1 = sqrt(1.0 + (z*z));
r2 = z;
if (z > DBL_EPSILON)
Eval = -RochePot
+ (1.0/r1)
+ (Mq * (1.0/r2 - 1.0))
+ ((Mq+1.0)/2.0);
else
{
if (Flags.debug[VERBOSE] == ON || Flags.debug[WARN] == ON)
fprintf(stderr, "**Warning**: Divide by zero, File %s, Line %d \n",
__FILE__, __LINE__);
Eval = -RochePot
+ (1.0/r1)
+ ((Mq+1.0)/2.0);
}
return(Eval);
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Potential
@tip Uses Global Variables Mq, RochePot, lambda, nu, F
@param (double) x Radius from origin
@return (double)
@heading Light Curve
****************************************************************************/
double RocheSurface(double x)
{
/* Uses Global Variable Mq (Mass Ratio) */
/* Uses Global Variable RochePot */
/* Uses Global Variable lambda (Coordinate) */
/* Uses Global Variable nu (Coordinate) */
/* Uses Global Variable F (nonsync ) */
/* input x is radius from origin */
double Eval; /* return value */
double sqrX; /* temp variable */
double lamX; /* temp variable */
sqrX = x * x;
lamX = x*lambda;
Eval = -RochePot
+ (1.0/x)
+ (Mq * (1.0/sqrt(1.0 + sqrX - 2*lamX) - lamX))
+ F*F* ((Mq+1.0)/2.0) * (sqrX * (1.0 - (nu * nu)));
return(Eval);
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Potential in x-y plane
@param (double) x x-vector
@param (double) y y-vector
@param (double) q mass ratio
@param (double) f async rotation
@return (double)
@heading Light Curve
****************************************************************************/
double RocheXYPrimary(double x, double y, double q, double f)
{
double Eval; /* return value */
double r1, r2; /* radii from Prim./Second. */
double sqrX; /* temp variable */
sqrX = x * x;
r1 = sqrt( sqrX +(y*y) );
r2 = sqrt( SQR(x-1.0) + (y*y) );
if (r1 > DBL_EPSILON && r2 > DBL_EPSILON)
Eval = 1.0/r1
+ q * ( 1.0/r2 - x )
+ ( (q+1.0)/2.0 ) * ( sqrX + (y*y) ) * (f*f);
else
{
if (Flags.debug[VERBOSE] == ON || Flags.debug[WARN] == ON)
fprintf(stderr, "**Warning**: Divide by zero, File %s, Line %d \n",
__FILE__, __LINE__);
Eval = 1.0;
}
return(Eval);
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Analytic polar radius
@param (double) q Mass ratio
@param (double) r Mean radius
@return (double)
@heading Light Curve
****************************************************************************/
double PolRad(double q, double r)
{
/* q is mass ratio, r mean radius */
/* these are the legendre functions for argument zero */
static double Pl2 = -1.0/2.0;
static double Pl4 = 3.0/8.0;
static double Pl6 = -5.0/16.0;
static double Pl8 = 35.0/128.0;
static double Pl10 = -315.0/1280.0;
register double powr; /* power of radius */
register double re; /* return value */
powr = r*r*r; /* r^3 */
re = powr*q*Pl2;
powr = powr*r*r; /* r^5 */
re = re + powr*q*Pl4;
powr = powr*r; /* r^6 */
re = re + powr*3.0*q*Pl2*q*Pl2;
powr = powr*r; /* r^7 */
re = re + powr*q*Pl6;
powr = powr*r; /* r^8 */
re = re + powr*8.0*q*q*Pl2*Pl4;
powr = powr*r; /* r^9 */
re = re + powr*q*Pl8 + powr*12.0*(q*Pl2)*(q*Pl2)*(q*Pl2);
powr = powr*r; /* r^10*/
re = re + powr*10.0*q*q*Pl2*Pl6 + powr*5.0*q*q*Pl4*Pl4;
powr = powr*r; /* r^11*/
re = re + powr*q*Pl10 + powr*55.0*q*q*q*Pl2*Pl2*Pl4;
re = r*re +r;
return(re);
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Returns (volume-scaledvolume) as function of r_0
@tip Zero if r_0 = RadMean(ScaledVolume)
@param (double) r_0 Radius
@return (double)
@heading Light Curve
****************************************************************************/
double VolOfR(double r_0)
{
/* uses global variable PhaseScaledVolume */
/* uses global variable Mq */
/* Volume analytic as by Kopal exact to O(r**11) */
register double NN; /* some temp variable */
register double RadPower; /* power of radius */
register double Vol1; /* volume */
double Volume; /* return value */
NN = (Mq + 1.0) / 2.0;
RadPower = r_0 * SQR(r_0); /* r^3 */
Vol1 = 1.0 + 2.0 * NN * RadPower;
RadPower = SQR(RadPower); /* r^6 */
Vol1 = Vol1 + (12.0/5.0)*(Mq*Mq)*RadPower
+ (32.0/5.0)*(NN*NN)*RadPower
+ (8.0/5.0)*NN*Mq*RadPower;
RadPower = RadPower*SQR(r_0); /* r^8 */
Vol1 = Vol1 + (15.0/7.0)*(Mq*Mq)*RadPower;
RadPower = RadPower*r_0; /* r^9 */
Vol1 = Vol1 + (22.0/7.0)*(Mq*Mq)*Mq*RadPower
+ (176.0/7.0)*(NN*NN)*NN*RadPower
+ (296.0/35.0)*NN*Mq
*(2.0*Mq + NN)*RadPower;
RadPower = RadPower*r_0; /* r^10*/
Vol1 = Vol1 + (18.0/9.0)*(Mq*Mq)*RadPower;
RadPower = RadPower*r_0; /* r^11*/
Vol1 = Vol1 + (157.0/7.0)*(Mq*Mq)*Mq*RadPower
+ (26.0/35.0)*NN*Mq
*(Mq + 3.0*NN)*RadPower;
Volume = Vol1
* (4.0*PI/3.0)*(r_0*r_0*r_0);
Volume = Volume - PhaseScaledVolume;
return(Volume);
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Zero at LagrangeOne
@param (double) x Position on x axix
@return (double)
@heading Light Curve
****************************************************************************/
double LagrangeOne(double x)
{
/* Uses Global Variable Mq (Mass Ratio) */
/* Eval is Zero at LagrangeOne (F=1) */
/* or at critical Radius (F != 1) */
double Eval = - 1.0/(x*x) + Mq * ( 1.0/SQR(x-1) - 1.0 ) + F*F*x*( Mq + 1);
return(Eval);
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Zero at LagrangeTwo
@param (double) x Position on x axix
@return (double)
@heading Light Curve
****************************************************************************/
double LagrangeTwo(double x)
{
/* Uses Global Variable Mq (Mass Ratio) */
/* Eval is Zero at LagrangeTwo (F=1) */
/* or at critical Radius (F != 1) */
double Eval = - 1.0/(x*x)
+ Mq * ( (1.0-x)/((x-1.0)*SQR(x-1.0)) - 1.0 )
+ F*F*x*( Mq + 1);
return(Eval);
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Zero at LagrangeThree
@param (double) x Position on x axix
@return (double)
@heading Light Curve
****************************************************************************/
double LagrangeThree(double x)
{
/* Uses Global Variable Mq (Mass Ratio) */
/* Eval is Zero at LagrangeThree (F=1) */
/* or at critical Radius (F != 1) */
double Eval = - x/((-x)*(-x)*(-x))
+ Mq * ( (1.0)/SQR(1.0-x) - 1.0 )
+ F*F*x*( Mq + 1);
return(Eval);
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Find root of a function in [Low,Up] with Tolerance
@long Uses Brents algorithm. f(Low) and f(Up) must have opposite signs,
thus indicating that they bracket a root of f().
If this is not the case, the program exits with an error message.
Adopted from procedure 'zero' in: Richard Brent, Algorithms for
minimization without derivatives, Prentice-Hall, Inc. (1973)
@param (double)(*Function)(double) Pointer to the function
@param (const double) Low Lower limit
@param (const double) Up Upper limit
@param (const double) Tolerance The convergence criterion
@param (const char*) caller The calling function
@param (int*) Err The exit status
@return (double)
@heading Light Curve
****************************************************************************/
double RootFind(double (*Function)(double),
const double Low, const double Up,
const double Tolerance, const char *caller, int *Err)
{
register unsigned long i = 0; /* loop variable */
double a = Low, b = Up, c = Up; /* storage */
double Delta = 0.0, e = 0.0; /* convergence testing */
double test1, test2; /* convergence testing */
double fa = (*Function)(a); /* function value */
double fb = (*Function)(b); /* function value */
double fc; /* function value */
double p, q, r, s; /* for quadratic interpol. */
double Tolerance1 = 0.0; /* convergence testing */
double Bisect, Eval; /* bisection */
static double laststore = 0.5; /* last root found */
register unsigned long j = 1; /* loop variable */
int is_ok = 0; /* brackets found */
*Err = 0;
/* -------- whether interval (a,b) brackets a root ------------------- */
if ( (fa >= FLT_EPSILON && fb >= FLT_EPSILON)
|| (fa <= (-FLT_EPSILON) && fb <= (-FLT_EPSILON) )) {
/* root not bracketed, try last value */
WARNING(_("Root not bracketed, try to recover"));
while (j < 500 && is_ok == 0) {
fa = (*Function)(laststore - j * 1.0e-4);
fb = (*Function)(laststore + j * 1.0e-4);
if ( (fa <= (-FLT_EPSILON) && fb >= FLT_EPSILON)
|| (fb <= (-FLT_EPSILON) && fa >= FLT_EPSILON)) {
is_ok = 1;
a = laststore - j * 1e-4;
b = laststore + j * 1e-4;
} else {
++j;
}
}
/* no success */
if (is_ok == 0) {
if (Flags.interactive == OFF) {
fprintf(stderr,
_("** ERROR **: %s: RootFind: Root not bracketed \n"), caller);
fprintf(stderr,
_("This is likely due to an extreme combination of\n"));
fprintf(stderr, _("mass ratio and fill factor\n"));
abort();
} else {
fprintf(stderr,
_("** ERROR **: %s: RootFind: Root not bracketed \n"), caller);
fprintf(stderr,
_("This is likely due to an extreme combination of\n"));
fprintf(stderr, _("mass ratio and fill factor\n"));
*Err = 1;
return(0.0);
}
}
}
/* -------------- initialize c -------------------------------------- */
fc = fb;
/* -------------- loop until convergence ------------------------------ */
for (i = 0; i < MAXITER; ++i) {
/* If (Root between a and b) then put it between b and c */
if ((fc >= FLT_EPSILON && fb >= FLT_EPSILON)
|| (fc <= (-FLT_EPSILON) && fb <= (-FLT_EPSILON) )) {
c = a;
fc = fa;
e = Delta = b-a;
}
/* make fb nearest to x-axis */
if (fabs(fc) <= fabs(fb)) {
a = b; fa = fb;
b = c; fb = fc;
c = a; fc = fa;
}
/* this is the bisection step */
Bisect = (c-b) / 2.0;
/* set tolerance for convergence testing */
/* in the following sum, the first term should ensure proper */
/* convergence testing in the case of a very steep function */
Tolerance1 = (2.0 * FLT_EPSILON * fabs(b)) + (0.5 * Tolerance);
/* Convergence Test */
if ((fabs(Bisect) <= Tolerance1) || (fabs(fb) <= (2.0*FLT_EPSILON))) {
laststore = b;
return(b);
}
/* Interpolation Step if Fast Enough */
if ((fabs(e) >= Tolerance1) && ((fabs(a) - fabs(b)) >= FLT_EPSILON )) {
s = fb/fa;
if (fabs(a - c) <= FLT_EPSILON) {
/* - only two function values, do linear interpolation ------ */
p = 2.0 * Bisect * s;
q = 1.0 - s;
} else {
/* - inverse quadratic interpolation ------- */
q = fa/fc;
r = fb/fc;
p = s * (2.0*Bisect*q*(q-r)-(b-a)*(r-1.0));
q = (q - 1.0) * ( r - 1.0) * (s - 1.0);
}
if (p <= 0.0) p = -p;
else q = -q;
s = e;
e = Delta;
/* Check Bounds */
test1 = (3.0 * Bisect * q) - fabs(Tolerance1 * q);
test2 = fabs( s * q );
if ( ((2.0*p) >= test1) || (p >= test2) ) {
/* we resort to bisection */
Delta = Bisect;
e = Delta;
} else {
/* we accept interpolation */
Delta = p/q;
}
/* Interpolation is creepingly slow, do Bisection */
} else {
Delta = Bisect;
e = Delta;
}
/* we reset a to the former value of b */
a = b; fa = fb;
/* and set b to its updated value */
/* we step at least by Tolerance1 */
if (fabs(Delta) >= Tolerance1) {
b = b + Delta;
} else {
if (Bisect >= DBL_EPSILON) {
Eval = fabs(Tolerance1);
} else {
Eval = -fabs(Tolerance1);
}
b = b + Eval;
}
fb = ((*Function)(b));
}
return(0.0);
}
syntax highlighted by Code2HTML, v. 0.9.1