cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
parse_set.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 /*ParseSet scan parameters off SET command */
4 #include "cddefines.h"
5 #include "physconst.h"
6 #include "geometry.h"
7 #include "input.h"
8 #include "prt.h"
9 #include "mole.h"
10 #include "rt.h"
11 #include "phycon.h"
12 #include "optimize.h"
13 #include "hcmap.h"
14 #include "hmi.h"
15 #include "dense.h"
16 #include "h2.h"
17 #include "iterations.h"
18 #include "conv.h"
19 #include "secondaries.h"
20 #include "rfield.h"
21 #include "ionbal.h"
22 #include "numderiv.h"
23 #include "dynamics.h"
24 #include "iso.h"
25 #include "punch.h"
26 #include "stopcalc.h"
27 #include "opacity.h"
28 #include "hydrogenic.h"
29 #include "peimbt.h"
30 #include "radius.h"
31 #include "atmdat.h"
32 #include "continuum.h"
33 #include "grains.h"
34 #include "grainvar.h"
35 #include "parse.h"
36 #include "lines.h"
37 #include "assertresults.h"
38 #include "thirdparty.h"
39 
40 void ParseSet(char *chCard )
41 {
42  bool lgEOL;
43  long int i,
44  ip;
45  char chString_quotes_lowercase[INPUT_LINE_LENGTH];
46  bool lgQuotesFound;
47 
48  DEBUG_ENTRY( "ParseSet()" );
49 
50  /* first check for any strings within quotes - this will get the string
51  * and blank it out, so that these are not confused with keywords. if
52  * double quotes not present routine returns unity, zero if found*/
53  lgQuotesFound = true;
54  if( GetQuote( chString_quotes_lowercase , chCard , false ) )
55  lgQuotesFound = false;
56 
57  /* commands to set certain variables, "SET XXX=" */
58  if( nMatch("ASSE",chCard) && nMatch("ABOR",chCard) )
59  {
60  /* set assert abort command, to tell the code to raise SIGABRT when an
61  * assert is thrown - this can be caught by the debugger. The keyword
62  * _FPE is accepted for backward compatibility */
63  cpu.setAssertAbort( true );
64  }
65 
66  else if( nMatch("ASSE",chCard) &&nMatch("SCIE",chCard) )
67  {
68  /* print asserts using scientific notation. Useful when an asserted value is
69  * less than 10^-4 of the normalization line. */
70  lgPrtSciNot = true;
71  }
72 
73  else if( nMatch("ATOM",chCard) && nMatch( "DATA", chCard ) )
74  {
75  /* set atomic data */
76  long int nelem , ion;
77  bool lgAs=false , lgCS=false , lgDR=false;
78  const char *chMole;
79 
80  /* As? */
81  if( nMatch( " AS " , chCard ) )
82  lgAs = true;
83  /* DR - dielectronic recombination */
84  else if( nMatch( " DR " , chCard ) )
85  lgDR = true;
86  /* lgCS collision strengths */
87  else if( nMatch( " CS " , chCard ) )
88  lgCS = true;
89  else
90  {
91  /* no keyword */
92  fprintf( ioQQQ, " One of the keywords AS DR or CS must appear to change"
93  " the transition probabilities, dielectronic recombination rate,"
94  " or collision strength.\n Sorry.\n" );
95  cdEXIT(EXIT_FAILURE);
96  }
97 
98  /* ion or molecule? */
99  if( nMatch(" ION",chCard) )
100  {
101  /* set atomic data - must have an element name and possibly an
102  * ionization stage - nelem is atomic number on C scale */
103  if( (nelem = GetElem(chCard))<0 )
104  {
105  fprintf( ioQQQ, " An element name must appear on this line\n Sorry.\n" );
106  cdEXIT(EXIT_FAILURE);
107  }
108  i = 5;
109  /* make sure that the ionization stage is ok if it is specified
110  * this is entered on physics scale, 1 is atom, 2 first ion */
111  ion = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
112  if( !lgEOL && (ion< 1 || ion>nelem+1) )
113  {
114  fprintf( ioQQQ, " The ionization stage is not valid\n Sorry.\n" );
115  cdEXIT(EXIT_FAILURE);
116  }
117 
118  /* two possible sets of transition probabilities for OII */
119  if( nelem==ipOXYGEN && ion == 2 && lgAs )
120  {
121  if( nMatch( " NEW" , chCard ) )
122  {
123  dense.lgAsChoose[nelem][ion-1] = true;
124  }
125  else if( nMatch( " OLD" , chCard ) )
126  {
127  /*>>chng 06 mar 30, this is the default */
128  dense.lgAsChoose[nelem][ion-1] = false;
129  }
130  else
131  {
132  fprintf( ioQQQ, " The keyword old or new must appear\n Sorry.\n" );
133  cdEXIT(EXIT_FAILURE);
134  }
135  }
136  else if( nelem==ipSULPHUR && lgDR )
137  {
138  /* change DR data set for S
139  * three cases for S DR - ionbal.nDR_S_guess
140  * 0, default larger of guess and Badnell
141  * 1, pure Badnell
142  * 3, scaled oxygen */
143  if( nMatch( " MIX" , chCard ) )
144  /* this is the default, larger of Badnell or guess */
145  ionbal.nDR_S_guess = 0;
146  else if( nMatch( "PURE" , chCard ) )
147  /* use only Badnell 1991 */
148  ionbal.nDR_S_guess = 1;
149  else if( nMatch( "SCAL" , chCard ) )
150  {
151  /* scaled oxygen */
152  ionbal.nDR_S_guess = 2;
153  i = 5;
154  ionbal.DR_S_scale[0] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
155  if( lgEOL )
156  NoNumb(chCard);
157  for( ion=1; ion<5; ++ion )
158  {
159  /* optionally get rest - if too few specified then use
160  * previously set value */
161  ionbal.DR_S_scale[ion] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
162  if( lgEOL )
163  ionbal.DR_S_scale[ion] = ionbal.DR_S_scale[ion-1];
164  }
165  }
166  else
167  {
168  /* disaster - no match */
169  fprintf(
170  ioQQQ, " One of the keywords MIX, PURE, or SCALe must"
171  "appear on the set atomic physics sulphur dr command.\n"
172  "Sorry.\n" );
173  cdEXIT(EXIT_FAILURE);
174  }
175  }
176  else
177  {
178  fprintf( ioQQQ, " None of the valid set atomic data atoms/ions were found\n Sorry.\n" );
179  cdEXIT(EXIT_FAILURE);
180  }
181  }
182  else
183  {
184  /* a molecule */
185  if( nMatch(" H2 ",chCard) )
186  {
187  /* H2 */
188  chMole = "H2";
189  }
190  else
191  {
192  /* nothing recognized */
193  fprintf( ioQQQ, " No molecule was on this SET ATOMIC DATA command.\n Sorry.\n" );
194  fprintf( ioQQQ, " Use SET ATOMIC DATA ION to change an ion.\n Sorry.\n" );
195  cdEXIT(EXIT_FAILURE);
196  }
197 
198  if( strcmp( chMole , "H2" )==0 )
199  {
200  if( nMatch(" H " , chCard ) && lgCS )
201  {
202  /* change the H - H2 collision data set */
203  i = 5;
204  ion = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
205  if( ion!=2 )
206  TotalInsanity();
207  long int nYear = (long)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
208  if( lgEOL )
209  NoNumb( chCard );
210  if( nYear == 1999 )
211  /* use the coefficients from
212  *>>refer H2 H collision Le Bourlot, J., Pineau des Forets,
213  *>>refercon G., & Flower, D.R. 1999, MNRAS, 305, 802 */
214  h2.lgH2_H_coll_07 = false;
215  else if( nYear == 2007 )
216  /* use the coefficients from
217  *>>refer H2 H collision Wrathmall, S. A., Gusdorf A.,
218  *>>refercon & Flower, D.R. 2007, MNRAS, 382, 133 */
219  h2.lgH2_H_coll_07 = true;
220  else
221  {
222  /* not an option */
223  fprintf(ioQQQ," the SET ATOMIC DATA MOLECULE H2"
224  " H CS command must have year 1999 or 2007.\n" );
225  cdEXIT(EXIT_FAILURE);
226  }
227  }
228  else if( nMatch(" He" , chCard ) && lgCS )
229  {
230  /* change the He - H2 collision data set */
231  if( nMatch("ORNL" , chCard ) )
232  {
233  /* use the coefficients from
234  *>>refer H2 He collision Lee et al in preparation */
235  mole.lgH2_He_ORNL = true;
236  }
237  else if( nMatch( "BOUR",chCard ) )
238  {
239  /* use the coefficients from
240  *>>refer H2 He collision Le Bourlot, J., Pineau des Forets,
241  *>>refercon G., & Flower, D.R. 1999, MNRAS, 305, 802*/
242  mole.lgH2_He_ORNL = false;
243  }
244  else
245  {
246  /* not an option */
247  fprintf(ioQQQ," the SET ATOMIC DATA MOLECULE H2"
248  "He CS command must have ORNL or Le BOURlot.\n" );
249  cdEXIT(EXIT_FAILURE);
250  }
251  }
252  }
253  else
254  TotalInsanity();
255  }
256  }
257 
258  else if( nMatch(" CHA",chCard) && !nMatch( "HO ", chCard ) )
259  {
260  /* set log of minimum charge transfer rate for high ions and H
261  * default of 1.92e-10 set in zero */
262  i = 5;
263  atmdat.HCTAlex = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
264  if( lgEOL )
265  {
266  NoNumb(chCard);
267  }
268  if( atmdat.HCTAlex < 0. )
269  {
270  atmdat.HCTAlex = pow(10.,atmdat.HCTAlex);
271  }
272  }
273 
274  else if( nMatch("CHEM",chCard ) && !nMatch( "HO ", chCard ) )
275  {
276  /* turn on Steve Federman's chemistry */
277  if( nMatch("FEDE",chCard ) )
278  {
279  if( nMatch( " ON " , chCard ) )
280  {
281  /* This turns on the diffuse cloud chemical rates of
282  * >>refer CO chemistry Zsargo, J. & Federman, S. R. 2003, ApJ, 589, 319*/
283  co.lgFederman = true;
284  }
285  else if( nMatch( " OFF" , chCard ) )
286  {
287  co.lgFederman = false;
288  }
289  else
290  {
291  /* this is the default when command used - true */
292  co.lgFederman = true;
293  }
294  }
295  /* >>chng 06 may 30 --NPA. Turn on non-equilibrium chemistry */
296  else if( nMatch(" NON",chCard ) && nMatch( "EQUI", chCard ))
297  {
298 
299  /* option to use effective temperature as defined in
300  * >>refer CO chemistry Zsargo, J. & Federman, S. R. 2003, ApJ, 589, 319
301  * By default, this is false - changed with set chemistry command */
302 
303  co.lgNonEquilChem = true;
304 
305  /* >>chng 06 jul 21 -- NPA. Option to include non-equilibrium
306  * effects for neutral reactions with a temperature dependent rate.
307  * Reasoning is that non-equilibrium chemistry is caused by MHD, and
308  * if magnetic field is only coupled to ions, then neutrals may not be
309  * affected.
310  * >>refer CO chemistry Zsargo, J. & Federman, S. R. 2003, ApJ, 589, 319
311  * By default, this is true - changed with set chemistry command */
312 
313  if( nMatch("NEUT",chCard ) )
314  {
315  if( nMatch( " ON " , chCard ) )
316  {
317  /* This turns on the diffuse cloud chemical rates of
318  * >>refer CO chemistry Zsargo, J. & Federman, S. R. 2003, ApJ, 589, 319*/
319  co.lgNeutrals = true;
320  }
321  else if( nMatch( " OFF" , chCard ) )
322  {
323  co.lgNeutrals = false;
324  }
325  else
326  {
327  /* this is the default when command used - true */
328  co.lgNeutrals = true;
329  }
330  }
331  }
332 
333  /* turn off proton elimination rates, which are of the form
334  *
335  *
336  * A + BH+ --> AB + H+ or
337  * AH + B+ --> AB + H+
338  *
339  * the following paper:
340  *
341  * >>refer CO chemistry Huntress, W. T., 1977, ApJS, 33, 495
342  * says reactions of these types are much less likely than
343  * identical reactions which leave the product AB ionized (AB+),
344  * leaving an H instead of H+ (this is called H atom elimination
345  * currently we only have one reaction of this type, it is
346  * C+ + OH -> CO + H+ */
347  else if( nMatch("PROT",chCard ) && nMatch( "ELIM", chCard ) )
348  {
349  if( nMatch( " ON " , chCard ) )
350  {
351  /* This turns on the diffuse cloud chemical rates of
352  * >>refer CO chemistry Zsargo, J. & Federman, S. R. 2003, ApJ, 589, 319*/
353  co.lgProtElim = true;
354  }
355  else if( nMatch( " OFF" , chCard ) )
356  {
357  co.lgProtElim = false;
358  }
359  else
360  {
361  /* this is the default when command used - true */
362  co.lgProtElim = true;
363  }
364  }
365 
366  else
367  {
368  /* should not have happened ... */
369  fprintf( ioQQQ, " There should have been an option on this SET CHEMISTRY command.\n" );
370  fprintf( ioQQQ, " consult Hazy to find valid options.\n Sorry.\n" );
371  cdEXIT(EXIT_FAILURE);
372  }
373  }
374 
375  /* set collision strength averaging ON / OFF */
376  else if( nMatch("COLL",chCard) && nMatch("STRE",chCard) && nMatch("AVER",chCard) )
377  {
378  if( nMatch(" OFF",chCard) )
379  {
380  iso.lgCollStrenThermAver = false;
381  }
382  else
383  {
384  /* this is default behavior of this command */
385  iso.lgCollStrenThermAver = true;
386  }
387  }
388 
389  else if( nMatch("COVE",chCard) )
390  {
391  iterations.lgConverge_set = true;
392  /* set coverage - limit number of iterations and zones */
393  if( nMatch("FAST",chCard) )
394  {
395  iterations.lim_zone = 1;
396  iterations.lim_iter = 0;
397  }
398  else
399  {
400  iterations.lim_zone = 10;
401  iterations.lim_iter = 1;
402  }
403  }
404 
405  else if( nMatch("CSUP",chCard) )
406  {
407  /* force secondary ionization rate, log entered */
408  i = 5;
410  secondaries.lgCSetOn = true;
411  if( lgEOL )
412  {
413  NoNumb(chCard);
414  }
416  }
417 
418  else if( nMatch(" D/H",chCard) )
419  {
420  /* set deuterium abundance, D to H ratio */
421  i = 5;
422  hydro.D2H_ratio = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
423  if( hydro.D2H_ratio <= 0. || nMatch( " LOG", chCard ) )
424  {
425  hydro.D2H_ratio = pow(10.,hydro.D2H_ratio);
426  }
427  if( lgEOL )
428  {
429  NoNumb(chCard);
430  }
431  }
432 
433  else if( nMatch( " HO " , chCard) && nMatch( "CHAR" , chCard) )
434  {
435  /* which solver will include H + O charge transfer? done in
436  * chemistry by default */
437  if( nMatch( "CHEM" , chCard ) )
438  {
439  /* this is the default */
440  ionbal.lgHO_ct_chem = true;
441  }
442  else if( nMatch( "IONI" , chCard ) )
443  {
444  /* do it in the ionization solver */
445  ionbal.lgHO_ct_chem = false;
446  }
447  else
448  {
449  fprintf( ioQQQ, " I did not recognize a subkey on this SET OH CHARge transfer line.\n" );
450  fprintf( ioQQQ, " Valid keys are CHEMistry and IONIzation.\n" );
451  cdEXIT(EXIT_FAILURE);
452  }
453  }
454 
455  else if( nMatch("12C1",chCard) )
456  {
457  /* set the 12C to 13C abundance ratio - default is 30 */
458  i = 5;
459  /* first two numbers on the line are 12 and 13 - we don't want them */
462 
463  /* now we can get the ratio */
465  if( lgEOL )
466  NoNumb(chCard);
467 
468  if( co.C12_C13_isotope_ratio <= 0. || nMatch( " LOG", chCard ) )
470  }
471 
472  /* set dynamics ... */
473  else if( nMatch("DYNA",chCard) )
474  {
475  /* set dynamics advection length */
476  if( nMatch("ADVE",chCard) && nMatch("LENG",chCard) )
477  {
478  /* <0 => relative fraction of length, + val in cm */
479  i = 5;
481  if( lgEOL )
482  NoNumb( chCard );
483  /* if fraction is present, then number was linear fraction, if not present
484  * then physical length in cm, log10 */
485  if( nMatch("FRAC",chCard) )
486  {
487  /* we scanned in the number, if it is a negative then it is the log of the fraction */
488  if( dynamics.AdvecLengthInit <= 0. )
490 
491  /* make neg sign as flag for this in dynamics.c */
492  dynamics.AdvecLengthInit *= -1.;
493  }
494  else
495  {
496  /* fraction did not occur, the number is the log of the length in cm -
497  * convert to linear cm */
499  }
500  }
501  else if( nMatch("PRES",chCard) && nMatch("MODE",chCard) )
502  {
503  dynamics.lgSetPresMode = true;
504  if( nMatch("SUBS",chCard) )
505  {
506  /* subsonic */
507  strcpy( dynamics.chPresMode , "subsonic" );
508  }
509  else if( nMatch("SUPE",chCard) )
510  {
511  /* supersonic */
512  strcpy( dynamics.chPresMode , "supersonic" );
513  }
514  else if( nMatch("STRO",chCard) )
515  {
516  /* strong d */
517  strcpy( dynamics.chPresMode , "strongd" );
518  }
519  else if( nMatch("ORIG",chCard) )
520  {
521  /* original */
522  strcpy( dynamics.chPresMode , "original" );
523  }
524  }
525  else if( nMatch("ANTI",chCard) && nMatch("DEPT",chCard) )
526  {
527  dynamics.lgSetPresMode = true;
528  strcpy( dynamics.chPresMode , "antishock" );
529  /* shock depth */
530  i = 5;
531  /* get log of shock depth in cm */
532  dynamics.ShockDepth = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
533  if( lgEOL )
534  NoNumb( chCard );
535  dynamics.ShockDepth = pow( 10., dynamics.ShockDepth );
536  }
537  else if( nMatch("ANTI",chCard) && nMatch("MACH",chCard) )
538  {
539  dynamics.lgSetPresMode = true;
540  strcpy( dynamics.chPresMode , "antishock-by-mach" );
541  /* Mach number */
542  i = 5;
543  /* get (isothermal) Mach number where we want antishock to occur */
544  dynamics.ShockMach = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
545  if( lgEOL )
546  NoNumb( chCard );
547  }
548  else if( nMatch("RELA",chCard) )
549  {
550  /* set how many iterations we will start with, before allowing
551  * changes. This allows the solution to relax to an equilibrium */
552  i = 5;
553  dynamics.n_initial_relax = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
554  if( lgEOL )
555  NoNumb( chCard );
556  else if( dynamics.n_initial_relax < 2 )
557  {
558  fprintf(ioQQQ," First iteration to relax dynamics must be > 1."
559  "It was %li. Sorry.\n",
561  cdEXIT(EXIT_FAILURE);
562  }
563  }
564  else if( nMatch("SHOC",chCard) && nMatch("DEPT",chCard) )
565  {
566  dynamics.lgSetPresMode = true;
567  strcpy( dynamics.chPresMode , "shock" );
568  /* shock depth */
569  i = 5;
570  /* get log of shock depth in cm */
571  dynamics.ShockDepth = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
572  if( lgEOL )
573  NoNumb( chCard );
574  dynamics.ShockDepth = pow( 10., dynamics.ShockDepth );
575  }
576  else
577  {
578  /* should not have happened ... */
579  fprintf( ioQQQ, " There should have been an option on this SET DYNAMICS command.\n" );
580  fprintf( ioQQQ, " consult Hazy to find valid options.\n Sorry.\n" );
581  cdEXIT(EXIT_FAILURE);
582  }
583  }
584 
585  else if( nMatch("DIDZ",chCard) )
586  {
587  /* set parameter to do with choice of dr;
588  * par is the largest optical depth to allow in the zone.
589  * >>chng 96 jan 08 had been two numbers - dtau1 removed */
590  i = 5;
591  radius.drChange = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
592  if( radius.drChange <= 0. )
593  {
595  }
596  if( lgEOL )
597  {
598  NoNumb(chCard);
599  }
600  }
601 
602  /* something to do with electron density */
603  else if( nMatch("EDEN",chCard) )
604  {
605  /* >>chng 02 apr 20, from ERROR to TOLERANCE to be parallel to set temp equivalent,
606  * >>chng 04 jun 03, but also accept error as keyword */
607  if( nMatch("CONV",chCard) || nMatch("ERRO",chCard))
608  {
609  /* keyword is eden convergence
610  * set tolerance in eden match */
611  i = 5;
612  conv.EdenErrorAllowed = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
613  if( lgEOL )
614  {
615  NoNumb(chCard);
616  }
617 
618  if( conv.EdenErrorAllowed < 0. )
619  {
621  }
622  }
623 
624  else if( nMatch("SOLV",chCard) )
625  {
626  /* the electron density solver */
627  if( nMatch("NEW",chCard) )
628  {
629  /* new method */
630  strcpy( conv.chSolverEden , "new" );
631  }
632  else
633  {
634  /* default method */
635  strcpy( conv.chSolverEden , "simple" );
636  }
637  }
638  else
639  {
640  /* no keyword, assume log of electron density */
641  i = 5;
642  /* set the electron density */
643  dense.EdenSet = (realnum)pow(10.,FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL));
644  /* warn that this model is meaningless */
645  phycon.lgPhysOK = false;
646  }
647  }
648 
649  else if( nMatch("FINE",chCard ) && nMatch("CONT",chCard ) )
650  {
651  /* set fine continuum resolution - an element name, used to get
652  * thermal width, and how many resolution elements to use to resolve
653  * a line of this element at 1e4 K
654  * first get an element name, nelem is atomic number on C scale
655  * default is iron */
656  if( (rfield.fine_opac_nelem = GetElem(chCard))<0 )
657  {
658  fprintf( ioQQQ, " An element name must appear on this line\n Sorry.\n" );
659  cdEXIT(EXIT_FAILURE);
660  }
661  i = 5;
662  /* set the number of resolution elements in HWHM at 1e4 K for turbulent
663  * velocity field, default is one element */
664  rfield.fine_opac_nresolv = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
665  if( rfield.fine_opac_nresolv< 1 )
666  {
667  fprintf( ioQQQ, " The number of resolution elements within FWHM of line must appear\n Sorry.\n" );
668  cdEXIT(EXIT_FAILURE);
669  }
670  }
671 
672  /* set grain command - but not set H2 grain command */
673  else if( nMatch("GRAI",chCard ) && !nMatch(" H2 ",chCard) )
674  {
675  if( nMatch("HEAT",chCard ) )
676  {
677  /* scale factor to change grain heating as per Allers et al. */
678  i = 5;
680  /* warn that this model is not the best we can do */
681  phycon.lgPhysOK = false;
682  if( lgEOL )
683  {
684  NoNumb(chCard);
685  }
686  }
687  else
688  {
689  fprintf( ioQQQ, " A keyword must appear on the SET GRAIN line - options are HEAT \n Sorry.\n" );
690  cdEXIT(EXIT_FAILURE);
691  }
692  }
693 
694  /* these are the "set Leiden hack" commands, used to turn off physics and
695  * sometimes replace with simple approximation */
696  else if( nMatch("LEID",chCard ) && nMatch("HACK",chCard ) )
697  {
698  if( nMatch( "H2* " , chCard ) && nMatch( " OFF" , chCard ) )
699  {
700  /* turn off reactions with H2* in the chemistry network */
701  hmi.lgLeiden_Keep_ipMH2s = false;
702  /* warn that this model is not the best we can do */
703  phycon.lgPhysOK = false;
704  }
705  else if( nMatch( "CR " , chCard ) && nMatch( " OFF" , chCard ) )
706  {
707  /* the CR Leiden hack option - turn off CR excitations of H2 */
708  hmi.lgLeidenCRHack = false;
709 
710  }
711  else if( nMatch("RATE",chCard ) && nMatch("UMIS",chCard ))
712  {
713  /* This command will use the rates given in the UMIST database, It
714  * will set to zero many reactions in hmole_step.c that are not
715  * in UMIST */
716  co.lgUMISTrates = false;
717  }
718  }
719 
720  /* set H2 ... */
721  else if( nMatch(" H2 ",chCard) )
722  {
723  i = 5;
724  ip = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
725  if( ip != 2 )
726  {
727  fprintf( ioQQQ, " The first number on this line must be the 2 in H2\n Sorry.\n" );
728  cdEXIT(EXIT_FAILURE);
729  }
730 
731  /* SET H2 SMALL MODEL or SOLOMON - which approximation to Solomon process */
732  if( nMatch("SOLO",chCard) || ( nMatch("SMAL",chCard)&&nMatch("MODE",chCard) ) )
733  {
734  if( nMatch("SOLO",chCard) )
735  {
736  /* warn that Solomon will not be used forever */
737  fprintf(ioQQQ,"PROBLEM - *set H2 Solomon* has been changed to *set H2 small model*."\
738  " This is OK for now but it may not work in a future version.\n");
739  }
740  if( nMatch("TH85", chCard ) )
741  {
742  /* rate from is eqn A8 of Tielens and Hollenbach 85a, */
744  }
745  else if( nMatch( " BHT" , chCard ) )
746  {
747  /* the improved H2 formalism given by
748  *>>refer H2 dissoc Burton, M.G., Hollenbach, D.J. & Tielens, A.G.G.M
749  >>refcon 1990, ApJ, 365, 620 */
751  }
752  else if( nMatch( "BD96" , chCard ) )
753  {
754  /* this rate is equation 23 of
755  *>>refer H2 dissoc Bertoldi, F., & Draine, B.T., 1996, 458, 222 */
756  /* this is the default */
758  }
759  else if( nMatch( "ELWE" , chCard ) )
760  {
761  /* this rate is equation 23 of
762  *>>refer H2 dissoc Elwert et al., in preparation */
763  /* this is the default */
765  }
766  else
767  {
768  fprintf( ioQQQ, " One of the keywords TH85, _BHT, BD96 or ELWErt must appear.\n Sorry.\n" );
769  cdEXIT(EXIT_FAILURE);
770  }
771  }
772 
773  /* series of commands that deal with grains */
774  /* which approximation to grain formation pumping */
775  if( nMatch("GRAI",chCard ) && nMatch("FORM",chCard ) && nMatch("PUMP",chCard ) )
776  {
777  if( nMatch( "DB96" , chCard ) )
778  {
779  /* Draine & Bertoldi 1996 */
780  hmi.chGrainFormPump = 'D';
781  }
782  else if( nMatch( "TAKA" , chCard ) )
783  {
784  /* Takahashi 2001 */
785  hmi.chGrainFormPump = 'T';
786  }
787  else if( nMatch( "THER" , chCard ) )
788  {
789  /* thermal distribution, upper right column of page 239 of
790  *>>refer H2 formation Le Bourlot, J, 1991, A&A, 242, 235 */
791  hmi.chGrainFormPump = 't';
792  }
793  else
794  {
795  fprintf( ioQQQ, " The grain form pump option is wrong.\n Sorry.\n" );
796  cdEXIT(EXIT_FAILURE);
797  }
798  }
799 
800  /* which approximation to Jura rate */
801  else if( nMatch("JURA",chCard) )
802  {
803  if( nMatch("TH85", chCard ) )
804  {
805  /* rate from is eqn A8 of Tielens and Hollenbach 85a*/
806  hmi.chJura = 'T';
807  }
808  else if( nMatch( "CT02" , chCard ) )
809  {
810  /* this rate is equation Cazeux & Tielens */
811  hmi.chJura = 'C';
812  }
813  else if( nMatch( "SN99" , chCard ) )
814  {
815  /* this rate is from Sternberg & Neufeld 99 */
816  hmi.chJura = 'S';
817  }
818  else if( nMatch( "RATE" , chCard ) )
819  {
820  /* set H2 rate - enters log of Jura rate - F for fixed
821  * no dependence on grain properties
822  * had been C, a bug since triggered Cazeux & Tielens
823  * >>chng 07 jan 10, bug caught by Robin Williams & Fixed by Nick Abel */
824  hmi.chJura = 'F';
825  hmi.rate_h2_form_grains_set = pow(10.,FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL) );
826  if( lgEOL )
827  {
828  /* no number on the line so use Jura's value, 3e-17
829  * >>refer H2 Jura Jura, M. 1975, ApJ, 197, 575 */
831  }
832  }
833  else if( nMatch( "SCAL" , chCard ) )
834  {
835  /* this is a scale factor to multiply the Jura rate */
836  hmi.ScaleJura = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
837  /* log or negative number means log was entered */
838  if( nMatch( " LOG" , chCard ) || hmi.ScaleJura<0. )
839  {
840  hmi.ScaleJura = (realnum)pow( 10., (double)hmi.ScaleJura );
841  }
842  if( lgEOL )
843  {
844  NoNumb( chCard );
845  }
846 
847  /* option to vary scale factor */
848  if( optimize.lgVarOn )
849  {
851  strcpy( optimize.chVarFmt[optimize.nparm], "SET H2 JURA SCALE %f" );
852 
853  /* pointer to where to write */
855 
856  /* log of Jura rate scale factor will be parameter */
858  optimize.vincr[optimize.nparm] = 0.3f;
859 
860  ++optimize.nparm;
861  }
862  }
863  else
864  {
865  fprintf( ioQQQ, " The Jura rate option is wrong.\n Sorry.\n" );
866  cdEXIT(EXIT_FAILURE);
867  }
868  }
869 
870  /* what temperature to use for binding energy, Tad in Le Bourlot, J., 2000, A&A, 360, 656-662 */
871  else if( nMatch(" TAD",chCard) )
872  {
873  hmi.Tad = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
874  if( lgEOL )
875  NoNumb(chCard);
876  /* log if <= 10. unless linear key appears too */
877  if( hmi.Tad <=10. && !nMatch("LINE",chCard) )
878  hmi.Tad = (realnum)pow((realnum)10.f,hmi.Tad);
879  }
880 
881  else if( nMatch("FRAC",chCard ) )
882  {
883  /* this is special option to force H2 abundance to value for testing
884  * this factor will multiply the hydrogen density to become the H2 density
885  * no attempt to conserve particles, or do the rest of the molecular equilibrium
886  * set consistently is made */
887  hmi.H2_frac_abund_set = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
888  if( lgEOL )
889  NoNumb( chCard );
890 
891  /* a number <= 0 is the log of the ratio */
892  if( hmi.H2_frac_abund_set <= 0. )
894  /* don't let it exceed 0.5 */
895  /* >>chng 03 jul 19, from 0.5 to 0.4999, do not want atomic density exactly zero */
897  }
898  }
899 
900  /* this is a scale factor that changes the n(H0)*1.7e-4 that is added to the
901  * electron density to account for collisions with atomic H. it is an order of
902  * magnitude guess, so this command provides ability to see whether it affects results */
903  else if( nMatch("HCOR",chCard) )
904  {
905  i = 5;
906  dense.HCorrFac = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
907  if( lgEOL )
908  NoNumb( chCard );
909  }
910 
911  else if( nMatch(" PAH",chCard) )
912  {
913  /* set one of several possible PAH abundance distribution functions */
914  if( nMatch(" H0 " , chCard ) )
915  {
916  /* the default, abundance depends on H0 / Htot */
917  strcpy( gv.chPAH_abundance_fcn , "H0" );
918  }
919  else if( nMatch("CONS" , chCard ) )
920  {
921  /* constant PAH abundance */
922  strcpy( gv.chPAH_abundance_fcn , "CON" );
923  }
924 
925  else if(nMatch("BAKE",chCard ) )
926  {
927  /* turn on simple PAH heating from Bakes & Tielens - this is very approximate */
928  /*>>>refer PAH heating Bakes, E.L.O., & Tielens, A.G.G.M. 1994, ApJ, 427, 822 */
929  gv.lgBakesPAH_heat = true;
930  /* warn that this model is not the best we can do */
931  phycon.lgPhysOK = false;
932  }
933  else
934  {
935  fprintf( ioQQQ, " one of the keywords H0, BAKES, or CONStant must appear.\n Sorry." );
936  cdEXIT(EXIT_FAILURE);
937  }
938  }
939 
940  else if( nMatch("PRES",chCard) && nMatch("IONI",chCard) )
941  {
942  /* set limit to number of calls from pressure to ionize solver,
943  * this set limit to conv.nPres2Ioniz */
944  i = 5;
945  conv.limPres2Ioniz = (long)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
946  if( lgEOL )
947  {
948  NoNumb(chCard);
949  }
950  else if( conv.limPres2Ioniz <= 0 )
951  {
952  fprintf( ioQQQ, " The limit must be greater than zero.\n Sorry." );
953  cdEXIT(EXIT_FAILURE);
954  }
955  }
956  /* something to do with pressure */
957  else if( nMatch("PRES",chCard) )
958  {
959  /* tolerance on pressure convergence */
960  if( nMatch("CONV",chCard) )
961  {
962  /* keyword is tolerance
963  * set tolerance or relative error in pressure match */
964  i = 5;
966  if( lgEOL )
967  NoNumb(chCard);
968 
969  if( conv.PressureErrorAllowed < 0. )
971  }
972 
973  else
974  {
975  /* >>chng 04 mar 02, printout had been wrong, saying TOLErange
976  * rather than CONVergence. Nick Abel */
977  fprintf( ioQQQ, " I didn\'t recognize a key on this SET PRESSURE line.\n" );
978  fprintf( ioQQQ, " The ones I know about are: CONVergence.\n" );
979  cdEXIT(EXIT_FAILURE);
980  }
981  }
982  else if( nMatch("RECO",chCard) && nMatch("MBIN",chCard) )
983  {
984  /* dielectronic recombination */
985  if( nMatch("DIEL",chCard) && nMatch("ECTR",chCard) )
986  {
987  /* change various factors for dielectronic recombination */
988  if( nMatch("KLUD",chCard) )
989  {
990  if( nMatch("BADN",chCard) )
991  {
992  /* set true to use Badnell mean in place of existing hacks */
994  }
995  else
996  {
997  int j;
998  /* option to turn off guess of dielectronic rec coefficient for 3rd row elements */
999  i = 3;
1000  /* this is first call, no number, lgEOL true but zero returned, the intended effect*/
1001  ionbal.GuessDiel[0] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
1002  /* >>chng 06 jan 30 also turn off guess if GuessDiel set zero */
1003  if( ionbal.GuessDiel[0]<=0. || nMatch(" OFF",chCard) )
1004  ionbal.lg_guess_coef = false;
1005 
1006  for( j=1; j<4; ++j )
1007  {
1008  ionbal.GuessDiel[j] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
1009  /* here lgEOL means to use previous values */
1010  if( lgEOL )
1011  ionbal.GuessDiel[j] = ionbal.GuessDiel[j-1];
1012  }
1013  }
1014  /* option to add log normal noise */
1015  if( nMatch("NOISE" , chCard ) )
1016  {
1017  i = 3;
1018  ionbal.guess_noise = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
1019  if( lgEOL )
1020  ionbal.guess_noise = 2.;
1021  /* Seed the random-number generator with current time so that
1022  * the numbers will be different every time we run. */
1023  init_genrand( (unsigned)time( NULL ) );
1024  }
1025  else
1026  {
1027  ionbal.guess_noise = 0.;
1028  }
1029  }
1030 
1031  else if( nMatch("BURG" , chCard) )
1032  {
1033  /* modify suppression of burgess dielectronic recombinations */
1034  if( nMatch(" ON ",chCard) )
1035  {
1036  /* leave suppression on - this is default state */
1037  ionbal.lgSupDie[0] = true;
1038  }
1039 
1040  else if( nMatch(" OFF",chCard) )
1041  {
1042  /* turn suppression off */
1043  ionbal.lgSupDie[0] = false;
1044  }
1045 
1046  else
1047  {
1048  fprintf( ioQQQ, " flag ON or OFF must appear.\n" );
1049  cdEXIT(EXIT_FAILURE);
1050  }
1051  }
1052 
1053  else if( nMatch("NUSS" , chCard) )
1054  {
1055  /* modify suppression of Nussbaumer and Storey dielectronic recombination */
1056  if( nMatch(" ON ",chCard) )
1057  {
1058  /* turn suppression on */
1059  ionbal.lgSupDie[1] = true;
1060  }
1061  else if( nMatch(" OFF",chCard) )
1062  {
1063  /* leave suppression off - this is default state */
1064  ionbal.lgSupDie[1] = false;
1065  }
1066  else
1067  {
1068  fprintf( ioQQQ, " flag ON or OFF must appear.\n" );
1069  cdEXIT(EXIT_FAILURE);
1070  }
1071  }
1072 
1073  else if( nMatch("BADN" , chCard ) )
1074  {
1075  /* use the new (in early 2006) Badnell DR rates */
1076  if( nMatch( " OFF" , chCard ) )
1077  {
1079  }
1080  else
1081  {
1083  }
1084 
1085  /* option to print rates then exit */
1086  if( nMatch("PRIN" , chCard ) )
1088  }
1089  else
1090  {
1091  fprintf( ioQQQ, " key KLUDge, BURGess, NUSSbaumer, or BADNell must appear.\n" );
1092  cdEXIT(EXIT_FAILURE);
1093  }
1094  }
1095  /* radiative recombination */
1096  else if( nMatch("RADI",chCard) && nMatch("ATIV",chCard) )
1097  {
1098  if( nMatch("BADN" , chCard ) )
1099  {
1100  /* use the new (in early 2006) Badnell DR rates */
1101  if( nMatch( " OFF" , chCard ) )
1102  {
1104  }
1105  else
1106  {
1108  }
1109 
1110  /* option to print rates then exit */
1111  if( nMatch("PRIN" , chCard ) )
1113  }
1114  else
1115  {
1116  fprintf( ioQQQ, " key BADNell must appear.\n" );
1117  cdEXIT(EXIT_FAILURE);
1118  }
1119  }
1120  else
1121  {
1122  fprintf( ioQQQ, " key RADIATIve or DIELECTRonic must appear on set recombination command.\n" );
1123  cdEXIT(EXIT_FAILURE);
1124  }
1125  }
1126 
1127  else if( nMatch("SPEC",chCard) )
1128  {
1129  /* "set spectrum" for optional parameters for punch spectrum command */
1130  if( nMatch("RANG" , chCard ) )
1131  {
1132  /* energy range option - argument is energy in Rydbergs */
1133  i = 5;
1134  /* >>chng 05 aug 25, rm pow(10., - was wrong */
1137  if( lgEOL )
1138  {
1139  NoNumb(chCard);
1140  }
1141  /* confirm that wavelengths are in correct order */
1143  {
1144  fprintf( ioQQQ, " The limits must be in increasing order.\n" );
1145  cdEXIT(EXIT_FAILURE);
1146  }
1147  }
1148 
1149  else if( nMatch("RESO" , chCard ) )
1150  {
1151  /* resolving power, default is zero, which means leave at native resolution */
1152  /* >>chng 05 aug 25, following had pow(number), so crashed for high resolving power
1153  * bug caught by Yan Changshou */
1154  i = 5;
1156  if( lgEOL )
1157  {
1158  NoNumb(chCard);
1159  }
1160  }
1161  }
1162 
1163  else if( nMatch(" DR ",chCard) )
1164  {
1165  /* set zone thickness by forcing drmax and drmin */
1166  i = 5;
1167  /* at this stage linear, but default is log */
1168  radius.sdrmax = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
1169  if( !nMatch("LINE",chCard) )
1170  {
1171  /* linear was not on command, so default is log */
1172  radius.sdrmax = pow( 10. , radius.sdrmax );
1173  }
1174  if( lgEOL )
1175  {
1176  NoNumb(chCard);
1177  }
1178 
1179  /* NB these being equal are tested in convinittemp to set dr */
1181  if( radius.sdrmax < DEPTH_OFFSET*1e4 )
1182  {
1183  fprintf( ioQQQ, "\n Thicknesses less than about %.0e will NOT give accurate results. If tricking the code\n",
1184  DEPTH_OFFSET*1e4 );
1185  fprintf( ioQQQ, " into computing emissivities instead of intensities, try to instead use a thickness of unity,\n" );
1186  fprintf( ioQQQ, " and then multiply (divide) the results by the necessary thickness (product of densities).\n\n" );
1187  cdEXIT(EXIT_FAILURE);
1188  }
1189  }
1190 
1191  else if( nMatch("DRMA",chCard) )
1192  {
1193  /* set maximum zone thickness */
1194  i = 5;
1195  radius.sdrmax = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
1196  if( lgEOL )
1197  NoNumb(chCard);
1198 
1199  /* log if less than 38 or if log keyword occurs */
1200  if( radius.sdrmax < log10(BIGFLOAT) || nMatch(" LOG",chCard) )
1201  radius.sdrmax = pow(10.,radius.sdrmax);
1202  }
1203 
1204  else if( nMatch("DRMI",chCard) )
1205  {
1206  /* option to set minimum zone thickness */
1207  i = 5;
1208  radius.sdrmin = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
1209  if( lgEOL )
1210  NoNumb(chCard);
1211 
1212  /* log if less than 38 or if log keyword occurs */
1213  if( radius.sdrmin < log10(BIGFLOAT ) || nMatch(" LOG",chCard) )
1214  radius.sdrmin = pow(10.,radius.sdrmin);
1215 
1216  radius.lgSMinON = true;
1217  }
1218 
1219  else if( nMatch("FLXF",chCard) )
1220  {
1221  /* faintest continuum flux to consider */
1222  i = 5;
1223  rfield.FluxFaint = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
1224  if( lgEOL )
1225  {
1226  NoNumb(chCard);
1227  }
1228  if( rfield.FluxFaint < 0. )
1229  {
1231  }
1232  }
1233 
1234  else if( nMatch("LINE",chCard) && nMatch("PREC",chCard) )
1235  {
1236  /* set line precision (number of significant figures)
1237  * this affects all aspects of reading and writing lines */
1238  i = 5;
1239  LineSave.sig_figs = (int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
1240  if( lgEOL )
1241  {
1242  NoNumb(chCard);
1243  }
1244  else if( LineSave.sig_figs > 6 )
1245  {
1246  fprintf( ioQQQ, " set line precision currently only works for up to 6 significant figures.\n" );
1247  cdEXIT(EXIT_FAILURE);
1248  }
1249  /* option to print main line array as a single column */
1250  /* prt.lgPrtLineArray = false; */
1251  }
1252 
1253  /* >>chng 00 dec 08, command added by Peter van Hoof */
1254  else if( nMatch("NFNU",chCard) )
1255  {
1256  /* set nFnu [incident_reflected] [incident_transmitted] [diffuse_inward] [diffuse_outward]
1257  * command for specifying what to include in the nFnu entries in LineSv */
1258  /* >>chng 01 nov 12, also accept form with space */
1259  /* "incident reflected" keyword */
1260  prt.lgSourceReflected = nMatch("NT R",chCard) || nMatch("NT_R",chCard);
1261  /* "incident transmitted" keyword */
1262  prt.lgSourceTransmitted = nMatch("NT_T",chCard) || nMatch("NT T",chCard);
1263  /* "diffuse inward" keyword */
1264  prt.lgDiffuseInward = nMatch("SE_I",chCard) || nMatch("SE I",chCard);
1265  /* "diffuse outward" keyword */
1266  prt.lgDiffuseOutward = nMatch("SE_O",chCard) || nMatch("SE O",chCard);
1267 
1268  /* at least one of these needs to be set ! */
1271  {
1272  fprintf( ioQQQ, " set nFnu expects one or more of the following keywords:\n" );
1273  fprintf( ioQQQ, " INCIDENT_REFLECTED, INCIDENT_TRANSMITTED, DIFFUSE_INWARD, DIFFUSE_OUTWARD\n" );
1274  cdEXIT(EXIT_FAILURE);
1275  }
1276  /* automatically print the nFnu items in the output - it is not necessary to also include
1277  * a print continua command if this is entered */
1278  prt.lgPrnDiff = true;
1279  }
1280 
1281  else if( nMatch("IND2",chCard) )
1282  {
1283  if( nMatch(" ON ",chCard) )
1284  {
1285  /* set flag saying to off or on induced two photon processes */
1286  iso.lgInd2nu_On = true;
1287  }
1288  else if( nMatch(" OFF",chCard) )
1289  {
1290  iso.lgInd2nu_On = false;
1291  }
1292  else
1293  {
1294  fprintf( ioQQQ, " set ind2 needs either ON or OFF.\n" );
1295  cdEXIT(EXIT_FAILURE);
1296  }
1297  }
1298 
1299  else if( nMatch("TEMP",chCard) )
1300  {
1301  /* set something to do with temperature, currently solver, tolerance, floor */
1302  if( nMatch("SOLV",chCard) )
1303  {
1304  /* solver */
1305  /* the electron density solver */
1306  if( nMatch("NEW",chCard) )
1307  {
1308  /* new method */
1309  strcpy( conv.chSolverTemp , "new" );
1310  }
1311  else
1312  {
1313  /* default method */
1314  strcpy( conv.chSolverTemp , "simple" );
1315  }
1316  }
1317 
1318  else if( nMatch("FLOO",chCard) )
1319  {
1320  i = 5;
1321  StopCalc.TeFloor = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
1322 
1323  /* linear option */
1324  if( StopCalc.TeFloor <= 10. && !nMatch("LINE",chCard) )
1325  {
1326  StopCalc.TeFloor = (realnum)pow(10.,StopCalc.TeFloor);
1327  }
1328 
1329  if( lgEOL )
1330  {
1331  NoNumb(chCard);
1332  }
1333 
1334  if( StopCalc.TeFloor < 2.8 )
1335  {
1336  fprintf( ioQQQ, " TE < 3K, reset to 2.8K.\n" );
1337  StopCalc.TeFloor = 2.8f;
1338  }
1339  }
1340 
1341  else if( nMatch("CONV",chCard) || nMatch("TOLE",chCard) )
1342  {
1343  /* error tolerance in heating cooling match, number is error/total */
1344  i = 5;
1346  if( lgEOL )
1347  {
1348  NoNumb(chCard);
1349  }
1350  if( conv.HeatCoolRelErrorAllowed <= 0. )
1351  {
1353  }
1354  }
1355 
1356  else
1357  {
1358  fprintf( ioQQQ, "\nI did not recognize a keyword on this SET TEPERATURE command.\n%s\n" , chCard);
1359  fprintf( ioQQQ, "The keywords are SOLVer, FLOOr, and CONVergence.\n" );
1360  cdEXIT(EXIT_FAILURE);
1361  }
1362  }
1363 
1364  else if( nMatch("TEST",chCard) )
1365  {
1366  /* set flag saying to turn on some test - this is in cddefines.h in the global namespace */
1367  lgTestCodeEnabled = true;
1368  }
1369 
1370  else if( nMatch("TRIM",chCard) )
1371  {
1372  /* set trim upper or lower, for ionization stage trimming
1373  * ion trimming ionization trimming
1374  * in routine TrimStage */
1375  i = 5;
1376  if( nMatch("UPPE",chCard) )
1377  {
1378  /* set trim upper - either set value or turn off */
1379  double save = ionbal.trimhi;
1380  ionbal.trimhi = pow(10.,FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL));
1381  if( lgEOL && nMatch(" OFF",chCard) )
1382  {
1383  /* turn off upward trimming */
1384  lgEOL = false;
1385  ionbal.lgTrimhiOn = false;
1386  /* reset high limit to proper value */
1387  ionbal.trimhi = save;
1388  }
1389  }
1390 
1391  else if( nMatch("LOWE",chCard) )
1392  {
1393  /* set trim lower */
1394  ionbal.trimlo = pow(10.,FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL));
1395  }
1396 
1397  /* turn off ionization stage trimming */
1398  else if( nMatch("SMAL",chCard) || nMatch(" OFF",chCard) )
1399  {
1400  /* set small limits to both upper and lower limit*/
1403  lgEOL = false;
1404  }
1405 
1406  else
1407  {
1408  /* set trim upper */
1409  ionbal.trimhi = pow(10.,FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL));
1410 
1411  /* set trim lower to same number */
1413  }
1414 
1415  if( lgEOL )
1416  {
1417  NoNumb(chCard);
1418  }
1419 
1420  if( ionbal.trimlo >= 1. || ionbal.trimhi >= 1. )
1421  {
1422  fprintf( ioQQQ, " number must be negative since log\n" );
1423  cdEXIT(EXIT_FAILURE);
1424  }
1425  }
1426 
1427  else if( nMatch("SKIP",chCard) )
1428  {
1429  /* skip punch command, for punching every n't point */
1430  i = 5;
1431  punch.ncPunchSkip = (long)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
1432  if( lgEOL )
1433  {
1434  NoNumb(chCard);
1435  }
1436  }
1437 
1438  else if( nMatch(" UTA",chCard) )
1439  {
1440  /* set UTA command, to determine which UTA data set to use
1441  * default is to use Gu et al. 2006 data set */
1442  if( nMatch("GU06" , chCard ) )
1443  {
1444  if( nMatch(" OFF" , chCard ) )
1445  {
1446  /* use Behar et al. 2001 */
1447  ionbal.lgInnerShell_Gu06 = false;
1448  }
1449  else
1450  {
1451  /* the default, use
1452  *>>refer Fe UTA Gu, M.F., Holczer, T., Behar, E. & Kahn, S.M.
1453  *>>refercon 2006, ApJ, 641, 1227 */
1454  ionbal.lgInnerShell_Gu06 = true;
1455  }
1456  }
1457  else if( nMatch("KISI" , chCard ) )
1458  {
1459  if( nMatch(" OFF" , chCard ) )
1460  {
1461  /* the default, do not use it */
1463  }
1464  else
1465  {
1466  /* use the data from Kisielius et al. 2003, MNRAS, 344, 696 */
1468  }
1469  }
1470  }
1471 
1472  else if( nMatch("EAKH",chCard) )
1473  {
1474  /* set WeakHeatCool, threshold on punch heating and cooling, default 0.05 */
1475  i = 5;
1476  punch.WeakHeatCool = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
1477 
1478  if( lgEOL )
1479  {
1480  NoNumb(chCard);
1481  }
1482 
1483  if( punch.WeakHeatCool < 0. )
1484  {
1486  }
1487  }
1488 
1489  else if( nMatch("KSHE",chCard) )
1490  {
1491  /* upper limit to opacities for k-shell ionization */
1492  i = 5;
1494  if( lgEOL )
1495  {
1496  NoNumb(chCard);
1497  }
1498 
1499  if( continuum.EnergyKshell == 0. )
1500  {
1501  /* arg of 0 causes upper limit to energy array */
1503  }
1504 
1505  else if( continuum.EnergyKshell < 194. )
1506  {
1507  fprintf( ioQQQ, " k-shell energy must be greater than 194 Ryd\n" );
1508  cdEXIT(EXIT_FAILURE);
1509  }
1510  }
1511 
1512  else if( nMatch("NCHR",chCard) )
1513  {
1514  /* option to set the number of charge states for grain model */
1515  long ii = 5;
1516  double val = FFmtRead(chCard,&ii,INPUT_LINE_LENGTH,&lgEOL);
1517 
1518  if( lgEOL )
1519  {
1520  NoNumb(chCard);
1521  }
1522  else
1523  {
1524  long nChrg = nint(val);
1525  if( nChrg < 2 || nChrg > NCHS )
1526  {
1527  fprintf( ioQQQ, " illegal value for number of charge states: %ld\n", nChrg );
1528  fprintf( ioQQQ, " choose a value between 2 and %d\n", NCHS );
1529  fprintf( ioQQQ, " or increase NCHS in grainvar.h and recompile\n" );
1530  cdEXIT(EXIT_FAILURE);
1531  }
1532  else
1533  {
1534  SetNChrgStates(nChrg);
1535  }
1536  }
1537  }
1538 
1539  else if( nMatch("NEGO",chCard) )
1540  {
1541  /* punch negative opacities if they occur, set negopac */
1542  opac.lgNegOpacIO = true;
1543  }
1544 
1545  else if( nMatch("NEND",chCard) )
1546  {
1547  /* default limit to number of zones to be computed
1548  * only do this if nend is NOT currently left at its default
1549  * nend is set to nEndDflt in routine zero
1550  * this command only has effect if stop zone not entered */
1551  if( geometry.lgEndDflt )
1552  {
1553  i = 5;
1554  /* this is default limit to number of zones, change it to this value */
1555  geometry.nEndDflt = (long)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
1556  geometry.lgEndDflt = false;
1557  if( lgEOL )
1558  {
1559  NoNumb(chCard);
1560  }
1561 
1562  /* now change all limits, for all iterations, to this value */
1563  for( i=0; i < iterations.iter_malloc; i++ )
1564  {
1566  }
1567  }
1568  }
1569 
1570  else if( nMatch("TSQD",chCard) )
1571  {
1572  /* upper limit for highest density considered in the
1573  * Peimbert-style t^2 section of the printout */
1574  i = 5;
1575  peimbt.tsqden = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
1576 
1577  if( lgEOL )
1578  {
1579  NoNumb(chCard);
1580  }
1581  peimbt.tsqden = (realnum)pow((realnum)10.f,peimbt.tsqden);
1582  }
1583 
1584  else if( nMatch("NMAP",chCard) )
1585  {
1586  /* how many steps in plot or punch of heating-cooling map */
1587  i = 5;
1588  hcmap.nMapStep = (long)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
1589  if( lgEOL )
1590  {
1591  NoNumb(chCard);
1592  }
1593  }
1594 
1595  else if( nMatch("NUME",chCard) && nMatch("DERI",chCard) )
1596  {
1597  /* this is an option to use numerical derivatives for heating and cooling */
1598  NumDeriv.lgNumDeriv = true;
1599  }
1600 
1601  else if( nMatch("PATH",chCard) )
1602  {
1603  fprintf( ioQQQ, " The SET PATH command is no longer supported.\n" );
1604  fprintf( ioQQQ, " Please set the correct path in the file path.h.\n" );
1605  fprintf( ioQQQ, " Or set the environment variable CLOUDY_DATA_PATH.\n Sorry.\n" );
1606  cdEXIT(EXIT_FAILURE);
1607  }
1608 
1609  else if( nMatch("PHFI",chCard) )
1610  {
1611  /* which version of PHFIT to use, 1995 or 1996 */
1612  i = 5;
1613  ip = (long)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
1614  if( lgEOL )
1615  {
1616  NoNumb(chCard);
1617  }
1618 
1619  if( ip == 1995 )
1620  {
1621  /* option to go back to old results, pre op */
1622  t_ADfA::Inst().set_version( PHFIT95 );
1623  }
1624  else if( ip == 1996 )
1625  {
1626  /* default is to use newer results, including opacity project */
1627  t_ADfA::Inst().set_version( PHFIT96 );
1628  }
1629  else
1630  {
1631  fprintf( ioQQQ, " Two possible values are 1995 and 1996.\n" );
1632  cdEXIT(EXIT_FAILURE);
1633  }
1634  }
1635 
1636  /* set punch xxx command */
1637  else if( nMatch("PUNC",chCard) )
1638  {
1639  if( nMatch("HASH",chCard) )
1640  {
1641  /* to use the hash option there must be double quotes on line - were there? */
1642  if( !lgQuotesFound )
1643  {
1644  fprintf( ioQQQ, " I didn\'t recognize a key on this SET PUNCH HASH line.\n" );
1645  cdEXIT(EXIT_FAILURE);
1646  }
1647  /* specify the terminator between punch output sets - normally a set of hash marks */
1648  /*
1649  * get any string within double quotes, and return it as
1650  * null terminated string
1651  * also sets name in OrgCard and chCard to blanks so that
1652  * do not trigger off it later */
1653  if( strcmp( chString_quotes_lowercase , "return" ) == 0 )
1654  {
1655  /* special case - return becomes new line */
1656  sprintf(punch.chHashString , "\n" );
1657  }
1658  else if( strcmp( chString_quotes_lowercase , "time" ) == 0 )
1659  {
1660  /* special case where output time between iterations */
1661  sprintf(punch.chHashString , "TIME_DEP" );
1662  }
1663  else
1664  {
1665  /* usual case, simply copy what is in quotes */
1666  sprintf(punch.chHashString , chString_quotes_lowercase );
1667  }
1668  }
1669 
1670  else if( nMatch("WIDT",chCard) )
1671  {
1672  /* set punch line width for contrast in continuum plots */
1673  i = 5;
1674  /* units are km/sec */
1675  punch.PunchLWidth = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
1676 
1677  /* lgEOL if no number hit, could have been c */
1678  if( lgEOL )
1679  {
1680  if( nMatch(" C " , chCard ) )
1681  {
1682  /* no number but _c_ entered, so enter speed of light in km/s */
1684  }
1685  else
1686  {
1687  NoNumb(chCard);
1688  }
1689  }
1690 
1691  if( punch.PunchLWidth <= 0. )
1692  {
1693  fprintf( ioQQQ, " line width must be greater than zero.\n" );
1694  cdEXIT(EXIT_FAILURE);
1695  }
1696 
1697  /* >>chng 06 dec 06, value of PunchLWidth is km/s but SPEEDLIGHT
1698  * is cm/s - add 1e5 before compare with SPEEDLIGHT */
1699  else if( punch.PunchLWidth*0.9999e5 > SPEEDLIGHT )
1700  {
1701  fprintf( ioQQQ, " line width must be entered in km/s and be less than c.\n" );
1702  cdEXIT(EXIT_FAILURE);
1703  }
1704  /* this is the factor that is used to add the lines to the continuum,
1705  * 15 converts to cm/s */
1707 
1708  }
1709 
1710  else if( nMatch("PREF",chCard) )
1711  {
1712  /* specify a prefix before all punch filenames */
1713  /*
1714  * get any string within double quotes, and return it as
1715  * null terminated string
1716  * also sets name in OrgCard and chCard to blanks so that
1717  * do not trigger off it later */
1718 
1719  strcpy( punch.chFilenamePrefix , chString_quotes_lowercase );
1720  }
1721 
1722  else if( nMatch("FLUS",chCard) )
1723  {
1724  /* flus the output after every iteration */
1725  punch.lgFLUSH = true;
1726  }
1727 
1728  else
1729  {
1730  fprintf( ioQQQ, " There should have been an option on this command.\n" );
1731  fprintf( ioQQQ, " Valid options are HASH and PREFIX.\n" );
1732  cdEXIT(EXIT_FAILURE);
1733  }
1734  }
1735 
1736  /* set continuum .... options */
1737  else if( nMatch("CONT" , chCard ) )
1738  {
1739  if( nMatch("RESO" , chCard ) )
1740  {
1741  /* set resolution, get factor that will multiply continuum resolution that
1742  * is contained in the file continuum_mesh.ini */
1743  i = 5;
1745  if( lgEOL )
1746  {
1747  NoNumb(chCard);
1748  }
1749  /* negative numbers were logs */
1750  if( continuum.ResolutionScaleFactor <= 0. )
1752  }
1753 
1754  else if( nMatch("SHIE" , chCard ) )
1755  {
1756  /* set continuum shielding function */
1757  /* these are all possible values of rt.nLineContShield,
1758  * first is default, these are set with set continuum shielding */
1759  /*#define LINE_CONT_SHIELD_PESC 1*/
1760  /*#define LINE_CONT_SHIELD_FEDERMAN 2*/
1761  /*#define LINE_CONT_SHIELD_FERLAND 3*/
1762  if( nMatch("PESC" , chCard ) )
1763  {
1764  /* this uses an inward looking escape probability */
1766  }
1767  else if( nMatch("FEDE" , chCard ) )
1768  {
1769  /* set continuum shielding Federman,
1770  * this is the default, and uses the appendix of
1771  * >>refer RT continuum shielding Federman, S.R., Glassgold, A.E., &
1772  * >>refercon Kwan, J. 1979, ApJ, 227, 466*/
1774  }
1775  else if( nMatch("FERL" , chCard ) )
1776  {
1778  }
1779  else if( nMatch("NONE" , chCard ) )
1780  {
1781  /* turn off continuum shielding */
1782  rt.nLineContShield = 0;
1783  }
1784  else
1785  {
1786  fprintf( ioQQQ, " I didn\'t recognize a key on this SET CONTINUUM SHIELDing line.\n" );
1787  fprintf( ioQQQ, " The ones I know about are: PESC, FEDErman, & FERLand.\n" );
1788  cdEXIT(EXIT_FAILURE);
1789  }
1790  }
1791 
1792  else
1793  {
1794  fprintf( ioQQQ, " I didn\'t recognize a key on this SET CONTINUUM line.\n" );
1795  fprintf( ioQQQ, " The ones I know about are: RESOlution and SHIEld.\n" );
1796  cdEXIT(EXIT_FAILURE);
1797  }
1798  }
1799 
1800  else
1801  {
1802  fprintf( ioQQQ, " I didn\'t recognize a key on this SET command.\n" );
1803  cdEXIT(EXIT_FAILURE);
1804  }
1805 
1806  return;
1807 }

Generated for cloudy by doxygen 1.8.1.1