cloudy  trunk
 All Data Structures Namespaces 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-2017 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 "input.h"
6 #include "prt.h"
7 #include "rt.h"
8 #include "phycon.h"
9 #include "optimize.h"
10 #include "hcmap.h"
11 #include "hmi.h"
12 #include "iterations.h"
13 #include "conv.h"
14 #include "secondaries.h"
15 #include "rfield.h"
16 #include "ionbal.h"
17 #include "numderiv.h"
18 #include "dynamics.h"
19 #include "iso.h"
20 #include "predcont.h"
21 #include "save.h"
22 #include "stopcalc.h"
23 #include "opacity.h"
24 #include "peimbt.h"
25 #include "radius.h"
26 #include "atmdat.h"
27 #include "continuum.h"
28 #include "grains.h"
29 #include "grainvar.h"
30 #include "parser.h"
31 #include "lines.h"
32 #include "monitor_results.h"
33 #include "mole.h"
34 #include "dense.h"
35 #include "hyperfine.h"
36 #include "service.h"
37 
38 void ParseSet(Parser &p)
39 {
40  long int ip;
41  string chString_quotes_lowercase;
42  DEBUG_ENTRY( "ParseSet()" );
43 
44  /* first check for any strings within quotes - this will get the string
45  * and blank it out, so that these are not confused with keywords. if
46  * double quotes not present routine returns unity, zero if found*/
47  bool lgQuotesFound = true;
48  if (p.GetQuote(chString_quotes_lowercase))
49  lgQuotesFound = false;
50 
51  /* commands to set certain variables, "SET XXX=" */
52 
53  if( p.nMatch("LYA") && p.nMatch("21CM") )
54  {
55  /* set the shape of the Lya source function at line center. Important for
56  * 21 cm pumping and the resulting spin temperature
57  */
58  if( p.nMatch("EXCI") )
59  {
60  /* the default, the 2p / 1s excitation temperature */
62  }
63  else if( p.nMatch("KINE") )
64  {
65  /* the gas kinetic temperature */
67  }
68  else if( p.nMatch("CONS") )
69  {
70  /* S_nu = constant, as assumed in stellar atmosphere calculations */
72  }
73  }
74 
75  else if (p.nMatch("ASSE") && p.nMatch("ABOR"))
76  {
77  /* set assert abort command, to tell the code to raise SIGABRT when an
78  * assert is thrown - this can be caught by the debugger. */
79  cpu.i().setAssertAbort(true);
80  }
81 
82  else if (p.nMatch("MONI") && p.nMatch("SCIE"))
83  {
84  /* print monitored values using scientific notation. Useful when an asserted value is
85  * less than 10^-4 of the normalization line. */
86  lgPrtSciNot = true;
87  }
88 
89  /* check energy conservation on every zone - relatively slow */
90  else if (p.nMatch("CHECK") && p.nMatch("ENERGY") && p.nMatch("EVERY")
91  && p.nMatch("ZONE"))
92  {
94  }
95 
96  else if (p.nMatch("UPDA") && p.nMatch("COUP") &&
97  p.nMatch("EVER") && p.nMatch("ION") )
98  {
99  conv.lgUpdateCouplings = true;
100  }
101 
102  /* enable Dere07 collisional ionization rate coeffs */
103  else if (p.nMatch(" COLL") && p.nMatch(" IONIZ"))
104  {
105  if (p.nMatch(" HYBRID"))
106  {
108  }
109  else if (p.nMatch(" DIMA"))
110  {
112  }
113  else
114  {
115  fprintf(ioQQQ,
116  " \nPROBLEM Unrecognized set collisional ionization data option.\n");
117  fprintf(ioQQQ, " Valid options are Dima or Hybrid.\n");
118  fprintf(ioQQQ, " See Hazy 1 for details.\n");
119  cdEXIT( EXIT_FAILURE );
120  }
121  }
122 
123  else if( p.nMatch(" H2 " ) && p.nMatch( "CONT" ) && p.nMatch( "DISS" ))
124  {
125  if( p.nMatch( "STAN" ) )
126  {
127  /* This uses on the H2 direct photodissociation cross sections
128  * calculated by Philip Stancil */
129  mole_global.lgStancil = true;
130  }
131 
132  else if( p.nMatch( "AD69" ) )
133  {
134  /* Use constant cross section from
135  >>refer H2 dissoc Allison, A.C. & Dalgarno, A. 1969,
136  Atomic Data, 1, 91 */
137 
138  mole_global.lgStancil = false;
139  }
140  else
141  {
142  /* should not have happened ... */
143  fprintf( ioQQQ, " There should have been an option on this SET H2 CONTinuum DISSociation command.\n" );
144  fprintf( ioQQQ, " consult Hazy to find valid options.\n Sorry.\n" );
146  }
147 
148  }
149 
150  else if( p.nMatch(" CHA") && !p.nMatch( "HO ") && !p.nMatch( "HHE") )
151  {
152  /* set log of minimum charge transfer rate for high ions and H
153  * default of 1.92e-10 set in zero */
154  atmdat.HCTAlex = p.FFmtRead();
155  if (p.lgEOL())
156  {
157  p.NoNumb("minimum charge transfer rate");
158  }
159  if (atmdat.HCTAlex < 0.)
160  {
162  }
163  }
164  else if( p.nMatch("CHEM") && !p.nMatch( "HO ") && !p.nMatch(" HHE") )
165  {
166  /* turn on Steve Federman's chemistry */
167  if (p.nMatch("FEDE"))
168  {
169  if (p.nMatch(" ON "))
170  {
171  /* This turns on the diffuse cloud chemical rates of
172  * >>refer CO chemistry Zsargo, J. & Federman, S. R. 2003, ApJ, 589, 319*/
173  mole_global.lgFederman = true;
174  }
175  else if( p.nMatch( " OFF" ) )
176  {
177  mole_global.lgFederman = false;
178  }
179  else
180  {
181  /* this is the default when command used - true */
182  mole_global.lgFederman = true;
183  }
184  }
185  /* >>chng 06 may 30 --NPA. Turn on non-equilibrium chemistry */
186  else if (p.nMatch(" NON") && p.nMatch("EQUI"))
187  {
188 
189  /* option to use effective temperature as defined in
190  * >>refer CO chemistry Zsargo, J. & Federman, S. R. 2003, ApJ, 589, 319
191  * By default, this is false - changed with set chemistry command */
192 
194 
195  /* >>chng 06 jul 21 -- NPA. Option to include non-equilibrium
196  * effects for neutral reactions with a temperature dependent rate.
197  * Reasoning is that non-equilibrium chemistry is caused by MHD, and
198  * if magnetic field is only coupled to ions, then neutrals may not be
199  * affected.
200  * >>refer CO chemistry Zsargo, J. & Federman, S. R. 2003, ApJ, 589, 319
201  * By default, this is true - changed with set chemistry command */
202 
203  if (p.nMatch("NEUT"))
204  {
205  if (p.nMatch(" ON "))
206  {
207  /* This turns on the diffuse cloud chemical rates of
208  * >>refer CO chemistry Zsargo, J. & Federman, S. R. 2003, ApJ, 589, 319*/
209  mole_global.lgNeutrals = true;
210  }
211  else if( p.nMatch( " OFF" ) )
212  {
213  mole_global.lgNeutrals = false;
214  }
215  else
216  {
217  /* this is the default when command used - true */
218  mole_global.lgNeutrals = true;
219  }
220  }
221  }
222 
223  /* turn off proton elimination rates, which are of the form
224  *
225  *
226  * A + BH+ --> AB + H+ or
227  * AH + B+ --> AB + H+
228  *
229  * the following paper:
230  *
231  * >>refer CO chemistry Huntress, W. T., 1977, ApJS, 33, 495
232  * says reactions of these types are much less likely than
233  * identical reactions which leave the product AB ionized (AB+),
234  * leaving an H instead of H+ (this is called H atom elimination
235  * currently we only have one reaction of this type, it is
236  * C+ + OH -> CO + H+ */
237  else if (p.nMatch("PROT") && p.nMatch("ELIM"))
238  {
239  if (p.nMatch(" ON "))
240  {
241  /* This turns on the diffuse cloud chemical rates of
242  * >>refer CO chemistry Zsargo, J. & Federman, S. R. 2003, ApJ, 589, 319*/
243  mole_global.lgProtElim = true;
244  }
245  else if( p.nMatch( " OFF" ) )
246  {
247  mole_global.lgProtElim = false;
248  }
249  else
250  {
251  /* this is the default when command used - true */
252  mole_global.lgProtElim = true;
253  }
254  }
255 
256  else
257  {
258  /* should not have happened ... */
259  fprintf(ioQQQ,
260  " There should have been an option on this SET CHEMISTRY command.\n");
261  fprintf(ioQQQ, " consult Hazy to find valid options.\n Sorry.\n");
263  }
264  }
265 
266  /* set collision strength averaging ON / OFF */
267  else if( p.nMatch("COLL") && p.nMatch("STRE") && p.nMatch("AVER") )
268  {
269  if( p.nMatch(" OFF") )
270  {
272  }
273  else
274  {
275  /* this is default behavior of this command */
277  }
278  }
279 
280  else if (p.nMatch("COVE"))
281  {
282  iterations.lgConverge_set = true;
283  /* set coverage - limit number of iterations and zones */
284  if (p.nMatch("FAST"))
285  {
286  iterations.lim_zone = 1;
287  iterations.lim_iter = 0;
288  }
289  else
290  {
291  iterations.lim_zone = 10;
292  iterations.lim_iter = 1;
293  }
294  }
295 
296  else if (p.nMatch("CSUP"))
297  {
298  /* force H^0 secondary ionization rate, log entered */
300  secondaries.lgCSetOn = true;
301  if (p.lgEOL())
302  {
303  p.NoNumb("secondary ionization rate");
304  }
306  }
307 
308  else if (p.nMatch("CUMU"))
309  {
310  if (p.nMatch("OFF"))
311  {
312  strcpy(rfield.chCumuType,"NONE");
313  }
314  else if (p.nMatch("MASS"))
315  {
316  strcpy(rfield.chCumuType,"MASS");
317  }
318  else if (p.nMatch("FLUX"))
319  {
320  strcpy(rfield.chCumuType,"FLUX");
321  }
322  else
323  {
324  fprintf(ioQQQ,
325  " Did not recognize option on 'set cumulative' command."
326  " It should be FLUX, MASS or OFF. Sorry.\n");
328  }
329  }
330 
331  else if (p.nMatch(" D/H"))
332  {
333  /* set deuterium abundance, D to H ratio */
334  double tmp = p.FFmtRead();
335  if (p.lgEOL())
336  {
337  fprintf(ioQQQ,
338  "The command 'set D/H' has been deprecated.\n"
339  "Please use 'element hydrogen isotopes' instead.\n");
340  }
341  else
342  {
343  if (tmp <= 0. || p.nMatch(" LOG"))
344  {
345  tmp = exp10(tmp);
346  }
347  fprintf(ioQQQ,
348  "The command 'set D/H' has been deprecated.\n"
349  "Please use 'element hydrogen isotopes (1, 1) (2, %g)' instead.\n",
350  tmp);
351  }
352  cdEXIT( EXIT_FAILURE );
353  }
354 
355  else if( p.nMatch("DENS") && p.nMatch("TOLE") )
356  {
357  /* set error in total gas-phase density of each element, including molecules */
359  if( p.lgEOL() )
360  {
361  p.NoNumb("density tolerance");
362  }
363  if( conv.GasPhaseAbundErrorAllowed <= 0. )
364  {
366  }
367  }
368 
369  else if( p.nMatch( " HO " ) && p.nMatch( "CHAR" ) )
370  {
371  fprintf(ioQQQ, " The SET HO CHAR command is no longer supported.\n" );
373  }
374 
375  else if( p.nMatch( " HHE" ) && p.nMatch( "CHAR" ) )
376  {
377  fprintf(ioQQQ, " The SET HHE CHAR command is no longer supported.\n" );
379  }
380 
381  else if( p.nMatch("12C1") )
382  {
383  /* set the 12C to 13C abundance ratio - default is 30 */
384  /* first two numbers on the line are 12 and 13 - we don't want them */
385  (void) p.FFmtRead();
386  (void) p.FFmtRead();
387 
388  /* now we can get the ratio */
389  realnum c12c13 = (realnum) p.FFmtRead();
390  if (p.lgEOL())
391  {
392  fprintf(ioQQQ, "Deprecated option. Please use command: \"element carbon isotope\"\n");
393  }
394  else
395  {
396  if (c12c13 <= 0. || p.nMatch(" LOG"))
397  c12c13 = exp10(c12c13);
398  fprintf(ioQQQ, "Deprecated option. Please use command: \"element carbon isotope (12, %.2f), (13, 1)\"\n",
399  c12c13);
400  }
402  }
403 
404  /* set dynamics ... */
405  else if (p.nMatch("DYNA"))
406  {
407  /* set dynamics advection length */
408  if (p.nMatch("ADVE") && p.nMatch("LENG"))
409  {
410  /* <0 => relative fraction of length, + val in cm */
412  if (p.lgEOL())
413  p.NoNumb("advection length");
414  /* if fraction is present, then number was linear fraction, if not present
415  * then physical length in cm, log10 */
416  if (p.nMatch("FRAC"))
417  {
418  /* we scanned in the number, if it is a negative then it is the log of the fraction */
419  if (dynamics.AdvecLengthInit <= 0.)
421 
422  /* make neg sign as flag for this in dynamics.c */
423  dynamics.AdvecLengthInit *= -1.;
424  }
425  else
426  {
427  /* fraction did not occur, the number is the log of the length in cm -
428  * convert to linear cm */
430  }
431  }
432  else if (p.nMatch("PRES") && p.nMatch("MODE"))
433  {
434  dynamics.lgSetPresMode = true;
435  if (p.nMatch("SUBS"))
436  {
437  /* subsonic */
438  strcpy(dynamics.chPresMode, "subsonic");
439  }
440  else if (p.nMatch("SUPE"))
441  {
442  /* supersonic */
443  strcpy(dynamics.chPresMode, "supersonic");
444  }
445  else if (p.nMatch("STRO"))
446  {
447  /* strong d */
448  strcpy(dynamics.chPresMode, "strongd");
449  }
450  else if (p.nMatch("ORIG"))
451  {
452  /* original */
453  strcpy(dynamics.chPresMode, "original");
454  }
455  }
456  else if (p.nMatch("ANTI") && p.nMatch("DEPT"))
457  {
458  dynamics.lgSetPresMode = true;
459  strcpy(dynamics.chPresMode, "antishock");
460  /* shock depth */
461  /* get log of shock depth in cm */
463  if (p.lgEOL())
464  p.NoNumb("antishock depth");
466  }
467  else if (p.nMatch("ANTI") && p.nMatch("MACH"))
468  {
469  dynamics.lgSetPresMode = true;
470  strcpy(dynamics.chPresMode, "antishock-by-mach");
471  /* Mach number */
472  /* get (isothermal) Mach number where we want antishock to occur */
474  if (p.lgEOL())
475  p.NoNumb("antishock-by-mach");
476  }
477  else if (p.nMatch("RELA"))
478  {
479  /* set how many iterations we will start with, before allowing
480  * changes. This allows the solution to relax to an equilibrium */
481  dynamics.n_initial_relax = (long int) p.FFmtRead();
482  if (p.lgEOL())
483  p.NoNumb("relaxation cycles before start of dynamics");
484  else if (dynamics.n_initial_relax < 2)
485  {
486  fprintf(ioQQQ,
487  " First iteration to relax dynamics must be > 1."
488  "It was %li. Sorry.\n", dynamics.n_initial_relax);
490  }
491  }
492  else if (p.nMatch("SHOC") && p.nMatch("DEPT"))
493  {
494  dynamics.lgSetPresMode = true;
495  strcpy(dynamics.chPresMode, "shock");
496  /* shock depth */
497  /* get log of shock depth in cm */
499  if (p.lgEOL())
500  p.NoNumb("shock depth");
502  }
503  else if (p.nMatch("POPU") && p.nMatch("EQUI"))
504  {
505  // set dynamics populations eqauilibrium
506  // force equilibrium populations
507  dynamics.lgEquilibrium = true;
508  }
509  else
510  {
511  /* should not have happened ... */
512  fprintf(ioQQQ,
513  " There should have been an option on this SET DYNAMICS command.\n");
514  fprintf(ioQQQ, " consult Hazy to find valid options.\n Sorry.\n");
516  }
517  }
518 
519  else if (p.nMatch("DIDZ"))
520  {
521  /* set parameter to do with choice of dr;
522  * par is the largest optical depth to allow in the zone.
523  * >>chng 96 jan 08 had been two numbers - dtau1 removed */
525  if (radius.drChange <= 0.)
526  {
528  }
529  if (p.lgEOL())
530  {
531  p.NoNumb("largest optical depth allowed in zone");
532  }
533  }
534 
535  /* something to do with electron density */
536  else if (p.nMatch("EDEN"))
537  {
538  /* set eden convergence sets convergence criterion, keyword parallel with
539  * set temperature convergence, but also accept (older) error as keyword */
540  if (p.nMatch("CONV") || p.nMatch("ERRO"))
541  {
542  /* keyword is eden convergence
543  * set tolerance in eden match */
545  if (p.lgEOL())
546  {
547  p.NoNumb("electron density error allowed");
548  }
549 
550  if (conv.EdenErrorAllowed < 0.)
551  {
553  }
554  }
555  else if(p.nMatch("FRACTION"))
556  {
557  /* set eden fraction sets the ratio eden/hden */
559  if (p.lgEOL())
560  {
561  p.NoNumb("electron density fraction");
562  }
563 
564  if (dense.EdenFraction <= 0. )
565  {
567  }
568  /* warn that this model is meaningless */
569  phycon.lgPhysOK = false;
570  }
571  else if (p.nMatch("SOLV"))
572  {
573  /* options are vWDB (default) and SECAnt */
574  if (p.nMatch("VWDB"))
575  {
576  strcpy( conv.chSolverEden, "vWDB" );
577  }
578  else if (p.nMatch("SECA"))
579  {
580  strcpy( conv.chSolverEden, "SECA" );
581  }
582  else
583  {
584  fprintf(ioQQQ,
585  "'set eden solver' options are 'vWDB' and 'SECA'.\n");
586  cdEXIT( EXIT_FAILURE );
587  }
588  }
589  else
590  {
591  /* no keyword, assume log of electron density */
592  /* set the electron density */
594  if (p.lgEOL())
595  {
596  p.NoNumb("electron density");
597  }
598 
599  /* warn that this model is meaningless */
600  phycon.lgPhysOK = false;
601  }
602  }
603 
604  /* Stop using gbar to fill in dBase transitions */
605  else if( p.nMatch("GBAR"))
606  {
607  if( p.nMatch(" OFF ") )
608  {
609  atmdat.lgGbarOn = false;
610  }
611  }
612 
613  /* something to do with electron density */
614  else if (p.nMatch("IONI") && p.nMatch("TOLE"))
615  {
616  /* keyword is ionization tolerance */
618  if (p.lgEOL())
619  {
620  p.NoNumb("ionization error allowed");
621  }
622 
623  if (conv.IonizErrorAllowed < 0.)
624  {
626  }
627  }
628 
629  else if (p.nMatch("ISOTOPE") || p.nMatch("ISOTOPO") )
630  {
631  if( p.nMatch(" ALL ") )
632  {
633  fprintf(ioQQQ,
634  "The command 'set isotopes all' has been deprecated.\n"
635  "Please use 'element isotopes all' instead.\n");
636  cdEXIT( EXIT_FAILURE );
637  }
638  }
639 
640 
641  else if (p.nMatch("FINE") && p.nMatch("CONT"))
642  {
643  /* set fine continuum resolution - an element name, used to get
644  * thermal width, and how many resolution elements to use to resolve
645  * a line of this element at 1e4 K
646  * first get an element name, nelem is atomic number on C scale
647  * default is iron */
648  if ((rfield.fine_opac_nelem = p.GetElem()) < 0)
649  {
650  fprintf(ioQQQ,
651  " An element name must appear on this line\n Sorry.\n");
653  }
654  /* set the number of resolution elements in HWHM at 1e4 K for turbulent
655  * velocity field, default is one element */
656  rfield.fine_opac_nresolv = (long int) p.FFmtRead();
657  if (rfield.fine_opac_nresolv < 1)
658  {
659  fprintf(ioQQQ,
660  " The number of resolution elements within FWHM of line must appear\n Sorry.\n");
662  }
663 
664  /* option to specify the lower and upper limits of the fine continuum
665  * the upper limit is always optional. */
666  if (!p.lgEOL())
667  {
668  realnum lower_limit = (realnum) p.FFmtRead();
669  if (lower_limit < 0)
670  lower_limit = exp10(lower_limit);
671 
672  if (lower_limit > 0.2f)
673  fprintf(
674  ioQQQ,
675  " The fine continuum lower limit is quite high (%f Ryd). Please check.\n",
676  lower_limit);
677 
678  /* reset to coarse limits if arguments are beyond them */
679  rfield.fine_ener_lo = MAX2(rfield.emm(), lower_limit);
680 
681  if (!p.lgEOL())
682  {
683  realnum upper_limit = (realnum) p.FFmtRead();
684  if (upper_limit < 0)
685  upper_limit = exp10(upper_limit);
686 
687  if (upper_limit < 10.f)
688  fprintf(
689  ioQQQ,
690  " The fine continuum upper limit is quite low (%f Ryd). Please check.\n",
691  upper_limit);
692 
693  /* reset to coarse limits if arguments are beyond them */
694  rfield.fine_ener_hi = MIN2(rfield.egamry(), upper_limit);
695  }
696  }
697  }
698 
699  /* set grain command - but not set H2 grain command */
700  else if (p.nMatch("GRAI") && !p.nMatch(" H2 "))
701  {
702  if (p.nMatch("HEAT"))
703  {
704  /* scale factor to change grain heating as per Allers et al. */
706  /* warn that this model is not the best we can do */
707  phycon.lgPhysOK = false;
708  if (p.lgEOL())
709  {
710  p.NoNumb("grain heating");
711  }
712  }
713  else
714  {
715  fprintf(ioQQQ,
716  " A keyword must appear on the SET GRAIN line - options are HEAT \n Sorry.\n");
718  }
719  }
720 
721  /* these are the "set Leiden hack" commands, used to turn off physics and
722  * sometimes replace with simple approximation */
723  else if (p.nMatch("LEID") && p.nMatch("HACK"))
724  {
725  if (p.nMatch("H2* ") && p.nMatch(" OFF"))
726  {
727  /* turn off reactions with H2* in the chemistry network */
728  hmi.lgLeiden_Keep_ipMH2s = false;
729  /* warn that this model is not the best we can do */
730  phycon.lgPhysOK = false;
731  }
732  else if (p.nMatch("CR ") && p.nMatch(" OFF"))
733  {
734  /* the CR Leiden hack option - turn off CR excitations of H2 */
735  hmi.lgLeidenCRHack = false;
736  }
737  else if( p.nMatch("RATE") && p.nMatch("UMIS"))
738  {
739  /* This command will use the rates given in the UMIST database,
740  * It will set to zero many reactions that are not in UMIST */
741  mole_global.lgLeidenHack = true;
742  }
743  }
744 
745  /* set H2 ... */
746  else if (p.nMatch(" H2 "))
747  {
748  ip = (long int) p.FFmtRead();
749  if (ip != 2)
750  {
751  fprintf(ioQQQ,
752  " The first number on this line must be the 2 in H2\n Sorry.\n");
754  }
755 
756  /* SET H2 SMALL MODEL or SOLOMON - which approximation to Solomon process */
757  if (p.nMatch("SOLO") || (p.nMatch("SMAL") && p.nMatch("MODE")))
758  {
759  if (p.nMatch("SOLO"))
760  {
761  /* warn that Solomon will not be used forever */
762  fprintf(ioQQQ,
763  "PROBLEM - *set H2 Solomon* has been changed to *set H2 small model*."
764  " This is OK for now but it may not work in a future version.\n");
765  }
766  if (p.nMatch("TH85"))
767  {
768  /* rate from is eqn A8 of Tielens and Hollenbach 85a, */
770  }
771  else if (p.nMatch(" BHT"))
772  {
773  /* the improved H2 formalism given by
774  *>>refer H2 dissoc Burton, M.G., Hollenbach, D.J. & Tielens, A.G.G.M
775  >>refcon 1990, ApJ, 365, 620 */
777  }
778  else if (p.nMatch("BD96"))
779  {
780  /* this rate is equation 23 of
781  *>>refer H2 dissoc Bertoldi, F., & Draine, B.T., 1996, 458, 222 */
782  /* this is the default */
784  }
785  else if (p.nMatch("ELWE"))
786  {
787  /* this rate is equation 23 of
788  *>>refer H2 dissoc Elwert et al., in preparation */
789  /* this is the default */
791  }
792  else
793  {
794  fprintf(ioQQQ,
795  " One of the keywords TH85, _BHT, BD96 or ELWErt must appear.\n Sorry.\n");
797  }
798  }
799 
800  /* series of commands that deal with grains */
801  /* which approximation to grain formation pumping */
802  if (p.nMatch("GRAI") && p.nMatch("FORM") && p.nMatch("PUMP"))
803  {
804  if (p.nMatch("DB96"))
805  {
806  /* Draine & Bertoldi 1996 */
807  hmi.chGrainFormPump = 'D';
808  }
809  else if (p.nMatch("TAKA"))
810  {
811  /* the default from
812  * >>refer H2 pump Takahashi, Junko, 2001, ApJ, 561, 254-263 */
813  hmi.chGrainFormPump = 'T';
814  }
815  else if (p.nMatch("THER"))
816  {
817  /* thermal distribution, upper right column of page 239 of
818  *>>refer H2 formation Le Bourlot, J, 1991, A&A, 242, 235 */
819  hmi.chGrainFormPump = 't';
820  }
821  else if (p.nMatch(" OFF"))
822  {
823  /* disable grain formation pumping */
824  hmi.chGrainFormPump = ' ';
825  }
826  else
827  {
828  fprintf(ioQQQ,
829  " The grain form pump option is wrong.\n Sorry.\n");
831  }
832  }
833 
834  /* which approximation to Jura rate */
835  else if (p.nMatch("JURA"))
836  {
837  if (p.nMatch("TH85"))
838  {
839  /* rate from is eqn A8 of Tielens and Hollenbach 85a*/
840  hmi.chJura = 'T';
841  }
842  else if (p.nMatch("CT02"))
843  {
844  /* this rate is equation Cazeux & Tielens */
845  hmi.chJura = 'C';
846  }
847  else if (p.nMatch("ELRD"))
848  {
849  /* this rate is equation C.1 from Rollig et al., and includes Eley-Rideal Effect*/
850  hmi.chJura = 'E';
851  }
852  else if (p.nMatch("SN99"))
853  {
854  /* this rate is from Sternberg & Neufeld 99 */
855  hmi.chJura = 'S';
856  }
857  else if (p.nMatch("RATE"))
858  {
859  /* set H2 rate - enters log of Jura rate - F for fixed
860  * no dependence on grain properties
861  * had been C, a bug since triggered Cazeux & Tielens
862  * >>chng 07 jan 10, bug caught by Robin Williams & Fixed by Nick Abel */
863  hmi.chJura = 'F';
865  if (p.lgEOL())
866  {
867  /* no number on the line so use Jura's value, 3e-17
868  * >>refer H2 Jura Jura, M. 1975, ApJ, 197, 575 */
870  }
871  }
872  else if (p.nMatch("SCAL"))
873  {
874  /* this is a scale factor to multiply the Jura rate */
875  hmi.ScaleJura = (realnum) p.FFmtRead();
876  /* log or negative number means log was entered */
877  if (p.nMatch(" LOG") || hmi.ScaleJura < 0.)
878  {
880  }
881  if (p.lgEOL())
882  p.NoNumb("scale for Jura rate");
883 
884  /* option to vary scale factor */
885  if (optimize.lgVarOn)
886  {
889  "SET H2 JURA SCALE %f LOG");
890 
891  /* pointer to where to write */
893 
894  /* log of Jura rate scale factor will be parameter */
895  optimize.vparm[0][optimize.nparm] = (realnum) log10(
896  hmi.ScaleJura);
897  optimize.vincr[optimize.nparm] = 0.3f;
898 
899  ++optimize.nparm;
900  }
901  }
902  else
903  {
904  fprintf(ioQQQ, " The Jura rate option is wrong.\n Sorry.\n");
906  }
907  }
908 
909  /* what temperature to use for binding energy, Tad in Le Bourlot, J., 2000, A&A, 360, 656-662 */
910  else if (p.nMatch(" TAD"))
911  {
912  hmi.Tad = (realnum) p.FFmtRead();
913  if (p.lgEOL())
914  p.NoNumb("temperature for binding energy");
915  /* log if <= 10. unless linear key appears too */
916  if (hmi.Tad <= 10. && !p.nMatch("LINE"))
917  hmi.Tad = exp10(hmi.Tad);
918  }
919 
920  else if (p.nMatch("FRAC"))
921  {
922  /* this is special option to force H2 abundance to value for testing
923  * this factor will multiply the hydrogen density to become the H2 density
924  * no attempt to conserve particles, or do the rest of the molecular equilibrium
925  * set consistently is made */
927  if (p.lgEOL())
928  p.NoNumb("H2 fractional abundance");
929 
930  /* a number <= 0 is the log of the ratio */
931  if (hmi.H2_frac_abund_set <= 0.)
933  /* don't let it exceed 0.5 */
934  /* >>chng 03 jul 19, from 0.5 to 0.4999, do not want atomic density exactly zero */
936  }
937 #if 0
938  else if( p.nMatch("FORM") && p.nMatch("SCAL") )
939  {
940  /* this is special option to scale H2 formation rate.
941  * In the fully molecular or fully ionized limits,
942  * this should supercede the above "FRAC" option because
943  * it allows the same thing without breaking the chemistry
944  * or any conservation checks */
945  hmi.H2_formation_scale = p.FFmtRead();
946  if (p.lgEOL())
947  p.NoNumb("H2 formation scale");
948 
949  /* a number <= 0 is the log of the ratio */
950  if (hmi.H2_formation_scale <= 0.)
951  hmi.H2_formation_scale = exp10(hmi.H2_frac_abund_set);
952  }
953 #endif
954  }
955 
956  /* this is a scale factor that changes the n(H0)*1.7e-4 that is added to the
957  * electron density to account for collisions with atomic H. it is an order of
958  * magnitude guess, so this command provides ability to see whether it affects results */
959  else if (p.nMatch("HCOR"))
960  {
961  dense.HCorrFac = (realnum) p.FFmtRead();
962  if (p.lgEOL())
963  p.NoNumb("scale for H0 correction to e- collision rate");
964  }
965 
966  else if (p.nMatch(" PAH"))
967  {
968  /* set one of several possible PAH abundance distribution functions */
969  if (lgQuotesFound)
970  {
971  /* abundance depends specified species as fraction of Htot */
972  gv.chPAH_abundance = chString_quotes_lowercase;
973  }
974  else if (p.nMatch("CONS"))
975  {
976  /* constant PAH abundance */
977  gv.chPAH_abundance = "CON";
978  }
979  else if (p.nMatch("BAKE"))
980  {
981  /* turn on simple PAH heating from Bakes & Tielens - this is very approximate */
982  /*>>>refer PAH heating Bakes, E.L.O., & Tielens, A.G.G.M. 1994, ApJ, 427, 822 */
983  gv.lgBakesPAH_heat = true;
984  /* warn that this model is not the best we can do: energy conservation is violated */
985  phycon.lgPhysOK = false;
986  }
987  else
988  {
989  fprintf(ioQQQ,
990  " a string, or one of the keywords BAKES, or CONStant must appear.\n Sorry.");
992  }
993  }
994 
995  /* something to do with pressure */
996  else if (p.nMatch("PRES"))
997  {
998  /* tolerance on pressure convergence */
999  if (p.nMatch("CONV"))
1000  {
1001  /* keyword is tolerance
1002  * set tolerance or relative error in pressure match */
1004  if (p.lgEOL())
1005  p.NoNumb("pressure convergence tolerance");
1006 
1007  if (conv.PressureErrorAllowed < 0.)
1009  }
1010  else if( p.nMatch("IONI") )
1011  {
1012  /* set limit to number of calls from pressure to ionize solver,
1013  * this set limit to conv.nPres2Ioniz */
1014  conv.limPres2Ioniz = (long) p.FFmtRead();
1015  if (p.lgEOL())
1016  {
1017  p.NoNumb("number of calls from pressure to ion solver");
1018  }
1019  else if (conv.limPres2Ioniz <= 0)
1020  {
1021  fprintf(ioQQQ, " The limit must be greater than zero.\n Sorry.");
1023  }
1024  }
1025 
1026  else
1027  {
1028  /* >>chng 04 mar 02, printout had been wrong, saying TOLErange
1029  * rather than CONVergence. Nick Abel */
1030  fprintf(ioQQQ,
1031  " I didn\'t recognize a key on this SET PRESSURE line.\n");
1032  fprintf(ioQQQ, " The ones I know about are: CONVergence and IONIze.\n");
1034  }
1035  }
1036  else if (p.nMatch("RECOMBIN"))
1037  {
1038  /* dielectronic recombination */
1039  if (p.nMatch("DIELECTR"))
1040  {
1041  if (p.nMatch("MEAN"))
1042  {
1043  /* change various factors for the dielectronic recombination means */
1044  if (p.nMatch(" OFF"))
1045  {
1046  for (int ion = 0; ion < LIMELM; ++ion)
1047  ionbal.DR_mean_scale[ion] = 0.;
1048  }
1049  /* multiply guess by scale factor */
1050  else if (p.nMatch("SCALE"))
1051  {
1052  // there must be at least one scale factor
1053  ionbal.DR_mean_scale[0] = (realnum) p.FFmtRead();
1054  if (p.lgEOL())
1055  {
1056  fprintf(
1057  ioQQQ,
1058  " There must be at least one scale factor on the SET RECOMBIANTION MEAN command.\n");
1060  }
1061 
1062  for (int ion = 1; ion < LIMELM; ++ion)
1063  {
1064  ionbal.DR_mean_scale[ion] = (realnum) p.FFmtRead();
1065  if (p.lgEOL())
1066  ionbal.DR_mean_scale[ion]
1067  = ionbal.DR_mean_scale[ion - 1];
1068  }
1069  for (int ion = 0; ion < LIMELM; ++ion)
1070  {
1071  if (ionbal.DR_mean_scale[ion] < 0)
1072  {
1073  fprintf(ioQQQ,
1074  " All scale factors on the SET RECOMBIANTION MEAN command must be >=0.\n");
1076  }
1077  }
1078  }
1079  /* option to add log normal noise */
1080  else if (p.nMatch("NOISE"))
1081  {
1083  if (p.lgEOL())
1084  ionbal.guess_noise = 2.;
1085  /* Seed the random-number generator with current time so that
1086  * the numbers will be different every time we run. */
1087  init_genrand((unsigned) time(NULL));
1088  }
1089  else
1090  {
1091  fprintf(ioQQQ, " key OFF or NOISE must appear.\n");
1093  }
1094  }
1095  else if( p.nMatch("SUPP") && p.nMatch(" OFF") )
1096  {
1097  ionbal.lgDRsup = false;
1098  }
1099  else
1100  {
1101  fprintf(ioQQQ, " key MEAN or SUPPression must appear.\n");
1103  }
1104  }
1105  else
1106  {
1107  fprintf(ioQQQ,
1108  " key DIELECTRonic must appear on set recombination command.\n");
1110  }
1111  }
1112 
1113  else if (p.nMatch(" DR "))
1114  {
1115  /* set zone thickness by forcing drmax and drmin */
1116  /* at this stage linear, but default is log */
1117  radius.sdrmax = p.FFmtRead();
1118  if (!p.nMatch("LINE"))
1119  {
1120  /* linear was not on command, so default is log */
1122  }
1123  if (p.lgEOL())
1124  {
1125  p.NoNumb("zone thickness");
1126  }
1127 
1128  radius.lgSdrmaxRel = p.nMatch("RELA");
1129 
1130  /* NB these being equal are tested in convinittemp to set dr */
1131  radius.lgFixed = true;
1134  if (!radius.lgSdrmaxRel && radius.sdrmax < DEPTH_OFFSET * 1e4)
1135  {
1136  fprintf(
1137  ioQQQ,
1138  "\n Thicknesses less than about %.0e will NOT give accurate results. If tricking the code\n",
1139  DEPTH_OFFSET * 1e4);
1140  fprintf(
1141  ioQQQ,
1142  " into computing emissivities instead of intensities, try to instead use a thickness of unity,\n");
1143  fprintf(
1144  ioQQQ,
1145  " and then multiply (divide) the results by the necessary thickness (product of densities).\n\n");
1147  }
1148  if (radius.lgSdrmaxRel && (radius.sdrmax <= 0. || radius.sdrmax >= 1.))
1149  {
1150  fprintf(
1151  ioQQQ,
1152  " When using a relative dr, a fraction between 0 and 1 must be entered. Found: %g\n",
1153  radius.sdrmax);
1155  }
1156  }
1157 
1158  else if (p.nMatch("DRMA"))
1159  {
1160  /* set maximum zone thickness" */
1161  radius.sdrmax = p.FFmtRead();
1162  if (p.lgEOL())
1163  p.NoNumb("maximum zone thickness");
1164 
1165  if ((radius.sdrmax < 38. || p.nMatch(" LOG")) && !p.nMatch("LINE"))
1167 
1168  // option to set min relative to current radius
1169  radius.lgSdrmaxRel = p.nMatch("RELA");
1170  if (radius.lgSdrmaxRel && (radius.sdrmax <= 0. || radius.sdrmax >= 1.))
1171  {
1172  fprintf(
1173  ioQQQ,
1174  " When using a relative drmax, a fraction between 0 and 1 must be entered. Found: %g\n",
1175  radius.sdrmax);
1177  }
1178  }
1179 
1180  else if (p.nMatch("DRMI") && p.nMatch("DEPT"))
1181  {
1182  /* option to set minimum zone thickness relative to depth */
1184  if (p.lgEOL())
1185  p.NoNumb("minimum zone thickness rel to depth");
1186 
1187  if ((radius.sdrmin_rel_depth < 38. || p.nMatch(" LOG")) && !p.nMatch("LINE"))
1189 
1190  if( radius.sdrmin_rel_depth <= 0. || radius.sdrmin_rel_depth >= 1.)
1191  {
1192  fprintf( ioQQQ, " When using a relative drmin, a fraction between 0 and 1 must be entered. Found: %g\n",
1195  }
1196  }
1197 
1198  else if (p.nMatch("DRMI"))
1199  {
1200  /* option to set minimum zone thickness */
1201  radius.sdrmin = p.FFmtRead();
1202  if (p.lgEOL())
1203  p.NoNumb("minimum zone thickness");
1204 
1205  if ((radius.sdrmin < 38. || p.nMatch(" LOG")) && !p.nMatch("LINE"))
1207 
1208  // option to set min relative to current radius
1209  radius.lgSdrminRel = p.nMatch("RELA");
1210  radius.lgSMinON = true;
1211  if (radius.lgSdrminRel && (radius.sdrmin <= 0. || radius.sdrmin >= 1.))
1212  {
1213  fprintf( ioQQQ, " When using a relative drmin, a fraction between 0 and 1 must be entered. Found: %g\n",
1214  radius.sdrmin);
1216  }
1217  }
1218 
1219  else if (p.nMatch("FLXF"))
1220  {
1221  /* faintest continuum flux to consider */
1222  rfield.FluxFaint = p.FFmtRead();
1223  if (p.lgEOL())
1224  {
1225  p.NoNumb("faintest continuum flux to consider");
1226  }
1227  if (rfield.FluxFaint < 0.)
1228  {
1230  }
1231  }
1232 
1233  else if (p.nMatch("LINE") && p.nMatch("PREC"))
1234  {
1235  fprintf(ioQQQ,"Sorry, this command changed to PRINT LINE PRECISION.\n");
1237  }
1238 
1239  /* >>chng 00 dec 08, command added by Peter van Hoof */
1240  else if (p.nMatch("NFNU"))
1241  {
1242  if (p.nMatch(" ADD"))
1243  {
1244  /* option to add a continuum point */
1245  double energy = p.FFmtRead();
1246  if (p.lgEOL())
1247  p.NoNumb("energy");
1248  t_PredCont::Inst().add(energy, p.StandardEnergyUnit());
1249  }
1250  else
1251  {
1252  /* set nFnu [incident_reflected] [incident_transmitted] [diffuse_inward] [diffuse_outward]
1253  * command for specifying what to include in the nFnu entries in LineSave.lines */
1254  /* >>chng 01 nov 12, also accept form with space */
1255  /* "incident reflected" keyword */
1256  prt.lgSourceReflected = p.nMatch("INCIDENT R") || p.nMatch(
1257  "INCIDENT_R");
1258  /* "incident transmitted" keyword */
1259  prt.lgSourceTransmitted = p.nMatch("INCIDENT_T") || p.nMatch(
1260  "INCIDENT T");
1261  /* "diffuse inward" keyword */
1262  prt.lgDiffuseInward = p.nMatch("DIFFUSE_I")
1263  || p.nMatch("DIFFUSE I");
1264  /* "diffuse outward" keyword */
1265  prt.lgDiffuseOutward = p.nMatch("DIFFUSE_O") || p.nMatch(
1266  "DIFFUSE O");
1267 
1268  /* at least one of these needs to be set ! */
1271  {
1272  fprintf(ioQQQ,
1273  " set nFnu expects one or more of the following keywords:\n");
1274  fprintf(ioQQQ, " INCIDENT_REFLECTED, INCIDENT_TRANSMITTED, "
1275  "DIFFUSE_INWARD, DIFFUSE_OUTWARD\n");
1277  }
1278  }
1279  }
1280 
1281  else if (p.nMatch("IND2"))
1282  {
1283  if (p.nMatch(" ON "))
1284  {
1285  /* set flag saying to off or on induced two photon processes */
1286  iso_ctrl.lgInd2nu_On = true;
1287  }
1288  else if( p.nMatch(" OFF") )
1289  {
1290  iso_ctrl.lgInd2nu_On = false;
1291  }
1292  else
1293  {
1294  fprintf( ioQQQ, " set ind2 needs either ON or OFF.\n" );
1296  }
1297  }
1298 
1299  else if (p.nMatch("SPECIES"))
1300  {
1301  /* Set the default collision strength for dBase transitions
1302  * without collision or radiative data */
1303  if (p.nMatch(" GBAR"))
1304  {
1305  atmdat.collstrDefault = (long)p.FFmtRead();
1306  if (p.lgEOL())
1307  {
1308  p.NoNumb("the default collision strengths when no collision or radiative data are available");
1309  }
1310  }
1311 
1312  /* Set the pseudo-continuum range and resolution for given species. */
1313  else if( p.nMatch("CONT") )
1314  {
1315  adjPseudoCont thisCont;
1316 
1317  if( lgQuotesFound )
1318  {
1319  string species = chString_quotes_lowercase;
1320  trimTrailingWhiteSpace( species );
1321  thisCont.speciesLabel = species;
1322  }
1323  else
1324  {
1325  fprintf( ioQQQ, "PROBLEM No species specified with the "
1326  "'set species continuum command'\n" );
1327  cdEXIT( EXIT_FAILURE );
1328  }
1329 
1330  /* set lower, upper wavelength range, and number of bins, for
1331  * save feii continuum command */
1332  thisCont.wlLo = (realnum) p.FFmtRead();
1333  thisCont.wlHi = (realnum) p.FFmtRead();
1334  thisCont.nBins = (long) p.FFmtRead();
1335 
1336  if (p.lgEOL())
1337  {
1338  fprintf(ioQQQ,
1339  "PROBLEM The 'set species continuum' command must have"
1340  " three numbers, the lower and upper wavelength range in Angstroms"
1341  " and the number of bins to divide this into.\n");
1343  }
1344 
1345  if( thisCont.wlLo >= thisCont.wlHi )
1346  {
1347  fprintf(ioQQQ,
1348  "PROBLEM The first two numbers on the 'set "
1349  "species continuum' command must be the lower and upper "
1350  "wavelength range in Angstroms and the first must be less "
1351  "than the second.\n");
1353  }
1354  if( thisCont.nBins < 2 )
1355  {
1356  fprintf(ioQQQ,
1357  "PROBLEM The third number on the 'set species continuum' "
1358  "command must be the number of bins to divide the range into and"
1359  " it must be greater than 1.\n");
1361  }
1362 
1363  save.setPseudoCont.push_back( thisCont );
1364  {
1365  enum { DEBUG_SPEC = false };
1366  if( DEBUG_SPEC )
1367  {
1368  fprintf( ioQQQ, "species :\t %s\n", thisCont.speciesLabel.c_str() );
1369  fprintf( ioQQQ, "conwlLo :\t %g\n", thisCont.wlLo );
1370  fprintf( ioQQQ, "conwlHi :\t %g\n", thisCont.wlHi );
1371  fprintf( ioQQQ, "ncon :\t %ld\n", thisCont.nBins );
1372  }
1373  }
1374  }
1375 
1376  else
1377  {
1378  fprintf(ioQQQ, " SET SPECIES takes options GBAR and CONTINUUM.\n");
1380  }
1381  }
1382 
1383  else if (p.nMatch("TEMP"))
1384  {
1385  /* set something to do with temperature, currently solver, tolerance, floor */
1386  if (p.nMatch("FLOO"))
1387  {
1388  StopCalc.TeFloor = (realnum) p.FFmtRead();
1389  if (p.lgEOL())
1390  p.NoNumb("temperature floor");
1391 
1392  /* linear option */
1393  if (p.nMatch(" LOG") || (StopCalc.TeFloor <= 10. && !p.nMatch(
1394  "LINE")))
1396 
1398  {
1399  fprintf(ioQQQ, " TE < %gK, reset to %gK.\n",
1402  }
1403 
1405  {
1406  fprintf(ioQQQ,
1407  " TE > %gK. Cloudy cannot handle this. Bailing out.\n",
1410  }
1411 
1412  /* option to vary scale factor */
1413  if (optimize.lgVarOn)
1414  {
1416  strcpy(optimize.chVarFmt[optimize.nparm],
1417  "SET TEMPERATURE FLOOR %f LOG");
1418 
1419  /* pointer to where to write */
1421 
1422  /* vary log of temperature */
1423  optimize.vparm[0][optimize.nparm] = (realnum) log10(
1424  StopCalc.TeFloor);
1425  optimize.varang[optimize.nparm][0] = (realnum) log10(
1426  1.00001 * phycon.TEMP_LIMIT_LOW);
1427  optimize.varang[optimize.nparm][1] = (realnum) log10(
1428  0.99999 * phycon.TEMP_LIMIT_HIGH);
1429  optimize.vincr[optimize.nparm] = 0.3f;
1430 
1431  ++optimize.nparm;
1432  }
1433  }
1434 
1435  else if (p.nMatch("CONV") || p.nMatch("TOLE"))
1436  {
1437  /* error tolerance in heating cooling match, number is error/total */
1439  if (p.lgEOL())
1440  {
1441  p.NoNumb("heating cooling tolerance");
1442  }
1443  if (conv.HeatCoolRelErrorAllowed <= 0.)
1444  {
1446  }
1447  }
1448 
1449  else
1450  {
1451  fprintf(ioQQQ,
1452  "\nI did not recognize a keyword on this SET TEMPERATURE command.\n");
1453  p.PrintLine(ioQQQ);
1454  fprintf(ioQQQ, "The keywords are FLOOr and CONVergence.\n");
1456  }
1457  }
1458 
1459  else if (p.nMatch("TEST"))
1460  {
1461  /* set flag saying to turn on some test - this is in cddefines.h in the global namespace */
1462  lgTestCodeEnabled = true;
1463  }
1464 
1465  else if (p.nMatch("TRIM"))
1466  {
1467  /* set trim upper or lower, for ionization stage trimming
1468  * ion trimming ionization trimming
1469  * in routine TrimStage */
1470  if (p.nMatch("UPPE"))
1471  {
1472  /* set trim upper - either set value or turn off */
1473  ionbal.trimhi = exp10(p.FFmtRead());
1474  ionbal.lgTrimhiOn = true;
1475  if (p.lgEOL() && p.nMatch(" OFF"))
1476  {
1477  /* turn off upward trimming */
1478  p.setEOL(false);
1479  ionbal.lgTrimhiOn = false;
1480  /* reset high limit to proper value */
1482  }
1483  }
1484 
1485  else if (p.nMatch("LOWE"))
1486  {
1487  /* set trim lower */
1488  ionbal.trimlo = exp10(p.FFmtRead());
1489  ionbal.lgTrimloOn = true;
1490  if (p.lgEOL() && p.nMatch(" OFF"))
1491  {
1492  /* turn off upward trimming */
1493  p.setEOL(false);
1494  ionbal.lgTrimloOn = false;
1495  /* reset low limit to proper value */
1497  }
1498  }
1499 
1500  /* turn off ionization stage trimming */
1501  else if (p.nMatch("SMAL") || p.nMatch(" OFF"))
1502  {
1503  /* set small limits to both upper and lower limit*/
1504  ionbal.lgTrimhiOn = false;
1505  ionbal.lgTrimloOn = false;
1508  p.setEOL(false);
1509  }
1510 
1511  else
1512  {
1513  /* set trim upper */
1514  ionbal.lgTrimhiOn = true;
1515  ionbal.trimhi = exp10(p.FFmtRead());
1516 
1517  /* set trim lower to same number */
1519  ionbal.lgTrimloOn = true;
1520  }
1521 
1522  if (p.lgEOL())
1523  {
1524  p.NoNumb("trimming parameter");
1525  }
1526 
1527  if (ionbal.trimlo >= 1. || ionbal.trimhi >= 1.)
1528  {
1529  fprintf(ioQQQ, " number must be negative since log\n");
1531  }
1532  }
1533 
1534  else if (p.nMatch("SKIP"))
1535  {
1536  /* skip save command, for saving every n't point */
1537  save.ncSaveSkip = (long) p.FFmtRead();
1538  if (p.lgEOL())
1539  {
1540  p.NoNumb("number of points to skip in save");
1541  }
1542  }
1543 
1544  else if (p.nMatch(" UTA"))
1545  {
1546  if (p.nMatch("KISI"))
1547  {
1548  if (p.nMatch(" OFF"))
1549  {
1550  /* the default, do not use it */
1552  }
1553  else if( p.nMatch( " ON " ) )
1554  {
1556  }
1557  else
1558  {
1559  fprintf( ioQQQ, "Error: SET UTA KISIELIUS expects ON or OFF\n" );
1560  cdEXIT( EXIT_FAILURE );
1561  }
1562  }
1563  else if( p.nMatch( " OFF" ) )
1564  {
1565  /* turn off ALL inner shell absorption ionization */
1566  atmdat.lgInnerShellLine_on = false;
1567  phycon.lgPhysOK = false;
1568  }
1569  }
1570 
1571  else if (p.nMatch("WEAKH"))
1572  {
1573  /* set WeakHeatCool, threshold on save heating and cooling, default 0.05 */
1575 
1576  if (p.lgEOL())
1577  {
1578  p.NoNumb("threshold on save heating and cooling");
1579  }
1580 
1581  if (save.WeakHeatCool < 0.)
1582  {
1584  = exp10(save.WeakHeatCool);
1585  }
1586  }
1587 
1588  else if (p.nMatch("KSHE"))
1589  {
1590  /* upper limit to opacities for k-shell ionization */
1592  if (p.lgEOL())
1593  {
1594  p.NoNumb("k-shell ionization opacity limit");
1595  }
1596 
1597  if (continuum.EnergyKshell == 0.)
1598  {
1599  /* arg of 0 causes upper limit to energy array */
1601  }
1602 
1603  else if (continuum.EnergyKshell < 194.)
1604  {
1605  fprintf(ioQQQ, " k-shell energy must be greater than 194 Ryd\n");
1607  }
1608  }
1609 
1610  else if (p.nMatch("NCHR"))
1611  {
1612  /* option to set the number of charge states for grain model */
1613  double val = p.FFmtRead();
1614 
1615  if (p.lgEOL())
1616  {
1617  p.NoNumb("number of charge states");
1618  }
1619  else
1620  {
1621  long nChrg = nint(val);
1622  if (nChrg < 2 || nChrg > NCHU)
1623  {
1624  fprintf(ioQQQ,
1625  " illegal value for number of charge states: %ld\n",
1626  nChrg);
1627  fprintf(ioQQQ, " choose a value between 2 and %d\n", NCHU);
1628  fprintf(ioQQQ,
1629  " or increase NCHS in grainvar.h and recompile\n");
1631  }
1632  else
1633  {
1634  SetNChrgStates(nChrg);
1635  }
1636  }
1637  }
1638 
1639  else if (p.nMatch("NEGO"))
1640  {
1641  /* save negative opacities if they occur, set negopac */
1642  opac.lgNegOpacIO = true;
1643  }
1644 
1645  else if (p.nMatch("NEND"))
1646  {
1647  /* default limit to number of zones to be computed
1648  * only do this if nend is NOT currently left at its default
1649  * nend is set to nEndDflt in routine zero
1650  * this command only has effect if stop zone not entered */
1651  if (iterations.lgEndDflt)
1652  {
1653  /* this is default limit to number of zones, change it to this value */
1654  iterations.nEndDflt = (long) p.FFmtRead();
1655  iterations.lgEndDflt = false;
1656  if (p.lgEOL())
1657  p.NoNumb("limit to zone number");
1658 
1659  /* now change all limits, for all iterations, to this value */
1660  for (long i = 0; i < iterations.iter_malloc; i++)
1662 
1663  if (iterations.nEndDflt > 2000)
1664  fprintf(
1665  ioQQQ,
1666  "CAUTION - it will take a lot of memory to save"
1667  " results for %li zones. Is this many zones really necessary?\n",
1669  }
1670  }
1671 
1672  else if (p.nMatch("TSQD"))
1673  {
1674  /* upper limit for highest density considered in the
1675  * Peimbert-style t^2 section of the printout */
1676  peimbt.tsqden = (realnum) p.FFmtRead();
1677 
1678  if (p.lgEOL())
1679  {
1680  p.NoNumb("highest density in t^2 section of printout");
1681  }
1683  }
1684 
1685  else if (p.nMatch("NMAP"))
1686  {
1687  /* how many steps in plot or save of heating-cooling map */
1688  hcmap.nMapStep = (long) p.FFmtRead();
1689  if (p.lgEOL())
1690  {
1691  p.NoNumb("steps in heating-cooling map");
1692  }
1693  }
1694 
1695  else if (p.nMatch("NUME") && p.nMatch("DERI"))
1696  {
1697  /* this is an option to use numerical derivatives for heating and cooling */
1698  NumDeriv.lgNumDeriv = true;
1699  }
1700 
1701  else if (p.nMatch("PATH"))
1702  {
1703  fprintf(ioQQQ, " The SET PATH command is no longer supported.\n");
1704  fprintf(ioQQQ, " Please set the correct path using the environment variable CLOUDY_DATA_PATH.\n");
1706  }
1707 
1708  else if (p.nMatch("PHFI"))
1709  {
1710  /* which version of PHFIT to use, 1995 or 1996 */
1711  ip = (long) p.FFmtRead();
1712  if (p.lgEOL())
1713  {
1714  p.NoNumb("version of PHFIT");
1715  }
1716 
1717  if (ip == 1995)
1718  {
1719  /* option to go back to old results, pre op */
1721  }
1722  else if (ip == 1996)
1723  {
1724  /* default is to use newer results, including opacity project */
1726  }
1727  else
1728  {
1729  fprintf(ioQQQ, " Two possible values are 1995 and 1996.\n");
1731  }
1732  }
1733 
1734  /* set save xxx command */
1735  else if (p.nMatch("PUNC") || p.nMatch("SAVE"))
1736  {
1737  if (p.nMatch("HASH"))
1738  {
1739  /* to use the hash option there must be double quotes on line - were there? */
1740  if (!lgQuotesFound)
1741  {
1742  fprintf(ioQQQ,
1743  " I didn\'t recognize a key on this SET SAVE HASH line.\n");
1745  }
1746  /* specify the terminator between save output sets - normally a set of hash marks */
1747  /*
1748  * get any string within double quotes, and return it as
1749  * null terminated string
1750  * also sets name in OrgCard and chCard to blanks so that
1751  * do not trigger off it later */
1752  if (strcmp(chString_quotes_lowercase.c_str(), "return") == 0)
1753  {
1754  /* special case - return becomes new line */
1755  sprintf(save.chHashString, "\n");
1756  }
1757  else if (strcmp(chString_quotes_lowercase.c_str(), "time") == 0)
1758  {
1759  /* special case where output time between iterations */
1760  sprintf(save.chHashString, "TIME_DEP");
1761  }
1762  else
1763  {
1764  /* usual case, simply copy what is in quotes */
1765  strcpy(save.chHashString, chString_quotes_lowercase.c_str());
1766  }
1767  }
1768 
1769  else if ((p.nMatch("LINE") && p.nMatch("WIDT")) || p.nMatch("RESO")
1770  || p.nMatch("LWID"))
1771  {
1772  /* set spectral resolution for contrast in continuum plots */
1773  if (p.nMatch(" C "))
1774  {
1775  fprintf(ioQQQ, " the keyword C is no longer necessary since"
1776  " energy conservation is now supported by default.\n");
1778  }
1779  else if (p.nMatch("SUPP"))
1780  /* suppress emission lines in save output */
1781  save.Resolution = realnum(0.);
1782  else
1783  {
1784  double number = p.FFmtRead();
1785  if (p.lgEOL())
1786  p.NoNumb("line width or resolution");
1787  if (number <= 0.)
1788  {
1789  fprintf(ioQQQ,
1790  " line width or resolution must be greater than zero.\n");
1792  }
1793 
1794  if (p.nMatch("RESO"))
1795  /* resolving power R = lambda/delta(lambda) */
1796  save.Resolution = realnum(number);
1797  else
1798  /* resolution FWHM in km/s */
1799  save.Resolution = realnum(SPEEDLIGHT / (number * 1.e5));
1800  }
1801 
1802  // option to do exactly the same thing with absorption lines
1803  // keyword is absorption
1804  if (p.nMatch("ABSO"))
1806  }
1807 
1808  else if (p.nMatch("PREF"))
1809  {
1810  if( save.nsave > 0 )
1811  {
1812  fprintf(ioQQQ, " The SET SAVE PREFIX command should precede all save commands.\n" );
1814  }
1815  // specify a prefix before all save filenames
1816  save.chFilenamePrefix = chString_quotes_lowercase;
1817  }
1818 
1819  else if (p.nMatch("FLUS"))
1820  {
1821  /* flush the output after every iteration */
1822  save.lgFLUSH = true;
1823  }
1824 
1825  else if (p.nMatch("LUMI") && p.nMatch(" OLD") )
1826  {
1827  /* use the old style of luminosity per inner cloud area in save continuum */
1828  save.lgLuminosityOld = true;
1829  }
1830 
1831  else
1832  {
1833  fprintf(ioQQQ,
1834  " There should have been an option on this command.\n");
1835  fprintf(ioQQQ,
1836  " Valid options for SET SAVE are summarized in Hazy 1 Miscellaneous Commands, and are:\n"
1837  " HASH, LINEWIDTH, RESOLUTION, PREFIX, FLUSH, LUMINOSITY OLD.\n");
1839  }
1840  }
1841 
1842  /* set continuum .... options */
1843  else if (p.nMatch("CONT"))
1844  {
1845  if (p.nMatch("RESO"))
1846  {
1847  /* set resolution, get factor that will multiply continuum resolution that
1848  * is contained in the file continuum_mesh.ini */
1849  (void)p.FFmtRead();
1850  if (p.lgEOL())
1851  {
1852  p.NoNumb("continuum resolution scale");
1853  }
1854  // no need to do anything further, the resolution scale factor was already set in cdRead()
1855  }
1856 
1857  else if (p.nMatch("SHIE"))
1858  {
1859  /* set continuum shielding function */
1860  /* these are all possible values of rt.nLineContShield,
1861  * first is default, these are set with set continuum shielding */
1862  /*#define LINE_CONT_SHIELD_PESC 1*/
1863  /*#define LINE_CONT_SHIELD_FEDERMAN 2*/
1864  /*#define LINE_CONT_SHIELD_FERLAND 3*/
1865  if (p.nMatch("PESC"))
1866  {
1867  /* this uses an inward looking escape probability */
1869  }
1870  else if (p.nMatch("FEDE"))
1871  {
1872  /* set continuum shielding Federman,
1873  * this is the default, and uses the appendix of
1874  * >>refer RT continuum shielding Federman, S.R., Glassgold, A.E., &
1875  * >>refercon Kwan, J. 1979, ApJ, 227, 466*/
1876  /* ...with a bug in the wing opacity fixed */
1878  }
1879  else if (p.nMatch("FBUG"))
1880  {
1881  /* set continuum shielding Federman,
1882  * this is the default, and uses the appendix of
1883  * >>refer RT continuum shielding Federman, S.R., Glassgold, A.E., &
1884  * >>refercon Kwan, J. 1979, ApJ, 227, 466*/
1886  }
1887  else if (p.nMatch("FERL"))
1888  {
1890  }
1891  else if (p.nMatch("RODG"))
1892  {
1893  /* set continuum shielding
1894  * >>refer Rodgers & Williams 1974JQSRT..14..319R
1895  * as in Draine & Bertoldi */
1897  }
1898  else if (p.nMatch("INTE"))
1899  {
1900  /* set continuum shielding, using best integral
1901  * formalism available -- may be slow... */
1903  }
1904  else if (p.nMatch("NONE"))
1905  {
1906  /* turn off continuum shielding */
1907  rt.nLineContShield = 0;
1908  }
1909  else
1910  {
1911  fprintf(ioQQQ,
1912  " I didn\'t recognize a key on this SET CONTINUUM SHIELDing line.\n");
1913  fprintf(ioQQQ,
1914  " The ones I know about are: PESC, FEDErman, FERLand, RODGers and INTEgral.\n");
1916  }
1917  }
1918 
1919  else
1920  {
1921  fprintf(ioQQQ,
1922  " I didn\'t recognize a key on this SET CONTINUUM line.\n");
1923  fprintf(ioQQQ,
1924  " The ones I know about are: RESOlution and SHIEld.\n");
1926  }
1927  }
1928 
1929  else
1930  {
1931  fprintf(ioQQQ, " I didn\'t recognize a key on this SET command.\n");
1933  }
1934 
1935  return;
1936 }
long int lim_iter
Definition: iterations.h:62
realnum EdenFraction
Definition: dense.h:220
void setEOL(bool val)
Definition: parser.h:107
double emm() const
Definition: mesh.h:84
bool nMatch(const char *chKey) const
Definition: parser.h:140
t_mole_global mole_global
Definition: mole.cpp:7
long int iter_malloc
Definition: iterations.h:40
bool lgBakesPAH_heat
Definition: grainvar.h:481
bool lgStancil
Definition: mole.h:337
char chGrainFormPump
Definition: hmi.h:182
bool lgTrimhiOn
Definition: ionbal.h:90
t_atmdat atmdat
Definition: atmdat.cpp:6
double EdenErrorAllowed
Definition: conv.h:265
long int fine_opac_nresolv
Definition: rfield.h:366
double FFmtRead(void)
Definition: parser.cpp:438
double exp10(double x)
Definition: cddefines.h:1383
t_input input
Definition: input.cpp:12
t_opac opac
Definition: opacity.cpp:5
realnum ResolutionAbs
Definition: save.h:478
t_hyperfine hyperfine
Definition: hyperfine.cpp:5
long int nvfpnt[LIMPAR]
Definition: optimize.h:198
string chFilenamePrefix
Definition: save.h:413
const realnum SMALLFLOAT
Definition: cpu.h:246
t_isoCTRL iso_ctrl
Definition: iso.cpp:9
t_cpu_i & i()
Definition: cpu.h:415
bool lgLeiden_Keep_ipMH2s
Definition: hmi.h:219
bool lgCheckEnergyEveryZone
Definition: continuum.h:102
vector< long int > nend
Definition: iterations.h:71
char chPresMode[20]
Definition: dynamics.h:120
realnum HCorrFac
Definition: dense.h:121
double ShockMach
Definition: dynamics.h:127
void ParseSet(Parser &p)
Definition: parse_set.cpp:38
t_StopCalc StopCalc
Definition: stopcalc.cpp:7
t_conv conv
Definition: conv.cpp:5
int GetQuote(string &chLabel)
Definition: parser.cpp:184
realnum drChange
Definition: radius.h:188
realnum varang[LIMPAR][2]
Definition: optimize.h:201
long int nRead
Definition: input.h:62
long int limPres2Ioniz
Definition: conv.h:154
t_phycon phycon
Definition: phycon.cpp:6
char chHashString[INPUT_LINE_LENGTH]
Definition: save.h:398
bool lgEndDflt
Definition: iterations.h:65
double trimhi
Definition: ionbal.h:83
realnum GasPhaseAbundErrorAllowed
Definition: conv.h:283
bool lgDiffuseInward
Definition: prt.h:209
void setAssertAbort(bool val)
Definition: cpu.h:364
double sdrmax
Definition: radius.h:158
CollIonRC CIRCData
Definition: atmdat.h:443
char chVarFmt[LIMPAR][FILENAME_PATH_LENGTH_2]
Definition: optimize.h:267
bool lgProtElim
Definition: mole.h:347
bool lgFLUSH
Definition: save.h:401
FILE * ioQQQ
Definition: cddefines.cpp:7
realnum vparm[LIMEXT][LIMPAR]
Definition: optimize.h:192
bool lgGbarOn
Definition: atmdat.h:421
double rate_h2_form_grains_set
Definition: hmi.h:192
t_dynamics dynamics
Definition: dynamics.cpp:42
bool lgConverge_set
Definition: iterations.h:60
char chJura
Definition: hmi.h:185
#define MIN2(a, b)
Definition: cddefines.h:807
Definition: parser.h:42
bool lgNegOpacIO
Definition: opacity.h:198
realnum Tad
Definition: hmi.h:135
realnum PressureErrorAllowed
Definition: conv.h:270
bool lgNumDeriv
Definition: numderiv.h:18
bool lgNonEquilChem
Definition: mole.h:342
void trimTrailingWhiteSpace(string &str)
Definition: service.cpp:155
long int nsave
Definition: save.h:303
bool lgVarOn
Definition: optimize.h:207
bool lgPhysOK
Definition: phycon.h:111
t_dense dense
Definition: global.cpp:15
static T & Inst()
Definition: cddefines.h:209
double sdrmin_rel_depth
Definition: radius.h:161
bool lgSourceReflected
Definition: prt.h:207
long int nMapStep
Definition: hcmap.h:26
const double TEMP_LIMIT_LOW
Definition: phycon.h:121
realnum SetCsupra
Definition: secondaries.h:45
long int fine_opac_nelem
Definition: rfield.h:363
t_ionbal ionbal
Definition: ionbal.cpp:8
double sdrmin
Definition: radius.h:157
realnum WeakHeatCool
Definition: save.h:470
long int nEndDflt
Definition: iterations.h:68
bool lgLeidenHack
Definition: mole.h:334
double H2_frac_abund_set
Definition: hmi.h:196
bool lgSdrminRel
Definition: radius.h:165
double energy(const genericState &gs)
bool lgDiffuseOutward
Definition: prt.h:210
realnum guess_noise
Definition: ionbal.h:226
t_continuum continuum
Definition: continuum.cpp:6
long int ncSaveSkip
Definition: save.h:466
long int nparm
Definition: optimize.h:204
const double TEMP_LIMIT_HIGH
Definition: phycon.h:123
double DR_mean_scale[LIMELM]
Definition: ionbal.h:215
t_rfield rfield
Definition: rfield.cpp:9
const int NCHU
Definition: grainvar.h:27
double collstrDefault
Definition: atmdat.h:426
realnum wlLo
Definition: save.h:195
float realnum
Definition: cddefines.h:124
const char * StandardEnergyUnit(void) const
Definition: parser.cpp:259
#define EXIT_FAILURE
Definition: cddefines.h:168
bool lgLuminosityOld
Definition: save.h:405
realnum IonizErrorAllowed
Definition: conv.h:278
bool lgSdrmaxRel
Definition: radius.h:166
bool lgEquilibrium
Definition: dynamics.h:180
double FluxFaint
Definition: rfield.h:58
#define cdEXIT(FAIL)
Definition: cddefines.h:484
NORETURN void NoNumb(const char *chDesc) const
Definition: parser.cpp:318
long int GetElem(void) const
Definition: parser.cpp:294
t_iterations iterations
Definition: iterations.cpp:6
realnum wlHi
Definition: save.h:195
t_optimize optimize
Definition: optimize.cpp:6
realnum HeatCoolRelErrorAllowed
Definition: conv.h:276
realnum EnergyKshell
Definition: continuum.h:99
t_radius radius
Definition: radius.cpp:5
realnum vincr[LIMPAR]
Definition: optimize.h:195
realnum tsqden
Definition: peimbt.h:18
long int lim_zone
Definition: iterations.h:61
long nBins
Definition: save.h:197
t_prt prt
Definition: prt.cpp:14
bool lgTestCodeEnabled
Definition: cddefines.cpp:12
bool lgLeidenCRHack
Definition: hmi.h:220
bool lgInd2nu_On
Definition: iso.h:375
realnum fine_ener_lo
Definition: rfield.h:383
bool lgPrtSciNot
t_peimbt peimbt
Definition: peimbt.cpp:5
realnum fine_ener_hi
Definition: rfield.h:383
bool lgNeutrals
Definition: mole.h:352
bool lgSetPresMode
Definition: dynamics.h:172
LyaSourceFunctionShape LyaSourceFunctionShape_assumed
Definition: hyperfine.h:72
bool lgInnerShell_Kisielius
Definition: atmdat.h:431
char chH2_small_model_type
Definition: hmi.h:179
const int LIMELM
Definition: cddefines.h:307
t_NumDeriv NumDeriv
Definition: numderiv.cpp:5
int nLineContShield
Definition: rt.h:190
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:729
bool lgSMinON
Definition: radius.h:169
double AdvecLengthInit
Definition: dynamics.h:114
realnum GrainHeatScaleFactor
Definition: grainvar.h:559
bool lgUpdateCouplings
Definition: conv.h:257
double egamry() const
Definition: mesh.h:88
bool lgInnerShellLine_on
Definition: atmdat.h:429
bool lgEOL(void) const
Definition: parser.h:103
#define MAX2(a, b)
Definition: cddefines.h:828
bool lgCollStrenThermAver
Definition: iso.h:371
double TeFloor
Definition: stopcalc.h:33
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1217
realnum ScaleJura
Definition: hmi.h:189
string chPAH_abundance
Definition: grainvar.h:504
bool lgSourceTransmitted
Definition: prt.h:208
double ShockDepth
Definition: dynamics.h:123
long int n_initial_relax
Definition: dynamics.h:132
double trimlo
Definition: ionbal.h:83
t_hcmap hcmap
Definition: hcmap.cpp:23
void init_genrand(unsigned long s)
int PrintLine(FILE *fp) const
Definition: parser.h:200
long nint(double x)
Definition: cddefines.h:762
realnum EdenSet
Definition: dense.h:214
bool lgFederman
Definition: mole.h:336
double HCTAlex
Definition: atmdat.h:321
bool lgTrimloOn
Definition: ionbal.h:93
GrainVar gv
Definition: grainvar.cpp:5
t_hmi hmi
Definition: hmi.cpp:5
t_secondaries secondaries
Definition: secondaries.cpp:5
static t_cpu cpu
Definition: cpu.h:423
double lgFixed
Definition: radius.h:159
t_save save
Definition: save.cpp:5
char chCumuType[5]
Definition: rfield.h:333
char chSolverEden[20]
Definition: conv.h:240
realnum Resolution
Definition: save.h:475
void SetNChrgStates(long nChrg)
Definition: grains.cpp:550
bool lgDRsup
Definition: ionbal.h:229
const double DEPTH_OFFSET
Definition: cddefines.h:321
void set_version(phfit_version val)
Definition: atmdat_adfa.h:50
long int nvarxt[LIMPAR]
Definition: optimize.h:198
vector< adjPseudoCont > setPseudoCont
Definition: save.h:492
string speciesLabel
Definition: save.h:194
t_rt rt
Definition: rt.cpp:5