/* NIGHTFALL Light Curve Synthesis Program */ /* Copyright (C) 2000 Rainer Wichmann */ /* */ /* This program is free software; you can redistribute it */ /* and/or modify */ /* it under the terms of the GNU General Public License as */ /* published by */ /* the Free Software Foundation; either version 2 of the License, or */ /* (at your option) any later version. */ /* */ /* This program is distributed in the hope that it will be useful, */ /* but WITHOUT ANY WARRANTY; without even the implied warranty of */ /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */ /* GNU General Public License for more details. */ /* */ /* You should have received a copy of the GNU General Public License */ /* along with this program; if not, write to the Free Software */ /* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */ #include #include #include #include #include #include "Light.h" char * MapPath = NULL; int MapBand = 2; /**************************************************************************** @package nightfall @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de) @version 1.0 @short Set master path to map @param (void) @return (void) @heading Surface Map ****************************************************************************/ int SetMapPath () { int len; char * p = getenv("NIGHTFALL_SMAP_PATH"); if (p == NULL) { return OFF; } len = strlen(p) + 1; MapPath = (char *) malloc (len); if (MapPath == NULL) { nf_error(_(errmsg[0])); return OFF; } strncpy (MapPath, p, len); p[len-1] = '\0'; return ON; } /**************************************************************************** @package nightfall @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de) @version 1.0 @short Set passband for map @param (void) @return (void) @heading Surface Map ****************************************************************************/ void SetMapBand () { char * p = getenv("NIGHTFALL_SMAP_BAND"); if (p == NULL) return; MapBand = atoi(p); MapBand = (MapBand >= NUM_MAG) ? (NUM_MAG-1) : MapBand; MapBand = (MapBand < 0) ? 0 : MapBand; return; } /**************************************************************************** @package nightfall @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de) @version 1.0 @short Write the map @param (int) thisPhaseStep The phase index @return (void) @heading Surface Map ****************************************************************************/ void PrintSurfaceMap (int thisPhaseStep) { double Trans[4][4]; /* transformation matrix */ double t10, t11, t12, t13; /* transformation matrix */ double t20, t21, t22, t23; /* transformation matrix */ double t13_C, t23_C; SurfaceElement *SurfPtr; /* pointer to surface */ long i; /* loop variable */ int Comp; /* the stellar component */ long CountElem; /* # of elements in X/Yplot */ double sinAlpha; /* sin inclination */ double cosAlpha; /* cos inclination */ double sinPhi; /* sin phase */ double cosPhi; /* cos phase */ double C; /* sign (for primary/secondary) */ float com; /* centre */ float com2; /* centre */ float Flux; /* flux */ float Xplot, Yplot; /* x, y position */ float Xplot_C, Yplot_C; /* x, y position */ OrbitType *OrbPtr = &Orbit; /* pointer to Orbit */ char * thisMapPath = NULL; /* File to hold surface map */ FILE * of; /* File handle for above */ double Cosgamma; /* LOS angle */ double Common; /* common part */ double OneMinCosgamma; /* 1.0 - Cosgamma */ if (NULL == MapPath) return; thisMapPath = (char *) malloc (strlen(MapPath) + 64); if (thisMapPath == NULL) { nf_error(_(errmsg[0])); return; } sprintf(thisMapPath, "%s.%02d.%03d", MapPath, MapBand, thisPhaseStep); of = fopen(thisMapPath, "w"); if (of == NULL) { nf_error(_(errmsg[0])); return; } if (Flags.debug[VERBOSE] == ON || Flags.debug[WARN] == ON) fprintf(stderr, _("Writing surface map %s.\n"), thisMapPath); /* loop over components */ CountElem = 0; sinAlpha = sin(OrbPtr->Inclination - 0.5*PI); cosAlpha = cos(OrbPtr->Inclination - 0.5*PI); for (Comp = 0; Comp < 2; ++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; com = 0.0; } else { C = -1.0; /* mirror y -> (-y) for secondary */ sinPhi = sin(-OrbPtr->Phase + PI); cosPhi = cos(-OrbPtr->Phase + PI); } /* the centre of the star */ com = 0.0; com2 = (((float)(Binary[Comp].Mq))/( 1.0 + (float)(Binary[Comp].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; */ t13_C = -com2 * sinPhi; t23_C = -com2 * cosPhi * sinAlpha; /* 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->CosGamma >= 0.0) { /* Observer looks from positive X-Axis */ /* this is the negative Y-Axis */ Xplot = - (SurfPtr->lambda * t10 + SurfPtr->mu * t11 * C + SurfPtr->nu * t12 + 1.0 * t13); /* and this the Z-Axis */ Yplot = (SurfPtr->lambda * t20 + SurfPtr->mu * t21 * C + SurfPtr->nu * t22 + 1.0 * t23); Xplot_C = Xplot - t13 + t13_C; Yplot_C = Yplot - t23 + t23_C; Cosgamma = SurfPtr->CosGamma; Common = SurfPtr->visibility * Cosgamma; OneMinCosgamma = (1.0 - Cosgamma); if (Flags.limb == 0) { Flux = SurfPtr->f_[MapBand] * Common * (1.0 - Limb[Comp][MapBand][0]*OneMinCosgamma); } else if ( Flags.limb == 1) { Flux = SurfPtr->f_[MapBand] * Common * (1.0 - Limb[Comp][MapBand][0]*OneMinCosgamma - Limb[Comp][MapBand][1]*(1.0- Cosgamma*Cosgamma)); } else if ( Flags.limb == 2) { Flux = SurfPtr->f_[MapBand] * Common * (1.0 - Limb[Comp][MapBand][0]*OneMinCosgamma - Limb[Comp][MapBand][1]*(1.0- sqrt(Cosgamma))); } else if (Flags.limb == 3 && Comp < Disk) { Flux = SurfPtr->f_[MapBand] * Common * SurfPtr->AreaType * (1.0 - SurfPtr->AreaType * Limb[Comp][MapBand][0]*OneMinCosgamma); } else if (Flags.limb == 3 && Comp >= Disk) { Flux = SurfPtr->f_[MapBand] * Common * (1.0 - Limb[Comp][MapBand][0]*OneMinCosgamma); } else Flux = 0.0; /* increase counter for elements in Xplot/Yplot */ fprintf(of, "%1d %5ld %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f %12.10f %g\n", Comp+1, i, Xplot, Yplot, Xplot_C, Yplot_C, Cosgamma, SurfPtr->temp, SurfPtr->grav, SurfPtr->area, Flux); ++CountElem; } ++SurfPtr; } fprintf(of, "\n"); } fclose(of); return; }