cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
stars.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 #include "cddefines.h"
4 #include "physconst.h"
5 #include "optimize.h"
6 #include "continuum.h"
7 #include "called.h"
8 #include "rfield.h"
9 #include "stars.h"
10 /*lint -e785 too few initializers */
11 /*lint -e801 use of go to depreciated */
12 
14 static const int NSB99 = 1250;
16 static const int MNTS = 200;
17 
19 static const int NRAUCH = 19951;
21 static const int NMODS_HCA = 66;
23 static const int NMODS_HNI = 51;
25 static const int NMODS_PG1159 = 71;
27 static const int NMODS_HYDR = 100;
29 static const int NMODS_HELIUM = 81;
31 static const int NMODS_HpHE = 117;
32 
33 /* set to 1 to turn on debug print statements in these routines */
34 #define DEBUGPRT 0
35 
36 #define FREE_CHECK(PTR) { ASSERT( PTR != NULL ); free( PTR ); PTR = NULL; }
37 #define FREE_SAFE(PTR) { if( PTR != NULL ) free( PTR ); PTR = NULL; }
38 
39 static const bool lgSILENT = false;
40 static const bool lgVERBOSE = true;
41 
42 static const bool lgLINEAR = false;
43 static const bool lgTAKELOG = true;
44 
45 typedef enum {
47 } IntStage;
48 
50 typedef struct
51 {
52  double par[MDIM];
53  int modid;
54  char chGrid;
55 } mpp;
56 
66 /* this is the structure of the binary atmosphere file (VERSION 20060612):
67  *
68  * ============================
69  * * int32 VERSION *
70  * * int32 MDIM *
71  * * int32 MNAM *
72  * * int32 ndim *
73  * * int32 npar *
74  * * int32 nmods *
75  * * int32 ngrid *
76  * * uint32 nOffset *
77  * * uint32 nBlocksize *
78  * * char names[MDIM][MNAM+1] *
79  * * mpp telg[nmods] *
80  * * realnum anu[ngrid] *
81  * * realnum mod1[ngrid] *
82  * * ... *
83  * * realnum modn[ngrid] *
84  * ============================
85  *
86  * nOffset == 7*sizeof(int32) + 2*sizeof(uint32) + MDIM*(MNAM+1)*sizeof(char) + nmods*sizeof(mpp)
87  * nBlocksize == ngrid*size(realnum) */
88 
90 typedef struct
91 {
93  string name;
97  FILE *ioIN;
100  const char *ident;
102  const char *command;
106  int32 ndim;
108  int32 npar;
110  int32 nmods;
112  int32 ngrid;
114  uint32 nOffset;
116  uint32 nBlocksize;
119  mpp *telg; /* telg[nmods] */
121  double **val; /* val[ndim][nval[n]] */
123  long *nval; /* nval[ndim] */
130  long *jlo; /* jlo(nval[0],...,nval[ndim-1]) */
131  long *jhi; /* jhi(nval[0],...,nval[ndim-1]) */
133  char names[MDIM][MNAM+1];
135  long *trackLen; /* trackLen[nTracks] */
137  long nTracks;
139  long *jval;
140 } stellar_grid;
141 
142 /* internal routines */
143 STATIC bool lgCompileAtmosphereCoStar(const char[],const char[],const realnum[],long,process_counter&);
144 STATIC void InterpolateGridCoStar(const stellar_grid*,const double[],double*,double*);
145 STATIC void FindHCoStar(const stellar_grid*,long,double,long,realnum*,long*,long*);
146 STATIC void FindVCoStar(const stellar_grid*,double,realnum*,long[]);
148 STATIC int RauchInitializeSub(const char[],const char[],const mpp[],long,long,
149  long,const double[],int);
150 STATIC bool lgCompileAtmosphere(const char[],const char[],const realnum[],long,process_counter&);
151 STATIC void InitGrid(stellar_grid*,bool);
153 STATIC bool lgValidAsciiFile(const char*,access_scheme);
155 STATIC void CheckVal(const stellar_grid*,double[],long*,long*);
156 STATIC void InterpolateRectGrid(const stellar_grid*,const double[],double*,double*);
158 STATIC void InterpolateModel(const stellar_grid*,const double[],double[],const long[],
159  const long[],long[],long,realnum[],IntStage);
160 STATIC void InterpolateModelCoStar(const stellar_grid*,const double[],double[],const long[],
161  const long[],long[],long,long,realnum[]);
162 STATIC void GetModel(const stellar_grid*,long,realnum[],bool,bool);
163 STATIC void SetLimits(const stellar_grid*,double,const long[],const long[],const long[],
164  const realnum[],double*,double*);
165 STATIC void SetLimitsSub(const stellar_grid*,double,const long[],const long[],long[],long,
166  double*,double*);
168 STATIC void FillJ(const stellar_grid*,long[],double[],long,bool);
169 STATIC long JIndex(const stellar_grid*,const long[]);
170 STATIC void SearchModel(const mpp[],long,const double[],long,long*,long*);
171 STATIC void FindIndex(const double[],long,double,long*,long*,bool*);
173 STATIC void ValidateGrid(const stellar_grid*,double);
174 STATIC bool lgValidModel(const realnum[],const realnum[],double,double);
175 STATIC void RebinAtmosphere(long,const realnum[],const realnum[],realnum[],long,const realnum[]);
176 STATIC realnum RebinSingleCell(realnum,realnum,const realnum[],const realnum[],const realnum[],long);
177 STATIC long RebinFind(const realnum[],long,realnum);
178 
179 
180 /* the version number for the ascii/binary atmosphere files */
181 static const long int VERSION_ASCII = 20060612L;
182 /* binary files are incompatible when floats are converted to doubles */
183 #ifdef FLT_IS_DBL
184 static const long int VERSION_BIN = 20071004L;
185 #else
186 static const long int VERSION_BIN = 20060612L;
187 #endif
188 
190 void AtmospheresAvail( void )
191 {
192  DEBUG_ENTRY( "AtmospheresAvail()" );
193 
194  /* This routine makes a list of all the stellar atmosphere grids that are valid,
195  * giving the parameters for use in the input script as well. It is simply a long
196  * list of if-statements, so if any grid is added to Cloudy, it should be added in
197  * this routine as well.
198  *
199  * NB NB NB -- test this routine regularly to see if the list is still complete! */
200 
201  fprintf( ioQQQ, "\n I will now list all stellar atmosphere grids that are ready to be used (if any).\n" );
202  fprintf( ioQQQ, " User-defined stellar atmosphere grids will not be included in this list.\n\n" );
203 
204  process_counter dum;
205 
206  /* we always look in the data directory regardless of where we are,
207  * it would be very confusing to the user if we did otherwise... */
209 
210  if( lgValidBinFile( "atlas_fp10k2.mod", dum, as ) )
211  fprintf( ioQQQ, " table star atlas Z+1.0 <Teff> [ <log(g)> ]\n" );
212  if( lgValidBinFile( "atlas_fp05k2.mod", dum, as ) )
213  fprintf( ioQQQ, " table star atlas Z+0.5 <Teff> [ <log(g)> ]\n" );
214  if( lgValidBinFile( "atlas_fp03k2.mod", dum, as ) )
215  fprintf( ioQQQ, " table star atlas Z+0.3 <Teff> [ <log(g)> ]\n" );
216  if( lgValidBinFile( "atlas_fp02k2.mod", dum, as ) )
217  fprintf( ioQQQ, " table star atlas Z+0.2 <Teff> [ <log(g)> ]\n" );
218  if( lgValidBinFile( "atlas_fp01k2.mod", dum, as ) )
219  fprintf( ioQQQ, " table star atlas Z+0.1 <Teff> [ <log(g)> ]\n" );
220  if( lgValidBinFile( "atlas_fp00k2.mod", dum, as ) )
221  fprintf( ioQQQ, " table star atlas Z+0.0 <Teff> [ <log(g)> ]\n" );
222  if( lgValidBinFile( "atlas_fm01k2.mod", dum, as ) )
223  fprintf( ioQQQ, " table star atlas Z-0.1 <Teff> [ <log(g)> ]\n" );
224  if( lgValidBinFile( "atlas_fm02k2.mod", dum, as ) )
225  fprintf( ioQQQ, " table star atlas Z-0.2 <Teff> [ <log(g)> ]\n" );
226  if( lgValidBinFile( "atlas_fm03k2.mod", dum, as ) )
227  fprintf( ioQQQ, " table star atlas Z-0.3 <Teff> [ <log(g)> ]\n" );
228  if( lgValidBinFile( "atlas_fm05k2.mod", dum, as ) )
229  fprintf( ioQQQ, " table star atlas Z-0.5 <Teff> [ <log(g)> ]\n" );
230  if( lgValidBinFile( "atlas_fm10k2.mod", dum, as ) )
231  fprintf( ioQQQ, " table star atlas Z-1.0 <Teff> [ <log(g)> ]\n" );
232  if( lgValidBinFile( "atlas_fm15k2.mod", dum, as ) )
233  fprintf( ioQQQ, " table star atlas Z-1.5 <Teff> [ <log(g)> ]\n" );
234  if( lgValidBinFile( "atlas_fm20k2.mod", dum, as ) )
235  fprintf( ioQQQ, " table star atlas Z-2.0 <Teff> [ <log(g)> ]\n" );
236  if( lgValidBinFile( "atlas_fm25k2.mod", dum, as ) )
237  fprintf( ioQQQ, " table star atlas Z-2.5 <Teff> [ <log(g)> ]\n" );
238  if( lgValidBinFile( "atlas_fm30k2.mod", dum, as ) )
239  fprintf( ioQQQ, " table star atlas Z-3.0 <Teff> [ <log(g)> ]\n" );
240  if( lgValidBinFile( "atlas_fm35k2.mod", dum, as ) )
241  fprintf( ioQQQ, " table star atlas Z-3.5 <Teff> [ <log(g)> ]\n" );
242  if( lgValidBinFile( "atlas_fm40k2.mod", dum, as ) )
243  fprintf( ioQQQ, " table star atlas Z-4.0 <Teff> [ <log(g)> ]\n" );
244  if( lgValidBinFile( "atlas_fm45k2.mod", dum, as ) )
245  fprintf( ioQQQ, " table star atlas Z-4.5 <Teff> [ <log(g)> ]\n" );
246  if( lgValidBinFile( "atlas_fm50k2.mod", dum, as ) )
247  fprintf( ioQQQ, " table star atlas Z-5.0 <Teff> [ <log(g)> ]\n" );
248 
249  if( lgValidBinFile( "atlas_fp05k2_odfnew.mod", dum, as ) )
250  fprintf( ioQQQ, " table star atlas odfnew Z+0.5 <Teff> [ <log(g)> ]\n" );
251  if( lgValidBinFile( "atlas_fp02k2_odfnew.mod", dum, as ) )
252  fprintf( ioQQQ, " table star atlas odfnew Z+0.2 <Teff> [ <log(g)> ]\n" );
253  if( lgValidBinFile( "atlas_fp00k2_odfnew.mod", dum, as ) )
254  fprintf( ioQQQ, " table star atlas odfnew Z+0.0 <Teff> [ <log(g)> ]\n" );
255  if( lgValidBinFile( "atlas_fm05k2_odfnew.mod", dum, as ) )
256  fprintf( ioQQQ, " table star atlas odfnew Z-0.5 <Teff> [ <log(g)> ]\n" );
257  if( lgValidBinFile( "atlas_fm10k2_odfnew.mod", dum, as ) )
258  fprintf( ioQQQ, " table star atlas odfnew Z-1.0 <Teff> [ <log(g)> ]\n" );
259  if( lgValidBinFile( "atlas_fm15k2_odfnew.mod", dum, as ) )
260  fprintf( ioQQQ, " table star atlas odfnew Z-1.5 <Teff> [ <log(g)> ]\n" );
261  if( lgValidBinFile( "atlas_fm20k2_odfnew.mod", dum, as ) )
262  fprintf( ioQQQ, " table star atlas odfnew Z-2.0 <Teff> [ <log(g)> ]\n" );
263  if( lgValidBinFile( "atlas_fm25k2_odfnew.mod", dum, as ) )
264  fprintf( ioQQQ, " table star atlas odfnew Z-2.5 <Teff> [ <log(g)> ]\n" );
265 
266  if( lgValidBinFile( "atlas_3d.mod", dum, as ) )
267  fprintf( ioQQQ, " table star atlas 3-dim <Teff> <log(g)> <log(Z)>\n" );
268 
269  if( lgValidBinFile( "atlas_3d_odfnew.mod", dum, as ) )
270  fprintf( ioQQQ, " table star atlas odfnew 3-dim <Teff> <log(g)> <log(Z)>\n" );
271 
272  if( lgValidBinFile( "Sc1_costar_solar.mod", dum, as ) )
273  fprintf( ioQQQ, " table star costar solar (see Hazy for parameters)\n" );
274  if( lgValidBinFile( "Sc1_costar_halo.mod", dum, as ) )
275  fprintf( ioQQQ, " table star costar halo (see Hazy for parameters)\n" );
276 
277  if( lgValidBinFile( "kurucz79.mod", dum, as ) )
278  fprintf( ioQQQ, " table star kurucz79 <Teff>\n" );
279 
280  if( lgValidBinFile( "mihalas.mod", dum, as ) )
281  fprintf( ioQQQ, " table star mihalas <Teff>\n" );
282 
283  if( lgValidBinFile( "rauch_h-ca_solar.mod", dum, as ) )
284  fprintf( ioQQQ, " table star rauch H-Ca solar <Teff> [ <log(g)> ]\n" );
285  if( lgValidBinFile( "rauch_h-ca_halo.mod", dum, as ) )
286  fprintf( ioQQQ, " table star rauch H-Ca halo <Teff> [ <log(g)> ]\n" );
287  if( lgValidBinFile( "rauch_h-ca_3d.mod", dum, as ) )
288  fprintf( ioQQQ, " table star rauch H-Ca 3-dim <Teff> <log(g)> <log(Z)>\n" );
289 
290  if( lgValidBinFile( "rauch_h-ni_solar.mod", dum, as ) )
291  fprintf( ioQQQ, " table star rauch H-Ni solar <Teff> [ <log(g)> ]\n" );
292  if( lgValidBinFile( "rauch_h-ni_halo.mod", dum, as ) )
293  fprintf( ioQQQ, " table star rauch H-Ni halo <Teff> [ <log(g)> ]\n" );
294  if( lgValidBinFile( "rauch_h-ni_3d.mod", dum, as ) )
295  fprintf( ioQQQ, " table star rauch H-Ni 3-dim <Teff> <log(g)> <log(Z)>\n" );
296 
297  if( lgValidBinFile( "rauch_pg1159.mod", dum, as ) )
298  fprintf( ioQQQ, " table star rauch pg1159 <Teff> [ <log(g)> ]\n" );
299 
300  if( lgValidBinFile( "rauch_hydr.mod", dum, as ) )
301  fprintf( ioQQQ, " table star rauch hydrogen <Teff> [ <log(g)> ]\n" );
302 
303  if( lgValidBinFile( "rauch_helium.mod", dum, as ) )
304  fprintf( ioQQQ, " table star rauch helium <Teff> [ <log(g)> ]\n" );
305 
306  if( lgValidBinFile( "rauch_h+he_3d.mod", dum, as ) )
307  fprintf( ioQQQ, " table star rauch H+He <Teff> <log(g)> <frac(He)>\n" );
308 
309  if( lgValidBinFile( "starburst99.mod", dum, as ) )
310  fprintf( ioQQQ, " table star \"starburst99.mod\" <age>\n" );
311 
312  if( lgValidBinFile( "bstar2006_p03.mod", dum, as ) )
313  fprintf( ioQQQ, " table star tlusty Bstar Z+0.3 <Teff> [ <log(g)> ]\n" );
314  if( lgValidBinFile( "bstar2006_p00.mod", dum, as ) )
315  fprintf( ioQQQ, " table star tlusty Bstar Z+0.0 <Teff> [ <log(g)> ]\n" );
316  if( lgValidBinFile( "bstar2006_m03.mod", dum, as ) )
317  fprintf( ioQQQ, " table star tlusty Bstar Z-0.3 <Teff> [ <log(g)> ]\n" );
318  if( lgValidBinFile( "bstar2006_m07.mod", dum, as ) )
319  fprintf( ioQQQ, " table star tlusty Bstar Z-0.7 <Teff> [ <log(g)> ]\n" );
320  if( lgValidBinFile( "bstar2006_m10.mod", dum, as ) )
321  fprintf( ioQQQ, " table star tlusty Bstar Z-1.0 <Teff> [ <log(g)> ]\n" );
322  if( lgValidBinFile( "bstar2006_m99.mod", dum, as ) )
323  fprintf( ioQQQ, " table star tlusty Bstar Z-inf <Teff> [ <log(g)> ]\n" );
324 
325  if( lgValidBinFile( "bstar2006_3d.mod", dum, as ) )
326  fprintf( ioQQQ, " table star tlusty Bstar 3-dim <Teff> <log(g)> <log(Z)>\n" );
327 
328  if( lgValidBinFile( "ostar2002_p03.mod", dum, as ) )
329  fprintf( ioQQQ, " table star tlusty Ostar Z+0.3 <Teff> [ <log(g)> ]\n" );
330  if( lgValidBinFile( "ostar2002_p00.mod", dum, as ) )
331  fprintf( ioQQQ, " table star tlusty Ostar Z+0.0 <Teff> [ <log(g)> ]\n" );
332  if( lgValidBinFile( "ostar2002_m03.mod", dum, as ) )
333  fprintf( ioQQQ, " table star tlusty Ostar Z-0.3 <Teff> [ <log(g)> ]\n" );
334  if( lgValidBinFile( "ostar2002_m07.mod", dum, as ) )
335  fprintf( ioQQQ, " table star tlusty Ostar Z-0.7 <Teff> [ <log(g)> ]\n" );
336  if( lgValidBinFile( "ostar2002_m10.mod", dum, as ) )
337  fprintf( ioQQQ, " table star tlusty Ostar Z-1.0 <Teff> [ <log(g)> ]\n" );
338  if( lgValidBinFile( "ostar2002_m15.mod", dum, as ) )
339  fprintf( ioQQQ, " table star tlusty Ostar Z-1.5 <Teff> [ <log(g)> ]\n" );
340  if( lgValidBinFile( "ostar2002_m17.mod", dum, as ) )
341  fprintf( ioQQQ, " table star tlusty Ostar Z-1.7 <Teff> [ <log(g)> ]\n" );
342  if( lgValidBinFile( "ostar2002_m20.mod", dum, as ) )
343  fprintf( ioQQQ, " table star tlusty Ostar Z-2.0 <Teff> [ <log(g)> ]\n" );
344  if( lgValidBinFile( "ostar2002_m30.mod", dum, as ) )
345  fprintf( ioQQQ, " table star tlusty Ostar Z-3.0 <Teff> [ <log(g)> ]\n" );
346  if( lgValidBinFile( "ostar2002_m99.mod", dum, as ) )
347  fprintf( ioQQQ, " table star tlusty Ostar Z-inf <Teff> [ <log(g)> ]\n" );
348 
349  if( lgValidBinFile( "ostar2002_3d.mod", dum, as ) )
350  fprintf( ioQQQ, " table star tlusty Ostar 3-dim <Teff> <log(g)> <log(Z)>\n" );
351 
352  if( lgValidBinFile( "kwerner.mod", dum, as ) )
353  fprintf( ioQQQ, " table star werner <Teff> [ <log(g)> ]\n" );
354 
355  if( lgValidBinFile( "wmbasic.mod", dum, as ) )
356  fprintf( ioQQQ, " table star wmbasic <Teff> <log(g)> <log(Z)>\n" );
357  return;
358 }
359 
360 /* AtlasCompile rebin Kurucz stellar models to match energy grid of code */
361 /* >>chng 05 nov 16, added return value to indicate success (0) or failure (1) */
363 {
364  /* these contain frequencies for the major absorption edges */
365  realnum Edges[3];
366 
367  bool lgFail = false;
368 
369  DEBUG_ENTRY( "AtlasCompile()" );
370 
371  /* This is a program to re-bin the Kurucz stellar models spectrum to match the
372  * CLOUDY grid. For wavelengths shorter than supplied in the Kurucz files,
373  * the flux will be set to zero. At long wavelengths a Rayleigh-Jeans
374  * extrapolation will be used. */
375 
376  /* This version uses power-law interpolation between the points of the stellar
377  * model.*/
378 
379  fprintf( ioQQQ, " AtlasCompile on the job.\n" );
380 
381  /* define the major absorption edges that require special attention during rebinning
382  *
383  * NB the frequencies should be chosen here such that they are somewhere in between
384  * the two frequency points that straddle the edge in the atmosphere model, the
385  * software in RebinAtmosphere will seek out the exact values of those two points
386  * e.g.: in the CoStar models the H I edge is straddled by wavelength points at
387  * 911.67 and 911.85 A, so Edges[0] should be chosen somewhere in between (e.g. at 911.76A).
388  *
389  * NB beware not to choose edges too close to one another (i.e. on the order of the
390  * resolution of the Cloudy frequency grid). E.g. the He II Balmer edge nearly coincides
391  * with the H I Ly edge, they should be treated as one edge. Trying to separate them will
392  * almost certainly lead to erroneous behaviour in RebinAtmosphere */
393  Edges[0] = (realnum)(RYDLAM/911.76);
394  Edges[1] = (realnum)(RYDLAM/504.26);
395  Edges[2] = (realnum)(RYDLAM/227.84);
396 
398 
399  /* >>chng 05 nov 19, add support for non-solar metalicities as well as odfnew models, PvH */
400  if( lgFileReadable( "atlas_fp10k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fp10k2.mod", pc, as ) )
401  lgFail = lgFail || lgCompileAtmosphere( "atlas_fp10k2.ascii", "atlas_fp10k2.mod", Edges, 3L, pc );
402  if( lgFileReadable( "atlas_fp05k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fp05k2.mod", pc, as ) )
403  lgFail = lgFail || lgCompileAtmosphere( "atlas_fp05k2.ascii", "atlas_fp05k2.mod", Edges, 3L, pc );
404  if( lgFileReadable( "atlas_fp03k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fp03k2.mod", pc, as ) )
405  lgFail = lgFail || lgCompileAtmosphere( "atlas_fp03k2.ascii", "atlas_fp03k2.mod", Edges, 3L, pc );
406  if( lgFileReadable( "atlas_fp02k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fp02k2.mod", pc, as ) )
407  lgFail = lgFail || lgCompileAtmosphere( "atlas_fp02k2.ascii", "atlas_fp02k2.mod", Edges, 3L, pc );
408  if( lgFileReadable( "atlas_fp01k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fp01k2.mod", pc, as ) )
409  lgFail = lgFail || lgCompileAtmosphere( "atlas_fp01k2.ascii", "atlas_fp01k2.mod", Edges, 3L, pc );
410  if( lgFileReadable( "atlas_fp00k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fp00k2.mod", pc, as ) )
411  lgFail = lgFail || lgCompileAtmosphere( "atlas_fp00k2.ascii", "atlas_fp00k2.mod", Edges, 3L, pc );
412  if( lgFileReadable( "atlas_fm01k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fm01k2.mod", pc, as ) )
413  lgFail = lgFail || lgCompileAtmosphere( "atlas_fm01k2.ascii", "atlas_fm01k2.mod", Edges, 3L, pc );
414  if( lgFileReadable( "atlas_fm02k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fm02k2.mod", pc, as ) )
415  lgFail = lgFail || lgCompileAtmosphere( "atlas_fm02k2.ascii", "atlas_fm02k2.mod", Edges, 3L, pc );
416  if( lgFileReadable( "atlas_fm03k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fm03k2.mod", pc, as ) )
417  lgFail = lgFail || lgCompileAtmosphere( "atlas_fm03k2.ascii", "atlas_fm03k2.mod", Edges, 3L, pc );
418  if( lgFileReadable( "atlas_fm05k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fm05k2.mod", pc, as ) )
419  lgFail = lgFail || lgCompileAtmosphere( "atlas_fm05k2.ascii", "atlas_fm05k2.mod", Edges, 3L, pc );
420  if( lgFileReadable( "atlas_fm10k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fm10k2.mod", pc, as ) )
421  lgFail = lgFail || lgCompileAtmosphere( "atlas_fm10k2.ascii", "atlas_fm10k2.mod", Edges, 3L, pc );
422  if( lgFileReadable( "atlas_fm15k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fm15k2.mod", pc, as ) )
423  lgFail = lgFail || lgCompileAtmosphere( "atlas_fm15k2.ascii", "atlas_fm15k2.mod", Edges, 3L, pc );
424  if( lgFileReadable( "atlas_fm20k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fm20k2.mod", pc, as ) )
425  lgFail = lgFail || lgCompileAtmosphere( "atlas_fm20k2.ascii", "atlas_fm20k2.mod", Edges, 3L, pc );
426  if( lgFileReadable( "atlas_fm25k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fm25k2.mod", pc, as ) )
427  lgFail = lgFail || lgCompileAtmosphere( "atlas_fm25k2.ascii", "atlas_fm25k2.mod", Edges, 3L, pc );
428  if( lgFileReadable( "atlas_fm30k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fm30k2.mod", pc, as ) )
429  lgFail = lgFail || lgCompileAtmosphere( "atlas_fm30k2.ascii", "atlas_fm30k2.mod", Edges, 3L, pc );
430  if( lgFileReadable( "atlas_fm35k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fm35k2.mod", pc, as ) )
431  lgFail = lgFail || lgCompileAtmosphere( "atlas_fm35k2.ascii", "atlas_fm35k2.mod", Edges, 3L, pc );
432  if( lgFileReadable( "atlas_fm40k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fm40k2.mod", pc, as ) )
433  lgFail = lgFail || lgCompileAtmosphere( "atlas_fm40k2.ascii", "atlas_fm40k2.mod", Edges, 3L, pc );
434  if( lgFileReadable( "atlas_fm45k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fm45k2.mod", pc, as ) )
435  lgFail = lgFail || lgCompileAtmosphere( "atlas_fm45k2.ascii", "atlas_fm45k2.mod", Edges, 3L, pc );
436  if( lgFileReadable( "atlas_fm50k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fm50k2.mod", pc, as ) )
437  lgFail = lgFail || lgCompileAtmosphere( "atlas_fm50k2.ascii", "atlas_fm50k2.mod", Edges, 3L, pc );
438 
439  if( lgFileReadable( "atlas_fp05k2_odfnew.ascii", pc, as ) &&
440  !lgValidBinFile( "atlas_fp05k2_odfnew.mod", pc, as ) )
441 
442  lgFail = lgFail || lgCompileAtmosphere( "atlas_fp05k2_odfnew.ascii", "atlas_fp05k2_odfnew.mod",
443  Edges, 3L, pc );
444  if( lgFileReadable( "atlas_fp02k2_odfnew.ascii", pc, as ) &&
445  !lgValidBinFile( "atlas_fp02k2_odfnew.mod", pc, as ) )
446  lgFail = lgFail || lgCompileAtmosphere( "atlas_fp02k2_odfnew.ascii", "atlas_fp02k2_odfnew.mod",
447  Edges, 3L, pc );
448  if( lgFileReadable( "atlas_fp00k2_odfnew.ascii", pc, as ) &&
449  !lgValidBinFile( "atlas_fp00k2_odfnew.mod", pc, as ) )
450  lgFail = lgFail || lgCompileAtmosphere( "atlas_fp00k2_odfnew.ascii", "atlas_fp00k2_odfnew.mod",
451  Edges, 3L, pc );
452  if( lgFileReadable( "atlas_fm05k2_odfnew.ascii", pc, as ) &&
453  !lgValidBinFile( "atlas_fm05k2_odfnew.mod", pc, as ) )
454  lgFail = lgFail || lgCompileAtmosphere( "atlas_fm05k2_odfnew.ascii", "atlas_fm05k2_odfnew.mod",
455  Edges, 3L, pc );
456  if( lgFileReadable( "atlas_fm10k2_odfnew.ascii", pc, as ) &&
457  !lgValidBinFile( "atlas_fm10k2_odfnew.mod", pc, as ) )
458  lgFail = lgFail || lgCompileAtmosphere( "atlas_fm10k2_odfnew.ascii", "atlas_fm10k2_odfnew.mod",
459  Edges, 3L, pc );
460  if( lgFileReadable( "atlas_fm15k2_odfnew.ascii", pc, as ) &&
461  !lgValidBinFile( "atlas_fm15k2_odfnew.mod", pc, as ) )
462  lgFail = lgFail || lgCompileAtmosphere( "atlas_fm15k2_odfnew.ascii", "atlas_fm15k2_odfnew.mod",
463  Edges, 3L, pc );
464  if( lgFileReadable( "atlas_fm20k2_odfnew.ascii", pc, as ) &&
465  !lgValidBinFile( "atlas_fm20k2_odfnew.mod", pc, as ) )
466  lgFail = lgFail || lgCompileAtmosphere( "atlas_fm20k2_odfnew.ascii", "atlas_fm20k2_odfnew.mod",
467  Edges, 3L, pc );
468  if( lgFileReadable( "atlas_fm25k2_odfnew.ascii", pc, as ) &&
469  !lgValidBinFile( "atlas_fm25k2_odfnew.mod", pc, as ) )
470  lgFail = lgFail || lgCompileAtmosphere( "atlas_fm25k2_odfnew.ascii", "atlas_fm25k2_odfnew.mod",
471  Edges, 3L, pc );
472 
473  if( lgFileReadable( "atlas_3d.ascii", pc, as ) && !lgValidBinFile( "atlas_3d.mod", pc, as ) )
474  lgFail = lgFail || lgCompileAtmosphere( "atlas_3d.ascii", "atlas_3d.mod", Edges, 3L, pc );
475 
476  if( lgFileReadable( "atlas_3d_odfnew.ascii", pc, as ) &&
477  !lgValidBinFile( "atlas_3d_odfnew.mod", pc, as ) )
478  lgFail = lgFail || lgCompileAtmosphere( "atlas_3d_odfnew.ascii", "atlas_3d_odfnew.mod", Edges, 3L, pc );
479  return lgFail;
480 }
481 
482 /* AtlasInterpolate read in and interpolate on Kurucz grid of atmospheres, originally by K Volk */
483 long AtlasInterpolate(double val[], /* val[nval] */
484  long *nval,
485  long *ndim,
486  const char *chMetalicity,
487  const char *chODFNew,
488  bool lgList,
489  double *Tlow,
490  double *Thigh)
491 {
492  char chIdent[13];
494 
495  DEBUG_ENTRY( "AtlasInterpolate()" );
496 
497  grid.name = "atlas_";
498  if( *ndim == 3 )
499  grid.name += "3d";
500  else
501  {
502  grid.name += "f";
503  grid.name += chMetalicity;
504  grid.name += "k2";
505  }
506  grid.name += chODFNew;
507  grid.name += ".mod";
508  grid.scheme = AS_DATA_OPTIONAL;
509  /* identification of this atmosphere set, used in
510  * the Cloudy output, *must* be 12 characters long */
511  if( *ndim == 3 )
512  {
513  strcpy( chIdent, "3-dim" );
514  }
515  else
516  {
517  strcpy( chIdent, "Z " );
518  strcat( chIdent, chMetalicity );
519  }
520  strcat( chIdent, ( strlen(chODFNew) == 0 ? " Kurucz" : " ODFNew" ) );
521  grid.ident = chIdent;
522  /* the Cloudy command needed to recompile the binary model file */
523  grid.command = "COMPILE STARS";
524 
525  InitGrid( &grid, lgList );
526 
527  CheckVal( &grid, val, nval, ndim );
528 
529  /* Note on the interpolation (solar abundance grid): 26 October 2000 (Peter van Hoof)
530  *
531  * I computed the effective temperature for a random sample of interpolated
532  * atmospheres by integrating the flux as shown above and compared the results
533  * with the expected effective temperature using DELTA = (COMP-EXPEC)/EXPEC.
534  *
535  * I found that the average discrepancy was:
536  *
537  * DELTA = -0.10% +/- 0.06% (sample size 5000)
538  *
539  * The most extreme discrepancies were
540  * -0.30% <= DELTA <= 0.21%
541  *
542  * The most negative discrepancies were for Teff = 36 - 39 kK, log(g) = 4.5 - 5
543  * The most positive discrepancies were for Teff = 3.5 - 4.0 kK, log(g) = 0 - 1
544  *
545  * The interpolation in the ATLAS grid is clearly very accurate */
546 
547  InterpolateRectGrid( &grid, val, Tlow, Thigh );
548 
549  FreeGrid( &grid );
550  return rfield.nupper;
551 }
552 
553 /* CoStarCompile rebin costar stellar models to match energy grid of code*/
555 {
556  realnum Edges[3];
557  bool lgFail = false;
558 
559  DEBUG_ENTRY( "CoStarCompile()" );
560 
561  fprintf( ioQQQ, " CoStarCompile on the job.\n" );
562 
563  /* define the major absorption edges that require special attention during rebinning
564  *
565  * NB the frequencies should be chosen here such that they are somewhere in between
566  * the two frequency points that straddle the edge in the atmosphere model, the
567  * software in RebinAtmosphere will seek out the exact values of those two points
568  * e.g.: in the CoStar models the H I edge is straddled by wavelength points at
569  * 911.67 and 911.85 A, so Edges[0] should be chosen somewhere in between (e.g. at 911.76A).
570  *
571  * NB beware not to choose edges too close to one another (i.e. on the order of the
572  * resolution of the Cloudy frequency grid). E.g. the He II Balmer edge nearly coincides
573  * with the H I Ly edge, they should be treated as one edge. Trying to separate them will
574  * almost certainly lead to erroneous behaviour in RebinAtmosphere */
575  Edges[0] = (realnum)(RYDLAM/911.76);
576  Edges[1] = (realnum)(RYDLAM/504.26);
577  Edges[2] = (realnum)(RYDLAM/227.84);
578 
580 
581  if( lgFileReadable( "Sc1_costar_z020_lb.fluxes", pc, as ) && !lgValidBinFile( "Sc1_costar_solar.mod", pc, as ) )
582  lgFail = lgFail || lgCompileAtmosphereCoStar( "Sc1_costar_z020_lb.fluxes", "Sc1_costar_solar.mod",
583  Edges, 3L, pc );
584  if( lgFileReadable( "Sc1_costar_z004_lb.fluxes", pc, as ) && !lgValidBinFile( "Sc1_costar_halo.mod", pc, as ) )
585  lgFail = lgFail || lgCompileAtmosphereCoStar( "Sc1_costar_z004_lb.fluxes", "Sc1_costar_halo.mod",
586  Edges, 3L, pc );
587  return lgFail;
588 }
589 
590 /* CoStarInterpolate read in and interpolate on CoStar grid of atmospheres */
591 long CoStarInterpolate(double val[], /* requested model parameters */
592  long *nval,
593  long *ndim,
594  IntMode imode, /* which interpolation mode is requested */
595  bool lgHalo, /* flag indicating whether solar (==0) or halo (==1) abundances */
596  bool lgList,
597  double *val0_lo,
598  double *val0_hi)
599 {
601 
602  DEBUG_ENTRY( "CoStarInterpolate()" );
603 
604  grid.name = ( lgHalo ? "Sc1_costar_halo.mod" : "Sc1_costar_solar.mod" );
605  grid.scheme = AS_DATA_OPTIONAL;
606  /* identification of this atmosphere set, used in
607  * the Cloudy output, *must* be 12 characters long */
608  grid.ident = " costar";
609  /* the Cloudy command needed to recompile the binary model file */
610  grid.command = "COMPILE STARS";
611 
612  /* listing the models in the grid is implemented in CoStarListModels() */
613  InitGrid( &grid, false );
614  /* now sort the models according to track */
615  InitGridCoStar( &grid );
616  /* override default interpolation mode */
617  grid.imode = imode;
618 
619  if( lgList )
620  {
621  CoStarListModels( &grid );
622  cdEXIT(EXIT_SUCCESS);
623  }
624 
625  CheckVal( &grid, val, nval, ndim );
626 
627  /* Note on the interpolation: 26 October 2000 (Peter van Hoof)
628  *
629  * I computed the effective temperature for a random sample of interpolated
630  * atmospheres by integrating the flux as shown above and compared the results
631  * with the expected effective temperature using DELTA = (COMP-EXPEC)/EXPEC.
632  *
633  * I found that the average discrepancy was:
634  *
635  * DELTA = -1.16% +/- 0.69% (SOLAR models, sample size 4590)
636  * DELTA = -1.17% +/- 0.70% (HALO models, sample size 4828)
637  *
638  * The most extreme discrepancies for the SOLAR models were
639  * -3.18% <= DELTA <= -0.16%
640  *
641  * The most negative discrepancies were for Teff = 35 kK, log(g) = 3.5
642  * The least negative discrepancies were for Teff = 50 kK, log(g) = 4.1
643  *
644  * The most extreme discrepancies for the HALO models were
645  * -2.90% <= DELTA <= -0.13%
646  *
647  * The most negative discrepancies were for Teff = 35 kK, log(g) = 3.5
648  * The least negative discrepancies were for Teff = 50 kK, log(g) = 4.1
649  *
650  * Since Cloudy checks the scaling elsewhere there is no need to re-scale
651  * things here, but this inaccuracy should be kept in mind since it could
652  * indicate problems with the flux distribution */
653 
654  InterpolateGridCoStar( &grid, val, val0_lo, val0_hi );
655 
656  FreeGrid( &grid );
657  return rfield.nupper;
658 }
659 
660 /* GridCompile rebin user supplied stellar models to match energy grid of code */
661 bool GridCompile(const char *InName)
662 {
663  bool lgFail = false;
664  realnum Edges[1];
665  string OutName( InName );
667 
668  DEBUG_ENTRY( "GridCompile()" );
669 
670  fprintf( ioQQQ, " GridCompile on the job.\n" );
671 
672  // replace filename extension with ".mod"
673  string::size_type ptr = OutName.find( '.' );
674  ASSERT( ptr != string::npos );
675  OutName.replace( ptr, string::npos, ".mod" );
676 
677  process_counter dum;
678  lgFail = lgCompileAtmosphere( InName, OutName.c_str(), Edges, 0L, dum );
679 
680  if( !lgFail )
681  {
682  /* the file must be local */
683  grid.name = OutName;
684  grid.scheme = AS_LOCAL_ONLY;
685  grid.ident = "bogus ident.";
686  grid.command = "bogus command.";
687 
688  InitGrid( &grid, false );
689 
690  /* check whether the models in the grid have the correct effective temperature */
691 
692  fprintf( ioQQQ, " GridCompile: checking effective temperatures...\n" );
693  ValidateGrid( &grid, 0.02 );
694 
695  FreeGrid( &grid );
696  }
697 
698  return lgFail;
699 }
700 
701 /* GridInterpolate read in and interpolate on user supplied grid of atmospheres */
702 long GridInterpolate(double val[], /* val[nval] */
703  long *nval,
704  long *ndim,
705  const char *FileName,
706  bool lgList,
707  double *Tlow,
708  double *Thigh)
709 {
710  char chIdent[13];
712 
713  DEBUG_ENTRY( "GridInterpolate()" );
714 
715  // make filename without extension
716  string chTruncName( FileName );
717  string::size_type ptr = chTruncName.find( '.' );
718  if( ptr != string::npos )
719  chTruncName.replace( ptr, string::npos, "" );
720 
721  grid.name = FileName;
722  grid.scheme = AS_DATA_OPTIONAL;
723  /* identification of this atmosphere set, used in
724  * the Cloudy output, *must* be 12 characters long */
725  sprintf( chIdent, "%12.12s", chTruncName.c_str() );
726  grid.ident = chIdent;
727  /* the Cloudy command needed to recompile the binary model file */
728  string chString( "COMPILE STARS \"" + chTruncName + ".ascii\"" );
729  grid.command = chString.c_str();
730 
731  InitGrid( &grid, lgList );
732 
733  CheckVal( &grid, val, nval, ndim );
734 
735  InterpolateRectGrid( &grid, val, Tlow, Thigh );
736 
737  FreeGrid( &grid );
738  return rfield.nupper;
739 }
740 
741 /* Kurucz79Compile rebin Kurucz 1979 stellar models to match energy grid of code */
743 {
744  realnum Edges[1];
745 
746  bool lgFail = false;
747 
748  DEBUG_ENTRY( "Kurucz79Compile()" );
749 
750  fprintf( ioQQQ, " Kurucz79Compile on the job.\n" );
751 
752  /* following atmospheres LTE from Kurucz 1979, Ap.J. Sup 40, 1. and
753  * Kurucz (1989) private communication, newer opacities */
754 
756 
757  if( lgFileReadable( "kurucz79.ascii", pc, as ) && !lgValidBinFile( "kurucz79.mod", pc, as ) )
758  lgFail = lgCompileAtmosphere( "kurucz79.ascii", "kurucz79.mod", Edges, 0L, pc );
759  return lgFail;
760 }
761 
762 /* Kurucz79Interpolate read in and interpolate on Kurucz79 grid of atmospheres */
763 long Kurucz79Interpolate(double val[], /* val[nval] */
764  long *nval,
765  long *ndim,
766  bool lgList,
767  double *Tlow,
768  double *Thigh)
769 {
771 
772  DEBUG_ENTRY( "Kurucz79Interpolate()" );
773 
774  grid.name = "kurucz79.mod";
775  grid.scheme = AS_DATA_OPTIONAL;
776  /* identification of this atmosphere set, used in
777  * the Cloudy output, *must* be 12 characters long */
778  grid.ident = " Kurucz 1979";
779  /* the Cloudy command needed to recompile the binary model file */
780  grid.command = "COMPILE STARS";
781 
782  InitGrid( &grid, lgList );
783 
784  CheckVal( &grid, val, nval, ndim );
785 
786  InterpolateRectGrid( &grid, val, Tlow, Thigh );
787 
788  FreeGrid( &grid );
789  return rfield.nupper;
790 }
791 
792 /* MihalasCompile rebin Mihalas stellar models to match energy grid of code */
794 {
795  realnum Edges[1];
796 
797  bool lgFail = false;
798 
799  DEBUG_ENTRY( "MihalasCompile()" );
800 
801  fprintf( ioQQQ, " MihalasCompile on the job.\n" );
802 
803  /* following atmospheres NLTE from Mihalas, NCAR-TN/STR-76 */
804 
806 
807  if( lgFileReadable( "mihalas.ascii", pc, as ) && !lgValidBinFile( "mihalas.mod", pc, as ) )
808  lgFail = lgCompileAtmosphere( "mihalas.ascii", "mihalas.mod", Edges, 0L, pc );
809  return lgFail;
810 }
811 
812 /* MihalasInterpolate read in and interpolate on Mihalas grid of atmospheres */
813 long MihalasInterpolate(double val[], /* val[nval] */
814  long *nval,
815  long *ndim,
816  bool lgList,
817  double *Tlow,
818  double *Thigh)
819 {
821 
822  DEBUG_ENTRY( "MihalasInterpolate()" );
823 
824  grid.name = "mihalas.mod";
825  grid.scheme = AS_DATA_OPTIONAL;
826  /* identification of this atmosphere set, used in
827  * the Cloudy output, *must* be 12 characters long */
828  grid.ident = " Mihalas";
829  /* the Cloudy command needed to recompile the binary model file */
830  grid.command = "COMPILE STARS";
831 
832  InitGrid( &grid, lgList );
833 
834  CheckVal( &grid, val, nval, ndim );
835 
836  InterpolateRectGrid( &grid, val, Tlow, Thigh );
837 
838  FreeGrid( &grid );
839  return rfield.nupper;
840 }
841 
842 /* RauchCompile create ascii and mod files for Rauch atmospheres */
844 {
845  bool lgFail = false;
846 
847  /* these contain frequencies for the major absorption edges */
848  realnum Edges[3];
849 
850  /* Before running this program issue the following command where the Rauch
851  * model atmosphere files are kept (0050000_50_solar_bin_0.1 and so on)
852  *
853  * ls *solar_bin_0.1 > rauchmods.list
854  *
855  * and check to see that there are 66 lines in the file.
856  */
857 
858  static const mpp telg1[NMODS_HCA] = {
859  {{50000.,5.0}}, {{50000.,6.0}}, {{50000.,7.0}}, {{50000.,8.0}},
860  {{60000.,5.0}}, {{60000.,6.0}}, {{60000.,7.0}}, {{60000.,8.0}},
861  {{70000.,5.0}}, {{70000.,6.0}}, {{70000.,7.0}}, {{70000.,8.0}},
862  {{80000.,5.0}}, {{80000.,6.0}}, {{80000.,7.0}}, {{80000.,8.0}},
863  {{90000.,5.0}}, {{90000.,6.0}}, {{90000.,7.0}}, {{90000.,8.0}},
864  {{100000.,5.0}},{{100000.,6.0}},{{100000.,7.0}},{{100000.,8.0}},
865  {{110000.,6.0}},{{110000.,7.0}},{{110000.,8.0}},
866  {{120000.,6.0}},{{120000.,7.0}},{{120000.,8.0}},
867  {{130000.,6.0}},{{130000.,7.0}},{{130000.,8.0}},
868  {{140000.,6.0}},{{140000.,7.0}},{{140000.,8.0}},
869  {{150000.,6.0}},{{150000.,7.0}},{{150000.,8.0}},
870  {{160000.,6.0}},{{160000.,7.0}},{{160000.,8.0}},
871  {{170000.,6.0}},{{170000.,7.0}},{{170000.,8.0}},
872  {{180000.,6.0}},{{180000.,7.0}},{{180000.,8.0}},
873  {{190000.,6.0}},{{190000.,7.0}},{{190000.,8.0}},
874  {{200000.,7.0}},{{200000.,8.0}},{{200000.,9.0}},
875  {{300000.,7.0}},{{300000.,8.0}},{{300000.,9.0}},
876  {{400000.,8.0}},{{400000.,9.0}},
877  {{500000.,8.0}},{{500000.,9.0}},
878  {{600000.,9.0}},
879  {{700000.,9.0}},
880  {{800000.,9.0}},
881  {{900000.,9.0}},
882  {{1000000.,9.0}}
883  };
884 
885  static const mpp telg2[NMODS_HNI] = {
886  {{50000.,5.0}}, {{50000.,6.0}}, {{50000.,7.0}}, {{50000.,8.0}},
887  {{60000.,5.0}}, {{60000.,6.0}}, {{60000.,7.0}}, {{60000.,8.0}},
888  {{70000.,5.0}}, {{70000.,6.0}}, {{70000.,7.0}}, {{70000.,8.0}},
889  {{80000.,5.0}}, {{80000.,6.0}}, {{80000.,7.0}}, {{80000.,8.0}},
890  {{90000.,5.0}}, {{90000.,6.0}}, {{90000.,7.0}}, {{90000.,8.0}},
891  {{100000.,5.0}},{{100000.,6.0}},{{100000.,7.0}},{{100000.,8.0}},
892  {{110000.,6.0}},{{110000.,7.0}},{{110000.,8.0}},
893  {{120000.,6.0}},{{120000.,7.0}},{{120000.,8.0}},
894  {{130000.,6.0}},{{130000.,7.0}},{{130000.,8.0}},
895  {{140000.,6.0}},{{140000.,7.0}},{{140000.,8.0}},
896  {{150000.,6.0}},{{150000.,7.0}},{{150000.,8.0}},
897  {{160000.,6.0}},{{160000.,7.0}},{{160000.,8.0}},
898  {{170000.,6.0}},{{170000.,7.0}},{{170000.,8.0}},
899  {{180000.,6.0}},{{180000.,7.0}},{{180000.,8.0}},
900  {{190000.,6.0}},{{190000.,7.0}},{{190000.,8.0}}
901  };
902 
903  static const mpp telg3[NMODS_PG1159] = {
904  {{40000.,5.0}}, {{40000.,6.0}}, {{40000.,7.0}}, {{40000.,8.0}}, {{40000.,9.0}},
905  {{50000.,5.0}}, {{50000.,6.0}}, {{50000.,7.0}}, {{50000.,8.0}}, {{50000.,9.0}},
906  {{60000.,5.0}}, {{60000.,6.0}}, {{60000.,7.0}}, {{60000.,8.0}}, {{60000.,9.0}},
907  {{70000.,5.0}}, {{70000.,6.0}}, {{70000.,7.0}}, {{70000.,8.0}}, {{70000.,9.0}},
908  {{80000.,5.0}}, {{80000.,6.0}}, {{80000.,7.0}}, {{80000.,8.0}}, {{80000.,9.0}},
909  {{90000.,5.0}}, {{90000.,6.0}}, {{90000.,7.0}}, {{90000.,8.0}}, {{90000.,9.0}},
910  {{100000.,5.0}},{{100000.,6.0}},{{100000.,7.0}},{{100000.,8.0}},{{100000.,9.0}},
911  {{110000.,6.0}},{{110000.,7.0}},{{110000.,8.0}},{{110000.,9.0}},
912  {{120000.,6.0}},{{120000.,7.0}},{{120000.,8.0}},{{120000.,9.0}},
913  {{130000.,6.0}},{{130000.,7.0}},{{130000.,8.0}},{{130000.,9.0}},
914  {{140000.,6.0}},{{140000.,7.0}},{{140000.,8.0}},{{140000.,9.0}},
915  {{150000.,6.0}},{{150000.,7.0}},{{150000.,8.0}},{{150000.,9.0}},
916  {{160000.,6.0}},{{160000.,7.0}},{{160000.,8.0}},{{160000.,9.0}},
917  {{170000.,6.0}},{{170000.,7.0}},{{170000.,8.0}},{{170000.,9.0}},
918  {{180000.,6.0}},{{180000.,7.0}},{{180000.,8.0}},{{180000.,9.0}},
919  {{190000.,6.0}},{{190000.,7.0}},{{190000.,8.0}},{{190000.,9.0}}
920  };
921 
922  static const mpp telg4[NMODS_HYDR] = {
923  {{20000.,4.0}}, {{20000.,5.0}}, {{20000.,6.0}}, {{20000.,7.0}}, {{20000.,8.0}}, {{20000.,9.0}},
924  {{30000.,4.0}}, {{30000.,5.0}}, {{30000.,6.0}}, {{30000.,7.0}}, {{30000.,8.0}}, {{30000.,9.0}},
925  {{40000.,4.0}}, {{40000.,5.0}}, {{40000.,6.0}}, {{40000.,7.0}}, {{40000.,8.0}}, {{40000.,9.0}},
926  {{50000.,4.0}}, {{50000.,5.0}}, {{50000.,6.0}}, {{50000.,7.0}}, {{50000.,8.0}}, {{50000.,9.0}},
927  {{60000.,5.0}}, {{60000.,6.0}}, {{60000.,7.0}}, {{60000.,8.0}}, {{60000.,9.0}},
928  {{70000.,5.0}}, {{70000.,6.0}}, {{70000.,7.0}}, {{70000.,8.0}}, {{70000.,9.0}},
929  {{80000.,5.0}}, {{80000.,6.0}}, {{80000.,7.0}}, {{80000.,8.0}}, {{80000.,9.0}},
930  {{90000.,5.0}}, {{90000.,6.0}}, {{90000.,7.0}}, {{90000.,8.0}}, {{90000.,9.0}},
931  {{100000.,5.0}},{{100000.,6.0}},{{100000.,7.0}},{{100000.,8.0}},{{100000.,9.0}},
932  {{110000.,6.0}},{{110000.,7.0}},{{110000.,8.0}},{{110000.,9.0}},
933  {{120000.,6.0}},{{120000.,7.0}},{{120000.,8.0}},{{120000.,9.0}},
934  {{130000.,6.0}},{{130000.,7.0}},{{130000.,8.0}},{{130000.,9.0}},
935  {{140000.,6.0}},{{140000.,7.0}},{{140000.,8.0}},{{140000.,9.0}},
936  {{150000.,6.0}},{{150000.,7.0}},{{150000.,8.0}},{{150000.,9.0}},
937  {{160000.,6.0}},{{160000.,7.0}},{{160000.,8.0}},{{160000.,9.0}},
938  {{170000.,6.0}},{{170000.,7.0}},{{170000.,8.0}},{{170000.,9.0}},
939  {{180000.,6.0}},{{180000.,7.0}},{{180000.,8.0}},{{180000.,9.0}},
940  {{190000.,6.0}},{{190000.,7.0}},{{190000.,8.0}},{{190000.,9.0}},
941  {{200000.,7.0}},{{200000.,8.0}},{{200000.,9.0}},
942  {{300000.,7.0}},{{300000.,8.0}},{{300000.,9.0}},
943  {{400000.,8.0}},{{400000.,9.0}},
944  {{500000.,8.0}},{{500000.,9.0}},
945  {{600000.,9.0}},
946  {{700000.,9.0}},
947  {{800000.,9.0}},
948  {{900000.,9.0}},
949  {{1000000.,9.0}}
950  };
951 
952  static const mpp telg5[NMODS_HELIUM] = {
953  {{50000.,5.0}}, {{50000.,6.0}}, {{50000.,7.0}}, {{50000.,8.0}}, {{50000.,9.0}},
954  {{60000.,5.0}}, {{60000.,6.0}}, {{60000.,7.0}}, {{60000.,8.0}}, {{60000.,9.0}},
955  {{70000.,5.0}}, {{70000.,6.0}}, {{70000.,7.0}}, {{70000.,8.0}}, {{70000.,9.0}},
956  {{80000.,5.0}}, {{80000.,6.0}}, {{80000.,7.0}}, {{80000.,8.0}}, {{80000.,9.0}},
957  {{90000.,5.0}}, {{90000.,6.0}}, {{90000.,7.0}}, {{90000.,8.0}}, {{90000.,9.0}},
958  {{100000.,5.0}},{{100000.,6.0}},{{100000.,7.0}},{{100000.,8.0}},{{100000.,9.0}},
959  {{110000.,6.0}},{{110000.,7.0}},{{110000.,8.0}},{{110000.,9.0}},
960  {{120000.,6.0}},{{120000.,7.0}},{{120000.,8.0}},{{120000.,9.0}},
961  {{130000.,6.0}},{{130000.,7.0}},{{130000.,8.0}},{{130000.,9.0}},
962  {{140000.,6.0}},{{140000.,7.0}},{{140000.,8.0}},{{140000.,9.0}},
963  {{150000.,6.0}},{{150000.,7.0}},{{150000.,8.0}},{{150000.,9.0}},
964  {{160000.,6.0}},{{160000.,7.0}},{{160000.,8.0}},{{160000.,9.0}},
965  {{170000.,6.0}},{{170000.,7.0}},{{170000.,8.0}},{{170000.,9.0}},
966  {{180000.,6.0}},{{180000.,7.0}},{{180000.,8.0}},{{180000.,9.0}},
967  {{190000.,6.0}},{{190000.,7.0}},{{190000.,8.0}},{{190000.,9.0}},
968  {{200000.,7.0}},{{200000.,8.0}},{{200000.,9.0}},
969  {{300000.,7.0}},{{300000.,8.0}},{{300000.,9.0}},
970  {{400000.,8.0}},{{400000.,9.0}},
971  {{500000.,8.0}},{{500000.,9.0}},
972  {{600000.,9.0}},
973  {{700000.,9.0}},
974  {{800000.,9.0}},
975  {{900000.,9.0}},
976  {{1000000.,9.0}}
977  };
978 
979  static const mpp telg6[NMODS_HpHE] = {
980  {{50000.,5.0}}, {{50000.,5.5}}, {{50000.,6.0}}, {{50000.,6.5}}, {{50000.,7.0}}, {{50000.,7.5}}, {{50000.,8.0}}, {{50000.,8.5}}, {{50000.,9.0}},
981  {{60000.,5.0}}, {{60000.,5.5}}, {{60000.,6.0}}, {{60000.,6.5}}, {{60000.,7.0}}, {{60000.,7.5}}, {{60000.,8.0}}, {{60000.,8.5}}, {{60000.,9.0}},
982  {{70000.,5.0}}, {{70000.,5.5}}, {{70000.,6.0}}, {{70000.,6.5}}, {{70000.,7.0}}, {{70000.,7.5}}, {{70000.,8.0}}, {{70000.,8.5}}, {{70000.,9.0}},
983  {{80000.,5.0}}, {{80000.,5.5}}, {{80000.,6.0}}, {{80000.,6.5}}, {{80000.,7.0}}, {{80000.,7.5}}, {{80000.,8.0}}, {{80000.,8.5}}, {{80000.,9.0}},
984  {{90000.,5.0}}, {{90000.,5.5}}, {{90000.,6.0}}, {{90000.,6.5}}, {{90000.,7.0}}, {{90000.,7.5}}, {{90000.,8.0}}, {{90000.,8.5}}, {{90000.,9.0}},
985  {{100000.,5.0}},{{100000.,5.5}},{{100000.,6.0}},{{100000.,6.5}},{{100000.,7.0}},{{100000.,7.5}},{{100000.,8.0}},{{100000.,8.5}},{{100000.,9.0}},
986  {{110000.,6.0}},{{110000.,6.5}},{{110000.,7.0}},{{110000.,7.5}},{{110000.,8.0}},{{110000.,8.5}},{{110000.,9.0}},
987  {{120000.,6.0}},{{120000.,6.5}},{{120000.,7.0}},{{120000.,7.5}},{{120000.,8.0}},{{120000.,8.5}},{{120000.,9.0}},
988  {{130000.,6.0}},{{130000.,6.5}},{{130000.,7.0}},{{130000.,7.5}},{{130000.,8.0}},{{130000.,8.5}},{{130000.,9.0}},
989  {{140000.,6.0}},{{140000.,6.5}},{{140000.,7.0}},{{140000.,7.5}},{{140000.,8.0}},{{140000.,8.5}},{{140000.,9.0}},
990  {{150000.,6.0}},{{150000.,6.5}},{{150000.,7.0}},{{150000.,7.5}},{{150000.,8.0}},{{150000.,8.5}},{{150000.,9.0}},
991  {{160000.,6.0}},{{160000.,6.5}},{{160000.,7.0}},{{160000.,7.5}},{{160000.,8.0}},{{160000.,8.5}},{{160000.,9.0}},
992  {{170000.,6.0}},{{170000.,6.5}},{{170000.,7.0}},{{170000.,7.5}},{{170000.,8.0}},{{170000.,8.5}},{{170000.,9.0}},
993  {{180000.,6.0}},{{180000.,6.5}},{{180000.,7.0}},{{180000.,7.5}},{{180000.,8.0}},{{180000.,8.5}},{{180000.,9.0}},
994  {{190000.,6.0}},{{190000.,6.5}},{{190000.,7.0}},{{190000.,7.5}},{{190000.,8.0}},{{190000.,8.5}},{{190000.,9.0}}
995  };
996 
997  /* metalicities of the solar and halo grid */
998  static const double par2[2] = { 0., -1. };
999 
1000  /* Helium fraction by mass */
1001  static const double par3[11] = { 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0 };
1002 
1003  DEBUG_ENTRY( "RauchCompile()" );
1004 
1005  fprintf( ioQQQ, " RauchCompile on the job.\n" );
1006 
1007  process_counter dum;
1009 
1010  /* this is the H-Ca grid */
1011  if( lgFileReadable( "0050000_50_solar_bin_0.1", dum, as ) && !lgValidAsciiFile( "rauch_h-ca_solar.ascii", as ) )
1012  {
1013  fprintf( ioQQQ, " Creating rauch_h-ca_solar.ascii....\n" );
1014  lgFail = lgFail || RauchInitializeSub( "rauch_h-ca_solar.ascii", "_solar_bin_0.1",
1015  telg1, NMODS_HCA, 1, 1, par2, 1 );
1016  }
1017 
1018  if( lgFileReadable( "0050000_50_halo__bin_0.1", dum, as ) && !lgValidAsciiFile( "rauch_h-ca_halo.ascii", as ) )
1019  {
1020  fprintf( ioQQQ, " Creating rauch_h-ca_halo.ascii....\n" );
1021  lgFail = lgFail || RauchInitializeSub( "rauch_h-ca_halo.ascii", "_halo__bin_0.1",
1022  telg1, NMODS_HCA, 1, 1, par2, 1 );
1023  }
1024 
1025  if( lgFileReadable( "0050000_50_solar_bin_0.1", dum, as ) &&
1026  lgFileReadable( "0050000_50_halo__bin_0.1", dum, as ) &&
1027  !lgValidAsciiFile( "rauch_h-ca_3d.ascii", as ) )
1028  {
1029  fprintf( ioQQQ, " Creating rauch_h-ca_3d.ascii....\n" );
1030  lgFail = lgFail || RauchInitializeSub( "rauch_h-ca_3d.ascii", "_solar_bin_0.1",
1031  telg1, NMODS_HCA, 1, 2, par2, 1 );
1032  lgFail = lgFail || RauchInitializeSub( "rauch_h-ca_3d.ascii", "_halo__bin_0.1",
1033  telg1, NMODS_HCA, 2, 2, par2, 1 );
1034  }
1035 
1036  /* this is the H-Ni grid */
1037  if( lgFileReadable( "0050000_50_solar_iron.bin_0.1", dum, as ) &&
1038  !lgValidAsciiFile( "rauch_h-ni_solar.ascii", as ) )
1039  {
1040  fprintf( ioQQQ, " Creating rauch_h-ni_solar.ascii....\n" );
1041  lgFail = lgFail || RauchInitializeSub( "rauch_h-ni_solar.ascii", "_solar_iron.bin_0.1",
1042  telg2, NMODS_HNI, 1, 1, par2, 1 );
1043  }
1044 
1045  if( lgFileReadable( "0050000_50_halo__iron.bin_0.1", dum, as ) &&
1046  !lgValidAsciiFile( "rauch_h-ni_halo.ascii", as ) )
1047  {
1048  fprintf( ioQQQ, " Creating rauch_h-ni_halo.ascii....\n" );
1049  lgFail = lgFail || RauchInitializeSub( "rauch_h-ni_halo.ascii", "_halo__iron.bin_0.1",
1050  telg2, NMODS_HNI, 1, 1, par2, 1 );
1051  }
1052 
1053  if( lgFileReadable( "0050000_50_solar_iron.bin_0.1", dum, as ) &&
1054  lgFileReadable( "0050000_50_halo__iron.bin_0.1", dum, as ) &&
1055  !lgValidAsciiFile( "rauch_h-ni_3d.ascii", as ) )
1056  {
1057  fprintf( ioQQQ, " Creating rauch_h-ni_3d.ascii....\n" );
1058  lgFail = lgFail || RauchInitializeSub( "rauch_h-ni_3d.ascii", "_solar_iron.bin_0.1",
1059  telg2, NMODS_HNI, 1, 2, par2, 1 );
1060  lgFail = lgFail || RauchInitializeSub( "rauch_h-ni_3d.ascii", "_halo__iron.bin_0.1",
1061  telg2, NMODS_HNI, 2, 2, par2, 1 );
1062  }
1063 
1064  /* this is the hydrogen deficient PG1159 grid */
1065  if( lgFileReadable( "0040000_5.00_33_50_02_15.bin_0.1", dum, as ) &&
1066  !lgValidAsciiFile( "rauch_pg1159.ascii", as ) )
1067  {
1068  fprintf( ioQQQ, " Creating rauch_pg1159.ascii....\n" );
1069  lgFail = lgFail || RauchInitializeSub( "rauch_pg1159.ascii", "_33_50_02_15.bin_0.1",
1070  telg3, NMODS_PG1159, 1, 1, par2, 2 );
1071  }
1072 
1073  /* this is the pure hydrogen grid */
1074  if( lgFileReadable( "0020000_4.00_H_00005-02000A.bin_0.1", dum, as ) &&
1075  !lgValidAsciiFile( "rauch_hydr.ascii", as ) )
1076  {
1077  fprintf( ioQQQ, " Creating rauch_hydr.ascii....\n" );
1078  lgFail = lgFail || RauchInitializeSub( "rauch_hydr.ascii", "_H_00005-02000A.bin_0.1",
1079  telg4, NMODS_HYDR, 1, 1, par2, 2 );
1080  }
1081 
1082  /* this is the pure helium grid */
1083  if( lgFileReadable( "0050000_5.00_He_00005-02000A.bin_0.1", dum, as ) &&
1084  !lgValidAsciiFile( "rauch_helium.ascii", as ) )
1085  {
1086  fprintf( ioQQQ, " Creating rauch_helium.ascii....\n" );
1087  lgFail = lgFail || RauchInitializeSub( "rauch_helium.ascii", "_He_00005-02000A.bin_0.1",
1088  telg5, NMODS_HELIUM, 1, 1, par2, 2 );
1089  }
1090 
1091  /* this is the 3D grid for arbitrary H+He mixtures */
1092  if( lgFileReadable( "0050000_5.00_H+He_1.000_0.000_00005-02000A.bin_0.1", dum, as ) &&
1093  !lgValidAsciiFile( "rauch_h+he_3d.ascii", as ) )
1094  {
1095  fprintf( ioQQQ, " Creating rauch_h+he_3d.ascii....\n" );
1096  lgFail = lgFail || RauchInitializeSub( "rauch_h+he_3d.ascii", "_H+He_1.000_0.000_00005-02000A.bin_0.1",
1097  telg6, NMODS_HpHE, 1, 11, par3, 2 );
1098  lgFail = lgFail || RauchInitializeSub( "rauch_h+he_3d.ascii", "_H+He_0.900_0.100_00005-02000A.bin_0.1",
1099  telg6, NMODS_HpHE, 2, 11, par3, 2 );
1100  lgFail = lgFail || RauchInitializeSub( "rauch_h+he_3d.ascii", "_H+He_0.800_0.200_00005-02000A.bin_0.1",
1101  telg6, NMODS_HpHE, 3, 11, par3, 2 );
1102  lgFail = lgFail || RauchInitializeSub( "rauch_h+he_3d.ascii", "_H+He_0.700_0.300_00005-02000A.bin_0.1",
1103  telg6, NMODS_HpHE, 4, 11, par3, 2 );
1104  lgFail = lgFail || RauchInitializeSub( "rauch_h+he_3d.ascii", "_H+He_0.600_0.400_00005-02000A.bin_0.1",
1105  telg6, NMODS_HpHE, 5, 11, par3, 2 );
1106  lgFail = lgFail || RauchInitializeSub( "rauch_h+he_3d.ascii", "_H+He_0.500_0.500_00005-02000A.bin_0.1",
1107  telg6, NMODS_HpHE, 6, 11, par3, 2 );
1108  lgFail = lgFail || RauchInitializeSub( "rauch_h+he_3d.ascii", "_H+He_0.400_0.600_00005-02000A.bin_0.1",
1109  telg6, NMODS_HpHE, 7, 11, par3, 2 );
1110  lgFail = lgFail || RauchInitializeSub( "rauch_h+he_3d.ascii", "_H+He_0.300_0.700_00005-02000A.bin_0.1",
1111  telg6, NMODS_HpHE, 8, 11, par3, 2 );
1112  lgFail = lgFail || RauchInitializeSub( "rauch_h+he_3d.ascii", "_H+He_0.200_0.800_00005-02000A.bin_0.1",
1113  telg6, NMODS_HpHE, 9, 11, par3, 2 );
1114  lgFail = lgFail || RauchInitializeSub( "rauch_h+he_3d.ascii", "_H+He_0.100_0.900_00005-02000A.bin_0.1",
1115  telg6, NMODS_HpHE, 10, 11, par3, 2 );
1116  lgFail = lgFail || RauchInitializeSub( "rauch_h+he_3d.ascii", "_H+He_0.000_1.000_00005-02000A.bin_0.1",
1117  telg6, NMODS_HpHE, 11, 11, par3, 2 );
1118  }
1119 
1120  /* define the major absorption edges that require special attention during rebinning
1121  *
1122  * NB the frequencies should be chosen here such that they are somewhere in between
1123  * the two frequency points that straddle the edge in the atmosphere model, the
1124  * software in RebinAtmosphere will seek out the exact values of those two points
1125  * e.g.: in the CoStar models the H I edge is straddled by wavelength points at
1126  * 911.67 and 911.85 A, so Edges[0] should be chosen somewhere in between (e.g. at 911.76A).
1127  *
1128  * NB beware not to choose edges too close to one another (i.e. on the order of the
1129  * resolution of the Cloudy frequency grid). E.g. the He II Balmer edge nearly coincides
1130  * with the H I Ly edge, they should be treated as one edge. Trying to separate them will
1131  * almost certainly lead to erroneous behaviour in RebinAtmosphere */
1132  Edges[0] = 0.99946789f;
1133  Edges[1] = 1.8071406f;
1134  Edges[2] = 3.9996377f;
1135 
1136  if( lgFileReadable( "rauch_h-ca_solar.ascii", pc, as ) && !lgValidBinFile( "rauch_h-ca_solar.mod", pc, as ) )
1137  lgFail = lgFail || lgCompileAtmosphere( "rauch_h-ca_solar.ascii", "rauch_h-ca_solar.mod",Edges,3L, pc );
1138  if( lgFileReadable( "rauch_h-ca_halo.ascii", pc, as ) && !lgValidBinFile( "rauch_h-ca_halo.mod", pc, as ) )
1139  lgFail = lgFail || lgCompileAtmosphere( "rauch_h-ca_halo.ascii", "rauch_h-ca_halo.mod", Edges, 3L, pc );
1140  if( lgFileReadable( "rauch_h-ca_3d.ascii", pc, as ) && !lgValidBinFile( "rauch_h-ca_3d.mod", pc, as ) )
1141  lgFail = lgFail || lgCompileAtmosphere( "rauch_h-ca_3d.ascii", "rauch_h-ca_3d.mod", Edges, 3L, pc );
1142 
1143  if( lgFileReadable( "rauch_h-ni_solar.ascii", pc, as ) && !lgValidBinFile( "rauch_h-ni_solar.mod", pc, as ) )
1144  lgFail = lgFail || lgCompileAtmosphere( "rauch_h-ni_solar.ascii", "rauch_h-ni_solar.mod",Edges,3L, pc );
1145  if( lgFileReadable( "rauch_h-ni_halo.ascii", pc, as ) && !lgValidBinFile( "rauch_h-ni_halo.mod", pc, as ) )
1146  lgFail = lgFail || lgCompileAtmosphere( "rauch_h-ni_halo.ascii", "rauch_h-ni_halo.mod", Edges, 3L, pc );
1147  if( lgFileReadable( "rauch_h-ni_3d.ascii", pc, as ) && !lgValidBinFile( "rauch_h-ni_3d.mod", pc, as ) )
1148  lgFail = lgFail || lgCompileAtmosphere( "rauch_h-ni_3d.ascii", "rauch_h-ni_3d.mod", Edges, 3L, pc );
1149 
1150  if( lgFileReadable( "rauch_pg1159.ascii", pc, as ) && !lgValidBinFile( "rauch_pg1159.mod", pc, as ) )
1151  lgFail = lgFail || lgCompileAtmosphere( "rauch_pg1159.ascii", "rauch_pg1159.mod", Edges, 3L, pc );
1152 
1153  if( lgFileReadable( "rauch_hydr.ascii", pc, as ) && !lgValidBinFile( "rauch_hydr.mod", pc, as ) )
1154  lgFail = lgFail || lgCompileAtmosphere( "rauch_hydr.ascii", "rauch_hydr.mod", Edges, 3L, pc );
1155 
1156  if( lgFileReadable( "rauch_helium.ascii", pc, as ) && !lgValidBinFile( "rauch_helium.mod", pc, as ) )
1157  lgFail = lgFail || lgCompileAtmosphere( "rauch_helium.ascii", "rauch_helium.mod", Edges, 3L, pc );
1158 
1159  if( lgFileReadable( "rauch_h+he_3d.ascii", pc, as ) && !lgValidBinFile( "rauch_h+he_3d.mod", pc, as ) )
1160  lgFail = lgFail || lgCompileAtmosphere( "rauch_h+he_3d.ascii", "rauch_h+he_3d.mod", Edges, 3L, pc );
1161  return lgFail;
1162 }
1163 
1164 /* RauchInterpolateHCa get one of the Rauch H-Ca model atmospheres, originally by K. Volk */
1165 long RauchInterpolateHCa(double val[], /* val[nval] */
1166  long *nval,
1167  long *ndim,
1168  bool lgHalo,
1169  bool lgList,
1170  double *Tlow,
1171  double *Thigh)
1172 {
1174 
1175  DEBUG_ENTRY( "RauchInterpolateHCa()" );
1176 
1177  if( *ndim == 3 )
1178  grid.name = "rauch_h-ca_3d.mod";
1179  else
1180  grid.name = ( lgHalo ? "rauch_h-ca_halo.mod" : "rauch_h-ca_solar.mod" );
1181  grid.scheme = AS_DATA_OPTIONAL;
1182  /* identification of this atmosphere set, used in
1183  * the Cloudy output, *must* be 12 characters long */
1184  grid.ident = " H-Ca Rauch";
1185  /* the Cloudy command needed to recompile the binary model file */
1186  grid.command = "COMPILE STARS";
1187 
1188  InitGrid( &grid, lgList );
1189 
1190  CheckVal( &grid, val, nval, ndim );
1191 
1192  InterpolateRectGrid( &grid, val, Tlow, Thigh );
1193 
1194  FreeGrid( &grid );
1195  return rfield.nupper;
1196 }
1197 
1198 /* RauchInterpolateHNi get one of the Rauch H-Ni model atmospheres */
1199 long RauchInterpolateHNi(double val[], /* val[nval] */
1200  long *nval,
1201  long *ndim,
1202  bool lgHalo,
1203  bool lgList,
1204  double *Tlow,
1205  double *Thigh)
1206 {
1208 
1209  DEBUG_ENTRY( "RauchInterpolateHNi()" );
1210 
1211  if( *ndim == 3 )
1212  grid.name = "rauch_h-ni_3d.mod";
1213  else
1214  grid.name = ( lgHalo ? "rauch_h-ni_halo.mod" : "rauch_h-ni_solar.mod" );
1215  grid.scheme = AS_DATA_OPTIONAL;
1216  /* identification of this atmosphere set, used in
1217  * the Cloudy output, *must* be 12 characters long */
1218  grid.ident = " H-Ni Rauch";
1219  /* the Cloudy command needed to recompile the binary model file */
1220  grid.command = "COMPILE STARS";
1221 
1222  InitGrid( &grid, lgList );
1223 
1224  CheckVal( &grid, val, nval, ndim );
1225 
1226  InterpolateRectGrid( &grid, val, Tlow, Thigh );
1227 
1228  FreeGrid( &grid );
1229  return rfield.nupper;
1230 }
1231 
1232 /* RauchInterpolatePG1159 get one of the Rauch PG1159 model atmospheres */
1233 long RauchInterpolatePG1159(double val[], /* val[nval] */
1234  long *nval,
1235  long *ndim,
1236  bool lgList,
1237  double *Tlow,
1238  double *Thigh)
1239 {
1241 
1242  DEBUG_ENTRY( "RauchInterpolatePG1159()" );
1243 
1244  grid.name = "rauch_pg1159.mod";
1245  grid.scheme = AS_DATA_OPTIONAL;
1246  /* identification of this atmosphere set, used in
1247  * the Cloudy output, *must* be 12 characters long */
1248  grid.ident = "PG1159 Rauch";
1249  /* the Cloudy command needed to recompile the binary model file */
1250  grid.command = "COMPILE STARS";
1251 
1252  InitGrid( &grid, lgList );
1253 
1254  CheckVal( &grid, val, nval, ndim );
1255 
1256  InterpolateRectGrid( &grid, val, Tlow, Thigh );
1257 
1258  FreeGrid( &grid );
1259  return rfield.nupper;
1260 }
1261 
1262 /* RauchInterpolateHydr get one of the Rauch pure hydrogen model atmospheres */
1263 long RauchInterpolateHydr(double val[], /* val[nval] */
1264  long *nval,
1265  long *ndim,
1266  bool lgList,
1267  double *Tlow,
1268  double *Thigh)
1269 {
1271 
1272  DEBUG_ENTRY( "RauchInterpolateHydr()" );
1273 
1274  grid.name = "rauch_hydr.mod";
1275  grid.scheme = AS_DATA_OPTIONAL;
1276  /* identification of this atmosphere set, used in
1277  * the Cloudy output, *must* be 12 characters long */
1278  grid.ident = " Hydr Rauch";
1279  /* the Cloudy command needed to recompile the binary model file */
1280  grid.command = "COMPILE STARS";
1281 
1282  InitGrid( &grid, lgList );
1283 
1284  CheckVal( &grid, val, nval, ndim );
1285 
1286  InterpolateRectGrid( &grid, val, Tlow, Thigh );
1287 
1288  FreeGrid( &grid );
1289  return rfield.nupper;
1290 }
1291 
1292 /* RauchInterpolateHelium get one of the Rauch pure helium model atmospheres */
1293 long RauchInterpolateHelium(double val[], /* val[nval] */
1294  long *nval,
1295  long *ndim,
1296  bool lgList,
1297  double *Tlow,
1298  double *Thigh)
1299 {
1301 
1302  DEBUG_ENTRY( "RauchInterpolateHelium()" );
1303 
1304  grid.name = "rauch_helium.mod";
1305  grid.scheme = AS_DATA_OPTIONAL;
1306  /* identification of this atmosphere set, used in
1307  * the Cloudy output, *must* be 12 characters long */
1308  grid.ident = "Helium Rauch";
1309  /* the Cloudy command needed to recompile the binary model file */
1310  grid.command = "COMPILE STARS";
1311 
1312  InitGrid( &grid, lgList );
1313 
1314  CheckVal( &grid, val, nval, ndim );
1315 
1316  InterpolateRectGrid( &grid, val, Tlow, Thigh );
1317 
1318  FreeGrid( &grid );
1319  return rfield.nupper;
1320 }
1321 
1322 /* RauchInterpolateHpHe get one of the Rauch hydrogen plus helium model atmospheres */
1323 long RauchInterpolateHpHe(double val[], /* val[nval] */
1324  long *nval,
1325  long *ndim,
1326  bool lgList,
1327  double *Tlow,
1328  double *Thigh)
1329 {
1331 
1332  DEBUG_ENTRY( "RauchInterpolateHpHe()" );
1333 
1334  grid.name = "rauch_h+he_3d.mod";
1335  grid.scheme = AS_DATA_OPTIONAL;
1336  /* identification of this atmosphere set, used in
1337  * the Cloudy output, *must* be 12 characters long */
1338  grid.ident = " H+He Rauch";
1339  /* the Cloudy command needed to recompile the binary model file */
1340  grid.command = "COMPILE STARS";
1341 
1342  InitGrid( &grid, lgList );
1343 
1344  CheckVal( &grid, val, nval, ndim );
1345 
1346  InterpolateRectGrid( &grid, val, Tlow, Thigh );
1347 
1348  FreeGrid( &grid );
1349  return rfield.nupper;
1350 }
1351 
1352 /* StarburstInitialize does the actual work of preparing the ascii file */
1353 bool StarburstInitialize(const char chInName[],
1354  const char chOutName[])
1355 {
1356  char chLine[INPUT_LINE_LENGTH]; /* used for getting input lines */
1357 
1358  bool lgHeader = true;
1359  long int i, j, nmods, ngp;
1360 
1361  size_t nsb_sz = (size_t)NSB99;
1362 
1363  double *wavl, *fluxes[MNTS], Age[MNTS], lwavl;
1364 
1365  FILE *ioOut, /* pointer to output file we came here to create*/
1366  *ioIn; /* pointer to input files we will read */
1367 
1368  DEBUG_ENTRY( "StarburstInitialize()" );
1369 
1370  for( i=0; i < MNTS; i++ )
1371  fluxes[i] = NULL;
1372 
1373  /* grab some space for the wavelengths and fluxes */
1374  wavl = (double *)MALLOC( sizeof(double)*nsb_sz);
1375 
1376  ioIn = open_data( chInName, "r", AS_LOCAL_ONLY );
1377 
1378  lwavl = 0.;
1379  nmods = 0;
1380  ngp = 0;
1381 
1382  while( read_whole_line( chLine, INPUT_LINE_LENGTH, ioIn ) != NULL )
1383  {
1384  if( !lgHeader )
1385  {
1386  double cage, cwavl, cfl1;
1387 
1388  /* format: age/yr wavl/Angstrom log10(flux_total) log10(flux_stellar) log10(flux_neb) */
1389  /* we are only interested in the total flux, so we ignore the remaining numbers */
1390  if( sscanf( chLine, " %le %le %le", &cage, &cwavl, &cfl1 ) != 3 )
1391  {
1392  fprintf( ioQQQ, "syntax error in data of Starburst grid.\n" );
1393  goto error;
1394  }
1395 
1396  if( cwavl < lwavl )
1397  {
1398  ++nmods;
1399  ngp = 0;
1400 
1401  if( nmods >= MNTS )
1402  {
1403  fprintf( ioQQQ, "too many time steps in Starburst grid.\n" );
1404  fprintf( ioQQQ, "please increase MNTS and recompile.\n" );
1405  goto error;
1406  }
1407  }
1408 
1409  if( ngp == 0 )
1410  {
1411  fluxes[nmods] = (double *)MALLOC( sizeof(double)*nsb_sz);
1412  Age[nmods] = cage;
1413  }
1414 
1415  if( ngp >= (long)nsb_sz )
1416  {
1417  /* this should only be needed when nmods == 0 */
1418  ASSERT( nmods == 0 );
1419 
1420  nsb_sz *= 2;
1421  fluxes[0] = (double *)REALLOC(fluxes[0],sizeof(double)*nsb_sz);
1422  wavl = (double *)REALLOC(wavl,sizeof(double)*nsb_sz);
1423  }
1424 
1425  if( fabs(Age[nmods]-cage) > 10.*DBL_EPSILON*Age[nmods] )
1426  {
1427  fprintf( ioQQQ, "age error in Starburst grid.\n" );
1428  goto error;
1429  }
1430 
1431  if( nmods == 0 )
1432  wavl[ngp] = cwavl;
1433  else
1434  {
1435  if( fabs(wavl[ngp]-cwavl) > 10.*DBL_EPSILON*wavl[ngp] )
1436  {
1437  fprintf( ioQQQ, "wavelength error in Starburst grid.\n" );
1438  goto error;
1439  }
1440  }
1441 
1442  /* arbitrarily renormalize to flux in erg/cm^2/s/A at 1kpc */
1443  /* constant is log10( 4*pi*(kpc/cm)^2 ) */
1444  fluxes[nmods][ngp] = pow( 10., cfl1 - 44.077911 );
1445 
1446  lwavl = cwavl;
1447  ++ngp;
1448  }
1449 
1450  if( lgHeader && strncmp( &chLine[1], "TIME [YR]", 9 ) == 0 )
1451  lgHeader = false;
1452  }
1453 
1454  if( lgHeader )
1455  {
1456  /* this happens when the "TIME [YR]" string was not found in column 1 of the file */
1457  fprintf( ioQQQ, "syntax error in header of Starburst grid.\n" );
1458  goto error;
1459  }
1460 
1461  ++nmods;
1462 
1463  /* finished - close the unit */
1464  fclose(ioIn);
1465 
1466  ioOut = open_data( chOutName, "w", AS_LOCAL_ONLY );
1467 
1468  fprintf( ioOut, " %ld\n", VERSION_ASCII );
1469  fprintf( ioOut, " %d\n", 1 );
1470  fprintf( ioOut, " %d\n", 1 );
1471  fprintf( ioOut, " Age\n" );
1472  fprintf( ioOut, " %ld\n", nmods );
1473  fprintf( ioOut, " %ld\n", ngp );
1474  /* Starburst99 models give the wavelength in Angstrom */
1475  fprintf( ioOut, " lambda\n" );
1476  /* conversion factor to Angstrom */
1477  fprintf( ioOut, " %.8e\n", 1. );
1478  /* Starburst99 models give the total flux F_lambda in erg/s/A, will be renormalized at 1 kpc */
1479  fprintf( ioOut, " F_lambda\n" );
1480  /* this factor is irrelevant since Teff check will not be carried out */
1481  fprintf( ioOut, " %.8e\n", 1. );
1482  /* write out the Ages */
1483  for( i=0; i < nmods; i++ )
1484  {
1485  fprintf( ioOut, " %.3e", Age[i] );
1486  if( ((i+1)%4) == 0 )
1487  fprintf( ioOut, "\n" );
1488  }
1489  if( (i%4) != 0 )
1490  fprintf( ioOut, "\n" );
1491 
1492  fprintf( ioQQQ, " Writing: " );
1493 
1494  /* write out the wavelength grid */
1495  for( j=0; j < ngp; j++ )
1496  {
1497  fprintf( ioOut, " %.4e", wavl[j] );
1498  /* want to have 5 numbers per line */
1499  if( ((j+1)%5) == 0 )
1500  fprintf( ioOut, "\n" );
1501  }
1502  /* need to throw a newline if we did not end on an exact line */
1503  if( (j%5) != 0 )
1504  fprintf( ioOut, "\n" );
1505 
1506  /* print to screen to show progress */
1507  fprintf( ioQQQ, "." );
1508  fflush( ioQQQ );
1509 
1510  for( i=0; i < nmods; i++ )
1511  {
1512  for( j=0; j < ngp; j++ )
1513  {
1514  fprintf( ioOut, " %.4e", fluxes[i][j] );
1515  /* want to have 5 numbers per line */
1516  if( ((j+1)%5) == 0 )
1517  fprintf( ioOut, "\n" );
1518  }
1519  /* need to throw a newline if we did not end on an exact line */
1520  if( (j%5) != 0 )
1521  fprintf( ioOut, "\n" );
1522 
1523  /* print to screen to show progress */
1524  fprintf( ioQQQ, "." );
1525  fflush( ioQQQ );
1526  }
1527 
1528  fprintf( ioQQQ, " done.\n" );
1529 
1530  fclose(ioOut);
1531 
1532  /* free the space grabbed above */
1533  for( i=0; i < MNTS; i++ )
1534  FREE_SAFE( fluxes[i] );
1535  FREE_CHECK( wavl );
1536  return false;
1537 
1538 error:
1539  for( i=0; i < MNTS; i++ )
1540  FREE_SAFE( fluxes[i] );
1541  FREE_CHECK( wavl );
1542  return true;
1543 }
1544 
1545 /* StarburstCompile, rebin Starburst99 model output to match energy grid of code */
1547 {
1548  realnum Edges[1];
1549 
1550  bool lgFail = false;
1551 
1552  DEBUG_ENTRY( "StarburstCompile()" );
1553 
1554  fprintf( ioQQQ, " StarburstCompile on the job.\n" );
1555 
1556  process_counter dum;
1558 
1559  if( lgFileReadable( "starburst99.stb99", dum, as ) && !lgValidAsciiFile( "starburst99.ascii", as ) )
1560  lgFail = lgFail || StarburstInitialize( "starburst99.stb99", "starburst99.ascii" );
1561  if( lgFileReadable( "starburst99.ascii", pc, as ) && !lgValidBinFile( "starburst99.mod", pc, as ) )
1562  lgFail = lgFail || lgCompileAtmosphere( "starburst99.ascii", "starburst99.mod", Edges, 0L, pc );
1563  return lgFail;
1564 }
1565 
1566 /* TlustyCompile rebin Tlusty BSTAR2006/OSTAR2002 stellar models to match energy grid of code */
1568 {
1569  /* these contain frequencies for the major absorption edges */
1570  realnum Edges[1];
1571 
1572  bool lgFail = false;
1573 
1574  DEBUG_ENTRY( "TlustyCompile()" );
1575 
1576  fprintf( ioQQQ, " TlustyCompile on the job.\n" );
1577 
1579 
1580  if( lgFileReadable( "bstar2006_p03.ascii", pc, as ) && !lgValidBinFile( "bstar2006_p03.mod", pc, as ) )
1581  lgFail = lgFail || lgCompileAtmosphere( "bstar2006_p03.ascii", "bstar2006_p03.mod", Edges, 0L, pc );
1582  if( lgFileReadable( "bstar2006_p00.ascii", pc, as ) && !lgValidBinFile( "bstar2006_p00.mod", pc, as ) )
1583  lgFail = lgFail || lgCompileAtmosphere( "bstar2006_p00.ascii", "bstar2006_p00.mod", Edges, 0L, pc );
1584  if( lgFileReadable( "bstar2006_m03.ascii", pc, as ) && !lgValidBinFile( "bstar2006_m03.mod", pc, as ) )
1585  lgFail = lgFail || lgCompileAtmosphere( "bstar2006_m03.ascii", "bstar2006_m03.mod", Edges, 0L, pc );
1586  if( lgFileReadable( "bstar2006_m07.ascii", pc, as ) && !lgValidBinFile( "bstar2006_m07.mod", pc, as ) )
1587  lgFail = lgFail || lgCompileAtmosphere( "bstar2006_m07.ascii", "bstar2006_m07.mod", Edges, 0L, pc );
1588  if( lgFileReadable( "bstar2006_m10.ascii", pc, as ) && !lgValidBinFile( "bstar2006_m10.mod", pc, as ) )
1589  lgFail = lgFail || lgCompileAtmosphere( "bstar2006_m10.ascii", "bstar2006_m10.mod", Edges, 0L, pc );
1590  if( lgFileReadable( "bstar2006_m99.ascii", pc, as ) && !lgValidBinFile( "bstar2006_m99.mod", pc, as ) )
1591  lgFail = lgFail || lgCompileAtmosphere( "bstar2006_m99.ascii", "bstar2006_m99.mod", Edges, 0L, pc );
1592 
1593  if( lgFileReadable( "bstar2006_3d.ascii", pc, as ) && !lgValidBinFile( "bstar2006_3d.mod", pc, as ) )
1594  lgFail = lgFail || lgCompileAtmosphere( "bstar2006_3d.ascii", "bstar2006_3d.mod", Edges, 0L, pc );
1595 
1596  if( lgFileReadable( "ostar2002_p03.ascii", pc, as ) && !lgValidBinFile( "ostar2002_p03.mod", pc, as ) )
1597  lgFail = lgFail || lgCompileAtmosphere( "ostar2002_p03.ascii", "ostar2002_p03.mod", Edges, 0L, pc );
1598  if( lgFileReadable( "ostar2002_p00.ascii", pc, as ) && !lgValidBinFile( "ostar2002_p00.mod", pc, as ) )
1599  lgFail = lgFail || lgCompileAtmosphere( "ostar2002_p00.ascii", "ostar2002_p00.mod", Edges, 0L, pc );
1600  if( lgFileReadable( "ostar2002_m03.ascii", pc, as ) && !lgValidBinFile( "ostar2002_m03.mod", pc, as ) )
1601  lgFail = lgFail || lgCompileAtmosphere( "ostar2002_m03.ascii", "ostar2002_m03.mod", Edges, 0L, pc );
1602  if( lgFileReadable( "ostar2002_m07.ascii", pc, as ) && !lgValidBinFile( "ostar2002_m07.mod", pc, as ) )
1603  lgFail = lgFail || lgCompileAtmosphere( "ostar2002_m07.ascii", "ostar2002_m07.mod", Edges, 0L, pc );
1604  if( lgFileReadable( "ostar2002_m10.ascii", pc, as ) && !lgValidBinFile( "ostar2002_m10.mod", pc, as ) )
1605  lgFail = lgFail || lgCompileAtmosphere( "ostar2002_m10.ascii", "ostar2002_m10.mod", Edges, 0L, pc );
1606  if( lgFileReadable( "ostar2002_m15.ascii", pc, as ) && !lgValidBinFile( "ostar2002_m15.mod", pc, as ) )
1607  lgFail = lgFail || lgCompileAtmosphere( "ostar2002_m15.ascii", "ostar2002_m15.mod", Edges, 0L, pc );
1608  if( lgFileReadable( "ostar2002_m17.ascii", pc, as ) && !lgValidBinFile( "ostar2002_m17.mod", pc, as ) )
1609  lgFail = lgFail || lgCompileAtmosphere( "ostar2002_m17.ascii", "ostar2002_m17.mod", Edges, 0L, pc );
1610  if( lgFileReadable( "ostar2002_m20.ascii", pc, as ) && !lgValidBinFile( "ostar2002_m20.mod", pc, as ) )
1611  lgFail = lgFail || lgCompileAtmosphere( "ostar2002_m20.ascii", "ostar2002_m20.mod", Edges, 0L, pc );
1612  if( lgFileReadable( "ostar2002_m30.ascii", pc, as ) && !lgValidBinFile( "ostar2002_m30.mod", pc, as ) )
1613  lgFail = lgFail || lgCompileAtmosphere( "ostar2002_m30.ascii", "ostar2002_m30.mod", Edges, 0L, pc );
1614  if( lgFileReadable( "ostar2002_m99.ascii", pc, as ) && !lgValidBinFile( "ostar2002_m99.mod", pc, as ) )
1615  lgFail = lgFail || lgCompileAtmosphere( "ostar2002_m99.ascii", "ostar2002_m99.mod", Edges, 0L, pc );
1616 
1617  if( lgFileReadable( "ostar2002_3d.ascii", pc, as ) && !lgValidBinFile( "ostar2002_3d.mod", pc, as ) )
1618  lgFail = lgFail || lgCompileAtmosphere( "ostar2002_3d.ascii", "ostar2002_3d.mod", Edges, 0L, pc );
1619  return lgFail;
1620 }
1621 
1622 /* TlustyInterpolate get one of the Tlusty BSTAR2006/OSTAR2002 model atmospheres */
1623 long TlustyInterpolate(double val[], /* val[nval] */
1624  long *nval,
1625  long *ndim,
1626  tl_grid tlg,
1627  const char *chMetalicity,
1628  bool lgList,
1629  double *Tlow,
1630  double *Thigh)
1631 {
1632  char chIdent[13];
1634 
1635  DEBUG_ENTRY( "TlustyInterpolate()" );
1636 
1637  if( tlg == TL_BSTAR )
1638  grid.name = "bstar2006_";
1639  else if( tlg == TL_OSTAR )
1640  grid.name = "ostar2002_";
1641  else
1642  TotalInsanity();
1643  if( *ndim == 3 )
1644  grid.name += "3d";
1645  else
1646  grid.name += chMetalicity;
1647  grid.name += ".mod";
1648  grid.scheme = AS_DATA_OPTIONAL;
1649  /* identification of this atmosphere set, used in
1650  * the Cloudy output, *must* be 12 characters long */
1651  if( *ndim == 3 )
1652  {
1653  strcpy( chIdent, "3-dim" );
1654  }
1655  else
1656  {
1657  strcpy( chIdent, "Z " );
1658  strcat( chIdent, chMetalicity );
1659  }
1660  if( tlg == TL_BSTAR )
1661  strcat( chIdent, " Bstr06" );
1662  else if( tlg == TL_OSTAR )
1663  strcat( chIdent, " Ostr02" );
1664  else
1665  TotalInsanity();
1666  grid.ident = chIdent;
1667  /* the Cloudy command needed to recompile the binary model file */
1668  grid.command = "COMPILE STARS";
1669 
1670  InitGrid( &grid, lgList );
1671 
1672  CheckVal( &grid, val, nval, ndim );
1673 
1674  InterpolateRectGrid( &grid, val, Tlow, Thigh );
1675 
1676  FreeGrid( &grid );
1677  return rfield.nupper;
1678 }
1679 
1680 /* WernerCompile rebin Werner stellar models to match energy grid of code */
1681 /* >>chng 05 nov 16, added return value to indicate success (0) or failure (1) */
1683 {
1684  /* these contain frequencies for the major absorption edges */
1685  realnum Edges[3];
1686 
1687  bool lgFail = false;
1688 
1689  DEBUG_ENTRY( "WernerCompile()" );
1690 
1691  fprintf( ioQQQ, " WernerCompile on the job.\n" );
1692 
1693  /* define the major absorption edges that require special attention during rebinning
1694  *
1695  * NB the frequencies should be chosen here such that they are somewhere in between
1696  * the two frequency points that straddle the edge in the atmosphere model, the
1697  * software in RebinAtmosphere will seek out the exact values of those two points
1698  * e.g.: in the CoStar models the H I edge is straddled by wavelength points at
1699  * 911.67 and 911.85 A, so Edges[0] should be chosen somewhere in between (e.g. at 911.76A).
1700  *
1701  * NB beware not to choose edges too close to one another (i.e. on the order of the
1702  * resolution of the Cloudy frequency grid). E.g. the He II Balmer edge nearly coincides
1703  * with the H I Ly edge, they should be treated as one edge. Trying to separate them will
1704  * almost certainly lead to erroneous behaviour in RebinAtmosphere */
1705  Edges[0] = 0.99946789f;
1706  Edges[1] = 1.8071406f;
1707  Edges[2] = 3.9996377f;
1708 
1709  /* The "kwerner.ascii" file is a modified ascii dump of the Klaus Werner
1710  * stellar model files which he gave to me in 1992. The first set of values
1711  * is the frequency grid (in Ryd) followed by the atmosphere models in order
1712  * of increasing temperature and log(g). The following comments are already
1713  * incorporated in the modified kwerner.ascii file that is supplied with Cloudy.
1714  *
1715  * >>chng 00 oct 18, The frequency grid was slightly tweaked compared to the
1716  * original values supplied by Klaus Werner to make it monotonically increasing;
1717  * this is due to there being fluxes above and below certain wavelengths where
1718  * the opacity changes (i.e. the Lyman and Balmer limits for example) which are
1719  * assigned the same wavelength in the original Klaus Werner files. PvH
1720  *
1721  * >>chng 00 oct 20, StarEner[172] is out of sequence. As per the Klaus Werner comment,
1722  * it should be omitted. The energy grid is very dense in this region and was most likely
1723  * intended to sample an absorption line which was not included in this particular grid.
1724  * StarFlux[172] is therefore always equal to the flux in neighbouring points (at least
1725  * those with slightly smaller energies). It is therefore safe to ignore this point. PvH
1726  *
1727  * >>chng 00 oct 20, As per the same comment, StarFlux[172] is also deleted. PvH */
1728 
1730 
1731  if( lgFileReadable( "kwerner.ascii", pc, as ) && !lgValidBinFile( "kwerner.mod", pc, as ) )
1732  lgFail = lgCompileAtmosphere( "kwerner.ascii", "kwerner.mod", Edges, 3L, pc );
1733  return lgFail;
1734 }
1735 
1736 /* WernerInterpolate read in and interpolate on Werner grid of PN atmospheres, originally by K Volk */
1737 long WernerInterpolate(double val[], /* val[nval] */
1738  long *nval,
1739  long *ndim,
1740  bool lgList,
1741  double *Tlow,
1742  double *Thigh)
1743 {
1745 
1746  DEBUG_ENTRY( "WernerInterpolate()" );
1747 
1748  /* This subroutine was added (28 dec 1992) to read from the set of
1749  * hot white dwarf model atmospheres from Klaus Werner at Kiel. The
1750  * values are read in (energy in Rydberg units, f_nu in cgs units)
1751  * for any of the 20 models. Each model had 513 points before rebinning.
1752  * The Rayleigh-Jeans tail was extrapolated. */
1753 
1754  grid.name = "kwerner.mod";
1755  grid.scheme = AS_DATA_OPTIONAL;
1756  /* identification of this atmosphere set, used in
1757  * the Cloudy output, *must* be 12 characters long */
1758  grid.ident = "Klaus Werner";
1759  /* the Cloudy command needed to recompile the binary model file */
1760  grid.command = "COMPILE STARS";
1761 
1762  InitGrid( &grid, lgList );
1763 
1764  CheckVal( &grid, val, nval, ndim );
1765 
1766  /* Note on the interpolation: 26 October 2000 (Peter van Hoof)
1767  *
1768  * I computed the effective temperature for a random sample of interpolated
1769  * atmospheres by integrating the flux as shown above and compared the results
1770  * with the expected effective temperature using DELTA = (COMP-EXPEC)/EXPEC.
1771  *
1772  * I found that the average discrepancy was:
1773  *
1774  * DELTA = -0.71% +/- 0.71% (sample size 5000)
1775  *
1776  * The most extreme discrepancies were
1777  * -4.37% <= DELTA <= 0.24%
1778  *
1779  * The most negative discrepancies were for Teff = 95 kK, log(g) = 5
1780  * The most positive discrepancies were for Teff = 160 kK, log(g) = 8
1781  *
1782  * Since Cloudy checks the scaling elsewhere there is no need to re-scale
1783  * things here, but this inaccuracy should be kept in mind since it could
1784  * indicate problems with the flux distribution */
1785 
1786  InterpolateRectGrid( &grid, val, Tlow, Thigh );
1787 
1788  FreeGrid( &grid );
1789  return rfield.nupper;
1790 }
1791 
1792 /* WMBASICCompile rebin WMBASIC stellar models to match energy grid of code */
1794 {
1795  /* these contain frequencies for the major absorption edges */
1796  realnum Edges[3];
1797 
1798  bool lgFail = false;
1799 
1800  DEBUG_ENTRY( "WMBASICCompile()" );
1801 
1802  fprintf( ioQQQ, " WMBASICCompile on the job.\n" );
1803 
1804  /* define the major absorption edges that require special attention during rebinning */
1805  Edges[0] = 0.99946789f;
1806  Edges[1] = 1.8071406f;
1807  Edges[2] = 3.9996377f;
1808 
1810 
1811  if( lgFileReadable( "wmbasic.ascii", pc, as ) && !lgValidBinFile( "wmbasic.mod", pc, as ) )
1812  lgFail = lgCompileAtmosphere( "wmbasic.ascii", "wmbasic.mod", Edges, 3L, pc );
1813  return lgFail;
1814 }
1815 
1816 /* WMBASICInterpolate read in and interpolate on WMBASIC grid of hot star atmospheres */
1817 long WMBASICInterpolate(double val[], /* val[nval] */
1818  long *nval,
1819  long *ndim,
1820  bool lgList,
1821  double *Tlow,
1822  double *Thigh)
1823 {
1825 
1826  DEBUG_ENTRY( "WMBASICInterpolate()" );
1827 
1828  grid.name = "wmbasic.mod";
1829  grid.scheme = AS_DATA_OPTIONAL;
1830  /* identification of this atmosphere set, used in
1831  * the Cloudy output, *must* be 12 characters long */
1832  grid.ident = " WMBASIC";
1833  /* the Cloudy command needed to recompile the binary model file */
1834  grid.command = "COMPILE STARS";
1835 
1836  InitGrid( &grid, lgList );
1837 
1838  CheckVal( &grid, val, nval, ndim );
1839 
1840  InterpolateRectGrid( &grid, val, Tlow, Thigh );
1841 
1842  FreeGrid( &grid );
1843  return rfield.nupper;
1844 }
1845 
1846 /* CompileAtmosphereCoStar rebin costar stellar atmospheres to match cloudy energy grid,
1847  * called by the compile stars command */
1848 STATIC bool lgCompileAtmosphereCoStar(const char chFNameIn[],
1849  const char chFNameOut[],
1850  const realnum Edges[], /* Edges[nedges] */
1851  long nedges,
1852  process_counter& pc)
1853 {
1854  char chLine[INPUT_LINE_LENGTH];
1855  char names[MDIM][MNAM+1];
1856  int32 val[7];
1857  uint32 uval[2];
1858  long int i, j, nskip, nModels, nWL;
1859 
1860  /* these will be malloced into large work arrays*/
1861  realnum *StarEner = NULL, *StarFlux = NULL, *CloudyFlux = NULL;
1862  /* this will hold all the model parameters */
1863  mpp *telg = NULL;
1864 
1865  FILE *ioIN; /* used for input */
1866  FILE *ioOUT; /* used for output */
1867 
1868  DEBUG_ENTRY( "CompileAtmosphereCoStar()" );
1869 
1870  /* This is a program to re-bin the costar stellar model spectra to match the
1871  * Cloudy grid. For short wavelengths I will use a power law extrapolation
1872  * of the model values (which should be falling rapidly) if needed. At long
1873  * wavelengths I will assume Rayleigh-Jeans from the last stellar model point
1874  * to extrapolate to 1 cm wavelength. */
1875 
1876  /* This version uses power-law interpolation between the points of the stellar model. */
1877 
1878  /* read the original data file obtained off the web,
1879  * open as read only */
1880  try
1881  {
1882  ioIN = open_data( chFNameIn, "r", AS_LOCAL_ONLY );
1883  }
1884  catch( cloudy_exit )
1885  {
1886  goto error;
1887  }
1888  fprintf( ioQQQ, " CompileAtmosphereCoStar got %s.\n", chFNameIn );
1889 
1890  /* get first line and see how many more to skip */
1891  if( read_whole_line( chLine, (int)sizeof(chLine), ioIN ) == NULL )
1892  {
1893  fprintf( ioQQQ, " CompileAtmosphereCoStar fails reading nskip.\n" );
1894  goto error;
1895  }
1896  sscanf( chLine, "%li", &nskip );
1897 
1898  /* now skip the header information */
1899  for( i=0; i < nskip; ++i )
1900  {
1901  if( read_whole_line( chLine, (int)sizeof(chLine), ioIN ) == NULL )
1902  {
1903  fprintf( ioQQQ, " CompileAtmosphereCoStar fails skipping header.\n" );
1904  goto error;
1905  }
1906  }
1907 
1908  /* now get number of models and number of wavelengths */
1909  if( read_whole_line( chLine, (int)sizeof(chLine), ioIN ) == NULL )
1910  {
1911  fprintf( ioQQQ, " CompileAtmosphereCoStar fails reading nModels, nWL.\n" );
1912  goto error;
1913  }
1914  sscanf( chLine, "%li%li", &nModels, &nWL );
1915 
1916  if( nModels <= 0 || nWL <= 0 )
1917  {
1918  fprintf( ioQQQ, " CompileAtmosphereCoStar scanned off impossible values for nModels=%li or nWL=%li\n",
1919  nModels, nWL );
1920  goto error;
1921  }
1922 
1923  /* allocate space for the stellar parameters */
1924  telg = (mpp *)CALLOC( (size_t)nModels, sizeof(mpp) );
1925 
1926  /* get all model parameters for the atmospheres */
1927  for( i=0; i < nModels; ++i )
1928  {
1929  if( read_whole_line( chLine, (int)sizeof(chLine), ioIN ) == NULL )
1930  {
1931  fprintf( ioQQQ, " CompileAtmosphereCoStar fails reading model parameters.\n" );
1932  goto error;
1933  }
1934  /* first letter on line is indicator of grid */
1935  telg[i].chGrid = chLine[0];
1936  /* get the model id number */
1937  sscanf( chLine+1, "%i", &telg[i].modid );
1938  /* get the temperature */
1939  sscanf( chLine+23, "%lg", &telg[i].par[0] );
1940  /* get the surface gravity */
1941  sscanf( chLine+31, "%lg", &telg[i].par[1] );
1942  /* get the ZAMS mass */
1943  sscanf( chLine+7, "%lg", &telg[i].par[2] );
1944  /* get the model age */
1945  sscanf( chLine+15, "%lg", &telg[i].par[3] );
1946 
1947  /* the code in parse_table.cpp implicitly depends on this! */
1948  ASSERT( telg[i].par[2] > 10. );
1949  ASSERT( telg[i].par[3] > 10. );
1950 
1951  /* convert ZAMS masses to logarithms */
1952  telg[i].par[2] = log10(telg[i].par[2]);
1953  }
1954 
1955  /* this will be the file we create, that will be read to compute models,
1956  * open to write binary */
1957  try
1958  {
1959  ioOUT = open_data( chFNameOut, "wb", AS_LOCAL_ONLY );
1960  }
1961  catch( cloudy_exit )
1962  {
1963  goto error;
1964  }
1965 
1966  /* >>chng 05 dec 01, use int32 instead of long so that binary file is compatible for 32/64-bit code */
1967  /* >>chng 05 dec 01, add several parameters to binary file, change sequence, PvH */
1968  /* >>chng 06 jun 10, add several parameters to binary file, change sequence, PvH */
1969  val[0] = (int32)VERSION_BIN;
1970  val[1] = (int32)MDIM;
1971  val[2] = (int32)MNAM;
1972  val[3] = (int32)2;
1973  val[4] = (int32)4;
1974  val[5] = (int32)nModels;
1975  val[6] = (int32)rfield.nupper;
1976  uval[0] = sizeof(val) + sizeof(uval) + sizeof(names) + nModels*sizeof(mpp); /* nOffset */
1977  uval[1] = rfield.nupper*sizeof(realnum); /* nBlocksize */
1978 
1979  strncpy( names[0], "Teff\0\0", MNAM+1 );
1980  strncpy( names[1], "log(g)", MNAM+1 );
1981  strncpy( names[2], "log(M)", MNAM+1 );
1982  strncpy( names[3], "Age\0\0\0", MNAM+1 );
1983 
1984  if( fwrite( val, sizeof(val), 1, ioOUT ) != 1 ||
1985  fwrite( uval, sizeof(uval), 1, ioOUT ) != 1 ||
1986  fwrite( names, sizeof(names), 1, ioOUT ) != 1 ||
1987  /* write out the array of {Teff,log(g)} pairs */
1988  fwrite( telg, sizeof(mpp), (size_t)nModels, ioOUT ) != (size_t)nModels ||
1989  /* write out the cloudy energy grid for later sanity checks */
1990  fwrite( rfield.AnuOrg, (size_t)uval[1], 1, ioOUT ) != 1 )
1991  {
1992  fprintf( ioQQQ, " CompileAtmosphereCoStar failed writing header of output file.\n" );
1993  goto error;
1994  }
1995 
1996  /* MALLOC some workspace */
1997  StarEner = (realnum *)MALLOC( sizeof(realnum)*nWL );
1998  StarFlux = (realnum *)MALLOC( sizeof(realnum)*nWL );
1999  CloudyFlux = (realnum *)MALLOC( (size_t)uval[1] );
2000 
2001  fprintf( ioQQQ, " Compiling: " );
2002 
2003  /* get the star data */
2004  for( i=0; i < nModels; ++i )
2005  {
2006  /* get number to skip */
2007  if( read_whole_line( chLine, (int)sizeof(chLine), ioIN ) == NULL )
2008  {
2009  fprintf( ioQQQ, " CompileAtmosphereCoStar fails reading the skip to next spectrum.\n" );
2010  goto error;
2011  }
2012  sscanf( chLine, "%li", &nskip );
2013 
2014  for( j=0; j < nskip; ++j )
2015  {
2016  if( read_whole_line( chLine, (int)sizeof(chLine), ioIN ) == NULL )
2017  {
2018  fprintf( ioQQQ, " CompileAtmosphereCoStar fails doing the skip.\n" );
2019  goto error;
2020  }
2021  }
2022 
2023  /* now read in the wavelength and flux for this star, read in
2024  * backwards since we want to be in increasing energy order rather
2025  * than wavelength */
2026  for( j=nWL-1; j >= 0; --j )
2027  {
2028  if( read_whole_line( chLine, (int)sizeof(chLine), ioIN ) == NULL )
2029  {
2030  fprintf( ioQQQ, " CompileAtmosphereCoStar fails reading the spectral data.\n" );
2031  goto error;
2032  }
2033  double help1, help2;
2034  sscanf( chLine, "%lg %lg", &help1, &help2 );
2035 
2036  /* continuum flux was log, convert to linear, also do
2037  * conversion from "astrophysical" flux to F_nu in cgs units */
2038  StarFlux[j] = (realnum)(PI*pow(10.,help2));
2039  /* StarEner was in Angstroms, convert to Ryd */
2040  StarEner[j] = (realnum)(RYDLAM/help1);
2041 
2042  /* sanity check */
2043  if( j < nWL-1 )
2044  ASSERT( StarEner[j] < StarEner[j+1] );
2045  }
2046 
2047  /* this will do the heavy lifting, and define arrays used below,
2048  * NB the lowest energy point in these grids appears to be bogus.
2049  * tell rebin about nWL-1 */
2050  RebinAtmosphere(nWL-1, StarEner+1, StarFlux+1, CloudyFlux, nedges, Edges );
2051 
2052  /* write the continuum out as a binary file */
2053  if( fwrite( CloudyFlux, (size_t)uval[1], 1, ioOUT ) != 1 )
2054  {
2055  fprintf( ioQQQ, " CompileAtmosphereCoStar failed writing star flux.\n" );
2056  goto error;
2057  }
2058 
2059  fprintf( ioQQQ, "." );
2060  fflush( ioQQQ );
2061  }
2062 
2063  fprintf( ioQQQ, " done.\n" );
2064 
2065  fclose( ioIN );
2066  fclose( ioOUT );
2067 
2068  FREE_CHECK( telg );
2069  FREE_CHECK( StarEner );
2070  FREE_CHECK( StarFlux );
2071  FREE_CHECK( CloudyFlux );
2072 
2073  fprintf( ioQQQ, "\n CompileAtmosphereCoStar completed ok.\n\n" );
2074  ++pc.nOK;
2075  return false;
2076 
2077 error:
2078  FREE_SAFE( telg );
2079  FREE_SAFE( StarEner );
2080  FREE_SAFE( StarFlux );
2081  FREE_SAFE( CloudyFlux );
2082  ++pc.nFail;
2083  return true;
2084 }
2085 
2086 /* InterpolateGridCoStar read in and interpolate on costar grid of windy O atmospheres */
2087 STATIC void InterpolateGridCoStar(const stellar_grid *grid, /* struct with all the grid parameters */
2088  const double val[], /* val[0]: Teff for imode = 1,2; M_ZAMS for imode = 3;
2089  * age for imode = 4 */
2090  /* val[1]: nmodid for imode = 1; log(g) for imode = 2;
2091  * age for imode = 3; M_ZAMS for imode = 4 */
2092  double *val0_lo,
2093  double *val0_hi)
2094 {
2095  long i, j, k, nmodid, off, ptr;
2096  long *indloTr, *indhiTr, useTr[2];
2097  long indlo[2], indhi[2], index[2];
2098  realnum *ValTr;
2099  double lval[2], aval[4];
2100 
2101  DEBUG_ENTRY( "InterpolateGridCoStar()" );
2102 
2103  switch( grid->imode )
2104  {
2105  case IM_COSTAR_TEFF_MODID:
2106  case IM_COSTAR_TEFF_LOGG:
2107  lval[0] = val[0];
2108  lval[1] = val[1];
2109  off = 0;
2110  break;
2111  case IM_COSTAR_MZAMS_AGE:
2112  lval[0] = log10(val[0]); /* use log10(M_ZAMS) internally */
2113  lval[1] = val[1];
2114  off = 2;
2115  break;
2116  case IM_COSTAR_AGE_MZAMS:
2117  /* swap parameters, hence mimic IM_COSTAR_MZAMS_AGE */
2118  lval[0] = log10(val[1]); /* use log10(M_ZAMS) internally */
2119  lval[1] = val[0];
2120  off = 2;
2121  break;
2122  default:
2123  fprintf( ioQQQ, " InterpolateGridCoStar called with insane value for imode: %d.\n", grid->imode );
2124  cdEXIT(EXIT_FAILURE);
2125  }
2126 
2127  nmodid = (long)(lval[1]+0.5);
2128 
2130 
2131  /* read in the saved cloudy energy scale so we can confirm this is a good image */
2133 
2134 # if DEBUGPRT
2135  /* check whether the models in the grid have the correct effective temperature */
2136  ValidateGrid( grid, 0.005 );
2137 # endif
2138 
2139  /* now allocate some temp workspace */
2140  ValTr = (realnum *)MALLOC( sizeof(realnum)*grid->nTracks );
2141  indloTr = (long *)MALLOC( sizeof(long)*grid->nTracks );
2142  indhiTr = (long *)MALLOC( sizeof(long)*grid->nTracks );
2143 
2144  /* first do horizontal search, i.e. search along individual tracks */
2145  for( j=0; j < grid->nTracks; j++ )
2146  {
2147  if( grid->imode == IM_COSTAR_TEFF_MODID )
2148  {
2149  if( grid->trackLen[j] >= nmodid ) {
2150  index[0] = nmodid - 1;
2151  index[1] = j;
2152  ptr = grid->jval[JIndex(grid,index)];
2153  indloTr[j] = ptr;
2154  indhiTr[j] = ptr;
2155  ValTr[j] = (realnum)grid->telg[ptr].par[off];
2156  }
2157  else
2158  {
2159  indloTr[j] = -2;
2160  indhiTr[j] = -2;
2161  ValTr[j] = -FLT_MAX;
2162  }
2163  }
2164  else
2165  {
2166  FindHCoStar( grid, j, lval[1], off, ValTr, indloTr, indhiTr );
2167  }
2168  }
2169 
2170 # if DEBUGPRT
2171  for( j=0; j < grid->nTracks; j++ )
2172  {
2173  if( indloTr[j] >= 0 )
2174  printf( "track %c: models %c%d, %c%d, val %g\n",
2175  (char)('A'+j), grid->telg[indloTr[j]].chGrid, grid->telg[indloTr[j]].modid,
2176  grid->telg[indhiTr[j]].chGrid, grid->telg[indhiTr[j]].modid, ValTr[j]);
2177  }
2178 # endif
2179 
2180  /* now do vertical search, i.e. interpolate between tracks */
2181  FindVCoStar( grid, lval[0], ValTr, useTr );
2182 
2183  /* This should only happen when InterpolateGridCoStar is called in non-optimizing mode,
2184  * when optimizing InterpolateGridCoStar should report back to optimize_func()...
2185  * The fact that FindVCoStar allows interpolation between non-adjoining tracks
2186  * should guarantee that this will not happen. */
2187  if( useTr[0] < 0 )
2188  {
2189  fprintf( ioQQQ, " The parameters for the requested CoStar model are out of range.\n" );
2190  cdEXIT(EXIT_FAILURE);
2191  }
2192 
2193  ASSERT( useTr[0] >= 0 && useTr[0] < grid->nTracks );
2194  ASSERT( useTr[1] >= 0 && useTr[1] < grid->nTracks );
2195  ASSERT( indloTr[useTr[0]] >= 0 && indloTr[useTr[0]] < (int)grid->nmods );
2196  ASSERT( indhiTr[useTr[0]] >= 0 && indhiTr[useTr[0]] < (int)grid->nmods );
2197  ASSERT( indloTr[useTr[1]] >= 0 && indloTr[useTr[1]] < (int)grid->nmods );
2198  ASSERT( indhiTr[useTr[1]] >= 0 && indhiTr[useTr[1]] < (int)grid->nmods );
2199 
2200 # if DEBUGPRT
2201  printf( "interpolate between tracks %c and %c\n", (char)('A'+useTr[0]), (char)('A'+useTr[1]) );
2202 # endif
2203 
2204  indlo[0] = indloTr[useTr[0]];
2205  indhi[0] = indhiTr[useTr[0]];
2206  indlo[1] = indloTr[useTr[1]];
2207  indhi[1] = indhiTr[useTr[1]];
2208 
2209  InterpolateModelCoStar( grid, lval, aval, indlo, indhi, index, 0, off, rfield.tslop[rfield.nspec] );
2210 
2211  for( i=0; i < rfield.nupper; i++ )
2212  {
2213  rfield.tslop[rfield.nspec][i] = (realnum)pow((realnum)10.f,rfield.tslop[rfield.nspec][i]);
2214  if( rfield.tslop[rfield.nspec][i] < 1e-37 )
2215  rfield.tslop[rfield.nspec][i] = 0.;
2216  }
2217 
2218  if( false )
2219  {
2220  FILE *ioBUG = fopen( "interpolated.txt", "w" );
2221  for( k=0; k < rfield.nupper; ++k )
2222  fprintf( ioBUG, "%e %e\n", rfield.tNuRyd[rfield.nspec][k], rfield.tslop[rfield.nspec][k] );
2223  fclose( ioBUG );
2224  }
2225 
2226  /* sanity check: see whether this model has the correct effective temperature */
2227  if( ! lgValidModel( rfield.tNuRyd[rfield.nspec], rfield.tslop[rfield.nspec], aval[0], 0.05 ) )
2228  TotalInsanity();
2229 
2230  /* set limits for optimizer */
2231  SetLimits( grid, lval[0], NULL, NULL, useTr, ValTr, val0_lo, val0_hi );
2232 
2233  /* now write some final info */
2234  if( called.lgTalk )
2235  {
2236  fprintf( ioQQQ, " * c<< FINAL: T_eff = %7.1f, ", aval[0] );
2237  fprintf( ioQQQ, "log(g) = %4.2f, M(ZAMS) = %5.1f, age = ", aval[1], pow(10.,aval[2]) );
2238  fprintf( ioQQQ, PrintEfmt("%8.2e",aval[3]) );
2239  fprintf( ioQQQ, " >>> *\n" );
2240  }
2241 
2242  FREE_CHECK( indhiTr );
2243  FREE_CHECK( indloTr );
2244  FREE_CHECK( ValTr );
2245  return;
2246 }
2247 
2248 /* find which models to use for interpolation along a given evolutionary track */
2250  long track,
2251  double par2, /* requested log(g) or age */
2252  long off, /* determines which parameter to match 0 -> log(g), 2 -> age */
2253  realnum *ValTr, /* ValTr[track]: Teff/log(M) value for interpolated model along track */
2254  long *indloTr, /* indloTr[track]: model number for first model used in interpolation */
2255  long *indhiTr) /* indhiTr[track]: model number for second model used in interpolation */
2256 {
2257  long index[2], j, mod1, mod2;
2258 
2259  DEBUG_ENTRY( "FindHCoStar()" );
2260 
2261  indloTr[track] = -2;
2262  indhiTr[track] = -2;
2263  ValTr[track] = -FLT_MAX;
2264 
2265  index[1] = track;
2266 
2267  for( j=0; j < grid->trackLen[track]; j++ )
2268  {
2269  index[0] = j;
2270  mod1 = grid->jval[JIndex(grid,index)];
2271 
2272  /* do we have an exact match ? */
2273  if( fabs(par2-grid->telg[mod1].par[off+1]) <= 10.*FLT_EPSILON*fabs(grid->telg[mod1].par[off+1]) )
2274  {
2275  indloTr[track] = mod1;
2276  indhiTr[track] = mod1;
2277  ValTr[track] = (realnum)grid->telg[mod1].par[off];
2278  return;
2279  }
2280  }
2281 
2282  for( j=0; j < grid->trackLen[track]-1; j++ )
2283  {
2284  index[0] = j;
2285  mod1 = grid->jval[JIndex(grid,index)];
2286  index[0] = j+1;
2287  mod2 = grid->jval[JIndex(grid,index)];
2288 
2289  /* do we interpolate ? */
2290  if( (par2 - grid->telg[mod1].par[off+1])*(par2 - grid->telg[mod2].par[off+1]) < 0. )
2291  {
2292  double frac;
2293 
2294  indloTr[track] = mod1;
2295  indhiTr[track] = mod2;
2296  frac = (par2 - grid->telg[mod2].par[off+1])/
2297  (grid->telg[mod1].par[off+1] - grid->telg[mod2].par[off+1]);
2298  ValTr[track] = (realnum)(frac*grid->telg[mod1].par[off] +
2299  (1.-frac)*grid->telg[mod2].par[off] );
2300  break;
2301  }
2302  }
2303  return;
2304 }
2305 
2306 /* find which tracks to use for interpolation in between tracks */
2308  double par1, /* requested Teff or ZAMS mass */
2309  realnum *ValTr, /* internal workspace */
2310  long useTr[]) /* useTr[0]: track number for first track to be used in interpolation
2311  * (i.e., 0 means 'A', etc.)
2312  * useTr[1]: track number for second track to be used in interpolation
2313  * NOTE: FindVCoStar raises a flag when interpolating between non-adjoining
2314  * tracks, i.e. when (useTr[1]-useTr[0]) > 1 */
2315 {
2316  long j;
2317 
2318  DEBUG_ENTRY( "FindVCoStar()" );
2319 
2320  useTr[0] = -1;
2321  useTr[1] = -1;
2322 
2323  for( j=0; j < grid->nTracks; j++ )
2324  {
2325  /* do we have an exact match ? */
2326  if( ValTr[j] != -FLT_MAX && fabs(par1-(double)ValTr[j]) <= 10.*FLT_EPSILON*fabs(ValTr[j]) )
2327  {
2328  useTr[0] = j;
2329  useTr[1] = j;
2330  break;
2331  }
2332  }
2333 
2334  if( useTr[0] >= 0 )
2335  {
2336  return;
2337  }
2338 
2339  for( j=0; j < grid->nTracks-1; j++ )
2340  {
2341  if( ValTr[j] != -FLT_MAX )
2342  {
2343  long int i,j2;
2344 
2345  /* find next valid track */
2346  j2 = 0;
2347  for( i = j+1; i < grid->nTracks; i++ )
2348  {
2349  if( ValTr[i] != -FLT_MAX )
2350  {
2351  j2 = i;
2352  break;
2353  }
2354  }
2355 
2356  /* do we interpolate ? */
2357  if( j2 > 0 && ((realnum)par1-ValTr[j])*((realnum)par1-ValTr[j2]) < 0.f )
2358  {
2359  useTr[0] = j;
2360  useTr[1] = j2;
2361  break;
2362  }
2363  }
2364  }
2365 
2366  /* raise caution when we interpolate between non-adjoining tracks */
2367  continuum.lgCoStarInterpolationCaution = ( useTr[1]-useTr[0] > 1 );
2368  return;
2369 }
2370 
2371 /* Make a listing of all the models in the CoStar grid */
2373 {
2374  long index[2], maxlen, n;
2375 
2376  DEBUG_ENTRY( "CoStarListModels()" );
2377 
2378  maxlen = 0;
2379  for( n=0; n < grid->nTracks; n++ )
2380  maxlen = MAX2( maxlen, grid->trackLen[n] );
2381 
2382  fprintf( ioQQQ, "\n" );
2383  fprintf( ioQQQ, " Track\\Index |" );
2384  for( n = 0; n < maxlen; n++ )
2385  fprintf( ioQQQ, " %5ld ", n+1 );
2386  fprintf( ioQQQ, "\n" );
2387  fprintf( ioQQQ, "--------------|" );
2388  for( n = 0; n < maxlen; n++ )
2389  fprintf( ioQQQ, "----------------" );
2390  fprintf( ioQQQ, "\n" );
2391 
2392  for( index[1]=0; index[1] < grid->nTracks; ++index[1] )
2393  {
2394  long ptr;
2395  double Teff, alogg, Mass;
2396 
2397  fprintf( ioQQQ, " %c", (char)('A'+index[1]) );
2398  index[0] = 0;
2399  ptr = grid->jval[JIndex(grid,index)];
2400  Mass = pow(10.,grid->telg[ptr].par[2]);
2401  fprintf( ioQQQ, " (%3.0f Msol) |", Mass );
2402 
2403  for( index[0]=0; index[0] < grid->trackLen[index[1]]; ++index[0] )
2404  {
2405  ptr = grid->jval[JIndex(grid,index)];
2406  Teff = grid->telg[ptr].par[0];
2407  alogg = grid->telg[ptr].par[1];
2408  fprintf( ioQQQ, " (%6.1f,%4.2f)", Teff, alogg );
2409  }
2410  fprintf( ioQQQ, "\n" );
2411  }
2412  return;
2413 }
2414 
2415 /* RauchInitializeSub does the actual work of preparing the ascii file */
2416 STATIC int RauchInitializeSub(const char chFName[],
2417  const char chSuff[],
2418  const mpp telg[],
2419  long nmods,
2420  long n,
2421  long ngrids,
2422  const double par2[], /* par2[ngrids] */
2423  int format)
2424 {
2425  char chLine[INPUT_LINE_LENGTH]; /* used for getting input lines */
2426 
2427  FILE *ioOut, /* pointer to output file we came here to create*/
2428  *ioIn; /* pointer to input files we will read */
2429 
2430  long int i, j;
2431 
2432  double *wavl, *fluxes;
2433 
2434  DEBUG_ENTRY( "RauchInitializeSub()" );
2435 
2436  /* grab some space for the wavelengths and fluxes */
2437  wavl = (double *)MALLOC( sizeof(double)*NRAUCH);
2438  fluxes = (double *)MALLOC( sizeof(double)*NRAUCH);
2439 
2440  try
2441  {
2442  if( n == 1 )
2443  ioOut = open_data( chFName, "w", AS_LOCAL_ONLY );
2444  else
2445  ioOut = open_data( chFName, "a", AS_LOCAL_ONLY );
2446  }
2447  catch( cloudy_exit )
2448  {
2449  goto error;
2450  }
2451 
2452  if( n == 1 )
2453  {
2454  fprintf( ioOut, " %ld\n", VERSION_ASCII );
2455  fprintf( ioOut, " %d\n", ( ngrids == 1 ? 2 : 3 ) );
2456  fprintf( ioOut, " %d\n", ( ngrids == 1 ? 2 : 3 ) );
2457  fprintf( ioOut, " Teff\n" );
2458  fprintf( ioOut, " log(g)\n" );
2459  if( ngrids == 2 )
2460  fprintf( ioOut, " log(Z)\n" );
2461  else if( ngrids == 11 )
2462  fprintf( ioOut, " f(He)\n" );
2463  /* NB - this is based on the assumption that each of the planes in the cubic grid is the same */
2464  fprintf( ioOut, " %ld\n", nmods*ngrids );
2465  fprintf( ioOut, " %d\n", NRAUCH );
2466  /* Rauch models give the wavelength in Angstrom */
2467  fprintf( ioOut, " lambda\n" );
2468  /* conversion factor to Angstrom */
2469  fprintf( ioOut, " %.8e\n", 1. );
2470  /* Rauch models give the "Astrophysical" flux F_lambda in erg/cm^2/s/cm */
2471  fprintf( ioOut, " F_lambda\n" );
2472  /* the factor PI*1e-8 is needed to convert to "regular" flux in erg/cm^2/s/Angstrom */
2473  fprintf( ioOut, " %.8e\n", PI*1.e-8 );
2474  /* NB - this is based on the assumption that each of the planes in the cubic grid is the same */
2475  for( j=0; j < ngrids; j++ )
2476  {
2477  /* write out the {Teff,log(g)} grid */
2478  for( i=0; i < nmods; i++ )
2479  {
2480  if( ngrids == 1 )
2481  fprintf( ioOut, " %.0f %.1f", telg[i].par[0], telg[i].par[1] );
2482  else
2483  fprintf( ioOut, " %.0f %.1f %.1f", telg[i].par[0], telg[i].par[1], par2[j] );
2484  if( ((i+1)%4) == 0 )
2485  fprintf( ioOut, "\n" );
2486  }
2487  if( (i%4) != 0 )
2488  fprintf( ioOut, "\n" );
2489  }
2490 
2491  fprintf( ioQQQ, " Writing: " );
2492  }
2493 
2494  for( i=0; i < nmods; i++ )
2495  {
2496  /* must create name of next stellar atmosphere */
2497  if( format == 1 )
2498  sprintf( chLine, "%7.7ld_%2ld", (long)(telg[i].par[0]+0.5), (long)(10.*telg[i].par[1]+0.5) );
2499  else if( format == 2 )
2500  sprintf( chLine, "%7.7ld_%.2f", (long)(telg[i].par[0]+0.5), telg[i].par[1] );
2501  else
2502  {
2503  fprintf( ioQQQ, " insanity in RauchInitializeSub\n" );
2504  ShowMe();
2505  cdEXIT(EXIT_FAILURE);
2506  }
2507  string chFileName( chLine );
2508  chFileName += chSuff;
2509  /* now open next stellar atmosphere for reading*/
2510  try
2511  {
2512  ioIn = open_data( chFileName.c_str(), "r", AS_LOCAL_ONLY );
2513  }
2514  catch( cloudy_exit )
2515  {
2516  goto error;
2517  }
2518 
2519  /* get first line */
2520  j = 0;
2521  if( read_whole_line( chLine, (int)sizeof(chLine), ioIn ) == NULL )
2522  {
2523  fprintf( ioQQQ, " RauchInitializeSub error in atmosphere file %4ld%4ld\n",
2524  i, j );
2525  goto error;
2526  }
2527  /* >>chng 02 nov 20, now keep reading them until don't hit the *
2528  * since number of comments may change */
2529  while( chLine[0] == '*' )
2530  {
2531  if( read_whole_line( chLine, (int)sizeof(chLine), ioIn ) == NULL )
2532  {
2533  fprintf( ioQQQ, " RauchInitializeSub error in atmosphere file %4ld%4ld\n",
2534  i, j );
2535  goto error;
2536  }
2537  ++j;
2538  }
2539 
2540  for( j=0; j < NRAUCH; j++ )
2541  {
2542  double ttemp, wl;
2543  /* get the input line */
2544  /* >>chng 02 nov 20, don't reread very first line image since we got it above */
2545  if( j > 0 )
2546  {
2547  if(read_whole_line( chLine, (int)sizeof(chLine), ioIn )==NULL )
2548  {
2549  fprintf( ioQQQ, " RauchInitializeSub error in atmosphere file %4ld%4ld\n",
2550  i, j );
2551  goto error;
2552  }
2553  }
2554 
2555  /* scan off wavelength and flux)*/
2556  if( sscanf( chLine, "%lf %le", &wl, &ttemp ) != 2 )
2557  {
2558  fprintf( ioQQQ, " RauchInitializeSub error in atmosphere file %4ld%4ld\n",
2559  i, j );
2560  goto error;
2561  }
2562 
2563  if( i == 0 )
2564  wavl[j] = wl;
2565  else
2566  {
2567  /* check if this model is on the same wavelength grid as the first */
2568  if( fabs(wavl[j]-wl) > 10.*DBL_EPSILON*wl )
2569  {
2570  fprintf( ioQQQ, " RauchInitializeSub error in atmosphere file %4ld%4ld\n",
2571  i, j );
2572  goto error;
2573  }
2574  }
2575  fluxes[j] = ttemp;
2576  }
2577 
2578  /* finished - close the unit */
2579  fclose(ioIn);
2580 
2581  /* now write to output file */
2582  if( i == 0 && n == 1 )
2583  {
2584  /* wavelength grid is the same for all models, so write only once */
2585  for( j=0; j < NRAUCH; j++ )
2586  {
2587  fprintf( ioOut, " %.4e", wavl[j] );
2588  /* want to have 5 numbers per line */
2589  if( ((j+1)%5) == 0 )
2590  fprintf( ioOut, "\n" );
2591  }
2592  /* need to throw a newline if we did not end on an exact line */
2593  if( (j%5) != 0 )
2594  fprintf( ioOut, "\n" );
2595  }
2596 
2597  for( j=0; j < NRAUCH; j++ )
2598  {
2599  fprintf( ioOut, " %.4e", fluxes[j] );
2600  /* want to have 5 numbers per line */
2601  if( ((j+1)%5) == 0 )
2602  fprintf( ioOut, "\n" );
2603  }
2604  /* need to throw a newline if we did not end on an exact line */
2605  if( (j%5) != 0 )
2606  fprintf( ioOut, "\n" );
2607 
2608  /* print to screen to show progress */
2609  fprintf( ioQQQ, "." );
2610  fflush( ioQQQ );
2611  }
2612 
2613  if( n == ngrids )
2614  fprintf( ioQQQ, " done.\n" );
2615 
2616  fclose(ioOut);
2617 
2618  /* free the space grabbed above */
2619  FREE_CHECK( fluxes );
2620  FREE_CHECK( wavl );
2621  return 0;
2622 
2623 error:
2624  FREE_CHECK( fluxes );
2625  FREE_CHECK( wavl );
2626  return 1;
2627 }
2628 
2629 /* lgCompileAtmosphere does the actual rebinning onto the Cloudy grid and writes the binary file */
2630 /* >>chng 01 feb 12, added return value to indicate success (0) or failure (1) */
2631 STATIC bool lgCompileAtmosphere(const char chFNameIn[],
2632  const char chFNameOut[],
2633  const realnum Edges[], /* Edges[nedges] */
2634  long nedges,
2635  process_counter& pc)
2636 {
2637  FILE *ioIN; /* used for input */
2638  FILE *ioOUT; /* used for output */
2639 
2640  char chDataType[11];
2641  char names[MDIM][MNAM+1];
2642 
2643  bool lgFreqX, lgFreqY, lgFlip;
2644  int32 val[7];
2645  uint32 uval[2];
2646  long int i, imod, version, nd, ndim, npar, nmods, ngrid;
2647 
2648  /* these will be malloced into large work arrays */
2649  realnum *StarEner = NULL, *StarFlux = NULL, *CloudyFlux = NULL, *scratch = NULL;
2650 
2651  double convert_wavl, convert_flux;
2652 
2653  mpp *telg = NULL;
2654 
2655  DEBUG_ENTRY( "lgCompileAtmosphere()" );
2656 
2657  try
2658  {
2659  ioIN = open_data( chFNameIn, "r", AS_LOCAL_ONLY );
2660  }
2661  catch( cloudy_exit )
2662  {
2663  goto error;
2664  }
2665  fprintf( ioQQQ, " lgCompileAtmosphere got %s.\n", chFNameIn );
2666 
2667  /* read version number */
2668  if( fscanf( ioIN, "%ld", &version ) != 1 )
2669  {
2670  fprintf( ioQQQ, " lgCompileAtmosphere failed reading VERSION.\n" );
2671  goto error;
2672  }
2673 
2674  if( version != VERSION_ASCII )
2675  {
2676  fprintf( ioQQQ, " lgCompileAtmosphere: there is a version number mismatch in"
2677  " the ascii atmosphere file: %s.\n", chFNameIn );
2678  fprintf( ioQQQ, " lgCompileAtmosphere: Please recreate this file or download the"
2679  " latest version following the instructions on the Cloudy website.\n" );
2680  goto error;
2681  }
2682 
2683  /* >>chng 06 jun 10, read the dimension of the grid, PvH */
2684  if( fscanf( ioIN, "%ld", &ndim ) != 1 || ndim <= 0 || ndim > MDIM )
2685  {
2686  fprintf( ioQQQ, " lgCompileAtmosphere failed reading valid dimension of grid.\n" );
2687  goto error;
2688  }
2689 
2690  /* >>chng 06 jun 12, read the number of model parameters, PvH */
2691  if( fscanf( ioIN, "%ld", &npar ) != 1 || npar <= 0 || npar < ndim || npar > MDIM )
2692  {
2693  fprintf( ioQQQ, " lgCompileAtmosphere failed reading valid no. of model parameters.\n" );
2694  goto error;
2695  }
2696 
2697  /* make sure valgrind doesn't trip on the binary write of this array */
2698  memset( names, '\0', MDIM*(MNAM+1) );
2699 
2700  for( nd=0; nd < npar; nd++ )
2701  {
2702  if( fscanf( ioIN, "%6s", names[nd] ) != 1 )
2703  {
2704  fprintf( ioQQQ, " lgCompileAtmosphere failed reading parameter label.\n" );
2705  goto error;
2706  }
2707  }
2708  if( strcmp( names[0], "Teff" ) != 0 && strcmp( names[0], "Age" ) != 0 )
2709  {
2710  fprintf( ioQQQ, " lgCompileAtmosphere: first parameter must be \"Teff\" or \"Age\".\n" );
2711  goto error;
2712  }
2713  if( ndim > 1 && strcmp( names[0], "Age" ) == 0 )
2714  {
2715  fprintf( ioQQQ, " First parameter is \"Age\", but ndim is not 1.\n" );
2716  goto error;
2717  }
2718  if( ndim >= 2 && strcmp( names[1], "log(g)" ) != 0 )
2719  {
2720  fprintf( ioQQQ, " lgCompileAtmosphere: second parameter must be \"log(g)\".\n" );
2721  goto error;
2722  }
2723 
2724  /* >>chng 05 nov 18, read the following extra parameters from the ascii file, PvH */
2725  if( fscanf( ioIN, "%ld", &nmods ) != 1 || nmods <= 0 )
2726  {
2727  fprintf( ioQQQ, " lgCompileAtmosphere failed reading valid number of models.\n" );
2728  goto error;
2729  }
2730 
2731  if( fscanf( ioIN, "%ld", &ngrid ) != 1 || ngrid <= 1 )
2732  {
2733  fprintf( ioQQQ, " lgCompileAtmosphere failed reading valid number of grid points.\n" );
2734  goto error;
2735  }
2736 
2737  /* read data type for wavelengths, allowed values are lambda, nu */
2738  if( fscanf( ioIN, "%10s", chDataType ) != 1 )
2739  {
2740  fprintf( ioQQQ, " lgCompileAtmosphere failed reading wavl DataType string.\n" );
2741  goto error;
2742  }
2743 
2744  if( strcmp( chDataType, "lambda" ) == 0 )
2745  lgFreqX = false;
2746  else if( strcmp( chDataType, "nu" ) == 0 )
2747  lgFreqX = true;
2748  else {
2749  fprintf( ioQQQ, " lgCompileAtmosphere found illegal wavl DataType: %s.\n", chDataType );
2750  goto error;
2751  }
2752 
2753  if( fscanf( ioIN, "%le", &convert_wavl ) != 1 || convert_wavl <= 0. )
2754  {
2755  fprintf( ioQQQ, " lgCompileAtmosphere failed reading valid wavl conversion factor.\n" );
2756  goto error;
2757  }
2758 
2759  /* read data type for flux, allowed values F_lambda, H_lambda, F_nu, H_nu */
2760  if( fscanf( ioIN, "%10s", chDataType ) != 1 )
2761  {
2762  fprintf( ioQQQ, " lgCompileAtmosphere failed reading flux DataType string.\n" );
2763  goto error;
2764  }
2765 
2766  if( strcmp( chDataType, "F_lambda" ) == 0 || strcmp( chDataType, "H_lambda" ) == 0 )
2767  lgFreqY = false;
2768  else if( strcmp( chDataType, "F_nu" ) == 0 || strcmp( chDataType, "H_nu" ) == 0 )
2769  lgFreqY = true;
2770  else {
2771  fprintf( ioQQQ, " lgCompileAtmosphere found illegal flux DataType: %s.\n", chDataType );
2772  goto error;
2773  }
2774 
2775  if( fscanf( ioIN, "%le", &convert_flux ) != 1 || convert_flux <= 0. )
2776  {
2777  fprintf( ioQQQ, " lgCompileAtmosphere failed reading valid flux conversion factor.\n" );
2778  goto error;
2779  }
2780 
2781  telg = (mpp *)CALLOC( (size_t)nmods, sizeof(mpp) );
2782 
2783  for( i=0; i < nmods; i++ )
2784  {
2785  for( nd=0; nd < npar; nd++ )
2786  {
2787  if( fscanf( ioIN, "%le", &telg[i].par[nd] ) != 1 )
2788  {
2789  fprintf( ioQQQ, " lgCompileAtmosphere failed reading valid model parameter.\n" );
2790  goto error;
2791  }
2792  }
2793  if( telg[i].par[0] <= 0. )
2794  {
2795  fprintf( ioQQQ, " lgCompileAtmosphere failed reading valid %s.\n", names[0] );
2796  goto error;
2797  }
2798  }
2799 
2800  try
2801  {
2802  ioOUT = open_data( chFNameOut, "wb", AS_LOCAL_ONLY );
2803  }
2804  catch( cloudy_exit )
2805  {
2806  goto error;
2807  }
2808 
2809  /* >>chng 05 nov 16, use int32 instead of long so that binary file is compatible for 32/64-bit code */
2810  /* >>chng 05 nov 18, add several parameters to binary file, change sequence, PvH */
2811  /* >>chng 06 jun 10, add several parameters to binary file, change sequence, PvH */
2812  val[0] = (int32)VERSION_BIN;
2813  val[1] = (int32)MDIM;
2814  val[2] = (int32)MNAM;
2815  val[3] = (int32)ndim;
2816  val[4] = (int32)npar;
2817  val[5] = (int32)nmods;
2818  val[6] = (int32)rfield.nupper;
2819  uval[0] = sizeof(val) + sizeof(uval) + sizeof(names) + nmods*sizeof(mpp); /* nOffset */
2820  uval[1] = rfield.nupper*sizeof(realnum); /* nBlocksize */
2821 
2822  if( fwrite( val, sizeof(val), 1, ioOUT ) != 1 ||
2823  fwrite( uval, sizeof(uval), 1, ioOUT ) != 1 ||
2824  fwrite( names, sizeof(names), 1, ioOUT ) != 1 ||
2825  /* write out the array of {Teff,log(g)} pairs */
2826  fwrite( telg, sizeof(mpp), (size_t)nmods, ioOUT ) != (size_t)nmods ||
2827  /* write out the cloudy energy grid for later sanity checks */
2828  fwrite( rfield.AnuOrg, (size_t)uval[1], 1, ioOUT ) != 1 )
2829  {
2830  fprintf( ioQQQ, " lgCompileAtmosphere failed writing header of output file.\n" );
2831  goto error;
2832  }
2833 
2834  /* MALLOC some workspace */
2835  StarEner = (realnum *)MALLOC( sizeof(realnum)*ngrid );
2836  scratch = (realnum *)MALLOC( sizeof(realnum)*ngrid );
2837  StarFlux = (realnum *)MALLOC( sizeof(realnum)*ngrid );
2838  CloudyFlux = (realnum *)MALLOC( (size_t)uval[1] );
2839 
2840  /* read wavelength grid */
2841  for( i=0; i < ngrid; i++ )
2842  {
2843  double help;
2844  if( fscanf( ioIN, "%lg", &help ) != 1 )
2845  {
2846  fprintf( ioQQQ, " lgCompileAtmosphere failed reading wavelength.\n" );
2847  goto error;
2848  }
2849  /* this conversion makes sure that scratch[i] is
2850  * either wavelength in Angstrom or frequency in Hz */
2851  scratch[i] = (realnum)(help*convert_wavl);
2852 
2853  ASSERT( scratch[i] > 0.f );
2854  }
2855 
2856  lgFlip = ( !lgFreqX && scratch[0] < scratch[1] ) || ( lgFreqX && scratch[0] > scratch[1] );
2857 
2858  /* convert continuum over to increasing frequency in Ryd */
2859  for( i=0; i < ngrid; i++ )
2860  {
2861  /* convert scratch[i] to frequency in Ryd */
2862  if( lgFreqX )
2863  scratch[i] /= (realnum)FR1RYD;
2864  else
2865  scratch[i] = (realnum)(RYDLAM/scratch[i]);
2866 
2867  if( lgFlip )
2868  StarEner[ngrid-i-1] = scratch[i];
2869  else
2870  StarEner[i] = scratch[i];
2871  }
2872 
2873  ASSERT( StarEner[0] > 0.f );
2874  /* make sure the array is in ascending order */
2875  for( i=1; i < ngrid; i++ )
2876  ASSERT( StarEner[i] > StarEner[i-1] );
2877 
2878  fprintf( ioQQQ, " Compiling: " );
2879 
2880  for( imod=0; imod < nmods; imod++ )
2881  {
2882  const realnum CONVERT_FNU = (realnum)(1.e8*SPEEDLIGHT/POW2(FR1RYD));
2883 
2884  /* now read the stellar fluxes */
2885  for( i=0; i < ngrid; i++ )
2886  {
2887  double help;
2888  if( fscanf( ioIN, "%lg", &help ) != 1 )
2889  {
2890  fprintf( ioQQQ, " lgCompileAtmosphere failed reading star flux.\n" );
2891  goto error;
2892  }
2893  /* this conversion makes sure that scratch[i] is either
2894  * F_nu in erg/cm^2/s/Hz or F_lambda in erg/cm^2/s/A */
2895  scratch[i] = (realnum)(help*convert_flux);
2896 
2897  /* this can underflow on the Wien tail */
2898  ASSERT( scratch[i] >= 0.f );
2899  }
2900 
2901  for( i=0; i < ngrid; i++ )
2902  {
2903  if( lgFlip )
2904  StarFlux[ngrid-i-1] = scratch[i];
2905  else
2906  StarFlux[i] = scratch[i];
2907  }
2908 
2909  for( i=0; i < ngrid; i++ )
2910  {
2911  /* this converts to F_nu in erg/cm^2/s/Hz */
2912  if( !lgFreqY )
2913  StarFlux[i] *= CONVERT_FNU/POW2(StarEner[i]);
2914  ASSERT( StarFlux[i] >= 0.f );
2915  }
2916 
2917  /* the re-binned values are returned in the "CloudyFlux" array */
2918  RebinAtmosphere( ngrid, StarEner, StarFlux, CloudyFlux, nedges, Edges );
2919 
2920  /* write the continuum out as a binary file */
2921  if( fwrite( CloudyFlux, (size_t)uval[1], 1, ioOUT ) != 1 )
2922  {
2923  fprintf( ioQQQ, " lgCompileAtmosphere failed writing star flux.\n" );
2924  goto error;
2925  }
2926 
2927  fprintf( ioQQQ, "." );
2928  fflush( ioQQQ );
2929  }
2930 
2931  fprintf( ioQQQ, " done.\n" );
2932 
2933  fclose(ioIN);
2934  fclose(ioOUT);
2935 
2936  /* now free up the memory we claimed */
2937  FREE_CHECK( CloudyFlux );
2938  FREE_CHECK( StarFlux );
2939  FREE_CHECK( StarEner );
2940  FREE_CHECK( scratch );
2941  FREE_CHECK( telg );
2942 
2943  fprintf( ioQQQ, " lgCompileAtmosphere completed ok.\n\n" );
2944  ++pc.nOK;
2945  return false;
2946 
2947 error:
2948  FREE_SAFE( CloudyFlux );
2949  FREE_SAFE( StarFlux );
2950  FREE_SAFE( StarEner );
2951  FREE_SAFE( scratch );
2952  FREE_SAFE( telg );
2953  ++pc.nFail;
2954  return true;
2955 }
2956 
2958  bool lgList)
2959 {
2960  long nd;
2961  int32 version, mdim, mnam;
2962 
2963  DEBUG_ENTRY( "InitGrid()" );
2964 
2965  try
2966  {
2967  grid->ioIN = open_data( grid->name.c_str(), "rb", grid->scheme );
2968  }
2969  catch( cloudy_exit )
2970  {
2971  /* something went wrong */
2972  /* NB NB - DO NOT CHANGE THE FOLLOWING ERROR MESSAGE! checkall.pl picks it up */
2973  fprintf( ioQQQ, " Error: stellar atmosphere file not found.\n" );
2974  fprintf(ioQQQ , "\n\n If the path is set then it is possible that the stellar atmosphere data files do not exist.\n");
2975  fprintf(ioQQQ , " Have the stellar data files been downloaded and compiled with the COMPILE STARS command?\n");
2976  fprintf(ioQQQ , " If you are simply running the test suite and do not need the stellar continua then you should simply ignore this failure\n");
2977  cdEXIT(EXIT_FAILURE);
2978  }
2979 
2980  /* >>chng 01 oct 17, add version and size to this array */
2981  if( fread( &version, sizeof(version), 1, grid->ioIN ) != 1 ||
2982  fread( &mdim, sizeof(mdim), 1, grid->ioIN ) != 1 ||
2983  fread( &mnam, sizeof(mnam), 1, grid->ioIN ) != 1 ||
2984  fread( &grid->ndim, sizeof(grid->ndim), 1, grid->ioIN ) != 1 ||
2985  fread( &grid->npar, sizeof(grid->npar), 1, grid->ioIN ) != 1 ||
2986  fread( &grid->nmods, sizeof(grid->nmods), 1, grid->ioIN ) != 1 ||
2987  fread( &grid->ngrid, sizeof(grid->ngrid), 1, grid->ioIN ) != 1 ||
2988  fread( &grid->nOffset, sizeof(grid->nOffset), 1, grid->ioIN ) != 1 ||
2989  fread( &grid->nBlocksize, sizeof(grid->nBlocksize), 1, grid->ioIN ) != 1 )
2990  {
2991  fprintf( ioQQQ, " InitGrid failed reading header.\n" );
2992  cdEXIT(EXIT_FAILURE);
2993  }
2994 
2995  /* do some sanity checks */
2996  if( version != VERSION_BIN )
2997  {
2998  fprintf( ioQQQ, " InitGrid: there is a version mismatch between"
2999  " the compiled atmospheres file I expected and the one I found.\n" );
3000  fprintf( ioQQQ, " InitGrid: Please recompile the stellar"
3001  " atmospheres file with the command: %s.\n", grid->command );
3002  cdEXIT(EXIT_FAILURE);
3003  }
3004 
3005  if( mdim != MDIM || mnam != MNAM )
3006  {
3007  fprintf( ioQQQ, " InitGrid: the compiled atmospheres file is produced"
3008  " with an incompatible version of Cloudy.\n" );
3009  fprintf( ioQQQ, " InitGrid: Please recompile the stellar"
3010  " atmospheres file with the command: %s.\n", grid->command );
3011  cdEXIT(EXIT_FAILURE);
3012  }
3013 
3014  ASSERT( grid->ndim > 0 && grid->ndim <= MDIM );
3015  ASSERT( grid->npar >= grid->ndim && grid->npar <= MDIM );
3016  ASSERT( grid->nmods > 0 );
3017  ASSERT( grid->ngrid > 0 );
3018  ASSERT( grid->nOffset > 0 );
3019  ASSERT( grid->nBlocksize > 0 );
3020 
3021  rfield.nupper = grid->ngrid;
3022 
3023  if( fread( &grid->names, sizeof(grid->names), 1, grid->ioIN ) != 1 )
3024  {
3025  fprintf( ioQQQ, " InitGrid failed reading names array.\n" );
3026  cdEXIT(EXIT_FAILURE);
3027  }
3028 
3029  grid->telg = (mpp *)MALLOC( sizeof(mpp)*grid->nmods );
3030  grid->val = (double **)MALLOC( sizeof(double*)*grid->ndim );
3031  for( nd=0; nd < grid->ndim; nd++ )
3032  {
3033  grid->val[nd] = (double *)MALLOC( sizeof(double)*grid->nmods );
3034  }
3035  grid->nval = (long *)MALLOC( sizeof(long)*grid->ndim );
3036 
3037  if( fread( grid->telg, sizeof(mpp), grid->nmods, grid->ioIN ) != (size_t)grid->nmods )
3038  {
3039  fprintf( ioQQQ, " InitGrid failed reading model parameter block.\n" );
3040  cdEXIT(EXIT_FAILURE);
3041  }
3042 
3043 # ifdef SEEK_END
3044  /* sanity check: does the file have the correct length ? */
3045  /* NOTE: this operation is not necessarily supported by all operating systems
3046  * but if the preprocessor symbol SEEK_END exists it is assumed to be supported */
3047  int res = fseek( grid->ioIN, 0, SEEK_END );
3048  if( res == 0 )
3049  {
3050  long End = ftell( grid->ioIN );
3051  long Expected = grid->nOffset + (grid->nmods+1)*grid->nBlocksize;
3052  if( End != Expected )
3053  {
3054  fprintf( ioQQQ, " InitGrid: Problem performing sanity check for size of binary file.\n" );
3055  fprintf( ioQQQ, " InitGrid: I expected to find %ld bytes, but actually found %ld bytes.\n",
3056  Expected, End );
3057  fprintf( ioQQQ, " InitGrid: Please recompile the stellar"
3058  " atmospheres file with the command: %s.\n", grid->command );
3059  cdEXIT(EXIT_FAILURE);
3060  }
3061  }
3062 # endif
3063 
3064  InitIndexArrays( grid, lgList );
3065 
3066  /* set default interpolation mode */
3067  grid->imode = IM_RECT_GRID;
3068  /* these are only used by CoStar grids */
3069  grid->trackLen = NULL;
3070  grid->nTracks = 0;
3071  grid->jval = NULL;
3072  return;
3073 }
3074 
3075 /* check whether a binary atmosphere exists and is up-to-date */
3076 STATIC bool lgValidBinFile(const char *binName, process_counter& pc, access_scheme scheme)
3077 {
3078  int32 version, mdim, mnam;
3080 
3081  DEBUG_ENTRY( "lgValidBinFile()" );
3082 
3083  grid.name = binName;
3084 
3085  if( (grid.ioIN = open_data( grid.name.c_str(), "rb", scheme )) == NULL )
3086  return false;
3087 
3088  if( fread( &version, sizeof(version), 1, grid.ioIN ) != 1 ||
3089  fread( &mdim, sizeof(mdim), 1, grid.ioIN ) != 1 ||
3090  fread( &mnam, sizeof(mnam), 1, grid.ioIN ) != 1 ||
3091  fread( &grid.ndim, sizeof(grid.ndim), 1, grid.ioIN ) != 1 ||
3092  fread( &grid.npar, sizeof(grid.npar), 1, grid.ioIN ) != 1 ||
3093  fread( &grid.nmods, sizeof(grid.nmods), 1, grid.ioIN ) != 1 ||
3094  fread( &grid.ngrid, sizeof(grid.ngrid), 1, grid.ioIN ) != 1 ||
3095  fread( &grid.nOffset, sizeof(grid.nOffset), 1, grid.ioIN ) != 1 ||
3096  fread( &grid.nBlocksize, sizeof(grid.nBlocksize), 1, grid.ioIN ) != 1 )
3097  {
3098  fclose( grid.ioIN );
3099  return false;
3100  }
3101 
3102  /* do some sanity checks */
3103  if( version != VERSION_BIN || mdim != MDIM || mnam != MNAM )
3104  {
3105  fclose( grid.ioIN );
3106  return false;
3107  }
3108 
3109 # ifdef SEEK_END
3110  /* sanity check: does the file have the correct length ? */
3111  /* NOTE: this operation is not necessarily supported by all operating systems
3112  * but if the preprocessor symbol SEEK_END exists it is assumed to be supported */
3113  int res = fseek( grid.ioIN, 0, SEEK_END );
3114  if( res == 0 )
3115  {
3116  long End = ftell( grid.ioIN );
3117  long Expected = grid.nOffset + (grid.nmods+1)*grid.nBlocksize;
3118  if( End != Expected )
3119  {
3120  fclose( grid.ioIN );
3121  return false;
3122  }
3123  }
3124 # endif
3125 
3126  fclose( grid.ioIN );
3127  ++pc.notProcessed; // the file is up-to-date -> no processing
3128  return true;
3129 }
3130 
3131 /* check whether a ascii atmosphere file exists and is up-to-date */
3132 STATIC bool lgValidAsciiFile(const char *ascName, access_scheme scheme)
3133 {
3134  long version;
3135  FILE *ioIN;
3136 
3137  DEBUG_ENTRY( "lgValidAsciiFile()" );
3138 
3139  /* can we read the file? */
3140  if( (ioIN = open_data( ascName, "r", scheme )) == NULL )
3141  return false;
3142 
3143  /* check version number */
3144  if( fscanf( ioIN, "%ld", &version ) != 1 || version != VERSION_ASCII )
3145  {
3146  fclose( ioIN );
3147  return false;
3148  }
3149 
3150  fclose( ioIN );
3151  return true;
3152 }
3153 
3154 /* sort CoStar models according to track and index number, store indices in grid->jval[] */
3155 STATIC void InitGridCoStar(stellar_grid *grid) /* the grid parameters */
3156 {
3157  char track;
3158  bool lgFound;
3159  long i, index[2];
3160 
3161  DEBUG_ENTRY( "InitGridCoStar()" );
3162 
3163  ASSERT( grid->ndim == 2 );
3164  ASSERT( grid->jlo != NULL && grid->jhi != NULL );
3165 
3166  grid->jval = grid->jlo;
3167  FREE_CHECK( grid->jhi );
3168  grid->jlo = grid->jhi = NULL;
3169 
3170  /* invalidate contents set by InitGrid first */
3171  memset( grid->jval, 0xff, (size_t)(grid->nval[0]*grid->nval[1]*sizeof(long)) );
3172 
3173  grid->trackLen = (long *)CALLOC( (size_t)grid->nmods, sizeof(long) );
3174 
3175  index[1] = 0;
3176  while( true )
3177  {
3178  index[0] = 0;
3179  track = (char)('A'+index[1]);
3180  do
3181  {
3182  lgFound = false;
3183  for( i=0; i < grid->nmods; i++ )
3184  {
3185  if( grid->telg[i].chGrid == track && grid->telg[i].modid == index[0]+1 )
3186  {
3187  grid->jval[JIndex(grid,index)] = i;
3188  ++index[0];
3189  lgFound = true;
3190  break;
3191  }
3192  }
3193  }
3194  while( lgFound );
3195 
3196  if( index[0] == 0 )
3197  break;
3198 
3199  grid->trackLen[index[1]] = index[0];
3200  ++index[1];
3201  }
3202 
3203  grid->nTracks = index[1];
3204  return;
3205 }
3206 
3208  double val[], /* val[ndim] */
3209  long *nval,
3210  long *ndim)
3211 {
3212  DEBUG_ENTRY( "CheckVal()" );
3213 
3214  if( *ndim == 0 )
3215  *ndim = (long)grid->ndim;
3216  if( *ndim == 2 && *nval == 1 )
3217  {
3218  /* default gravity is maximum gravity */
3219  val[*nval] = grid->val[1][grid->nval[1]-1];
3220  ++(*nval);
3221  }
3222  if( *ndim != (long)grid->ndim )
3223  {
3224  fprintf( ioQQQ, " A %ld-dim grid was requested, but a %ld-dim grid was found.\n",
3225  *ndim, (long)grid->ndim );
3226  cdEXIT(EXIT_FAILURE);
3227  }
3228  if( *nval < *ndim )
3229  {
3230  fprintf( ioQQQ, " A %ld-dim grid was requested, but only %ld parameters were entered.\n",
3231  *ndim, *nval );
3232  cdEXIT(EXIT_FAILURE);
3233  }
3234 }
3235 
3237  const double val[], /* val[ndim] */
3238  double *Tlow,
3239  double *Thigh)
3240 {
3241  bool lgInvalid;
3242  long int i,
3243  *indlo,
3244  *indhi,
3245  *index,
3246  k,
3247  nd;
3248  double *aval;
3249 
3250  DEBUG_ENTRY( "InterpolateRectGrid()" );
3251 
3252  /* create some space */
3253  indlo = (long *)MALLOC((size_t)(grid->ndim*sizeof(long)) );
3254  indhi = (long *)MALLOC((size_t)(grid->ndim*sizeof(long)) );
3255  index = (long *)MALLOC((size_t)(grid->ndim*sizeof(long)) );
3256  aval = (double *)MALLOC((size_t)(grid->npar*sizeof(double)) );
3257 
3259  ASSERT( grid->nBlocksize == rfield.nupper*sizeof(realnum) );
3260 
3261  /* save energy scale for check against code's in conorm (scale not yet defined when this routine called) */
3263 
3264 # if DEBUGPRT
3265  /* check whether the models have the correct effective temperature, for debugging only */
3266  ValidateGrid( grid, 0.02 );
3267 # endif
3268 
3269  /* now generate pointers for models to use */
3270  for( nd=0; nd < grid->ndim; nd++ )
3271  {
3272  FindIndex( grid->val[nd], grid->nval[nd], val[nd], &indlo[nd], &indhi[nd], &lgInvalid );
3273  if( lgInvalid )
3274  {
3275  fprintf( ioQQQ,
3276  " Requested parameter %s = %.2f is not within the range %.2f to %.2f\n",
3277  grid->names[nd], val[nd], grid->val[nd][0], grid->val[nd][grid->nval[nd]-1] );
3278  cdEXIT(EXIT_FAILURE);
3279  }
3280  }
3281 
3282  InterpolateModel( grid, val, aval, indlo, indhi, index, grid->ndim, rfield.tslop[rfield.nspec], IS_UNDEFINED );
3283 
3284  /* print the parameters of the interpolated model */
3285  if( called.lgTalk )
3286  {
3287  if( grid->npar == 1 )
3288  fprintf( ioQQQ,
3289  " * c<< FINAL: %6s = %13.2f"
3290  " >>> *\n",
3291  grid->names[0], aval[0] );
3292  else if( grid->npar == 2 )
3293  fprintf( ioQQQ,
3294  " * c<< FINAL: %6s = %10.2f"
3295  " %6s = %8.5f >>> *\n",
3296  grid->names[0], aval[0], grid->names[1], aval[1] );
3297  else if( grid->npar == 3 )
3298  fprintf( ioQQQ,
3299  " * c<< FINAL: %6s = %7.0f"
3300  " %6s = %5.2f %6s = %5.2f >>> *\n",
3301  grid->names[0], aval[0], grid->names[1], aval[1],
3302  grid->names[2], aval[2] );
3303  else if( grid->npar >= 4 )
3304  {
3305  fprintf( ioQQQ,
3306  " * c<< FINAL: %4s = %7.0f"
3307  " %6s = %4.2f %6s = %5.2f %6s = ",
3308  grid->names[0], aval[0], grid->names[1], aval[1],
3309  grid->names[2], aval[2], grid->names[3] );
3310  fprintf( ioQQQ, PrintEfmt( "%9.2e", aval[3] ) );
3311  fprintf( ioQQQ, " >>> *\n" );
3312  }
3313  }
3314 
3315  for( i=0; i < rfield.nupper; i++ )
3316  {
3317  rfield.tslop[rfield.nspec][i] = (realnum)pow((realnum)10.f,rfield.tslop[rfield.nspec][i]);
3318  if( rfield.tslop[rfield.nspec][i] < 1e-37 )
3319  rfield.tslop[rfield.nspec][i] = 0.;
3320  }
3321 
3322  if( false )
3323  {
3324  FILE *ioBUG = fopen( "interpolated.txt", "w" );
3325  for( k=0; k < rfield.nupper; ++k )
3326  fprintf( ioBUG, "%e %e\n", rfield.tNuRyd[rfield.nspec][k], rfield.tslop[rfield.nspec][k] );
3327  fclose( ioBUG );
3328  }
3329 
3330  if( strcmp( grid->names[0], "Teff" ) == 0 )
3331  {
3332  if( ! lgValidModel( rfield.tNuRyd[rfield.nspec], rfield.tslop[rfield.nspec], val[0], 0.10 ) )
3333  TotalInsanity();
3334  }
3335 
3336  /* set limits for optimizer */
3337  SetLimits( grid, val[0], indlo, indhi, NULL, NULL, Tlow, Thigh );
3338 
3339  FREE_CHECK( aval );
3340  FREE_CHECK( index );
3341  FREE_CHECK( indhi );
3342  FREE_CHECK( indlo );
3343  return;
3344 }
3345 
3347 {
3348  long i;
3349 
3350  DEBUG_ENTRY( "FreeGrid()" );
3351 
3352  /* this was opened/allocated in InitGrid and subsidiaries,
3353  * this should become a destructor in C++ */
3354  fclose( grid->ioIN );
3355  FREE_CHECK( grid->telg );
3356  for( i = 0; i < grid->ndim; i++ )
3357  FREE_CHECK( grid->val[i] );
3358  FREE_CHECK( grid->val );
3359  FREE_CHECK( grid->nval );
3360  FREE_SAFE( grid->jlo );
3361  FREE_SAFE( grid->jhi );
3362  FREE_SAFE( grid->trackLen )
3363  FREE_SAFE( grid->jval );
3364  return;
3365 }
3366 
3368  const double val[],
3369  double aval[],
3370  const long indlo[],
3371  const long indhi[],
3372  long index[],
3373  long nd,
3374  realnum flux1[],
3375  IntStage stage)
3376 {
3377  bool lgDryRun;
3378  long i, ind, j;
3379 
3380  DEBUG_ENTRY( "InterpolateModel()" );
3381 
3382  --nd;
3383 
3384  lgDryRun = ( flux1 == NULL );
3385 
3386  if( nd < 0 )
3387  {
3388  long n = JIndex(grid,index);
3389  if( stage == IS_FIRST )
3390  ind = ( grid->jlo[n] >= 0 ) ? grid->jlo[n] : grid->jhi[n];
3391  else if( stage == IS_SECOND )
3392  ind = ( grid->jhi[n] >= 0 ) ? grid->jhi[n] : grid->jlo[n];
3393  else if( grid->ndim == 1 )
3394  /* in this case grid->jlo[n] and grid->jhi[n] should be identical */
3395  ind = grid->jlo[n];
3396  else
3397  TotalInsanity();
3398 
3399  if( ind < 0 )
3400  {
3401  fprintf( ioQQQ, " The requested interpolation could not be completed, sorry.\n" );
3402  fprintf( ioQQQ, " No suitable match was found for a model with" );
3403  for( i=0; i < grid->ndim; i++ )
3404  fprintf( ioQQQ, " %s=%.6g ", grid->names[i], grid->val[i][index[i]] );
3405  fprintf( ioQQQ, "\n" );
3406  cdEXIT(EXIT_FAILURE);
3407  }
3408 
3409  for( i=0; i < grid->npar; i++ )
3410  aval[i] = grid->telg[ind].par[i];
3411 
3412  if( !lgDryRun )
3413  {
3414  for( i=0; i < grid->ndim && called.lgTalk; i++ )
3415  {
3416  if( fabs(grid->val[i][index[i]]-aval[i]) > 10.*DBL_EPSILON*fabs(aval[i]) )
3417  {
3418  fprintf( ioQQQ, " No exact match was found for a model with" );
3419  for( j=0; j < grid->ndim; j++ )
3420  fprintf( ioQQQ, " %s=%.6g ", grid->names[j], grid->val[j][index[j]] );
3421  fprintf( ioQQQ, "- using the following model instead:\n" );
3422  break;
3423  }
3424  }
3425 
3426  GetModel( grid, ind, flux1, lgVERBOSE, lgTAKELOG );
3427  }
3428  }
3429  else
3430  {
3431  realnum *flux2;
3432  double *aval2;
3433 
3434 # ifndef NDEBUG
3435  const realnum SECURE = 10.f*FLT_EPSILON;
3436 # endif
3437 
3438  flux2 = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
3439  aval2 = (double*)MALLOC((size_t)(grid->npar*sizeof(double)) );
3440 
3441  /* Interpolation is carried out first in Teff (i.e., nd == 0), then the
3442  * parameter with nd == 1 (log(g)), etc. One or two atmosphere models
3443  * are read depending on whether the parameter was matched exactly or not.
3444  * If needed, logarithmic interpolation is done.
3445  */
3446 
3447  if( nd == 1 )
3448  stage = IS_FIRST;
3449 
3450  index[nd] = indlo[nd];
3451  InterpolateModel( grid, val, aval, indlo, indhi, index, nd, flux1, stage );
3452 
3453  if( nd == 1 )
3454  stage = IS_SECOND;
3455 
3456  index[nd] = indhi[nd];
3457  InterpolateModel( grid, val, aval2, indlo, indhi, index, nd, NULL, stage );
3458 
3459  /* the test MUST be "greater than", otherwise it will fail for aval[nd] == aval2[nd] == 0 */
3460  if( fabs(aval2[nd]-aval[nd]) > 10.*DBL_EPSILON*fabs(aval[nd]) )
3461  {
3462  double fr1, fr2, fc1 = 0., fc2 = 0.;
3463 
3464  if( !lgDryRun )
3465  InterpolateModel( grid, val, aval2, indlo, indhi, index, nd, flux2, stage );
3466 
3467  fr1 = (aval2[nd]-val[nd])/(aval2[nd]-aval[nd]);
3468  /* when interpolating in log(g) it can happen that fr1 is outside the range 0 .. 1
3469  * this can be the case when the requested log(g) was not present in the grid
3470  * and it had to be approximated by another model. In this case do not extrapolate */
3471  if( nd == 1 )
3472  fr1 = MIN2( MAX2( fr1, 0. ), 1. );
3473  fr2 = 1. - fr1;
3474 
3475  ASSERT( 0.-SECURE <= fr1 && fr1 <= 1.+SECURE );
3476 
3477  if( !lgDryRun )
3478  {
3479 # if DEBUGPRT
3480  fprintf( ioQQQ, "interpolation nd=%ld fr1=%g\n", nd, fr1 );
3481 # endif
3482 
3483  /* special treatment for high-temperature Rauch models */
3484  if( nd == 0 && strcmp( grid->names[nd], "Teff" ) == 0 )
3485  {
3486  /* The following is an approximate scaling to use for the range of
3487  * temperatures above 200000 K in the H-Ca Rauch models where the
3488  * temperature steps are large and thus the interpolations are over
3489  * large ranges. For the lower temperatures I assume that there is
3490  * no need for this.
3491  *
3492  * It should be remembered that this interpolation is not exact, and
3493  * the possible error at high temperatures might be large enough to
3494  * matter. (Kevin Volk)
3495  */
3496  fc1 = ( val[nd] > 200000. ) ? log10(val[nd]/grid->val[nd][indlo[nd]])*4. : 0.;
3497  fc2 = ( val[nd] > 200000. ) ? log10(val[nd]/grid->val[nd][indhi[nd]])*4. : 0.;
3498  }
3499 
3500  for( i=0; i < rfield.nupper; ++i )
3501  flux1[i] = (realnum)(fr1*(flux1[i]+fc1) + fr2*(flux2[i]+fc2));
3502  }
3503 
3504  for( i=0; i < grid->npar; i++ )
3505  aval[i] = fr1*aval[i] + fr2*aval2[i];
3506  }
3507 
3508  FREE_CHECK( aval2 );
3509  FREE_CHECK( flux2 );
3510  }
3511  return;
3512 }
3513 
3515  const double val[],
3516  double aval[],
3517  const long indlo[],
3518  const long indhi[],
3519  long index[],
3520  long nd,
3521  long off,
3522  realnum flux1[])
3523 {
3524  long i, ind;
3525 
3526  DEBUG_ENTRY( "InterpolateModelCoStar()" );
3527 
3528  if( nd == 2 )
3529  {
3530  ind = ( index[1] == 0 ) ? indlo[index[0]] : indhi[index[0]];
3531 
3532  GetModel( grid, ind, flux1, lgVERBOSE, lgTAKELOG );
3533 
3534  for( i=0; i < grid->npar; i++ )
3535  aval[i] = grid->telg[ind].par[i];
3536  }
3537  else
3538  {
3539  bool lgSkip;
3540 # ifndef NDEBUG
3541  const realnum SECURE = 10.f*FLT_EPSILON;
3542 # endif
3543 
3544  /* Interpolation is carried out first along evolutionary tracks, then
3545  * in between evolutionary tracks. Between 1 and 4 atmosphere models are read
3546  * depending on whether the parameter/track was matched exactly or not.
3547  */
3548 
3549  index[nd] = 0;
3550  InterpolateModelCoStar( grid, val, aval, indlo, indhi, index, nd+1, off, flux1 );
3551 
3552  lgSkip = ( nd == 1 ) ? ( indhi[index[0]] == indlo[index[0]] ) :
3553  ( indlo[0] == indlo[1] && indhi[0] == indhi[1] );
3554 
3555  if( ! lgSkip )
3556  {
3557  realnum *flux2;
3558  double fr1, fr2, *aval2;
3559 
3560  flux2 = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
3561  aval2 = (double*)MALLOC((size_t)(grid->npar*sizeof(double)) );
3562 
3563  index[nd] = 1;
3564  InterpolateModelCoStar( grid, val, aval2, indlo, indhi, index, nd+1, off, flux2 );
3565 
3566  fr1 = (aval2[nd+off]-val[nd])/(aval2[nd+off]-aval[nd+off]);
3567  fr2 = 1. - fr1;
3568 
3569 # if DEBUGPRT
3570  fprintf( ioQQQ, "interpolation nd=%ld fr1=%g\n", nd, fr1 );
3571 # endif
3572 
3573  ASSERT( 0.-SECURE <= fr1 && fr1 <= 1.+SECURE );
3574 
3575  for( i=0; i < rfield.nupper; ++i )
3576  flux1[i] = (realnum)(fr1*flux1[i] + fr2*flux2[i]);
3577 
3578  for( i=0; i < grid->npar; i++ )
3579  aval[i] = fr1*aval[i] + fr2*aval2[i];
3580 
3581  FREE_CHECK( aval2 );
3582  FREE_CHECK( flux2 );
3583  }
3584  }
3585  return;
3586 }
3587 
3589  long ind,
3590  realnum flux[],
3591  bool lgTalk,
3592  bool lgTakeLog)
3593 {
3594  long i;
3595 
3596  DEBUG_ENTRY( "GetModel()" );
3597 
3598  /* add 1 to account for frequency grid that is stored in front of all the atmospheres */
3599  ind++;
3600 
3601  /* make sure ident is exactly 12 characters long, otherwise output won't fit */
3602  ASSERT( strlen(grid->ident) == 12 );
3603  /* ind == 0 is the frequency grid, ind == 1 .. nmods are the atmosphere models */
3604  ASSERT( ind >= 0 && ind <= grid->nmods );
3605 
3606  /* skip over ind stars */
3607  /* >>chng 01 oct 18, add nOffset */
3608  if( fseek( grid->ioIN, (long)(ind*grid->nBlocksize+grid->nOffset), SEEK_SET ) != 0 )
3609  {
3610  fprintf( ioQQQ, " Error seeking atmosphere %ld\n", ind );
3611  cdEXIT(EXIT_FAILURE);
3612  }
3613 
3614  if( fread( flux, 1, grid->nBlocksize, grid->ioIN ) != grid->nBlocksize )
3615  {
3616  fprintf( ioQQQ, " Error trying to read atmosphere %ld\n", ind );
3617  cdEXIT(EXIT_FAILURE);
3618  }
3619 
3620  /* print the parameters of the atmosphere model */
3621  if( called.lgTalk && lgTalk )
3622  {
3623  /* ind-1 below since telg doesn't have the entry for the frequency grid */
3624  if( grid->npar == 1 )
3625  fprintf( ioQQQ,
3626  " * c<< %s model%5ld read. "
3627  " %6s = %13.2f >>> *\n",
3628  grid->ident, ind, grid->names[0], grid->telg[ind-1].par[0] );
3629  else if( grid->npar == 2 )
3630  fprintf( ioQQQ,
3631  " * c<< %s model%5ld read. "
3632  " %6s = %10.2f %6s = %8.5f >>> *\n",
3633  grid->ident, ind, grid->names[0], grid->telg[ind-1].par[0],
3634  grid->names[1], grid->telg[ind-1].par[1] );
3635  else if( grid->npar == 3 )
3636  fprintf( ioQQQ,
3637  " * c<< %s model%5ld read. "
3638  " %6s=%7.0f %6s=%5.2f %6s=%5.2f >>> *\n",
3639  grid->ident, ind, grid->names[0], grid->telg[ind-1].par[0],
3640  grid->names[1], grid->telg[ind-1].par[1],
3641  grid->names[2], grid->telg[ind-1].par[2] );
3642  else if( grid->npar >= 4 )
3643  {
3644  fprintf( ioQQQ,
3645  " * c< %s mdl%4ld"
3646  " %4s=%5.0f %6s=%4.2f %6s=%5.2f %6s=",
3647  grid->ident, ind, grid->names[0], grid->telg[ind-1].par[0],
3648  grid->names[1], grid->telg[ind-1].par[1],
3649  grid->names[2], grid->telg[ind-1].par[2], grid->names[3] );
3650  fprintf( ioQQQ, PrintEfmt( "%9.2e", grid->telg[ind-1].par[3] ) );
3651  fprintf( ioQQQ, " >> *\n" );
3652  }
3653  }
3654 
3655  if( lgTakeLog )
3656  {
3657  /* convert to logs since we will interpolate in log flux */
3658  for( i=0; i < rfield.nupper; ++i )
3659  flux[i] = (realnum)log10( MAX2( 1e-37, (double)flux[i] ) );
3660  }
3661  return;
3662 }
3663 
3665  double val,
3666  const long indlo[],
3667  const long indhi[],
3668  const long useTr[],
3669  const realnum ValTr[],
3670  double *loLim,
3671  double *hiLim)
3672 {
3673  DEBUG_ENTRY( "SetLimits()" );
3674 
3675  if( optimize.lgVarOn )
3676  {
3677  int ptr0,ptr1;
3678  long index[MDIM], j;
3679  const double SECURE = (1. + 20.*(double)FLT_EPSILON);
3680 
3681  *loLim = +DBL_MAX;
3682  *hiLim = -DBL_MAX;
3683 
3684  switch( grid->imode )
3685  {
3686  case IM_RECT_GRID:
3687  *loLim = -DBL_MAX;
3688  *hiLim = +DBL_MAX;
3689  SetLimitsSub( grid, val, indlo, indhi, index, grid->ndim, loLim, hiLim );
3690  break;
3691  case IM_COSTAR_TEFF_MODID:
3692  case IM_COSTAR_TEFF_LOGG:
3693  case IM_COSTAR_MZAMS_AGE:
3694  for( j=0; j < grid->nTracks; j++ )
3695  {
3696  if( ValTr[j] != -FLT_MAX )
3697  {
3698  /* M_ZAMS is already logarithm, Teff is linear */
3699  double temp = ( grid->imode == IM_COSTAR_MZAMS_AGE ) ?
3700  pow(10.,(double)ValTr[j]) : ValTr[j];
3701  *loLim = MIN2(*loLim,temp);
3702  *hiLim = MAX2(*hiLim,temp);
3703  }
3704  }
3705  break;
3706  case IM_COSTAR_AGE_MZAMS:
3707  index[0] = 0;
3708  index[1] = useTr[0];
3709  ptr0 = grid->jval[JIndex(grid,index)];
3710  index[1] = useTr[1];
3711  ptr1 = grid->jval[JIndex(grid,index)];
3712  *loLim = MAX2(grid->telg[ptr0].par[3],grid->telg[ptr1].par[3]);
3713 # if DEBUGPRT
3714  printf( "set limit 0: (models %d, %d) %f %f\n",
3715  ptr0+1, ptr1+1, grid->telg[ptr0].par[3], grid->telg[ptr1].par[3] );
3716 # endif
3717  index[0] = grid->trackLen[useTr[0]]-1;
3718  index[1] = useTr[0];
3719  ptr0 = grid->jval[JIndex(grid,index)];
3720  index[0] = grid->trackLen[useTr[1]]-1;
3721  index[1] = useTr[1];
3722  ptr1 = grid->jval[JIndex(grid,index)];
3723  *hiLim = MIN2(grid->telg[ptr0].par[3],grid->telg[ptr1].par[3]);
3724 # if DEBUGPRT
3725  printf( "set limit 1: (models %d, %d) %f %f\n",
3726  ptr0+1, ptr1+1, grid->telg[ptr0].par[3], grid->telg[ptr1].par[3] );
3727 # endif
3728  break;
3729  default:
3730  fprintf( ioQQQ, " SetLimits called with insane value for imode: %d.\n", grid->imode );
3731  cdEXIT(EXIT_FAILURE);
3732  }
3733 
3734  ASSERT( fabs(*loLim) < DBL_MAX && fabs(*hiLim) < DBL_MAX );
3735 
3736  /* check sanity of optimization limits */
3737  if( *hiLim <= *loLim )
3738  {
3739  fprintf( ioQQQ, " no room to optimize: lower limit %.4f, upper limit %.4f.\n",
3740  *loLim,*hiLim );
3741  cdEXIT(EXIT_FAILURE);
3742  }
3743 
3744  /* make a bit of room for round-off errors */
3745  *loLim *= SECURE;
3746  *hiLim /= SECURE;
3747 
3748 # if DEBUGPRT
3749  printf("set limits: %g %g\n",*loLim,*hiLim);
3750 # endif
3751  }
3752  else
3753  {
3754  *loLim = 0.;
3755  *hiLim = 0.;
3756  }
3757  return;
3758 }
3759 
3761  double val,
3762  const long indlo[],
3763  const long indhi[],
3764  long index[],
3765  long nd,
3766  double *loLim,
3767  double *hiLim)
3768 {
3769  long n;
3770 
3771  DEBUG_ENTRY( "SetLimitsSub()" );
3772 
3773  --nd;
3774 
3775  if( nd < 2 )
3776  {
3777  double loLoc = +DBL_MAX;
3778  double hiLoc = -DBL_MAX;
3779 
3780  index[1] = 0;
3781  for( index[0]=0; index[0] < grid->nval[0]; ++index[0] )
3782  {
3783  /* grid->val[0][i] is the array of Teff/Age values in the grid, which
3784  * it is sorted in strict monotonically increasing order. This routine
3785  * searches for the largest range [loLoc,hiLoc] in Teff/Age such that
3786  * loLoc <= val <= hiLoc (but loLoc < hiLoc) and at least one model
3787  * exists for each Teff/Age value in this range. This assures that
3788  * interpolation is safe and the optimizer will not trip... */
3789  n = JIndex(grid,index);
3790  if( grid->jhi[n] < 0 )
3791  {
3792  /* there are no models with this value of Teff */
3793  /* this value of Teff should be outside of allowed range */
3794  if( grid->val[0][index[0]] < val )
3795  loLoc = DBL_MAX;
3796  if( grid->val[0][index[0]] > val )
3797  break;
3798  }
3799  else
3800  {
3801  /* there are models with this value of Teff */
3802  /* update range to include this value of Teff */
3803  if( grid->val[0][index[0]] <= val )
3804  {
3805  /* remember lowest legal value of loLoc
3806  * -> only update if previous value was illegal */
3807  if( loLoc == DBL_MAX )
3808  loLoc = grid->val[0][index[0]];
3809  }
3810  if( grid->val[0][index[0]] >= val )
3811  {
3812  /* remember highest legal value of hiLoc
3813  * -> always update */
3814  hiLoc = grid->val[0][index[0]];
3815  }
3816  }
3817  }
3818 
3819  ASSERT( fabs(loLoc) < DBL_MAX && fabs(hiLoc) < DBL_MAX && loLoc < hiLoc );
3820 
3821  *loLim = MAX2(*loLim,loLoc);
3822  *hiLim = MIN2(*hiLim,hiLoc);
3823  }
3824  else
3825  {
3826  index[nd] = indlo[nd];
3827  SetLimitsSub( grid, val, indlo, indhi, index, nd, loLim, hiLim );
3828 
3829  if( indhi[nd] != indlo[nd] )
3830  {
3831  index[nd] = indhi[nd];
3832  SetLimitsSub( grid, val, indlo, indhi, index, nd, loLim, hiLim );
3833  }
3834  }
3835  return;
3836 }
3837 
3839  bool lgList)
3840 {
3841  long i, *index, j, jsize, nd;
3842  double *val;
3843 
3844  DEBUG_ENTRY( "InitIndexArrays()" );
3845 
3846  ASSERT( grid->telg != NULL );
3847  ASSERT( grid->nmods > 0 );
3848 
3849  jsize = 1;
3850 
3851  /* this loop creates a list of all unique model parameter values in increasing order */
3852  for( nd=0; nd < grid->ndim; nd++ )
3853  {
3854  double pval = grid->telg[0].par[nd];
3855  grid->val[nd][0] = pval;
3856  grid->nval[nd] = 1;
3857 
3858  for( i=1; i < grid->nmods; i++ )
3859  {
3860  bool lgOutOfRange;
3861  long i1, i2;
3862 
3863  pval = grid->telg[i].par[nd];
3864  FindIndex( grid->val[nd], grid->nval[nd], pval, &i1, &i2, &lgOutOfRange );
3865  /* if i1 < i2, the new parameter value was not present yet and
3866  * it needs to be inserted in between i1 and i2 --> first move
3867  * all entries from i2 to grid->nval[nd]-1 one slot upward and
3868  * then insert the new value at i2; this also works correctly
3869  * if lgOutOfRange is set, hence no special check is needed */
3870  if( i1 < i2 )
3871  {
3872  /* val[nd] has grid->nmods entries, so cannot overflow */
3873  for( j = grid->nval[nd]-1; j >= i2; j-- )
3874  grid->val[nd][j+1] = grid->val[nd][j];
3875  grid->val[nd][i2] = pval;
3876  grid->nval[nd]++;
3877  }
3878  }
3879 
3880  jsize *= grid->nval[nd];
3881 
3882 # if DEBUGPRT
3883  printf( "%s[%ld]:", grid->names[nd], grid->nval[nd] );
3884  for( i=0; i < grid->nval[nd]; i++ )
3885  printf( " %g", grid->val[nd][i] );
3886  printf( "\n" );
3887 # endif
3888  }
3889 
3890  index = (long *)MALLOC( sizeof(long)*grid->ndim );
3891  val = (double *)MALLOC( sizeof(double)*grid->ndim );
3892 
3893  /* this memory will be freed in the calling function */
3894  grid->jlo = (long *)MALLOC( sizeof(long)*jsize );
3895  grid->jhi = (long *)MALLOC( sizeof(long)*jsize );
3896 
3897  /* set up square array of model indices; this will be used to
3898  * choose the correct models for the interpolation process */
3899  FillJ( grid, index, val, grid->ndim, lgList );
3900 
3901  FREE_CHECK( val );
3902  FREE_CHECK( index );
3903 
3904  if( lgList )
3905  cdEXIT(EXIT_SUCCESS);
3906  return;
3907 }
3908 
3910  long index[], /* index[grid->ndim] */
3911  double val[], /* val[grid->ndim] */
3912  long nd,
3913  bool lgList)
3914 {
3915  DEBUG_ENTRY( "FillJ()" );
3916 
3917  --nd;
3918 
3919  if( nd < 0 )
3920  {
3921  long n = JIndex(grid,index);
3922  SearchModel( grid->telg, grid->nmods, val, grid->ndim, &grid->jlo[n], &grid->jhi[n] );
3923  }
3924  else
3925  {
3926  for( index[nd]=0; index[nd] < grid->nval[nd]; index[nd]++ )
3927  {
3928  val[nd] = grid->val[nd][index[nd]];
3929  FillJ( grid, index, val, nd, lgList );
3930  }
3931  }
3932 
3933  if( lgList && nd == MIN2(grid->ndim-1,1) )
3934  {
3935  long n;
3936  bool lgAge = ( strcmp( grid->names[0], "Age" ) == 0 );
3937 
3938  fprintf( ioQQQ, "\n" );
3939  if( grid->ndim > 2 )
3940  {
3941  fprintf( ioQQQ, "subgrid for" );
3942  for( n = nd+1; n < grid->ndim; n++ )
3943  fprintf( ioQQQ, " %s=%g", grid->names[n], val[n] );
3944  fprintf( ioQQQ, ":\n\n" );
3945  }
3946  if( grid->ndim > 1 )
3947  {
3948  fprintf( ioQQQ, "Teff\\lg g|" );
3949  for( n = 0; n < grid->nval[1]; n++ )
3950  fprintf( ioQQQ, " %5.2f", grid->val[1][n] );
3951  fprintf( ioQQQ, "\n" );
3952  fprintf( ioQQQ, "---------|" );
3953  for( n = 0; n < grid->nval[1]; n++ )
3954  fprintf( ioQQQ, "------" );
3955  }
3956  else
3957  {
3958  if( lgAge )
3959  fprintf( ioQQQ, " Age |\n" );
3960  else
3961  fprintf( ioQQQ, " Teff |\n" );
3962  fprintf( ioQQQ, "---------|------" );
3963  }
3964  fprintf( ioQQQ, "\n" );
3965  for( index[0]=0; index[0] < grid->nval[0]; index[0]++ )
3966  {
3967  if( lgAge )
3968  fprintf( ioQQQ, "%8.2e |", grid->val[0][index[0]] );
3969  else
3970  fprintf( ioQQQ, "%8.0f |", grid->val[0][index[0]] );
3971  if( grid->ndim > 1 )
3972  {
3973  for( index[1]=0; index[1] < grid->nval[1]; index[1]++ )
3974  if( grid->jlo[JIndex(grid,index)] == grid->jhi[JIndex(grid,index)] &&
3975  grid->jlo[JIndex(grid,index)] >= 0 )
3976  fprintf( ioQQQ, " %5ld", grid->jlo[JIndex(grid,index)]+1 );
3977  else
3978  fprintf( ioQQQ, " --" );
3979  }
3980  else
3981  {
3982  fprintf( ioQQQ, " %5ld", grid->jlo[JIndex(grid,index)]+1 );
3983  }
3984  fprintf( ioQQQ, "\n" );
3985  }
3986  }
3987  return;
3988 }
3989 
3991  const long index[]) /* index[grid->ndim] */
3992 {
3993  long i, ind, mul;
3994 
3995  DEBUG_ENTRY( "JIndex()" );
3996 
3997  ind = 0;
3998  mul = 1;
3999  for( i=0; i < grid->ndim; i++ )
4000  {
4001  ind += index[i]*mul;
4002  mul *= grid->nval[i];
4003  }
4004  return ind;
4005 }
4006 
4007 STATIC void SearchModel(const mpp telg[], /* telg[nmods] */
4008  long nmods,
4009  const double val[], /* val[ndim] */
4010  long ndim,
4011  long *index_low,
4012  long *index_high)
4013 {
4014  long i, nd;
4015  double alogg_low = -DBL_MAX, alogg_high = DBL_MAX;
4016 
4017  DEBUG_ENTRY( "SearchModel()" );
4018 
4019  /* given values for the model parameters, this routine searches for
4020  * the atmosphere model that is the best match. If all parameters
4021  * can be matched simultaneously the choice is obvious, but this
4022  * cannot always be achieved (typically for high Teff, the low
4023  * log(g) models will be missing). The rule is that all parameters
4024  * except log(g) must always be matched (such a model is not always
4025  * guaranteed to exist). If all requested parameters can be matched
4026  * exactly, both index_low and index_high will point to that model.
4027  * If all parameters except log(g) can be matched exactly, it will
4028  * return the model with the lowest log(g) value larger than the
4029  * requested value in index_high, and the model with the highest
4030  * log(g) value lower than the requested value in index_low. If
4031  * either requirement cannot be fulfilled, -2 will be returned. */
4032 
4033  *index_low = *index_high = -2;
4034  for( i=0; i < nmods; i++ )
4035  {
4036  bool lgNext = false;
4037  /* ignore models with different parameters */
4038  for( nd=0; nd < ndim; nd++ )
4039  {
4040  if( nd != 1 && fabs(telg[i].par[nd]-val[nd]) > 10.*DBL_EPSILON*fabs(val[nd]) )
4041  {
4042  lgNext = true;
4043  break;
4044  }
4045  }
4046  if( lgNext )
4047  continue;
4048 
4049  /* an exact match is found */
4050  /* this needs to be <= otherwise the test will fail for values equal to 0 */
4051  if( ndim == 1 || fabs(telg[i].par[1]-val[1]) <= 10.*DBL_EPSILON*fabs(val[1]) )
4052  {
4053  *index_low = i;
4054  *index_high = i;
4055  return;
4056  }
4057  /* keep a record of the highest log(g) model smaller than alogg */
4058  if( telg[i].par[1] < val[1] && telg[i].par[1] > alogg_low )
4059  {
4060  *index_low = i;
4061  alogg_low = telg[i].par[1];
4062  }
4063  /* also keep a record of the lowest log(g) model greater than alogg */
4064  if( telg[i].par[1] > val[1] && telg[i].par[1] < alogg_high )
4065  {
4066  *index_high = i;
4067  alogg_high = telg[i].par[1];
4068  }
4069  }
4070  return;
4071 }
4072 
4073 STATIC void FindIndex(const double xval[], /* xval[NVAL] */
4074  long NVAL,
4075  double x,
4076  long *ind1,
4077  long *ind2,
4078  bool *lgInvalid)
4079 {
4080  bool lgOutLo, lgOutHi;
4081  long i;
4082 
4083  DEBUG_ENTRY( "FindIndex()" );
4084 
4085  /* this routine searches for indices ind1, ind2 such that
4086  * xval[ind1] < x < xval[ind2]
4087  * if x is equal to one of the values in xval, then
4088  * ind1 == ind2 and xval[ind1] == x
4089  *
4090  * if x is outside the range xval[0] ... xval[NVAL-1]
4091  * then lgInvalid will be set to true
4092  *
4093  * NB NB -- this routine implicitly assumes that xval is
4094  * strictly monotonically increasing!
4095  */
4096 
4097  ASSERT( NVAL > 0 );
4098 
4099  /* is x outside of range xval[0] ... xval[NVAL-1]? */
4100  lgOutLo = ( x-xval[0] < -10.*DBL_EPSILON*fabs(xval[0]) );
4101  lgOutHi = ( x-xval[NVAL-1] > 10.*DBL_EPSILON*fabs(xval[NVAL-1]) );
4102 
4103  if( lgOutLo || lgOutHi )
4104  {
4105  /* pretend there are two fictitious array elements
4106  * xval[-1] = -Inf and xval[NVAL] = +Inf,
4107  * and return ind1 and ind2 accordingly. This behavior
4108  * is needed for InitIndexArrays() to work correctly */
4109  *ind1 = lgOutLo ? -1 : NVAL-1;
4110  *ind2 = lgOutLo ? 0 : NVAL;
4111  *lgInvalid = true;
4112  return;
4113  }
4114 
4115  *lgInvalid = false;
4116 
4117  /* there are more efficient ways of doing this, e.g. a binary search.
4118  * However, the xval arrays typically only have 1 or 2 dozen elements,
4119  * so the overhead is negligible and the clarity of this code is preferred */
4120 
4121  /* first look for an "exact" match */
4122  for( i=0; i < NVAL; i++ )
4123  {
4124  /* this needs to be <= otherwise the test will fail for x == 0 */
4125  if( fabs(xval[i]-x) <= 10.*DBL_EPSILON*fabs(xval[i]) )
4126  {
4127  *ind1 = i;
4128  *ind2 = i;
4129  return;
4130  }
4131  }
4132 
4133  /* no match was found -> bracket the x value */
4134  for( i=0; i < NVAL-1; i++ )
4135  {
4136  if( xval[i] < x && x < xval[i+1] )
4137  {
4138  *ind1 = i;
4139  *ind2 = i+1;
4140  return;
4141  }
4142  }
4143 
4144  /* this should never be reached ! */
4145  fprintf( ioQQQ, " insanity in FindIndex\n" );
4146  ShowMe();
4147  cdEXIT(EXIT_FAILURE);
4148 }
4149 
4150 STATIC bool lgFileReadable(const char *chFnam, process_counter& pc, access_scheme scheme)
4151 {
4152  DEBUG_ENTRY( "lgFileReadable()" );
4153 
4154  FILE *ioIN;
4155 
4156  ioIN = open_data( chFnam, "r", scheme );
4157  if( ioIN != NULL )
4158  {
4159  fclose( ioIN );
4160  ++pc.nFound;
4161  return true;
4162  }
4163  else
4164  {
4165  return false;
4166  }
4167 }
4168 
4169 /*ValidateGrid: check each model in the grid to see if it has the correct Teff */
4171  double toler)
4172 {
4173  long i, k, nd;
4174  realnum *anu, *flux;
4175 
4176  DEBUG_ENTRY( "ValidateGrid()" );
4177 
4178  if( strcmp( grid->names[0], "Teff" ) != 0 )
4179  {
4180  return;
4181  }
4182 
4183  anu = (realnum *)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
4184  flux = (realnum *)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
4185 
4186  GetModel( grid, -1, anu, lgSILENT, lgLINEAR );
4187 
4188  for( i=0; i < grid->nmods; i++ )
4189  {
4190  fprintf( ioQQQ, "testing model %ld ", i+1 );
4191  for( nd=0; nd < grid->npar; nd++ )
4192  fprintf( ioQQQ, " %s %g", grid->names[nd], grid->telg[i].par[nd] );
4193 
4194  GetModel( grid, i, flux, lgSILENT, lgLINEAR );
4195 
4196  if( lgValidModel( anu, flux, grid->telg[i].par[0], toler ) )
4197  fprintf( ioQQQ, " OK\n" );
4198 
4199  if( false )
4200  {
4201  FILE *ioBUG = fopen( "atmosphere_dump.txt", ( i == 0 ) ? "w" : "a" );
4202 
4203  fprintf( ioBUG, "######## MODEL %ld", i+1 );
4204  for( nd=0; nd < grid->npar; nd++ )
4205  fprintf( ioBUG, " %s %g", grid->names[nd], grid->telg[i].par[nd] );
4206  fprintf( ioBUG, " ####################\n" );
4207 
4208  for( k=0; k < rfield.nupper; ++k )
4209  fprintf( ioBUG, "%e %e\n", anu[k], flux[k] );
4210 
4211  fclose( ioBUG );
4212  }
4213  }
4214 
4215  FREE_CHECK( flux );
4216  FREE_CHECK( anu );
4217  return;
4218 }
4219 
4220 STATIC bool lgValidModel(const realnum anu[],
4221  const realnum flux[],
4222  double Teff,
4223  double toler)
4224 {
4225  bool lgPassed = true;
4226  long k;
4227  double chk, lumi;
4228 
4229  DEBUG_ENTRY( "lgValidModel()" );
4230 
4231  ASSERT( Teff > 0. );
4232 
4233  lumi = 0.;
4234  /* rebinned models are in cgs F_nu units */
4235  for( k=1; k < rfield.nupper; k++ )
4236  lumi += (anu[k] - anu[k-1])*(flux[k] + flux[k-1])/2.;
4237 
4238  /* now convert luminosity to effective temperature */
4239  chk = pow(lumi*FR1RYD/STEFAN_BOLTZ,0.25);
4240  /* the allowed tolerance is set by the caller in toler */
4241  if( fabs(Teff - chk) > toler*Teff ) {
4242  fprintf( ioQQQ, "\n*** WARNING, Teff discrepancy for this model, expected Teff %.2f, ", Teff);
4243  fprintf( ioQQQ, "integration yielded Teff %.2f, delta %.2f%%\n", chk, (chk/Teff-1.)*100. );
4244  lgPassed = false;
4245  }
4246  return lgPassed;
4247 }
4248 
4249 /*RebinAtmosphere: generic routine for rebinning atmospheres onto Cloudy grid */
4250 STATIC void RebinAtmosphere(long nCont, /* the number of points in the incident continuum*/
4251  const realnum StarEner[], /* StarEner[nCont], the freq grid for the model, in Ryd*/
4252  const realnum StarFlux[], /* StarFlux[nCont], the original model flux */
4253  realnum CloudyFlux[], /* CloudyFlux[NC_ELL], the model flux on the cloudy grid */
4254  long nEdge, /* the number of bound-free continuum edges in AbsorbEdge */
4255  const realnum AbsorbEdge[]) /* AbsorbEdge[nEdge], energies of the edges */
4256 {
4257  bool lgDone;
4258  long int ind,
4259  j,
4260  k;
4261  /* >>chng 00 jun 02, demoted next two to realnum, PvH */
4262  realnum BinHigh,
4263  BinLow,
4264  BinMid,
4265  BinNext,
4266  *EdgeLow=NULL,
4267  *EdgeHigh=NULL,
4268  *StarPower;
4269 
4270  DEBUG_ENTRY( "RebinAtmosphere()" );
4271 
4272  if( nEdge > 0 )
4273  {
4274  EdgeLow = (realnum*)MALLOC( sizeof(realnum)*(unsigned)nEdge );
4275  EdgeHigh = (realnum*)MALLOC( sizeof(realnum)*(unsigned)nEdge );
4276  }
4277 
4278  /* this loop should be before the next loop, otherwise models with a
4279  * very strong He II edge (i.e. no flux beyond that edge) will fail */
4280  for( j=0; j < nEdge; j++ )
4281  {
4282  ind = RebinFind(StarEner,nCont,AbsorbEdge[j]);
4283 
4284  /* sanity check */
4285  ASSERT( ind >= 0 && ind+1 < nCont );
4286 
4287  EdgeLow[j] = StarEner[ind];
4288  EdgeHigh[j] = StarEner[ind+1];
4289  }
4290 
4291  /* cut off that part of the Wien tail that evaluated to zero */
4292  /* >> chng 05 nov 22, inverted loop, slightly faster PvH */
4293  /*for( j=nCont-1; j >= 0; j-- )*/
4294  for( j=0; j < nCont; j++ )
4295  {
4296  if( StarFlux[j] == 0.f )
4297  {
4298  nCont = j;
4299  break;
4300  }
4301  }
4302  ASSERT( nCont > 0 );
4303 
4304  StarPower = (realnum *)MALLOC( sizeof(realnum)*(unsigned)(nCont-1) );
4305 
4306  for( j=0; j < nCont-1; j++ )
4307  {
4308  double ratio_x, ratio_y;
4309 
4310  /* >>chng 05 nov 22, add sanity check to prevent invalid fp operations */
4311  ASSERT( StarEner[j+1] > StarEner[j] );
4312 
4313  /* >>chng 06 aug 11, on some systems (e.g., macbook pro) y/x can get evaluated as y*(1/x);
4314  * this causes overflows if x is a denormalized number, hence we force a cast to double, PvH */
4315  ratio_x = (double)StarEner[j+1]/(double)StarEner[j];
4316  ratio_y = (double)StarFlux[j+1]/(double)StarFlux[j];
4317  StarPower[j] = (realnum)(log(ratio_y)/log(ratio_x));
4318  }
4319 
4320  for( j=0; j < rfield.nupper; j++ )
4321  {
4322  /* >>chng 05 nov 22, modified BinLow, BinHigh, BinNext to make boundaries match exactly, PvH */
4323  /* BinLow is lower bound of this continuum cell */
4324  BinLow = ( j > 0 ) ?
4325  (realnum)sqrt(rfield.anu[j-1]*rfield.anu[j]) : (realnum)sqrt(POW3(rfield.anu[0])/rfield.anu[1]);
4326 
4327  /* BinHigh is upper bound of this continuum cell */
4328  BinHigh = ( j+1 < rfield.nupper ) ?
4329  (realnum)sqrt(rfield.anu[j]*rfield.anu[j+1]) : rfield.anu[rfield.nupper-1];
4330 
4331  /* BinNext is upper bound of next continuum cell */
4332  BinNext = ( j+2 < rfield.nupper ) ?
4333  (realnum)sqrt(rfield.anu[j+1]*rfield.anu[j+2]) : rfield.anu[rfield.nupper-1];
4334 
4335  lgDone = false;
4336 
4337  /* >>chng 00 aug 14, take special care not to interpolate over major edges,
4338  * the region in between EdgeLow and EdgeHigh should be avoided,
4339  * the spectrum is extremely steep there, leading to significant roundoff error, PvH */
4340  for( k=0; k < nEdge; k++ )
4341  {
4342  if( BinLow < EdgeLow[k] && BinNext > EdgeHigh[k] )
4343  {
4344  BinMid = 0.99999f*EdgeLow[k];
4345  CloudyFlux[j] = RebinSingleCell(BinLow,BinMid,StarEner,StarFlux,StarPower,nCont);
4346  j++;
4347 
4348  /* sanity check */
4349  ASSERT( j < rfield.nupper );
4350 
4351  BinMid = 1.00001f*EdgeHigh[k];
4352  CloudyFlux[j] = RebinSingleCell(BinMid,BinNext,StarEner,StarFlux,StarPower,nCont);
4353  lgDone = true;
4354  break;
4355  }
4356  }
4357 
4358  /* default case when we are not close to an edge */
4359  if( !lgDone )
4360  {
4361  CloudyFlux[j] = RebinSingleCell(BinLow,BinHigh,StarEner,StarFlux,StarPower,nCont);
4362  }
4363  }
4364 
4365  FREE_CHECK( StarPower );
4366  FREE_SAFE( EdgeHigh );
4367  FREE_SAFE( EdgeLow );
4368  return;
4369 }
4370 
4372  realnum BinHigh,
4373  const realnum StarEner[], /* StarEner[nCont] */
4374  const realnum StarFlux[], /* StarFlux[nCont] */
4375  const realnum StarPower[], /* StarPower[nCont-1] */
4376  long nCont)
4377 {
4378  long int i,
4379  ipHi,
4380  ipLo;
4381  double anu,
4382  retval,
4383  widflx;
4384  double sum,
4385  v1,
4386  val,
4387  x1,
4388  x2;
4389 
4390  DEBUG_ENTRY( "RebinSingleCell()" );
4391 
4392  /* >>chng 05 nov 22, use geometric mean instead of arithmetic mean, PvH */
4393  anu = sqrt(BinLow*BinHigh);
4394  /* >>chng 05 nov 22, reduce widflx if cell sticks out above highest frequency in model, PvH */
4395  widflx = MIN2(BinHigh,StarEner[nCont-1])-BinLow;
4396 
4397  if( BinLow < StarEner[0] )
4398  {
4399  /* this is case where Cloudy's continuum is below stellar continuum,
4400  * (at least for part of the cell), so we do Rayleigh Jeans extrapolation */
4401  retval = (realnum)(StarFlux[0]*pow(anu/StarEner[0],2.));
4402  }
4403  else if( BinLow > StarEner[nCont-1] )
4404  {
4405  /* case where cloudy continuum is entirely above highest stellar point */
4406  retval = 0.0e00;
4407  }
4408  else
4409  {
4410  /* now go through stellar continuum to find bins corresponding to
4411  * this cloudy cell, stellar continuum defined through nCont cells */
4412  ipLo = RebinFind(StarEner,nCont,BinLow);
4413  ipHi = RebinFind(StarEner,nCont,BinHigh);
4414  /* sanity check */
4415  ASSERT( ipLo >= 0 && ipLo < nCont-1 && ipHi >= ipLo );
4416 
4417  if( ipHi == ipLo )
4418  {
4419  /* Do the case where the cloudy cell and its edges are between
4420  * two adjacent stellar model points: do power-law interpolation */
4421  retval = (realnum)(StarFlux[ipLo]*pow(anu/StarEner[ipLo],(double)StarPower[ipLo]));
4422  }
4423  else
4424  {
4425  /* Do the case where the cloudy cell and its edges span two or more
4426  * stellar model cells: add segments with power-law interpolation up to
4427  * do the averaging.*/
4428 
4429  sum = 0.;
4430 
4431  /* ipHi points to stellar point at high end of cloudy continuum cell,
4432  * if the Cloudy cell extends beyond the stellar grid, ipHi == nCont-1
4433  * and the MIN2(ipHi,nCont-2) prevents access beyond allocated memory
4434  * ipLo points to low end, above we asserted that 0 <= ipLo < nCont-1 */
4435  for( i=ipLo; i <= MIN2(ipHi,nCont-2); i++ )
4436  {
4437  double pp1 = StarPower[i] + 1.;
4438 
4439  if( i == ipLo )
4440  {
4441  x1 = BinLow;
4442  x2 = StarEner[i+1];
4443  v1 = StarFlux[i]*pow(x1/StarEner[i],(double)StarPower[i]);
4444  /*v2 = StarFlux[i+1];*/
4445  }
4446 
4447  else if( i == ipHi )
4448  {
4449  x2 = BinHigh;
4450  x1 = StarEner[i];
4451  /*v2 = StarFlux[i]*pow(x2/StarEner[i],StarPower[i]);*/
4452  v1 = StarFlux[i];
4453  }
4454 
4455  /*if( i > ipLo && i < ipHi )*/
4456  else
4457  {
4458  x1 = StarEner[i];
4459  x2 = StarEner[i+1];
4460  v1 = StarFlux[i];
4461  /*v2 = StarFlux[i+1];*/
4462  }
4463 
4464  if( fabs(pp1) < 0.001 )
4465  {
4466  val = x1*v1*log(x2/x1);
4467  }
4468  else
4469  {
4470  val = pow(x2/x1,pp1) - 1.;
4471  val = val*x1*v1/pp1;
4472  }
4473  sum += val;
4474  }
4475 
4476  retval = sum/widflx;
4477  }
4478  }
4479  return (realnum)retval;
4480 }
4481 
4482 STATIC long RebinFind(const realnum array[], /* array[nArr] */
4483  long nArr,
4484  realnum val)
4485 {
4486  long i1,
4487  i2,
4488  i3,
4489  ind = -2,
4490  sgn;
4491 
4492  DEBUG_ENTRY( "RebinFind()" );
4493 
4494  /* sanity check */
4495  ASSERT( nArr > 1 );
4496 
4497  /* return ind(val) such that array[ind] <= val <= array[ind+1],
4498  *
4499  * NB NB: this routine assumes that array[] increases monotonically !
4500  *
4501  * the first two clauses indicate out-of-bounds conditions and
4502  * guarantee that when val1 <= val2, also ind(val1) <= ind(val2) */
4503 
4504  if( val < array[0] )
4505  {
4506  ind = -1;
4507  }
4508  else if( val > array[nArr-1] )
4509  {
4510  ind = nArr-1;
4511  }
4512  else
4513  {
4514  /* do a binary search for ind */
4515  i1 = 0;
4516  i3 = nArr-1;
4517  while( i3-i1 > 1 )
4518  {
4519  i2 = (i1+i3)/2;
4520  sgn = sign3(val-array[i2]);
4521 
4522  switch(sgn)
4523  {
4524  case -1:
4525  i3 = i2;
4526  break;
4527  case 0:
4528  ind = i2;
4529  return( ind );
4530  case 1:
4531  i1 = i2;
4532  break;
4533  }
4534  }
4535  ind = i1;
4536  }
4537 
4538  /* sanity check */
4539  ASSERT( ind > -2 );
4540  return ind;
4541 }
4542 /*lint +e785 too few initializers */
4543 /*lint +e801 use of go to depreciated */

Generated for cloudy by doxygen 1.8.3.1