/* 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 <string.h>
#include <float.h>
#include "Light.h"
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Blackbody flux for an element
@param (SurfaceElement *) SurfPtr the surface element
@return (void)
@heading Flux Computation
****************************************************************************/
static void BBFlux(SurfaceElement *SurfPtr, int Comp)
{
OrbitType *OrbPtr = &Orbit; /* orbit pointer */
int band; /* the passband */
double Temp; /* temperature */
double ScaledArea; /* scaled surface area */
static double CC1 = 3.7418498e+15;/* 10e(5*4) erg cm / sec */
static double CC2 = 14388.3361; /* 10e4 cm K */
/* WaveL in 10e4 cm */
/* we compute I(1) -- Intensity normal to surface */
/* in blackbody approximation */
/* B = CC1/WaveL^5 * 1/(exp(CC2/(T*WaveL))-1) WaveL(ength) in _cm_ */
/* ----------- surface rescaling --------------------------------- */
ScaledArea = SurfPtr->area * SQR(OrbPtr->Dist);
/* ----------- flux computation --------------------------------- */
Temp = SurfPtr->temp;
for (band = 0; band < NUM_MAG; ++band)
{
SurfPtr->f_[band] = ScaledArea
* CC1/
(WaveL[Comp][band] * WaveL[Comp][band] * WaveL[Comp][band]
* WaveL[Comp][band] * WaveL[Comp][band])
* 1.0/(exp(CC2/(WaveL[Comp][band] *
Temp))-1.0);
/*
if (myrank == (numprocs - 1) && band == 2)
{
fprintf(stderr, "%3d: ", band);
fprintf(stderr, "%8.4g %8.4g %8.4g %8.4g %8.4g %8.4g %8.4g %8.4g",
SurfPtr->f_[band], ScaledArea, SurfPtr->area, OrbPtr->Dist,
WaveL[Comp][band], CC1, CC2, Temp);
fprintf(stderr, "\n");
}
*/
}
return;
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Computation of blackbody fluxes
@param (int) Comp The stellar component
@return (void)
@heading Flux Computation
****************************************************************************/
void BlackbodyFlux(int Comp)
{
SurfaceElement *SurfPtr; /* pointer to surface */
long i; /* loop variable */
SurfPtr = Surface[Comp];
for (i = 0; i < Binary[Comp].NumElem; ++i)
{
BBFlux(SurfPtr, Comp);
++SurfPtr;
}
}
/****************************************************************************
@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 a centered 4-value subarray
@heading Flux Computation
****************************************************************************/
static int search_key (float key, float * r, int n )
{
static int i = 0;
int high, low, k;
/* speedup (?)
*/
if (key >= r[i] && key < r[i+1])
{
k = MAX ( (i-1), 0);
k = MIN ( k, (n-4));
return ( k );
}
/* 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;
k = MAX ( (i-1), 0);
k = MIN ( k, (n-4));
return ( k );
}
/* Kurucz ATLAS Intensities from
*
* http://kurucz.harvard.edu/grids/gridP00/ip00k2.pck19
* (file dated 10-Mar-98 18:50)
*/
#include "FluxKurucz_35.h"
#include "FluxKurucz_45.h"
#include "FluxKurucz_40.h"
#include "FluxKurucz_50.h"
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.1
@short Return upper limits for atmosphere models
@param (float) log_g The surface gravity
@return (float)
@heading Flux Computation
****************************************************************************/
float ModelLimits(float log_g)
{
/*********************************************
if (fabs(log_g - 3.5) < FLT_EPSILON)
return k_35_high;
else if (fabs(log_g - 4.0) < FLT_EPSILON)
return k_40_high;
else if (fabs(log_g - 4.5) < FLT_EPSILON)
return k_45_high;
else if (fabs(log_g - 5.0) < FLT_EPSILON)
return k_50_high;
else
return k_35_high;
***********************************************/
/* Changed to return BB limit, will cause automatic fallback on BB */
return LIM_TB_H;
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Computation of atmospheric model fluxes
@param (int) Comp The stellar component
@return (void)
@heading Flux Computation
****************************************************************************/
void ModelFlux(int Comp)
{
/* PHOENIX Intensities from Peter Hauschildt
* Models below 3000K are 'dusty'
*/
#include "FluxHauschildt_35.h"
#include "FluxHauschildt_40.h"
#include "FluxHauschildt_45.h"
#include "FluxHauschildt_50.h"
/* Scaling factors between models and/or BB fluxes
*/
#include "FluxScales_50.h"
#include "FluxScales_45.h"
#include "FluxScales_40.h"
#include "FluxScales_35.h"
int maxnum; /* number of models */
int maxnum_h; /* number of Hauschildt models */
int maxnum_k; /* number of Kurucz models */
float * tarr; /* the model temperatures */
float * tarr_k; /* the model temperatures (Ku) */
float * tarr_h; /* the model temperatures (Ha) */
float * farr[NUM_MAG]; /* the model fluxes */
float * farr_k[NUM_MAG]; /* the model fluxes (Ku) */
float * farr_h[NUM_MAG]; /* the model fluxes (Ha) */
float * adjust; /* scale factors to relate models */
float * adjust_k; /* model scale factors (Ku) */
float * adjust_bb_lo; /* model scale factors (BB lower) */
float * adjust_bb_hi; /* model scale factors (BB upper) */
int last; /* last switch */
SurfaceElement * SurfPtr; /* pointer to surface */
OrbitType *OrbPtr = &Orbit; /* orbit pointer */
register int j; /* loop variable */
register long i; /* loop variable */
int band; /* passband */
int hkey; /* hash key */
double T; /* temperature */
double k_high; /* max. temperature Kurucz */
double h_low; /* min. temperature Hauschildt */
double Hi, Mid, Lo; /* bracketing temperatures in model */
double FHi, FMid, FLo; /* corresponding fluxes */
double FLLo, LLo; /* more bracketing temperatures */
double Phi; /* temperature */
double Cubic; /* cubic interpolation */
double ScaledArea; /* scaled surface area */
double ScaleFactor; /* scale factor for surface area */
double Parabolic; /* parabolic interpolation */
/* unit scaling factors for PHOENIX
*/
static float adjust_unity[NUM_MAG] = {
1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
1.0, 1.0, 1.0, 1.0, 1.0, 1.0
};
/* ----------- initialize data --------------------------------- */
SurfPtr = Surface[Comp];
/* silence compiler warning */
(void) k_35_low; (void) k_40_low; (void) k_45_low; (void) k_50_low;
(void) h_35_high; (void) h_40_high; (void) h_45_high; (void) h_50_high;
/*
* scale factor for surface area
*/
ScaleFactor = SQR(OrbPtr->Dist);
/*
* model atmospheres - log_g switching logic outside loop
*/
if (fabs(Binary[Comp].log_g - 3.5) < FLT_EPSILON) {
maxnum_h = h_35_max; maxnum_k = k_35_max;
tarr_h = h_temp; tarr_k = k_temp_35;
h_low = h_35_low; k_high = k_35_high;
adjust_bb_lo = adjust_bb_lo_35;
adjust_bb_hi = adjust_bb_hi_35;
adjust_k = adjust_flux_35;
for (j = 0; j < NUM_MAG; ++j)
{
farr_k[j] = k_flux_35[j];
farr_h[j] = h_flux_35[j];
}
} else if (fabs(Binary[Comp].log_g - 4.0) < FLT_EPSILON) {
maxnum_h = h_40_max; maxnum_k = k_40_max;
tarr_h = h_temp; tarr_k = k_temp_40;
h_low = h_40_low; k_high = k_40_high;
adjust_bb_lo = adjust_bb_lo_40;
adjust_bb_hi = adjust_bb_hi_40;
adjust_k = adjust_flux_40;
for (j = 0; j < NUM_MAG; ++j)
{
farr_k[j] = k_flux_40[j];
farr_h[j] = h_flux_40[j];
}
} else if (fabs(Binary[Comp].log_g - 4.5) < FLT_EPSILON) {
maxnum_h = h_45_max; maxnum_k = k_45_max;
tarr_h = h_temp; tarr_k = k_temp_45;
h_low = h_45_low; k_high = k_45_high;
adjust_bb_lo = adjust_bb_lo_45;
adjust_bb_hi = adjust_bb_hi_45;
adjust_k = adjust_flux_45;
for (j = 0; j < NUM_MAG; ++j)
{
farr_k[j] = k_flux_45[j];
farr_h[j] = h_flux_45[j];
}
} else if (fabs(Binary[Comp].log_g - 5.0) < FLT_EPSILON) {
maxnum_h = h_50_max; maxnum_k = k_50_max;
tarr_h = h_temp; tarr_k = k_temp_50;
h_low = h_50_low; k_high = k_50_high;
adjust_bb_lo = adjust_bb_lo_50;
adjust_bb_hi = adjust_bb_hi_50;
adjust_k = adjust_flux_50;
for (j = 0; j < NUM_MAG; ++j)
{
farr_k[j] = k_flux_50[j];
farr_h[j] = h_flux_50[j];
}
} else {
WARNING(_("log g out of range"));
return;
}
last = 0;
/* silence compiler warning */
maxnum = maxnum_h;
tarr = tarr_h;
adjust = adjust_unity;
/* ----------- main loop --------------------------------------- */
for (i = 0; i < Binary[Comp].NumElem; ++i) {
T = SurfPtr->temp;
if (T >= h_low && T <= SWITCH_TEMP)
{
if (last == 1)
goto IsOK;
maxnum = maxnum_h;
tarr = tarr_h;
adjust = adjust_unity;
for (j = 0; j < NUM_MAG; ++j)
farr[j] = farr_h[j];
last = 1;
}
else if (T > SWITCH_TEMP && T <= k_high )
{
if (last == 2)
goto IsOK;
maxnum = maxnum_k;
tarr = tarr_k;
adjust = adjust_k;
for (j = 0; j < NUM_MAG; ++j)
farr[j] = farr_k[j];
last = 2;
}
else
{
BBFlux(SurfPtr, Comp);
if (T > SWITCH_TEMP)
for (j = 0; j < NUM_MAG; ++j)
{
Cubic = SurfPtr->f_[j] / adjust_bb_hi[j];
Cubic *= adjust_k[j];
SurfPtr->f_[j] = Cubic;
}
else
for (j = 0; j < NUM_MAG; ++j)
SurfPtr->f_[j] /= adjust_bb_lo[j];
++SurfPtr;
continue;
}
IsOK:
ScaledArea = SurfPtr->area * ScaleFactor;
/* ----------- get hashkey --------------------------------------- */
hkey = search_key (T, tarr, maxnum);
if (hkey == 0 ) { /* parabolic interpolation */
for (band=0; band < NUM_MAG; ++band) {
FLo = farr[band][hkey];
FMid = farr[band][hkey+1];
FHi = farr[band][hkey+2];
Phi = T;
Lo = (double) tarr[hkey];
Mid = (double) tarr[hkey+1];
Hi = (double) tarr[hkey+2];
/* ------------- do a parabolic intepolation -----------------
* (uses Lagrange interpolation polynome)
*/
Parabolic = FLo * ((Phi-Mid)*(Phi-Hi )) / ((Lo -Mid)*(Lo -Hi ))
+ FMid * ((Phi-Lo )*(Phi-Hi )) / ((Mid-Lo )*(Mid-Hi ))
+ FHi * ((Phi-Lo )*(Phi-Mid)) / ((Hi -Lo )*(Hi -Mid));
SurfPtr->f_[band] = (ScaledArea * Parabolic) * adjust[band];
}
} else { /* cubic interpolation */
for (band=0; band < NUM_MAG; ++band) {
FLLo = farr[band][hkey];
FLo = farr[band][hkey+1];
FMid = farr[band][hkey+2];
FHi = farr[band][hkey+3];
Phi = T;
LLo = (double) tarr[hkey];
Lo = (double) tarr[hkey+1];
Mid = (double) tarr[hkey+2];
Hi = (double) tarr[hkey+3];
/* do a cubic intepolation
* (uses Lagrange interpolation polynome)
*/
Cubic =
FLLo* ((Phi-Lo )*(Phi-Mid)*(Phi-Hi )) /
((LLo- Lo)*(LLo-Mid)*(LLo-Hi ))
+ FLo * ((Phi-LLo)*(Phi-Mid)*(Phi-Hi )) /
((Lo -LLo)*(Lo -Mid)*(Lo -Hi ))
+ FMid* ((Phi-LLo)*(Phi-Lo )*(Phi-Hi )) /
((Mid-LLo)*(Mid-Lo )*(Mid-Hi ))
+ FHi * ((Phi-LLo)*(Phi-Lo )*(Phi-Mid)) /
((Hi -LLo)*(Hi -Lo )*(Hi -Mid));
SurfPtr->f_[band] = (ScaledArea * Cubic) * adjust[band];
}
}
++SurfPtr;
}
return;
}
syntax highlighted by Code2HTML, v. 0.9.1