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