/* 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. */ /* Define this if you want to print coordinates */ /* #define PRINT_COORDS 1 */ #include #include #include #include #include #include "Light.h" #ifdef _WITH_PGPLOT #ifndef _WITH_GNUPLOT #include "cpgplot.h" #endif #endif #ifdef _WITH_PGPLOT /****************************************************************** @package nightfall @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de) @version 1.0 @short Real-time animation of lightcurve computation @param (int) j The phase index @return (void) @heading Animation *******************************************************************/ void Animate(int j) { long i; /* loop variable */ int test; /* test variable */ int Comp; /* the stellar component */ int MaxComp = 2; /* the maximum # of components */ long CountElem; /* # of elements in X/Yplot */ float *Xplot; /* plot array */ float *Yplot; /* plot array */ float *Zplot; /* plot array */ double sinAlpha; /* sin inclination */ double cosAlpha; /* cos inclination */ double sinPhi; /* sin phase */ double cosPhi; /* cos phase */ double Trans[4][4]; /* transformation matrix */ double t10, t11, t12, t13; /* transformation matrix */ double t20, t21, t22, t23; /* transformation matrix */ double Calib; /* normalization for flux */ double C; /* sign (for primary/secondary) */ float x1, x2, y1, y2; /* plot window borders */ float MinMag, MaxMag; /* plot scaling */ float MinVel, MaxVel; /* plot scaling */ float Upper; /* plot window borders */ float Xtick; /* tickmark distance */ PhotoOut *FluxPtr = FluxOut; /* pointer to flux */ SurfaceElement *SurfPtr; /* pointer to surface */ OrbitType *OrbPtr = &Orbit; /* pointer to Orbit */ char titlestring[MAX_CFG_INLINE+1]; /* the title string */ float com; /* centre of mass */ int js; /* j minus invalid */ #ifdef HAVE_DISK static float init_phase; /* initial phase */ #endif #if defined(PRINT_COORDS) static FILE * outfile; double sx[3], sy[3]; double ssin, scos, posw, hyp; #endif /* >>>>>>>>>>> Initialize <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */ Xplot = malloc(3 * MAXELEMENTS * sizeof(float)); Yplot = malloc(3 * MAXELEMENTS * sizeof(float)); Zplot = malloc(PHASESTEPS * sizeof(float)); if (Xplot == NULL || Yplot == NULL || Zplot == NULL) #ifdef _WITH_GTK { if (Flags.interactive == ON) { make_dialog (_(errmsg[0])); if (Xplot != NULL) free(Xplot); if (Yplot != NULL) free(Yplot); if (Zplot != NULL) free(Zplot); return; } else nf_error(_(errmsg[5])); } #else nf_error(_(errmsg[5])); #endif sinAlpha = sin(OrbPtr->Inclination - 0.5*PI); cosAlpha = cos(OrbPtr->Inclination - 0.5*PI); /* >>>>>>>>>>> open device <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */ #ifdef _WITH_GNUPLOT if (j == 0) { if (cpgopen("/XSERVE") != 1) #ifdef _WITH_GTK { if (Flags.interactive == ON) { make_dialog (_(errmsg[2])); free(Xplot); free(Yplot); free(Zplot); return; } else nf_error(_(errmsg[2])); } #else nf_error(_(errmsg[2])); #endif } gnu_start(); #else if(cpgopen("/XSERVE") != 1) #ifdef _WITH_GTK { if (Flags.interactive == ON) { make_dialog (_(errmsg[2])); free(Xplot); free(Yplot); free(Zplot); return; } else nf_error(_(errmsg[2])); } #else nf_error(_(errmsg[2])); #endif #endif /* plot box */ if (Orbit.Name[0] != '\0') strncpy(titlestring, Orbit.Name, MAX_CFG_INLINE); else strcpy(titlestring, "Nightfall"); cpgscf(2); cpgslw(1); cpgsch(0.9); cpgsvp(0.05,0.6,0.05,0.9); cpgwnad(-0.98, 0.980, -0.9, 0.9); /* cpgbox("BCNST", 0,0, "BCNST", 0,0); */ cpgsch(2.0); cpgmtxt("T", 1, 0.8, 0.5, titlestring); cpgsch(0.9); /* >>>>>>>>>>>> Light Curve <<<<<<<<<<<<<<<<<<<<<<<<<< */ /* photometric zeropoint */ Calib = FluxOut[0].Mag[Vmag]; /* get photometric data and scale */ js = 0; MinMag = 10.; MaxMag = -10.; FluxPtr = FluxOut; for (i = 0; i <= j; ++i) { if (FluxPtr->Mag[Vmag] != 0) { Xplot[js] = (FluxPtr->Phase/(PI+PI)) - 0.5; Yplot[js] = 2.5 * log10(FluxPtr->Mag[Vmag]/Calib); MinMag = MIN(MinMag,Yplot[js]); MaxMag = MAX(MaxMag,Yplot[js]); ++js; } ++FluxPtr; } /* shift phase if necessary */ test = ON; do { if (Xplot[0] > 0.0) { for (i = 0; i <= js; ++i) { Xplot[i] = Xplot[i] - 1.0; } test = ON; /* there was something to do */ } else { if (Xplot[0] < (-1.0)) { for (i = 0; i <= js; ++i) { Xplot[i] = Xplot[i] + 1.0; } test = ON; /* there was something to do */ } else { test = OFF; } } } while (test == ON); if (fabs(MinMag-MaxMag) < 0.1) MinMag = -0.1; /* get lower window border */ cpgqvp(0, &x1, &x2, &y1, &y2); /* store upper window border */ Upper = y2; /* plot box */ Xtick = (int)(100.*(MaxMag-MinMag+0.1)/3.)/100.; cpgscf(2); cpgslw(1); cpgsch(0.9); cpgsvp(0.6,0.95,y1,y1+(y2-y1)/2.); /* was commented out */ cpgswin(-0.4, 0.9, MinMag-0.05, MaxMag+0.05); cpgbox("BCNST", 0.5, 5, "BCMST", Xtick, 4); cpgmtxt("B", 3, 0.5, 0.5, "Phase"); /* reset point size and plot lightcurve */ cpgsch(0.1); cpgslw(1); cpgsci(2); cpgline(js, Xplot, Yplot); cpgsci(1); /* >>>>>>>>>>>> Radial Velocities <<<<<<<<<<<<<<<<<< */ FluxPtr = FluxOut; MinVel = (FluxPtr->RV[Primary]) /1000.; MaxVel = (FluxPtr->RV[Secondary])/1000.; js = 0; if (Flags.elliptic == ON) { for (i = 0; i <= j; ++i) { if (FluxPtr->Mag[Vmag] != 0) { Yplot[js] = (FluxPtr->RV[Primary]) /1000.; Zplot[js] = (FluxPtr->RV[Secondary])/1000.; MinVel = MIN(MinVel,Yplot[js]); MaxVel = MAX(MaxVel,Yplot[js]); MinVel = MIN(MinVel,Zplot[js]); MaxVel = MAX(MaxVel,Zplot[js]); ++js; } ++FluxPtr; } } else { for (i = 0; i <= j; ++i) { if (FluxPtr->Mag[Vmag] != 0) { Yplot[js] = - (FluxPtr->RV[Primary]) /1000.; Zplot[js] = - (FluxPtr->RV[Secondary])/1000.; MinVel = MIN(MinVel,Yplot[js]); MaxVel = MAX(MaxVel,Yplot[js]); MinVel = MIN(MinVel,Zplot[js]); MaxVel = MAX(MaxVel,Zplot[js]); ++js; } ++FluxPtr; } } if (fabs(MinVel-MaxVel) < 5.) { MinVel = -5.01; MaxVel = 5.01; } /* get lower window border */ cpgqvp(0, &x1, &x2, &y1, &y2); /* plot box */ cpgscf(2); cpgslw(1); cpgsch(0.9); cpgsvp(0.6,0.95,y2,Upper); cpgswin(-0.4, 0.9, MinVel-1., MaxVel+1.); #ifndef _WITH_GNUPLOT cpgbox("BCST", 0.5, 5, "BCMST", (int)((2.+MaxVel-MinVel)/3), 4); #else cpgbox("BCNST", 0.5, 5, "BCMST", (int)((2.+MaxVel-MinVel)/3), 4); #endif /* reset point size and plot lightcurve */ cpgsch(0.1); cpgslw(1); #ifndef _WITH_GNUPLOT /* primary */ cpgsci(5); cpgline(js, Xplot, Yplot); cpgsci(1); /* secondary */ cpgsci(4); cpgline(js, Xplot, Zplot); cpgsci(1); #else cpgline2(js, js, Xplot, Yplot, Xplot, Zplot); #endif #ifndef _WITH_GNUPLOT /* plot the non-keplerian disturbance */ /* zoomed by factor two */ js = 0; FluxPtr = FluxOut; if (Flags.elliptic == OFF){ for (i = 0; i <= j; ++i) { if (FluxPtr->Mag[Vmag] != 0) { Yplot[js] = -(FluxPtr->D_RV[Primary]) /500.; Zplot[js] = -(FluxPtr->D_RV[Secondary])/500.; ++js; } ++FluxPtr; } } else { for (i = 0; i <= j; ++i) { if (FluxPtr->Mag[Vmag] != 0) { Yplot[js] = (FluxPtr->D_RV[Primary]) /500.; Zplot[js] = (FluxPtr->D_RV[Secondary])/500.; ++js; } ++FluxPtr; } } cpgsls(4); /* primary */ cpgsci(5); cpgline(js, Xplot, Yplot); cpgsci(1); /* secondary */ cpgsci(4); cpgline(js, Xplot, Zplot); cpgsci(1); cpgsls(1); #endif /* >>>>>>>>>>>>> the stars <<<<<<<<<<<<<<<<<<<<<<< */ /* plot box */ cpgscf(2); cpgslw(1); cpgsch(0.9); cpgsvp(0.05,0.6,0.05,0.9); cpgwnad(-0.98, 0.980, -0.9, 0.9); cpgsch(2.0); cpgmtxt("T", 1, 0.8, 0.5, titlestring); cpgsch(0.9); /* reset size */ cpgsch(0.1); cpgslw(1); #ifdef HAVE_DISK if (j == 0) { /* die scheibe muss raumfest sein, z.B. wg. tilt */ init_phase = OrbPtr->Phase; } #endif /* loop over components */ #if defined(PRINT_COORDS) if (j == 0) { outfile = fopen ("NightfallCoordinates.txt", "w"); if (!outfile) perror("Error opening NightfallCoordinates.txt"); } #endif CountElem = 0; #ifdef HAVE_DISK if (Flags.disk == ON) MaxComp = 3; #endif for (Comp = 0; Comp < MaxComp; ++Comp) { /* now we need the phase */ /* we get it from Orbit.Phase, which is set either */ /* in KeplerSolve or in LightCurve */ /* so this routine must be called after LightCurve */ if (Comp == Primary) { sinPhi = sin(-OrbPtr->Phase); cosPhi = cos(-OrbPtr->Phase); C = 1.0; } else if (Comp == Secondary) { C = -1.0; /* mirror y -> (-y) for secondary */ sinPhi = sin(-OrbPtr->Phase + PI); cosPhi = cos(-OrbPtr->Phase + PI); } else { /* (Comp == Disk) */ C = -1.0; /* mirror y -> (-y) for secondary */ sinPhi = sin(-OrbPtr->Phase + PI); cosPhi = cos(-OrbPtr->Phase + PI); } /* the centre of mass */ if (Comp == Primary || Comp == Secondary) com = (((float)(Binary[Comp].Mq))/( 1.0 + (float)(Binary[Comp].Mq))); else com = (((float)(Binary[Secondary].Mq))/( 1.0 + (float)(Binary[Secondary].Mq))); /* the transformation matrix */ /* first column, then row */ /* Trans[0][0] = cosPhi * cosAlpha; */ Trans[1][0] = sinPhi; Trans[2][0] = cosPhi * sinAlpha; /* Trans[3][0] = 0.0; */ /* Trans[0][1] = -sinPhi * cosAlpha;*/ Trans[1][1] = cosPhi; Trans[2][1] = -sinPhi * sinAlpha; /* Trans[3][1] = 0.0; */ /* Trans[0][2] = -sinAlpha; */ Trans[1][2] = 0.0; Trans[2][2] = cosAlpha; /* Trans[3][2] = 0.0; */ /* Trans[0][3] = com * (1.0 - cosPhi * cosAlpha); */ Trans[1][3] = -com * sinPhi; Trans[2][3] = -com * cosPhi * sinAlpha; /* Trans[3][3] = 1.0; */ /* Trans[0][3] = Binary[Comp].RLag1 * (1.0 - cosPhi * cosAlpha); Trans[1][3] = -Binary[Comp].RLag1 * sinPhi; Trans[2][3] = -Binary[Comp].RLag1 * cosPhi * sinAlpha; Trans[3][3] = 1.0; */ t10 = Trans[1][0]; t11 = Trans[1][1]; t12 = Trans[1][2]; t13 = Trans[1][3]; t20 = Trans[2][0]; t21 = Trans[2][1]; t22 = Trans[2][2]; t23 = Trans[2][3]; /* do the transformation */ SurfPtr = Surface[Comp]; for (i = 0; i < Binary[Comp].NumElem; ++i) { if (SurfPtr->EclipseFlag == 0 && SurfPtr->SpotNum == 0) { /* Observer looks from positive X-Axis */ /* this is the negative Y-Axis */ Xplot[CountElem] = - ( SurfPtr->lambda * t10 + SurfPtr->mu * t11 * C + SurfPtr->nu * t12 + t13); /* + 1.0 * t13); */ /* and this the Z-Axis */ Yplot[CountElem] = ( SurfPtr->lambda * t20 + SurfPtr->mu * t21 * C + SurfPtr->nu * t22 + t23); /* + 1.0 * t23); */ /* increase counter for elements in Xplot/Yplot */ ++CountElem; } ++SurfPtr; } } #if defined(PRINT_COORDS) hyp = sqrt(((sx[0] - sx[1])*(sx[0] - sx[1]))+((sy[0] - sy[1])*(sy[0] - sy[1]))); ssin = (sy[0] - sy[1])/hyp; scos = (sx[0] - sx[1])/hyp; posw = atan((sy[0] - sy[1])/(sx[0] - sx[1]))*(180/M_PI); if (ssin < 0 && scos > 0) posw += 180; if (ssin > 0 && scos > 0) posw += 180; if (ssin > 0 && scos < 0) posw += 360; if (outfile) { fprintf (outfile, "%3d %12.6f %12.6f %12.6f\n", j, Orbit.Dist*hyp, posw /* , ssin, scos */, (FluxOut[j].Phase/(PI+PI)) - 0.5); fflush(outfile); } #endif /* plot the data */ #ifdef _WITH_GNUPLOT cpgsvp(0.05,0.6,0.05,0.9); cpgwnad(-0.98, 0.980, -0.9, 0.9); cpgbox("BCST", 0, 0, "BCST", 0, 0); cpgmtxt("T", 1, 0.8, 0.5, titlestring); #endif cpgpt(CountElem, Xplot, Yplot, 17); /* plot finished */ #ifndef _WITH_GNUPLOT cpgend(); #else gnu_end(); if (j == (PhaseSteps - 1)) cpgend(); #endif free(Xplot); free(Yplot); free(Zplot); return; } #endif