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