cloudy  trunk
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 */
28 STATIC void ReadTable(const char * fnam);
29 
30 static string chLINE_LIST;
31 
32 /* tables of various built in continue */
33 /* Davidson 1985 version of crab nebula SED */
34 static const int NCRAB = 10;
35 static double tnucrb[NCRAB],
36  fnucrb[NCRAB];
37 
38 /* Bob Rubin's corrected theta 1 Ori C continuum */
39 static const int NRUBIN = 56;
40 static 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 */
53 static const int NCFL = 40;
54 static double cfle[NCFL],
55  cflf[NCFL];
56 
57 /* Brad Peterson's AKN 120 continuum */
58 static const int NKN120 = 11;
59 static double tnuakn[NKN120],
60  fnuakn[NKN120];
61 
62 /* Black's ISM continuum, with He hole filled in */
63 static const int NISM = 23;
64 static double tnuism[NISM],
65  fnuism[NISM];
66 
67 /* z=2 background,
68  * >>refer continuum background Haardt, Francesco, & Madau, Piero, 1996,
69  * >>refercon ApJ, 461, 20 */
70 static const int NHM96 = 14;
71 /* log energy in Ryd */
72 static 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,
75 0.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})*/
78 static 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 */
82 static const int NAGN = 8;
83 static double tnuagn[NAGN],
84  tslagn[NAGN];
85 
86 /* table Draine ISM continuum */
87 static const int NDRAINE = 15;
88 static double tnudrn[NDRAINE] , tsldrn[NDRAINE];
89 
90 /* routine that stores values for above vectors */
91 STATIC 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 */
98 STATIC 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 */
133 STATIC 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 );
174  cdEXIT(EXIT_FAILURE);
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 )
205  TotalInsanity();
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");
223  cdEXIT(EXIT_FAILURE);
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 )
241  TotalInsanity();
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");
283  cdEXIT(EXIT_FAILURE);
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");
309  cdEXIT(EXIT_FAILURE);
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 
340 void ParseTable(Parser &p)
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 );
379  cdEXIT(EXIT_FAILURE);
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" );
428  cdEXIT(EXIT_FAILURE);
429  }
430  }
431 
432  if( ConBreak == 0. )
433  {
434  fprintf( ioQQQ, " The break must be greater than 0.2 Ryd.\n Sorry.\n" );
435  cdEXIT(EXIT_FAILURE);
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" );
453  cdEXIT(EXIT_FAILURE);
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() );
460  cdEXIT(EXIT_FAILURE);
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() );
467  cdEXIT(EXIT_FAILURE);
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 );
605  cdEXIT(EXIT_FAILURE);
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" );
615  cdEXIT(EXIT_FAILURE);
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");
677  cdEXIT(EXIT_FAILURE);
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 );
716  cdEXIT(EXIT_FAILURE);
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" );
726  cdEXIT(EXIT_FAILURE);
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 );
813  cdEXIT(EXIT_FAILURE);
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" );
823  cdEXIT(EXIT_FAILURE);
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 );
887  cdEXIT(EXIT_FAILURE);
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" );
897  cdEXIT(EXIT_FAILURE);
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");
948  cdEXIT(EXIT_FAILURE);
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");
971  cdEXIT(EXIT_FAILURE);
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 );
1080  cdEXIT(EXIT_FAILURE);
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");
1107  cdEXIT(EXIT_FAILURE);
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" );
1129  cdEXIT(EXIT_FAILURE);
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  {
1228  AtmospheresAvail();
1229  cdEXIT(EXIT_SUCCESS);
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" );
1243  cdEXIT(EXIT_FAILURE);
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");
1260  cdEXIT(EXIT_FAILURE );
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" );
1314  cdEXIT(EXIT_FAILURE);
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" );
1323  cdEXIT(EXIT_FAILURE);
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" );
1347  cdEXIT(EXIT_FAILURE);
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" );
1368  cdEXIT(EXIT_FAILURE);
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" );
1430  cdEXIT(EXIT_FAILURE);
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" );
1452  cdEXIT(EXIT_FAILURE);
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" );
1667  cdEXIT(EXIT_FAILURE);
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 */
1775 STATIC void ZeroContin(void)
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 */
2042 int lines_table()
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 */
2092 STATIC 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 
2277 error:
2278  fprintf( ioQQQ, " please recreate this file using the SAVE TRANSMITTED CONTINUUM command.\n" );
2279  cdEXIT(EXIT_FAILURE);
2280 }
2281 
2282 #if defined(__HP_aCC)
2283 #pragma OPTIMIZE OFF
2284 #pragma OPTIMIZE ON
2285 #endif
RauchInterpolateHydr
long RauchInterpolateHydr(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
Definition: stars.cpp:1201
tnucrb
static double tnucrb[NCRAB]
Definition: parse_table.cpp:35
TL_BSTAR
@ TL_BSTAR
Definition: stars.h:22
FR1RYD
const UNUSED double FR1RYD
Definition: physconst.h:195
t_rfield::totpow
double totpow[LIMSPC]
Definition: rfield.h:300
Parser::nMatch
bool nMatch(const char *chKey) const
Definition: parser.h:135
ParseTable
void ParseTable(Parser &p)
Definition: parse_table.cpp:340
Kurucz79Interpolate
long Kurucz79Interpolate(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
Definition: stars.cpp:799
open_data
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition: cpu.cpp:625
lines.h
t_optimize::vincr
realnum vincr[LIMPAR]
Definition: optimize.h:191
Parser::FFmtRead
double FFmtRead(void)
Definition: parser.cpp:353
RauchInterpolateHCa
long RauchInterpolateHCa(double val[], long *nval, long *ndim, bool lgHalo, bool lgList, double *Tlow, double *Thigh)
Definition: stars.cpp:1073
t_optimize::nparm
long int nparm
Definition: optimize.h:201
cfle
static double cfle[NCFL]
Definition: parse_table.cpp:54
t_input::nRead
long int nRead
Definition: input.h:49
NISM
static const int NISM
Definition: parse_table.cpp:63
rfield
t_rfield rfield
Definition: rfield.cpp:8
ioQQQ
FILE * ioQQQ
Definition: cddefines.cpp:7
t_rfield::tNu
vector< Energy > tNu[LIMSPC]
Definition: rfield.h:330
IM_COSTAR_TEFF_LOGG
@ IM_COSTAR_TEFF_LOGG
Definition: stars.h:18
FFmtRead
double FFmtRead(const char *chCard, long int *ipnt, long int last, bool *lgEOL)
Definition: service.cpp:381
realnum
float realnum
Definition: cddefines.h:103
rfield.h
tslagn
static double tslagn[NAGN]
Definition: parse_table.cpp:84
STATIC
#define STATIC
Definition: cddefines.h:97
AS_LOCAL_DATA_TRY
@ AS_LOCAL_DATA_TRY
Definition: cpu.h:207
NCFL
static const int NCFL
Definition: parse_table.cpp:53
HIONPOT
const UNUSED double HIONPOT
Definition: physconst.h:119
AS_LOCAL_ONLY
@ AS_LOCAL_ONLY
Definition: cpu.h:208
fnuism
static double fnuism[NISM]
Definition: parse_table.cpp:65
t_optimize::lgVarOn
bool lgVarOn
Definition: optimize.h:203
thirdparty.h
fnuHM96
static const double fnuHM96[NHM96]
Definition: parse_table.cpp:78
LIMEXT
const long LIMEXT
Definition: optimize.h:60
cpu
static t_cpu cpu
Definition: cpu.h:355
t_rfield::RSFCheck
double RSFCheck[LIMSPC]
Definition: rfield.h:339
trace.h
t_rfield::chSpType
char chSpType[LIMSPC][6]
Definition: rfield.h:353
stars.h
WMBASICInterpolate
long WMBASICInterpolate(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
Definition: stars.cpp:1788
LIMSPC
const int LIMSPC
Definition: rfield.h:18
SDIV
sys_float SDIV(sys_float x)
Definition: cddefines.h:952
t_rfield::lgContMalloc
bool lgContMalloc[LIMSPC]
Definition: rfield.h:343
MDIM
#define MDIM
Definition: stars.h:9
Parser::GetQuote
int GetQuote(char *chLabel, bool lgABORT)
Definition: parser.h:209
NHM96
static const int NHM96
Definition: parse_table.cpp:70
cflf
static double cflf[NCFL]
Definition: parse_table.cpp:55
input
t_input input
Definition: input.cpp:12
ASSERT
#define ASSERT(exp)
Definition: cddefines.h:578
PI4
const UNUSED double PI4
Definition: physconst.h:35
t_continuum::mesh_md5sum
string mesh_md5sum
Definition: continuum.h:127
MihalasInterpolate
long MihalasInterpolate(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
Definition: stars.cpp:849
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
EXIT_SUCCESS
#define EXIT_SUCCESS
Definition: cddefines.h:138
t_rfield::tslop
vector< realnum > tslop[LIMSPC]
Definition: rfield.h:331
NCRAB
static const int NCRAB
Definition: parse_table.cpp:34
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
t_rfield::lgMustBlockHIon
bool lgMustBlockHIon
Definition: rfield.h:110
NMD5
static const unsigned int NMD5
Definition: thirdparty.h:384
cddrive.h
t_radius::lgRadiusKnown
bool lgRadiusKnown
Definition: radius.h:116
t_trace::lgConBug
bool lgConBug
Definition: trace.h:100
ZeroContin
STATIC void ZeroContin(void)
Definition: parse_table.cpp:1775
radius
t_radius radius
Definition: radius.cpp:5
IM_COSTAR_TEFF_MODID
@ IM_COSTAR_TEFF_MODID
Definition: stars.h:17
tsldrn
static double tsldrn[NDRAINE]
Definition: parse_table.cpp:88
optimize
t_optimize optimize
Definition: optimize.cpp:5
t_cpu::i
t_cpu_i & i()
Definition: cpu.h:347
t_optimize::vparm
realnum vparm[LIMEXT][LIMPAR]
Definition: optimize.h:188
t_optimize::chVarFmt
char chVarFmt[LIMPAR][FILENAME_PATH_LENGTH_2]
Definition: optimize.h:263
EXIT_FAILURE
#define EXIT_FAILURE
Definition: cddefines.h:140
t_rfield::egamry
realnum egamry
Definition: rfield.h:52
t_rfield::chRSpec
char chRSpec[LIMSPC][5]
Definition: rfield.h:352
t_cpu_i::big_endian
bool big_endian() const
Definition: cpu.h:287
Parser
Definition: parser.h:31
lines_table
int lines_table()
Definition: parse_table.cpp:2042
tnuakn
static double tnuakn[NKN120]
Definition: parse_table.cpp:59
trace
t_trace trace
Definition: trace.cpp:5
cddefines.h
IM_COSTAR_AGE_MZAMS
@ IM_COSTAR_AGE_MZAMS
Definition: stars.h:18
lgCalled
static bool lgCalled
Definition: cddrive.cpp:425
t_radius::Radius
double Radius
Definition: radius.h:25
ReadTable
STATIC void ReadTable(const char *fnam)
Definition: parse_table.cpp:2092
Illuminate::SYMMETRIC
@ SYMMETRIC
Definition: rfield.h:30
t_rfield::tFluxLog
vector< realnum > tFluxLog[LIMSPC]
Definition: rfield.h:333
optimize.h
TotalInsanity
NORETURN void TotalInsanity(void)
Definition: service.cpp:886
t_radius::rdfalt
double rdfalt
Definition: radius.h:124
tnuHM96
static const double tnuHM96[NHM96]
Definition: parse_table.cpp:72
t_optimize::nvarxt
long int nvarxt[LIMPAR]
Definition: optimize.h:194
t_rfield::Illumination
Illuminate::IlluminationType Illumination[LIMSPC]
Definition: rfield.h:316
t_rfield::ncont
long ncont[LIMSPC]
Definition: rfield.h:335
radius.h
Parser::nMatchErase
bool nMatchErase(const char *chKey)
Definition: parser.h:158
MALLOC
#define MALLOC(exp)
Definition: cddefines.h:501
tnudrn
static double tnudrn[NDRAINE]
Definition: parse_table.cpp:88
cdGetLineList
long int cdGetLineList(const char chFile[], vector< char * > &chLabels, vector< realnum > &wl)
Definition: cdgetlinelist.cpp:10
NCELL
const int NCELL
Definition: rfield.h:21
RauchInterpolateHpHe
long RauchInterpolateHpHe(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
Definition: stars.cpp:1261
tnurbn
static double tnurbn[NRUBIN]
Definition: parse_table.cpp:40
RYDLAM
const UNUSED double RYDLAM
Definition: physconst.h:176
RauchInterpolateCOWD
long RauchInterpolateCOWD(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
Definition: stars.cpp:1171
tnuagn
static double tnuagn[NAGN]
Definition: parse_table.cpp:83
cdEXIT
#define cdEXIT(FAIL)
Definition: cddefines.h:434
Parser::lgEOL
bool lgEOL(void) const
Definition: parser.h:98
BIGFLOAT
const UNUSED realnum BIGFLOAT
Definition: cpu.h:189
cdLine
long int cdLine(const char *chLabel, realnum wavelength, double *relint, double *absint)
Definition: cddrive.cpp:1228
save.h
t_rfield::nupper
long int nupper
Definition: rfield.h:46
TL_OBSTAR
@ TL_OBSTAR
Definition: stars.h:22
read_hm05
STATIC void read_hm05(const char chFile[], double **tnuHM, double **fnuHM, long int *nhm, double redshift)
Definition: parse_table.cpp:133
IM_COSTAR_MZAMS_AGE
@ IM_COSTAR_MZAMS_AGE
Definition: stars.h:18
NRUBIN
static const int NRUBIN
Definition: parse_table.cpp:39
prt.h
FILENAME_PATH_LENGTH_2
const int FILENAME_PATH_LENGTH_2
Definition: cddefines.h:249
t_optimize::varang
realnum varang[LIMPAR][2]
Definition: optimize.h:198
INPUT_LINE_LENGTH
const int INPUT_LINE_LENGTH
Definition: cddefines.h:254
t_rfield::lgTimeVary
bool lgTimeVary[LIMSPC]
Definition: rfield.h:306
t_rfield::lgBeamed
bool lgBeamed[LIMSPC]
Definition: rfield.h:310
t_rfield::range
double range[LIMSPC][2]
Definition: rfield.h:347
TL_OSTAR
@ TL_OSTAR
Definition: stars.h:22
fp_equal
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition: cddefines.h:812
parser.h
physconst.h
tnuism
static double tnuism[NISM]
Definition: parse_table.cpp:64
fnuakn
static double fnuakn[NKN120]
Definition: parse_table.cpp:60
NAGN
static const int NAGN
Definition: parse_table.cpp:82
t_optimize::nvfpnt
long int nvfpnt[LIMPAR]
Definition: optimize.h:195
AtmospheresAvail
void AtmospheresAvail(void)
Definition: stars.cpp:202
NDRAINE
static const int NDRAINE
Definition: parse_table.cpp:87
VERSION_TRNCON
static const long VERSION_TRNCON
Definition: save.h:15
fnucrb
static double fnucrb[NCRAB]
Definition: parse_table.cpp:36
prt_wl
void prt_wl(FILE *ioOUT, realnum wl)
Definition: prt.cpp:13
Parser::m_nqh
long int m_nqh
Definition: parser.h:41
read_whole_line
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
Definition: service.cpp:70
NKN120
static const int NKN120
Definition: parse_table.cpp:58
t_rfield::nShape
long int nShape
Definition: rfield.h:322
RauchInterpolateHelium
long RauchInterpolateHelium(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
Definition: stars.cpp:1231
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
RauchInterpolatePG1159
long RauchInterpolatePG1159(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
Definition: stars.cpp:1141
continuum
t_continuum continuum
Definition: continuum.cpp:5
fnurbn
static double fnurbn[NRUBIN]
Definition: parse_table.cpp:45
t_rfield::emm
realnum emm
Definition: rfield.h:49
t_rfield::chSpNorm
char chSpNorm[LIMSPC][5]
Definition: rfield.h:351
resetBltin
STATIC void resetBltin(double *tnu, double *fluxlog, bool lgLog)
Definition: parse_table.cpp:98
IntMode
IntMode
Definition: stars.h:16
IM_ILLEGAL_MODE
@ IM_ILLEGAL_MODE
Definition: stars.h:17
continuum.h
WernerInterpolate
long WernerInterpolate(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
Definition: stars.cpp:1708
EVRYD
const UNUSED double EVRYD
Definition: physconst.h:189
AS_LOCAL_DATA
@ AS_LOCAL_DATA
Definition: cpu.h:208
chLINE_LIST
static string chLINE_LIST
Definition: parse_table.cpp:30
input.h
GridInterpolate
long GridInterpolate(double val[], long *nval, long *ndim, const char *FileName, bool lgList, double *Tlow, double *Thigh)
Definition: stars.cpp:738
DEBUG_ENTRY
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:684
t_trace::lgTrace
bool lgTrace
Definition: trace.h:12
SMALLFLOAT
const realnum SMALLFLOAT
Definition: cpu.h:191
strchr_s
const char * strchr_s(const char *s, int c)
Definition: cddefines.h:1439
RauchInterpolateHNi
long RauchInterpolateHNi(double val[], long *nval, long *ndim, bool lgHalo, bool lgList, double *Tlow, double *Thigh)
Definition: stars.cpp:1107
Generated by   doxygen 1.8.17