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-2022 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  /* use new trimming approach */
1512  else if (p.nMatch("NEW"))
1513  {
1514  ionbal.lgNewTrim = true;
1515  p.setEOL(false);
1516  }
1517 
1518  /* use old trimming approach */
1519  else if (p.nMatch("NEW"))
1520  {
1521  ionbal.lgNewTrim = false;
1522  p.setEOL(false);
1523  }
1524 
1525  else
1526  {
1527  /* set trim upper */
1528  ionbal.lgTrimhiOn = true;
1529  ionbal.trimhi = exp10(p.FFmtRead());
1530 
1531  /* set trim lower to same number */
1533  ionbal.lgTrimloOn = true;
1534  }
1535 
1536  if (p.lgEOL())
1537  {
1538  p.NoNumb("trimming parameter");
1539  }
1540 
1541  if (ionbal.trimlo >= 1. || ionbal.trimhi >= 1.)
1542  {
1543  fprintf(ioQQQ, " number must be negative since log\n");
1545  }
1546  }
1547 
1548  else if (p.nMatch("SKIP"))
1549  {
1550  /* skip save command, for saving every n't point */
1551  save.ncSaveSkip = (long) p.FFmtRead();
1552  if (p.lgEOL())
1553  {
1554  p.NoNumb("number of points to skip in save");
1555  }
1556  }
1557 
1558  else if (p.nMatch(" UTA"))
1559  {
1560  if (p.nMatch("KISI"))
1561  {
1562  if (p.nMatch(" OFF"))
1563  {
1564  /* the default, do not use it */
1566  }
1567  else if( p.nMatch( " ON " ) )
1568  {
1570  }
1571  else
1572  {
1573  fprintf( ioQQQ, "Error: SET UTA KISIELIUS expects ON or OFF\n" );
1574  cdEXIT( EXIT_FAILURE );
1575  }
1576  }
1577  else if( p.nMatch( " OFF" ) )
1578  {
1579  /* turn off ALL inner shell absorption ionization */
1580  atmdat.lgInnerShellLine_on = false;
1581  phycon.lgPhysOK = false;
1582  }
1583  }
1584 
1585  else if (p.nMatch("WEAKH"))
1586  {
1587  /* set WeakHeatCool, threshold on save heating and cooling, default 0.05 */
1589 
1590  if (p.lgEOL())
1591  {
1592  p.NoNumb("threshold on save heating and cooling");
1593  }
1594 
1595  if (save.WeakHeatCool < 0.)
1596  {
1598  = exp10(save.WeakHeatCool);
1599  }
1600  }
1601 
1602  else if (p.nMatch("KSHE"))
1603  {
1604  /* upper limit to opacities for k-shell ionization */
1606  if (p.lgEOL())
1607  {
1608  p.NoNumb("k-shell ionization opacity limit");
1609  }
1610 
1611  if (continuum.EnergyKshell == 0.)
1612  {
1613  /* arg of 0 causes upper limit to energy array */
1615  }
1616 
1617  else if (continuum.EnergyKshell < 194.)
1618  {
1619  fprintf(ioQQQ, " k-shell energy must be greater than 194 Ryd\n");
1621  }
1622  }
1623 
1624  else if (p.nMatch("NCHR"))
1625  {
1626  /* option to set the number of charge states for grain model */
1627  double val = p.FFmtRead();
1628 
1629  if (p.lgEOL())
1630  {
1631  p.NoNumb("number of charge states");
1632  }
1633  else
1634  {
1635  long nChrg = nint(val);
1636  if (nChrg < 2 || nChrg > NCHU)
1637  {
1638  fprintf(ioQQQ,
1639  " illegal value for number of charge states: %ld\n",
1640  nChrg);
1641  fprintf(ioQQQ, " choose a value between 2 and %d\n", NCHU);
1642  fprintf(ioQQQ,
1643  " or increase NCHS in grainvar.h and recompile\n");
1645  }
1646  else
1647  {
1648  SetNChrgStates(nChrg);
1649  }
1650  }
1651  }
1652 
1653  else if (p.nMatch("NEGO"))
1654  {
1655  /* save negative opacities if they occur, set negopac */
1656  opac.lgNegOpacIO = true;
1657  }
1658 
1659  else if (p.nMatch("NEND"))
1660  {
1661  /* default limit to number of zones to be computed
1662  * only do this if nend is NOT currently left at its default
1663  * nend is set to nEndDflt in routine zero
1664  * this command only has effect if stop zone not entered */
1665  if (iterations.lgEndDflt)
1666  {
1667  /* this is default limit to number of zones, change it to this value */
1668  iterations.nEndDflt = (long) p.FFmtRead();
1669  iterations.lgEndDflt = false;
1670  if (p.lgEOL())
1671  p.NoNumb("limit to zone number");
1672 
1673  /* now change all limits, for all iterations, to this value */
1674  for (long i = 0; i < iterations.iter_malloc; i++)
1676 
1677  if (iterations.nEndDflt > 2000)
1678  fprintf(
1679  ioQQQ,
1680  "CAUTION - it will take a lot of memory to save"
1681  " results for %li zones. Is this many zones really necessary?\n",
1683  }
1684  }
1685 
1686  else if (p.nMatch("TSQD"))
1687  {
1688  /* upper limit for highest density considered in the
1689  * Peimbert-style t^2 section of the printout */
1690  peimbt.tsqden = (realnum) p.FFmtRead();
1691 
1692  if (p.lgEOL())
1693  {
1694  p.NoNumb("highest density in t^2 section of printout");
1695  }
1697  }
1698 
1699  else if (p.nMatch("NMAP"))
1700  {
1701  /* how many steps in plot or save of heating-cooling map */
1702  hcmap.nMapStep = (long) p.FFmtRead();
1703  if (p.lgEOL())
1704  {
1705  p.NoNumb("steps in heating-cooling map");
1706  }
1707  }
1708 
1709  else if (p.nMatch("NUME") && p.nMatch("DERI"))
1710  {
1711  /* this is an option to use numerical derivatives for heating and cooling */
1712  NumDeriv.lgNumDeriv = true;
1713  }
1714 
1715  else if (p.nMatch("PATH"))
1716  {
1717  fprintf(ioQQQ, " The SET PATH command is no longer supported.\n");
1718  fprintf(ioQQQ, " Please set the correct path using the environment variable CLOUDY_DATA_PATH.\n");
1720  }
1721 
1722  else if (p.nMatch("PHFI"))
1723  {
1724  /* which version of PHFIT to use, 1995 or 1996 */
1725  ip = (long) p.FFmtRead();
1726  if (p.lgEOL())
1727  {
1728  p.NoNumb("version of PHFIT");
1729  }
1730 
1731  if (ip == 1995)
1732  {
1733  /* option to go back to old results, pre op */
1735  }
1736  else if (ip == 1996)
1737  {
1738  /* default is to use newer results, including opacity project */
1740  }
1741  else
1742  {
1743  fprintf(ioQQQ, " Two possible values are 1995 and 1996.\n");
1745  }
1746  }
1747 
1748  /* set save xxx command */
1749  else if (p.nMatch("PUNC") || p.nMatch("SAVE"))
1750  {
1751  if (p.nMatch("HASH"))
1752  {
1753  /* to use the hash option there must be double quotes on line - were there? */
1754  if (!lgQuotesFound)
1755  {
1756  fprintf(ioQQQ,
1757  " I didn\'t recognize a key on this SET SAVE HASH line.\n");
1759  }
1760  /* specify the terminator between save output sets - normally a set of hash marks */
1761  /*
1762  * get any string within double quotes, and return it as
1763  * null terminated string
1764  * also sets name in OrgCard and chCard to blanks so that
1765  * do not trigger off it later */
1766  if (chString_quotes_lowercase == "return")
1767  {
1768  /* special case - return becomes new line */
1769  sprintf(save.chHashString, "\n");
1770  }
1771  else if (chString_quotes_lowercase == "time")
1772  {
1773  /* special case where output time between iterations */
1774  sprintf(save.chHashString, "TIME_DEP");
1775  }
1776  else
1777  {
1778  /* usual case, simply copy what is in quotes */
1779  strcpy(save.chHashString, chString_quotes_lowercase.c_str());
1780  }
1781  }
1782 
1783  else if ((p.nMatch("LINE") && p.nMatch("WIDT")) || p.nMatch("RESO")
1784  || p.nMatch("LWID"))
1785  {
1786  /* set spectral resolution for contrast in continuum plots */
1787  if (p.nMatch(" C "))
1788  {
1789  fprintf(ioQQQ, " the keyword C is no longer necessary since"
1790  " energy conservation is now supported by default.\n");
1792  }
1793  else if (p.nMatch("SUPP"))
1794  /* suppress emission lines in save output */
1795  save.Resolution = realnum(0.);
1796  else
1797  {
1798  double number = p.FFmtRead();
1799  if (p.lgEOL())
1800  p.NoNumb("line width or resolution");
1801  if (number <= 0.)
1802  {
1803  fprintf(ioQQQ,
1804  " line width or resolution must be greater than zero.\n");
1806  }
1807 
1808  if (p.nMatch("RESO"))
1809  /* resolving power R = lambda/delta(lambda) */
1810  save.Resolution = realnum(number);
1811  else
1812  /* resolution FWHM in km/s */
1813  save.Resolution = realnum(SPEEDLIGHT / (number * 1.e5));
1814  }
1815 
1816  // option to do exactly the same thing with absorption lines
1817  // keyword is absorption
1818  if (p.nMatch("ABSO"))
1820  }
1821 
1822  else if (p.nMatch("PREF"))
1823  {
1824  if( save.nsave > 0 )
1825  {
1826  fprintf(ioQQQ, " The SET SAVE PREFIX command should precede all save commands.\n" );
1828  }
1829  if( lgQuotesFound )
1830  {
1831  // specify a prefix before all save filenames
1832  save.chFilenamePrefix = chString_quotes_lowercase;
1833  }
1834  else
1835  {
1836  fprintf( ioQQQ, " PROBLEM No prefix between double quotes was found on this line.\n" );
1837  cdEXIT( EXIT_FAILURE );
1838  }
1839  }
1840 
1841  else if (p.nMatch("FLUS"))
1842  {
1843  /* flush the output after every iteration */
1844  save.lgFLUSH = true;
1845  }
1846 
1847  else if (p.nMatch("LUMI") && p.nMatch(" OLD") )
1848  {
1849  /* use the old style of luminosity per inner cloud area in save continuum */
1850  save.lgLuminosityOld = true;
1851  }
1852 
1853  else
1854  {
1855  fprintf(ioQQQ,
1856  " There should have been an option on this command.\n");
1857  fprintf(ioQQQ,
1858  " Valid options for SET SAVE are summarized in Hazy 1 Miscellaneous Commands, and are:\n"
1859  " HASH, LINEWIDTH, RESOLUTION, PREFIX, FLUSH, LUMINOSITY OLD.\n");
1861  }
1862  }
1863 
1864  /* set continuum .... options */
1865  else if (p.nMatch("CONT"))
1866  {
1867  if (p.nMatch("RESO"))
1868  {
1869  /* set resolution, get factor that will multiply continuum resolution that
1870  * is contained in the file continuum_mesh.ini */
1871  (void)p.FFmtRead();
1872  if (p.lgEOL())
1873  {
1874  p.NoNumb("continuum resolution scale");
1875  }
1876  // no need to do anything further, the resolution scale factor was already set in cdRead()
1877  }
1878 
1879  else if (p.nMatch("SHIE"))
1880  {
1881  /* set continuum shielding function */
1882  /* these are all possible values of rt.nLineContShield,
1883  * first is default, these are set with set continuum shielding */
1884  /*#define LINE_CONT_SHIELD_PESC 1*/
1885  /*#define LINE_CONT_SHIELD_FEDERMAN 2*/
1886  /*#define LINE_CONT_SHIELD_FERLAND 3*/
1887  if (p.nMatch("PESC"))
1888  {
1889  /* this uses an inward looking escape probability */
1891  }
1892  else if (p.nMatch("FEDE"))
1893  {
1894  /* set continuum shielding Federman,
1895  * this is the default, and uses the appendix of
1896  * >>refer RT continuum shielding Federman, S.R., Glassgold, A.E., &
1897  * >>refercon Kwan, J. 1979, ApJ, 227, 466*/
1898  /* ...with a bug in the wing opacity fixed */
1900  }
1901  else if (p.nMatch("FBUG"))
1902  {
1903  /* set continuum shielding Federman,
1904  * this is the default, and uses the appendix of
1905  * >>refer RT continuum shielding Federman, S.R., Glassgold, A.E., &
1906  * >>refercon Kwan, J. 1979, ApJ, 227, 466*/
1908  }
1909  else if (p.nMatch("FERL"))
1910  {
1912  }
1913  else if (p.nMatch("RODG"))
1914  {
1915  /* set continuum shielding
1916  * >>refer Rodgers & Williams 1974JQSRT..14..319R
1917  * as in Draine & Bertoldi */
1919  }
1920  else if (p.nMatch("INTE"))
1921  {
1922  /* set continuum shielding, using best integral
1923  * formalism available -- may be slow... */
1925  }
1926  else if (p.nMatch("NONE"))
1927  {
1928  /* turn off continuum shielding */
1929  rt.nLineContShield = 0;
1930  }
1931  else
1932  {
1933  fprintf(ioQQQ,
1934  " I didn\'t recognize a key on this SET CONTINUUM SHIELDing line.\n");
1935  fprintf(ioQQQ,
1936  " The ones I know about are: PESC, FEDErman, FERLand, RODGers and INTEgral.\n");
1938  }
1939  }
1940 
1941  else
1942  {
1943  fprintf(ioQQQ,
1944  " I didn\'t recognize a key on this SET CONTINUUM line.\n");
1945  fprintf(ioQQQ,
1946  " The ones I know about are: RESOlution and SHIEld.\n");
1948  }
1949  }
1950 
1951  else
1952  {
1953  fprintf(ioQQQ, " I didn\'t recognize a key on this SET command.\n");
1955  }
1956 
1957  return;
1958 }
long int lim_iter
Definition: iterations.h:62
realnum EdenFraction
Definition: dense.h:220
void setEOL(bool val)
Definition: parser.h:117
double emm() const
Definition: mesh.h:93
bool nMatch(const char *chKey) const
Definition: parser.h:150
t_mole_global mole_global
Definition: mole.cpp:7
long int iter_malloc
Definition: iterations.h:40
bool lgBakesPAH_heat
Definition: grainvar.h:479
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
bool lgNewTrim
Definition: ionbal.h:96
double EdenErrorAllowed
Definition: conv.h:263
long int fine_opac_nresolv
Definition: rfield.h:368
double FFmtRead(void)
Definition: parser.cpp:472
double exp10(double x)
Definition: cddefines.h:1368
t_input input
Definition: input.cpp:12
t_opac opac
Definition: opacity.cpp:5
realnum ResolutionAbs
Definition: save.h:496
t_hyperfine hyperfine
Definition: hyperfine.cpp:5
long int nvfpnt[LIMPAR]
Definition: optimize.h:198
string chFilenamePrefix
Definition: save.h:431
const realnum SMALLFLOAT
Definition: cpu.h:246
t_isoCTRL iso_ctrl
Definition: iso.cpp:9
t_cpu_i & i()
Definition: cpu.h:419
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:213
realnum drChange
Definition: radius.h:189
realnum varang[LIMPAR][2]
Definition: optimize.h:201
long int nRead
Definition: input.h:105
long int limPres2Ioniz
Definition: conv.h:154
t_phycon phycon
Definition: phycon.cpp:6
char chHashString[INPUT_LINE_LENGTH]
Definition: save.h:416
bool lgEndDflt
Definition: iterations.h:65
double trimhi
Definition: ionbal.h:83
realnum GasPhaseAbundErrorAllowed
Definition: conv.h:281
bool lgDiffuseInward
Definition: prt.h:209
void setAssertAbort(bool val)
Definition: cpu.h:367
double sdrmax
Definition: radius.h:159
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:419
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:803
Definition: parser.h:43
bool lgNegOpacIO
Definition: opacity.h:199
realnum Tad
Definition: hmi.h:135
realnum PressureErrorAllowed
Definition: conv.h:268
bool lgNumDeriv
Definition: numderiv.h:18
bool lgNonEquilChem
Definition: mole.h:342
void trimTrailingWhiteSpace(string &str)
Definition: service.cpp:153
long int nsave
Definition: save.h:318
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:162
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:365
t_ionbal ionbal
Definition: ionbal.cpp:8
double sdrmin
Definition: radius.h:158
realnum WeakHeatCool
Definition: save.h:488
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:166
double energy(const genericState &gs)
bool lgDiffuseOutward
Definition: prt.h:210
realnum guess_noise
Definition: ionbal.h:229
t_continuum continuum
Definition: continuum.cpp:6
long int ncSaveSkip
Definition: save.h:484
long int nparm
Definition: optimize.h:204
const double TEMP_LIMIT_HIGH
Definition: phycon.h:123
double DR_mean_scale[LIMELM]
Definition: ionbal.h:218
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:202
float realnum
Definition: cddefines.h:124
const char * StandardEnergyUnit(void) const
Definition: parser.cpp:286
#define EXIT_FAILURE
Definition: cddefines.h:168
bool lgLuminosityOld
Definition: save.h:423
realnum IonizErrorAllowed
Definition: conv.h:276
bool lgSdrmaxRel
Definition: radius.h:167
bool lgEquilibrium
Definition: dynamics.h:180
double FluxFaint
Definition: rfield.h:56
#define cdEXIT(FAIL)
Definition: cddefines.h:482
NORETURN void NoNumb(const char *chDesc) const
Definition: parser.cpp:345
long int GetElem(void) const
Definition: parser.cpp:321
t_iterations iterations
Definition: iterations.cpp:6
realnum wlHi
Definition: save.h:202
t_optimize optimize
Definition: optimize.cpp:6
realnum HeatCoolRelErrorAllowed
Definition: conv.h:274
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:204
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:385
bool lgPrtSciNot
t_peimbt peimbt
Definition: peimbt.cpp:5
realnum fine_ener_hi
Definition: rfield.h:385
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:308
t_NumDeriv NumDeriv
Definition: numderiv.cpp:5
int nLineContShield
Definition: rt.h:190
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:723
bool lgSMinON
Definition: radius.h:170
double AdvecLengthInit
Definition: dynamics.h:114
realnum GrainHeatScaleFactor
Definition: grainvar.h:557
bool lgUpdateCouplings
Definition: conv.h:255
double egamry() const
Definition: mesh.h:97
bool lgInnerShellLine_on
Definition: atmdat.h:429
bool lgEOL(void) const
Definition: parser.h:113
#define MAX2(a, b)
Definition: cddefines.h:824
bool lgCollStrenThermAver
Definition: iso.h:371
double TeFloor
Definition: stopcalc.h:33
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1121
realnum ScaleJura
Definition: hmi.h:189
string chPAH_abundance
Definition: grainvar.h:502
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:206
long nint(double x)
Definition: cddefines.h:758
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:427
double lgFixed
Definition: radius.h:160
t_save save
Definition: save.cpp:5
char chCumuType[5]
Definition: rfield.h:335
char chSolverEden[20]
Definition: conv.h:238
realnum Resolution
Definition: save.h:493
void SetNChrgStates(long nChrg)
Definition: grains.cpp:545
bool lgDRsup
Definition: ionbal.h:232
const double DEPTH_OFFSET
Definition: cddefines.h:322
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:510
string speciesLabel
Definition: save.h:201
t_rt rt
Definition: rt.cpp:5