/* 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 <stdlib.h>
#include <ctype.h>
#include <float.h>
#include "Light.h"
/******************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Parse the command line
@tip This is the most awkward function in the whole program ...
@param (int) argc The number of arguments
@param (char *[]) argv The arguments
@return (void)
@heading Main
*******************************************************************/
/* Get The Command Line Arguments */
void GetArguments(int *argc, char *argv[])
{
int i; /* loop variable */
int Ichar;
int CheckIn, mandatory = 0;
char infile[256];
char ErrMsg[256];
char *debug;
double G_SI = 6.6726e-11; /* gravitational constant */
int FitDim = 0;
char Fitopt[] = "0123456789ABCDEFGHIJKLMNOPabcdefghijkl";
/* these are globals, we set the default values */
if (PHASESTEPS > 80) PhaseSteps = 80;
else PhaseSteps = PHASESTEPS;
StepsPerPi = STEPS_PER_PI;
MaximumElements = MAXELEMENTS;
/* optional arguments */
while ((*argc) > 1) {
if ((argv[1][0] == '-') && (argv[1][1] != '-')) {
switch (argv[1][1]) {
case 'h':
UsageInfo();
MPI_Abort(MPI_COMM_WORLD, EXIT_FAILURE);
exit (0);
break;
case 'D':
debug = &argv[1][2];
while (debug != NULL && isalpha (*debug)) {
if ( ! isalpha (*debug) ) {
WARNERR(_("Invalid character in debug string encountered, skipped"),
argv[1]);
} else {
if ( (*debug) == 'w') Flags.debug[WARN] = 1;
else if ( (*debug) == 'v') Flags.debug[VERBOSE] = 1;
else if ( (*debug) == 'b') Flags.debug[BUSY] = 1;
else if ( (*debug) == 'g') Flags.debug[GEO] = 1;
else if ( (*debug) == 'i') Flags.debug[INFILE] = 1;
else if ( (*debug) == 'm') Flags.debug[MINFIND] = 1;
else if ( (*debug) == 'q') Flags.debug[CHISQR] = 1;
else if ( (*debug) == 'o') Flags.debug[OPTIMIZE] = 1;
else if ( (*debug) == 's') Flags.debug[SURFACE] = 1;
}
debug++;
}
break;
case 'U':
Flags.interactive = ON;
break;
case 'A':
Flags.animate = ON;
break;
#ifdef _WITH_OPENGL
case 'g':
Flags.use_opengl = OFF;
break;
#endif
case 'V':
if (argv[1][2] == 'i') Flags.visualize = 2;
else if (argv[1][2] == 'c') Flags.visualize = 3;
else if (argv[1][2] == 'a') Flags.visualize = 4;
else Flags.visualize = ON;
break;
case 'G':
if (argv[1][2] == 'P') Flags.plot = 3;
else if (argv[1][2] == '2') Flags.plot = 2;
else if (argv[1][2] == 'S') Flags.plot = 4;
else Flags.plot = ON;
break;
case 'H':
Flags.eps = ON;
break;
case 'N':
CheckIn = sscanf(&argv[1][2], "%d", &PhaseSteps);
if (CheckIn != 1) {
WARNERR(_("Steps for Lightcurve (maybe a 'blank' after 'N'), Default used"),
argv[1] );
PhaseSteps = IMIN(PHASESTEPS, 80);
}
else if (PhaseSteps > PHASESTEPS) {
sprintf(ErrMsg, _("Lightcurve Steps > Maximum (%d), Default used"),
PHASESTEPS);
WARNERR(ErrMsg, argv[1] );
PhaseSteps = IMIN(PHASESTEPS, 80);
}
else if (PhaseSteps < 0) {
WARNERR(_("Lightcurve Steps < Zero, Default used"), argv[1] );
PhaseSteps = IMIN(PHASESTEPS, 80);
}
break;
case 'P':
CheckIn = sscanf(&argv[1][2], "%lf", &Orbit.LambdaZero);
if (CheckIn != 1) {
WARNERR(_("Line Wavelength missing (maybe a 'blank' after 'P'), default used"),
argv[1] );
Orbit.LambdaZero = LAMDAZERO_DEFAULT;
}
else if ( (Orbit.LambdaZero <= (LIM_PRF_L-FLT_EPSILON))
|| (Orbit.LambdaZero >= (LIM_PRF_H+FLT_EPSILON*LIM_PRF_H))){
sprintf(ErrMsg,
_("Line Profile Rest Wavelength must be in [%f, %f]"),
LIM_PRF_L, LIM_PRF_H );
WARNERR(ErrMsg, argv[1]);
Orbit.LambdaZero = LAMDAZERO_DEFAULT;
}
Flags.lineprofile = ON;
break;
case 'X':
debug = &argv[1][2];
if ( ! isalnum (*debug) ) {
sprintf(ErrMsg, _("Bad Option for Simplex Fit %c, skipped"),
*debug);
WARNERR(ErrMsg, argv[1]);
}
while ( *debug != '\0' ) {
for (Ichar = 0; Ichar < (int) strlen(Fitopt); ++Ichar) {
if (*debug == Fitopt[Ichar]) Flags.simplex[Ichar] = ON;
}
debug++;
}
if (argv[2] != NULL) CheckIn = sscanf(argv[2],"%f", &Flags.SimplexTol);
else CheckIn = 0;
if (CheckIn != 1) {
WARNERR(_("Simplex Fit: No Tolerance specified, default used"),
argv[1]);
Flags.SimplexTol = SIMPLEX_FIT_TOLERANCE;
} else {
++argv;
--(*argc);
}
Flags.WantFit = ON;
if (fabs (Flags.SimplexTol - 0.001) <= FLT_EPSILON) Flags.anneal = ON;
break;
case 'Y':
debug = &argv[1][2];
if ( ! isalnum (*debug) ) {
sprintf(ErrMsg,_("Bad Option for ChiSquare Map %c, skipped"),
*debug);
WARNERR(ErrMsg, argv[1]);
}
while ( *debug != '\0' ) {
for (Ichar = 0; Ichar < (int) strlen(Fitopt); ++Ichar) {
if (*debug == Fitopt[Ichar]) Flags.simplex[Ichar] = ON;
}
debug++;
}
++argv;
--(*argc);
if (argv[1] != NULL) CheckIn = sscanf(argv[1],"%f", &Flags.Step[0]);
else CheckIn = 0;
if (CheckIn != 1) INERR(_("Chi Square Map: No Step1 specified"),
argv[1])
++argv;
--(*argc);
if (argv[1] != NULL) CheckIn = sscanf(argv[1],"%f", &Flags.Step[1]);
else CheckIn = 0;
if (CheckIn != 1) INERR(_("Chi Square Map: No Step2 specified"),
argv[1])
Flags.WantMap = ON;
break;
case 'M':
Flags.blackbody = OFF;
break;
case 'F':
Flags.fractional = ON;
break;
case 'L':
(void) sscanf(&argv[1][2], "%d", &Flags.limb);
if (Flags.limb > LD_LAW_MAX || Flags.limb < 0) {
WARNERR(_("Invalid limb darkening option (maybe 'blank' after 'L'), linear used"),
argv[1]);
Flags.limb = 0;
}
break;
case 'R':
(void) sscanf(&argv[1][2], "%d", &Flags.reflect);
if (Flags.reflect > 9 || Flags.reflect < 1) {
WARNERR(_("Invalid Reflection Iteration Count (maybe 'blank' after 'R'), 3 used"),
argv[1]);
Flags.reflect = 3;
}
break;
case 'B':
switch (argv[1][2]) {
case 'U':
Flags.PlotBand = Umag;
break;
case 'B':
Flags.PlotBand = Bmag;
break;
case 'V':
Flags.PlotBand = Vmag;
break;
case 'R':
Flags.PlotBand = Rmag;
break;
case 'I':
Flags.PlotBand = Imag;
break;
case 'u':
Flags.PlotBand = umag;
break;
case 'v':
Flags.PlotBand = vmag;
break;
case 'b':
Flags.PlotBand = bmag;
break;
case 'y':
Flags.PlotBand = ymag;
break;
case 'J':
Flags.PlotBand = Jmag;
break;
case 'H':
Flags.PlotBand = Hmag;
break;
case 'K':
Flags.PlotBand = Kmag;
break;
case '1':
Flags.PlotBand = mag1;
break;
case '2':
Flags.PlotBand = mag2;
break;
default:
WARNERR(_("Bad Option For Bandpass to Graph, V Band used"), argv[1])
break;
}
break;
case '3':
switch (argv[1][2]) {
case 'U':
CheckIn = sscanf(&argv[1][3], "%lf", &Orbit.Third[Umag]);
if (CheckIn != 1)
WARNERR(_("No Third Light given (misplaced 'blank' ?)"), argv[1])
break;
case 'B':
CheckIn = sscanf(&argv[1][3], "%lf", &Orbit.Third[Bmag]);
if (CheckIn != 1)
WARNERR(_("No Third Light given (misplaced 'blank' ?)"), argv[1])
break;
case 'V':
CheckIn = sscanf(&argv[1][3], "%lf", &Orbit.Third[Vmag]);
if (CheckIn != 1)
WARNERR(_("No Third Light given (misplaced 'blank' ?)"), argv[1])
break;
case 'R':
CheckIn = sscanf(&argv[1][3], "%lf", &Orbit.Third[Rmag]);
if (CheckIn != 1)
WARNERR(_("No Third Light given (misplaced 'blank' ?)"), argv[1])
break;
case 'I':
CheckIn = sscanf(&argv[1][3], "%lf", &Orbit.Third[Imag]);
if (CheckIn != 1)
WARNERR(_("No Third Light given (misplaced 'blank' ?)"), argv[1])
break;
case 'u':
CheckIn = sscanf(&argv[1][3], "%lf", &Orbit.Third[umag]);
if (CheckIn != 1)
WARNERR(_("No Third Light given (misplaced 'blank' ?)"), argv[1])
break;
case 'v':
CheckIn = sscanf(&argv[1][3], "%lf", &Orbit.Third[vmag]);
if (CheckIn != 1)
WARNERR(_("No Third Light given (misplaced 'blank' ?)"), argv[1])
break;
case 'b':
CheckIn = sscanf(&argv[1][3], "%lf", &Orbit.Third[bmag]);
if (CheckIn != 1)
WARNERR(_("No Third Light given (misplaced 'blank' ?)"), argv[1])
break;
case 'y':
CheckIn = sscanf(&argv[1][3], "%lf", &Orbit.Third[ymag]);
if (CheckIn != 1)
WARNERR(_("No Third Light given (misplaced 'blank' ?)"), argv[1])
break;
case 'J':
CheckIn = sscanf(&argv[1][3], "%lf", &Orbit.Third[Jmag]);
if (CheckIn != 1)
WARNERR(_("No Third Light given (misplaced 'blank' ?)"), argv[1])
break;
case 'H':
CheckIn = sscanf(&argv[1][3], "%lf", &Orbit.Third[Hmag]);
if (CheckIn != 1)
WARNERR(_("No Third Light given (misplaced 'blank' ?)"), argv[1])
break;
case 'K':
CheckIn = sscanf(&argv[1][3], "%lf", &Orbit.Third[Kmag]);
if (CheckIn != 1)
WARNERR(_("No Third Light given (misplaced 'blank' ?)"), argv[1])
break;
default:
WARNERR(_("Bad Option For Third Light"), argv[1])
break;
}
break;
case 'C':
if ((*argc) < 3) {
UsageInfo();
MPI_Abort(MPI_COMM_WORLD, EXIT_FAILURE);
exit(EXIT_FAILURE);
}
else if (strlen(argv[2]) < (sizeof(infile)-1)) {
(void) sscanf(argv[2], "%s", infile);
ParseConfig(infile, &mandatory);
++argv; --(*argc);
}
else {
WARNERR(_("Path to Config File exceeds Maximum Length (255), skipped"),
argv[2]);
++argv; --(*argc);
}
break;
case 'I':
if ((*argc) < 3) {
UsageInfo();
MPI_Abort(MPI_COMM_WORLD, EXIT_FAILURE);
exit(EXIT_FAILURE);
}
else if (strlen(argv[2]) < (sizeof(infile)-1)) {
(void) sscanf(argv[2], "%s", infile);
Read(infile);
++argv; --(*argc);
}
else {
WARNING(_("Path to Input File exceeds Maximum Length (255), skipped"));
++argv; --(*argc);
}
break;
case 's':
if ((*argc) < 5) {
UsageInfo();
MPI_Abort(MPI_COMM_WORLD, EXIT_FAILURE);
exit(EXIT_FAILURE);
}
switch (argv[1][2]) {
case 'P':
++argv; --(*argc);
CheckIn = sscanf(argv[1],
"%lf", &Spot[Primary][Flags.Spots1].longitude);
if (CheckIn != 1) INERR(_("No Longitude for Spot"), argv[1] );
if( (Spot[Primary][Flags.Spots1].longitude-LIM_SLO_L)
<=(-FLT_EPSILON) ||
(Spot[Primary][Flags.Spots1].longitude-LIM_SLO_H)>=FLT_EPSILON)
INERR(_("Spot longitude must be in [-360,360]"), argv[1]);
++argv; --(*argc);
CheckIn = sscanf(argv[1],
"%lf", &Spot[Primary][Flags.Spots1].latitude);
if (CheckIn != 1) INERR(_("No Latitude for Spot"), argv[1]);
if((Spot[Primary][Flags.Spots1].latitude-LIM_SLA_L)
<=(-FLT_EPSILON) ||
(Spot[Primary][Flags.Spots1].latitude-LIM_SLA_H)>=FLT_EPSILON)
INERR(_("Spot latitude must be in [-90,90]"), argv[1]);
++argv; --(*argc);
CheckIn = sscanf(argv[1],"%lf",&Spot[Primary][Flags.Spots1].radius);
if (CheckIn != 1) INERR(_("No Radius for Spot"), argv[1]);
if((Spot[Primary][Flags.Spots1].radius-LIM_SRA_L)<=(-FLT_EPSILON)||
(Spot[Primary][Flags.Spots1].radius-LIM_SRA_H)>=FLT_EPSILON)
INERR(_("Spot radius must be in [1,180]"), argv[1] );
++argv; --(*argc);
CheckIn = sscanf(argv[1],
"%lf", &Spot[Primary][Flags.Spots1].dimfactor);
if (CheckIn != 1) INERR(_("No Dim Factor for Spot"), argv[1]);
if((Spot[Primary][Flags.Spots1].dimfactor-LIM_SDF_L)
<=(-FLT_EPSILON))
INERR(_("Spot Dim Factor must be > 0.5"), argv[1]);
if((Spot[Primary][Flags.Spots1].dimfactor-LIM_SDF_H)>=FLT_EPSILON)
INERR(_("Spot Dim Factor must be < 2. "), argv[1]);
++Flags.Spots1;
break;
case 'S':
++argv; --(*argc);
CheckIn = sscanf(argv[1],
"%lf", &Spot[Secondary][Flags.Spots2].longitude);
if (CheckIn != 1) INERR(_("No Longitude for Spot"), argv[1] );
if((Spot[Secondary][Flags.Spots2].longitude-LIM_SLO_L)
<=(-FLT_EPSILON) ||
(Spot[Secondary][Flags.Spots2].longitude-LIM_SLO_H)>=FLT_EPSILON)
INERR(_("Spot longitude must be in [-360,360]"), argv[1]);
++argv; --(*argc);
CheckIn = sscanf(argv[1],
"%lf", &Spot[Secondary][Flags.Spots2].latitude);
if (CheckIn != 1) INERR(_("No Latitude for Spot"), argv[1]);
if((Spot[Secondary][Flags.Spots2].latitude-LIM_SLA_L)
<=(-FLT_EPSILON) ||
(Spot[Secondary][Flags.Spots2].latitude-LIM_SLA_H)>=FLT_EPSILON)
INERR(_("Spot latitude must be in [-90,90]"), argv[1]);
++argv; --(*argc);
CheckIn = sscanf(argv[1],
"%lf", &Spot[Secondary][Flags.Spots2].radius);
if (CheckIn != 1) INERR(_("No Radius for Spot"), argv[1]);
if((Spot[Secondary][Flags.Spots2].radius-LIM_SRA_L)
<=(-FLT_EPSILON) ||
(Spot[Secondary][Flags.Spots2].radius-LIM_SRA_H)>=FLT_EPSILON)
INERR(_("Spot radius must be in [1,180]"), argv[1] );
++argv; --(*argc);
CheckIn = sscanf(argv[1],
"%lf", &Spot[Secondary][Flags.Spots2].dimfactor);
if (CheckIn != 1) INERR(_("No Dim Factor for Spot"), argv[1]);
if((Spot[Secondary][Flags.Spots2].dimfactor-LIM_SDF_L)
<=(-FLT_EPSILON))
INERR(_("Spot Dim Factor must be > 0.5"), argv[1]);
if((Spot[Secondary][Flags.Spots2].dimfactor-LIM_SDF_H)
>=FLT_EPSILON)
INERR(_("Spot Dim Factor must be < 2."), argv[1]);
++Flags.Spots2;
break;
default:
WARNERR(_("Bad option for spot, skipped"), argv[1]);
break;
}
break;
case 't':
if ((*argc) < 3) {
UsageInfo();
MPI_Abort(MPI_COMM_WORLD, EXIT_FAILURE);
exit(EXIT_FAILURE);
}
switch (argv[1][2]) {
case 'P':
++argv;
--(*argc);
CheckIn = sscanf(argv[1], "%lf", &Orbit.TruePeriod);
if (CheckIn != 1)
WARNERR(_("No Period given (maybe a 'blank' after 'P')"), argv[1]);
if (Orbit.TruePeriod <= LIM_PERI_L)
INERR(_("Absolute Period"), argv[1]);
/* convert days to seconds */
Orbit.TruePeriod = Orbit.TruePeriod * 86400.;
break;
case 'M':
++argv;
--(*argc);
CheckIn = sscanf(argv[1], "%lf", &Orbit.TrueMass);
if (CheckIn != 1)
WARNERR(_("No Mass given (maybe a 'blank' after 'M')"), argv[1]);
if (Orbit.TrueMass <= LIM_MASS_L) INERR(_("Absolute Mass"), argv[1]);
/* convert solar masses to KG */
Orbit.TrueMass = Orbit.TrueMass * 1.989E30;
break;
case 'D':
++argv;
--(*argc);
CheckIn = sscanf(argv[1], "%lf", &Orbit.TrueDistance);
if (CheckIn != 1)
WARNERR(_("No Separation given (maybe a 'blank' after 'D')"),
argv[1]);
if (Orbit.TrueDistance <= LIM_DIST_L)
INERR(_("Absolute Separation"), argv[1]);
/* convert solar radius to meter */
Orbit.TrueDistance = Orbit.TrueDistance * 6.960E8;
break;
default:
WARNERR(_("Bad option for absolute values, default used"), argv[1]);
break;
}
break;
case 'e':
if ((*argc) < 4) {
UsageInfo();
MPI_Abort(MPI_COMM_WORLD, EXIT_FAILURE);
exit(EXIT_FAILURE);
}
++argv;
--(*argc);
CheckIn = sscanf(argv[1],"%lf", &Orbit.Excentricity);
if (CheckIn != 1) INERR(_("No Eccentricity found"), argv[1]);
if ( (Orbit.Excentricity - LIM_EX_L) <= (-FLT_EPSILON)
|| (Orbit.Excentricity - LIM_EX_H) >= FLT_EPSILON )
INERR(_("Eccentricity must be in (0., 0.95)"), argv[1]);
++argv;
--(*argc);
CheckIn = sscanf(argv[1],"%lf", &Orbit.Omega);
if (CheckIn != 1) INERR(_("No Periastron Length found"), argv[1]);
if ((Orbit.Omega - LIM_PA_L) <= (-FLT_EPSILON)
|| (Orbit.Omega - LIM_PA_H) >= FLT_EPSILON )
INERR(_("Periastron Length must be in [-360, 360]"), argv[1]);
Flags.elliptic = ON;
break;
case 'f':
if ((*argc) < 3) {
UsageInfo();
MPI_Abort(MPI_COMM_WORLD, EXIT_FAILURE);
exit(EXIT_FAILURE);
}
switch (argv[1][2]) {
case 'P':
++argv;
--(*argc);
CheckIn = sscanf(argv[1],"%lf", &Binary[Primary].Fratio);
if (CheckIn != 1) INERR(_("No Rotation Rate found"), argv[1]);
if( (Binary[Primary].Fratio - LIM_FR_L) <= (-FLT_EPSILON) ||
(Binary[Primary].Fratio - LIM_FR_H) >= FLT_EPSILON )
INERR(_("Rotation Rate must be in (0., 10.]"), argv[1]);
Flags.asynchron1 = ON;
break;
case 'S':
++argv;
--(*argc);
CheckIn = sscanf(argv[1],"%lf", &Binary[Secondary].Fratio);
if (CheckIn != 1) INERR(_("No Rotation Rate found"), argv[1]);
if( (Binary[Secondary].Fratio - LIM_FR_L) <= (-FLT_EPSILON) ||
(Binary[Primary].Fratio - LIM_FR_H) >= FLT_EPSILON )
INERR(_("Rotation Rate must be in (0., 10.]"), argv[1]);
Flags.asynchron2 = ON;
break;
default:
WARNERR(_("Bad Option (use 'nightfall <RETURN> to get a list of options): "),
argv[1]);
break;
}
break;
case 'O':
if ((*argc) < 3) {
UsageInfo();
MPI_Abort(MPI_COMM_WORLD, EXIT_FAILURE);
exit(EXIT_FAILURE);
}
switch (argv[1][2]) {
case 'P':
++argv;
--(*argc);
CheckIn = sscanf(argv[1],"%lf", &Binary[Primary].log_g);
if (CheckIn != 1) INERR(_("No log g found"), argv[1]);
if( (Binary[Primary].log_g - LIM_LOGG_L) <= (-FLT_EPSILON) ||
(Binary[Primary].log_g - LIM_LOGG_H) >= FLT_EPSILON )
INERR(_("Log g must be in [3.5, 5.0]"), argv[1]);
Binary[Primary].log_g = 0.5 * ROUND(2.0 * Binary[Primary].log_g);
break;
case 'S':
++argv;
--(*argc);
CheckIn = sscanf(argv[1],"%lf", &Binary[Secondary].log_g);
if (CheckIn != 1) INERR(_("No log_g found"), argv[1]);
if( (Binary[Secondary].log_g - LIM_LOGG_L) <= (-FLT_EPSILON) ||
(Binary[Primary].log_g - LIM_LOGG_H) >= FLT_EPSILON )
INERR(_("Log g must be in [3.5, 5.0]"), argv[1]);
Binary[Secondary].log_g = 0.5 * ROUND(2.0 * Binary[Secondary].log_g);
break;
default:
WARNERR(_("Bad Option (use 'nightfall <RETURN> to get a list of options): "),
argv[1]);
break;
}
break;
case 'a':
if ((*argc) < 3) {
UsageInfo();
MPI_Abort(MPI_COMM_WORLD, EXIT_FAILURE);
exit(EXIT_FAILURE);
}
switch (argv[1][2]) {
case 'P':
++argv;
--(*argc);
SetAlbedo(Primary, argv[1]);
break;
case 'S':
++argv;
--(*argc);
SetAlbedo(Secondary, argv[1]);
break;
default:
WARNERR(_("Bad Option (use 'nightfall <RETURN> to get a list of options): "),
argv[1]);
break;
}
break;
default:
WARNERR(_("Bad Option (use 'nightfall <RETURN> to get a list of options): "),
argv[1]);
}
++argv;
--(*argc);
} else if ((argv[1][0] == '-') && (argv[1][1] == '-')) {
/* --------------- ignore ----------------------------------------- */
++argv;
--(*argc);
if ( (*argc) > 1 ) {
if (argv[1][0] != '-') {
++argv;
--(*argc);
}
}
} else {
/* -------------- Argument did not start with a '-' or '--' ----- */
break;
}
}
/* >>>>>>>>>>>>>>>>>>>>>>>> mandatory arguments <<<<<<<<<<<<<<<<<<< */
if ( ((*argc) < 7) && (mandatory < 6) && (Flags.interactive == OFF))
nf_error(_("One of the mandatory Arguments is missing. Type 'nightfall' to get more info. (One possible reason might be a missing argument for a preceding option)"));
/*
+ now we scan the rest of the CL. if there are gtk arguments at end,
+ they will be passed. if they are in front of the mandatory arguments,
+ a read error will occur (hopefully)
*/
if ((*argc) >= 7 ) {
(void) sscanf(argv[1], "%lf", &Binary[Primary].Mq);
if((Binary[Primary].Mq - LIM_MQ_L) <= (-FLT_EPSILON) ||
(Binary[Primary].Mq - LIM_MQ_H) >= FLT_EPSILON )
INERR(_("mass ratio must be in [0.001,50.]"), argv[1] );
Binary[Secondary].Mq = 1.0/Binary[Primary].Mq;
++argv;
--(*argc);
(void) sscanf(argv[1], "%lf", &Orbit.Inclination);
if((Orbit.Inclination-LIM_IN_L) <= (-FLT_EPSILON)
|| (Orbit.Inclination - RTOD*LIM_IN_H) >= FLT_EPSILON)
INERR(_("Inclination must be in [0,90]"), argv[1]);
Orbit.Inclination = DTOR * Orbit.Inclination;
++argv;
--(*argc);
(void) sscanf(argv[1], "%lf", &Binary[Primary].RocheFill);
if((Binary[Primary].RocheFill - LIM_RF_L) <= (-FLT_EPSILON) ||
(Binary[Primary].RocheFill - LIM_RO_H) >= FLT_EPSILON)
INERR(_("Roche lobe filling must be in [0.001,1.3]"),argv[1] );
++argv;
--(*argc);
(void) sscanf(argv[1], "%lf", &Binary[Secondary].RocheFill);
if((Binary[Secondary].RocheFill - LIM_RF_L) <= (-FLT_EPSILON) ||
(Binary[Secondary].RocheFill - LIM_RO_H) >= FLT_EPSILON)
INERR(_("Roche lobe filling must be in [0.001,1.3]"),argv[1] );
++argv;
--(*argc);
(void) sscanf(argv[1], "%lf", &Binary[Primary].Temperature);
if(Flags.blackbody == OFF) {
if((Binary[Primary].Temperature - LIM_TM_L ) <= (-FLT_EPSILON) ||
(Binary[Primary].Temperature - LIM_TM_H(Binary[Primary].log_g)) >= FLT_EPSILON)
INERR(_("Temperature must be in [3000,35000]"), argv[1]);
} else {
if((Binary[Primary].Temperature - LIM_TB_L) <= (-FLT_EPSILON) ||
(Binary[Primary].Temperature - LIM_TB_H )>= FLT_EPSILON)
INERR(_("Temperature must be in [350,350000]"), argv[1]);
}
++argv;
--(*argc);
(void) sscanf(argv[1], "%lf", &Binary[Secondary].Temperature);
if(Flags.blackbody == OFF) {
if((Binary[Secondary].Temperature - LIM_TM_L) <= (-FLT_EPSILON) ||
(Binary[Secondary].Temperature - LIM_TM_H(Binary[Secondary].log_g)) >= FLT_EPSILON)
INERR(_("Temperature must be in [3000,35000]"), argv[4]);
} else {
if((Binary[Secondary].Temperature - LIM_TB_L) <= (-FLT_EPSILON) ||
(Binary[Secondary].Temperature - LIM_TB_H) >= FLT_EPSILON)
INERR(_("Temperature must be in [350,350000]"), argv[1]);
}
++argv;
--(*argc);
}
/* >>>>>>>>>>>>>>>>>>> ADJUST INPUT <<<<<<<<<<<<<<<<<<<<<<<< */
for (i = 0; i < 32; ++i) {
if (Flags.simplex[i] == ON) ++FitDim;
}
if (FitDim < 2 && Flags.WantFit == ON) {
Flags.WantFit = OFF;
WARNING(_("Less than two Parameters to fit, discard Fit Option"));
}
/* no overcontact if elliptic */
if ( (Flags.elliptic == ON)
&& (Binary[Primary].RocheFill >= (1+FLT_EPSILON) ) ) {
Binary[Primary].RocheFill = 1.0;
WARNING(_("Eccentric Orbit --> Decrease Roche Fill Factor to One"));
}
if ( (Flags.elliptic == ON)
&& (Binary[Secondary].RocheFill >= (1+FLT_EPSILON) ) ) {
Binary[Secondary].RocheFill = 1.0;
WARNING(_("Eccentric Orbit --> Decrease Roche Fill Factor to One"));
}
/* no overcontact if disk */
/* no disk if elliptic */
#ifdef HAVE_DISK
if ( (Flags.elliptic == ON) && (Flags.disk == ON) )
{
Flags.disk = OFF;
WARNING(_("Eccentric Orbit --> No accretion disk possible"));
}
if ( (Flags.disk == ON)
&& (Binary[Primary].RocheFill >= (1+FLT_EPSILON) ) ) {
Binary[Primary].RocheFill = 1.0;
WARNING(_("Eccentric Orbit --> Decrease Roche Fill Factor to One"));
}
if ( (Flags.disk == ON)
&& (Binary[Secondary].RocheFill >= (1+FLT_EPSILON) ) ) {
Binary[Secondary].RocheFill = 1.0;
WARNING(_("Eccentric Orbit --> Decrease Roche Fill Factor to One"));
}
#endif
/* common fillfactor if overcontact */
if ( Binary[Primary].RocheFill >= (1+FLT_EPSILON)
|| Binary[Secondary].RocheFill >= (1+FLT_EPSILON) ){
if ( fabs(Binary[Primary].RocheFill-Binary[Secondary].RocheFill)
>= FLT_EPSILON)
WARNING(_("Overcontact --> Set Common Roche Fill Factor to Maximum of Both"));
Binary[Primary].RocheFill =
MAX(Binary[Primary].RocheFill,Binary[Secondary].RocheFill);
Binary[Secondary].RocheFill = Binary[Primary].RocheFill;
Flags.fill = ON;
} else {
Flags.fill = OFF;
}
/* Limits for Mass Ratio */
if (Flags.fill == ON) {
if (Binary[Primary].Mq >= LIM_MQO_H) {
Binary[Primary].Mq = LIM_MQO_H;
Binary[Secondary].Mq = 1.0/Binary[Primary].Mq;
WARNING(_("Overcontact --> Reduce Mass Ratio to limit"));
}
} else {
if (Binary[Primary].Mq <= LIM_MQ_L ) {
Binary[Primary].Mq = LIM_MQ_L;
Binary[Secondary].Mq = 1.0/Binary[Primary].Mq;
WARNING(_("Avoid Numerical Problems --> Increase Mass Ratio to limit"));
}
}
/* no nonsynchroneous rotation if overcontact */
if (Flags.fill == ON) {
if ( (fabs(Binary[Primary].Fratio - 1.0) >= FLT_EPSILON) ||
(fabs(Binary[Secondary].Fratio - 1.0) >= FLT_EPSILON) ) {
WARNING(_("Overcontact --> Reset NonSync Rotation to Unity"));
Binary[Primary].Fratio = 1.0;
Binary[Secondary].Fratio = 1.0;
Flags.asynchron1 = OFF;
Flags.asynchron2 = OFF;
}
}
/* ----------- set default band to plot -------------------- */
if (Flags.PlotBand == -1) Flags.PlotBand = Vmag;
ComputeGravDark();
/* >>>>>>>>>>> consistency of mass,period, and distance <<<<<<<<<<<<<< */
if (fabs(Orbit.TruePeriod) >= FLT_EPSILON
&& fabs(Orbit.TrueDistance ) >= FLT_EPSILON ) {
Orbit.TrueMass =
(Orbit.TrueDistance * Orbit.TrueDistance) * Orbit.TrueDistance
* 4.0 * PI * PI
/ (G_SI * (Orbit.TruePeriod * Orbit.TruePeriod));
}
else if (fabs(Orbit.TruePeriod ) >= FLT_EPSILON
&& fabs(Orbit.TrueMass ) >= FLT_EPSILON ) {
Orbit.TrueDistance =
pow((Orbit.TruePeriod*Orbit.TruePeriod)
* G_SI * Orbit.TrueMass / (4.0 * (PI*PI)),
(1.0/3.0));
}
else if (fabs(Orbit.TrueDistance) >= FLT_EPSILON
&& fabs(Orbit.TrueMass ) >= FLT_EPSILON ) {
Orbit.TruePeriod =
sqrt((Orbit.TrueDistance*Orbit.TrueDistance)
* Orbit.TrueDistance * 4.0 * (PI*PI)
/ (G_SI * Orbit.TrueMass));
}
else {
/* -------------------- default -------------------------------- */
Orbit.TruePeriod = 2.2226262e7; /* 0.7 year */
Orbit.TrueDistance = 1.496e11; /* one astronomical unit */
Orbit.TrueMass = 3.978e30; /* two solar mass */
/* M*G/(4*PI*PI) = D*D*D/P*P */
/* (Keplers Third Law) */
}
}
/******************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Print usage info
@param (void)
@return (void)
@heading Main
*******************************************************************/
void UsageInfo()
{
/* -------- PACKAGE, VERSION are declared in the makefile --- */
printf("\n");
/*@i@*/ printf(_(" Welcome to the %10s (version %5s)"), PACKAGE, VERSION);
printf(_(" Light Curve Synthesis Program\n"));
printf(_("\tCopyright (C) 1998 Rainer Wichmann\n"));
printf("\n");
printf(_(" Usage: nightfall <options (many)> <Q I RF1 RF2 T1 T2>\n"));
printf(_("\t (example: nightfall -G 0.8 85 0.8 1.0 5500 5800)\n"));
printf(_("\t (interactive: nightfall -U 0.8 85 0.8 1.0 5500 5800)\n"));
printf(_("\t (config_file: nightfall -U -C ty_boo.cfg)\n"));
printf("\n");
printf(_(" Basic System Parameters:\n"));
printf(_("\tQ = Mass(Secondary)/Mass(Primary),\n"));
printf(_("\tI = Inclination angle of orbital plane (degree)\n"));
printf(_("\tRF1, T1 = Primary: Roche_fill_factor, Temperature\n"));
printf(_("\tRF2, T2 = Secondary: Roche_fill_factor, Temperature\n"));
printf("\n");
printf(_(" Options:\n"));
printf(_("\t<-U> Interactive (GUI)\n"));
printf(_("\t<-A> Animated view\n"));
#ifdef _WITH_OPENGL
printf(_("\t<-g> Don't use OpenGL for animation\n"));
#endif
printf(_("\t<-C config_file> Input configuration file \n"));
printf(_("\t<-I datafile> Input observed data \n"));
printf(_("\t<-G[P,S,1,2]> Graph of lightcurve (default: 1)\n"));
printf(_("\t\tP,S: close-up of Primary/Secondary eclipse\n"));
printf(_("\t\t1,2: display 1/2 orbital cycles\n"));
printf(_("\t<-B[U/B/V...]> Bandpass to display in graph \n"));
printf("\t (default: V)\n");
printf(_("\t<-V[v,i,c,a]> Visualize geometry (default: v)\n"));
printf(_("\t\tv: view of stars, \n"));
printf(_("\t\ti: image of potential,\n"));
printf(_("\t\tc: contourmap of potential\n"));
printf(_("\t\ta: all of the above\n"));
printf(_("\t<-H> Hardcopy (postscript plot)\n"));
printf("\n");
printf(_(" Advanced System Parameters:\n"));
printf(_("\t<-O[P/S] value> log g (surface gravity)\n"));
printf(_("\t<-f[P/S] value> asynchroneous rotation (Omega/Omega_System)\n"));
printf(_("\t<-s[P/S] longitude latitude radius dimfactor> Spot on Primary/Secondary\n"));
printf(_("\t<-e eccentricity periastron_length> eccentric orbit parameter\n"));
printf(_("\t<-t[P/M/D] value> period/total mass/distance in days/solar masses/radii\n"));
printf("\n");
printf(_(" Debugging Options:\n"));
printf(_("\t<-D[vwb]> Debug [verbose,warn,busy]\n"));
/**********************************************
printf(" <-D[vwbdgflmrsecxiqoa]> Debug [verbose,warn,busy,\n");
printf(" divide,geometry,fractional\n");
printf(" lagrange,minfind,rootfind,spot,\n");
printf(" eclipsetest,computeflux,\n");
printf(" xcentric,infile\n");
printf(" chisquare,optimize,atmosphere]\n");
*********************************************/
printf("\n");
printf(_(" Data Fitting Options:\n"));
printf(_("\t<-X[fit_params] Tolerance>\n"));
printf(_("\t\t where 'fit_params' is a string for selecting\n"));
printf(_("\t\t the parameters to fit, as follows:\n"));
printf("\t\t0123456789: Q, i, RF1, RF2, T1, T2, e, w, F1, F2,\n");
printf(_("\t\tA-HI-P: 2 Spots (Primary), 2 Spots (Secondary),\n"));
printf(_("\t\tQR: Mass, Distance,\n"));
printf(_("\t\ta-l: Third Light (UBVRIJHKuvby) \n"));
printf(_("\t\t(set Tolerance = 0.001 to switch on Simulated Annealing)\n"));
printf(_("\t<-Y[as above] Step1 Step2> Chi Square Map (2 Parameters)\n"));
printf("\n");
printf(_(" Computation Options:\n"));
printf(_("\t<-Plamda_zero> Profile of absorption line at\n"));
printf(_("\t rest wavelength lamda_zero (nm)\n"));
printf(_("\t<-Nnn> nn steps for lightcurve \n"));
printf("\t (default 80)\n");
printf(_("\t<-M> use Model atmosphere \n"));
printf(_("\t<-F> compute Fractional visibility \n"));
printf(_("\t<-L[0-2]> Limb darkening method \n"));
printf("\t (default: linear == 0)\n");
printf(_("\t\t0: linear, 1: quadratic, 2: square root, 3: detailed\n"));
printf(_("\t<-R[1-9]> Reflection treatment \n"));
printf(_("\t (default: sloppy == 0)\n"));
printf(_("\t\t0: point source, 1-9: exact, up to 9 iterations\n"));
printf(_("\t<-a[P/S] albedo> Albedo des Primary/Secondary\n"));
printf(_("\t<-3CM> Third light, C: colour code, M: relative brightness \n"));
printf("\n");
printf(_("Compiled with:\n"));
printf(_("\t Support for interactive usage: "));
#ifdef _WITH_GTK
printf(_("yes\n"));
#else
printf(_("no\n"));
#endif
printf(_("\t Support for graphics: \t\t"));
#ifdef _WITH_PGPLOT
#ifdef _WITH_GNUPLOT
printf(_("gnuplot\n"));
#else
printf(_("pgplot\n"));
#endif
#else
printf(_("none\n"));
#endif
printf(_("\t Support for OPENGL: \t\t"));
#ifdef _WITH_OPENGL
printf(_("yes\n"));
#else
printf(_("no\n"));
#endif
printf(_("\t Support for MPI: \t\t"));
#ifdef _WITH_MPI
printf(_("yes\n"));
#else
printf(_("no\n"));
#endif
}
/******************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Initialize
@param (void)
@return (void)
@heading Main
*******************************************************************/
void InitFlags()
{
int j;
Flags.Step[0] = 0.01;
Flags.Step[1] = 0.01;
#ifdef _WITH_OPENGL
Flags.use_opengl = ON;
Flags.labels = OFF;
#endif
Flags.visualize = OFF;
for (j = 0; j < 128; ++j) Flags.debug[j] = OFF;
Flags.interactive = OFF;
Flags.animate = OFF;
/* added disk Flag MK 13.05.01 */
Flags.disk = OFF;
Flags.warp = OFF;
Flags.plot = OFF;
Flags.PlotBand = -1; /* signals that nothing read yet */
for (j = 0; j < (NUM_MAG+2); ++j) Flags.Passbands[j] = OFF;
for (j = 0; j < 64; ++j) Flags.simplex[j] = OFF;
Flags.fractional = OFF;
Flags.lineprofile = OFF;
Flags.WantFit = OFF;
Flags.SimplexTol = SIMPLEX_FIT_TOLERANCE;
Flags.anneal = OFF;
Flags.limb = 0;
Flags.reflect = 0;
Flags.Spots1 = 0;
Flags.Spots2 = 0;
Flags.fill = OFF;
Flags.elliptic = OFF;
Flags.asynchron1 = OFF;
Flags.asynchron2 = OFF;
Flags.blackbody = ON;
Flags.GL = OFF;
#ifdef _WITH_OPENGL
Flags.axes = OFF;
Flags.movie = OFF;
Flags.frame = 0;
Flags.wireframe = ON;
Flags.points = OFF;
Flags.texture = OFF;
Flags.textype = IMAGE;
Flags.GLdistancez = 1.5;
Flags.GLdistancex = 0.0;
#endif
Flags.first_pass = ON;
Flags.ProComputed = OFF;
Flags.ConfFile[0] = '\0';
strcpy (Orbit.Name, "(void)");
Orbit.Dist = 1.0;
Orbit.MinDist = 1.0;
Orbit.Excentricity = 0.0;
Orbit.Omega = 90.;
Orbit.TruePeriod = 0.0;
Orbit.TrueDistance = 0.0;
Orbit.TrueMass = 0.0;
Orbit.LambdaZero = LAMDAZERO_DEFAULT; /* should be nanometer */
for (j = 0; j < NUM_MAG; ++j) Orbit.Third[j] = 0.0;
Binary[Primary].Fratio = 1.0;
Binary[Secondary].Fratio = 1.0;
/* see Light.h for more details */
#ifdef HAVE_DISK
Binary[Disk].Tilt = 0.0;
Binary[Disk].Warp = 0.0;
Binary[Disk].Thick = 0.01;
Binary[Disk].HR = 0.0; /* 0 for a flat disk */
Binary[Disk].Rout = 0.8;
Binary[Disk].Rin = 0.2;
Binary[Disk].Temperature = 8000.0;
Binary[Disk].tempHS = Binary[Disk].Temperature;
Binary[Disk].longitudeHS = 0.0;
Binary[Disk].extentHS = 0.0;
Binary[Disk].depthHS = 0.0;
#endif
Flags.plotOpened = OFF;
strcpy(Out_Plot_File, "Nightfall.ps");
if (getenv("NIGHTFALL_PLOTFILE") != NULL)
{
strncpy( Out_Plot_File,
getenv("NIGHTFALL_PLOTFILE"),
(size_t) ((long) sizeof(Out_Plot_File) - 1) );
Out_Plot_File[sizeof(Out_Plot_File) - 1] = '\0';
}
return;
}
syntax highlighted by Code2HTML, v. 0.9.1