/* 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>
#include "Light.h"
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Search index for temperature in array
@param (float) key temperature
@param (float *) array temperature array
@param (int) n size of array
@return (int) the left index of the two bracketing values
@heading Light Curve
****************************************************************************/
static int search_lkey (float key, float * r, int n )
{
static int i = 0;
static int high, low;
/* speedup (?)
*/
if (key >= r[i] && key < r[i+1])
{
return ( i );
}
/* return 0 if < r[0];
*/
for ( low = (-1), high = n; high-low > 1; )
{
i = (high+low) / 2;
if ( key < r[i])
high = i;
else
low = i;
}
/* make binary search stable by explicitely selecting the low value
*/
i = (low < 0) ? 0 : low;
return ( i );
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichmann@hs.uni-hamburg.de)
@version 1.0
@short Compute limb darkening coefficients
@long Limb darkening coefficients
from Claret A. 2000 Astron. Astrophys. 363, 1081
@tip log(g) = 3.5, 4.0, 4.5, 5.0
@param (int) Comp The stellar component
@return (void)
@heading Light Curve
****************************************************************************/
static void FindbestLC ();
#include "LDC_PHOENIX_3.5.h"
#include "LDC_PHOENIX_4.0.h"
#include "LDC_PHOENIX_4.5.h"
#include "LDC_PHOENIX_5.0.h"
/* For surface gravities not extending to 50000 K, the last value
is repeated to fill the tables. */
#include "LDC_ATLAS_3.5.h"
#include "LDC_ATLAS_4.0.h"
#include "LDC_ATLAS_4.5.h"
#include "LDC_ATLAS_5.0.h"
#define hMAX 35
#define kMAX 76
/* stellar temperature grid for limb darkening coefficients */
static
float TH[hMAX] ={
2000., 2200., 2400., 2600., 2800., 3000., 3200., 3400., 3600., 3800.,
4000., 4200., 4400., 4600., 4800., 5000., 5200., 5400., 5600., 5800.,
7000., 7200., 7400., 7600., 7800., 8000., 8200., 8400., 8600., 8800.,
9000., 9200., 9400., 9600., 9800.};
static
float TK[kMAX] ={
3500., 3750., 4000., 4250., 4500., 4750., 5000., 5250., 5500.,
5750., 6000., 6250., 6500., 6750., 7000., 7250., 7500., 7750.,
8000., 8250., 8500., 8750., 9000., 9250., 9500., 9750., 10000.,
10250., 10500., 10750., 11000., 11250., 11500., 11750., 12000., 12250.,
12500., 12750., 13000., 14000., 15000., 16000., 17000., 18000., 19000.,
20000., 21000., 22000., 23000., 24000., 25000., 26000., 27000., 28000.,
29000., 30000., 31000., 32000., 33000., 34000., 35000., 36000., 37000.,
38000., 39000., 40000., 41000., 42000., 43000., 44000., 45000., 46000.,
47000., 48000., 49000., 50000
};
void LightLimbDark(int Comp)
{
float Tn; /* temperature */
register long j; /* loop variable */
int band; /* passband */
double log_g = Binary[Comp].log_g; /* log g */
SurfaceElement *SurfPtr; /* pointer to surface */
int tnum; /* size of LDC arrays */
float *T, *u, *a, *b, *c, *d; /* pointers to T, LDC arrays */
/* ----------- initialize ------------------------------------------ */
/* array[p][q] = ARRVAL(array, max(q), p, q)
*/
#define ARRVAL(array, ncolumns, i, j) array[(i) * (ncolumns) + (j)]
SurfPtr = Surface[Comp];
for (band = 0; band < NUM_MAG; ++band) {
for (j = 0; j < 2; ++j) {
Limb[Comp][band][j] = 0.0;
}
}
/* we approximate T = Binary[Comp].Temperature for all */
/* surface elements */
Tn = (float) Binary[Comp].Temperature;
if (Tn > SWITCH_TEMP) {
tnum = kMAX;
T = TK;
if (log_g <= 3.5) {
u = &kL_35[0][0]; a = &kQ1_35[0][0]; b = &kQ2_35[0][0];
c = &kS1_35[0][0]; d = &kS2_35[0][0];
}
else if (fabs(log_g - 4.0) < FLT_EPSILON) {
u = &kL_40[0][0]; a = &kQ1_40[0][0]; b = &kQ2_40[0][0];
c = &kS1_40[0][0]; d = &kS2_40[0][0];
}
else if (fabs(log_g - 4.5) < FLT_EPSILON) {
u = &kL_45[0][0]; a = &kQ1_45[0][0]; b = &kQ2_45[0][0];
c = &kS1_45[0][0]; d = &kS2_45[0][0];
}
else {
u = &kL_50[0][0]; a = &kQ1_50[0][0]; b = &kQ2_50[0][0];
c = &kS1_50[0][0]; d = &kS2_50[0][0];
}
}
else {
tnum = hMAX;
T = TH;
if (log_g <= 3.5) {
u = &hL_35[0][0]; a = &hQ1_35[0][0]; b = &hQ2_35[0][0];
c = &hS1_35[0][0]; d = &hS2_35[0][0];
}
else if (fabs(log_g - 4.0) < FLT_EPSILON) {
u = &hL_40[0][0]; a = &hQ1_40[0][0]; b = &hQ2_40[0][0];
c = &hS1_40[0][0]; d = &hS2_40[0][0];
}
else if (fabs(log_g - 4.5) < FLT_EPSILON) {
u = &hL_45[0][0]; a = &hQ1_45[0][0]; b = &hQ2_45[0][0];
c = &hS1_45[0][0]; d = &hS2_45[0][0];
}
else {
u = &hL_50[0][0]; a = &hQ1_50[0][0]; b = &hQ2_50[0][0];
c = &hS1_50[0][0]; d = &hS2_50[0][0];
}
}
/* ----------- linear law ----------------------------------------- */
if (Flags.limb == 0) {
j = search_lkey (Tn, T, tnum);
for (band = 0; band < NUM_MAG; ++band) {
/* Limb[Comp][band][0] = u[band][j]; */
Limb[Comp][band][0] = ARRVAL(u, tnum, band, j);
}
}
/* ----------- quadratic law--------------------------------------- */
else if (Flags.limb == 1) {
j = search_lkey (Tn, T, tnum);
for (band = 0; band < NUM_MAG; ++band) {
/* Limb[Comp][band][0] = a[band][j]; */
/* Limb[Comp][band][1] = b[band][j]; */
Limb[Comp][band][0] = ARRVAL(a, tnum, band, j);
Limb[Comp][band][1] = ARRVAL(b, tnum, band, j);
}
}
/* ----------- square root law------------------------------------- */
else if (Flags.limb == 2) {
j = search_lkey (Tn, T, tnum);
for (band = 0; band < NUM_MAG; ++band) {
/* Limb[Comp][band][0] = c[band][j]; */
/* Limb[Comp][band][1] = d[band][j]; */
Limb[Comp][band][0] = ARRVAL(c, tnum, band, j);
Limb[Comp][band][1] = ARRVAL(d, tnum, band, j);
}
}
/* ----------- inverse ------------------------------------- */
else if (Flags.limb == 3) {
j = search_lkey (Tn, T, tnum);
for (band = 0; band < NUM_MAG; ++band) {
/* Limb[Comp][band][0] = u[band][j]; */
Limb[Comp][band][0] = ARRVAL(u, tnum, band, j);
}
}
/* ----------- Modify, if using monochrome wavelengths ------------ */
if (Flags.monochrome == ON)
FindbestLC ();
/* ----------- Verbose output ------------------------------------- */
if (Flags.debug[VERBOSE] == ON) {
printf(_("Limb Darkening Coefficients: \n"));
if(Comp == Primary) {
printf(_(" Primary: \n"));
printf(" U: %6.3f %6.3f B: %6.3f %6.3f V: %6.3f %6.3f\n",
Limb[Primary][Umag][0],Limb[Primary][Umag][1],
Limb[Primary][Bmag][0],Limb[Primary][Bmag][1],
Limb[Primary][Vmag][0],Limb[Primary][Vmag][1]);
printf(" u: %6.3f %6.3f v: %6.3f %6.3f \n",
Limb[Primary][umag][0],Limb[Primary][umag][1],
Limb[Primary][vmag][0],Limb[Primary][vmag][1]);
printf(" b: %6.3f %6.3f y: %6.3f %6.3f\n",
Limb[Primary][bmag][0],Limb[Primary][bmag][1],
Limb[Primary][ymag][0],Limb[Primary][ymag][1]);
printf(" R: %6.3f %6.3f I: %6.3f %6.3f\n",
Limb[Primary][Rmag][0],Limb[Primary][Rmag][1],
Limb[Primary][Imag][0],Limb[Primary][Imag][1]);
printf(" J: %6.3f %6.3f H: %6.3f %6.3f K: %6.3f %6.3f\n",
Limb[Primary][Jmag][0],Limb[Primary][Jmag][1],
Limb[Primary][Hmag][0],Limb[Primary][Hmag][1],
Limb[Primary][Kmag][0],Limb[Primary][Kmag][1]);
}
if(Comp == Secondary) {
printf(_(" Secondary: \n"));
printf(" U: %6.3f %6.3f B: %6.3f %6.3f V: %6.3f %6.3f\n",
Limb[Secondary][Umag][0],Limb[Secondary][Umag][1],
Limb[Secondary][Bmag][0],Limb[Secondary][Bmag][1],
Limb[Secondary][Vmag][0],Limb[Secondary][Vmag][1]);
printf(" u: %6.3f %6.3f v: %6.3f %6.3f \n",
Limb[Secondary][umag][0],Limb[Secondary][umag][1],
Limb[Secondary][vmag][0],Limb[Secondary][vmag][1]);
printf(" b: %6.3f %6.3f y: %6.3f %6.3f\n",
Limb[Secondary][bmag][0],Limb[Secondary][bmag][1],
Limb[Secondary][ymag][0],Limb[Secondary][ymag][1]);
printf(" R: %6.3f %6.3f I: %6.3f %6.3f\n",
Limb[Secondary][Rmag][0],Limb[Secondary][Rmag][1],
Limb[Secondary][Imag][0],Limb[Secondary][Imag][1]);
printf(" J: %6.3f %6.3f H: %6.3f %6.3f K: %6.3f %6.3f\n",
Limb[Secondary][Jmag][0],Limb[Secondary][Jmag][1],
Limb[Secondary][Hmag][0],Limb[Secondary][Hmag][1],
Limb[Secondary][Kmag][0],Limb[Secondary][Kmag][1]);
}
#ifdef HAVE_DISK
if(Comp == Disk) {
printf(_(" Disk: \n"));
printf(" U: %6.3f %6.3f B: %6.3f %6.3f V: %6.3f %6.3f\n",
Limb[Disk][Umag][0],Limb[Disk][Umag][1],
Limb[Disk][Bmag][0],Limb[Disk][Bmag][1],
Limb[Disk][Vmag][0],Limb[Disk][Vmag][1]);
printf(" u: %6.3f %6.3f v: %6.3f %6.3f \n",
Limb[Disk][umag][0],Limb[Disk][umag][1],
Limb[Disk][vmag][0],Limb[Disk][vmag][1]);
printf(" b: %6.3f %6.3f y: %6.3f %6.3f\n",
Limb[Disk][bmag][0],Limb[Disk][bmag][1],
Limb[Disk][ymag][0],Limb[Disk][ymag][1]);
printf(" R: %6.3f %6.3f I: %6.3f %6.3f\n",
Limb[Disk][Rmag][0],Limb[Disk][Rmag][1],
Limb[Disk][Imag][0],Limb[Disk][Imag][1]);
printf(" J: %6.3f %6.3f H: %6.3f %6.3f K: %6.3f %6.3f\n",
Limb[Disk][Jmag][0],Limb[Disk][Jmag][1],
Limb[Disk][Hmag][0],Limb[Disk][Hmag][1],
Limb[Disk][Kmag][0],Limb[Disk][Kmag][1]);
}
#endif
}
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Verbose output of wavelengths
@param (int) Comp The stellar component
@return (void)
@heading Light Curve
****************************************************************************/
static
void VerboseMessage (int Comp)
{
if (Flags.debug[VERBOSE] == ON) {
printf("\n");
printf(_("Effective Wavelengths:\n"));
if(Comp == Primary) {
printf(_(" Primary: \n"));
printf(" U: %6.4f B: %6.4f V: %6.4f R: %6.4f I: %6.4f\n",
WaveL[Primary][Umag], WaveL[Primary][Bmag],
WaveL[Primary][Vmag], WaveL[Primary][Rmag],
WaveL[Primary][Imag]);
printf(" u: %6.4f v: %6.4f b: %6.4f y: %6.4f\n",
WaveL[Primary][umag], WaveL[Primary][vmag],
WaveL[Primary][bmag], WaveL[Primary][ymag]);
printf(" J: %6.4f H: %6.4f K: %6.4f\n",
WaveL[Primary][Jmag], WaveL[Primary][Hmag],
WaveL[Primary][Kmag]);
}
if(Comp == Secondary) {
printf(_(" Secondary: \n"));
printf(" U: %6.4f B: %6.4f V: %6.4f R: %6.4f I: %6.4f\n",
WaveL[Secondary][Umag], WaveL[Secondary][Bmag],
WaveL[Secondary][Vmag], WaveL[Secondary][Rmag],
WaveL[Secondary][Imag]);
printf(" u: %6.4f v: %6.4f b: %6.4f y: %6.4f\n",
WaveL[Secondary][umag], WaveL[Secondary][vmag],
WaveL[Secondary][bmag], WaveL[Secondary][ymag]);
printf(" J: %6.4f H: %6.4f K: %6.4f\n",
WaveL[Secondary][Jmag], WaveL[Secondary][Hmag],
WaveL[Secondary][Kmag]);
}
#ifdef HAVE_DISK
if(Comp == Disk) {
printf(_(" Disk: \n"));
printf(" U: %6.4f B: %6.4f V: %6.4f R: %6.4f I: %6.4f\n",
WaveL[Disk][Umag], WaveL[Disk][Bmag],
WaveL[Disk][Vmag], WaveL[Disk][Rmag],
WaveL[Disk][Imag]);
printf(" u: %6.4f v: %6.4f b: %6.4f y: %6.4f\n",
WaveL[Disk][umag], WaveL[Disk][vmag],
WaveL[Disk][bmag], WaveL[Disk][ymag]);
printf(" J: %6.4f H: %6.4f K: %6.4f\n",
WaveL[Disk][Jmag], WaveL[Disk][Hmag],
WaveL[Disk][Kmag]);
}
#endif
}
return;
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Set monochromatic wavelengths for blackbody
@param (int) Comp The stellar component
@return (void)
@heading Light Curve
****************************************************************************/
float mono_zero[NUM_MAG] = {
0.3585, 0.4383, 0.5502,
0.6400, 0.7900,
1.25, 1.62, 2.2,
0.3500, 0.4110, 0.4670, 0.5470
};
static
void SetMonoWavelengths(int Comp)
{
int band;
for (band = 0; band < NUM_MAG; ++band)
{
WaveL[Comp][band] = mono_zero[band];
}
/* Verbose output */
VerboseMessage(Comp);
return;
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Define monochromatic wavelengths for blackbody
@param (void)
@return (int) Flag ON/OFF
@heading Light Curve
****************************************************************************/
int DefineMonoWavelengths ()
{
int band;
int len;
float fin;
char * pp;
char * qq;
if (NULL == (pp = getenv("NIGHTFALL_MONO_WAVE")))
return OFF;
len = strlen(pp) + 1;
qq = (char *) malloc (len);
if (NULL == qq)
{
nf_error(_(errmsg[0]));
return OFF;
}
strncpy (qq, pp, len);
qq[len-1] = '\0';
band = 0;
pp = strtok(qq, ",");
if (pp != NULL)
{
fin = atof(pp);
mono_zero[band] = (fin > 0) ? fin : mono_zero[band];
++band;
}
else
return OFF;
while (pp != NULL && band < NUM_MAG)
{
pp = strtok(NULL, ",");
if (pp != NULL)
{
fin = atof(pp);
mono_zero[band] = (fin > 0) ? fin : mono_zero[band];
++band;
}
}
return ON;
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Search closest value in an array
@param (float) fin A float value
@param (float *) farr A float array
@param (int) arrsiz Size of farr
@return (int) The index of closest match
@heading Light Curve
****************************************************************************/
static
int BestIndex (float fin, float * farr, int arrsiz)
{
register int k;
int i = 0;
float tmp;
float best = fabs (fin - farr[0]);
for (k = 1; k < arrsiz; ++k)
{
tmp = fabs (fin - farr[k]);
if (tmp < best)
{
best = tmp;
i = k;
}
}
return i;
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Find best limb darkening coefficient for monochrome wavelengths
@param (void)
@return (void)
@heading Light Curve
****************************************************************************/
static double tmp_limb[2][NUM_MAG][2];
/* effective wavelength */
/*
U, B, V,
R, I,
J, H, K,
u, v, b, y
*/
static float l_zero[NUM_MAG] = {
0.3585, 0.4383, 0.5502,
0.6400, 0.7900,
1.25, 1.62, 2.2,
0.3500, 0.4110, 0.4670, 0.5470
};
static
void FindbestLC ()
{
register int i, j, k;
for (k = 0; k < NUM_MAG; ++k)
{
i = BestIndex (mono_zero[k], l_zero, NUM_MAG);
tmp_limb[0][k][0] = Limb[0][i][0];
tmp_limb[1][k][0] = Limb[1][i][0];
tmp_limb[0][k][1] = Limb[0][i][1];
tmp_limb[1][k][1] = Limb[1][i][1];
}
/* copy from temporary array
*/
for (i = 0; i < 2; ++i)
for (j = 0; j < 2; ++j)
for (k = 0; k < NUM_MAG; ++k)
Limb[i][k][j] = tmp_limb[i][k][j];
return;
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Compute temperature-dependent effective wavelengths for blackbody
@long Uses second moments of passbands following the prescription by:
Young A.T. 1992, Astronomy and Astrophysics 257, 366 (Equation 12)
Effective wavelengths
calculated following Equation 3.30 in:
Budding E. 1993, An Introduction to Astronomic Photometry,
Cambridge University Press
@param (int) Comp The stellar component
@return (void)
@heading Light Curve
****************************************************************************/
void EffectiveWavelength(int Comp)
{
double l_max; /* blackbody peak */
int band; /* passband */
double Wien = 2897.756; /* constant for Wien's law */
/* second moments of passbands */
/* approximated following the prescription by */
/* Young A.T. 1992, Astronomy and Astrophysics 257, 366; Equation 12 */
static float M2ND[NUM_MAG] = {
0.000448, 0.001236, 0.001405,
4.556e-3, 1.845e-3,
4.648e-3, 5.971e-3, 8.264e-3,
1.8061e-5, 1.3514e-5, 1.6357e-5, 2.8619e-5
};
/* --------------------- monochromatic wavelength -------------------- */
if (Flags.monochrome == ON)
{
SetMonoWavelengths(Comp);
return;
}
/* --------------------- effective wavelength ------------------------ */
/* effective wavelengths */
/* calculated following Equation 3.30 in */
/* Budding E. 1993, An Introduction to Astronomic Photometry */
/* Cambridge University Press */
/* blackbody peak */
l_max = Wien/Binary[Comp].Temperature;
for (band = 0; band < NUM_MAG; ++band) {
WaveL[Comp][band] = l_zero[band]
+ 5*(l_max-l_zero[band])*M2ND[band]/SQR(l_zero[band]);
}
/* Verbose output */
VerboseMessage(Comp);
return;
}
syntax highlighted by Code2HTML, v. 0.9.1