/* 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 #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; }