/* 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 <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <float.h>
#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;
}
syntax highlighted by Code2HTML, v. 0.9.1