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