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

Generated for cloudy by doxygen 1.8.1.1