/* 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 #ifdef HAVE_SYS_TIMERS_H #include #endif #include #include "Light.h" #ifdef _WITH_GTK #include #endif /* for _exit() */ #ifdef HAVE_UNISTD_H #include #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