/* 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 <stdlib.h>
#include <string.h>
#include <float.h>
#ifdef HAVE_SYS_TIMERS_H
#include <sys/timers.h>
#endif
#include <signal.h>
#include "Light.h"
#ifdef _WITH_GTK
#include <gtk/gtk.h>
#endif
/* for _exit()
*/
#ifdef HAVE_UNISTD_H
#include <unistd.h>
#endif
/*********************************/
/* GLOBAL VARIABLES */
/*********************************/
/* include error messages */
#include "LightMsg.h"
char * program_invocation_name = NULL;
/* Flags */
FlagsHandle Flags;
/* Input Data */
Photo3 *Observed[NUM_MAG+2]; /* [colour][datum] */
/*@null@*/ FileType *FileList = NULL; /* data file list */
/* the profile */
/*@out@*/ /* float **ProData; */
float ProData[PHASESTEPS][PROFILE_ARRAY+1];
/* Flux */
/*@out@*/ PhotoOut *FluxOut; /* Output Flux */
double L90[3][NUM_MAG]; /* Luminosities */
double F90[NUM_MAG]; /* Normalization Variable */
int Umag = 0, Bmag = 1, Vmag = 2;
int Rmag = 3, Imag = 4;
int Jmag = 5, Hmag = 6, Kmag = 7;
int umag = 8, vmag = 9, bmag =10, ymag =11;
int mag1 = 12, mag2 = 13;
char Filters[NUM_MAG+2][24] =
{ "U Band", "B Band", "V Band", "R Band",
"I Band", "J Band", "H Band", "K Band", "u Band",
"v Band", "b Band", "y Band",
"Primary RV", "Secondary RV" };
/* Limb Darkening Coefficients */
double Limb[3][NUM_MAG][2]; /* component, passband, number */
double WaveL[3][NUM_MAG]; /* wavelengths */
/* Components */
const int Primary = 0; /* literal definition */
const int Secondary = 1; /* literal definition */
const int Disk = 2; /* literal definition */
BinaryComponent Binary[NUM_COMP]; /* properties of the components */
SurfaceElement *Surface[NUM_COMP]; /* surface elements of stars+Disk */
SpotLight *Spot[2]; /* surface spots */
OrbitType Orbit; /* orbit data */
/* These are Global Variables that are used (mainly) to pass */
/* Arguments to Functions that are called by RootFind/MinFind */
double F; /* nonsync rotation w/wsystem */
double Mq; /* mass ratio, used in Roche functions */
double RochePot; /* Roche Potential, used in Roche functions */
double PhaseScaledVolume; /* scaled volume */
double lambda, mu, nu; /* coordinates */
/* else */
char Out_Plot_File[256]; /* output plot file */
int StepsPerPi; /* No of steps per Pi (surface) */
double MaximumElements; /* max. no of surface elements */
int PhaseSteps; /* phasesteps for lightcurve */
double GArgSquare; /* for the SQR macro */
/* MPI */
int myrank = 0;
int numprocs = 1;
char processor_name[MPI_MAX_PROCESSOR_NAME] = "Unknown Processor";
char error_name[MPI_MAX_ERROR_STRING] = "No error";
void NF_Mpi_wait(void);
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short The signal handler
@param (int) sig The signal
@heading Main
****************************************************************************/
#if defined(_WITH_MPI) || defined(_WITH_MPI_COARSE)
void SIGhandler(int sig)
{
signal(sig, SIG_IGN);
#ifdef NF_MPI_DEBUG
fprintf(stderr, "MPI: Process %d of %d on %s: Killed by signal %d\n",
myrank, numprocs, processor_name, sig);
#endif
NF_MPI_Abort();
exit(EXIT_FAILURE);
}
#endif
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short The 'main' program
@param (int) argc The number of aguments
@param (char) *argv[] The argument list
@return (int) status The exit status
@heading Main
****************************************************************************/
int main(int argc, char *argv[])
{
double Merit = 0; /* the merit for (observed-computed) */
int i; /* a loop variable */
#ifdef _WITH_GTK
char * args[24];
int argcount = 0;
#endif
#if defined(_WITH_MPI) || defined(_WITH_MPI_COARSE)
struct sigaction act, oldact; /* signal handling */
double starttime, endtime;
MPI_Init(&argc,&argv);
/* set up a signal handler so we can clean up MPI if killed */
act.sa_handler = &SIGhandler; /* signal action */
sigemptyset( &act.sa_mask ); /* set an empty mask */
act.sa_flags = 0; /* init sa_flags */
#ifdef SIGHUP
sigaction(SIGHUP, &act, &oldact);
#endif
#ifdef SIGINT
sigaction(SIGINT, &act, &oldact);
#endif
#ifdef SIGQUIT
sigaction(SIGQUIT, &act, &oldact);
#endif
#ifdef SIGILL
sigaction(SIGILL, &act, &oldact);
#endif
#ifdef SIGABRT
sigaction(SIGABRT, &act, &oldact);
#endif
#ifdef SIGSEGV
sigaction(SIGSEGV, &act, &oldact);
#endif
#ifdef SIGPIPE
sigaction(SIGPIPE, &act, &oldact);
#endif
#ifdef SIGTRAP
sigaction(SIGTRAP, &act, &oldact);
#endif
#ifdef SIGIOT
sigaction(SIGIOT, &act, &oldact);
#endif
#ifdef SIGTERM
sigaction(SIGTERM, &act, &oldact);
#endif
#ifdef SIGBUS
sigaction(SIGBUS, &act, &oldact);
#endif
#ifdef SIGIO
sigaction(SIGIO, &act, &oldact);
#endif
#ifdef SIGPOLL
sigaction(SIGPOLL, &act, &oldact);
#endif
#ifdef SIGXCPU
sigaction(SIGXCPU, &act, &oldact);
#endif
#ifdef SIGPWR
sigaction(SIGPWR, &act, &oldact);
#endif
MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
MPI_Get_processor_name(processor_name, &i);
starttime = MPI_Wtime();
#endif
#ifdef NF_MPI_DEBUG
fprintf(stderr, "MPI: Process %d of %d on %s: Started\n",
myrank, numprocs, processor_name);
#endif
#ifdef ENABLE_NLS
/* Need this because setlocale (LC_NUMERIC, "C"); does not work, at least
* not with glibc 2.2 on linux.
*/
#ifdef HAVE_PUTENV
putenv("LC_NUMERIC=C");
#endif
#ifdef LC_ALL
/* for some stupid reason, german umlauts do not work
* with setlocale (LC_MESSAGES, "");
* but setlocale (LC_ALL, "");
* is fine (http://www.debianforum.de/forum/viewtopic.php?p=166242)
*/
/* setlocale (LC_MESSAGES, ""); */
setlocale (LC_ALL, "");
#else
#error NLS enabled, but LC_ALL undefined
#error you may want to use --with-included-gettext or --disable-nls
#endif
/* Set this to 'C' in order not to run into problems with the datafiles.
* UPDATE: Contrary to the documentation on the 'setlocale' function,
* this does not work at all, forcing us to fall back to the crude
* 'putenv' hack above.
*/
#ifdef LC_NUMERIC
setlocale (LC_NUMERIC, "C");
#else
#error NLS enabled, but LC_NUMERIC undefined
#error you may want to use --with-included-gettext or --disable-nls
#endif
bindtextdomain (PACKAGE,
(getenv("NIGHTFALL_LOCALE_DIR") == NULL) ?
LOCALEDIR : getenv("NIGHTFALL_LOCALE_DIR"));
textdomain (PACKAGE);
#endif
/* >>>>>>>>>>>>>>>> ALLOCATE MEMORY <<<<<<<<<<<<<<<<<<<<< */
#if 0
program_invocation_name = malloc (1 + strlen(argv[0]));
strcpy(program_invocation_name, argv[0]);
#endif
/* should be freed by OS on program exit (?) */
Surface[0] = malloc(MAXELEMENTS * sizeof(SurfaceElement));
Surface[1] = malloc(MAXELEMENTS * sizeof(SurfaceElement));
#ifdef HAVE_DISK
Surface[2] = malloc(MAXELEMENTS * sizeof(SurfaceElement));
#endif
Spot[0] = malloc(N_SPOT * sizeof(SpotLight));
Spot[1] = malloc(N_SPOT * sizeof(SpotLight));
FluxOut = malloc(PHASESTEPS * sizeof(PhotoOut));
/* Is a fixed-size array */
/* ProData = matrix(PHASESTEPS,PROFILE_ARRAY+1); */
if (Surface[0] == NULL || Surface[1] == NULL
#ifdef HAVE_DISK
|| Surface[2]== NULL
#endif
|| Spot[0] == NULL || Spot[1] == NULL
|| FluxOut == NULL
/* || ProData == NULL */)
nf_error(_(errmsg[0]));
for (i=0; i < NUM_MAG+2; ++i) {
Observed[i] = malloc(MAXOBS * sizeof(Photo3));
if (Observed[i] == NULL)
nf_error(_(errmsg[0]));
}
#if defined(_WITH_MPI_COARSE)
if (myrank > 0) {
NF_Mpi_wait();
}
#endif
/* >>>>>>>>>>>>>>>> USAGE INFO IF NO INPUT <<<<<<<<<<<<<<<<<<<<< */
if (argc == 1) {
if (myrank == 0) {
UsageInfo();
MPI_Abort(MPI_COMM_WORLD, EXIT_FAILURE);
exit(EXIT_FAILURE);
} else {
MPI_Finalize();
exit(EXIT_FAILURE);
}
}
#ifdef _WITH_GTK
args[0] = argv[0]; argcount = 1;
for (i = 1; i < argc; ++i) {
if (argcount < 24) {
if ((argv[i][0] == '-') && (argv[i][1] == '-') ) {
args[argcount] = argv[i];
++argcount;
if ((i+1) < argc) {
args[argcount] = argv[i+1];
++argcount;
}
}
}
}
#endif
#ifdef _WITH_PGPLOT
#ifndef _WITH_GNUPLOT
SetPgplotWindow(); /* set the PGPLOT window size */
#endif
#endif
/* >>>>>>>>>>>>>>>> INITIALIZE <<<<<<<<<<<<<<<<<<<<< */
InitFlags();
Flags.smap = SetMapPath ();
if (Flags.smap == ON)
SetMapBand ();
Flags.monochrome = DefineMonoWavelengths ();
/* >>>>>>>>>>>>>>>> GET ARGUMENTS <<<<<<<<<<<<<<<<<<<<<< */
Flags.parseCL = ON; /* we parse the command line */
GetArguments(&argc, argv);
Flags.parseCL = OFF; /* finished parsing cl */
/* >>>>>>>>>>>>>>>> MAIN SECTION <<<<<<<<<<<<<<<<<<<<< */
if (Flags.interactive == OFF && myrank == 0) {
if (Flags.WantMap == ON) (void) ChiMap();
else if (Flags.WantFit == ON) (void) Simplex();
else (void) MainLoop(&Merit);
WriteOutput();
} else if (myrank == 0) {
#ifdef _WITH_GTK
/******************************
fprintf(stderr,
"%f %f %f %f %f %f\n",
Orbit.Inclination,
Binary[Primary].Mq,
Binary[Primary].Temperature,
Binary[Secondary].Temperature,
Binary[Primary].RocheFill,
Binary[Secondary].RocheFill
);
*******************************/
/* some reasonable defaults */
if (Orbit.Inclination <= FLT_EPSILON)
Orbit.Inclination = DTOR * 80.0;
if (Binary[Primary].Mq <=
(Flags.fill == ON ? LIM_MQO_L : LIM_MQ_L))
{
Binary[Primary].Mq = 1.0;
Binary[Secondary].Mq = 1.0;
}
if (Binary[Primary].Temperature <= FLT_EPSILON)
Binary[Primary].Temperature = 5700.0;
if (Binary[Secondary].Temperature <= FLT_EPSILON)
Binary[Secondary].Temperature = 5700.0;
if (Binary[Primary].RocheFill <= FLT_EPSILON)
Binary[Primary].RocheFill = 0.6;
if (Binary[Secondary].RocheFill <= FLT_EPSILON)
Binary[Secondary].RocheFill = 0.6;
if (Binary[Primary].log_g <= FLT_EPSILON)
Binary[Primary].log_g = 4.0;
if (Binary[Secondary].log_g <= FLT_EPSILON)
Binary[Secondary].log_g = 4.0;
/* start the GUI */
the_gui(argcount, args);
#else
nf_error (_(errmsg[20]));
#endif
}
else if (myrank > 0) {
NF_Mpi_wait();
}
/* >>>>>>>>>>>>>>>> END <<<<<<<<<<<<<<<<<<<<<< */
/* successful exit */
#ifdef _WITH_GNUPLOT
if (Flags.plotOpened == ON && myrank == 0)
{
Flags.plotOpened = OFF; cpgend();
}
#endif
#if defined(_WITH_MPI) || defined(_WITH_MPI_COARSE)
endtime = MPI_Wtime();
if (Flags.debug[BUSY] == ON || Flags.debug[VERBOSE] == ON)
printf("Elapsed time: %f sec\n", endtime-starttime);
#endif
NF_Mpi_call(NF_TAG_END);
MPI_Finalize();
#ifdef NF_MPI_DEBUG
fprintf(stderr, "MPI: Process %d of %d on %s: Exiting\n",
myrank, numprocs, processor_name);
#endif
#ifdef _WITH_GTK
if (Flags.interactive == ON && myrank == 0)
{
#ifdef _WITH_OPENGL
/* 2002-05-03 rwichmann: very dirty workaround for segfault on exit(),
* probably in some atexit() cleanup function.
*/
_exit (0);
#else
gtk_exit(0);
#endif
}
#endif
return (0);
}
/*********************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.37
@short Loop for computation of lightcurve
@param (double) *Merit The merit for (observed-computed)
@return (int) status The exit status
@heading Main Loop
**********************************************************************/
int MainLoop(/*@out@*/ double *Merit)
{
int PhaseIndex; /* the index for the phase */
int ErrCode; /* exit status of subroutines */
long i; /* loop variable */
SurfaceElement *SurfPtrP; /* pointer to surface primary */
SurfaceElement *SurfPtrS; /* pointer to surface secondary */
#ifdef HAVE_DISK /* Changed by MPR */
SurfaceElement *SurfPtrD; /* pointer to surface disk */
#endif
#ifdef _WITH_GTK
float pvalue; /* progress bar increment */
char message[64]; /* statusbar message */
#endif
/* >>>>>>>>>>>>> Initialize <<<<<<<<<<<<<< */
#if defined(_WITH_MPI_COARSE)
int loop_alloc[PHASESTEPS];
int k = 0;
if (myrank == 0) {
NF_Mpi_call(NF_TAG_MAINLOOP);
}
NF_Mpi_distribute();
for (i = 0; i < PhaseSteps; ++i)
{
loop_alloc[i] = k; ++k;
if (k == numprocs)
k = 0;
}
#endif
SurfPtrP = Surface[Primary];
SurfPtrS = Surface[Secondary];
#ifdef HAVE_DISK /* Changed by MPR */
SurfPtrD = Surface[Disk];
#endif
#ifdef _WITH_GTK
pvalue = 0.0;
#endif
*Merit = 0.0;
for (i = 0; i < MAXELEMENTS; ++i) {
SurfPtrP->SpotNum = 0;
SurfPtrS->SpotNum = 0;
SurfPtrP->SumDimFactor = 0;
SurfPtrS->SumDimFactor = 0;
++SurfPtrP;
++SurfPtrS;
#ifdef HAVE_DISK /* Changed by MPR */
SurfPtrD->SpotNum = 0;
SurfPtrD->SumDimFactor = 0;
++SurfPtrD;
#endif
}
/* >>>>>>>>>>>>> LightGeometry <<<<<<<<<<<<<< */
#ifdef _WITH_GTK
if (Flags.interactive == ON)
my_appbar_push(_("Setup geometry") );
#endif
/* first pass (initial setup of stars) */
Flags.first_pass = ON;
#ifdef _WITH_OPENGL
/* reset frame counter for GL movie */
if (Flags.movie == ON && Flags.use_opengl == ON) {
Flags.frame=0;
}
#endif
/* store Roche fill factor -- needed for eccentric orbit */
/* (at final phase not equal to initial value) */
/* we don't use eccentric orbits plus a disk at the moment */
/* -> nothing to do here for the disk */
Binary[Primary].RocheStore = Binary[Primary].RocheFill;
Binary[Secondary].RocheStore = Binary[Secondary].RocheFill;
/* elliptic orbit */
if (Flags.elliptic == ON) FindPeriastron();
/* Primary at x=y=z=0, Secondary at x=1, y=z=0 */
if (Flags.fill == OFF) {
/* -- no overcontact -- */
/* THE PRIMARY */
ErrCode = DefineParam(Primary);
if (ErrCode == 8) return(8);
if (Flags.debug[VERBOSE] == ON && myrank == 0)
{
printf("\n");
printf(_(" -- Primary defined --\n") );
}
/* THE SECONDARY */
ErrCode = DefineParam(Secondary);
if (ErrCode == 8) return(8);
if (Flags.debug[VERBOSE] == ON && myrank == 0)
{
printf("\n");
printf(_(" -- Secondary defined --\n") );
}
/* THE DISK */
#ifdef HAVE_DISK
if (Flags.disk == ON) {
/* new code for the disk parameters TBD */
ErrCode = DefineParamDisk();
if (ErrCode == 8) return(8);
if (Flags.debug[VERBOSE] == ON && myrank == 0)
{
printf("\n");
printf(_(" -- Disk defined --\n") );
}
}
#endif
} else if (Flags.disk == ON) {
/* Disk + overcontact is not possible */
#ifdef _WITH_GTK
if (Flags.interactive == OFF) {
nf_error(_(" -- Overcontact+Disk is not possible --"));
} else {
make_dialog(_(" -- Overcontact+Disk is not possible --"));
return (8);
}
#else
nf_error(_(" -- Overcontact+Disk is not possible --"));
#endif
}
else {
/* -- Overcontact system -- */
ErrCode = DefineParamOver();
if (ErrCode == 8) return(8);
if (Flags.debug[VERBOSE] == ON && myrank == 0)
{
printf("\n");
printf(_(" -- Overcontact Geometry defined --\n"));
}
}
/* >>>>>>>>>>>>> LightLimbDarkCoeff <<<<<<<<<<<<<<<<<<< */
/* Compute effective wavelength for filters */
EffectiveWavelength(Primary);
EffectiveWavelength(Secondary);
/* works fine for disk, too. checked MK */
#ifdef HAVE_DISK
if (Flags.disk == ON) EffectiveWavelength(Disk);
#endif
/* >>>>>>>>>>>>> LightDivide <<<<<<<<<<<<<<<<<<< */
/* ----------- DO SECONDARY SECOND !!! ----------------------- */
/* set up the surface grid */
/* THE PRIMARY */
ErrCode = DivideSurface(Primary);
if (ErrCode == 8) return(8);
if (Flags.debug[VERBOSE] == ON && myrank == 0)
{ printf("\n"); printf(_(" Primary divided\n") );}
/* THE SECONDARY */
ErrCode = DivideSurface(Secondary);
if (ErrCode == 8) return(8);
if (Flags.debug[VERBOSE] == ON && myrank == 0)
{ printf("\n"); printf(_(" Secondary divided\n") );}
#ifdef HAVE_DISK /* Changed by MPR */
/* THE DISK */
if (Flags.disk == ON) {
ErrCode = DivideDisk(Disk);
if (ErrCode == 8) return(8);
if (Flags.debug[VERBOSE] == ON && myrank == 0)
{ printf("\n"); printf(_(" Disk divided\n") );}
}
#endif
/* >>>>>>>>>>>>> LightLimbDarkCoeff <<<<<<<<<<<<<<<<<<< */
/* compute limb darkening coefficients */
LightLimbDark(Primary);
LightLimbDark(Secondary);
#ifdef HAVE_DISK
if (Flags.disk == ON) {
LightLimbDark(Disk);
}
#endif
/* >>>>>>>>>>>>> LightSynthesize <<<<<<<<<<<<<<<<<< */
/* Calculate various quantities for geometric eclipse tests */
/* if Asynch Rotation < 1.0, eventually the stars intersect */
/* we test for this and abort */
if (LightSetupTests() > 0) {
#ifdef _WITH_GTK
if (Flags.interactive == OFF) {
nf_error(_(errmsg[5]));
} else {
make_dialog(_(errmsg[5]));
return (8);
}
#else
nf_error(_(errmsg[5]));
#endif
}
Flags.first_pass = OFF;
/* --------------- LightVisualize ------------------------- */
#ifdef _WITH_PGPLOT
if (Flags.visualize > OFF && myrank == 0) PlotGeometry();
#endif
/* ------------- LightOutput ------------------------ */
/* Initialize flux array (total flux summed over all components)
*/
InitOutFlux();
/* >>>>>>>>>>>>> LightSynthesize <<<<<<<<<<<<<<<<<< */
if (Flags.debug[VERBOSE] == ON && myrank == 0)
{
printf("\n");
printf(_(" Starting Lightcurve Synthesis\n") );
}
/* ************** LOOP OVER PHASE ***************************** */
if (Flags.debug[BUSY] == ON && myrank == 0) {
printf(" 0 per cent");
(void) fflush(stdout);
}
#ifdef _WITH_GTK
if (Flags.interactive == ON)
my_appbar_push(_("Start lightcurve computation"));
#endif
for (PhaseIndex = 0; PhaseIndex < PhaseSteps; ++PhaseIndex) {
#if defined(_WITH_MPI_COARSE)
if (loop_alloc[PhaseIndex] != myrank)
continue;
#endif
#ifdef _WITH_GTK
/* update the progress bar */
pvalue = (float)(PhaseIndex)/(float)(PhaseSteps);
pvalue = (pvalue > 1.0 ? 1.0 : pvalue);
if (Flags.interactive == ON) {
#ifdef USING_GTK2
gtk_progress_bar_set_fraction (GTK_PROGRESS_BAR (progress_bar), pvalue);
#else
gtk_progress_bar_update (GTK_PROGRESS_BAR (progress_bar), pvalue);
#endif
while (gtk_events_pending()) gtk_main_iteration();
}
#endif
/* a quick hack for the MakeSpots program after shifting the */
/* call to LightDivide */
/* which was necessary to include spots correctly in the */
/* reflection treatment */
Orbit.PhaseIndex = PhaseIndex;
if (Flags.debug[BUSY] == ON) {
printf("\b\b\b\b\b\b\b\b\b\b\b\b%3d per cent",
(int)((100.0*PhaseIndex)/PhaseSteps));
(void) fflush(stdout);
}
/* ------------- LightElliptic ------------------------ */
if (Flags.elliptic == ON) {
/* update (nearly) everything */
/* disk and elliptic are mutually exclusive */
SolveKepler(PhaseIndex);
ErrCode = UpdateGeometry(Primary);
if (ErrCode == 8) return(8);
ErrCode = UpdateGeometry(Secondary);
if (ErrCode == 8) return(8);
/* DO SECONDARY SECOND !!! */
ErrCode = DivideSurface(Primary);
if (ErrCode == 8) return(8);
ErrCode = DivideSurface(Secondary);
if (ErrCode == 8) return(8);
ErrCode = LightSetupTests();
if (ErrCode > 0) return(8);
}
/* ----------- move spots if asynchroneous rotation ------------ */
if ( (Flags.Spots1 > OFF)
&& Flags.asynchron1 == ON
&& Flags.elliptic == OFF )
{
SurfPtrP = Surface[Primary];
for (i = 0; i < Binary[Primary].NumElem; ++i) {
/* Gravitational Darkening */
SurfPtrP->temp =
Binary[Primary].Temperature
* (exp(Binary[Primary].GravDarkCoeff
* log(SurfPtrP->grav))
/ Binary[Primary].DarkeningGrav );
++SurfPtrP;
}
MakeSpots(Primary, Orbit.PhaseIndex);
}
if ( (Flags.Spots2 > OFF)
&& Flags.asynchron1 == ON
&& Flags.elliptic == OFF )
{
SurfPtrS = Surface[Secondary];
for (i = 0; i < Binary[Secondary].NumElem; ++i) {
/* Gravitational Darkening */
SurfPtrS->temp =
Binary[Secondary].Temperature
* (exp(Binary[Secondary].GravDarkCoeff
* log(SurfPtrS->grav))
/ Binary[Secondary].DarkeningGrav );
++SurfPtrP;
}
MakeSpots(Secondary, Orbit.PhaseIndex);
}
/* >>>>>>>>>>>>> LightReflect <<<<<<<<<<<<<<<< */
/* TBD reflection on disk */
if (Flags.elliptic == OFF &&
( (Flags.Spots2 > OFF
&& Flags.asynchron1 == ON) ||
(Flags.Spots1 > OFF
&& Flags.asynchron2 == ON) ) ){
if (Flags.reflect > 0) {
LightReflect();
}
if (Flags.reflect == 0) {
SimpleReflect(Primary);
SimpleReflect(Secondary);
}
}
/* >>>>>>>>>>>>> LightSynthesize <<<<<<<<<<<<<<<< */
if (Flags.fill == ON) LightCopyThroat();
/* eclipse checking */
LightCurve(PhaseIndex);
#ifdef HAVE_DISK
if (Flags.disk == ON)
LightCurveDisk(PhaseIndex);
#endif
/* ------------- LightFractional --------------------- */
/* fractional visibility */
if (Flags.fractional == ON) {
if (Flags.InEclipse1 == ON) {
LightFractionalPot(Primary);
}
if (Flags.InEclipse2 == ON) {
LightFractionalPot(Secondary);
}
/* code for disk TBD */
}
/* >>>>>>>>>>>>> LightSynthesize <<<<<<<<<<<<<<<< */
if (Flags.fill == ON) LightCopyThroatBack();
/* add support for disk TBD */
RadVel(PhaseIndex);
/* light curve */
LightFlux(Primary, PhaseIndex);
/*
* Now check if we are at normalization phase,
* then fill L90[n][bands]
*/
LightLuminosity(Primary, PhaseIndex);
LightFlux(Secondary, PhaseIndex);
LightLuminosity(Secondary, PhaseIndex);
#ifdef HAVE_DISK
/* code for disk TBD */
if (Flags.disk == ON) {
LightFlux(Disk, PhaseIndex);
LightLuminosity(Disk, PhaseIndex);
}
#endif
/* >>>>>>>>>>>>> LightAnimate <<<<<<<<<<<<<<<<< */
/* if we have OpenGL, use it if supported by the display and if the
* user really wants it, otherwise use the conventional animation
*/
#ifdef _WITH_OPENGL
#ifdef _WITH_PGPLOT
if ((Flags.animate == ON)&&(Flags.use_opengl == OFF || Flags.GL == OFF)) {
Animate(PhaseIndex);
} else {
if (Flags.animate == ON)
AnimateGL(PhaseIndex);
}
#else
if ((Flags.animate == ON) && (Flags.GL == ON))
AnimateGL(PhaseIndex);
#endif
/* no OpenGL, use the conventional animation
*/
#else
#ifdef _WITH_PGPLOT
if (Flags.animate == ON)
Animate(PhaseIndex);
#endif
#endif
/* >>>>>>>>>>>>> LightMap <<<<<<<<<<<<<<<<< */
if (Flags.smap == ON)
/* add disk geometry TBD */
PrintSurfaceMap (PhaseIndex);
} /* END LOOP OVER PHASE */
#if defined(_WITH_MPI_COARSE)
NF_Mpi_collect();
if (myrank != 0) {
return (0);
} else {
NF_Mpi_endcall();
}
#endif
/* >>>>>>>>>>>>>>> LightLib <<<<<<<<<<<<<<<<< */
if (Flags.elliptic == ON) PhaseShift();
/* check for disk geometry */
NormalizeFlux();
/* >>>>>>>>>>>>>>> LightOptimize <<<<<<<<<<<<<<<<< */
/* compute the merit for (observed-computed) */
*Merit = MeritFunction();
if (Flags.debug[BUSY] == ON) {
if (*Merit >= FLT_EPSILON) {
printf("\b\b\b\b\b\b\b\b\b\b\b\b100 per cent, ChiSquare = %15.5f\n",
*Merit);
(void) fflush(stdout);
} else {
printf("\b\b\b\b\b\b\b\b\b\b\b\b100 per cent\n");
(void) fflush(stdout);
}
}
Flags.IsComputed = ON;
#ifdef _WITH_GTK
if (Flags.interactive == ON) {
sprintf(message, _("Chi Square = %15.5f"), *Merit);
my_appbar_push(message);
#ifdef USING_GTK2
gtk_progress_bar_set_fraction (GTK_PROGRESS_BAR (progress_bar), 1.0);
#else
gtk_progress_bar_update (GTK_PROGRESS_BAR (progress_bar), 1.0);
#endif
while (gtk_events_pending()) gtk_main_iteration();
}
#endif
#ifdef _WITH_PGPLOT
if (Flags.plot > OFF) PlotOutput();
#endif
/* get back Roche fill factor */
Binary[Primary].RocheFill = Binary[Primary].RocheStore;
/* use fill factor for disk, too ??? TBD */
Binary[Secondary].RocheFill = Binary[Secondary].RocheStore;
return(0);
}
/****************************************************************
*
*
* M P I S T U F F
*
****************************************************************/
#if defined(_WITH_MPI) || defined(_WITH_MPI_COARSE)
static int nf_mpi_called = 0;
void NF_MPI_Abort()
{
if (myrank == 0 && nf_mpi_called == 0)
{
NF_Mpi_call(NF_TAG_END);
/* Hangs in MPI_Finalize(); then aborted b/o
* p0_4640: (15.850686) net_recv failed for fd = 5
* p0_4640: p4_error: net_recv read, errno = : 104
*/
MPI_Finalize();
}
else
{
MPI_Abort(MPI_COMM_WORLD, EXIT_FAILURE);
}
return;
}
void NF_Mpi_endcall()
{
nf_mpi_called = 0;
return;
}
void NF_Mpi_call(int tag)
{
int i;
for (i = 1; i < numprocs; ++i)
MPI_Send(&i, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
nf_mpi_called = 1;
return;
}
void NF_Mpi_wait() {
MPI_Request request;
int flag;
MPI_Status status;
int buf;
int tag;
double Merit; /* dummy, only evaluated by process 0 */
do {
MPI_Irecv(&buf, 1, MPI_INT, 0 /* source */, MPI_ANY_TAG /* tag */,
MPI_COMM_WORLD, &request);
do {
MPI_Test(&request, &flag, &status);
} while (!flag);
tag = status.MPI_TAG;
switch (tag) {
case NF_TAG_END:
#ifdef NF_MPI_DEBUG
fprintf(stderr, "MPI: Process %d of %d on %s: Exiting\n",
myrank, numprocs, processor_name);
#endif
MPI_Finalize();
exit(EXIT_SUCCESS);
case NF_TAG_MAINLOOP:
MainLoop(&Merit);
break;
default:
#ifdef NF_MPI_DEBUG
fprintf(stderr,
"MPI: Process %d of %d on %s: Unexpected message - aborting\n",
myrank, numprocs, processor_name);
#endif
MPI_Abort(MPI_COMM_WORLD, EXIT_FAILURE);
exit(EXIT_FAILURE);
}
} while (1 == 1);
return;
}
void NF_Mpi_collect ()
{
int i, j = 0, k = 0, size;
double * msg_comp;
int loop_alloc[PHASESTEPS];
MPI_Status status;
if (numprocs < 2)
return;
for (i = 0; i < PhaseSteps; ++i)
{
loop_alloc[i] = k; ++k;
if (k == numprocs)
k = 0;
}
size = PhaseSteps * (NUM_MAG + 3 + 3 + 1);
msg_comp = malloc(size * sizeof(double));
if (myrank != (numprocs - 1))
{
MPI_Recv(msg_comp, size, MPI_DOUBLE, myrank + 1,
NF_MSG_TAG + 4, MPI_COMM_WORLD, &status);
}
j = 0;
if (myrank != 0)
{
for (i = 0; i < PhaseSteps; ++i)
{
if (loop_alloc[i] == myrank)
{
msg_comp[j++] = FluxOut[i].Phase;
for (k = 0; k < NUM_MAG; ++k)
msg_comp[j++] = FluxOut[i].Mag[k];
for (k = 0; k < 3; ++k)
msg_comp[j++] = FluxOut[i].RV[k];
for (k = 0; k < 3; ++k)
msg_comp[j++] = FluxOut[i].D_RV[k];
}
else
{
j += (NUM_MAG + 7);
}
}
MPI_Send(msg_comp, size, MPI_DOUBLE, myrank-1,
NF_MSG_TAG + 4, MPI_COMM_WORLD);
}
else
{
for (i = 0; i < PhaseSteps; ++i)
{
if (loop_alloc[i] != 0)
{
FluxOut[i].Phase = msg_comp[j++];
for (k = 0; k < NUM_MAG; ++k)
FluxOut[i].Mag[k] = msg_comp[j++];
for (k = 0; k < 3; ++k)
FluxOut[i].RV[k] = msg_comp[j++];
for (k = 0; k < 3; ++k)
FluxOut[i].D_RV[k] = msg_comp[j++];
}
else
{
j += (NUM_MAG + 7);
}
}
}
return;
}
void NF_Mpi_distribute ()
{
int i, j = 0, k, size;
double * msg_comp;
int * msg_int;
extern int MapBand;
extern char * MapPath;
extern float mono_zero[NUM_MAG];
int MapPathSize = 0;
if (MapPath != NULL)
MapPathSize = strlen(MapPath) + 1;
/* ---------- the integer values -------------------
*/
size = (256 + 12 + 24);
msg_int = malloc(size * sizeof(int));
if (!msg_int) {
nf_error(_(errmsg[0]));
}
j = 0;
for (k = 0; k < 256; ++k)
msg_int[j++] = (int) Flags.debug[k];
msg_int[j++] = PhaseSteps;
msg_int[j++] = Flags.lineprofile;
msg_int[j++] = Flags.blackbody;
msg_int[j++] = Flags.fractional;
msg_int[j++] = Flags.asynchron1;
msg_int[j++] = Flags.asynchron2;
msg_int[j++] = Flags.elliptic;
msg_int[j++] = Flags.limb;
msg_int[j++] = Flags.reflect;
msg_int[j++] = Flags.Spots1;
msg_int[j++] = Flags.Spots2;
msg_int[j++] = StepsPerPi;
msg_int[j++] = Flags.monochrome;
msg_int[j++] = Flags.smap;
msg_int[j++] = MapBand;
msg_int[j++] = MapPathSize;
MPI_Bcast(msg_int, size, MPI_INT, 0, MPI_COMM_WORLD );
j = 0;
for (k = 0; k < 256; ++k)
Flags.debug[k] = msg_int[j++];
PhaseSteps = msg_int[j++];
Flags.lineprofile = msg_int[j++];
Flags.blackbody = msg_int[j++];
Flags.fractional = msg_int[j++];
Flags.asynchron1 = msg_int[j++];
Flags.asynchron2 = msg_int[j++];
Flags.elliptic = msg_int[j++];
Flags.limb = msg_int[j++];
Flags.reflect = msg_int[j++];
Flags.Spots1 = msg_int[j++];
Flags.Spots2 = msg_int[j++];
StepsPerPi = msg_int[j++];
Flags.monochrome = msg_int[j++];
Flags.smap = msg_int[j++];
MapBand = msg_int[j++];
MapPathSize = msg_int[j++];
if (myrank != 0)
{
Flags.WantMap = OFF;
Flags.WantFit = OFF;
Flags.anneal = OFF;
Flags.interactive = OFF;
Flags.animate = OFF;
Flags.visualize = OFF;
Flags.plot = OFF;
#ifdef _WITH_OPENGL
Flags.use_opengl = OFF;
#endif
}
/* ---------- the char arrays ------------------------
*/
if (MapPathSize != 0)
{
if (myrank == 0)
{
MPI_Bcast(MapPath, MapPathSize, MPI_BYTE, 0, MPI_COMM_WORLD );
}
else
{
MapPath = (char *) malloc (MapPathSize);
if (MapPath == NULL)
{
nf_error(_(errmsg[0]));
exit(EXIT_FAILURE);
}
MPI_Bcast(MapPath, MapPathSize, MPI_BYTE, 0, MPI_COMM_WORLD );
}
}
/* ---------- the floating-point values ------------
*/
size = (12 + (4*2*N_SPOT) + 11 + (2 * NUM_MAG) + 24);
msg_comp = malloc(size * sizeof(double));
if (!msg_comp) {
nf_error(_(errmsg[0]));
}
j = 0;
for (i = 0; i < 2; ++i)
{
msg_comp[j++] = Binary[i].Mq;
msg_comp[j++] = Binary[i].RocheFill;
msg_comp[j++] = Binary[i].Temperature;
msg_comp[j++] = Binary[i].Albedo;
msg_comp[j++] = Binary[i].GravDarkCoeff;
msg_comp[j++] = Binary[i].Fratio;
for (k = 0; k < N_SPOT; ++k)
{
msg_comp[j++] = Spot[i][k].longitude;
msg_comp[j++] = Spot[i][k].latitude;
msg_comp[j++] = Spot[i][k].radius;
msg_comp[j++] = Spot[i][k].dimfactor;
}
}
msg_comp[j++] = Orbit.Inclination;
msg_comp[j++] = Orbit.TrueDistance;
msg_comp[j++] = Orbit.TruePeriod;
msg_comp[j++] = Orbit.TrueMass;
msg_comp[j++] = Orbit.Excentricity;
msg_comp[j++] = Orbit.Omega;
msg_comp[j++] = Orbit.LambdaZero;
msg_comp[j++] = Orbit.Dist;
msg_comp[j++] = Flags.Step[0];
msg_comp[j++] = Flags.Step[1];
msg_comp[j++] = MaximumElements;
for (k = 0; k < NUM_MAG; ++k)
msg_comp[j++] = Orbit.Third[k];
for (k = 0; k < NUM_MAG; ++k)
msg_comp[j++] = mono_zero[k];
MPI_Bcast(msg_comp, size, MPI_DOUBLE, 0, MPI_COMM_WORLD );
j = 0;
for (i = 0; i < 2; ++i)
{
Binary[i].Mq = msg_comp[j++];
Binary[i].RocheFill = msg_comp[j++];
Binary[i].Temperature = msg_comp[j++];
Binary[i].Albedo = msg_comp[j++];
Binary[i].GravDarkCoeff = msg_comp[j++];
Binary[i].Fratio = msg_comp[j++];
for (k = 0; k < N_SPOT; ++k)
{
Spot[i][k].longitude = msg_comp[j++];
Spot[i][k].latitude = msg_comp[j++];
Spot[i][k].radius = msg_comp[j++];
Spot[i][k].dimfactor = msg_comp[j++];
}
}
Orbit.Inclination = msg_comp[j++];
Orbit.TrueDistance = msg_comp[j++];
Orbit.TruePeriod = msg_comp[j++];
Orbit.TrueMass = msg_comp[j++];
Orbit.Excentricity = msg_comp[j++];
Orbit.Omega = msg_comp[j++];
Orbit.LambdaZero = msg_comp[j++];
Orbit.Dist = msg_comp[j++];
Flags.Step[0] = msg_comp[j++];
Flags.Step[1] = msg_comp[j++];
MaximumElements = msg_comp[j++];
for (k = 0; k < NUM_MAG; ++k)
Orbit.Third[k] = msg_comp[j++];
for (k = 0; k < NUM_MAG; ++k)
mono_zero[k] = msg_comp[j++];
return;
}
#else
/* ---------- NO MPI -----------
*/
void NF_Mpi_wait() { return; }
void NF_Mpi_call(int tag) { return; }
void NF_Mpi_endcall() { return; }
void NF_MPI_Abort() { return; }
#endif
syntax highlighted by Code2HTML, v. 0.9.1