/* NIGHTFALL Light Curve Synthesis Program */ /* Copyright (C) 1998 Rainer Wichmann */ /* */ /* This program is free software; you can redistribute it */ /* and/or modify */ /* it under the terms of the GNU General Public License as */ /* published by */ /* the Free Software Foundation; either version 2 of the License, or */ /* (at your option) any later version. */ /* */ /* This program is distributed in the hope that it will be useful, */ /* but WITHOUT ANY WARRANTY; without even the implied warranty of */ /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */ /* GNU General Public License for more details. */ /* */ /* You should have received a copy of the GNU General Public License */ /* along with this program; if not, write to the Free Software */ /* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */ #include #include #include #include #include #include #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); }