cloudy trunk
Loading...
Searching...
No Matches
parse_commands.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/*ParseCommands main command line parser, decode command, then call other routines to read */
4#include "cddefines.h"
5#include "physconst.h"
6#include "optimize.h"
7#include "doppvel.h"
8#include "stopcalc.h"
9#include "abund.h"
10#include "geometry.h"
11#include "dense.h"
12#include "plot.h"
13#include "grid.h"
14#include "rfield.h"
15#include "grainvar.h"
16#include "dynamics.h"
17#include "magnetic.h"
18#include "trace.h"
19#include "atmdat.h"
20#include "h2.h"
21#include "mole.h"
22#include "hmi.h"
23#include "rt.h"
24#include "thermal.h"
25#include "opacity.h"
26#include "atomfeii.h"
27#include "called.h"
28#include "wind.h"
29#include "hextra.h"
30#include "iterations.h"
31#include "radius.h"
32#include "input.h"
33#include "monitor_results.h"
34#include "phycon.h"
35#include "fudgec.h"
36#include "version.h"
37#include "conv.h"
38#include "parse.h"
39#include "cosmology.h"
40#include "pressure.h"
41#include "parser.h"
42#include "dark_matter.h"
43
44void ParseAperture(Parser &p);
45void ParseAtom(Parser &p);
47void ParseCExtra(Parser &p);
48void ParseCMBOuter(Parser &p);
49void ParseCosm(Parser &p);
50void ParseCovering(Parser &p);
51void ParseCylinder(Parser &p);
52void ParseDarkMatter(Parser &p);
54void ParseDiffuse(Parser &p);
55void ParseDistance(Parser &p);
57void ParseEden(Parser &p);
58void ParseEnergy(Parser &p);
59void ParseFail(Parser &p);
60void ParseFill(Parser &p);
63void ParseFudge(Parser &p);
64void ParsePGrains(Parser &);
65void ParseGravity(Parser &p);
66void ParseHeLike(Parser &);
67void ParseHelp(Parser &);
68void ParseHExtra(Parser &p);
70void ParseHydrogen(Parser &);
71void ParseInitCount(Parser &p);
72void ParseIntensity(Parser &p);
73void ParseIterations(Parser &p);
74void ParseL_nu(Parser &p);
75void ParseLaser(Parser &p);
76void ParseLuminosity(Parser &p);
77void ParseNeutrons(Parser &p);
78void ParseNuF_nu(Parser &p);
79void ParseNuL_nu(Parser &p);
80void ParsePhi(Parser &p);
81void ParseQH(Parser &p);
82void ParseRoberto(Parser &);
83void ParseSpecial(Parser &);
84void ParseTauMin(Parser &p);
85void ParseTitle(Parser &);
87void ParseVLaw(Parser &p);
88void ParseTurbulence(Parser &p);
89
90void ParseCommands(void)
91{
92 bool lgStop ,
93 lgStop_not_enough_info;
94
95 DEBUG_ENTRY( "ParseCommands()" );
96
97 /* following says abundances are solar */
98 abund.lgAbnSolar = true;
99
100 /* there are no plots desired yet */
101 plotCom.lgPlotON = false;
102 plotCom.nplot = 0;
103
104 /* this flag remembers whether grains have ever been turned on. It is passed
105 * to routine ParseAbundances - there grains will be turned on with commands
106 * such as abundances ism, unless grains were previously set
107 * with a grains command. this way abundances will not clobber explicitly set
108 * grains. */
109
110 radius.Radius = 0.;
111 radius.rinner = 0.;
112 rfield.nShape = 0;
113
114 /* will be set true if no buffering command entered */
115 input.lgSetNoBuffering = false;
116
117 /* initialize some code variables in case assert command used in input stream */
119
120 for( long int i=0; i < LIMSPC; i++ )
121 {
122 strcpy( rfield.chRSpec[i], "UNKN" );
123 }
124 optimize.nparm = 0;
125
126 /* this is an option to turn on trace printout on the nth
127 * call from the optimizer */
128 if( optimize.lgTrOpt )
129 {
130 /* nTrOpt was set with "optimize trace" command,
131 * is iteration to turn on trace */
132 if( optimize.nTrOpt == optimize.nOptimiz )
133 {
134 trace.lgTrace = true;
135 called.lgTalk = cpu.i().lgMPI_talk();
136 trace.lgTrOvrd = true;
137 fprintf( ioQQQ, " READR turns on trace from optimize option.\n" );
138 }
139 }
140
141 /* say this is a beta version if we are talking and it is the truth */
142 if( t_version::Inst().nBetaVer > 0 && called.lgTalk )
143 {
144 fprintf( ioQQQ,
145 "\n This is a beta release of Cloudy, and is intended for testing only.\n" );
146 fprintf( ioQQQ,
147 "Please help make Cloudy better by posing problems or suggestions on http://tech.groups.yahoo.com/group/cloudy_simulations/.\n\n" );
148 }
149
150 if( called.lgTalk )
151 {
152 /* this code prints pretty lines at top of output box */
153 int indent = (int)((122 - strlen(t_version::Inst().chVersion))/2);
154 fprintf( ioQQQ, "%*cCloudy %s\n", indent, ' ', t_version::Inst().chVersion );
155
156 fprintf( ioQQQ, "%57cwww.nublado.org\n\n", ' ' );
157
158 /* now print box and date of version, before printing commands */
159 fprintf( ioQQQ, "%23c", ' ' );
160 fprintf( ioQQQ, "**************************************");
161 fprintf( ioQQQ, "%7.7s", t_version::Inst().chDate);
162 fprintf( ioQQQ, "**************************************\n");
163
164 fprintf( ioQQQ, "%23c*%81c*\n", ' ', ' ' );
165 }
166
167 /* read in commands and print them */
168
169 /* following signals which read is in progress,
170 * 1 is forward read, of input commands
171 * -1 is reverse read from bottom of line array, of cloudy.ini file */
172 input.iReadWay = 1;
173
174 /* initialize array reader, this sub does nothing but set
175 * initial value of a variable, depending on value of iReadWay */
176 input.init();
177
178 static const CloudyCommand commands[] = {
179 {"ABSO",ParseAbsMag},
180 /* enter luminosity in absolute magnitudes, in reads2 */
181 {"AGE",ParseAge},
182 /* enter age of cloud so we can check for time-steady reads2 */
183 {"AGN ",ParseAgn},
184 /* enter generic style AGN continuum, in reads2 */
186 {"APER",ParseAperture},
187 {"ATOM",ParseAtom},
188 {"BACK", ParseBackgrd},
189 /* cosmic background, in parse_backgrd */
190 {"BLAC", ParseBlackbody},
191 /* black body, in reads2 */
192 {"BREM", ParseBremsstrahlung},
193 {"CASE",ParseCaseB},
194 /* do Case A or Case B (usually Case B since most realistic) */
195 {"CEXT",ParseCExtra},
196 /* CMB command */
197 {"CLUM",ParseFill},
198 {"CMB", ParseCMBOuter},
199 /* cosmic thermal background radiation, argument is redshift */
200 /* if no number on line then (realnum)FFmtRead returns z=0; i.e., now */
201 {"COMP", ParseCompile},
202 /* compile ascii version of stellar atmosphere continua in volk */
203 {"CONS",ParseConstant},
204 /* constant temperature, pressure, density, or gas pressure
205 * in readsun */
206 {"CORO",ParseCoronal},
207 /* coronal equilibrium; set constant temperature to number on line
208 * in readsun */
209 {"COSMOLOGY",ParseCosmology},
210 {"COSMI", ParseCosmicRays},
211 {"COSM",ParseCosm}, // Backstop for ambiguity
212 {"COVE",ParseCovering},
213 {"CRAS",ParseCrashDo},
214 /* any of several tests to cause the code to crash */
215 {"CYLI",ParseCylinder},
216 {"DARK",ParseDarkMatter},
217 {"DIEL",ParseDielectronic},
218 {"DIFF",ParseDiffuse},
219 {"DIST",ParseDistance},
220 {"DLAW",ParseDLaw},
221 /* either use dense_fabden routine, or read in table of points */
222 {"DOUB",ParseDoubleTau},
223 /* double optical depth scale after each iteration */
224 {"DRIV",ParseDrive},
225 /* call one of several drivers, readsun */
226 {"EDEN",ParseEden},
227 /* option to add "extra" electrons, to test Compton limit
228 * for very low T(star) - option is log of eden */
229 {"ELEM",ParseElement},
230 /* element command;
231 * scale or abundance options, to change abundance of specific element
232 * read option to change order of elements
233 * in reads2.f */
234 {"ENER",ParseEnergy},
235 {"EXTI",ParseExtinguish},
236 /* extinguish ionizing continuum by absorbing column AFTER
237 * setting luminosity or Q(H). First number is the column
238 * density (log), second number is leakage (def=0%)
239 * last number is lowest energy (ryd), last two may be omitted
240 * from right to left
241 *
242 * extinction is actually done in extin, which is called by ContSetIntensity */
243 {"FAIL",ParseFail},
244 /* reset number of temp failures allowed, default=20 */
245 {"FILL",ParseFill},
246 {"FLUC",ParseFluc},
247 /* rapid density fluctuations, in readsun */
248 {"F(NU",ParseF_nuSpecific},
249 /* this is the specific flux at nu
250 * following says F(nu) not nuF(nu) */
251 {"FORC",ParseForceTemperature},
252 /* force temperature of first zone, but don't keep constant
253 * allow to then go to nearest equilibrium
254 * log if < 10 */
255 {"FUDG",ParseFudge},
256 {"GLOB",ParseGlobule},
257 /* globule with density increasing inward
258 * parameters are outer density, radius of globule, and density power */
259 {"GRAI",ParseGrain},
260 /* read parameters dealing with grains, in reads2 */
261 {"PGRA",ParsePGrains},
262 {"GRAV",ParseGravity},
263 /* (self-)Gravity forces: Yago Ascasibar (UAM, Spring 2009) */
264 {"GRID",ParseGrid},
265 /* option to run grid of models by varying certain parameters
266 * in reads2 */
267 {"HDEN",ParseHDEN},
268 /* parse the hden command to set the hydrogen density, in reads2 */
269 {"HELI",ParseHeLike},
270 {"HELP",ParseHelp},
271 {"HEXT",ParseHExtra},
272 /* "extra" heating rate, so that first= log(erg(cm-3, s-1},
273 * second optional number is scale radius, so that HXTOT = TurbHeat*SEXP(DEPTH/SCALE)
274 * if missing then constant heating.
275 * third option is depth from shielded face, to mimic irradiation from both sides*/
276 {"HIGH",ParseConvHighT},
277 /* approach equilibrium from high te */
278 {"HYDROGEN",ParseHydrogen},
279 {"ILLU",ParseIlluminate},
280 // illuminate command
281 {"INIT",ParseInitCount},
282 {"INTEN",ParseIntensity},
283 {"INTER",ParseInterp},
284 /* interpolate on input tables for continuum, set of power laws used
285 * input ordered pairs nu( ryd or log(Hz},, log(fnu)
286 * additional lines begin CONTINUE
287 * first check that this is the one and only INTERP command
288 * in readsun */
289 {"IONI",ParseIonParI},
290 /* inter ionization parameter U=Q/12 R*R N C;
291 * defined per hydrogen, not per electron (as before)
292 * radius must also be entered if spherical, not needed if plane */
293 {"ITER",ParseIterations},
294 {"L(NU",ParseL_nu},
295 {"LASE",ParseLaser},
296 {"LUMI",ParseLuminosity},
297 {"MAGN", ParseMagnet},
298 /* parse the magnetic field command, routine in magnetic.c */
299 {"MAP ",ParseMap},
300 /* do cooling space map for specified zones, in reads2 */
301 {"META",ParseMetal},
302 /* read depletion for metals, all elements heavier than He
303 * in reads2 */
304 {"MONI",ParseMonitorResults},
305 /* monitor that code predicts certain results, in monitor_results.h */
306 {"NEUT",ParseNeutrons},
307 {"NO ", ParseDont},
308 /* don't do something, in readsun */
309 {"NORM",ParseNorm},
310 /* normalize lines to this rather than h-b, sec number is scale factor */
311 {"NUF(",ParseNuF_nu},
312 {"NUL(",ParseNuL_nu},
313 {"OPTI",ParseOptimize},
314 /* option to optimize the model by varying certain parameters
315 * in reads2 */
316 {"PHI(",ParsePhi},
317 {"PLOT", ParsePlot},
318 /* make plot of log(nu.f(nu}, vs log(nu) or opacity
319 * in file plot */
320 {"POWE", ParsePowerlawContinuum},
321 /* power law with cutoff, in reads2 */
322 {"PRIN",ParsePrint},
323 /* adjust print schedule, in readsun */
324 {"PUNC",ParseSave},
325 /* save something, in save */
326 {"Q(H)",ParseQH},
327 {"RATI",ParseRatio},
328 /* enter a continuum luminosity as a ratio of
329 * nuFnu for this continuum relative to a previous continuum
330 * format; first number is ratio of second to first continuum
331 * second number is energy for this ratio
332 * if third number on line, then 2nd number is energy of
333 * first continuum, while 3rd number is energy of second continuum
334 * in reads2 */
335 {"RADI", ParseRadius},
336 /* log of inner and outer radii, default second=infinity,
337 * if R2<R1 then R2=R1+R2
338 * there is an optional keyword, "PARSEC" on the line
339 * to use PC as units, reads2 */
340 {"ROBE",ParseRoberto},
341 {"SAVE",ParseSave},
342 {"SET ",ParseSet},
343 /* set something, in reads2 */
344 {"SPEC", ParseSpecial},
345 {"SPHE", ParseSphere},
346 /* compute a spherical model, diffuse field from other side in
347 * in reads2 */
348 {"STAT",ParseState},
349 /* either get or put the code's state as a file on disk */
350 {"STOP",ParseStop},
351 /* stop model at desired zone, temperature, column density or tau-912
352 * in readsun */
353 {"TABL",ParseTable},
354 /* interpolate on input tables for continuum, set of power laws used
355 * input stored in big BLOCK data
356 * first check that this is the one and only INTERP command
357 * in readsun */
358 {"TAUM",ParseTauMin},
359 {"TEST",ParseTest},
360 /* parse the test command and its options */
361 {"TIME",ParseDynaTime},
362 /* parse the time dependent command, in dynamics.c */
363 {"TITL",ParseTitle},
364 {"TLAW", ParseTLaw},
365 /* some temperature vs depth law */
366 {"TOLE", ParseTolerance},
367 {"TRAC", ParseTrace},
368 /* turn on trace output, in reads2 */
369 {"VLAW",ParseVLaw},
370 {"TURB",ParseTurbulence},
371 {"WIND",ParseDynaWind},
372 /* NB - advection and wind commands are now a single command */
373 /* parse the wind command, in dynamics.c */
374 {"XI",ParseIonParX},
375 {NULL,NULL}}; // {NULL,NULL} sentinel must come last
376
377 Parser p(commands);
378
379 p.m_nqh = 0;
380 p.m_nInitFile=0;/* used to count number of init files read in */
381 p.m_lgDSet = false;
382 p.m_lgEOF = false;
383
384 /* read until eof or blank line, then return control back to main program */
385
386 while (p.getline())
387 {
388 /* line starting with blank is taken as end of commands */
389 if ( p.last( ) )
390 break;
391
392 /* echo the line but only if it does not contain the keyword HIDE */
393 p.echo( );
394
395 /* check whether VARY is on line */
396 if( p.nMatch("VARY") )
397 {
398 optimize.lgVarOn = true;
399 if( optimize.nparm + 1 > LIMPAR )
400 {
401 fprintf( ioQQQ, " Too many VARY lines entered; the limit is%4ld\n",
402 LIMPAR );
404 }
405 }
406 else
407 {
408 optimize.lgVarOn = false;
409 }
410
411 if( p.isCommandComment() )
412 {
413 ((void)0); // do nothing for comments
414 }
415 else if( p.isVar() )
416 {
417 p.doSetVar();
418 }
419 else
420 {
421 long int i;
422 for (i=0; commands[i].name != NULL; ++i)
423 {
424 if (p.Command(commands[i].name,commands[i].action))
425 break;
426 }
427 if (commands[i].name == NULL)
428 {
429 p.CommandError(); // No command was recognized
430 }
431 }
432 }
433
434 /*------------------------------------------------------------------- */
435 /* fall through - hit lgEOF or blank line */
436
437 if( called.lgTalk )
438 {
439 fprintf( ioQQQ, "%23c*%81c*\n", ' ', ' ' );
440 fprintf( ioQQQ, "%23c***********************************************************************************\n\n\n\n", ' ' );
441 }
442
443 /* this is an option to turn on trace printout on the nth
444 * call from the optimizer - only allow trace if
445 * this is the case and nOptimiz 1 below nTrOpt */
446 if( optimize.lgTrOpt )
447 {
448 /* nTrOpt was set with "optimize trace" command,
449 * is iteration to turn on trace */
450 if( optimize.nTrOpt != optimize.nOptimiz )
451 {
452 trace.lgTrace = false;
453 /* following overrides turning on trace elsewhere */
454 trace.lgTrOvrd = false;
455 }
456 else
457 {
458 trace.lgTrace = true;
459 called.lgTalk = cpu.i().lgMPI_talk();
460 trace.lgTrOvrd = true;
461 fprintf( ioQQQ, " READR turns on trace from optimize option.\n" );
462 }
463 }
464
465 /* set density from DLAW command, must be done here since it may depend on later input commands */
466 if( strcmp(dense.chDenseLaw,"DLW1") == 0 )
467 {
468 dense.SetGasPhaseDensity( ipHYDROGEN, (realnum)dense_fabden(radius.Radius,radius.depth) );
469
470 if( dense.gas_phase[ipHYDROGEN] <= 0. )
471 {
472 fprintf( ioQQQ, " PROBLEM DISASTER Hydrogen density set by DLAW must be > 0.\n" );
474 }
475 }
476 else if( strcmp(dense.chDenseLaw,"DLW2") == 0 )
477 {
478 dense.SetGasPhaseDensity( ipHYDROGEN, (realnum)dense_tabden(radius.Radius,radius.depth) );
479
480 if( dense.gas_phase[ipHYDROGEN] <= 0. )
481 {
482 fprintf( ioQQQ, " PROBLEM DISASTER Hydrogen density set by DLAW must be > 0.\n" );
484 }
485 }
486 else if( strcmp(dense.chDenseLaw,"DLW3") == 0 )
487 {
488 dense.SetGasPhaseDensity( ipHYDROGEN, (realnum)dense_parametric_wind(radius.Radius) );
489
490 if( dense.gas_phase[ipHYDROGEN] <= 0. )
491 {
492 fprintf( ioQQQ, " PROBLEM DISASTER Hydrogen density set by DLAW must be > 0.\n" );
494 }
495 }
496
497 /* start checks on parameters set properly - this begins with same line saying start .. */
498
499 /* lgStop_not_enough_info says that not enough info for model, so stop
500 * set true in following tests if anything missing */
501 lgStop_not_enough_info = false;
502 lgStop = false;
503
504 /* check whether hydrogen density has been set - this value was set to 0 in zero */
505 if( dense.gas_phase[ipHYDROGEN] <= 0. )
506 {
507 fprintf( ioQQQ, " PROBLEM DISASTER Hydrogen density MUST be specified.\n" );
508 lgStop_not_enough_info = true;
509 lgStop = true;
510 /* need to set something since used below - will abort
511 * since lgStop is set */
512 dense.SetGasPhaseDensity( ipHYDROGEN, 1. );
513 }
514
515 /* the SAVE XSPEC command cannot be combined with negative increments on the GRID command */
516 if( grid.lgSaveXspec && grid.lgNegativeIncrements )
517 {
518 if( called.lgTalk )
519 {
520 fprintf( ioQQQ, " PROBLEM DISASTER The SAVE XSPEC command cannot be combined with negative grid increments.\n" );
521 fprintf( ioQQQ, " PROBLEM DISASTER Please check your GRID commands.\n\n\n" );
523 }
524 }
525
526 if( geometry.covaper < 0.f || geometry.iEmissPower == 2 )
527 geometry.covaper = geometry.covgeo;
528
529 radius.rinner = radius.Radius;
530
531 /* mass flux for wind model - used for mass conservation */
532 wind.emdot = dense.gas_phase[ipHYDROGEN]*wind.windv0;
533
534 /* set converge criteria - limit number of iterations and zones */
535 if( iterations.lgConverge_set)
536 {
537 iterations.itermx = MIN2( iterations.itermx , iterations.lim_iter );
538 for( long int j=0; j < iterations.iter_malloc; j++ )
539 {
540 geometry.nend[j] = MIN2( geometry.nend[j] , iterations.lim_zone );
541 }
542 }
543
544 /* NSAVE is number of lines saved, =0 if no commands entered */
545 if( input.nSave < 0 )
546 {
547 fprintf( ioQQQ, " PROBLEM DISASTER No commands were entered - whats up?\n" );
549 }
550
551 /* iterate to convergence and wind models are mutually exclusive */
552 if( wind.lgBallistic() && conv.lgAutoIt )
553 {
554 if( called.lgTalk )
555 {
556 fprintf( ioQQQ, " NOTE PROBLEM Due to the nature of the Sobolev approximation, it makes no sense to converge a windy model.\n" );
557 fprintf( ioQQQ, " NOTE Iterate to convergence is turned off\n\n\n" );
558 }
559 conv.lgAutoIt = false;
560 iterations.itermx = 0;
561 }
562
563 /* iterate to convergence and case b do not make sense together */
564 /* WJH 22 May 2004: unless we are doing i-front dynamics (negative wind.windv0) */
565 if( opac.lgCaseB && conv.lgAutoIt && (wind.lgBallistic() || wind.lgStatic()) )
566 {
567 if( called.lgTalk )
568 {
569 fprintf( ioQQQ, " NOTE Case B is an artificial test, it makes no sense to converge this model.\n" );
570 fprintf( ioQQQ, " NOTE Iterate to convergence is turned off.\n\n\n" );
571 }
572 conv.lgAutoIt = false;
573 iterations.itermx = 0;
574 }
575
576 /* specifying a density power and constant pressure makes no sense */
577 if( dense.DensityPower!=0. && strcmp( dense.chDenseLaw, "CPRE" )==0 )
578 {
579 if( called.lgTalk )
580 {
581 fprintf( ioQQQ, " NOTE Specifying both a density power law and constant pressure is impossible.\n" );
582 }
583 lgStop = true;
584 }
585
586 if( !rfield.lgIonizReevaluate && strcmp( dense.chDenseLaw, "CDEN" )!=0 )
587 {
588 if( called.lgTalk )
589 {
590 fprintf( ioQQQ, " NOTE NO REEVALUATE IONIZATION can only be used with constant density.\n" );
591 fprintf( ioQQQ, " NOTE Resetting to reevaluate ionization.\n\n" );
592 }
593 rfield.lgIonizReevaluate = true;
594 }
595
596 if( !rfield.lgOpacityReevaluate && strcmp( dense.chDenseLaw, "CDEN" )!=0 )
597 {
598 if( called.lgTalk )
599 {
600 fprintf( ioQQQ, " NOTE NO REEVALUATE OPACITY can only be used with constant density.\n" );
601 fprintf( ioQQQ, " NOTE Resetting to reevaluate opacity.\n\n" );
602 }
603 rfield.lgOpacityReevaluate = true;
604 }
605
606 /* check that a symmetry is specified if gravity from an external mass has been added */
607 if( pressure.external_mass[0].size()>0 && pressure.gravity_symmetry==-1 )
608 {
609 if( called.lgTalk )
610 {
611 fprintf( ioQQQ, " NOTE Gravity from an external mass has been added, but no symmetry (spherical/mid-plane) was specified.\n" );
612 fprintf( ioQQQ, " NOTE It will be ignored.\n\n\n" );
613 }
614 }
615
616 /* check if the combination of stopping column density and H density are physically plausible */
617 double thickness = min( MIN3( StopCalc.tauend, StopCalc.colpls, StopCalc.colnut ),
618 MIN3( StopCalc.col_h2, StopCalc.col_h2_nut, StopCalc.HColStop ) );
619 if( thickness < COLUMN_INIT )
620 {
621 /* a stop column density was specified - check on physical thickness this corresponds to */
622 thickness /= (dense.gas_phase[ipHYDROGEN]*geometry.FillFac);
623 /* don't complain if outer radius set small with `stop thickness' command. */
624 if( thickness > 1e25 && radius.StopThickness[0] > 1e25 )
625 {
626 fprintf( ioQQQ,
627 "NOTE The specified column density and hydrogen density correspond to a thickness of %.2e cm.\n",
628 thickness);
629 fprintf( ioQQQ,
630 "NOTE This seems large to me.\n");
631 fprintf(ioQQQ,"NOTE a very large radius may cause overflow.\n\n");
632 }
633 }
634
635 if( gv.lgDColOn && thermal.ConstGrainTemp>0 && called.lgTalk )
636 {
637 /* warn if constant grain temperature but gas-grain thermal effects
638 * are still included */
639 fprintf( ioQQQ,
640 "NOTE The grain temperatures are set to a constant value with the "
641 "CONSTANT GRAIN TEMPERATURE command, but "
642 "energy exchange \n");
643 fprintf( ioQQQ,
644 "NOTE is still included. The grain-gas heating-cooling will be incorrect. "
645 "Consider turning off gas-grain collisional energy\n");
646 fprintf( ioQQQ,
647 "NOTE exchange with the NO GRAIN GAS COLLISIONAL ENERGY EXCHANGE command.\n\n\n");
648 }
649
650 if( !rfield.lgDoLineTrans && rfield.lgOpacityFine )
651 {
652 if( called.lgTalk )
653 {
654 fprintf( ioQQQ, " NOTE NO LINE TRANSER set but fine opacities still computed.\n" );
655 fprintf( ioQQQ, " NOTE Turning off fine opacities.\n\n" );
656 }
657 rfield.lgOpacityFine = false;
658 }
659
660 if( h2.lgEnabled && (!rfield.lgDoLineTrans || !rfield.lgOpacityFine) )
661 {
662 if( called.lgTalk )
663 {
664 fprintf( ioQQQ, " NOTE Large H2 molecule turned on but line transfer and fine opacities are not.\n" );
665 fprintf( ioQQQ, " NOTE Turning on line transfer and fine opacities.\n\n" );
666 }
667 rfield.lgOpacityFine = true;
668 rfield.lgDoLineTrans = true;
669 }
670
671 if( rfield.lgMustBlockHIon && !rfield.lgBlockHIon )
672 {
673 /* one of the input continua needs to have H-ionizing radiation
674 * blocked with extinguish command, but it was not done */
675 fprintf( ioQQQ, "\n NOTE\n"
676 " NOTE One of the incident continuum is a form used when no H-ionizing radiation is produced.\n" );
677 fprintf( ioQQQ, " NOTE You must also include the EXTINGUISH command to make sure this is done.\n" );
678 fprintf( ioQQQ, " NOTE The EXTINGUISH command was not included.\n" );
679 fprintf( ioQQQ, " NOTE YOU MAY BE MAKING A BIG MISTAKE!!\n NOTE\n\n\n\n" );
680 }
681
682 /* if stop temp set below default then we are going into cold and possibly molecular
683 * gas - check some parameters in this case */
684 if( called.lgTalk && (StopCalc.TempLoStopZone < phycon.TEMP_STOP_DEFAULT ||
685 /* thermal.ConstTemp def is zero, set pos when constant temperature is set */
686 (thermal.ConstTemp > 0. && thermal.ConstTemp < phycon.TEMP_STOP_DEFAULT ) ) )
687 {
688 /* print warning if temperature set below default but cosmic rays not turned on
689 * do not print if molecules are off */
690 if( (hextra.cryden == 0.) && !mole_global.lgNoMole )
691 {
692 fprintf( ioQQQ, "\n NOTE\n"
693 " NOTE The simulation is going into neutral gas but cosmic rays are not included.\n" );
694 fprintf( ioQQQ, " NOTE Ion-molecule chemistry will not occur without a source of ionization.\n" );
695 fprintf( ioQQQ, " NOTE The chemistry network may collapse deep in molecular regions.\n" );
696 fprintf( ioQQQ, " NOTE Consider adding galactic background cosmic rays with the COSMIC RAYS BACKGROUND command.\n" );
697 fprintf( ioQQQ, " NOTE You may be making a BIG mistake.\n NOTE\n\n\n\n" );
698 }
699 }
700
701 /* dense.gas_phase[ipHYDROGEN] is linear hydrogen density (cm-3) */
702 /* test for hydrogen density properly set has already been done above */
703 if( called.lgTalk && dense.gas_phase[ipHYDROGEN] < 1e-4 )
704 {
705 fprintf( ioQQQ, " NOTE Is the entered value of the hydrogen density (%.2e) reasonable?\n",
706 dense.gas_phase[ipHYDROGEN]);
707 fprintf( ioQQQ, " NOTE It seems pretty low to me.\n\n\n" );
708 }
709 else if( called.lgTalk && dense.gas_phase[ipHYDROGEN] > 1e15 )
710 {
711 fprintf( ioQQQ, " NOTE Is this value of the hydrogen density reasonable?\n" );
712 fprintf( ioQQQ, " NOTE It seems pretty high to me.\n\n\n" );
713 }
714
715 /* is the model going to crash because of extreme density? */
716 if( called.lgTalk && !lgStop && !lgStop_not_enough_info )
717 {
718 if( dense.gas_phase[ipHYDROGEN] < 1e-6 || dense.gas_phase[ipHYDROGEN] > 1e19 )
719 {
720 fprintf( ioQQQ, " NOTE Simulation may crash because of extreme "
721 "density. The value was %.2e\n\n" ,
722 dense.gas_phase[ipHYDROGEN] );
723 }
724 }
725
726 if( rfield.nShape == 0 )
727 {
728 fprintf( ioQQQ, " PROBLEM DISASTER No incident radiation field was specified - "
729 "at least put in the CMB.\n" );
730 lgStop = true;
731 lgStop_not_enough_info = true;
732 }
733
734 if( (p.m_nqh) == 0 )
735 {
736 fprintf( ioQQQ, " PROBLEM DISASTER Luminosity of continuum MUST be specified.\n" );
737 lgStop = true;
738 lgStop_not_enough_info = true;
739 }
740
741 /* Testing on 0 is safe since this it is initialized that way
742 * only print comment if continuum has been specified,
743 * if continuum not given then we are aborting anyway */
744 if( radius.Radius == 0. && rfield.nShape> 0)
745 {
746 fprintf( ioQQQ, " PROBLEM DISASTER Starting radius MUST be specified.\n" );
747 lgStop = true;
748 lgStop_not_enough_info = true;
749 }
750
751 if( rfield.nShape != (p.m_nqh) )
752 {
753 fprintf( ioQQQ, " PROBLEM DISASTER There were not the same number of continuum shapes and luminosities entered.\n" );
754 lgStop = true;
755 }
756
757 /* we only want to do this test on the first call to the command
758 * parser - it will be called many more times but with no grid command
759 * during the actual grid calculation */
760 static bool lgFirstPass = true;
761
762 /* the NO VARY option sets this flag, and can be used to turn off
763 * the grid command, as well as the optimizer */
764 if( optimize.lgNoVary && grid.lgGrid )
765 {
766 /* ignore grids */
767 grid.lgGrid = false;
768 optimize.nparm = 0;
769 grid.nGridCommands = 0;
770 }
771
772 if( lgFirstPass && grid.lgGrid && (optimize.nparm!=grid.nGridCommands) )
773 {
774 /* number of grid vary options do match */
775 fprintf( ioQQQ, " PROBLEM DISASTER The GRID command was entered "
776 "but there were %li GRID commands and %li commands with a VARY option.\n" ,
777 grid.nGridCommands , optimize.nparm);
778 fprintf( ioQQQ, " There must be the same number of GRIDs and VARY.\n" );
779 lgStop = true;
780 }
781 lgFirstPass = false;
782
783 if( lgStop_not_enough_info )
784 {
785 fprintf( ioQQQ, " PROBLEM DISASTER I do not have enough information to do the simulation, I cannot go on.\n" );
786 fprintf( ioQQQ, "\n\n Sorry.\n\n\n" );
788 }
789
790 {
791 bool lgParserTest = false;
792 if ( lgParserTest )
793 {
794 // Quit after parse phase, to allow quick verification that
795 // there are no immediate syntax handling failures.
796 fprintf(ioQQQ,"Parser phase PASSED\n");
797 exit(0);
798 }
799 }
800
801 if( lgStop )
802 {
804 }
805
806 /* end checks on parameters set properly - this begins with same line saying start .. */
807 return;
808}
809
811{
812 /* chemical abundances, readsun, p.m_lgDSet is set true if grains command
813 * ever entered, so that abundances command does not override grains command */
815 /* abundances no longer solar */
816 abund.lgAbnSolar = false;
817}
818
820{
821 DEBUG_ENTRY("ParseAperture()");
822 /* aperture command to simulate pencil beam or long slit */
823
824 /* if the "BEAM" or "SLIT" option is specified then only part
825 * of the geometry is observed, and intensities
826 * should not be weighted by r^2. There are two limiting cases, SLIT in which
827 * the slit is longer than the diameter of the nebula and the contribution to the
828 * detected luminosity goes as r^1, and BEAM when the contribution is r^0,
829 * the same as plane parallel
830 *
831 * default value of geometry.iEmissPower is 2 (set in zero.c) for full geometry
832 */
833 if( p.nMatch("SLIT") )
834 {
835 /* long slit is case where slit is longer than diameter, so emissivity contributes
836 * r^1 to the observed luminosity */
837 geometry.iEmissPower = 1;
838 }
839 else if( p.nMatch("BEAM") )
840 {
841 /* pencil beam is case where we view the nebula through a narrow square
842 * centered on the central source, this gives r^0 dependence */
843 geometry.iEmissPower = 0;
844 }
845 else if( p.nMatch("SIZE") )
846 {
847 /* set the aperture size: slit width or suface area of the pencil beam */
848 /* units are arcsec for slit width, or arcsec^2 for pencil beam */
849 geometry.size = realnum(p.FFmtRead());
850 if( p.lgEOL() )
851 {
852 p.NoNumb("aperture size");
853 }
854 if( geometry.size <= 0.f )
855 {
856 fprintf( ioQQQ, " The aperture size must be positive. Sorry.\n" );
858 }
859 geometry.lgSizeSet = true;
860 }
861 else if( p.nMatch("COVE") )
862 {
863 /* set the aperture covering factor, see Hazy 1 for a discussion
864 * this is a dimensionless fraction between 0 and 1 */
865 geometry.covaper = realnum(p.FFmtRead());
866 if( p.lgEOL() )
867 {
868 p.NoNumb("aperture covering factor");
869 }
870 if( geometry.covaper <= 0.f || geometry.covaper > 1.f )
871 {
872 fprintf( ioQQQ, " The aperture covering factor must be > 0 and <= 1. Sorry.\n" );
874 }
875 }
876 else
877 {
878 fprintf( ioQQQ, " One of the keywords SLIT, BEAM, SIZE or COVEring factor must appear.\n" );
879 fprintf( ioQQQ, " Sorry.\n" );
881 }
882}
884{
885 DEBUG_ENTRY("ParseAtom()");
886 /* accept both forms of feii */
887 if( p.nMatch("FEII") || p.nMatch("FE II") )
888 {
889 /* parse the atom FeII command - routine in atom_feii.c */
890 ParseAtomFeII(p);
891 }
892
893 else if( p.nMatch("H-LI") )
894 {
895 /* parse the atom h-like command */
897 }
898
899 else if( p.nMatch("HE-L") )
900 {
901 /* parse the atom he-like command */
903 }
904
905 else if( p.nMatch("ROTO") || p.nMatch(" CO ") )
906 {
907 /* ATOM ROTOR same as ATOM CO */
908 fprintf(ioQQQ," The old CO models no longer exist, and this command is no longer supported.\n" );
909 fprintf( ioQQQ, " Sorry.\n" );
911 }
912
913 else if( p.nMatch(" H2 ") )
914 {
915 ParseAtomH2(p);
916 }
917
918 /* enable models from CHIANTI database */
919 else if (p.nMatch("CHIA"))
920 {
921 char chString_quotes_lowercase[INPUT_LINE_LENGTH];
922 bool lgQuotesFound = true;
923 if (p.GetQuote(chString_quotes_lowercase, false))
924 lgQuotesFound = false;
925
926 // option to specify different CloudyChianti.ini file, was initialized
927 // with string CloudyChianti.ini
928 if (lgQuotesFound == true)
929 strcpy(atmdat.chCloudyChiantiFile, chString_quotes_lowercase);
930
931 atmdat.lgChiantiOn = true;
932 if (p.nMatch(" OFF"))
933 {
934 atmdat.lgChiantiOn = false;
935 atmdat.lgChiantiHybrid = false;
936 }
937
938 // hybrid, chianti with OP for higher energy lines
939 // Turn off hybrid, use Chianti only
940 if (p.nMatch(" NO HYBR"))
941 atmdat.lgChiantiHybrid = false;
942
943 // Print which species are being used in output and # of levels
944 if (p.nMatch(" PR"))
945 atmdat.lgChiantiPrint = true;
946
947 // Use Experimental energies exclusively. Default use experimental.
948 if (p.nMatch(" THEO"))
949 atmdat.lgChiantiExp = false;
950
951 // Input the maximum number of Chianti levels to use
952 if (p.nMatch(" LEV"))
953 {
954 if (p.nMatch(" MAX"))
955 {
956 atmdat.nChiantiMaxLevels = LONG_MAX;
957 atmdat.nChiantiMaxLevelsFe = LONG_MAX;
958 }
959 else
960 {
961 atmdat.nChiantiMaxLevelsFe = (long)p.FFmtRead();
962 atmdat.nChiantiMaxLevels = (long)p.FFmtRead();
963 if (p.lgEOL())
964 {
965 p.NoNumb("two numbers, the maximum number of levels in Fe, and in other elements,");
966 }
967 if( atmdat.nChiantiMaxLevelsFe < 2 || atmdat.nChiantiMaxLevels < 2 )
968 {
969 fprintf(ioQQQ,
970 " \nPROBLEM The maximum number of chianti levels should be two or greater.\n");
971 fprintf(ioQQQ, " To turn off the Chianti data use \"Set Chianti off\" instead.\n");
972 fprintf(ioQQQ, " See Hazy 1 for details.\n");
974 }
975 }
976 atmdat.lgChiantiLevelsSet = true;
977 }
978 }
979
980 /* enable models from STOUT database */
981 else if (p.nMatch("STOUT"))
982 {
983 char chString_quotes_lowercase[INPUT_LINE_LENGTH];
984 bool lgQuotesFound = true;
985 if (p.GetQuote(chString_quotes_lowercase, false))
986 lgQuotesFound = false;
987
988 // option to specify different Stout.ini file, was initialized
989 // with string Stout.ini
990 if (lgQuotesFound == true)
991 strcpy(atmdat.chStoutFile, chString_quotes_lowercase);
992
993 atmdat.lgStoutOn = true;
994 // Print which species are being used in output and # of levels
995 if (p.nMatch(" PR"))
996 atmdat.lgStoutPrint = true;
997
998 if (p.nMatch(" OFF"))
999 {
1000 atmdat.lgStoutOn = false;
1001 atmdat.lgStoutHybrid = false;
1002 }
1003
1004 // hybrid, Stout with OP for higher energy lines
1005 // Turn off hybrid, use Stout only
1006 if (p.nMatch(" NO HYBR"))
1007 atmdat.lgStoutHybrid = false;
1008
1009 // Input the maximum number of Stout levels to use
1010 if (p.nMatch(" LEV"))
1011 {
1012 if (p.nMatch(" MAX"))
1013 {
1014 atmdat.nStoutMaxLevels = LONG_MAX;
1015 }
1016 else
1017 {
1018 atmdat.nStoutMaxLevels = (long)p.FFmtRead();
1019 if (p.lgEOL())
1020 {
1021 p.NoNumb("the maximum number of Stout levels,");
1022 }
1023 if( atmdat.nStoutMaxLevels < 2 )
1024 {
1025 fprintf(ioQQQ,
1026 " \nPROBLEM The maximum number of Stout levels should be two or greater.\n");
1027 fprintf(ioQQQ, " To turn off the Stout data use \"Set Stout off\" instead.\n");
1028 fprintf(ioQQQ, " See Hazy 1 for details.\n");
1030 }
1031 }
1032 }
1033
1034 }
1035
1036 /* enable models from LAMDA (Leiden Atomic and Molecular Database) */
1037 else if (p.nMatch("LAMD"))
1038 {
1039 if (p.nMatch(" OFF"))
1040 atmdat.lgLamdaOn = false;
1041 else if (p.nMatch(" ON "))
1042 atmdat.lgLamdaOn = true;
1043 else
1044 {
1045 /* should not have happened ... */
1046 fprintf(ioQQQ, " There should have been an option on this SET LAMDA command.\n");
1047 fprintf(ioQQQ, " consult Hazy to find valid options.\n Sorry.\n");
1049 }
1050 }
1051 else
1052 {
1053 fprintf( ioQQQ, " I could not recognize a keyword on this atom command.\n");
1054 fprintf( ioQQQ, " The available keys are FeII, H-Like, He-like, rotor and H2.\n");
1055 fprintf( ioQQQ, " Sorry.\n" );
1057 }
1058}
1060{
1061 DEBUG_ENTRY("ParseBremsstrahlung()");
1062 /* bremsstrahlung continuum from central object */
1063 strcpy( rfield.chSpType[rfield.nShape], "BREMS" );
1064 rfield.slope[rfield.nShape] =
1065 (realnum)p.FFmtRead();
1066 if( p.lgEOL() )
1067 {
1068 p.NoNumb("temperature");
1069 }
1070
1071 /* temperature is interpreted as log if <= 10, or if keyword is present */
1072 if( rfield.slope[rfield.nShape] <= 10. || p.nMatch(" LOG") )
1073 {
1074 rfield.slope[rfield.nShape] =
1075 pow(10.,rfield.slope[rfield.nShape]);
1076 }
1077 rfield.cutoff[rfield.nShape][0] = 0.;
1078
1079 /* option for vary keyword */
1080 if( optimize.lgVarOn )
1081 {
1082 /* only one parameter */
1083 optimize.nvarxt[optimize.nparm] = 1;
1084 strcpy( optimize.chVarFmt[optimize.nparm], "BREMS, T=%f LOG" );
1085 /* pointer to where to write */
1086 optimize.nvfpnt[optimize.nparm] = input.nRead;
1087 /* log of temp will be pointer */
1088 optimize.vparm[0][optimize.nparm] = (realnum)log10(rfield.slope[rfield.nShape]);
1089 optimize.vincr[optimize.nparm] = 0.5;
1090 ++optimize.nparm;
1091 }
1092 ++rfield.nShape;
1093 if( rfield.nShape >= LIMSPC )
1094 {
1095 /* too many continua were entered */
1096 fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" );
1098 }
1099}
1101{
1102 /* add "extra" cooling, power law temp dependence */
1103 thermal.lgCExtraOn = true;
1104 thermal.CoolExtra = (realnum)pow(10.,p.FFmtRead());
1105 if( p.lgEOL() )
1106 {
1107 p.NoNumb("extra cooling");
1108 }
1109 thermal.cextpw = (realnum)p.FFmtRead();
1110}
1112{
1113 cosmology.redshift_start = (realnum)p.FFmtRead();
1114 cosmology.redshift_current = cosmology.redshift_start;
1115
1116 /* >>chng 06 mar 22, add time option to vary only some continua with time */
1117 if( p.nMatch( "TIME" ) )
1118 rfield.lgTimeVary[(p.m_nqh)] = true;
1119
1120 ParseCMB(cosmology.redshift_start,&(p.m_nqh));
1121
1122 /* option to also set the hydrogen density at given redshift. */
1123 if( p.nMatch("DENS") && !p.lgEOL() )
1124 {
1125 char chStuff[INPUT_LINE_LENGTH];
1126
1127 /* hydrogen density */
1128 sprintf( chStuff , "HDEN %.2e LINEAR", GetDensity( cosmology.redshift_start ) );
1129 p.setline(chStuff);
1130 p.set_point(4);
1131 ParseHDEN(p);
1132 }
1133}
1135{
1136 DEBUG_ENTRY("ParseCosm()");
1137 fprintf(ioQQQ,
1138 "This command is now ambiguous -- please specify either COSMIC RAYS or COSMOLOGY.\nSorry.\n");
1140}
1142{
1143 DEBUG_ENTRY("ParseCovering()");
1144 /* covering factor for gas */
1145 /* The geometric covering factor accounts for how much of 4\pi is
1146 * covered by gas, and so linearly multiplies the predicted intensities */
1147 geometry.covgeo = (realnum)p.FFmtRead();
1148
1149 if( p.lgEOL() )
1150 {
1151 p.NoNumb("covering factor");
1152 }
1153
1154 /* if negative then log, convert to linear */
1155 if( geometry.covgeo <= 0. )
1156 {
1157 geometry.covgeo = (realnum)pow((realnum)10.f,geometry.covgeo);
1158 }
1159
1160 if( geometry.covgeo > 1. )
1161 {
1162 fprintf( ioQQQ, " A covering factor greater than 1 makes no physical sense. Sorry.\n" );
1164 }
1165
1166 /* radiative transfer covering factor will be equal to the geometric one */
1167 geometry.covrt = geometry.covgeo;
1168}
1170{
1171 /* height of cylinder in cm */
1172 radius.lgCylnOn = true;
1173 radius.CylindHigh = pow(10.,p.FFmtRead());
1174 if( p.lgEOL() )
1175 {
1176 p.NoNumb("height");
1177 }
1178}
1180{
1181 DEBUG_ENTRY( "ParseDarkMatter()" );
1182
1183 if( p.nMatch(" NFW") )
1184 {
1185 /* specify an NFW profile */
1186 /* >>refer dark matter Navarro, Frenk, & White, 1996, ApJ, 462, 563 */
1187
1188 dark.r_200 = (realnum)p.getNumberCheckAlwaysLog("NFW r_200");
1189 /* Let r_s have a default corresponding to c_200 = 10. */
1190 dark.r_s = (realnum)p.getNumberDefaultAlwaysLog("NFW r_s", log10(dark.r_200)-1. );
1191 dark.lgNFW_Set = true;
1192
1193 /* option for vary keyword */
1194 if( optimize.lgVarOn )
1195 {
1196 /* only one parameter */
1197 optimize.nvarxt[optimize.nparm] = 1;
1198 strcpy( optimize.chVarFmt[optimize.nparm], "DARK NFW %f" );
1199 /* pointer to where to write */
1200 optimize.nvfpnt[optimize.nparm] = input.nRead;
1201 /* log of temp will be pointer */
1202 optimize.vparm[0][optimize.nparm] = (realnum)log10(dark.r_200);
1203 optimize.vincr[optimize.nparm] = 0.5;
1204 ++optimize.nparm;
1205 }
1206 }
1207 else
1208 {
1209 fprintf( ioQQQ, " Did not recognize a valid option for this DARK command.\nSorry.\n\n" );
1211 }
1212
1213 return;
1214}
1216{
1217 DEBUG_ENTRY("ParseDielectronic()");
1218 /* >>chng 05 dec 21, change dielectronic command to set dielectronic recombination command */
1219 fprintf( ioQQQ, " The DIELectronic command has been replaced with the SET DIELectronic recombination command.\n" );
1220 fprintf( ioQQQ, " Please have a look at Hazy.\n Sorry.\n\n" );
1222}
1224{
1225 DEBUG_ENTRY("ParseDiffuse()");
1226 /* set method of transferring diffuse fields,
1227 * default is outward only, cdDffTrns label "OU2", set in zero.c */
1228 if( p.nMatch(" OTS") )
1229 {
1230 if( p.nMatch("SIMP") )
1231 {
1232 /* this is simple ots, which never allows photons to escape */
1233 strcpy( rfield.chDffTrns, "OSS" );
1234 }
1235 else
1236 {
1237 /* this is regular ots, which allows photons to escape */
1238 strcpy( rfield.chDffTrns, "OTS" );
1239 }
1240 rfield.lgOutOnly = false;
1241 }
1242 else if( p.nMatch(" OUT") )
1243 {
1244 rfield.lgOutOnly = true;
1245 long int j = (long int)p.FFmtRead();
1246 if( p.lgEOL() )
1247 {
1248 /* this is the default set in zero */
1249 strcpy( rfield.chDffTrns, "OU2" );
1250 }
1251 else
1252 {
1253 if( j > 0 && j < 10 )
1254 {
1255 sprintf( rfield.chDffTrns, "OU%1ld", j );
1256 }
1257 else
1258 {
1259 fprintf( ioQQQ, " must be between 1 and 9 \n" );
1261 }
1262 }
1263 }
1264
1265 else
1266 {
1267 fprintf( ioQQQ, " There should have been OUTward or OTS on this line. Sorry.\n" );
1269 }
1270}
1272{
1273 /* distance to the object */
1274 radius.distance = p.FFmtRead();
1275 if( p.lgEOL() )
1276 {
1277 p.NoNumb("distance");
1278 }
1279 /* default is for quantity to be log of distance, linear keyword
1280 * overrides this - is LINEAR is not on line then exp */
1281 if( !p.nMatch("LINE" ) )
1282 {
1283 radius.distance = pow(10., radius.distance );
1284 }
1285 /* default is radius in cm - if parsecs appears then must
1286 * convert to cm */
1287 if( p.nMatch("PARS") )
1288 {
1289 radius.distance *= PARSEC;
1290 }
1291
1292 /* vary option */
1293 if( optimize.lgVarOn )
1294 {
1295 strcpy( optimize.chVarFmt[optimize.nparm], "DISTANCE %f LOG" );
1296 optimize.nvfpnt[optimize.nparm] = input.nRead;
1297 optimize.vparm[0][optimize.nparm] = realnum(log10(radius.distance));
1298 optimize.vincr[optimize.nparm] = 0.3f;
1299 optimize.nvarxt[optimize.nparm] = 1;
1300 ++optimize.nparm;
1301 }
1302}
1304{
1305 rt.DoubleTau = 2.;
1306}
1308{
1309 dense.EdenExtra = (realnum)pow(10.,p.FFmtRead());
1310 if( p.lgEOL() )
1311 {
1312 p.NoNumb("electron density");
1313 }
1314 /* warn that this model is meaningless */
1315 phycon.lgPhysOK = false;
1316}
1318{
1319 DEBUG_ENTRY("ParseEnergy()");
1320 /* energy density (degrees K) of this continuum source */
1321 if( (p.m_nqh) >= LIMSPC )
1322 {
1323 /* too many continua were entered */
1324 fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" );
1326 }
1327 /* silly, but calms down the lint */
1328 ASSERT( (p.m_nqh) < LIMSPC );
1329 strcpy( rfield.chRSpec[(p.m_nqh)], "SQCM" );
1330 realnum teset = (realnum)p.FFmtRead();
1331 if( p.lgEOL() )
1332 {
1333 p.NoNumb("energy density");
1334 }
1335 /* set radius to very large value if not already set */
1336 /* >>chng 01 jul 24, from Radius == 0 to this, as per PvH comments */
1337 if( !radius.lgRadiusKnown )
1338 {
1339 /* set radius to large value */
1340 radius.Radius = pow(10.,radius.rdfalt);
1341 }
1342
1343 /* convert temp to log, recognize linear option */
1344 if( !p.nMatch(" LOG") && (p.nMatch("LINE") || teset > 10.) )
1345 {
1346 /* option for linear temperature, must store log */
1347 teset = (realnum)log10(teset);
1348 }
1349
1350 if( teset > 5. )
1351 {
1352 fprintf( ioQQQ, " This intensity may be too large. The code may crash due to overflow. Was log intended?\n" );
1353 }
1354
1355 /* teset is now log of temp, now get log of total luminosity */
1356 strcpy( rfield.chSpNorm[(p.m_nqh)], "LUMI" );
1357
1358 /* full range of continuum will be used */
1359 rfield.range[(p.m_nqh)][0] = rfield.emm;
1360 rfield.range[(p.m_nqh)][1] = rfield.egamry;
1361 rfield.totpow[(p.m_nqh)] = (4.*(teset) - 4.2464476 + 0.60206);
1362
1363 /* >>chng 06 mar 22, add time option to vary only some continua with time */
1364 if( p.nMatch( "TIME" ) )
1365 rfield.lgTimeVary[(p.m_nqh)] = true;
1366
1367 /* vary option */
1368 if( optimize.lgVarOn )
1369 {
1370 strcpy( optimize.chVarFmt[optimize.nparm], "ENERGY DENSITY %f LOG" );
1371 if( rfield.lgTimeVary[(p.m_nqh)] )
1372 strcat( optimize.chVarFmt[optimize.nparm], " TIME" );
1373 /* pointer to where to write */
1374 optimize.nvfpnt[optimize.nparm] = input.nRead;
1375 optimize.vparm[0][optimize.nparm] = teset;
1376 optimize.vincr[optimize.nparm] = 0.1f;
1377 optimize.nvarxt[optimize.nparm] = 1;
1378 ++optimize.nparm;
1379 }
1380
1381 /* finally increment number of continua */
1382 ++(p.m_nqh);
1383}
1385{
1386 /* save previous value */
1387 long int j = conv.LimFail;
1388 conv.LimFail = (long int)p.FFmtRead();
1389 if( p.lgEOL() )
1390 {
1391 p.NoNumb("limit");
1392 }
1393
1394 /* >>chng 01 mar 14, switch logic on maps */
1395 /* ' map' flag, make map when failure, default is no map,
1396 * second check is so that no map does not trigger a map */
1397 if( p.nMatch(" MAP") && !p.nMatch(" NO ") )
1398 {
1399 conv.lgMap = true;
1400 }
1401
1402 /* complain if failures was increased above default */
1403 if( conv.LimFail > j )
1404 {
1405 fprintf( ioQQQ, " This command should not be necessary.\n" );
1406 fprintf( ioQQQ, " Please show this input stream to Gary Ferland if this command is really needed for this simulation.\n" );
1407 }
1408}
1410{
1411 /* filling factor, power law exponent (default=1., 0.) */
1412 realnum a = (realnum)p.FFmtRead();
1413 if( p.lgEOL() )
1414 {
1415 p.NoNumb("filling factor");
1416 }
1417
1418 if( a <= 0. || p.nMatch(" LOG") )
1419 {
1420 /* number less than or equal to 0, must have been entered as log */
1421 geometry.FillFac = (realnum)pow((realnum)10.f,a);
1422 }
1423 else
1424 {
1425 /* number greater than zero, must have been the real thing */
1426 geometry.FillFac = a;
1427 }
1428 if( geometry.FillFac > 1.0 )
1429 {
1430 if( called.lgTalk )
1431 fprintf( ioQQQ, " Filling factor > 1, reset to 1\n" );
1432 geometry.FillFac = 1.;
1433 }
1434 geometry.fiscal = geometry.FillFac;
1435
1436 /* option to have power law dependence, filpow set to 0 in zerologic */
1437 geometry.filpow = (realnum)p.FFmtRead();
1438
1439 /* vary option */
1440 if( optimize.lgVarOn )
1441 {
1442 strcpy( optimize.chVarFmt[optimize.nparm], "FILLING FACTOR= %f LOG power= %f" );
1443
1444 /* pointer to where to write */
1445 optimize.nvfpnt[optimize.nparm] = input.nRead;
1446 optimize.vparm[0][optimize.nparm] = (realnum)log10(geometry.FillFac);
1447 optimize.vincr[optimize.nparm] = 0.5f;
1448
1449 /* power law dependence here, but cannot be varied */
1450 optimize.vparm[1][optimize.nparm] = geometry.filpow;
1451 optimize.nvarxt[optimize.nparm] = 2;
1452
1453 /* do not allow filling factor to go positive */
1454 optimize.varang[optimize.nparm][0] = -FLT_MAX;
1455 optimize.varang[optimize.nparm][1] = 0.f;
1456 ++optimize.nparm;
1457 }
1458}
1460{
1461 bool lgNu2 = false;
1462 /* in reads2 */
1463 ParseF_nu(p,"SQCM",lgNu2);
1464}
1466{
1467 thermal.ConstTemp = (realnum)p.FFmtRead();
1468 if( p.lgEOL() )
1469 {
1470 p.NoNumb("temperature");
1471 }
1472
1473 if( p.nMatch(" LOG") || (thermal.ConstTemp <= 10. &&
1474 !p.nMatch("LINE")) )
1475 {
1476 thermal.ConstTemp = (realnum)pow((realnum)10.f,thermal.ConstTemp);
1477 }
1478
1479 /* low energy bounds of continuum array do not permit too-low a temp */
1480 if( thermal.ConstTemp < 3. )
1481 {
1482 fprintf( ioQQQ, " TE reset to 3K: entered number too small.\n" );
1483 thermal.ConstTemp = 3.;
1484 }
1485}
1487{
1488 /* enter a fudge factor */
1489 fudgec.nfudge = 0;
1490 for( long int j=0; j < NFUDGC; j++ )
1491 {
1492 fudgec.fudgea[j] = (realnum)p.FFmtRead();
1493 /* fudgec.nfudge is number of factors on the line */
1494 if( !p.lgEOL() )
1495 fudgec.nfudge = j+1;
1496 }
1497 if( fudgec.nfudge == 0 )
1498 p.NoNumb("fudge factor");
1499
1500 /* vary option */
1501 if( optimize.lgVarOn )
1502 {
1503 /* no luminosity options on vary */
1504 optimize.nvarxt[optimize.nparm] = 1;
1505 strcpy( optimize.chVarFmt[optimize.nparm], "FUDGE= %f" );
1506 /* enter remaining parameters as constants */
1507 char chHold[1000];
1508 for( long int j=1; j<fudgec.nfudge; ++j )
1509 {
1510 sprintf( chHold , " %f" , fudgec.fudgea[j] );
1511 strcat( optimize.chVarFmt[optimize.nparm] , chHold );
1512 }
1513 optimize.lgOptimizeAsLinear[optimize.nparm] = true;
1514
1515 /* pointer to where to write */
1516 optimize.nvfpnt[optimize.nparm] = input.nRead;
1517 /* fudge factor stored here */
1518 optimize.vparm[0][optimize.nparm] = fudgec.fudgea[0];
1519 /* the increment in the first steps away from the original value */
1520 optimize.vincr[optimize.nparm] = abs(0.2f*fudgec.fudgea[0]);
1521 /* we have no clue what to use when initial estimate is 0... */
1522 if( optimize.vincr[optimize.nparm] == 0.f )
1523 optimize.vincr[optimize.nparm] = 1.f;
1524 ++optimize.nparm;
1525 }
1526}
1528{
1529 DEBUG_ENTRY("ParsePGrains()");
1530 fprintf(ioQQQ," Sorry, this command is obsolete, you can now use the normal GRAINS command.\n");
1532}
1534{
1535 DEBUG_ENTRY("ParseGravity()");
1536
1537 if( p.nMatch("EXTE") )
1538 {
1539 /* external mass: M/Msun or Sigma/(Msun/pc^2) depending on symmetry */
1540 /* if no number is read, FFmtRead returns 0 */
1541 /* default is linear, unless "log" keyword is present */
1542 double M_i = p.FFmtRead();
1543
1544 if( !p.lgEOL() && p.nMatch("LOG") )
1545 M_i = pow( 10., M_i );
1546 pressure.external_mass[0].push_back( M_i );
1547
1548 /* extent of the mass distribution, in pc */
1549 /* default is linear, unless "log" keyword is present */
1550 double x_i = p.FFmtRead();
1551
1552 if( !p.lgEOL() && p.nMatch("LOG") )
1553 x_i = pow( 10., x_i );
1554 pressure.external_mass[1].push_back( x_i * PARSEC );
1555
1556 /* exponential index of the mass distribution */
1557 pressure.external_mass[2].push_back( p.FFmtRead() );
1558 }
1559 else
1560 {
1561 /* choose spherical or mid-plane symmetry for the gas mass distribution */
1562 if( p.nMatch("SPHE") )
1563 {
1564 pressure.gravity_symmetry = 0;
1565 }
1566 else if( p.nMatch("PLAN") )
1567 {
1568 pressure.gravity_symmetry = 1;
1569 }
1570 else
1571 {
1572 fprintf( ioQQQ, " The symmetry of the gravitational mass must be specified explicitly. Sorry.\n" );
1574 }
1575
1576 /* external mass, proportional to the gas density */
1577 /* default is linear, unless "log" keyword is present */
1578 pressure.self_mass_factor = p.FFmtRead();
1579
1580 if( p.lgEOL() )
1581 pressure.self_mass_factor = 1.;
1582 else if( p.nMatch("LOG") )
1583 {
1584 pressure.self_mass_factor = pow( 10., pressure.self_mass_factor );
1585 }
1586 }
1587}
1589{
1590 DEBUG_ENTRY("ParseHeLike()");
1591 fprintf(ioQQQ,"Sorry, this command is replaced with ATOM HE-LIKE\n");
1593}
1595{
1596 DEBUG_ENTRY("ParseHelp()");
1597 p.help(ioQQQ);
1598}
1600{
1601 DEBUG_ENTRY(" ParseHExtra()");
1602 hextra.TurbHeat = (realnum)pow(10.,p.FFmtRead());
1603 if( p.lgEOL() )
1604 p.NoNumb("extra heating first parameter" );
1605
1606 /* save initial value in case reset with time dependent command */
1607 hextra.TurbHeatSave = hextra.TurbHeat;
1608
1609 /* remember which type of scale dependence we will use */
1610 const char *chHextraScale;
1611 /* these are optional numbers on depth or density option */
1612 realnum par1=0. , par2=0.;
1613 long int nPar;
1614 if( p.nMatch("DEPT") )
1615 {
1616 /* depend on depth */
1617 hextra.lgHextraDepth = true;
1618 chHextraScale = "DEPTH";
1619 /* optional scale radius */
1620 realnum a = (realnum)p.FFmtRead();
1621 if( p.lgEOL() )
1622 p.NoNumb("depth");
1623 else
1624 hextra.turrad = (realnum)pow((realnum)10.f,a);
1625
1626 /* depth from shielded face, to mimic illumination from both sides
1627 * may not be specified */
1628 realnum b = (realnum)p.FFmtRead();
1629 if( p.lgEOL() )
1630 {
1631 hextra.turback = 0.;
1632 nPar = 2;
1633 }
1634 else
1635 {
1636 hextra.turback = (realnum)pow((realnum)10.f,b);
1637 nPar = 3;
1638 }
1639 par1 = a;
1640 par2 = b;
1641 }
1642 else if( p.nMatch("DENS") )
1643 {
1644 /* depends on density */
1645 chHextraScale = "DENSITY";
1646 hextra.lgHextraDensity = true;
1647
1648 /* the log of the density */
1649 realnum a = (realnum)p.FFmtRead();
1650 /* if number not entered then unity is used, heating
1651 * proportional to density */
1652 hextra.HextraScaleDensity = (realnum)pow((realnum)10.f,a);
1653 par1 = a;
1654 par2 = 0;
1655 nPar = 2;
1656 }
1657 else if( p.nMatch("SS") )
1658 {
1659 /* alpha disk model, specify alpha, mass of black hole, and distance */
1660 chHextraScale = "SS";
1661 hextra.lgHextraSS = true;
1662
1663 /* the parameter alpha of alpha-disk model */
1664 hextra.HextraSSalpha = hextra.TurbHeat;
1665
1666 /* mass in solar masses of center black hole */
1667 realnum a = (realnum)p.FFmtRead();
1668 if( p.lgEOL() )
1669 p.NoNumb("hextraSS Mass");
1670 hextra.HextraSS_M = a*SOLAR_MASS;
1671
1672 /* radius (cm) from center */
1673 realnum b = (realnum)p.FFmtRead();
1674 if( p.lgEOL() )
1675 p.NoNumb("hextraSS radius");
1676 hextra.HextraSSradius = b;
1677 par1 = a;
1678 par2 = 0;
1679 nPar = 2;
1680 }
1681 else
1682 {
1683 /* constant heating */
1684 chHextraScale = "";
1685 nPar = 1;
1686 }
1687
1688 /* option to have heating vary with time */
1689 if( p.nMatch( "TIME" ) )
1690 hextra.lgTurbHeatVaryTime = true;
1691
1692 /* vary option */
1693 if( optimize.lgVarOn )
1694 {
1695 if( hextra.lgHextraSS )
1696 {
1697 fprintf(ioQQQ,"Sorry, HEXTRA SS command does not now support vary option.\n");
1699 }
1700 /* 1 to 3 options on hextra vary */
1701 optimize.nvarxt[optimize.nparm] = nPar;
1702 optimize.vparm[0][optimize.nparm] = log10(hextra.TurbHeat);
1703 optimize.vparm[1][optimize.nparm] = par1;
1704 optimize.vparm[2][optimize.nparm] = par2;
1705
1706 /* the keyword LOG is not used above, but is checked elsewhere */
1707 strcpy( optimize.chVarFmt[optimize.nparm], "HEXTra %f LOG " );
1708 strcat( optimize.chVarFmt[optimize.nparm], chHextraScale );
1709 while( nPar > 1 )
1710 {
1711 /* add extra spots to write parameters */
1712 --nPar;
1713 strcat( optimize.chVarFmt[optimize.nparm], " %f" );
1714 }
1715 if( hextra.lgTurbHeatVaryTime )
1716 strcat( optimize.chVarFmt[optimize.nparm], " TIME" );
1717 /* pointer to where to write */
1718 optimize.nvfpnt[optimize.nparm] = input.nRead;
1719 /* the increment in the first steps away from the original value */
1720 optimize.vincr[optimize.nparm] = 0.1f;
1721 ++optimize.nparm;
1722 }
1723}
1724
1726{
1727 thermal.lgTeHigh = true;
1728}
1730{
1731 DEBUG_ENTRY("ParseHydrogen()");
1732 fprintf(ioQQQ," Sorry, this command has been replaced with the ATOM H-LIKE command.\n");
1734}
1736{
1737 DEBUG_ENTRY("ParseInitCount()");
1738 /* read cloudy.ini initialization file
1739 * need to read in file every time, since some vars
1740 * get reset in zero - would be unsafe not to read in again */
1741 ParseInit(p);
1742
1743 /* check that only one init file was in the input stream -
1744 * we cannot now read more than one */
1745 ++p.m_nInitFile;
1746 if( p.m_nInitFile > 1 )
1747 {
1748 fprintf( ioQQQ,
1749 " This is the second init file, I can only handle one.\nSorry.\n" );
1751 }
1752
1753 /* we will put the data for the ini file at the end of the array of
1754 * line images, this is the increment for stuffing the lines in - negative */
1755 input.iReadWay = -1;
1756
1757 /* initialize array reader, this sub does nothing but set
1758 * initial value of a variable, depending on value of iReadWay */
1759 input.init();
1760}
1762{
1763 DEBUG_ENTRY("ParseIntensity()");
1764 /* intensity of this continuum source */
1765 if( (p.m_nqh) >= LIMSPC )
1766 {
1767 /* too many continua were entered */
1768 fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" );
1770 }
1771
1772 /* silly, but calms down the lint */
1773 ASSERT( (p.m_nqh) < LIMSPC );
1774 strcpy( rfield.chRSpec[(p.m_nqh)], "SQCM" );
1775 rfield.totpow[(p.m_nqh)] = p.FFmtRead();
1776 if( p.lgEOL() )
1777 p.NoNumb("intensity");
1778
1779 /* set radius to very large value if not already set */
1780 /* >>chng 01 jul 24, from Radius == 0 to this, as per PvH comments */
1781 if( !radius.lgRadiusKnown )
1782 {
1783 /* set radius to large value */
1784 radius.Radius = pow(10.,radius.rdfalt);
1785 }
1786 if( p.nMatch("LINE") )
1787 {
1788 /* silly, but calms down the lint */
1789 ASSERT( (p.m_nqh) < LIMSPC );
1790 /* option for linear input parameter */
1791 rfield.totpow[(p.m_nqh)] = log10(rfield.totpow[(p.m_nqh)]);
1792 }
1793 strcpy( rfield.chSpNorm[(p.m_nqh)], "LUMI" );
1794 /* ParseRangeOption in readsun */
1796
1797 /* >>chng 06 mar 22, add time option to vary only some continua with time */
1798 if( p.nMatch( "TIME" ) )
1799 rfield.lgTimeVary[(p.m_nqh)] = true;
1800
1801 /* vary option */
1802 if( optimize.lgVarOn )
1803 {
1804 /* create the format we will use to write the command */
1805 strcpy( optimize.chVarFmt[optimize.nparm], "INTENSITY %f LOG range %f %f" );
1806 if( rfield.lgTimeVary[(p.m_nqh)] )
1807 strcat( optimize.chVarFmt[optimize.nparm], " TIME" );
1808 /* array index for this command */
1809 optimize.nvfpnt[optimize.nparm] = input.nRead;
1810 optimize.vparm[0][optimize.nparm] = (realnum)rfield.totpow[(p.m_nqh)];
1811 optimize.vincr[optimize.nparm] = 0.5;
1812 /* range option, but cannot be varied */
1813 optimize.vparm[1][optimize.nparm] = (realnum)log10(rfield.range[(p.m_nqh)][0]);
1814 optimize.vparm[2][optimize.nparm] = (realnum)log10(rfield.range[(p.m_nqh)][1]);
1815 optimize.nvarxt[optimize.nparm] = 3;
1816 ++optimize.nparm;
1817 }
1818 ++(p.m_nqh);
1819}
1821{
1822 /* iterate to get optical depths and diffuse field properly */
1823 iterations.itermx = (long int)p.FFmtRead() - 1;
1824 iterations.itermx = MAX2(iterations.itermx,1);
1825 /* >>chng 06 mar 19, malloc itrdim arrays
1826 * iterations.iter_malloc is -1 before space allocated and set to
1827 * itermx if not previously set */
1828 if( iterations.itermx > iterations.iter_malloc - 1 )
1829 {
1830 long int iter_malloc_save = iterations.iter_malloc;
1831 /* add one for mismatch between array dim and iterations defn,
1832 * another few for safety */
1833 iterations.iter_malloc = iterations.itermx+3;
1834 iterations.IterPrnt = (long int*)REALLOC(iterations.IterPrnt ,
1835 (size_t)iterations.iter_malloc*sizeof(long int) );
1836 geometry.nend = (long int*)REALLOC(geometry.nend ,
1837 (size_t)iterations.iter_malloc*sizeof(long int) );
1838 radius.StopThickness = (double*)REALLOC(radius.StopThickness ,
1839 (size_t)iterations.iter_malloc*sizeof(double) );
1840 /* >>chng 06 jun 07, need to set IterPrnt, thickness, and nend */
1841 for(long int j=iter_malloc_save; j<iterations.iter_malloc; ++j )
1842 {
1843 iterations.IterPrnt[j] = iterations.IterPrnt[iter_malloc_save-1];
1844 geometry.nend[j] = geometry.nend[iter_malloc_save-1];
1845 radius.StopThickness[j] = radius.StopThickness[iter_malloc_save-1];
1846 }
1847 }
1848
1849 if( p.nMatch("CONV") )
1850 {
1851 /* option to keep iterating until it converges, checks on convergence
1852 * are done in update, and checked again in prtcomment */
1853 conv.lgAutoIt = true;
1854 /* above would have been legitimate setting of ITERMX, set to default 10 */
1855 if( p.lgEOL() )
1856 {
1857 iterations.itermx = 10 - 1;
1858 }
1859 realnum a = (realnum)p.FFmtRead();
1860 /* change convergence criteria, preset in zero */
1861 if( !p.lgEOL() )
1862 {
1863 conv.autocv = a;
1864 }
1865 }
1866}
1868{
1869 /* this is the specific luminosity at nu
1870 * following says n(nu) not nuF(nu) */
1871 bool lgNu2 = false;
1872 /* in reads2 */
1873 ParseF_nu(p,"4 PI",lgNu2);
1874}
1876{
1877 DEBUG_ENTRY("ParseLaser()");
1878 /* mostly one frequency (a laser) to check gamma's,*/
1879 /* say the continuum type */
1880 strcpy( rfield.chSpType[rfield.nShape], "LASER" );
1881
1882 /* scan off the laser's energy ion Rydbergs */
1883 rfield.slope[rfield.nShape] = p.FFmtRead();
1884
1885 /* negative energies are logs */
1886 if( rfield.slope[rfield.nShape] <= 0. )
1887 {
1888 rfield.slope[rfield.nShape] =
1889 pow(10.,rfield.slope[rfield.nShape]);
1890 }
1891 if( p.lgEOL() )
1892 {
1893 p.NoNumb("energy");
1894 }
1895
1896 /* next number is optional relative width of laser */
1897 rfield.cutoff[rfield.nShape][0] = p.FFmtRead();
1898 if( p.lgEOL() )
1899 {
1900 /* default width is 0.05 */
1901 rfield.cutoff[rfield.nShape][0] = 0.05;
1902 }
1903
1904 /* increment counter making sure we still have space enough */
1905 ++rfield.nShape;
1906 if( rfield.nShape >= LIMSPC )
1907 {
1908 /* too many continua were entered */
1909 fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" );
1911 }
1912}
1914{
1915 DEBUG_ENTRY("ParseLuminosity()");
1916 /* luminosity of this continuum source */
1917 if( (p.m_nqh) >= LIMSPC )
1918 {
1919 /* too many continua were entered */
1920 fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" );
1922 }
1923
1924 strcpy( rfield.chRSpec[(p.m_nqh)], "4 PI" );
1925 rfield.totpow[(p.m_nqh)] = p.FFmtRead();
1926 if( p.lgEOL() )
1927 {
1928 p.NoNumb("luminosity");
1929 }
1930 if( p.nMatch("LINE") )
1931 {
1932 /* option for linear input parameter */
1933 rfield.totpow[(p.m_nqh)] = log10(rfield.totpow[(p.m_nqh)]);
1934 }
1935
1936 strcpy( rfield.chSpNorm[(p.m_nqh)], "LUMI" );
1937
1938 /* the solar option - number is log total solar luminosity */
1939 if( p.nMatch("SOLA") )
1940 {
1941 /* option to use log of total luminosity in solar units */
1942 rfield.range[(p.m_nqh)][0] = rfield.emm;
1943 rfield.range[(p.m_nqh)][1] = rfield.egamry;
1944 /* log of solar luminosity in erg s-1 */
1945 rfield.totpow[(p.m_nqh)] += 33.5827f;
1946 }
1947 else
1948 {
1949 /* check if range is present and parse it if it is - ParseRangeOption in readsun
1950 * this includes TOTAL keyword for total luminosity */
1952 }
1953
1954 /* >>chng 06 mar 22, add time option to vary only some continua with time */
1955 if( p.nMatch( "TIME" ) )
1956 rfield.lgTimeVary[(p.m_nqh)] = true;
1957
1958 /* vary option */
1959 if( optimize.lgVarOn )
1960 {
1961 strcpy( optimize.chVarFmt[optimize.nparm], "LUMINOSITY %f LOG range %f %f" );
1962 if( rfield.lgTimeVary[(p.m_nqh)] )
1963 strcat( optimize.chVarFmt[optimize.nparm], " TIME" );
1964 /* pointer to where to write */
1965 optimize.nvfpnt[optimize.nparm] = input.nRead;
1966 optimize.vparm[0][optimize.nparm] = (realnum)rfield.totpow[(p.m_nqh)];
1967 optimize.vincr[optimize.nparm] = 0.5;
1968 /* range option in, but cannot be varied */
1969 optimize.vparm[1][optimize.nparm] = (realnum)log10(rfield.range[(p.m_nqh)][0]);
1970 optimize.vparm[2][optimize.nparm] = (realnum)log10(rfield.range[(p.m_nqh)][1]);
1971 optimize.nvarxt[optimize.nparm] = 3;
1972 ++optimize.nparm;
1973 }
1974 ++(p.m_nqh);
1975}
1977{
1978 /* heating and ionization due to fast neutrons */
1979 hextra.lgNeutrnHeatOn = true;
1980
1981 /* first number is neutron luminosity
1982 * relative to bolometric luminosity */
1983 hextra.frcneu = (realnum)p.FFmtRead();
1984 if( p.lgEOL() )
1985 p.NoNumb("neutron luminosity");
1986
1987 /* save as log of fraction */
1988 if( hextra.frcneu > 0. )
1989 hextra.frcneu = (realnum)log10(hextra.frcneu);
1990
1991 /* second number is efficiency */
1992 hextra.effneu = (realnum)p.FFmtRead();
1993 if( p.lgEOL() )
1994 {
1995 hextra.effneu = 1.0;
1996 }
1997 else
1998 {
1999 if( hextra.effneu <= 0. )
2000 hextra.effneu = (realnum)pow((realnum)10.f,hextra.effneu);
2001 }
2002}
2004{
2005 /* flux density of this continuum source, at optional frequency
2006 * expressed as product nu*f_nu */
2007 bool lgNu2 = true;
2008 /* in reads2 */
2009 ParseF_nu(p,"SQCM",lgNu2);
2010}
2012{
2013 /* specific luminosity density of this continuum source, at opt freq
2014 * expressed as product nu*f_nu */
2015 bool lgNu2 = true;
2016 /* in reads2 */
2017 ParseF_nu(p,"4 PI",lgNu2);
2018}
2020{
2021 DEBUG_ENTRY("ParsePhi()");
2022 /* enter phi(h), the number of h-ionizing photons/cm2 */
2023 if( (p.m_nqh) >= LIMSPC )
2024 {
2025 /* too many continua were entered */
2026 fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" );
2028 }
2029 /* silly, but calms down the lint */
2030 ASSERT( (p.m_nqh) < LIMSPC );
2031 strcpy( rfield.chRSpec[(p.m_nqh)], "SQCM" );
2032 strcpy( rfield.chSpNorm[(p.m_nqh)], "PHI " );
2033 rfield.totpow[(p.m_nqh)] = p.FFmtRead();
2034 if( p.lgEOL() )
2035 {
2036 p.NoNumb("number of h-ionizing photons");
2037 }
2038 /* set radius to very large value if not already set */
2039 /* >>chng 01 jul 24, from Radius == 0 to this, as per PvH comments */
2040 if( !radius.lgRadiusKnown )
2041 {
2042 /* set radius to large value */
2043 radius.Radius = pow(10.,radius.rdfalt);
2044 }
2045 /* check if continuum so intense that crash is likely */
2046 if( rfield.totpow[(p.m_nqh)] > 29. )
2047 {
2048 fprintf( ioQQQ, " Is the flux for this continuum correct?\n" );
2049 fprintf( ioQQQ, " It appears too bright to me.\n" );
2050 }
2051 /* ParseRangeOption in readsun xx parse_rangeoption*/
2053
2054 /* >>chng 06 mar 22, add time option to vary only some continua with time */
2055 if( p.nMatch( "TIME" ) )
2056 rfield.lgTimeVary[(p.m_nqh)] = true;
2057
2058 /* vary option */
2059 if( optimize.lgVarOn )
2060 {
2061 strcpy( optimize.chVarFmt[optimize.nparm], "phi(h) %f LOG range %f %f" );
2062 if( rfield.lgTimeVary[(p.m_nqh)] )
2063 strcat( optimize.chVarFmt[optimize.nparm], " TIME" );
2064 /* pointer to where to write */
2065 optimize.nvfpnt[optimize.nparm] = input.nRead;
2066 optimize.vparm[0][optimize.nparm] = (realnum)rfield.totpow[(p.m_nqh)];
2067 optimize.vincr[optimize.nparm] = 0.5;
2068 /* range option in, but cannot be varied */
2069 optimize.vparm[1][optimize.nparm] = (realnum)log10(rfield.range[(p.m_nqh)][0]);
2070 optimize.vparm[2][optimize.nparm] = (realnum)log10(rfield.range[(p.m_nqh)][1]);
2071 optimize.nvarxt[optimize.nparm] = 3;
2072 ++optimize.nparm;
2073 }
2074 ++(p.m_nqh);
2075}
2077{
2078 DEBUG_ENTRY("ParseQH()");
2079 /* log of number of ionizing photons */
2080 if( (p.m_nqh) >= LIMSPC )
2081 {
2082 /* too many continua were entered */
2083 fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" );
2085 }
2086
2087 /* silly, but calms down the lint */
2088 ASSERT( (p.m_nqh) < LIMSPC );
2089 strcpy( rfield.chRSpec[(p.m_nqh)], "4 PI" );
2090 strcpy( rfield.chSpNorm[(p.m_nqh)], "Q(H)" );
2091 rfield.totpow[(p.m_nqh)] = p.FFmtRead();
2092 if( rfield.totpow[(p.m_nqh)] > 100. && called.lgTalk )
2093 {
2094 fprintf( ioQQQ, " Is this reasonable?\n" );
2095 }
2096 if( p.lgEOL() )
2097 {
2098 p.NoNumb("number of ionizing photons");
2099 }
2100 /* ParseRangeOption in readsun */
2102
2103 /* >>chng 06 mar 22, add time option to vary only some continua with time */
2104 if( p.nMatch( "TIME" ) )
2105 rfield.lgTimeVary[(p.m_nqh)] = true;
2106
2107 /* vary option */
2108 if( optimize.lgVarOn )
2109 {
2110 strcpy( optimize.chVarFmt[optimize.nparm], "Q(H) %f LOG range %f %f" );
2111 if( rfield.lgTimeVary[(p.m_nqh)] )
2112 strcat( optimize.chVarFmt[optimize.nparm], " TIME" );
2113 /* pointer to where to write */
2114 optimize.nvfpnt[optimize.nparm] = input.nRead;
2115 optimize.vparm[0][optimize.nparm] = (realnum)rfield.totpow[(p.m_nqh)];
2116 optimize.vincr[optimize.nparm] = 0.5;
2117 /* range option in, but cannot be varied */
2118 optimize.vparm[1][optimize.nparm] = (realnum)log10(rfield.range[(p.m_nqh)][0]);
2119 optimize.vparm[2][optimize.nparm] = (realnum)log10(rfield.range[(p.m_nqh)][1]);
2120 optimize.nvarxt[optimize.nparm] = 3;
2121 ++optimize.nparm;
2122 }
2123 /* increment number of luminosity sources specified */
2124 ++(p.m_nqh);
2125}
2127{
2128 /* this is the Roberto Terlivich command, no telling if it still works */
2129 radius.dRadSign = -1.;
2130}
2132{
2133 DEBUG_ENTRY("ParseSpecial()");
2134 /* special key, can do anything */
2136}
2138{
2139 /* taumin command minimum optical depths for lines dafault 1e-20 */
2140 opac.taumin = (realnum)pow(10.,p.FFmtRead());
2141 if( p.lgEOL() )
2142 p.NoNumb("minimum optical depth");
2143}
2145{
2146 /* read in title of model starting in col 5 -- prefer to get string
2147 in quotes, but for the moment if it's not present take what you
2148 can get */
2149 if (p.GetQuote(input.chTitle,false) != 0)
2150 strcpy( input.chTitle , p.getRawTail().c_str()+1 );
2151}
2153{
2154 DEBUG_ENTRY("ParseTolerance()");
2155 fprintf(ioQQQ,
2156 "Sorry, this command has been replaced with the SET TEMPERATURE TOLERANCE command.\n");
2158}
2160{
2161 /* velocity power law as a function of radius. */
2162 DoppVel.TurbVelLaw = (realnum)p.FFmtRead();
2163
2164 DoppVel.lgTurbLawOn = true;
2166 /* for now, insist upon non-positive power,
2167 * so that velocity decreases with increasing radius. */
2168 ASSERT( DoppVel.TurbVelLaw <= 0.f );
2169}
2171{
2172 DEBUG_ENTRY("ParseTurbulence()");
2173 string ExtraPars;
2174
2175 if( p.nMatch("EQUIPART") )
2176 {
2177 /* turbulence equipartition option */
2178 DoppVel.lgTurbEquiMag = true;
2179
2180 /* this is optional F parameter from equation 34 of
2181 *>>refer pressure turb Heiles, C. & Troland, T.H. 2005, ApJ, 624, 773
2182 * turbulent energy density will be rho F v^2 / 2 */
2183 DoppVel.Heiles_Troland_F = (realnum)p.FFmtRead();
2184 if( p.lgEOL() )
2185 {
2186 /* this is the default value of 3 for isotropic turbulence */
2187 DoppVel.Heiles_Troland_F = 3.f;
2188 }
2189 }
2190 else
2191 {
2192 /* line has turbulent velocity in km/s */
2193 DoppVel.lgTurbEquiMag = false;
2194 DoppVel.TurbVel = (realnum)p.FFmtRead();
2195 if( p.lgEOL() )
2196 p.NoNumb("microturbulent velocity");
2197
2198 if( p.nMatch(" LOG") )
2199 {
2200 if( DoppVel.TurbVel > 32. )
2201 {
2202 fprintf( ioQQQ, "PROBLEM the log of the turbulence is "
2203 "%.2e - I cannot handle a number this big.\n",
2204 DoppVel.TurbVel );
2205 fprintf( ioQQQ, " The line image was\n" );
2206 p.PrintLine(ioQQQ);
2207 fprintf( ioQQQ, " Sorry.\n" );
2209 }
2210 DoppVel.TurbVel = (realnum)pow((realnum)10.f,DoppVel.TurbVel);
2211 }
2212 /* now convert from km/s to cm/s */
2213 DoppVel.TurbVel *= 1e5;
2214
2215 if( DoppVel.TurbVel < 0. )
2216 {
2217 fprintf( ioQQQ, " PROBLEM: the turbulent velocity needs to be > 0, but this was entered: %e\n",
2218 DoppVel.TurbVel );
2219 fprintf( ioQQQ, " Bailing out. Sorry.\n" );
2221 }
2222 else if( DoppVel.TurbVel >= SPEEDLIGHT )
2223 {
2224 fprintf( ioQQQ, " PROBLEM: A turbulent velocity greater than speed of light is not allowed, this was entered: %e\n",
2225 DoppVel.TurbVel );
2226 fprintf( ioQQQ, " Bailing out. Sorry.\n" );
2228 }
2229
2230 /* this is optional F parameter from equation 34 of
2231 *>>refer pressure turb Heiles, C. & Troland, T.H. 2005, ApJ, 624, 773
2232 * turbulent energy density will be rho F v^2 / 2 */
2233 DoppVel.Heiles_Troland_F = (realnum)p.FFmtRead();
2234 if( p.lgEOL() )
2235 {
2236 /* this is the default value of 3 for isotropic turbulence */
2237 DoppVel.Heiles_Troland_F = 3.f;
2238 }
2239
2240 /* option to have turbulence be dissipative, keyword is dissipate,
2241 * argument is log of scale length in cm */
2242 if( p.nMatch("DISS") )
2243 {
2244 DoppVel.DispScale = (realnum)pow(10., p.FFmtRead() );
2245 if( p.lgEOL() )
2246 p.NoNumb("turbulence dissipation scale");
2247 ExtraPars += " DISSIPATE %f";
2248 }
2249 }
2250
2251 /* include turbulent pressure in equation of state?
2252 * >>chng 06 mar 24, include turbulent pressure in gas equation of state by default,
2253 * but not if NO PRESSURE occurs */
2254 if( p.nMatch(" NO ") && p.nMatch("PRES") )
2255 {
2256 DoppVel.lgTurb_pressure = false;
2257 ExtraPars += " NO PRESSURE";
2258 }
2259 else
2260 {
2261 /* default is to include turbulent pressure in equation of state */
2262 DoppVel.lgTurb_pressure = true;
2263 }
2264
2265 /* vary option */
2266 if( optimize.lgVarOn && !p.nMatch("EQUIPART") )
2267 {
2268 /* only one parameter to vary */
2269 optimize.nvarxt[optimize.nparm] = 2;
2270 // the line image
2271 strcpy( optimize.chVarFmt[optimize.nparm], "TURBULENCE= %f LOG %f" );
2272 strcat( optimize.chVarFmt[optimize.nparm], ExtraPars.c_str() );
2273
2274 /* pointer to where to write */
2275 optimize.nvfpnt[optimize.nparm] = input.nRead;
2276 /* turbulent velocity */
2277 optimize.vparm[0][optimize.nparm] = log10(DoppVel.TurbVel/1e5f);
2278 optimize.vparm[1][optimize.nparm] = DoppVel.Heiles_Troland_F;
2279 if( p.nMatch("DISS") )
2280 {
2281 optimize.nvarxt[optimize.nparm] = 3;
2282 optimize.vparm[2][optimize.nparm] = log10(DoppVel.DispScale);
2283 }
2284 optimize.vincr[optimize.nparm] = 0.1f;
2285 ++optimize.nparm;
2286 }
2287
2288 /* set initial turbulence */
2289 DoppVel.TurbVelZero = DoppVel.TurbVel;
2290}
t_abund abund
Definition abund.cpp:5
t_atmdat atmdat
Definition atmdat.cpp:6
void ParseAtomFeII(Parser &p)
t_called called
Definition called.cpp:5
FILE * ioQQQ
Definition cddefines.cpp:7
#define ASSERT(exp)
Definition cddefines.h:578
#define REALLOC
Definition cddefines.h:519
#define MIN2
Definition cddefines.h:761
#define MIN3(a, b, c)
Definition cddefines.h:766
const int INPUT_LINE_LENGTH
Definition cddefines.h:254
#define EXIT_FAILURE
Definition cddefines.h:140
#define cdEXIT(FAIL)
Definition cddefines.h:434
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
const int ipHYDROGEN
Definition cddefines.h:305
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
long min(int a, long b)
Definition cddefines.h:723
bool getline(void)
Definition parser.cpp:164
double FFmtRead(void)
Definition parser.cpp:353
bool nMatch(const char *chKey) const
Definition parser.h:135
void setline(const char *const card)
Definition parser.h:69
bool Command(const char *name, OptionParser doOpts)
Definition parser.h:187
void set_point(long int ipnt)
Definition parser.h:77
double getNumberCheckAlwaysLog(const char *chDesc)
Definition parser.cpp:308
bool isCommandComment(void) const
Definition parser.cpp:97
bool m_lgDSet
Definition parser.h:42
bool last(void) const
Definition parser.h:200
bool lgEOL(void) const
Definition parser.h:98
NORETURN void NoNumb(const char *chDesc) const
Definition parser.cpp:233
void echo(void) const
Definition parser.cpp:147
long int m_nInitFile
Definition parser.h:41
bool m_lgEOF
Definition parser.h:42
int GetQuote(char *chLabel, bool lgABORT)
Definition parser.h:209
double getNumberDefaultAlwaysLog(const char *chDesc, double fdef)
Definition parser.cpp:327
long int m_nqh
Definition parser.h:41
string getRawTail()
Definition parser.h:220
void doSetVar(void)
Definition parser.cpp:119
NORETURN void CommandError(void) const
Definition parser.cpp:154
bool isVar(void) const
Definition parser.cpp:102
int PrintLine(FILE *fp) const
Definition parser.h:204
void help(FILE *fp) const
Definition parser.cpp:182
static t_version & Inst()
Definition cddefines.h:175
t_conv conv
Definition conv.cpp:5
t_cosmology cosmology
Definition cosmology.cpp:11
realnum GetDensity(realnum z)
Definition cosmology.cpp:31
static t_cpu cpu
Definition cpu.h:355
t_dark_matter dark
t_dense dense
Definition dense.cpp:24
double dense_tabden(double r0, double depth)
double dense_fabden(double radius, double depth)
double dense_parametric_wind(double rad)
t_DoppVel DoppVel
Definition doppvel.cpp:5
void ParseDynaWind(Parser &p)
void ParseDynaTime(Parser &p)
t_fudgec fudgec
Definition fudgec.cpp:5
#define NFUDGC
Definition fudgec.h:9
t_geometry geometry
Definition geometry.cpp:5
GrainVar gv
Definition grainvar.cpp:5
t_grid grid
Definition grid.cpp:5
diatomics h2("h2", 4100., &hmi.H2_total, Yan_H2_CS)
t_hextra hextra
Definition hextra.cpp:5
t_input input
Definition input.cpp:12
const int ipHE_LIKE
Definition iso.h:63
const int ipH_LIKE
Definition iso.h:62
t_iterations iterations
Definition iterations.cpp:5
void ParseMagnet(Parser &p)
Definition magnetic.cpp:156
t_mole_global mole_global
Definition mole.cpp:6
void InitMonitorResults(void)
void ParseMonitorResults(Parser &p)
t_opac opac
Definition opacity.cpp:5
t_optimize optimize
Definition optimize.cpp:5
const long LIMPAR
Definition optimize.h:61
void ParseCMB(double z, long int *nqh)
Definition parse_CMB.cpp:11
void ParseAbsMag(Parser &p)
void ParseAbundances(Parser &p)
void ParseAge(Parser &p)
Definition parse_age.cpp:39
void ParseAgn(Parser &p)
Definition parse_agn.cpp:10
void ParseAtomH2(Parser &p)
void ParseAtomISO(long ipISO, Parser &p)
void ParseBackgrd(Parser &p)
void ParseBlackbody(Parser &p)
void ParseCaseB(Parser &p)
void ParseDarkMatter(Parser &p)
void ParseQH(Parser &p)
void ParseIntensity(Parser &p)
void ParseHydrogen(Parser &)
void ParseCommands(void)
void ParseSpecial(Parser &)
void ParseNuL_nu(Parser &p)
void ParseHelp(Parser &)
void ParseLaser(Parser &p)
void ParseConvHighT(Parser &)
void ParseTurbulence(Parser &p)
void ParseFail(Parser &p)
void ParseVLaw(Parser &p)
void ParseLuminosity(Parser &p)
void ParseDielectronic(Parser &)
void ParseBremsstrahlung(Parser &p)
void ParseRoberto(Parser &)
void ParseCMBOuter(Parser &p)
void ParseHExtra(Parser &p)
void ParseInitCount(Parser &p)
void ParseEnergy(Parser &p)
void ParseNuF_nu(Parser &p)
void ParseDoubleTau(Parser &)
void ParseHeLike(Parser &)
void ParseCosm(Parser &p)
void ParseFudge(Parser &p)
void ParseL_nu(Parser &p)
void ParseGravity(Parser &p)
void ParseAbundancesNonSolar(Parser &p)
void ParseDistance(Parser &p)
void ParseForceTemperature(Parser &p)
void ParseAtom(Parser &p)
void ParseF_nuSpecific(Parser &p)
void ParseCovering(Parser &p)
void ParseCylinder(Parser &p)
void ParseCExtra(Parser &p)
void ParseTauMin(Parser &p)
void ParseNeutrons(Parser &p)
void ParseTitle(Parser &)
void ParseAperture(Parser &p)
void ParseDiffuse(Parser &p)
void ParsePGrains(Parser &)
void ParsePhi(Parser &p)
void ParseFill(Parser &p)
void ParseTolerance(Parser &)
void ParseIterations(Parser &p)
void ParseEden(Parser &p)
void ParseCompile(Parser &p)
void ParseConstant(Parser &p)
void ParseCoronal(Parser &p)
void ParseCosmicRays(Parser &p)
void ParseCosmology(Parser &p)
void ParseCrashDo(Parser &p)
void ParseDLaw(Parser &p)
void ParseDont(Parser &p)
void ParseDrive(Parser &p)
void ParseElement(Parser &p)
void ParseExtinguish(Parser &p)
void ParseF_nu(Parser &p, const char *chType, bool lgNU2)
void ParseFluc(Parser &p)
Definition parse_fluc.cpp:8
void ParseGlobule(Parser &p)
void ParseGrain(Parser &p)
void ParseGrid(Parser &p)
void ParseHDEN(Parser &p)
void ParseIlluminate(Parser &p)
void ParseInit(Parser &p)
Definition parse_init.cpp:9
void ParseInterp(Parser &p)
void ParseIonParI(Parser &p)
void ParseIonParX(Parser &p)
void ParseMap(Parser &p)
Definition parse_map.cpp:9
void ParseMetal(Parser &p)
void ParseNorm(Parser &p)
void ParseOptimize(Parser &p)
void ParsePlot(Parser &p)
void ParsePowerlawContinuum(Parser &p)
void ParsePrint(Parser &p)
void ParseRadius(Parser &p)
void ParseRangeOption(Parser &p)
void ParseRatio(Parser &p)
void ParseSave(Parser &p)
void ParseSet(Parser &p)
Definition parse_set.cpp:44
void ParseSphere(Parser &p)
void ParseState(Parser &p)
void ParseStop(Parser &p)
void ParseTable(Parser &p)
void ParseTest(Parser &p)
void ParseTLaw(Parser &p)
void ParseTrace(Parser &p)
t_phycon phycon
Definition phycon.cpp:6
UNUSED const double SOLAR_MASS
Definition physconst.h:71
UNUSED const double SPEEDLIGHT
Definition physconst.h:100
UNUSED const double PARSEC
Definition physconst.h:138
t_plotCom plotCom
Definition plot.cpp:19
t_pressure pressure
Definition pressure.cpp:5
t_radius radius
Definition radius.cpp:5
t_rfield rfield
Definition rfield.cpp:8
const int LIMSPC
Definition rfield.h:18
t_rt rt
Definition rt.cpp:5
t_StopCalc StopCalc
Definition stopcalc.cpp:5
const realnum COLUMN_INIT
Definition stopcalc.h:14
OptionParser action
Definition parser.h:25
const char * name
Definition parser.h:24
t_thermal thermal
Definition thermal.cpp:5
t_trace trace
Definition trace.cpp:5
Wind wind
Definition wind.cpp:5