cloudy trunk
Loading...
Searching...
No Matches
parse_table.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/*ParseTable parse the table read command */
4/*lines_table invoked by table lines command, check if we can find all lines in a given list */
5/*read_hm05 read in the data files and interpolate the continuum to
6 * the correct redshift */
7#include "cddefines.h"
8#include "cddrive.h"
9#include "physconst.h"
10#include "optimize.h"
11#include "rfield.h"
12#include "trace.h"
13#include "radius.h"
14#include "input.h"
15#include "stars.h"
16#include "lines.h"
17#include "prt.h"
18#include "parser.h"
19#include "save.h"
20#include "thirdparty.h"
21#include "continuum.h"
22/* HP cc cannot compile this routine with any optimization */
23#if defined(__HP_aCC)
24# pragma OPT_LEVEL 1
25#endif
26
27/*ReadTable called by TABLE READ to read in continuum from PUNCH TRANSMITTED CONTINUUM */
28STATIC void ReadTable(const char * fnam);
29
30static string chLINE_LIST;
31
32/* tables of various built in continue */
33/* Davidson 1985 version of crab nebula SED */
34static const int NCRAB = 10;
35static double tnucrb[NCRAB],
37
38/* Bob Rubin's corrected theta 1 Ori C continuum */
39static const int NRUBIN = 56;
40static double tnurbn[NRUBIN] = {1.05E-08,1.05E-07,1.05E-06,1.04E-05,1.00E-04,1.00E-03,1.00E-02,3.01E-02,1.00E-01,
41 1.50E-01,2.50E-01,4.01E-01,6.01E-01,9.8E-01,9.96E-01,1.00E+00,1.02445,1.07266,1.12563,1.18411,1.23881,
42 1.29328,1.35881,1.42463,1.48981,1.55326,1.6166,1.68845,1.76698,1.8019,1.808,1.84567,1.9317,2.04891,2.14533,
43 2.19702,2.27941,2.37438,2.43137,2.49798,2.56113,2.59762,2.66597,2.80543,2.95069,3.02911,3.11182,3.22178,
44 3.3155,3.42768,3.50678,3.56157,3.61811,3.75211,3.89643,4.},
45 fnurbn[NRUBIN] = {1.56E-20,1.55E-17,1.54E-14,1.53E-11,1.35E-08,1.34E-05,1.35E-02,3.62E-01,1.27E+01,
46 3.90E+01,1.48E+02,4.52E+02,1.02E+03,2.27E+03,8.69E+02,8.04E+02,6.58E+02,6.13E+02,6.49E+02,6.06E+02,
47 6.30E+02,5.53E+02,5.55E+02,5.24E+02,4.86E+02,4.49E+02,4.42E+02,3.82E+02,3.91E+02,2.91E+02,2.61E+02,
48 1.32E+02,1.20E+02,1.16E+02,9.56E+01,9.94E+01,9.10E+01,4.85E+01,7.53E+01,9.53E+01,8.52E+01,5.76E+01,
49 6.72E+01,5.20E+01,8.09E+00,3.75E+00,5.57E+00,3.80E+00,2.73E+00,2.22E+00,3.16E-01,1.26E-01,1.39E-01,
50 6.15E-02,3.22E-02,7.98E-03};
51
52/* stores numbers for table cooling flow */
53static const int NCFL = 40;
54static double cfle[NCFL],
56
57/* Brad Peterson's AKN 120 continuum */
58static const int NKN120 = 11;
59static double tnuakn[NKN120],
61
62/* Black's ISM continuum, with He hole filled in */
63static const int NISM = 23;
64static double tnuism[NISM],
66
67/* z=2 background,
68 * >>refer continuum background Haardt, Francesco, & Madau, Piero, 1996,
69 * >>refercon ApJ, 461, 20 */
70static const int NHM96 = 14;
71/* log energy in Ryd */
72static const double tnuHM96[NHM96]={-8,-1.722735683,-0.351545683,-0.222905683,-0.133385683,
73/* changeg these two energies to prevent degeneracy */
74-0.127655683,-0.004575683,0.297544317,0.476753,0.476756,0.588704317,
750.661374317,1.500814317,2.245164317};
76/*-0.127655683,-0.004575683,0.297544317,0.476754317,0.476754317,0.588704317,*/
77/*log J in the units of (erg cm^{-2} s^{-1} Hz^{-1} sr^{-1})*/
78static const double fnuHM96[NHM96]={-32.53342863,-19.9789,-20.4204,-20.4443,-20.5756,-20.7546,
79-21.2796,-21.6256,-21.8404,-21.4823,-22.2102,-22.9263,-23.32,-24.2865};
80
81/* Mathews and Ferland generic AGN continuum */
82static const int NAGN = 8;
83static double tnuagn[NAGN],
85
86/* table Draine ISM continuum */
87static const int NDRAINE = 15;
88static double tnudrn[NDRAINE] , tsldrn[NDRAINE];
89
90/* routine that stores values for above vectors */
91STATIC void ZeroContin(void);
92
93/* this allows the low energy point of any built in array to be reset to the
94 * current low energy point in the code - nothing need be done if this is reset
95 * tnu is array of energies, [0] is first, and we want it to be lower
96 * fluxlog is flux at tnu, and may or may not be log
97 * lgLog says whether it is */
98STATIC void resetBltin( double *tnu , double *fluxlog , bool lgLog )
99{
100 /* this will multiply low-energy bounds of code and go into element[0]
101 * ensures that energy range is fully covered */
102 const double RESETFACTOR = 0.98;
103 double power;
104 /* this makes sure we are called after emm is defined */
105 ASSERT( rfield.emm > 0. );
106
107 if( lgLog )
108 {
109 /* continuum comes in as log of flux */
110 /* this is current power-law slope of low-energy continuum */
111 power = (fluxlog[1] - fluxlog[0] ) / log10( tnu[1]/tnu[0] );
112 /* this will be new low energy bounds to this continuum */
113 tnu[0] = rfield.emm*RESETFACTOR;
114 fluxlog[0] = fluxlog[1] + power * log10( tnu[0] /tnu[1] );
115 }
116 else
117 {
118 /* continuum comes in as linear flux */
119 /* this is current power-law slope of low-energy continuum */
120 power = log10( fluxlog[1]/fluxlog[0]) / log10( tnu[1]/tnu[0] );
121 /* this will be new low energy bounds to this continuum */
122 tnu[0] = rfield.emm*RESETFACTOR;
123 fluxlog[0] = log10(fluxlog[1]) + power * log10( tnu[0] /tnu[1] );
124 /* flux is not really log, we want linear */
125 fluxlog[0] = pow(10. , fluxlog[0]);
126 }
127 /*fprintf(ioQQQ," power is %f lgLog is %i\n", power, lgLog );*/
128 return;
129}
130
131/*read_hm05 read in the data files and interpolate the Haardt & Madau continuum to
132 * the correct redshift */
133STATIC void read_hm05( const char chFile[] , double **tnuHM , double **fnuHM ,
134 long int *nhm , double redshift )
135{
136
137 FILE *ioFILE;
138 double *table_redshifts = NULL,
139 *table_wl = NULL ,
140 **table_fn=NULL,
141 frac_hi;
142 char chLine[1000];
143 long int nRedshift , i , n , nz , ipLo , ipHi;
144 bool lgEOL;
145
146 DEBUG_ENTRY( "read_hm05()" );
147
148 ioFILE = open_data( chFile, "r", AS_LOCAL_DATA );
149
150 /* this will be the number of continuum points in their table */
151 *nhm = 0;
152 /* the number of redshifts in their table */
153 nRedshift = -1;
154 /* first read past comments and get magic number */
155 while( read_whole_line( chLine , (int)sizeof(chLine) , ioFILE ) != NULL )
156 {
157 /* we want to count the lines that do not start with #
158 * since these contain data */
159 if( chLine[0] != '#')
160 {
161 ++*nhm;
162 /* check magic number if first valid line */
163 if( *nhm==1 )
164 {
165 long mag_read;
166 /* this is the magic number - read it in */
167 i = 1;
168 mag_read = (long)FFmtRead(chLine,&i,(int)sizeof(chLine),&lgEOL);
169 if( mag_read != 50808 )
170 {
171 fprintf(ioQQQ,
172 " Magic number in Haardt & Madau file is not correct, I expected 50808 and found %li\n",
173 mag_read );
175 }
176 }
177 /* second valid line count number of redshifts */
178 else if( *nhm==2 )
179 {
180 double a;
181 i = 2;
182 lgEOL = false;
183 nRedshift = 0;
184 while( !lgEOL )
185 {
186 ++nRedshift;
187 a = FFmtRead(chLine,&i,(int)sizeof(chLine),&lgEOL);
188 ASSERT( a >= 0. );
189 }
190 /* have over counted by one - back up one */
191 --nRedshift;
192 /* highest redshift data missing in file i received */
193 --nRedshift;
194 /*fprintf(ioQQQ," number of z is %li\n", nRedshift);*/
195 /* malloc some space */
196 table_redshifts = (double*)MALLOC( (size_t)nRedshift*sizeof(double) );
197 /* now read in the redshifts */
198 i = 2;
199 for( n=0; n<nRedshift; ++n )
200 {
201 table_redshifts[n] = FFmtRead(chLine,&i,(int)sizeof(chLine),&lgEOL);
202 /*fprintf(ioQQQ,"DEBUG %li z %.3f\n", n , table_redshifts[n] );*/
203 }
204 if( lgEOL )
206 }
207 }
208 }
209
210 /* malloc space for the arrays first wavelength array */
211 table_wl = (double*)MALLOC( (size_t)*nhm*sizeof(double) );
212 /* the spectrum array table_fn[z][nu] */
213 table_fn = (double**)MALLOC( (size_t)nRedshift*sizeof(double*) );
214 for(n=0; n<nRedshift; ++n )
215 {
216 table_fn[n] = (double*)MALLOC( (size_t)(*nhm)*sizeof(double) );
217 }
218
219 /* rewind the file so we can read it a second time*/
220 if( fseek( ioFILE , 0 , SEEK_SET ) != 0 )
221 {
222 fprintf( ioQQQ, " read_hm05 could not rewind HM05 date file.\n");
224 }
225 n = 0;
226 *nhm = 0;
227 while( read_whole_line( chLine , (int)sizeof(chLine) , ioFILE ) != NULL )
228 {
229 /* we want to count the lines that do not start with #
230 * since these contain data */
231 if( chLine[0] != '#')
232 {
233 /* this is number of non-comment lines - will skip magic number
234 * and line giving redshift */
235 ++n;
236 if( n>2 )
237 {
238 i = 1;
239 table_wl[*nhm] = FFmtRead(chLine,&i,(int)sizeof(chLine),&lgEOL);
240 if( lgEOL )
242 for( nz=0; nz<nRedshift; ++nz )
243 {
244 /* >>chng 07 feb 18, PvH change from last branck to first */
245 table_fn[nz][*nhm] = FFmtRead(chLine,&i,(int)sizeof(chLine),&lgEOL);
246 }
247 ++*nhm;
248 }
249 }
250 }
251
252 fclose( ioFILE );
253
254 {
255 /* change following to true to print their original table */
256 enum {DEBUG_LOC=false};
257 if( DEBUG_LOC)
258 {
259 fprintf(ioQQQ,"wavelength/z");
260 for(nz=0; nz<nRedshift; ++nz )
261 fprintf(ioQQQ,"\t%.3f", table_redshifts[nz] );
262 fprintf(ioQQQ,"\n");
263 for( i=0; i<*nhm; ++i )
264 {
265 fprintf(ioQQQ,"%.3e", table_wl[i] );
266 for( nz=0; nz<nRedshift; ++nz )
267 fprintf(ioQQQ,"\t%.3e", table_fn[nz][i] );
268 fprintf(ioQQQ,"\n");
269 }
270 }
271 }
272
273 /* this is just to shut up the lint and must not be ASSERT */
274 assert( table_redshifts!=NULL );
275
276 /* first check that desired redshift is within range of table */
277 if( redshift < table_redshifts[0] ||
278 redshift > table_redshifts[nRedshift-1] )
279 {
280 fprintf(ioQQQ," The redshift requested on table HM05 is %g but is not within the range of the table, which goes from %g to %g.\n",
281 redshift, table_redshifts[0] , table_redshifts[nRedshift-1] );
282 fprintf(ioQQQ," Sorry.\n");
284 }
285
286 ipLo = -1;
287 ipHi = -1;
288 /* find which redshift bin we need */
289 for( nz=0; nz<nRedshift-1; ++nz )
290 {
291 if( redshift >= table_redshifts[nz] &&
292 redshift <= table_redshifts[nz+1] )
293 {
294 ipLo = nz;
295 ipHi = nz+1;
296 break;
297 }
298 }
299 ASSERT( ipLo>=0 && ipHi >=0 );
300
301 /* make sure that the wavelengths are in increasing order - they were not in the
302 * original data table, but had repeated values near important ionization edges */
303 for( i=0; i<*nhm-1; ++i )
304 {
305 if( table_wl[i]>=table_wl[i+1] )
306 {
307 fprintf(ioQQQ," The wavelengths in the table HM05 data table are not in increasing order. This is required.\n");
308 fprintf(ioQQQ," Sorry.\n");
310 }
311 }
312
313 /* malloc the space for the returned arrays */
314 *fnuHM = (double *)MALLOC( (size_t)(*nhm)*sizeof(double ) );
315 *tnuHM = (double *)MALLOC( (size_t)(*nhm)*sizeof(double ) );
316
317 /* now fill in the arrays with the interpolated table
318 * we will interpolate linearly in redshift
319 * in general the redshift will be between the tabulated redshift */
320 frac_hi = (redshift-table_redshifts[ipLo]) / (table_redshifts[ipHi]-table_redshifts[ipLo]);
321 for( i=0; i<*nhm; ++i )
322 {
323 /* we have wavelengths in angstroms and want log Ryd
324 * original table was in decreasing energy order, must also
325 * flip it around since need increasing energy order */
326 (*tnuHM)[*nhm-1-i] = RYDLAM / table_wl[i];
327 /* original table in correct units so no conversion needed */
328 (*fnuHM)[*nhm-1-i] = table_fn[ipLo][i]*(1.-frac_hi) + table_fn[ipHi][i]*frac_hi;
329 /*fprintf(ioQQQ,"DEBUG1\t%.3e\t%.3e\n",(*tnuHM)[*nhm-1-i] , (*fnuHM)[*nhm-1-i] );*/
330 }
331
332 for( n=0; n<nRedshift; ++n )
333 free( table_fn[n] );
334 free( table_fn );
335 free( table_wl );
336 free( table_redshifts );
337 return;
338}
339
341{
342 char chFile[FILENAME_PATH_LENGTH_2]; /*file name for table read */
343
344 IntMode imode = IM_ILLEGAL_MODE;
345 bool lgHit,
346 lgLogSet;
347 static bool lgCalled=false;
348
349 long int i,
350 j,
351 nstar,
352 ipNORM;
353
354 double alpha,
355 brakmm,
356 brakxr,
357 ConBreak,
358 fac,
359 scale,
360 slopir,
361 slopxr;
362
363 bool lgNoContinuum = false,
364 lgQuoteFound;
365
366 DEBUG_ENTRY( "ParseTable()" );
367
368 /* if first call then set up values for table */
369 if( !lgCalled )
370 {
371 ZeroContin();
372 lgCalled = true;
373 }
374
375 if( rfield.nShape >= LIMSPC )
376 {
377 fprintf( ioQQQ, " %ld is too many spectra entered. Increase LIMSPC\n Sorry.\n",
378 rfield.nShape );
380 }
381
382 /* three commands, tables line, read, and star, have quotes on the
383 * lines giving file names. must get quotes first so that filename
384 * does not confuse parser */
385 lgQuoteFound = true;
386 if( p.GetQuote( chFile , false ) )
387 lgQuoteFound = false;
388
389 /* set flag telling interpolate */
390 strcpy( rfield.chSpType[rfield.nShape], "INTER" );
391 /* >>chng 06 jul 16, fix memory leak when optimizing, PvH */
392 if( !rfield.lgContMalloc[rfield.nShape] )
393 {
394 rfield.tNu[rfield.nShape].resize( NCELL );
395 rfield.tslop[rfield.nShape].resize( NCELL );
396 rfield.tFluxLog[rfield.nShape].resize( NCELL );
397 rfield.ncont[rfield.nShape] = 0;
398 rfield.lgContMalloc[rfield.nShape] = true;
399 }
400
401 /* NB when adding more keys also change the comment at the end */
402 if( p.nMatch(" AGN") )
403 {
404 /* do Mathews and Ferland generic AGN continuum */
405 ASSERT( NAGN < NCELL);
406 for( i=0; i < NAGN; i++ )
407 {
408 rfield.tNu[rfield.nShape][i].set(tnuagn[i]);
409 rfield.tslop[rfield.nShape][i] = (realnum)tslagn[i];
410 }
411 rfield.ncont[rfield.nShape] = NAGN;
412
413 /* optional keyword break, to adjust IR cutoff */
414 if( p.nMatch("BREA") )
415 {
416 ConBreak = p.FFmtRead();
417
418 if( p.lgEOL() )
419 {
420 /* no break, set to low energy limit of code */
421 if( p.nMatch(" NO ") )
422 {
423 ConBreak = rfield.emm*1.01f;
424 }
425 else
426 {
427 fprintf( ioQQQ, " There must be a number for the break.\n Sorry.\n" );
429 }
430 }
431
432 if( ConBreak == 0. )
433 {
434 fprintf( ioQQQ, " The break must be greater than 0.2 Ryd.\n Sorry.\n" );
436 }
437
438 if( p.nMatch("MICR") )
439 {
440 /* optional keyword, ``microns'', convert to Rydbergs */
441 ConBreak = 0.0912/ConBreak;
442 }
443
444 if( ConBreak < 0. )
445 {
446 /* option to enter break as LOG10 */
447 ConBreak = pow(10.,ConBreak);
448 }
449
450 else if( ConBreak == 0. )
451 {
452 fprintf( ioQQQ, " An energy of 0 is not allowed.\n Sorry.\n" );
454 }
455
456 if( ConBreak >= rfield.tNu[rfield.nShape][2].Ryd() )
457 {
458 fprintf( ioQQQ, " The energy of the break cannot be greater than%10.2e Ryd.\n Sorry.\n",
459 rfield.tNu[rfield.nShape][2].Ryd() );
461 }
462
463 else if( ConBreak <= rfield.tNu[rfield.nShape][0].Ryd() )
464 {
465 fprintf( ioQQQ, " The energy of the break cannot be less than%10.2e Ryd.\n Sorry.\n",
466 rfield.tNu[rfield.nShape][0].Ryd() );
468 }
469
470 rfield.tNu[rfield.nShape][1].set(ConBreak);
471
472 rfield.tslop[rfield.nShape][1] =
473 (realnum)(rfield.tslop[rfield.nShape][2] +
474 log10(rfield.tNu[rfield.nShape][2].Ryd()/rfield.tNu[rfield.nShape][1].Ryd()));
475
476 rfield.tslop[rfield.nShape][0] =
477 (realnum)(rfield.tslop[rfield.nShape][1] -
478 2.5*log10(rfield.tNu[rfield.nShape][1].Ryd()/rfield.tNu[rfield.nShape][0].Ryd()));
479 }
480 }
481
482 else if( p.nMatchErase("AKN1") )
483 {
484 /* AKN120 continuum from Brad Peterson */
485 ASSERT( NKN120 < NCELL );
486 for( i=0; i < NKN120; i++ )
487 {
488 rfield.tNu[rfield.nShape][i].set(tnuakn[i]);
489 rfield.tslop[rfield.nShape][i] = (realnum)log10(fnuakn[i]);
490 }
491 rfield.ncont[rfield.nShape] = NKN120;
492 }
493
494 else if( p.nMatch("CRAB") )
495 {
496 if( p.nMatch("DAVIDSON") )
497 {
498 ASSERT( NCRAB < NCELL );
499 /* Crab Nebula SED from Davidson & Fesen 1985 Ann Rev Ast Ap */
500 for( i=0; i < NCRAB; i++ )
501 {
502 rfield.tNu[rfield.nShape][i].set(tnucrb[i]);
503 rfield.tslop[rfield.nShape][i] = (realnum)log10(fnucrb[i]);
504 }
505 rfield.ncont[rfield.nShape] = NCRAB;
506 }
507 else
508 {
509 /* Crab Nebula SED from
510 * *>>refer continuum CrabNebula Hester, J, 2008, ARA&A, 46, 127-155
511 * log Hz, log nuLnu, digitized from figure*/
512 const int NCRAB08 = 14;
513 ASSERT( NCRAB08 < NCELL );
514 static const double tnucrbHz08[NCRAB08] = {
515 10.0114,
516 12.0499,
517 13.8657,
518 14.8718,
519 15.6126,
520 16.3815,
521 18.0873,
522 18.6467,
523 19.6535,
524 20.9825,
525 21.8082,
526 22.5783,
527 23.293,
528 23.644 };
529 static const double crbnuLnu08[NCRAB08] = {
530 34.4381,
531 35.855,
532 36.7703,
533 37.0142,
534 37.0441,
535 37.0297,
536 36.8609,
537 36.7579,
538 36.5962,
539 36.0585,
540 35.5279,
541 34.8277,
542 33.862,
543 33.0141 };
544 static double tnu[NCRAB08] , tflux[NCRAB08];
545 for( i=0; i < NCRAB08; i++ )
546 {
547 // Hester plot shows log Hz - convert to linear Rydbergs
548 tnu[i] = pow(10. , tnucrbHz08[i] ) / FR1RYD;
549 // plot shows log nu L_nu - convert to linear L_nu
550 tflux[i] = crbnuLnu08[i] - tnucrbHz08[i];
551 }
552 resetBltin( tnu , tflux , false );
553 for( i=0; i < NCRAB08; i++ )
554 {
555 rfield.tNu[rfield.nShape][i].set(tnu[i]);
556 // plot shows log nu L_nu - convert to linear L_nu
557 rfield.tslop[rfield.nShape][i] = (realnum)tflux[i];
558 }
559 rfield.ncont[rfield.nShape] = NCRAB08;
560 }
561 }
562
563 else if( p.nMatch("COOL") )
564 {
565 ASSERT( NCFL < NCELL );
566 /* cooling flow from Johnstone et al. 1992, MN 255, 431. */
567 for( i=0; i < NCFL; i++ )
568 {
569 rfield.tNu[rfield.nShape][i].set(cfle[i]);
570 rfield.tslop[rfield.nShape][i] = (realnum)cflf[i];
571 }
572 rfield.ncont[rfield.nShape] = NCFL;
573 }
574
575 else if( p.nMatchErase("HM96") )
576 {
577 /* this is the old Haardt & Madau continuum, one set of points
578 * with only the quasars
579 * this command does not include the CMB - do that separately with the CMB command */
580 /* set flag telling interpolate */
581 strcpy( rfield.chSpType[rfield.nShape], "INTER" );
582
583 ASSERT( NHM96 < NCELL );
584 /* z=2 background,
585 * >>refer continuum background Haardt, Francesco, & Madau, Piero, 1996, ApJ, 461, 20 */
586 for( j=0; j < NHM96; j++ )
587 {
588 /* frequency was stored as log of ryd */
589 rfield.tNu[rfield.nShape][j].set(pow( 10. , tnuHM96[j] ));
590 rfield.tslop[rfield.nShape][j] = (realnum)fnuHM96[j];
591 }
592 rfield.ncont[rfield.nShape] = NHM96;
593
594 /* optional scale factor to change default intensity from their value
595 * assumed to be log if negative, and linear otherwise */
596 scale = p.FFmtRead();
597 if( scale > 0. )
598 scale = log10(scale);
599
600 /* this also sets continuum intensity*/
601 if( p.m_nqh >= LIMSPC )
602 {
603 fprintf( ioQQQ, " %ld is too many continua entered. Increase LIMSPC\n Sorry.\n",
604 p.m_nqh );
606 }
607
608 /* check that stack of shape and luminosity specifications
609 * is parallel, stop if not - this happens is background comes
610 * BETWEEN another set of shape and luminosity commands */
611 if( rfield.nShape != p.m_nqh )
612 {
613 fprintf( ioQQQ, " This command has come between a previous ordered pair of continuum shape and luminosity commands.\n Reorder the commands to complete each continuum specification before starting another.\n" );
614 fprintf( ioQQQ, " Sorry.\n" );
616 }
617
618 strcpy( rfield.chRSpec[p.m_nqh], "SQCM" );
619 strcpy( rfield.chSpNorm[p.m_nqh], "FLUX" );
620
621 /* this is an isotropic radiation field */
622 rfield.lgBeamed[p.m_nqh] = false;
623 rfield.Illumination[p.m_nqh] = Illuminate::SYMMETRIC;
624
625 /* this will be flux density at some frequency on the table. the numbers
626 * are per Hz and sr so must multiply by 4 pi
627 * [2] is not special, could have been any within array*/
628 rfield.range[p.m_nqh][0] = pow(10. , tnuHM96[2] )*1.0001;
629
630 /* convert intensity HM96 give to current units of code */
631 rfield.totpow[p.m_nqh] = (fnuHM96[2] + log10(PI4) + scale);
632
633 /* set radius to very large value if not already set */
634 /* >>chng 01 jul 24, from Radius == 0 to this, as per PvH comments */
635 if( !radius.lgRadiusKnown )
636 {
637 radius.Radius = pow(10.,radius.rdfalt);
638 }
639 ++p.m_nqh;
640
641 /* set radius to very large value if not already set */
642 /* >>chng 01 jul 24, from Radius == 0 to this, as per PvH comments */
643 if( !radius.lgRadiusKnown )
644 {
645 radius.Radius = pow(10.,radius.rdfalt);
646 }
647 }
648
649 else if( p.nMatchErase("HM05") )
650 {
651 double *tnuHM , *fnuHM;
652 double redshift;
653 long int nhm;
654 /* the Haardt & Madau 2005 continuum, read in a table and interpolate
655 * for any redshift, background includes both quasars and galaxies
656 * >>refer continuum background Haardt, Francesco, & Madau, Piero, 2005, in preparation
657 * this command does not include the CMB - do that separately with the CMB command
658 * set flag telling interpolate */
659 strcpy( rfield.chSpType[rfield.nShape], "INTER" );
660 if( p.nMatch("QUAS") )
661 {
662 /* quasar-only continuum */
663 strcpy( chFile , "haardt_madau_quasar.dat" );
664 }
665 else
666 {
667 /* the default, quasar plus galaxy continuum */
668 strcpy( chFile , "haardt_madau_galaxy.dat" );
669 }
670
671 /* find the redshift - must be within bounds of table, which we
672 * do not know yet */
673 redshift = p.FFmtRead();
674 if( p.lgEOL() )
675 {
676 fprintf(ioQQQ," The redshift MUST be specified on the table HM05 command. I did not find one.\n Sorry.\n");
678 }
679
680 /* read in the data files and interpolate the continuum to
681 * the correct redshift */
682 read_hm05( chFile , &tnuHM , &fnuHM , &nhm , redshift );
683 /* now find a point near 1 Ryd where continuum may not be far too faint */
684 ipNORM = -1;
685 for( j=0; j < nhm-1; j++ )
686 {
687 /* looking for continuum frequency near 1 Ryd
688 * does not need to be precise */
689 if( tnuHM[j]<=1. && 1. <= tnuHM[j+1] )
690 {
691 ipNORM = j;
692 break;
693 }
694 }
695 ASSERT( ipNORM>=0 );
696
697 /* optional scale factor to change default intensity from their value
698 * assumed to be log if negative, and linear otherwise
699 * increment i to move past the 96 in the keyword */
700 scale = p.FFmtRead();
701 if( scale > 0. )
702 scale = log10(scale);
703
704 /* this will be flux density at some frequency on the table. the numbers
705 * are per Hz and sr so must multiply by 4 pi */
706 rfield.range[p.m_nqh][0] = tnuHM[ipNORM];
707
708 /* convert intensity HM96 give to current units of code */
709 rfield.totpow[p.m_nqh] = log10(fnuHM[ipNORM]) + log10(PI4) + scale;
710
711 /* this also sets continuum intensity*/
712 if( p.m_nqh >= LIMSPC )
713 {
714 fprintf( ioQQQ, " %ld is too many continua entered. Increase LIMSPC\n Sorry.\n",
715 p.m_nqh );
717 }
718
719 /* check that stack of shape and luminosity specifications
720 * is parallel, stop if not - this happens is background comes
721 * BETWEEN another set of shape and luminosity commands */
722 if( rfield.nShape != p.m_nqh )
723 {
724 fprintf( ioQQQ, " This command has come between a previous ordered pair of continuum shape and luminosity commands.\n Reorder the commands to complete each continuum specification before starting another.\n" );
725 fprintf( ioQQQ, " Sorry.\n" );
727 }
728
729 strcpy( rfield.chRSpec[p.m_nqh], "SQCM" );
730 strcpy( rfield.chSpNorm[p.m_nqh], "FLUX" );
731
732 /* this is an isotropic radiation field */
733 rfield.lgBeamed[p.m_nqh] = false;
734 rfield.Illumination[p.m_nqh] = Illuminate::SYMMETRIC;
735
736 /* many of their values are outside the range of a 32 bit cpu and we don't
737 * want to fill array with zeros - so rescale to cont near 1 ryd */
738 scale = SDIV( fnuHM[ipNORM] );
739 /* their table does not extend to the low-energy limit of the code
740 * usually does not matter since CMB will take over, but there is
741 * an obvious gap at low, z = 0, redshift. assume Rayleigh-Jeans tail
742 * and add a first continuum point */
743 rfield.tNu[rfield.nShape][0].set(rfield.emm);
744 rfield.tslop[rfield.nShape][0] =
745 (realnum)log10(fnuHM[0]*pow((double)(rfield.emm/tnuHM[0]),2.)/scale);
746 /*fprintf(ioQQQ,"DEBUG2\t%.3e\t%.3e\n",
747 rfield.tNuRyd[rfield.nShape][0] ,
748 rfield.tslop[rfield.nShape][0] );*/
749 for( j=0; j < nhm; j++ )
750 {
751 /* frequency is already Ryd */
752 rfield.tNu[rfield.nShape][j+1].set(tnuHM[j]);
753 rfield.tslop[rfield.nShape][j+1] = (realnum)log10(fnuHM[j]/scale);
754 /*fprintf(ioQQQ,"DEBUG2\t%.3e\t%.3e\n",
755 rfield.tNuRyd[rfield.nShape][j+1] ,
756 rfield.tslop[rfield.nShape][j+1] );*/
757 }
758 ++nhm;
759 rfield.ncont[rfield.nShape] = nhm;
760
761 /* now make sure they are in increasing order */
762 for( j=0; j < nhm-1; j++ )
763 ASSERT( rfield.tNu[rfield.nShape][j].Ryd() < rfield.tNu[rfield.nShape][j+1].Ryd() );
764
765 /* set radius to very large value if not already set */
766 /* >>chng 01 jul 24, from Radius == 0 to this, as per PvH comments */
767 if( !radius.lgRadiusKnown )
768 {
769 radius.Radius = pow(10.,radius.rdfalt);
770 }
771 ++p.m_nqh;
772
773 /* set radius to very large value if not already set */
774 /* >>chng 01 jul 24, from Radius == 0 to this, as per PvH comments */
775 if( !radius.lgRadiusKnown )
776 {
777 radius.Radius = pow(10.,radius.rdfalt);
778 }
779 free( tnuHM );
780 free( fnuHM );
781
782 }
783 else if( p.nMatch(" ISM") )
784 {
785 ASSERT( NISM < NCELL );
786 /* local ISM radiation field from Black 1987, Interstellar Processes */
787 /* >>chng 04 mar 16, rm CMB from field so that it can be used at
788 * any redshift */
789 rfield.tNu[rfield.nShape][0].set(6.);
790 rfield.tslop[rfield.nShape][0] = -21.21f - 6.f;
791 for( i=6; i < NISM; i++ )
792 {
793 rfield.tNu[rfield.nShape][i-5].set(tnuism[i]);
794 rfield.tslop[rfield.nShape][i-5] = (realnum)(fnuism[i] -
795 tnuism[i]);
796 }
797
798 rfield.ncont[rfield.nShape] = NISM -5;
799
800 /* optional scale factor to change default luminosity
801 * from observed value
802 * want final number to be log
803 * assumed to be log if negative, and linear otherwise unless log option is present */
804 scale = p.FFmtRead();
805 if( scale > 0. && !p.nMatch(" LOG") )
806 scale = log10(scale);
807
808 /* this also sets continuum intensity*/
809 if( p.m_nqh >= LIMSPC )
810 {
811 fprintf( ioQQQ, " %4ld is too many continua entered. Increase LIMSPC\n Sorry.\n",
812 p.m_nqh );
814 }
815
816 /* check that stack of shape and luminosity specifications
817 * is parallel, stop if not - this happens is background comes
818 * BETWEEN another set of shape and luminosity commands */
819 if( rfield.nShape != p.m_nqh )
820 {
821 fprintf( ioQQQ, " This command has come between a previous ordered pair of continuum shape and luminosity commands.\n Reorder the commands to complete each continuum specification before starting another.\n" );
822 fprintf( ioQQQ, " Sorry.\n" );
824 }
825
826 strcpy( rfield.chRSpec[p.m_nqh], "SQCM" );
827 strcpy( rfield.chSpNorm[p.m_nqh], "FLUX" );
828
829 /* this is an isotropic radiation field */
830 rfield.lgBeamed[p.m_nqh] = false;
831 rfield.Illumination[p.m_nqh] = Illuminate::SYMMETRIC;
832
833 /* this will be flux density at 1 Ryd
834 * >>chng 96 dec 18, from 1 Ryd to H mass Rydberg
835 * >>chng 97 jan 10, had HLevNIonRyd but not defined yet */
836 rfield.range[p.m_nqh][0] = HIONPOT;
837
838 /* interpolated from Black 1987 */
839 rfield.totpow[p.m_nqh] = (-18.517 + scale);
840
841 /* set radius to very large value if not already set */
842 /* >>chng 01 jul 24, from Radius == 0 to this, as per PvH comments */
843 if( !radius.lgRadiusKnown )
844 {
845 radius.Radius = pow(10.,radius.rdfalt);
846 }
847 ++p.m_nqh;
848
849 if( optimize.lgVarOn )
850 {
851 optimize.nvarxt[optimize.nparm] = 1;
852 strcpy( optimize.chVarFmt[optimize.nparm], "TABLE ISM LOG %f");
853 /* pointer to where to write */
854 optimize.nvfpnt[optimize.nparm] = input.nRead;
855 /* the scale factor */
856 optimize.vparm[0][optimize.nparm] = (realnum)scale;
857 optimize.vincr[optimize.nparm] = 0.2f;
858 ++optimize.nparm;
859 }
860 }
861 else if( p.nMatch("DRAI") )
862 {
863 ASSERT( NDRAINE < NCELL );
864 rfield.lgMustBlockHIon = true;
865 /* local ISM radiation field from equation 23
866 *>>refer ISM continuum Draine & Bertoldi 1996 */
867 for( i=0; i < NDRAINE; i++ )
868 {
869 rfield.tNu[rfield.nShape][i].set(tnudrn[i]);
870 rfield.tslop[rfield.nShape][i] = (realnum)tsldrn[i];
871 }
872
873 rfield.ncont[rfield.nShape] = NDRAINE;
874
875 /* optional scale factor to change default luminosity
876 * from observed value
877 * assumed to be log if negative, and linear otherwise unless log option is present */
878 scale = p.FFmtRead();
879 if( scale > 0. && !p.nMatch(" LOG") )
880 scale = log10(scale);
881
882 /* this also sets continuum intensity*/
883 if( p.m_nqh >= LIMSPC )
884 {
885 fprintf( ioQQQ, " %4ld is too many continua entered. Increase LIMSPC\n Sorry.\n",
886 p.m_nqh );
888 }
889
890 /* check that stack of shape and luminosity specifications
891 * is parallel, stop if not - this happens is background comes
892 * BETWEEN another set of shape and luminosity commands */
893 if( rfield.nShape != p.m_nqh )
894 {
895 fprintf( ioQQQ, " This command has come between a previous ordered pair of continuum shape and luminosity commands.\n Reorder the commands to complete each continuum specification before starting another.\n" );
896 fprintf( ioQQQ, " Sorry.\n" );
898 }
899
900 strcpy( rfield.chRSpec[p.m_nqh], "SQCM" );
901 strcpy( rfield.chSpNorm[p.m_nqh], "FLUX" );
902
903 /* this is an isotropic radiation field */
904 rfield.lgBeamed[p.m_nqh] = false;
905 rfield.Illumination[p.m_nqh] = Illuminate::SYMMETRIC;
906
907 /* continuum normalization given by flux density at first point,
908 * must set energy a bit higher to make sure it is within energy bounds
909 * that results from float arithmetic */
910 rfield.range[p.m_nqh][0] = tnudrn[0]*1.01;
911
912 /* this is f_nu at this first point */
913 rfield.totpow[p.m_nqh] = tsldrn[0] + scale;
914
915 if( optimize.lgVarOn )
916 {
917 optimize.nvarxt[optimize.nparm] = 1;
918 strcpy( optimize.chVarFmt[optimize.nparm], "TABLE DRAINE LOG %f");
919 /* pointer to where to write */
920 optimize.nvfpnt[optimize.nparm] = input.nRead;
921 /* the scale factor */
922 optimize.vparm[0][optimize.nparm] = (realnum)scale;
923 optimize.vincr[optimize.nparm] = 0.2f;
924 ++optimize.nparm;
925 }
926
927 /* set radius to very large value if not already set */
928 /* >>chng 01 jul 24, from Radius == 0 to this, as per PvH comments */
929 if( !radius.lgRadiusKnown )
930 {
931 radius.Radius = pow(10.,radius.rdfalt);
932 }
933 ++p.m_nqh;
934 }
935
936 else if( p.nMatch("LINE") )
937 {
938 /* table lines command - way to check that lines within a data
939 * file are still valid */
940
941 /* say that this is not a continuum command, so don't try to work with unallocated space */
942 /* this is not a continuum source - it is to read a table of lines */
943 lgNoContinuum = true;
944
945 if( chLINE_LIST.size() > 0 )
946 {
947 fprintf(ioQQQ," sorry, only one table line per input stream\n");
949 }
950
951 /* get file name within double quotes, if not present will use default
952 * return value of 1 indicates did not find double quotes on line */
953 if( lgQuoteFound && strlen(chFile) > 0 )
954 chLINE_LIST = chFile;
955 else
956 chLINE_LIST = "LineList_BLR.dat";
957
958 /* this flag says this is not a continuum source - nearly all table XXX
959 * commands deal with continuum sources */
960 rfield.ncont[rfield.nShape] = -9;
961
962 // check if the file exists
963 FILE* ioData = open_data( chLINE_LIST.c_str(), "r", AS_LOCAL_DATA_TRY );
964 if( ioData == NULL )
965 {
966 /* did not find file, abort */
967 fprintf(ioQQQ,"\n DISASTER PROBLEM ParseTable could not find "
968 "line list file %s\n", chLINE_LIST.c_str() );
969 fprintf(ioQQQ," Please check the spelling of the file name and that it "
970 "is in either the local or data directory.\n\n");
972 }
973 else
974 {
975 fclose(ioData);
976 }
977 /* actually reading the data is done in lines_table() */
978 }
979
980 else if( p.nMatch("POWE") )
981 {
982 /* simple power law continuum between 10 micron and 50 keV
983 * option to read in any slope for the intermediate continuum */
984 alpha = p.FFmtRead();
985
986 /* default (no number on line) is f_nu proportional nu^-1 */
987 if( p.lgEOL() )
988 alpha = -1.;
989
990 /* this is low energy for code */
991 rfield.tNu[rfield.nShape][0].set(rfield.emm);
992 /* and the value of the flux at this point (f_nu units)*/
993 rfield.tslop[rfield.nShape][0] = -5.f;
994
995 lgLogSet = false;
996
997 /* option to adjust sub-millimeter break */
998 brakmm = p.FFmtRead();
999
1000 /* default is 10 microns */
1001 if( p.lgEOL() )
1002 {
1003 lgLogSet = true;
1004 brakmm = 9.115e-3;
1005 }
1006
1007 else if( brakmm == 0. )
1008 {
1009 /* if second number on line is zero then set lower limit to
1010 * low-energy limit of the code. Also set linear mode,
1011 * so that last number will also be linear. */
1012 lgLogSet = false;
1013 brakmm = rfield.tNu[rfield.nShape][0].Ryd()*1.001;
1014 }
1015
1016 else if( brakmm < 0. )
1017 {
1018 /* if number is negative then this and next are logs */
1019 lgLogSet = true;
1020 brakmm = pow(10.,brakmm);
1021 }
1022
1023 /* optional microns keyword - convert to Rydbergs */
1024 if( p.nMatch("MICR") )
1025 brakmm = RYDLAM / (1e4*brakmm);
1026
1027 rfield.tNu[rfield.nShape][1].set(brakmm);
1028
1029 /* check whether this is a reasonable mm break */
1030 if( brakmm > 1. )
1031 fprintf(ioQQQ,
1032 " Check the order of parameters on this table power law command - the low-energy break of %f Ryd seems high to me.\n",
1033 brakmm );
1034
1035 /* this is spectral slope, in F_nu units, between the low energy limit
1036 * and the break that may have been adjusted above
1037 * this is the slope appropriate for self-absorbed synchrotron, see eq 6.54, p.190
1038 *>>refer continuum synchrotron Rybicki, G. B., & Lightman, A.P. 1979,
1039 *>>refercon Radiative Processes in Astrophysics (New York: Wiley)*/
1040 slopir = 5./2.;
1041
1042 /* now extrapolate a flux at this energy using the flux entered for
1043 * the first point, and this slope */
1044 rfield.tslop[rfield.nShape][1] =
1045 (realnum)(rfield.tslop[rfield.nShape][0] +
1046 slopir*log10(rfield.tNu[rfield.nShape][1].Ryd()/rfield.tNu[rfield.nShape][0].Ryd()));
1047
1048 /* this is energy of the high-energy limit to code */
1049 rfield.tNu[rfield.nShape][3].set(rfield.egamry);
1050
1051 /* option to adjust hard X-ray break */
1052 brakxr = p.FFmtRead();
1053
1054 /* default is 50 keV */
1055 if( p.lgEOL() )
1056 {
1057 brakxr = 3676.;
1058 }
1059
1060 else if( brakxr == 0. )
1061 {
1062 brakxr = rfield.tNu[rfield.nShape][3].Ryd()/1.001;
1063 }
1064
1065 else if( lgLogSet )
1066 {
1067 /* first number was negative this is a logs */
1068 brakxr = pow(10.,brakxr);
1069 }
1070
1071 /* note that this second cutoff does not have the micron keyword */
1072 rfield.tNu[rfield.nShape][2].set(brakxr);
1073
1074 /* >>chng 03 jul 19, check that upper energy is greater than lower energy,
1075 * quit if this is not the case */
1076 if( brakmm >= brakxr )
1077 {
1078 fprintf( ioQQQ, " HELP!! The lower energy for the power law, %f, is greater than the upper energy, %f. This is not possible.\n Sorry.\n",
1079 brakmm , brakxr );
1081 }
1082
1083 /* alpha was first option on line, is slope of mid-range */
1084 rfield.tslop[rfield.nShape][2] =
1085 (realnum)(rfield.tslop[rfield.nShape][1] +
1086 alpha*log10(rfield.tNu[rfield.nShape][2].Ryd()/rfield.tNu[rfield.nShape][1].Ryd()));
1087
1088 /* high energy range is nu^-2 */
1089 slopxr = -2.;
1090
1091 rfield.tslop[rfield.nShape][3] =
1092 (realnum)(rfield.tslop[rfield.nShape][2] +
1093 slopxr*log10(rfield.tNu[rfield.nShape][3].Ryd()/rfield.tNu[rfield.nShape][2].Ryd()));
1094
1095 /* following is number of portions of continuum */
1096 rfield.ncont[rfield.nShape] = 4;
1097 }
1098
1099 else if( p.nMatch("READ") )
1100 {
1101 /* set up eventual read of table of points previously punched by code
1102 * get file name within double quotes, return as null terminated string
1103 * also blank out original, chCard version of name so do not trigger */
1104 if( !lgQuoteFound )
1105 {
1106 fprintf( ioQQQ, " Name of file must appear on TABLE READ.\n");
1108 }
1109
1110 /* set flag saying really just read in continuum exactly as punched */
1111 strcpy( rfield.chSpType[rfield.nShape], "READ " );
1112
1113 ReadTable(chFile);
1114
1115 /* this is flag saying continuum has been read in here */
1116 rfield.tslop[rfield.nShape][0] = 1.;
1117 rfield.tslop[rfield.nShape][1] = 0.;
1118
1119 /* number of spectra shapes that have been specified */
1120 ++rfield.nShape;
1121 return;
1122 }
1123
1124 else if( p.nMatch("TLUSTY") && !p.nMatch("STAR") )
1125 {
1126 /* >>chng 04 nov 30, retired TABLE TLUSTY command, PvH */
1127 fprintf( ioQQQ, " The TABLE TLUSTY command is no longer supported.\n" );
1128 fprintf( ioQQQ, " Please use TABLE STAR TLUSTY instead. See Hazy for details.\n" );
1130 }
1131
1132 else if( p.nMatch("RUBI") )
1133 {
1134 /* do Rubin modified theta 1 Ori c */
1135 for( i=0; i < NRUBIN; i++ )
1136 {
1137 rfield.tNu[rfield.nShape][i].set(tnurbn[i]);
1138 rfield.tslop[rfield.nShape][i] = (realnum)log10(fnurbn[i] /tnurbn[i] );
1139 }
1140 rfield.ncont[rfield.nShape] = NRUBIN;
1141 }
1142
1143 /* >>chng 06 jul 10, retired TABLE STARBURST command, PvH */
1144
1145 else if( p.nMatch("STAR") )
1146 {
1147 char chMetalicity[5] = "", chODFNew[8], chVaryFlag[7] = "";
1148 bool lgHCa = false, lgHNi = false;
1149 long nval, ndim=0;
1150 double Tlow = -1., Thigh = -1.;
1151 double val[MDIM];
1152
1153 /* >>chng 06 jun 22, add support for 3 and 4-dimensional grids, PvH */
1154 if( p.nMatchErase("1-DI") )
1155 ndim = 1;
1156 else if( p.nMatchErase("2-DI") )
1157 ndim = 2;
1158 else if( p.nMatchErase("3-DI") )
1159 ndim = 3;
1160 else if( p.nMatchErase("4-DI") )
1161 ndim = 4;
1162
1163 if( ndim != 0 )
1164 {
1165 /* remember keyword for possible use in optimization command */
1166 sprintf(chVaryFlag,"%1ld-DIM",ndim);
1167 }
1168
1169 /* time option to vary only some continua with time */
1170 rfield.lgTimeVary[p.m_nqh] = p.nMatch( "TIME" );
1171
1172 static const char table[][2][10] = {
1173 {"Z+1.0 ", "p10"},
1174 {"Z+0.75", "p075"},
1175 {"Z+0.5 ", "p05"},
1176 {"Z+0.4 ", "p04"},
1177 {"Z+0.3 ", "p03"},
1178 {"Z+0.25", "p025"},
1179 {"Z+0.2 ", "p02"},
1180 {"Z+0.1 ", "p01"},
1181 {"Z+0.0 ", "p00"},
1182 {"Z-0.1 ", "m01"},
1183 {"Z-0.2 ", "m02"},
1184 {"Z-0.25", "m025"},
1185 {"Z-0.3 ", "m03"},
1186 {"Z-0.4 ", "m04"},
1187 {"Z-0.5 ", "m05"},
1188 {"Z-0.7 ", "m07"},
1189 {"Z-0.75", "m075"},
1190 {"Z-1.0 ", "m10"},
1191 {"Z-1.3 ", "m13"},
1192 {"Z-1.5 ", "m15"},
1193 {"Z-1.7 ", "m17"},
1194 {"Z-2.0 ", "m20"},
1195 {"Z-2.5 ", "m25"},
1196 {"Z-3.0 ", "m30"},
1197 {"Z-3.5 ", "m35"},
1198 {"Z-4.0 ", "m40"},
1199 {"Z-4.5 ", "m45"},
1200 {"Z-5.0 ", "m50"},
1201 {"Z-INF ", "m99"}
1202 };
1203
1204 strncpy( chMetalicity, "p00", sizeof(chMetalicity) ); // default
1205 for (i=0; i < (long)(sizeof(table)/sizeof(*table)); ++i)
1206 {
1207 if (p.nMatchErase(table[i][0]))
1208 {
1209 strncpy( chVaryFlag, table[i][0], sizeof(chVaryFlag) );
1210 strncpy( chMetalicity, table[i][1], sizeof(chMetalicity) );
1211 break;
1212 }
1213 }
1214
1215
1216 /* there are two abundance sets (solar and halo) for CoStar and Rauch H-Ca/H-Ni models.
1217 * If halo keyword appears use halo, else use solar unless other metalicity was requested */
1218 bool lgHalo = p.nMatch("HALO");
1219 bool lgSolar = ( !lgHalo && strcmp( chMetalicity, "p00" ) == 0 );
1220
1221 /* >>chng 05 oct 27, added support for PG1159 Rauch models, PvH */
1222 bool lgPG1159 = p.nMatchErase("PG1159");
1223
1224 bool lgList = p.nMatch("LIST");
1225
1226 if( p.nMatch("AVAI") )
1227 {
1230 }
1231
1232 /* now that all the keywords are out of the way, scan numbers from line image */
1233 for( nval=0; nval < MDIM; nval++ )
1234 {
1235 val[nval] = p.FFmtRead();
1236 if( p.lgEOL() )
1237 break;
1238 }
1239
1240 if( nval == 0 && !lgList )
1241 {
1242 fprintf( ioQQQ, " The stellar temperature MUST be entered.\n" );
1244 }
1245
1246 /* option for log if less than or equal to 10 */
1247 /* Caution: For CoStar models this implicitly assumes that
1248 * the minimum ZAMS mass and youngest age is greater than 10.
1249 * As of June 1999 this is the case. However, this is not
1250 * guaranteed for possible future upgrades. */
1251 /* match on "LOG " rather than " LOG" to avoid confusion with "LOG(G)" */
1252 if( ( val[0] <= 10. && !p.nMatch("LINE") ) || p.nMatch("LOG ") )
1253 {
1254 if( val[0] < log10(BIGFLOAT) )
1255 val[0] = pow(10.,val[0]);
1256 else
1257 {
1258 fprintf(ioQQQ," DISASTER the log of the temperature was specified, "
1259 "but the numerical value is too large.\n Sorry.\n\n");
1261 }
1262 }
1263
1264 if( lgQuoteFound )
1265 {
1266 nstar = GridInterpolate(val,&nval,&ndim,chFile,lgList,&Tlow,&Thigh);
1267 }
1268
1269 else if( p.nMatch("ATLA") )
1270 {
1271 /* this sub-branch added by Kevin Volk, to read in large
1272 * grid of Kurucz atmospheres */
1273 /* >>chng 05 nov 19, option TRACE is no longer supported, PvH */
1274
1275 if( p.nMatch("ODFN") )
1276 strncpy( chODFNew, "_odfnew", sizeof(chODFNew) );
1277 else
1278 strncpy( chODFNew, "", sizeof(chODFNew) );
1279
1280 /* >>chng 05 nov 19, add support for non-solar metalicities, PvH */
1281 nstar = AtlasInterpolate(val,&nval,&ndim,chMetalicity,chODFNew,lgList,&Tlow,&Thigh);
1282 }
1283
1284 else if( p.nMatch("COST") )
1285 {
1286 /* >>chng 99 apr 30, added CoStar stellar atmospheres */
1287 /* this subbranch reads in CoStar atmospheres, no log(g),
1288 * but second parameter is age sequence, a number between 1 and 7,
1289 * default is 1 */
1290 if( p.nMatch("INDE") )
1291 {
1292 imode = IM_COSTAR_TEFF_MODID;
1293 if( nval == 1 )
1294 {
1295 val[1] = 1.;
1296 nval++;
1297 }
1298
1299 /* now make sure that second parameter is within allowed range -
1300 * we do not have enough information at this time to verify temperature */
1301 if( val[1] < 1. || val[1] > 7. )
1302 {
1303 fprintf( ioQQQ, " The second number must be the id sequence number, 1 to 7.\n" );
1304 fprintf( ioQQQ, " reset to 1.\n" );
1305 val[1] = 1.;
1306 }
1307 }
1308 else if( p.nMatch("ZAMS") )
1309 {
1310 imode = IM_COSTAR_MZAMS_AGE;
1311 if( nval == 1 )
1312 {
1313 fprintf( ioQQQ, " A second number, the age of the star, must be present.\n" );
1315 }
1316 }
1317 else if( p.nMatch(" AGE") )
1318 {
1319 imode = IM_COSTAR_AGE_MZAMS;
1320 if( nval == 1 )
1321 {
1322 fprintf( ioQQQ, " A second number, the ZAMS mass of the star, must be present.\n" );
1324 }
1325 }
1326 else
1327 {
1328 if( nval == 1 )
1329 {
1330 imode = IM_COSTAR_TEFF_MODID;
1331 /* default is to use ZAMS models, i.e. use index 1 */
1332 val[1] = 1.;
1333 nval++;
1334 }
1335 else
1336 {
1337 /* Caution: the code in CoStarInterpolate implicitly
1338 * assumes that the dependence of log(g) with age is
1339 * strictly monotonic. As of June 1999 this is the case. */
1340 imode = IM_COSTAR_TEFF_LOGG;
1341 }
1342 }
1343
1344 if( ! ( lgSolar || lgHalo ) )
1345 {
1346 fprintf( ioQQQ, " Please choose SOLAR or HALO abundances.\n" );
1348 }
1349
1350 nstar = CoStarInterpolate(val,&nval,&ndim,imode,lgHalo,lgList,&Tlow,&Thigh);
1351 }
1352
1353 else if( p.nMatch("KURU") )
1354 {
1355 nstar = Kurucz79Interpolate(val,&nval,&ndim,lgList,&Tlow,&Thigh);
1356 }
1357
1358 else if( p.nMatch("MIHA") )
1359 {
1360 nstar = MihalasInterpolate(val,&nval,&ndim,lgList,&Tlow,&Thigh);
1361 }
1362
1363 else if( p.nMatch("RAUC") )
1364 {
1365 if( ! ( lgSolar || lgHalo ) )
1366 {
1367 fprintf( ioQQQ, " Please choose SOLAR or HALO abundances.\n" );
1369 }
1370
1371 /* >>chng 97 aug 23, K Volk added Rauch stellar atmospheres */
1372 /* >>chng 05 oct 27, added support for PG1159 Rauch models, PvH */
1373 /* >>chng 06 jun 26, added support for H, He, H+He Rauch models, PvH */
1374 if( p.nMatch("H-CA") || p.nMatch(" OLD") )
1375 {
1376 lgHCa = true;
1377 nstar = RauchInterpolateHCa(val,&nval,&ndim,lgHalo,lgList,&Tlow,&Thigh);
1378 }
1379 else if( p.nMatch("HYDR") )
1380 {
1381 nstar = RauchInterpolateHydr(val,&nval,&ndim,lgList,&Tlow,&Thigh);
1382 }
1383 else if( p.nMatch("HELI") )
1384 {
1385 nstar = RauchInterpolateHelium(val,&nval,&ndim,lgList,&Tlow,&Thigh);
1386 }
1387 else if( p.nMatch("H+HE") )
1388 {
1389 nstar = RauchInterpolateHpHe(val,&nval,&ndim,lgList,&Tlow,&Thigh);
1390 }
1391 else if( lgPG1159 ) /* the keyword PG1159 was already matched and erased above */
1392 {
1393 nstar = RauchInterpolatePG1159(val,&nval,&ndim,lgList,&Tlow,&Thigh);
1394 }
1395 else if( p.nMatch("CO W") )
1396 {
1397 nstar = RauchInterpolateCOWD(val,&nval,&ndim,lgList,&Tlow,&Thigh);
1398 }
1399 else
1400 {
1401 lgHNi = true;
1402 nstar = RauchInterpolateHNi(val,&nval,&ndim,lgHalo,lgList,&Tlow,&Thigh);
1403 }
1404 }
1405
1406 else if( p.nMatch("TLUS") )
1407 {
1408 if( p.nMatch("OBST") )
1409 {
1410 /* >>chng 09 dec 24, this sub-branch added to read in
1411 * merged Tlusty O-star & B-star atmospheres, PvH */
1412 nstar = TlustyInterpolate(val,&nval,&ndim,TL_OBSTAR,chMetalicity,lgList,&Tlow,&Thigh);
1413 }
1414 else if( p.nMatch("BSTA") )
1415 {
1416 /* >>chng 06 oct 19, this sub-branch added to read in
1417 * large 2006 grid of Tlusty B-star atmospheres, PvH */
1418 nstar = TlustyInterpolate(val,&nval,&ndim,TL_BSTAR,chMetalicity,lgList,&Tlow,&Thigh);
1419 }
1420 else if( p.nMatch("OSTA") )
1421 {
1422 /* >>chng 05 nov 21, this sub-branch added to read in
1423 * large 2002 grid of Tlusty O-star atmospheres, PvH */
1424 nstar = TlustyInterpolate(val,&nval,&ndim,TL_OSTAR,chMetalicity,lgList,&Tlow,&Thigh);
1425 }
1426 else
1427 {
1428 fprintf( ioQQQ, " There must be a third key on TABLE STAR TLUSTY;" );
1429 fprintf( ioQQQ, " the options are OBSTAR, BSTAR, OSTAR.\n" );
1431 }
1432 }
1433
1434 else if( p.nMatch("WERN") )
1435 {
1436 /* this subbranch added by Kevin Volk, to read in large
1437 * grid of hot PN atmospheres computed by Klaus Werner */
1438 nstar = WernerInterpolate(val,&nval,&ndim,lgList,&Tlow,&Thigh);
1439 }
1440
1441 else if( p.nMatch("WMBA") )
1442 {
1443 /* >>chng 06 jun 27, this subbranch added to read in
1444 * grid of hot atmospheres computed by Pauldrach */
1445 nstar = WMBASICInterpolate(val,&nval,&ndim,lgList,&Tlow,&Thigh);
1446 }
1447
1448 else
1449 {
1450 fprintf( ioQQQ, " There must be a second key on TABLE STAR;" );
1451 fprintf( ioQQQ, " the options are ATLAS, KURUCZ, MIHALAS, RAUCH, WERNER, and WMBASIC.\n" );
1453 }
1454
1455 /* set flag saying really just read in continuum exactly as punched */
1456 strcpy( rfield.chSpType[rfield.nShape], "VOLK " );
1457
1458 rfield.ncont[rfield.nShape] = nstar;
1459
1460 /* vary option */
1461 if( optimize.lgVarOn )
1462 {
1463 optimize.vincr[optimize.nparm] = (realnum)0.1;
1464
1465 if( lgQuoteFound )
1466 {
1467 /* this is number of parameters to feed onto the input line */
1468 optimize.nvarxt[optimize.nparm] = nval;
1469 sprintf( optimize.chVarFmt[optimize.nparm], "TABLE STAR \"%s\"", chFile );
1470 strcat( optimize.chVarFmt[optimize.nparm], " %f LOG" );
1471 for( i=1; i < nval; i++ )
1472 strcat( optimize.chVarFmt[optimize.nparm], " %f" );
1473 }
1474
1475 if( p.nMatch("ATLA") )
1476 {
1477 /* this is number of parameters to feed onto the input line */
1478 optimize.nvarxt[optimize.nparm] = ndim;
1479 strcpy( optimize.chVarFmt[optimize.nparm], "TABLE STAR ATLAS " );
1480 strcat( optimize.chVarFmt[optimize.nparm], chVaryFlag );
1481 if( p.nMatch("ODFN") )
1482 strcat( optimize.chVarFmt[optimize.nparm], " ODFNEW" );
1483
1484 strcat( optimize.chVarFmt[optimize.nparm], " %f LOG %f" );
1485
1486 if( ndim == 3 )
1487 strcat( optimize.chVarFmt[optimize.nparm], " %f" );
1488
1489 }
1490
1491 else if( p.nMatch("COST") )
1492 {
1493 /* this is number of parameters to feed onto the input line */
1494 optimize.nvarxt[optimize.nparm] = 2;
1495
1496 strcpy( optimize.chVarFmt[optimize.nparm], "TABLE STAR COSTAR " );
1497 if( lgHalo )
1498 strcat( optimize.chVarFmt[optimize.nparm], "HALO " );
1499
1500 if( imode == IM_COSTAR_TEFF_MODID )
1501 {
1502 strcat( optimize.chVarFmt[optimize.nparm], "%f LOG , INDEX %f" );
1503 }
1504 else if( imode == IM_COSTAR_TEFF_LOGG )
1505 {
1506 strcat( optimize.chVarFmt[optimize.nparm], "%f LOG %f" );
1507 }
1508 else if( imode == IM_COSTAR_MZAMS_AGE )
1509 {
1510 strcat( optimize.chVarFmt[optimize.nparm], "ZAMS %f LOG %f" );
1511 }
1512 else if( imode == IM_COSTAR_AGE_MZAMS )
1513 {
1514 strcat( optimize.chVarFmt[optimize.nparm], "AGE %f LOG %f" );
1515 optimize.vincr[optimize.nparm] = (realnum)0.5;
1516 }
1517 }
1518
1519 else if( p.nMatch("KURU") )
1520 {
1521 /* this is number of parameters to feed onto the input line */
1522 optimize.nvarxt[optimize.nparm] = 1;
1523 strcpy( optimize.chVarFmt[optimize.nparm],
1524 "TABLE STAR KURUCZ %f LOG" );
1525 }
1526
1527 else if( p.nMatch("MIHA") )
1528 {
1529 /* this is number of parameters to feed onto the input line */
1530 optimize.nvarxt[optimize.nparm] = 1;
1531 strcpy( optimize.chVarFmt[optimize.nparm],
1532 "TABLE STAR MIHALAS %f LOG" );
1533 }
1534
1535 else if( p.nMatch("RAUC") )
1536 {
1537 /* this is number of parameters to feed onto the input line */
1538 optimize.nvarxt[optimize.nparm] = ndim;
1539
1540 strcpy( optimize.chVarFmt[optimize.nparm], "TABLE STAR RAUCH " );
1541
1542 if( p.nMatch("HYDR") )
1543 strcat( optimize.chVarFmt[optimize.nparm], "HYDR " );
1544 else if( p.nMatch("HELI") )
1545 strcat( optimize.chVarFmt[optimize.nparm], "HELIUM " );
1546 else if( p.nMatch("H+HE") )
1547 strcat( optimize.chVarFmt[optimize.nparm], "H+HE " );
1548 else if( lgPG1159 )
1549 strcat( optimize.chVarFmt[optimize.nparm], "PG1159 " );
1550 else if( p.nMatch("CO W") )
1551 strcat( optimize.chVarFmt[optimize.nparm], "CO WD " );
1552 else if( lgHCa )
1553 strcat( optimize.chVarFmt[optimize.nparm], "H-CA " );
1554
1555 if( ( lgHCa || lgHNi ) && ndim == 2 )
1556 {
1557 if( lgHalo )
1558 strcat( optimize.chVarFmt[optimize.nparm], "HALO " );
1559 else
1560 strcat( optimize.chVarFmt[optimize.nparm], "SOLAR " );
1561 }
1562
1563 strcat( optimize.chVarFmt[optimize.nparm], "%f LOG %f" );
1564
1565 if( ndim == 3 )
1566 {
1567 if( p.nMatch("H+HE") )
1568 strcat( optimize.chVarFmt[optimize.nparm], " %f" );
1569 else
1570 strcat( optimize.chVarFmt[optimize.nparm], " %f 3-DIM" );
1571 }
1572 }
1573
1574 if( p.nMatch("TLUS") )
1575 {
1576 /* this is number of parameters to feed onto the input line */
1577 optimize.nvarxt[optimize.nparm] = ndim;
1578 strcpy( optimize.chVarFmt[optimize.nparm], "TABLE STAR TLUSTY " );
1579 if( p.nMatch("OBST") )
1580 strcat( optimize.chVarFmt[optimize.nparm], "OBSTAR " );
1581 else if( p.nMatch("BSTA") )
1582 strcat( optimize.chVarFmt[optimize.nparm], "BSTAR " );
1583 else if( p.nMatch("OSTA") )
1584 strcat( optimize.chVarFmt[optimize.nparm], "OSTAR " );
1585 else
1586 TotalInsanity();
1587 strcat( optimize.chVarFmt[optimize.nparm], chVaryFlag );
1588 strcat( optimize.chVarFmt[optimize.nparm], " %f LOG %f" );
1589
1590 if( ndim == 3 )
1591 strcat( optimize.chVarFmt[optimize.nparm], " %f" );
1592 }
1593
1594 else if( p.nMatch("WERN") )
1595 {
1596 /* this is number of parameters to feed onto the input line */
1597 optimize.nvarxt[optimize.nparm] = 2;
1598 strcpy( optimize.chVarFmt[optimize.nparm],
1599 "TABLE STAR WERNER %f LOG %f" );
1600 }
1601
1602 else if( p.nMatch("WMBA") )
1603 {
1604 /* this is number of parameters to feed onto the input line */
1605 optimize.nvarxt[optimize.nparm] = 3;
1606 strcpy( optimize.chVarFmt[optimize.nparm],
1607 "TABLE STAR WMBASIC %f LOG %f %f" );
1608 }
1609
1610 if( rfield.lgTimeVary[p.m_nqh] )
1611 strcat( optimize.chVarFmt[optimize.nparm], " TIME" );
1612
1613 /* pointer to where to write */
1614 optimize.nvfpnt[optimize.nparm] = input.nRead;
1615
1616 ASSERT( nval <= LIMEXT );
1617
1618 /* log of temp will be pointer */
1619 optimize.vparm[0][optimize.nparm] = (realnum)log10(val[0]);
1620 for( i=1; i < nval; i++ )
1621 optimize.vparm[i][optimize.nparm] = (realnum)val[i];
1622
1623 /* following are upper and lower limits to temperature range */
1624 optimize.varang[optimize.nparm][0] = (realnum)log10(Tlow);
1625 optimize.varang[optimize.nparm][1] = (realnum)log10(Thigh);
1626
1627 /* finally increment this counter */
1628 ++optimize.nparm;
1629 }
1630 }
1631
1632 else if( p.nMatch(" XDR") )
1633 {
1634 /* The XDR SED described by
1635 //>>ref XDR Maloney, P.R. and Hollenbach, D.J. \& Tielens, A. G. G. M., 1996, ApJ, 466, 561
1636 */
1637 // low energy limit of code rfield.emm
1638 i = 0;
1639 rfield.tNu[rfield.nShape][i].set(rfield.emm);
1640 rfield.tslop[rfield.nShape][i] = (realnum)log10(SMALLFLOAT);
1641 // their radiation field is only defined from 1 to 100 keV
1642 ++i;
1643 rfield.tNu[rfield.nShape][i].set(1e3 / EVRYD * 0.95);
1644 rfield.tslop[rfield.nShape][i] = (realnum)log10(SMALLFLOAT);
1645 // _nu propto nu^-0.7 1 - 100 keV
1646 ++i;
1647 rfield.tNu[rfield.nShape][i].set(1e3 / EVRYD );
1648 rfield.tslop[rfield.nShape][i] = (realnum)log10(1.);
1649 ++i;
1650 rfield.tNu[rfield.nShape][i].set(1e5 / EVRYD );
1651 rfield.tslop[rfield.nShape][i] = (realnum)log10(pow(100.,-0.7));
1652 ++i;
1653 rfield.tNu[rfield.nShape][i].set(1e5 / EVRYD * 1.05);
1654 rfield.tslop[rfield.nShape][i] = (realnum)log10(SMALLFLOAT);
1655 ++i;
1656 rfield.tNu[rfield.nShape][i].set(rfield.egamry);
1657 rfield.tslop[rfield.nShape][i] = (realnum)log10(SMALLFLOAT);
1658
1659 rfield.ncont[rfield.nShape] = i+1;
1660 }
1661
1662 else
1663 {
1664 fprintf( ioQQQ, " There MUST be a keyword on this line. The keys are:"
1665 "AGN, AKN120, CRAB, COOL, DRAINE, HM96, HM05, _ISM, LINE, POWERlaw, "
1666 "READ, RUBIN, STAR, and XDR.\n Sorry.\n" );
1668 }
1669
1670 /* table star Werner and table star atlas are special
1671 * cases put in by Kevin Volk - they are not really tables
1672 * at all, so return if chSpType is Volk */
1673 if( strcmp(rfield.chSpType[rfield.nShape],"VOLK ") == 0 )
1674 {
1675 ++rfield.nShape;
1676 return;
1677 }
1678
1679 /* this flag set true if we did not parse a continuum source, thus creating
1680 * the arrays that are needed - return at this point */
1681 if( lgNoContinuum )
1682 {
1683 return;
1684 }
1685
1686 /* ncont only positive when real continuum entered
1687 * convert from log(Hz) to Ryd if first nu>5 */
1688 if( rfield.tNu[rfield.nShape][0].Ryd() >= 5. )
1689 {
1690 for( i=0; i < rfield.ncont[rfield.nShape]; i++ )
1691 {
1692 rfield.tNu[rfield.nShape][i].set(
1693 pow(10.,rfield.tNu[rfield.nShape][i].Ryd() - 15.51718));
1694 }
1695 }
1696
1697 else if( rfield.tNu[rfield.nShape][0].Ryd() < 0. )
1698 {
1699 /* energies entered as logs */
1700 for( i=0; i < rfield.ncont[rfield.nShape]; i++ )
1701 {
1702 rfield.tNu[rfield.nShape][i].set(
1703 (realnum)pow(10.,rfield.tNu[rfield.nShape][i].Ryd()));
1704 }
1705 }
1706
1707 /* tFluxLog will be log(fnu) at that spot, read into tslop */
1708 for( i=0; i < rfield.ncont[rfield.nShape]; i++ )
1709 {
1710 rfield.tFluxLog[rfield.nShape][i] = rfield.tslop[rfield.nShape][i];
1711 }
1712
1713 for( i=0; i < rfield.ncont[rfield.nShape]-1; i++ )
1714 {
1715 rfield.tslop[rfield.nShape][i] =
1716 (realnum)((rfield.tslop[rfield.nShape][i+1] -
1717 rfield.tslop[rfield.nShape][i])/
1718 log10(rfield.tNu[rfield.nShape][i+1].Ryd()/
1719 rfield.tNu[rfield.nShape][i].Ryd()));
1720 }
1721
1722 if( rfield.ncont[rfield.nShape] > 1 && (rfield.ncont[rfield.nShape] + 1 < rfield.nupper) )
1723 {
1724 /* zero out remainder of array */
1725 for( i=rfield.ncont[rfield.nShape]; i < rfield.nupper; i++ )
1726 {
1727 rfield.tNu[rfield.nShape][i].set(0.);
1728 }
1729 }
1730
1731 if( trace.lgConBug && trace.lgTrace )
1732 {
1733 fprintf( ioQQQ, " Table for this continuum; TNU, TFAC, TSLOP\n" );
1734 for( i=0; i < rfield.ncont[rfield.nShape]; i++ )
1735 {
1736 fprintf( ioQQQ, "%12.4e %12.4e %12.4e\n", rfield.tNu[rfield.nShape][i].Ryd(),
1737 rfield.tFluxLog[rfield.nShape][i], rfield.tslop[rfield.nShape][i] );
1738 }
1739 }
1740
1741 /* renormalize continuum so that quantity passes through unity at 1 Ryd
1742 * lgHit will be set false when we get a hit in following loop,
1743 * then test made to make sure it happened*/
1744 lgHit = false;
1745 /*following will be reset when hit occurs*/
1746 fac = -DBL_MAX;
1747 /* >>chng 04 mar 16, chng loop so breaks when hit, previously had cycled
1748 * until last point reached, so last point always used */
1749 for( i=0; i < rfield.ncont[rfield.nShape] && !lgHit; i++ )
1750 {
1751 if( rfield.tNu[rfield.nShape][i].Ryd() > 1. )
1752 {
1753 fac = rfield.tFluxLog[rfield.nShape][i];
1754 lgHit = true;
1755 }
1756 }
1757
1758 if( rfield.ncont[rfield.nShape] > 1 && lgHit )
1759 {
1760 /* do the renormalization */
1761 for( i=0; i < rfield.ncont[rfield.nShape] ; i++ )
1762 {
1763 rfield.tFluxLog[rfield.nShape][i] -= (realnum)fac;
1764 }
1765 }
1766
1767 if( rfield.ncont[rfield.nShape] > 1 )
1768 ++rfield.nShape;
1769 return;
1770}
1771
1772
1773
1774/*ZeroContin store sets of built in continua */
1776{
1777
1778 DEBUG_ENTRY( "ZeroContin()" );
1779
1780 /* Draine 1978 ISM continuum */
1781 /* freq in ryd */
1782 tnudrn[0] = 0.3676;
1783 tnudrn[1] = 0.4144;
1784 tnudrn[2] = 0.4558;
1785 tnudrn[3] = 0.5064;
1786 tnudrn[4] = 0.5698;
1787 tnudrn[5] = 0.6511;
1788 tnudrn[6] = 0.7012;
1789 tnudrn[7] = 0.7597;
1790 tnudrn[8] = 0.8220;
1791 tnudrn[9] = 0.9116;
1792 tnudrn[10] = 0.9120;
1793 tnudrn[11] = 0.9306;
1794 tnudrn[12] = 0.9600;
1795 tnudrn[13] = 0.9806;
1796 /* >>chng 05 aug 03, this energy is so close to 1 ryd that it spills over
1797 * into the H-ionizing continuum - move down to 0.99 */
1798 /* >>chng 05 aug 08, this destabilized pdr_leiden_hack_v4 (!) so put back to
1799 * original value and include extinguish command */
1800 tnudrn[14] = 0.9999;
1801 /*tnudrn[14] = 0.99;*/
1802
1803 /* f_nu from equation 23 of
1804 * >>refer ism field Draine, B.T. & Bertoldi, F. 1996, ApJ, 468, 269 */
1805 int i;
1806 i= 0;
1807 tsldrn[i] = -17.8063;
1808 ++i;tsldrn[i] = -17.7575;
1809 ++i;tsldrn[i] = -17.7268;
1810 ++i;tsldrn[i] = -17.7036;
1811 ++i;tsldrn[i] = -17.6953;
1812 ++i;tsldrn[i] = -17.7182;
1813 ++i;tsldrn[i] = -17.7524;
1814 ++i;tsldrn[i] = -17.8154;
1815 ++i;tsldrn[i] = -17.9176;
1816 ++i;tsldrn[i] = -18.1675;
1817 ++i;tsldrn[i] = -18.1690;
1818 ++i;tsldrn[i] = -18.2477;
1819 ++i;tsldrn[i] = -18.4075;
1820 ++i;tsldrn[i] = -18.5624;
1821 ++i;tsldrn[i] = -18.7722;
1822
1823 /* generic AGN continuum taken from
1824 * >>refer AGN cont Mathews and Ferland ApJ Dec15 '87
1825 * except self absorption at 10 microns, nu**-2.5 below that */
1826 tnuagn[0] = 1e-5;
1827 tnuagn[1] = 9.12e-3;
1828 tnuagn[2] = .206;
1829 tnuagn[3] = 1.743;
1830 tnuagn[4] = 4.13;
1831 tnuagn[5] = 26.84;
1832 tnuagn[6] = 7353.;
1833 tnuagn[7] = 7.4e6;
1834
1835 tslagn[0] = -3.388;
1836 tslagn[1] = 4.0115;
1837 tslagn[2] = 2.6576;
1838 tslagn[3] = 2.194;
1839 tslagn[4] = 1.819;
1840 tslagn[5] = -.6192;
1841 tslagn[6] = -2.326;
1842 tslagn[7] = -7.34;
1843 resetBltin( tnuagn , tslagn , true );
1844
1845 /* Crab Nebula continuum from
1846 *>>refer continuum CrabNebula Davidson, K., & Fesen, 1985, ARAA,
1847 * second vector is F_nu as seen from Earth */
1848 tnucrb[0] = 1.0e-5;
1849 tnucrb[1] = 5.2e-4;
1850 tnucrb[2] = 1.5e-3;
1851 tnucrb[3] = 0.11;
1852 tnucrb[4] = 0.73;
1853 tnucrb[5] = 7.3;
1854 tnucrb[6] = 73.;
1855 tnucrb[7] = 7300.;
1856 tnucrb[8] = 1.5e6;
1857 tnucrb[9] = 7.4e6;
1858
1859 fnucrb[0] = 3.77e-21;
1860 fnucrb[1] = 1.38e-21;
1861 fnucrb[2] = 2.10e-21;
1862 fnucrb[3] = 4.92e-23;
1863 fnucrb[4] = 1.90e-23;
1864 fnucrb[5] = 2.24e-24;
1865 fnucrb[6] = 6.42e-26;
1866 fnucrb[7] = 4.02e-28;
1867 fnucrb[8] = 2.08e-31;
1868 fnucrb[9] = 1.66e-32;
1869 resetBltin( tnucrb , fnucrb , false );
1870
1871 /* Bob Rubin's theta 1 Ori C continuum, modified from Kurucz to
1872 * get NeIII right */
1873 /* >>chng 02 may 02, revise continuum while working with Bob Rubin on NII revisit */
1874 resetBltin( tnurbn , fnurbn , false );
1875
1876 /* cooling flow continuum from Johnstone et al. MNRAS 255, 431. */
1877 cfle[0] = 0.0000100;
1878 cflf[0] = -0.8046910;
1879 cfle[1] = 0.7354023;
1880 cflf[1] = -0.8046910;
1881 cfle[2] = 1.4708046;
1882 cflf[2] = -0.7436830;
1883 cfle[3] = 2.2062068;
1884 cflf[3] = -0.6818757;
1885 cfle[4] = 2.9416091;
1886 cflf[4] = -0.7168990;
1887 cfle[5] = 3.6770115;
1888 cflf[5] = -0.8068384;
1889 cfle[6] = 4.4124136;
1890 cflf[6] = -0.6722584;
1891 cfle[7] = 5.1478162;
1892 cflf[7] = -0.7626385;
1893 cfle[8] = 5.8832183;
1894 cflf[8] = -1.0396487;
1895 cfle[9] = 6.6186204;
1896 cflf[9] = -0.7972314;
1897 cfle[10] = 7.3540230;
1898 cflf[10] = -0.9883416;
1899 cfle[11] = 14.7080460;
1900 cflf[11] = -1.1675659;
1901 cfle[12] = 22.0620689;
1902 cflf[12] = -1.1985949;
1903 cfle[13] = 29.4160919;
1904 cflf[13] = -1.2263466;
1905 cfle[14] = 36.7701149;
1906 cflf[14] = -1.2918345;
1907 cfle[15] = 44.1241379;
1908 cflf[15] = -1.3510833;
1909 cfle[16] = 51.4781609;
1910 cflf[16] = -1.2715496;
1911 cfle[17] = 58.8321838;
1912 cflf[17] = -1.1098027;
1913 cfle[18] = 66.1862030;
1914 cflf[18] = -1.4315782;
1915 cfle[19] = 73.5402298;
1916 cflf[19] = -1.1327956;
1917 cfle[20] = 147.080459;
1918 cflf[20] = -1.6869649;
1919 cfle[21] = 220.620681;
1920 cflf[21] = -2.0239367;
1921 cfle[22] = 294.160919;
1922 cflf[22] = -2.2130392;
1923 cfle[23] = 367.701141;
1924 cflf[23] = -2.3773901;
1925 cfle[24] = 441.241363;
1926 cflf[24] = -2.5326197;
1927 cfle[25] = 514.7816162;
1928 cflf[25] = -2.5292389;
1929 cfle[26] = 588.3218384;
1930 cflf[26] = -2.8230250;
1931 cfle[27] = 661.8620605;
1932 cflf[27] = -2.9502323;
1933 cfle[28] = 735.4022827;
1934 cflf[28] = -3.0774822;
1935 cfle[29] = 1470.8045654;
1936 cflf[29] = -4.2239799;
1937 cfle[30] = 2206.2067871;
1938 cflf[30] = -5.2547927;
1939 cfle[31] = 2941.6091309;
1940 cflf[31] = -6.2353640;
1941 cfle[32] = 3677.0114746;
1942 cflf[32] = -7.1898708;
1943 cfle[33] = 4412.4135742;
1944 cflf[33] = -8.1292381;
1945 cfle[34] = 5147.8159180;
1946 cflf[34] = -9.0594845;
1947 cfle[35] = 5883.2182617;
1948 cflf[35] = -9.9830370;
1949 cfle[36] = 6618.6206055;
1950 cflf[36] = -10.9028034;
1951 cfle[37] = 7354.0229492;
1952 cflf[37] = -11.8188877;
1953 cfle[38] = 7400.0000000;
1954 cflf[38] = -30.0000000;
1955 cfle[39] = 10000000.0000000;
1956 cflf[39] = -30.0000000;
1957 resetBltin( cfle , cflf , true );
1958
1959 /* AKN120 continuum from Brad Peterson's plot
1960 * second vector is F_nu*1E10 as seen from Earth */
1961 tnuakn[0] = 1e-5;
1962 tnuakn[1] = 1.9e-5;
1963 tnuakn[2] = 3.0e-4;
1964 tnuakn[3] = 2.4e-2;
1965 tnuakn[4] = 0.15;
1966 tnuakn[5] = 0.30;
1967 tnuakn[6] = 0.76;
1968 tnuakn[7] = 2.0;
1969 tnuakn[8] = 76.0;
1970 tnuakn[9] = 760.;
1971 tnuakn[10] = 7.4e6;
1972 fnuakn[0] = 1.5e-16;
1973 fnuakn[1] = 1.6e-16;
1974 fnuakn[2] = 1.4e-13;
1975 fnuakn[3] = 8.0e-15;
1976 fnuakn[4] = 1.6e-15;
1977 fnuakn[5] = 1.8e-15;
1978 fnuakn[6] = 7.1e-16;
1979 fnuakn[7] = 7.9e-17;
1980 fnuakn[8] = 1.1e-18;
1981 fnuakn[9] = 7.1e-20;
1982 fnuakn[10] = 1.3e-24;
1983 resetBltin( tnuakn , fnuakn , false );
1984
1985 /* interstellar radiation field from Black 1987, "Interstellar Processes"
1986 * table of log nu, log nu*fnu taken from his figure 2 */
1987 /* >>chng 99 jun 14 energy range lowered to 1e-8 ryd */
1988 tnuism[0] = 6.00;
1989 /*tnuism[0] = 9.00;*/
1990 tnuism[1] = 10.72;
1991 tnuism[2] = 11.00;
1992 tnuism[3] = 11.23;
1993 tnuism[4] = 11.47;
1994 tnuism[5] = 11.55;
1995 tnuism[6] = 11.85;
1996 tnuism[7] = 12.26;
1997 tnuism[8] = 12.54;
1998 tnuism[9] = 12.71;
1999 tnuism[10] = 13.10;
2000 tnuism[11] = 13.64;
2001 tnuism[12] = 14.14;
2002 tnuism[13] = 14.38;
2003 tnuism[14] = 14.63;
2004 tnuism[15] = 14.93;
2005 tnuism[16] = 15.08;
2006 tnuism[17] = 15.36;
2007 tnuism[18] = 15.43;
2008 tnuism[19] = 16.25;
2009 tnuism[20] = 17.09;
2010 tnuism[21] = 18.00;
2011 tnuism[22] = 23.00;
2012 /* these are log nu*Fnu */
2013 fnuism[0] = -16.708;
2014 /*fnuism[0] = -7.97;*/
2015 fnuism[1] = -2.96;
2016 fnuism[2] = -2.47;
2017 fnuism[3] = -2.09;
2018 fnuism[4] = -2.11;
2019 fnuism[5] = -2.34;
2020 fnuism[6] = -3.66;
2021 fnuism[7] = -2.72;
2022 fnuism[8] = -2.45;
2023 fnuism[9] = -2.57;
2024 fnuism[10] = -3.85;
2025 fnuism[11] = -3.34;
2026 fnuism[12] = -2.30;
2027 fnuism[13] = -1.79;
2028 fnuism[14] = -1.79;
2029 fnuism[15] = -2.34;
2030 fnuism[16] = -2.72;
2031 fnuism[17] = -2.55;
2032 fnuism[18] = -2.62;
2033 fnuism[19] = -5.68;
2034 fnuism[20] = -6.45;
2035 fnuism[21] = -6.30;
2036 fnuism[22] = -11.3;
2037
2038 return;
2039}
2040
2041/*lines_table invoked by table lines command, check if we can find all lines in a given list */
2043{
2044 DEBUG_ENTRY( "lines_table()" );
2045
2046 if( chLINE_LIST.empty() )
2047 return 0;
2048
2049 vector<char*> chLabel;
2050 vector<realnum> wl;
2051
2052 long nLINE_TABLE = cdGetLineList( chLINE_LIST.c_str(), chLabel, wl );
2053
2054 // the check if the file exists has already been done by the parser
2055 if( nLINE_TABLE == 0 )
2056 return 0;
2057
2058 fprintf( ioQQQ , "lines_table checking lines within data table %s\n", chLINE_LIST.c_str() );
2059 long miss = 0;
2060
2061 /* \todo 2 DOS carriage return on linux screws this up. Can we overlook the CR? Skip in cdgetlinelist? */
2062 for( long n=0; n < nLINE_TABLE; ++n )
2063 {
2064 double relative , absolute;
2065 if( (cdLine( chLabel[n], wl[n] , &relative , &absolute )) <= 0 )
2066 {
2067 ++miss;
2068 fprintf(ioQQQ,"lines_table in parse_table.cpp did not find line %4s ",chLabel[n]);
2069 prt_wl(ioQQQ,wl[n]);
2070 fprintf(ioQQQ,"\n");
2071 }
2072 }
2073 if( miss )
2074 {
2075 /* this is so that we pick up problem automatically */
2076 fprintf( ioQQQ , " BOTCHED MONITORS!!! Botched Monitors!!! lines_table could not find a total of %li lines\n\n", miss );
2077 }
2078 else
2079 {
2080 fprintf( ioQQQ , "lines_table found all lines\n\n" );
2081 }
2082
2083 // deallocate the memory allocated in cdGetLineList()
2084 for( unsigned j=0; j < chLabel.size(); ++j )
2085 delete[] chLabel[j];
2086 chLabel.clear();
2087
2088 return miss;
2089}
2090
2091/*ReadTable called by TABLE READ to read in continuum from SAVE TRANSMITTED CONTINUUM */
2092STATIC void ReadTable(const char *fnam)
2093{
2094 char chLine[INPUT_LINE_LENGTH];
2095 long int i;
2096 FILE *io;
2097
2098 DEBUG_ENTRY( "ReadTable()" );
2099
2100 /* make unused array has valid zeros */
2101 for( i=0; i < NCELL; i++ )
2102 {
2103 rfield.tFluxLog[rfield.nShape][i] = -70.;
2104 }
2105
2106 io = open_data( fnam, "r", AS_LOCAL_ONLY );
2107
2108 string unit;
2109 char *last;
2110
2111 /* read in first line of header and parse for units, if present */
2112 if( read_whole_line( chLine , (int)sizeof(chLine) , io )==NULL )
2113 {
2114 fprintf( ioQQQ, " the input continuum file has been truncated.\n" );
2115 goto error;
2116 }
2117
2118 unit = "Ryd"; // default
2119 last = strchr_s(chLine,'\t');
2120 if (last)
2121 {
2122 *last = '\0';
2123 char *first = strrchr(chLine,'/');
2124 if (first)
2125 {
2126 unit = string(first+1);
2127 }
2128 *last = '\t';
2129 };
2130
2131 /* line 2: empty comment line */
2132 if( read_whole_line( chLine , (int)sizeof(chLine) , io )==NULL )
2133 {
2134 fprintf( ioQQQ, " the input continuum file has been truncated.\n" );
2135 goto error;
2136 }
2137
2138 /* line 3: the version number */
2139 if( read_whole_line( chLine , (int)sizeof(chLine) , io )!=NULL )
2140 {
2141 long version;
2142 sscanf( chLine, "%ld", &version );
2143 if( version != VERSION_TRNCON )
2144 {
2145 fprintf( ioQQQ,
2146 " the input continuum file has the wrong version number, I expected %li and found %li.\n",
2147 VERSION_TRNCON, version);
2148 goto error;
2149 }
2150 }
2151 else
2152 {
2153 fprintf( ioQQQ, " the input continuum file has been truncated.\n" );
2154 goto error;
2155 }
2156
2157 char md5sum[NMD5];
2158 long nflux;
2159 double mesh_lo, mesh_hi;
2160 union {
2161 double x;
2162 uint32 i[2];
2163 } u;
2164
2165 /* line 4: the MD5 checksum */
2166 if( read_whole_line( chLine , (int)sizeof(chLine) , io )!=NULL )
2167 {
2168 strncpy( md5sum, chLine, NMD5 );
2169 }
2170 else
2171 {
2172 fprintf( ioQQQ, " the input continuum file has been truncated.\n" );
2173 goto error;
2174 }
2175
2176 /* line 5: the lower limit of the energy grid */
2177 if( read_whole_line( chLine , (int)sizeof(chLine) , io )!=NULL )
2178 {
2179 if( cpu.i().big_endian() )
2180 sscanf( chLine, "%x %x", &u.i[0], &u.i[1] );
2181 else
2182 sscanf( chLine, "%x %x", &u.i[1], &u.i[0] );
2183 mesh_lo = u.x;
2184 }
2185 else
2186 {
2187 fprintf( ioQQQ, " the input continuum file has been truncated.\n" );
2188 goto error;
2189 }
2190
2191 /* line 6: the upper limit of the energy grid */
2192 if( read_whole_line( chLine , (int)sizeof(chLine) , io )!=NULL )
2193 {
2194 if( cpu.i().big_endian() )
2195 sscanf( chLine, "%x %x", &u.i[0], &u.i[1] );
2196 else
2197 sscanf( chLine, "%x %x", &u.i[1], &u.i[0] );
2198 mesh_hi = u.x;
2199 }
2200 else
2201 {
2202 fprintf( ioQQQ, " the input continuum file has been truncated.\n" );
2203 goto error;
2204 }
2205
2206 if( strncmp( md5sum, continuum.mesh_md5sum.c_str(), NMD5 ) != 0 ||
2207 !fp_equal( mesh_lo, double(rfield.emm) ) ||
2208 !fp_equal( mesh_hi, double(rfield.egamry) ) )
2209 {
2210 fprintf( ioQQQ, " the input continuum file has an incompatible energy grid.\n" );
2211 goto error;
2212 }
2213
2214 /* line 7: the energy grid resolution scale factor */
2215 if( read_whole_line( chLine , (int)sizeof(chLine) , io )!=NULL )
2216 {
2217 if( cpu.i().big_endian() )
2218 sscanf( chLine, "%x %x", &u.i[0], &u.i[1] );
2219 else
2220 sscanf( chLine, "%x %x", &u.i[1], &u.i[0] );
2221 rfield.RSFCheck[rfield.nShape] = u.x;
2222 }
2223 else
2224 {
2225 fprintf( ioQQQ, " the input continuum file has been truncated.\n" );
2226 goto error;
2227 }
2228
2229 /* line 8: the number of frequency grid points contained in the file */
2230 if( read_whole_line( chLine , (int)sizeof(chLine) , io )!=NULL )
2231 {
2232 sscanf( chLine, "%ld", &nflux );
2233 }
2234 else
2235 {
2236 fprintf( ioQQQ, " the input continuum file has been truncated.\n" );
2237 goto error;
2238 }
2239
2240 /* line 9: empty comment line */
2241 if( read_whole_line( chLine , (int)sizeof(chLine) , io )==NULL )
2242 {
2243 fprintf( ioQQQ, " the input continuum file has been truncated.\n" );
2244 goto error;
2245 }
2246
2247 /* now read in the file of numbers */
2248 i = 0;
2249 /* keep reading until we hit eol or run out of room in the array */
2250 while( (read_whole_line(chLine, (int)sizeof(chLine),io)!=NULL) && (i<NCELL) )
2251 {
2252 double help[2];
2253 sscanf( chLine, "%lf%lf ", &help[0], &help[1] );
2254 rfield.tNu[rfield.nShape][i].set(help[0],unit.c_str());
2255 // Convert to standard FluxLog units
2256 if (help[1] > 0.0)
2257 {
2258 rfield.tFluxLog[rfield.nShape][i] = (realnum) log10(help[1]/
2259 rfield.tNu[rfield.nShape][i].Ryd());
2260 }
2261 ++i;
2262 }
2263 /* put pointer at last good value */
2264 rfield.ncont[rfield.nShape] = i;
2265
2266 /* check that sane number of values entered */
2267 if( rfield.ncont[rfield.nShape] != nflux )
2268 {
2269 fprintf( ioQQQ, " the input continuum file has the wrong number of points: %ld\n",
2270 rfield.ncont[rfield.nShape] );
2271 goto error;
2272 }
2273
2274 fclose(io);
2275 return;
2276
2277error:
2278 fprintf( ioQQQ, " please recreate this file using the SAVE TRANSMITTED CONTINUUM command.\n" );
2280}
2281
2282#if defined(__HP_aCC)
2283#pragma OPTIMIZE OFF
2284#pragma OPTIMIZE ON
2285#endif
FILE * ioQQQ
Definition cddefines.cpp:7
const int FILENAME_PATH_LENGTH_2
Definition cddefines.h:249
#define ASSERT(exp)
Definition cddefines.h:578
#define STATIC
Definition cddefines.h:97
#define MALLOC(exp)
Definition cddefines.h:501
const int INPUT_LINE_LENGTH
Definition cddefines.h:254
#define EXIT_SUCCESS
Definition cddefines.h:138
#define EXIT_FAILURE
Definition cddefines.h:140
double FFmtRead(const char *chCard, long int *ipnt, long int last, bool *lgEOL)
Definition service.cpp:381
#define cdEXIT(FAIL)
Definition cddefines.h:434
float realnum
Definition cddefines.h:103
const char * strchr_s(const char *s, int c)
Definition cddefines.h:1439
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition cddefines.h:812
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
Definition service.cpp:70
NORETURN void TotalInsanity(void)
Definition service.cpp:886
sys_float SDIV(sys_float x)
Definition cddefines.h:952
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
static bool lgCalled
Definition cddrive.cpp:425
long int cdLine(const char *chLabel, realnum wavelength, double *relint, double *absint)
Definition cddrive.cpp:1228
long int cdGetLineList(const char chFile[], vector< char * > &chLabels, vector< realnum > &wl)
double FFmtRead(void)
Definition parser.cpp:353
bool nMatch(const char *chKey) const
Definition parser.h:135
bool nMatchErase(const char *chKey)
Definition parser.h:158
bool lgEOL(void) const
Definition parser.h:98
int GetQuote(char *chLabel, bool lgABORT)
Definition parser.h:209
long int m_nqh
Definition parser.h:41
t_continuum continuum
Definition continuum.cpp:5
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition cpu.cpp:625
static t_cpu cpu
Definition cpu.h:355
UNUSED const realnum BIGFLOAT
Definition cpu.h:189
@ AS_LOCAL_ONLY
Definition cpu.h:208
@ AS_LOCAL_DATA
Definition cpu.h:208
@ AS_LOCAL_DATA_TRY
Definition cpu.h:207
const realnum SMALLFLOAT
Definition cpu.h:191
t_input input
Definition input.cpp:12
@ SYMMETRIC
Definition rfield.h:30
t_optimize optimize
Definition optimize.cpp:5
const long LIMEXT
Definition optimize.h:60
static const double tnuHM96[NHM96]
static double cfle[NCFL]
static const int NDRAINE
static double tnurbn[NRUBIN]
STATIC void resetBltin(double *tnu, double *fluxlog, bool lgLog)
static double fnuism[NISM]
static double tnuagn[NAGN]
static double tnuism[NISM]
static const int NCFL
static double fnuakn[NKN120]
static string chLINE_LIST
static const double fnuHM96[NHM96]
static const int NISM
int lines_table()
static double tsldrn[NDRAINE]
STATIC void ReadTable(const char *fnam)
static double tslagn[NAGN]
static double tnuakn[NKN120]
void ParseTable(Parser &p)
static const int NKN120
static const int NCRAB
STATIC void ZeroContin(void)
STATIC void read_hm05(const char chFile[], double **tnuHM, double **fnuHM, long int *nhm, double redshift)
static const int NRUBIN
static const int NHM96
static double tnudrn[NDRAINE]
static double cflf[NCFL]
static const int NAGN
static double fnucrb[NCRAB]
static double tnucrb[NCRAB]
static double fnurbn[NRUBIN]
UNUSED const double FR1RYD
Definition physconst.h:195
UNUSED const double PI4
Definition physconst.h:35
UNUSED const double HIONPOT
Definition physconst.h:119
UNUSED const double EVRYD
Definition physconst.h:189
UNUSED const double RYDLAM
Definition physconst.h:176
void prt_wl(FILE *ioOUT, realnum wl)
Definition prt.cpp:13
t_radius radius
Definition radius.cpp:5
t_rfield rfield
Definition rfield.cpp:8
const int NCELL
Definition rfield.h:21
const int LIMSPC
Definition rfield.h:18
static const long VERSION_TRNCON
Definition save.h:15
long RauchInterpolateHydr(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
Definition stars.cpp:1201
long Kurucz79Interpolate(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
Definition stars.cpp:799
long WMBASICInterpolate(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
Definition stars.cpp:1788
long AtlasInterpolate(double val[], long *nval, long *ndim, const char *chMetalicity, const char *chODFNew, bool lgList, double *Tlow, double *Thigh)
Definition stars.cpp:515
long WernerInterpolate(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
Definition stars.cpp:1708
long RauchInterpolateHNi(double val[], long *nval, long *ndim, bool lgHalo, bool lgList, double *Tlow, double *Thigh)
Definition stars.cpp:1107
long CoStarInterpolate(double val[], long *nval, long *ndim, IntMode imode, bool lgHalo, bool lgList, double *val0_lo, double *val0_hi)
Definition stars.cpp:623
long GridInterpolate(double val[], long *nval, long *ndim, const char *FileName, bool lgList, double *Tlow, double *Thigh)
Definition stars.cpp:738
long RauchInterpolateHpHe(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
Definition stars.cpp:1261
long RauchInterpolateHelium(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
Definition stars.cpp:1231
long TlustyInterpolate(double val[], long *nval, long *ndim, tl_grid tlg, const char *chMetalicity, bool lgList, double *Tlow, double *Thigh)
Definition stars.cpp:1590
long RauchInterpolateHCa(double val[], long *nval, long *ndim, bool lgHalo, bool lgList, double *Tlow, double *Thigh)
Definition stars.cpp:1073
void AtmospheresAvail(void)
Definition stars.cpp:202
long RauchInterpolateCOWD(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
Definition stars.cpp:1171
long MihalasInterpolate(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
Definition stars.cpp:849
long RauchInterpolatePG1159(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
Definition stars.cpp:1141
IntMode
Definition stars.h:16
@ IM_COSTAR_MZAMS_AGE
Definition stars.h:18
@ IM_COSTAR_AGE_MZAMS
Definition stars.h:18
@ IM_COSTAR_TEFF_LOGG
Definition stars.h:18
@ IM_ILLEGAL_MODE
Definition stars.h:17
@ IM_COSTAR_TEFF_MODID
Definition stars.h:17
@ TL_OBSTAR
Definition stars.h:22
@ TL_OSTAR
Definition stars.h:22
@ TL_BSTAR
Definition stars.h:22
#define MDIM
Definition stars.h:9
static const unsigned int NMD5
Definition thirdparty.h:384
t_trace trace
Definition trace.cpp:5