cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
parse_commands.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 /*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 "h2.h"
20 #include "mole.h"
21 #include "hmi.h"
22 #include "rt.h"
23 #include "thermal.h"
24 #include "opacity.h"
25 #include "atomfeii.h"
26 #include "called.h"
27 #include "wind.h"
28 #include "hextra.h"
29 #include "iterations.h"
30 #include "radius.h"
31 #include "input.h"
32 #include "assertresults.h"
33 #include "phycon.h"
34 #include "fudgec.h"
35 #include "version.h"
36 #include "conv.h"
37 #include "parse.h"
38 
39 void ParseCommands(void)
40 {
41  char chCard[INPUT_LINE_LENGTH],
42  chKey2[3],
43  chKey3[4],
44  chKey4[5],
45  chKey5[6];
46 
47  bool lgDSet,
48  lgEOF,
49  lgEOL,
50  lgNu2;
51  bool lgStop ,
52  lgStop_not_enough_info;
53 
54  long int i,
55  j,
56  nqh;
57  long nInitFile=0;/* used to count number of init files read in */
58 
59  realnum a,
60  /* \todo 1 ar1 saves initial radius but does not appear
61  * to be used for anything. all uses are on left side of =
62  * this could be removed across the code. But first check
63  * that it is indeed not ever used */
64  ar1,
65  teset,
66  z;
67  double thickness;
68 
69  DEBUG_ENTRY( "ParseCommands()" );
70 
71  /* following says abundances are solar */
72  abund.lgAbnSolar = true;
73 
74  /* there are no plots desired yet */
75  plotCom.lgPlotON = false;
76  plotCom.nplot = 0;
77 
78  /* this flag remembers whether grains have ever been turned on. It is passed
79  * to routine ParseAbundances - there grains will be turned on with commands
80  * such as abundances ism, unless grains were previously set
81  * with a grains command. this way abundances will not clobber explicitly set
82  * grains. */
83  lgDSet = false;
84 
85  radius.Radius = 0.;
86  radius.rinner = 0.;
87  nqh = 0;
88  rfield.nspec = 0;
89 
90  /* will be set true if no buffering command entered */
91  input.lgSetNoBuffering = false;
92 
93  /* initialize some code variables in case assert command used in input stream */
95 
96  for( i=0; i < LIMSPC; i++ )
97  {
98  strcpy( rfield.chRSpec[i], "UNKN" );
99  }
100  optimize.nparm = 0;
101 
102  /* this is an option to turn on trace printout on the nth
103  * call from the optimizer */
104  if( optimize.lgTrOpt )
105  {
106  /* nTrOpt was set with "optimize trace" command,
107  * is iteration to turn on trace */
109  {
110  trace.lgTrace = true;
111  called.lgTalk = true;
112  trace.lgTrOvrd = true;
113  fprintf( ioQQQ, " READR turns on trace from optimize option.\n" );
114  }
115  }
116 
117  /* say this is a beta version if we are talking and it is the truth */
118  if( t_version::Inst().nBetaVer > 0 && called.lgTalk )
119  {
120  fprintf( ioQQQ,
121  "\n This is a beta release of Cloudy, and is intended for testing only.\n\n" );
122  }
123 
124  if( called.lgTalk )
125  {
126  /* this code prints pretty lines at top of output box */
127  int indent = (int)((122 - strlen(t_version::Inst().chVersion))/2);
128  fprintf( ioQQQ, "%*cCloudy %s\n", indent, ' ', t_version::Inst().chVersion );
129 
130  fprintf( ioQQQ, "%57cwww.nublado.org\n\n", ' ' );
131 
132  /* now print box and date of version, before printing commands */
133  fprintf( ioQQQ, "%23c", ' ' );
134  fprintf( ioQQQ, "**************************************");
135  fprintf( ioQQQ, "%7.7s", t_version::Inst().chDate);
136  fprintf( ioQQQ, "**************************************\n");
137 
138  fprintf( ioQQQ, "%23c*%81c*\n", ' ', ' ' );
139  }
140 
141  /* read in commands and print them */
142 
143  /* following signals which read is in progress,
144  * 1 is forward read, of input commands
145  * -1 is reverse read from bottom of line array, of cloudy.ini file */
146  input.iReadWay = 1;
147 
148  /* initialize array reader, this sub does nothing but set
149  * initial value of a variable, depending on value of iReadWay */
150  input_init();
151 
152  /* read until eof or blank line, then return control back to main program */
153 
154  /*
155  * input_readarray puts the next stored line image into chCard
156  * it will be <=INPUT_LINE_LENGTH char long, with eol at end
157  * it will be in mixed upper and lower case
158  * lgEOF is set true if we hit eof and no more lines present
159  * code moving over to info contained in structures,
160  *
161  * input.chOrgCard == original image of command line
162  * input.chCARDCAPS == original image converted to caps
163  */
164  input_readarray(chCard,&lgEOF);
165 
166  /* line starting with blank is eof, this checks we are not at eof */
167  while( !lgEOF && chCard[0] != ' ' )
168  {
169  /* echo the line before we cap it, but only if it does
170  * not contain the keyword HIDE */
171  /* >>chng 04 jan 21, add HIDE option, mostly for print quiet command */
172  if( called.lgTalk && !nMatch("HIDE",input.chCARDCAPS) )
173  fprintf( ioQQQ, "%23c* %-80s*\n", ' ', chCard );
174 
175  /* change chCard to all caps */
176  caps( chCard );
177 
178  /* now set up several partial keys */
179  /* first two characters into chKey2 */
180  strncpy( chKey2 , chCard , 2 );
181  chKey2[2] = '\0';
182 
183  /* first three characters into chKey3 */
184  strncpy( chKey3 , chCard , 3 );
185  chKey3[3] = '\0';
186 
187  /* first four characters into chKey4 */
188  strncpy( chKey4 , chCard , 4 );
189  chKey4[4] = '\0';
190 
191  /* first four characters into chKey4 */
192  strncpy( chKey5 , chCard , 5 );
193  chKey5[5] = '\0';
194 
195  /* check whether VARY is on line */
196  if( nMatch("VARY",chCard) )
197  {
198  optimize.lgVarOn = true;
199  if( optimize.nparm + 1 > LIMPAR )
200  {
201  fprintf( ioQQQ, " Too many VARY lines entered; the limit is%4ld\n",
202  LIMPAR );
203  cdEXIT(EXIT_FAILURE);
204  }
205  }
206 
207  else
208  {
209  optimize.lgVarOn = false;
210  }
211 
212  /* look for "C " comment - this form of test should not be necessary since
213  * input line is padded with extra space - but works and would continue to
214  * work if padding were removed */
215  if( chCard[0]=='C' && (chCard[1]==' ' || chCard[1]== 0) )
216  {
217  /* this is null test since lines starting with "C " are comments,
218  * the char after the c can be a space or a newline */
219  i = 1;
220  }
221 
222  /* start to look for keywords */
223  else if( strcmp(chKey4,"ABSO") == 0 )
224  {
225  /* enter luminosity in absolute magnitudes, in reads2 */
226  ParseAbsMag(chCard,&nqh);
227  }
228 
229  else if( strcmp(chKey3,"AGE") == 0 )
230  {
231  /* enter age of cloud so we can check for time-steady reads2 */
232  ParseAge(chCard);
233  }
234 
235  else if( strcmp(chKey4,"AGN ") == 0 )
236  {
237  /* enter generic style AGN continuum, in reads2 */
238  ParseAgn(chCard);
239  }
240 
241  else if( strcmp(chKey4,"ABUN") == 0 )
242  {
243  /* chemical abundances, readsun, lgDSet is set true if grains command
244  * ever entered, so that abundances command does not override grains command */
245  ParseAbundances(chCard,lgDSet);
246  /* abundances no longer solar */
247  abund.lgAbnSolar = false;
248  }
249 
250  else if( strcmp(chKey4,"APER") == 0 )
251  {
252  /* aperture command to simulate pencil beam or long slit */
253 
254  /* if the "BEAM" or "SLIT" option is specified then only part
255  * of the geometry is observed, and intensities
256  * should not be weighted by r^2. There are two limiting cases, SLIT in which
257  * the slit is longer than the diameter of the nebula and the contribution to the
258  * detected luminosity goes as r^1, and BEAM when the contribution is r^0,
259  * the same as plane parallel
260  *
261  * default value of geometry.iEmissPower is 2 (set in zero.c) for full geometry
262  */
263  if( nMatch("SLIT",chCard) )
264  {
265  /* long slit is case where slit is longer than diameter, so emissivity contributes
266  * r^1 to the observed luminosity */
267  geometry.iEmissPower = 1;
268  }
269  else if( nMatch("BEAM",chCard) )
270  {
271  /* long slit is case where slit is longer than diameter, so emissivity contributes
272  * r^1 to the observed luminosity */
273  /* this is the default or SHORT case, where r^0 is dependence */
274  geometry.iEmissPower = 0;
275  }
276  else
277  {
278  fprintf( ioQQQ, " One of the keywords SLIT or BEAM must appear.\n" );
279  fprintf( ioQQQ, " Sorry.\n" );
280  cdEXIT(EXIT_FAILURE);
281  }
282  }
283 
284  else if( strcmp(chKey4,"ASSE") == 0 )
285  {
286  /* assert that code predicts certain results, in assertresults.h */
288  }
289 
290  else if( strcmp(chKey4,"ATOM") == 0 )
291  {
292  /* accept both forms of feii */
293  if( nMatch("FEII",chCard) || nMatch("FE II",chCard) )
294  {
295  /* parse the atom FeII command - routine in atom_feii.c */
296  ParseAtomFeII(chCard);
297  }
298 
299  else if( nMatch("H-LI",chCard) )
300  {
301  /* parse the atom h-like command */
302  ParseAtomISO(ipH_LIKE, chCard);
303  }
304 
305  else if( nMatch("HE-L",chCard) )
306  {
307  /* parse the atom he-like command */
308  ParseAtomISO(ipHE_LIKE, chCard);
309  }
310 
311  else if( nMatch("ROTO",chCard) )
312  {
313  /* ATOM ROTOR same as ATOM CO */
314  fprintf(ioQQQ,"PROBLEM - the atom rotor command is now the ATOM CO command. "
315  "Please use that instead. I will accept this for now but may not in future versions.\n");
316  ParseAtomCO(chCard);
317  }
318 
319  else if( nMatch(" CO ",chCard) )
320  {
321  /* parameters for the rigid rotor CO molecules (yes, its not an atom)
322  * command is ATOM CO */
323  ParseAtomCO(chCard);
324  }
325 
326  else if( nMatch(" H2 ",chCard) )
327  {
328  /* parameters for the rigid rotor CO molecules (yes, its not an atom)
329  * command is ATOM ROTOR, routine is stored in parse_atomh2 */
330  ParseAtomH2(chCard);
331  }
332 
333  else
334  {
335  fprintf( ioQQQ, " I could not recognize a keyword on this atom command.\n");
336  fprintf( ioQQQ, " The available keys are FeII, H-Like, He-like, rotor and H2.\n");
337  fprintf( ioQQQ, " Sorry.\n" );
338  cdEXIT(EXIT_FAILURE);
339  }
340  }
341 
342  else if( strcmp(chKey4,"BACK") == 0 )
343  {
344  /* cosmic background, in parse_backgrd */
345  ParseBackgrd(&nqh,chCard,&ar1);
346  }
347 
348  else if( strcmp(chKey4,"BLAC") == 0 )
349  {
350  /* black body, in reads2 */
351  ParseBlackbody(chCard,&nqh,&ar1);
352 
353  /* vary option */
354  if( optimize.lgVarOn )
355  {
356  /* no luminosity options on vary */
358  strcpy( optimize.chVarFmt[optimize.nparm], "BLACKbody=%f" );
359  /* pointer to where to write */
361  /* log of temp stored here */
362  /* >>chng 02 feb 11, log was not here, don't know what happened to it */
364  /* the increment in the first steps away from the original value */
365  optimize.vincr[optimize.nparm] = 0.5;
366  optimize.nparm += 1;
367  }
368  }
369 
370  else if( strcmp(chKey4,"BREM") == 0 )
371  {
372  /* bremsstrahlung continuum from central object */
373  strcpy( rfield.chSpType[rfield.nspec], "BREMS" );
374  i = 5;
376  (realnum)FFmtRead(chCard,&i, INPUT_LINE_LENGTH,&lgEOL);
377  if( lgEOL )
378  {
379  NoNumb(chCard);
380  }
381 
382  /* temperature is interpreted as log if <=10, linear otherwise*/
383  if( rfield.slope[rfield.nspec] <= 10. )
384  {
386  pow(10.,rfield.slope[rfield.nspec]);
387  }
388  rfield.cutoff[rfield.nspec][0] = 0.;
389 
390  /* option for vary keyword */
391  if( optimize.lgVarOn )
392  {
393  /* only one parameter */
395  strcpy( optimize.chVarFmt[optimize.nparm], "BREMS, T=%f" );
396  /* pointer to where to write */
398  /* log of temp will be pointer */
400  optimize.vincr[optimize.nparm] = 0.5;
401  ++optimize.nparm;
402  }
403  ++rfield.nspec;
404  if( rfield.nspec >= LIMSPC )
405  {
406  /* too many continua were entered */
407  fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" );
408  cdEXIT(EXIT_FAILURE);
409  }
410  }
411 
412  else if( strcmp(chKey4,"CASE") == 0 )
413  {
414  /* do hydrogen case b */
415  ParseCaseB( chCard );
416  }
417 
418  else if( strcmp(chKey4,"CEXT") == 0 )
419  {
420  /* add "extra" cooling, power law temp dependence */
421  thermal.lgCExtraOn = true;
422  i = 5;
423  thermal.CoolExtra = (realnum)pow(10.,FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL));
424  if( lgEOL )
425  {
426  NoNumb(chCard);
427  }
428  thermal.cextpw = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
429  }
430 
431  /* replace with CMB command, both will continue to be parsed */
432  else if( (strcmp(chKey3,"CMB") == 0) || (strcmp(chKey4,"FIRE") == 0) )
433  {
434  /* cosmic thermal background radiation, argument is redshift */
435  i = 5;
436  /* if no number on line then (realnum)FFmtRead returns z=0; i.e., now */
437  z = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
438  /* in readsun */
439  ParseCMB(z,&nqh,&ar1);
440  }
441 
442  else if( strcmp(chKey4,"COMP") == 0 )
443  {
444  /* compile ascii version of stellar atmosphere continua in volk */
445  ParseCompile(chCard);
446  }
447 
448  else if( strcmp(chKey4,"CONS") == 0 )
449  {
450  /* constant temperature, pressure, density, or gas pressure
451  * in readsun */
452  ParseConstant(chCard);
453  }
454 
455  else if( strcmp(chKey4,"CORO") == 0 )
456  {
457  /* coronal equilibrium; set constant temperature to number on line
458  * in readsun */
459  ParseCoronal( chCard,&nqh,&ar1);
460  }
461 
462  else if( strcmp(chKey4,"COSM") == 0 )
463  {
464  /* cosmic ray density, log of rate relative to background, or density of rel. electrons,
465  * quantity is log unless keyword linear appears */
466  ParseCosmicRays( chCard );
467  }
468 
469  else if( strcmp(chKey4,"COVE") == 0 )
470  {
471  /* covering factor for gas */
472  i = 5;
473  /* The geometric covering factor accounts for how much of 4\pi is
474  * covered by gas, and so linearly multiplies the predicted intensities */
475  geometry.covgeo = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
476 
477  if( lgEOL )
478  {
479  NoNumb(chCard);
480  }
481 
482  /* if negative then log, convert to linear */
483  if( geometry.covgeo <= 0. )
484  {
486  }
487 
488  if( geometry.covgeo > 1. )
489  {
490  fprintf( ioQQQ, " A covering factor greater than 1 makes no physical sense. Sorry.\n" );
491  cdEXIT(EXIT_FAILURE);
492  }
493 
494  /* radiative transfer covering factor will be equal to the geometric one */
496  }
497 
498  else if( strcmp(chKey4,"CRAS") == 0 )
499  {
500  /* any of several tests to cause the code to crash */
501  ParseCrashDo(chCard);
502  }
503 
504  else if( strcmp(chKey4,"CYLI") == 0 )
505  {
506  /* height of cylinder in cm */
507  i = 5;
508  radius.lgCylnOn = true;
509  radius.CylindHigh = pow(10.,FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL));
510  if( lgEOL )
511  {
512  NoNumb(chCard);
513  }
514  }
515 
516  else if( strcmp(chKey4,"DIEL") == 0 )
517  {
518  /* >>chng 05 dec 21, change dielectronic command to set dielectronic recombination command */
519  fprintf( ioQQQ, " The DIELectronic command has been replaced with the SET DIELectronic recombination command.\n" );
520  fprintf( ioQQQ, " Please have a look at Hazy.\n Sorry.\n\n" );
521  cdEXIT(EXIT_FAILURE);
522  }
523 
524  else if( strcmp(chKey4,"DIFF") == 0 )
525  {
526  /* set method of transferring diffuse fields,
527  * default is outward only, cdDffTrns label "OU2", set in zero.c */
528  if( nMatch(" OTS",chCard) )
529  {
530  if( nMatch("SIMP",chCard) )
531  {
532  /* this is simple ots, which never allows photons to escape */
533  strcpy( rfield.chDffTrns, "OSS" );
534  }
535  else
536  {
537  /* this is regular ots, which allows photons to escape */
538  strcpy( rfield.chDffTrns, "OTS" );
539  }
540  rfield.lgOutOnly = false;
541  }
542  else if( nMatch(" OUT",chCard) )
543  {
544  rfield.lgOutOnly = true;
545  i = 4;
546  j = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
547  if( lgEOL )
548  {
549  /* this is the default set in zero */
550  strcpy( rfield.chDffTrns, "OU2" );
551  }
552  else
553  {
554  if( j > 0 && j < 10 )
555  {
556  sprintf( rfield.chDffTrns, "OU%1ld", j );
557  }
558  else
559  {
560  fprintf( ioQQQ, " must be between 1 and 9 \n" );
561  cdEXIT(EXIT_FAILURE);
562  }
563  }
564  }
565 
566  else
567  {
568  fprintf( ioQQQ, " There should have been OUTward or OTS on this line. Sorry.\n" );
569  cdEXIT(EXIT_FAILURE);
570  }
571  }
572 
573  else if( strcmp(chKey4,"DIST") == 0 )
574  {
575  /* distance to the object */
576  i = 5;
577  radius.distance = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
578  if( lgEOL )
579  {
580  NoNumb(chCard);
581  }
582  /* default is for quantity to be log of distance, linear keyword
583  * overrides this - is LINEAR is not on line then exp */
584  if( !nMatch("LINE",chCard ) )
585  {
586  radius.distance = pow(10., radius.distance );
587  }
588  /* default is radius in cm - if parsecs appears then must
589  * convert to cm */
590  if( nMatch("PARS",chCard ) )
591  {
593  }
594  }
595 
596  else if( strcmp(chKey4,"DLAW") == 0 )
597  {
598  /* either use dense_fabden routine, or read in table of points */
599  ParseDLaw(chCard);
600  }
601 
602  else if( strcmp(chKey4,"DOUB") == 0 )
603  {
604  /* double optical depth scale after each iteration */
605  rt.DoubleTau = 2.;
606  }
607 
608  else if( strcmp(chKey4,"DRIV") == 0 )
609  {
610  /* call one of several drivers, readsun */
611  ParseDrive(chCard);
612  }
613 
614  else if( strcmp(chKey4,"EDEN") == 0 )
615  {
616  /* option to add "extra" electrons, to test Compton limit
617  * for very low T(star) - option is log of eden */
618  i = 5;
619  dense.EdenExtra = (realnum)pow(10.,FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL));
620  if( lgEOL )
621  {
622  NoNumb(chCard);
623  }
624  /* warn that this model is meaningless */
625  phycon.lgPhysOK = false;
626  }
627 
628  else if( strcmp(chKey4,"ELEM") == 0 )
629  {
630  /* element command;
631  * scale or abundance options, to change abundance of specific element
632  * read option to change order of elements
633  * in reads2.f */
634  ParseElement(chCard);
635  }
636 
637  else if( strcmp(chKey4,"ENER") == 0 )
638  {
639  /* energy density (degrees K) of this continuum source */
640  i = 5;
641  if( nqh >= LIMSPC )
642  {
643  /* too many continua were entered */
644  fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" );
645  cdEXIT(EXIT_FAILURE);
646  }
647  /* silly, but calms down the lint */
648  ASSERT( nqh < LIMSPC );
649  strcpy( rfield.chRSpec[nqh], "SQCM" );
650  teset = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
651  if( lgEOL )
652  {
653  NoNumb(chCard);
654  }
655  /* set radius to very large value if not already set */
656  /* >>chng 01 jul 24, from Radius == 0 to this, as per PvH comments */
657  if( !radius.lgRadiusKnown )
658  {
659  /* set radius to large value */
660  ar1 = (realnum)radius.rdfalt;
661  radius.Radius = pow(10.,radius.rdfalt);
662  }
663 
664  /* convert temp to log, recognize linear option */
665  if( nMatch("LINE",chCard) || teset > 10. )
666  {
667  /* option for linear temperature, must store log */
668  teset = (realnum)log10(teset);
669  }
670 
671  if( teset > 5. )
672  {
673  fprintf( ioQQQ, " This intensity may be too large. The code may crash due to overflow. Was log intended?\n" );
674  }
675 
676  /* teset is not log of temp, now get log of total luminosity */
677  strcpy( rfield.chSpNorm[nqh], "LUMI" );
678 
679  /* full range of continuum will be used */
680  rfield.range[nqh][0] = rfield.emm;
681  rfield.range[nqh][1] = rfield.egamry;
682  rfield.totpow[nqh] = (4.*teset - 4.2464476 + 0.60206);
683 
684  /* >>chng 06 mar 22, add time option to vary only some continua with time */
685  if( nMatch( "TIME" , chCard ) )
686  rfield.lgTimeVary[nqh] = true;
687 
688  /* vary option */
689  if( optimize.lgVarOn )
690  {
691  strcpy( optimize.chVarFmt[optimize.nparm], "ENERGY DENSITY %f log " );
692  /* pointer to where to write */
694  optimize.vparm[0][optimize.nparm] = (realnum)log10(rfield.totpow[nqh]);
695  optimize.vincr[optimize.nparm] = 0.1f;
697  optimize.nparm += 1;
698  }
699 
700  /* finally increment number of continua */
701  ++nqh;
702  }
703 
704  else if( strcmp(chKey4,"EXTI") == 0 )
705  {
706  /* extinguish ionizing continuum by absorbing column AFTER
707  * setting luminosity or Q(H). First number is the column
708  * density (log), second number is leakage (def=0%)
709  * last number is lowest energy (ryd), last two may be omitted
710  * from right to left
711  *
712  * extinction is actually done in extin, which is called by ContSetIntensity */
713  ParseExtinguish( chCard );
714  }
715 
716  else if( strcmp(chKey4,"FAIL") == 0 )
717  {
718  /* reset number of temp failures allowed, default=20 */
719  i = 5;
720 
721  /* save previous value */
722  j = conv.LimFail;
723  conv.LimFail = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
724  if( lgEOL )
725  {
726  NoNumb(chCard);
727  }
728 
729  /* >>chng 01 mar 14, switch logic on maps */
730  /* ' map' flag, make map when failure, default is no map,
731  * second check is so that no map does not trigger a map */
732  if( nMatch(" MAP",chCard) && !nMatch(" NO ",chCard) )
733  {
734  conv.lgMap = true;
735  }
736 
737  /* complain if failures was increased above default */
738  if( conv.LimFail > j )
739  {
740  fprintf( ioQQQ, " This command should not be necessary.\n" );
741  fprintf( ioQQQ, " Please show this input stream to Gary Ferland if this command is really needed for this simulation.\n" );
742  }
743  }
744 
745  else if( strcmp(chKey4,"FILL") == 0 )
746  {
747  /* filling factor, power law exponent (default=1., 0.) */
748  i = 5;
749  a = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
750  if( lgEOL )
751  {
752  NoNumb(chCard);
753  }
754 
755  if( a <= 0. )
756  {
757  /* number less than or equal to 0, must have been entered as log */
758  geometry.FillFac = (realnum)pow((realnum)10.f,a);
759  }
760  else
761  {
762  /* number greater than zero, must have been the real thing */
763  geometry.FillFac = a;
764  if( geometry.FillFac > 1.0 )
765  {
766  if( called.lgTalk )
767  {
768  fprintf( ioQQQ, " Filling factor > 1, reset to 1\n" );
769  }
770  geometry.FillFac = 1.;
771  }
772  }
774 
775  /* option to have power law dependence, filpow set to 0 in zerologic */
776  geometry.filpow = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
777 
778  /* vary option */
779  if( optimize.lgVarOn )
780  {
781  strcpy( optimize.chVarFmt[optimize.nparm], "FILLING FACTOR= %f power=%f" );
782 
783  /* pointer to where to write */
786  optimize.vincr[optimize.nparm] = 0.5;
787 
788  /* power law dependence here, but cannot be varied */
791 
792  /* do not allow filling factor to go positive */
793  optimize.varang[optimize.nparm][0] = -1e38f;
794  optimize.varang[optimize.nparm][1] = 0.f;
795  optimize.nparm += 1;
796  }
797  }
798 
799  else if( strcmp(chKey4,"FLUC") == 0 )
800  {
801  /* rapid density fluctuations, in readsun */
802  ParseFluc(chCard);
803  }
804 
805  else if( strcmp(chKey4,"F(NU") == 0 )
806  {
807  /* this is the specific flux at nu
808  * following says F(nu) not nuF(nu) */
809  lgNu2 = false;
810  /* in reads2 */
811  ParseF_nu(chCard,&nqh,&ar1,"SQCM",lgNu2);
812  }
813 
814  else if( strcmp(chKey4,"FORC") == 0 )
815  {
816  /* force temperature of first zone, but don't keep constant
817  * allow to then go to nearest equilibrium
818  * log if < 10 */
819  i = 5;
820  thermal.ConstTemp = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
821  if( lgEOL )
822  {
823  NoNumb(chCard);
824  }
825 
826  if( nMatch(" LOG",chCard) || (thermal.ConstTemp <= 10. &&
827  !nMatch("LINE",chCard)) )
828  {
830  }
831 
832  /* low energy bounds of continuum array do not permit too-low a temp */
833  if( thermal.ConstTemp < 3. )
834  {
835  fprintf( ioQQQ, " TE reset to 3K: entered number too small.\n" );
836  thermal.ConstTemp = 3.;
837  }
838  }
839 
840  else if( strcmp(chKey4,"FUDG") == 0 )
841  {
842  /* enter a fudge factor */
843  i = 5;
844  for( j=0; j < NFUDGC; j++ )
845  {
846  fudgec.fudgea[j] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
847  /* fudgec.nfudge is number of factors on the line */
848  if( !lgEOL )
849  fudgec.nfudge = j+1;
850  }
851 
852  /* vary option */
853  if( optimize.lgVarOn )
854  {
855  /* no luminosity options on vary */
857  strcpy( optimize.chVarFmt[optimize.nparm], "FUDGE=%f" );
858  /* enter remaining parameters as constants */
859  char chHold[1000];
860  for( j=1; j<fudgec.nfudge; ++j )
861  {
862  sprintf( chHold , " %f" , fudgec.fudgea[j] );
863  strcat( optimize.chVarFmt[optimize.nparm] , chHold );
864  }
865 
866  /* pointer to where to write */
868  /* fudge factor stored here */
870  /* the increment in the first steps away from the original value */
871  optimize.vincr[optimize.nparm] = 0.5;
872  optimize.nparm += 1;
873  }
874  }
875 
876  else if( strcmp(chKey4,"GLOB") == 0 )
877  {
878  /* globule with density increasing inward
879  * parameters are outer density, radius of globule, and density power */
880  ParseGlobule(chCard);
881  }
882 
883  else if( (strcmp(chKey4,"GRAI") == 0 )||( strcmp(chKey4,"PGRA") == 0 ) )
884  {
885  /* read parameters dealing with grains, in reads2 */
886  ParseGrain(chCard,&lgDSet);
887  }
888 
889  else if( strcmp(chKey4,"GRID") == 0 )
890  {
891  /* option to run grid of models by varying certain parameters
892  * in reads2 */
893  ParseGrid(chCard);
894  }
895 
896  else if( strcmp(chKey4,"HDEN") == 0 )
897  {
898  /* parse the hden command to set the hydrogen density, in reads2 */
899  ParseHDEN(chCard);
900  }
901 
902  else if( strcmp(chKey4,"HELI") == 0 )
903  {
904  fprintf(ioQQQ,"Sorry, this command is replaced with ATOM HE-LIKE\n");
905  cdEXIT(EXIT_FAILURE);
906  }
907 
908  else if( strcmp(chKey4,"HEXT") == 0 )
909  {
910  /* "extra" heating rate, so that first= log(erg(cm-3, s-1))
911  * second optional number is scale radius, so that HXTOT = TurbHeat*SEXP(DEPTH/SCALE)
912  * if missing then constant heating.
913  * third option is depth from shielded face, to mimic irradiation from both sides*/
914  i = 5;
915  hextra.TurbHeat = (realnum)pow(10.,FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL));
916  if( lgEOL )
917  NoNumb( chCard );
918  /* save initial value in case reset with time dependent command */
920 
921  /* remember which type of scale dependence we will use */
922  const char *chHextraScale;
923  /* these are optional numbers on depth or density option */
924  realnum par1=0. , par2=0.;
925  long int nPar;
926  if( nMatch("DEPT" , chCard ) )
927  {
928  /* depend on depth */
929  hextra.lgHextraDepth = true;
930  chHextraScale = "DEPTH";
931  /* optional scale radius */
932  a = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
933  if( lgEOL )
934  NoNumb(chCard);
935  else
936  hextra.turrad = (realnum)pow((realnum)10.f,a);
937 
938  /* depth from shielded face, to mimic illumination from both sides
939  * may not be specified */
940  realnum b = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
941  if( lgEOL )
942  {
943  hextra.turback = 0.;
944  nPar = 2;
945  }
946  else
947  {
948  hextra.turback = (realnum)pow((realnum)10.f,b);
949  nPar = 3;
950  }
951  par1 = a;
952  par2 = b;
953  }
954  else if( nMatch("DENS" , chCard ) )
955  {
956  /* depends on density */
957  chHextraScale = "DENSITY";
958  hextra.lgHextraDensity = true;
959 
960  /* the log of the density */
961  a = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
962  /* if number not entered then unity is used, heating
963  * proportional to density */
964  hextra.HextraScaleDensity = (realnum)pow((realnum)10.f,a);
965  par1 = a;
966  par2 = 0;
967  nPar = 2;
968  }
969  else
970  {
971  /* constant heating */
972  chHextraScale = "";
973  nPar = 1;
974  }
975 
976  /* option to have heating vary with time */
977  if( nMatch( "TIME" , chCard ) )
978  hextra.lgTurbHeatVaryTime = true;
979 
980  /* vary option */
981  if( optimize.lgVarOn )
982  {
983  /* 1 to 3 options on hextra vary */
984  optimize.nvarxt[optimize.nparm] = nPar;
986  optimize.vparm[1][optimize.nparm] = par1;
987  optimize.vparm[2][optimize.nparm] = par2;
988 
989  strcpy( optimize.chVarFmt[optimize.nparm], "HEXTra %f " );
990  strcat( optimize.chVarFmt[optimize.nparm], chHextraScale );
991  while( nPar > 1 )
992  {
993  /* add extra spots to write parameters */
994  --nPar;
995  strcat( optimize.chVarFmt[optimize.nparm], " %f" );
996  }
998  strcat( optimize.chVarFmt[optimize.nparm], " TIME" );
999  /* pointer to where to write */
1001  /* the increment in the first steps away from the original value */
1002  optimize.vincr[optimize.nparm] = 0.1f;
1003  optimize.nparm += 1;
1004  }
1005  }
1006 
1007  else if( strcmp(chKey4,"HIGH") == 0 )
1008  {
1009  /* approach equilibrium from high te */
1010  thermal.lgTeHigh = true;
1011  }
1012 
1013  else if( strncmp( chCard ,"HYDROGEN",8) == 0 )
1014  {
1015  fprintf(ioQQQ," Sorry, this command has been replaced with the ATOM H-LIKE command.\n");
1016  cdEXIT(EXIT_FAILURE);
1017  }
1018 
1019  else if( strcmp(chKey4,"ILLU") == 0 )
1020  {
1021  /* illumination angle */
1022  i = 5;
1024  if( lgEOL )
1025  {
1026  NoNumb(chCard);
1027  }
1028  /* default is degrees, but radian key also there */
1029  if( nMatch("RADI",chCard) )
1030  {
1031  /* entered as radians, convert to degrees first */
1032  geometry.AngleIllumRadian /= (realnum)(PI/180.);
1033  }
1034  /* we now have degrees, make sure between 0 and 90 */
1035  if( geometry.AngleIllumRadian < 0. || geometry.AngleIllumRadian >= 90. )
1036  {
1037  fprintf( ioQQQ, " Angle of illumination must be between 0 and 90 degrees "
1038  "or 0 and pi/2 radians.\n" );
1039  cdEXIT(EXIT_FAILURE);
1040  }
1041  /* we really want 1. / COS( theta ) with theta in radians
1042  * first convert angle to radians */
1044 
1045  /* this is effective path - dTau_eff = dTau_normal * geometry.DirectionalCosin */
1047 
1048  /* vary option */
1049  if( optimize.lgVarOn )
1050  {
1051  /* no luminosity options on vary */
1053  strcpy( optimize.chVarFmt[optimize.nparm], "ILLUminate %f radians " );
1054  /* pointer to where to write */
1057  /* the increment in the first steps away from the original value */
1058  optimize.vincr[optimize.nparm] = 0.1f;
1059  optimize.nparm += 1;
1060  }
1061  }
1062 
1063  else if( strcmp(chKey4,"INIT") == 0 )
1064  {
1065  /* read cloudy.ini initialization file
1066  * need to read in file every time, since some vars
1067  * get reset in zero - would be unsafe not to read in again */
1068  ParseInit(chCard);
1069 
1070  /* check that only one init file was in the input stream -
1071  * we cannot now read more than one */
1072  ++nInitFile;
1073  if( nInitFile > 1 )
1074  {
1075  fprintf( ioQQQ,
1076  " This is the second init file, I can only handle one.\nSorry.\n" );
1077  cdEXIT(EXIT_FAILURE);
1078  }
1079 
1080  /* we will put the data for the ini file at the end of the array of
1081  * line images, this is the increment for stuffing the lines in - negative */
1082  input.iReadWay = -1;
1083 
1084  /* initialize array reader, this sub does nothing but set
1085  * initial value of a variable, depending on value of iReadWay */
1086  input_init();
1087  }
1088 
1089  else if( strcmp(chKey5,"INTEN") == 0 )
1090  {
1091  /* intensity of this continuum source */
1092  i = 5;
1093  if( nqh >= LIMSPC )
1094  {
1095  /* too many continua were entered */
1096  fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" );
1097  cdEXIT(EXIT_FAILURE);
1098  }
1099 
1100  /* silly, but calms down the lint */
1101  ASSERT( nqh < LIMSPC );
1102  strcpy( rfield.chRSpec[nqh], "SQCM" );
1103  rfield.totpow[nqh] = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
1104  if( lgEOL )
1105  {
1106  NoNumb(chCard);
1107  }
1108 
1109  /* set radius to very large value if not already set */
1110  /* >>chng 01 jul 24, from Radius == 0 to this, as per PvH comments */
1111  if( !radius.lgRadiusKnown )
1112  {
1113  /* set radius to large value */
1114  ar1 = (realnum)radius.rdfalt;
1115  radius.Radius = pow(10.,radius.rdfalt);
1116  }
1117  if( nMatch("LINE",chCard) )
1118  {
1119  /* silly, but calms down the lint */
1120  ASSERT( nqh < LIMSPC );
1121  /* option for linear input parameter */
1122  rfield.totpow[nqh] = log10(rfield.totpow[nqh]);
1123  }
1124  strcpy( rfield.chSpNorm[nqh], "LUMI" );
1125  /* ParseRangeOption in readsun */
1126  ParseRangeOption(nqh,chCard);
1127 
1128  /* >>chng 06 mar 22, add time option to vary only some continua with time */
1129  if( nMatch( "TIME" , chCard ) )
1130  rfield.lgTimeVary[nqh] = true;
1131 
1132  /* vary option */
1133  if( optimize.lgVarOn )
1134  {
1135  /* create the format we will use to write the command */
1136  /* >>chng 03 apr 30, had used log for range option, but if one or other
1137  * was 1 ryd, then became 0 in the log, which was misinterpreted as
1138  * one of the continuum bounds */
1139  /*strcpy( optimize.chVarFmt[optimize.nparm], "INTENSITY%f log range %f %f" );*/
1140  strcpy( optimize.chVarFmt[optimize.nparm], "INTENSITY %f range %f %f " );
1141  /* array index for this command */
1144  optimize.vincr[optimize.nparm] = 0.5;
1145  /* range option, but cannot be varied */
1146  /*optimize.vparm[1][optimize.nparm] = (realnum)log10(rfield.range[nqh][0]);
1147  optimize.vparm[2][optimize.nparm] = (realnum)log10(rfield.range[nqh][1]);*/
1148  /* >>change to linear option */
1149  optimize.vparm[1][optimize.nparm] = (realnum)rfield.range[nqh][0];
1150  optimize.vparm[2][optimize.nparm] = (realnum)rfield.range[nqh][1];
1152  ++optimize.nparm;
1153  }
1154  ++nqh;
1155  }
1156 
1157  else if( strcmp(chKey5,"INTER") == 0 )
1158  {
1159  /* interpolate on input tables for continuum, set of power laws used
1160  * input ordered pairs nu( ryd or log(Hz) ), log(fnu)
1161  * additional lines begin CONTINUE
1162  * first check that this is the one and only INTERP command
1163  * in readsun */
1164  ParseInterp(chCard,&lgEOF);
1165  }
1166 
1167  else if( strcmp(chKey4,"IONI") == 0 )
1168  {
1169  /* inter ionization parameter U=Q/12 R*R N C;
1170  * defined per hydrogen, not per electron (as before)
1171  * radius must also be entered if spherical, but not needed if plane */
1172  ParseIonPar(&nqh,chCard,&ar1);
1173  }
1174 
1175  else if( strcmp(chKey4,"ITER") == 0 )
1176  {
1177  /* iterate to get optical depths and diffuse field properly */
1178  i = 5;
1179  iterations.itermx = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL) - 1;
1181  /* >>chng 06 mar 19, malloc itrdim arrays
1182  * iterations.iter_malloc is -1 before space allocated and set to
1183  * itermx if not previously set */
1185  {
1186  long int iter_malloc_save = iterations.iter_malloc;
1187  /* add one for mismatch between array dim and iterations defn,
1188  * another few for safety */
1191  (size_t)iterations.iter_malloc*sizeof(long int) );
1192  geometry.nend = (long int*)REALLOC(geometry.nend ,
1193  (size_t)iterations.iter_malloc*sizeof(long int) );
1194  radius.router = (double*)REALLOC(radius.router ,
1195  (size_t)iterations.iter_malloc*sizeof(double) );
1196  /* >>chng 06 jun 07, need to set IterPrnt, router, and nend */
1197  for(j=iter_malloc_save; j<iterations.iter_malloc; ++j )
1198  {
1199  iterations.IterPrnt[j] = iterations.IterPrnt[iter_malloc_save-1];
1200  geometry.nend[j] = geometry.nend[iter_malloc_save-1];
1201  radius.router[j] = radius.router[iter_malloc_save-1];
1202  }
1203  }
1204 
1205  if( nMatch("CONV",chCard) )
1206  {
1207  /* option to keep iterating until it converges, checks on convergence
1208  * are done in update, and checked again in prtcomment */
1209  conv.lgAutoIt = true;
1210  /* above would have been legitimate setting of ITERMX, set to default 10 */
1211  if( lgEOL )
1212  {
1213  iterations.itermx = 10 - 1;
1214  }
1215  a = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
1216  /* change convergence criteria, preset in zero */
1217  if( !lgEOL )
1218  {
1219  conv.autocv = a;
1220  }
1221  }
1222  }
1223 
1224  else if( strcmp(chKey4,"L(NU") == 0 )
1225  {
1226  /* this is the specific luminosity at nu
1227  * following says n(nu) not nuF(nu) */
1228  lgNu2 = false;
1229  /* in reads2 */
1230  ParseF_nu(chCard,&nqh,&ar1,"4 PI",lgNu2);
1231  }
1232 
1233  else if( strcmp(chKey4,"LASE") == 0 )
1234  {
1235  /* mostly one frequency (a laser) to check gamma's,*/
1236  /* say the continuum type */
1237  strcpy( rfield.chSpType[rfield.nspec], "LASER" );
1238 
1239  /* scan off the laser's energy */
1240  i = 5;
1241  rfield.slope[rfield.nspec] = FFmtRead(chCard,&i, INPUT_LINE_LENGTH,&lgEOL);
1242 
1243  /* negative energies are logs */
1244  if( rfield.slope[rfield.nspec] <= 0. )
1245  {
1246  rfield.slope[rfield.nspec] =
1247  pow(10.,rfield.slope[rfield.nspec]);
1248  }
1249  if( lgEOL )
1250  {
1251  NoNumb(chCard);
1252  }
1253 
1254  /* next number is optional relative width of laser */
1255  rfield.cutoff[rfield.nspec][0] = FFmtRead(chCard,&i, INPUT_LINE_LENGTH,&lgEOL);
1256  if( lgEOL )
1257  {
1258  /* default width is 0.05 */
1259  rfield.cutoff[rfield.nspec][0] = 0.05;
1260  }
1261 
1262  /* increment counter making sure we still have space enough */
1263  ++rfield.nspec;
1264  if( rfield.nspec >= LIMSPC )
1265  {
1266  /* too many continua were entered */
1267  fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" );
1268  cdEXIT(EXIT_FAILURE);
1269  }
1270  }
1271 
1272  else if( strcmp(chKey4,"LUMI") == 0 )
1273  {
1274  /* luminosity of this continuum source */
1275  i = 5;
1276  if( nqh >= LIMSPC )
1277  {
1278  /* too many continua were entered */
1279  fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" );
1280  cdEXIT(EXIT_FAILURE);
1281  }
1282 
1283  strcpy( rfield.chRSpec[nqh], "4 PI" );
1284  rfield.totpow[nqh] = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
1285  if( lgEOL )
1286  {
1287  NoNumb(chCard);
1288  }
1289  if( nMatch("LINE",chCard) )
1290  {
1291  /* option for linear input parameter */
1292  rfield.totpow[nqh] = log10(rfield.totpow[nqh]);
1293  }
1294 
1295  strcpy( rfield.chSpNorm[nqh], "LUMI" );
1296 
1297  /* the solar option - number is log total solar luminosity */
1298  if( nMatch("SOLA",chCard) )
1299  {
1300  /* option to use log of total luminosity in solar units */
1301  rfield.range[nqh][0] = rfield.emm;
1302  rfield.range[nqh][1] = rfield.egamry;
1303  /* log of solar luminosity in erg s-1 */
1304  rfield.totpow[nqh] += 33.5827f;
1305  }
1306  else
1307  {
1308  /* check if range is present and parse it if it is - ParseRangeOption in readsun
1309  * this includes TOTAL keyword for total luminosity */
1310  ParseRangeOption(nqh,chCard);
1311  }
1312 
1313  /* >>chng 06 mar 22, add time option to vary only some continua with time */
1314  if( nMatch( "TIME" , chCard ) )
1315  rfield.lgTimeVary[nqh] = true;
1316 
1317  /* vary option */
1318  if( optimize.lgVarOn )
1319  {
1320  /* >>chng 03 apr 30, had used log for range option, but if one or other
1321  * was 1 ryd, then became 0 in the log, which was misinterpreted as
1322  * one of the continuum bounds */
1323  strcpy( optimize.chVarFmt[optimize.nparm], "LUMINOSITY %f range %f %f " );
1324  /* pointer to where to write */
1327  optimize.vincr[optimize.nparm] = 0.5;
1328  /* range option in, but cannot be varied */
1329  optimize.vparm[1][optimize.nparm] = (realnum)rfield.range[nqh][0];
1330  optimize.vparm[2][optimize.nparm] = (realnum)rfield.range[nqh][1];
1332  optimize.nparm += 1;
1333  }
1334  ++nqh;
1335  }
1336 
1337  else if( strcmp(chKey4,"MAGN") == 0 )
1338  {
1339  /* parse the magnetic field command, routine in magnetic.c */
1340  ParseMagnet( chCard );
1341  }
1342 
1343  else if( strcmp(chKey4,"MAP ") == 0 )
1344  {
1345  /* do cooling space map for specified zones, in reads2 */
1346  ParseMap(chCard);
1347  }
1348 
1349  else if( strcmp(chKey4,"META") == 0 )
1350  {
1351  /* read depletion for metals, all elements heavier than He
1352  * in reads2 */
1353  ParseMetal(chCard);
1354  /* abundances no longer solar */
1355  abund.lgAbnSolar = false;
1356  }
1357 
1358  else if( strcmp(chKey4,"NEUT") == 0 )
1359  {
1360  /* heating and ionization due to fast neutrons */
1361  hextra.lgNeutrnHeatOn = true;
1362 
1363  /* first number is neutron luminosity
1364  * relative to bolometric luminosity */
1365  i = 5;
1366  hextra.frcneu = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
1367  if( lgEOL )
1368  NoNumb(chCard);
1369 
1370  /* save as log of fraction */
1371  if( hextra.frcneu > 0. )
1372  hextra.frcneu = (realnum)log10(hextra.frcneu);
1373 
1374  /* second number is efficiency */
1375  hextra.effneu = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
1376  if( lgEOL )
1377  {
1378  hextra.effneu = 1.0;
1379  }
1380  else
1381  {
1382  if( hextra.effneu <= 0. )
1383  hextra.effneu = (realnum)pow((realnum)10.f,hextra.effneu);
1384  }
1385  }
1386 
1387  else if( strcmp(chKey3,"NO ") == 0 )
1388  {
1389  /* don't do something, in readsun */
1390  ParseDont(chCard);
1391  }
1392 
1393  else if( strcmp(chKey4,"NORM") == 0 )
1394  {
1395  /* normalize lines to this rather than h-b, sec number is scale factor */
1396  ParseNorm(chCard);
1397  }
1398 
1399  else if( strcmp(chKey4,"NUF(") == 0 )
1400  {
1401  /* flux density of this continuum source, at optional frequency
1402  * expressed as product nu*f_nu */
1403  lgNu2 = true;
1404  /* in reads2 */
1405  ParseF_nu(chCard,&nqh,&ar1,"SQCM",lgNu2);
1406  }
1407 
1408  else if( strcmp(chKey4,"NUL(") == 0 )
1409  {
1410  /* specific luminosity density of this continuum source, at opt freq
1411  * expressed as product nu*f_nu */
1412  lgNu2 = true;
1413  /* in reads2 */
1414  ParseF_nu(chCard,&nqh,&ar1,"4 PI",lgNu2);
1415  }
1416 
1417  else if( strcmp(chKey4,"OPTI") == 0 )
1418  {
1419  /* option to optimize the model by varying certain parameters
1420  * in reads2 */
1421  ParseOptimize(chCard);
1422  }
1423 
1424  else if( strcmp(chKey4,"PHI(") == 0 )
1425  {
1426  /* enter phi(h), the number of h-ionizing photons/cm2 */
1427  i = 5;
1428  if( nqh >= LIMSPC )
1429  {
1430  /* too many continua were entered */
1431  fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" );
1432  cdEXIT(EXIT_FAILURE);
1433  }
1434  /* silly, but calms down the lint */
1435  ASSERT( nqh < LIMSPC );
1436  strcpy( rfield.chRSpec[nqh], "SQCM" );
1437  strcpy( rfield.chSpNorm[nqh], "PHI " );
1438  rfield.totpow[nqh] = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
1439  if( lgEOL )
1440  {
1441  NoNumb(chCard);
1442  }
1443  /* set radius to very large value if not already set */
1444  /* >>chng 01 jul 24, from Radius == 0 to this, as per PvH comments */
1445  if( !radius.lgRadiusKnown )
1446  {
1447  /* set radius to large value */
1448  ar1 = (realnum)radius.rdfalt;
1449  radius.Radius = pow(10.,radius.rdfalt);
1450  }
1451  /* check if continuum so intense that crash is likely */
1452  if( rfield.totpow[nqh] > 29. )
1453  {
1454  fprintf( ioQQQ, " Is the flux for this continuum correct?\n" );
1455  fprintf( ioQQQ, " It appears too bright to me.\n" );
1456  }
1457  /* ParseRangeOption in readsun xx parse_rangeoption*/
1458  ParseRangeOption(nqh,chCard);
1459 
1460  /* >>chng 06 mar 22, add time option to vary only some continua with time */
1461  if( nMatch( "TIME" , chCard ) )
1462  rfield.lgTimeVary[nqh] = true;
1463 
1464  /* vary option */
1465  if( optimize.lgVarOn )
1466  {
1467  /* >>chng 03 apr 30, had used log for range option, but if one or other
1468  * was 1 ryd, then became 0 in the log, which was misinterpreted as
1469  * one of the continuum bounds */
1470  strcpy( optimize.chVarFmt[optimize.nparm], "phi(h) %f range %f %f" );
1471  /* pointer to where to write */
1474  optimize.vincr[optimize.nparm] = 0.5;
1475  /* range option in, but cannot be varied */
1476  optimize.vparm[1][optimize.nparm] = (realnum)rfield.range[nqh][0];
1477  optimize.vparm[2][optimize.nparm] = (realnum)rfield.range[nqh][1];
1479  optimize.nparm += 1;
1480  }
1481  ++nqh;
1482  }
1483 
1484  else if( strcmp(chKey4,"PLOT") == 0 )
1485  {
1486  /* make plot of log(nu.f(nu)) vs log(nu) or opacity
1487  * in file plot */
1488  ParsePlot(chCard);
1489  }
1490 
1491  else if( strcmp(chKey4,"POWE") == 0 )
1492  {
1493  /* power law with cutoff, in reads2 */
1494  ParsePowerlawContinuum(chCard);
1495  }
1496 
1497  else if( strcmp(chKey4,"PRIN") == 0 )
1498  {
1499  /* adjust print schedule, in readsun */
1500  ParsePrint(chCard);
1501  }
1502 
1503  else if( strcmp(chKey4,"PUNC") == 0 )
1504  {
1505  /* punch something, in punch */
1506  ParsePunch(chCard);
1507  }
1508 
1509  else if( strcmp(chKey4,"Q(H)") == 0 )
1510  {
1511  /* log of number of ionizing photons */
1512  i = 5;
1513  if( nqh >= LIMSPC )
1514  {
1515  /* too many continua were entered */
1516  fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" );
1517  cdEXIT(EXIT_FAILURE);
1518  }
1519 
1520  /* silly, but calms down the lint */
1521  ASSERT( nqh < LIMSPC );
1522  strcpy( rfield.chRSpec[nqh], "4 PI" );
1523  strcpy( rfield.chSpNorm[nqh], "Q(H)" );
1524  rfield.totpow[nqh] = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
1525  if( rfield.totpow[nqh] > 100. && called.lgTalk )
1526  {
1527  fprintf( ioQQQ, " Is this reasonable?\n" );
1528  }
1529  if( lgEOL )
1530  {
1531  NoNumb(chCard);
1532  }
1533  /* ParseRangeOption in readsun */
1534  ParseRangeOption(nqh,chCard);
1535 
1536  /* >>chng 06 mar 22, add time option to vary only some continua with time */
1537  if( nMatch( "TIME" , chCard ) )
1538  rfield.lgTimeVary[nqh] = true;
1539 
1540  /* vary option */
1541  if( optimize.lgVarOn )
1542  {
1543  /* >>chng 03 apr 30, had used log for range option, but if one or other
1544  * was 1 ryd, then became 0 in the log, which was misinterpreted as
1545  * one of the continuum bounds */
1546  strcpy( optimize.chVarFmt[optimize.nparm], "Q(H) %f range %f %f" );
1547  /* pointer to where to write */
1550  optimize.vincr[optimize.nparm] = 0.5;
1551  /* range option in, but cannot be varied */
1552  optimize.vparm[1][optimize.nparm] = (realnum)rfield.range[nqh][0];
1553  optimize.vparm[2][optimize.nparm] = (realnum)rfield.range[nqh][1];
1555  ++optimize.nparm;
1556  }
1557  /* increment number of luminosity sources specified */
1558  ++nqh;
1559  }
1560 
1561  else if( strcmp(chKey4,"RATI") == 0 )
1562  {
1563  /* enter a continuum luminosity as a ratio of
1564  * nuFnu for this continuum relative to a previous continuum
1565  * format; first number is ratio of second to first continuum
1566  * second number is energy for this ratio
1567  * if third number on line, then 2nd number is energy of
1568  * first continuum, while 3rd number is energy of second continuum
1569  * in reads2 */
1570  ParseRatio(chCard,&nqh);
1571  }
1572 
1573  else if( strcmp(chKey4,"RADI") == 0 )
1574  {
1575  /* log of inner and outer radii, default second=infinity,
1576  * if R2<R1 then R2=R1+R2
1577  * there is an optional keyword, "PARSEC" on the line
1578  * to use PC as units, reads2 */
1579  ParseRadius(chCard,&ar1);
1580  }
1581 
1582  else if( strcmp(chKey4,"ROBE") == 0 )
1583  {
1584  /* this is the Roberto Terlivich command, no telling if it still works */
1585  radius.dRadSign = -1.;
1586  }
1587 
1588  else if( strcmp(chKey4,"SET ") == 0 )
1589  {
1590  /* set something, in reads2 */
1591  ParseSet(chCard);
1592  }
1593 
1594  else if( strcmp(chKey4,"SPEC") == 0 )
1595  {
1596  /* special key, can do anything */
1597  cdEXIT(EXIT_FAILURE);
1598  }
1599 
1600  else if( strcmp(chKey4,"SPHE") == 0 )
1601  {
1602  /* compute a spherical model, diffuse field from other side in
1603  * in reads2 */
1604  ParseSphere(chCard);
1605  }
1606 
1607  else if( strcmp(chKey4,"STAT") == 0 )
1608  {
1609  /* either get or put the code's state as a file on disk */
1610  ParseState(chCard);
1611  }
1612 
1613  else if( strcmp(chKey4,"STOP") == 0 )
1614  {
1615  /* stop model at desired zone, temperature, column density or tau-912
1616  * in readsun */
1617  ParseStop(chCard);
1618  }
1619 
1620  else if( strcmp(chKey4,"TABL") == 0 )
1621  {
1622  /* interpolate on input tables for continuum, set of power laws used
1623  * input stored in big BLOCK data
1624  * first check that this is the one and only INTERP command
1625  * in readsun */
1626  ParseTable(&nqh,chCard,&ar1);
1627  }
1628 
1629  else if( strcmp(chKey4,"TAUM") == 0 )
1630  {
1631  /* minimum optical depths for lines (normally 1e-20)
1632  * (used in STARTER to set lines up) */
1633  i = 5;
1634  opac.taumin = (realnum)pow(10.,FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL));
1635  if( lgEOL )
1636  {
1637  NoNumb(chCard);
1638  }
1639  }
1640 
1641  else if( strcmp(chKey4 , "TEST" ) == 0 )
1642  {
1643  /* parse the test command and its options */
1644  ParseTest(chCard,&nqh , &ar1, lgDSet);
1645  }
1646 
1647  else if( strcmp(chKey4,"TIME") == 0 )
1648  {
1649  /* parse the time dependent command, in dynamics.c */
1650  ParseDynaTime(chCard);
1651  }
1652 
1653  else if( strcmp(chKey4,"TITL") == 0 )
1654  {
1655  /* read in title of model starting in col 5 */
1656  strcpy( input.chTitle , &input.chOrgCard[5] );
1657  }
1658 
1659  else if( strcmp(chKey4,"TLAW") == 0 )
1660  {
1661  /* some temperature vs depth law */
1662  ParseTLaw(chCard);
1663  }
1664 
1665  else if( strcmp(chKey4,"TOLE") == 0 )
1666  {
1667  fprintf(ioQQQ,
1668  "Sorry, this command has been replaced with the SET TEMPERATURE TOLERANCE command.\n");
1669  cdEXIT(EXIT_FAILURE);
1670  }
1671 
1672  else if( strcmp(chKey4,"TRAC") == 0 )
1673  {
1674  /* turn on trace output, in reads2 */
1675  ParseTrace(chCard);
1676  }
1677 
1678  else if( strcmp(chKey4,"VLAW") == 0 )
1679  {
1680 
1681  /* velocity power law as a function of radius. */
1682  i = 5;
1683  DoppVel.TurbVelLaw = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
1684 
1685  DoppVel.lgTurbLawOn = true;
1687  /* for now, insist upon non-positive power,
1688  * so that velocity decreases with increasing radius. */
1689  ASSERT( DoppVel.TurbVelLaw <= 0.f );
1690  }
1691 
1692  else if( strcmp(chKey4,"TURB") == 0 )
1693  {
1694  /* line has turbulent velocity in km/s */
1695  i = 5;
1696  DoppVel.TurbVel = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
1697 
1698  /* include turbulent pressure in equation of state?
1699  * >>chng 06 mar 24, include turbulent pressure in gas equation of state by default,
1700  * but not if NO PRESSURE occurs */
1701  if( nMatch(" NO ",chCard) && nMatch("PRES",chCard) )
1702  {
1703  DoppVel.lgTurb_pressure = false;
1704  }
1705  else
1706  {
1707  /* default is to include turbulent pressure in equation of state */
1708  DoppVel.lgTurb_pressure = true;
1709  }
1710 
1711  /* true if no number on line - check for equipartition */
1712  DoppVel.lgTurbEquiMag = false;
1713 
1714  if( nMatch("EQUI",chCard) && nMatch("PART",chCard) )
1715  {
1716  /* turbulence equipartition option */
1717  DoppVel.lgTurbEquiMag = true;
1718  }
1719 
1720  if( nMatch(" LOG",chCard) )
1721  {
1722  DoppVel.TurbVel = (realnum)pow((realnum)10.f,DoppVel.TurbVel);
1723  }
1724  /* now convert from km/s to cm/s */
1725  DoppVel.TurbVel *= 1e5;
1726 
1727  /* this is optional F parameter from equation 34 of
1728  *>>refer pressure turb Heiles, C. & Troland, T.H. 2005, ApJ, 624, 773
1729  * turbulent energy density will be rho F v^2 / 2 */
1731  if( lgEOL )
1732  {
1733  /* this is the default value of 3 for isotropic turbulence */
1734  DoppVel.Heiles_Troland_F = 3.f;
1735  }
1736 
1737  /* option to have turbulence be dissipative, keyword is dissipate,
1738  * argument is log of scale length in cm */
1739  if( nMatch("DISS",chCard) )
1740  {
1741  DoppVel.DispScale = (realnum)pow(10., FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL) );
1742  if( lgEOL )
1743  {
1744  NoNumb(chCard);
1745  }
1746  }
1747 
1748  /* vary option */
1749  if( optimize.lgVarOn )
1750  {
1751  /* only one parameter to vary */
1753  strcpy( optimize.chVarFmt[optimize.nparm], "TURBULENCE= %f " );
1754  /* pointer to where to write */
1756  /* turbulent velocity */
1759  ++optimize.nparm;
1760  }
1761 
1762  /* set initial turbulence */
1764  }
1765 
1766  else if( strcmp(chKey4,"WIND") == 0 )
1767  {
1768  /* NB - advection and wind commands are now a single command */
1769  /* parse the wind command, in dynamics.c */
1770  ParseDynaWind(chCard);
1771  }
1772 
1773  else if( strcmp(chKey2,"XI") == 0 )
1774  {
1775  /* inter x-ray ionization parameter xi=L/r^2 n;
1776  * defined per hydrogen
1777  * radius must also be entered if spherical, but not needed if plane */
1778  ParseIonPar(&nqh,chCard,&ar1);
1779  }
1780 
1781  /* a line starting with these characters is a comment and so ok
1782  * >>chng 06 sep 04 use routine to check for comments */
1783  else if( !lgInputComment(chCard) /*(chK1 != '#' && chK1 != '*') && chK1 != '%'*/ )
1784  {
1785  /* end of keys - if we get here key is unrecognized---------- */
1786  fprintf( ioQQQ, " Unrecognized command. Key=\"%4.4s\". This is routine ParseCommands.\n",
1787  chKey4 );
1788  fprintf( ioQQQ, " The line image was==%s==\n",
1789  chCard );
1790  fprintf( ioQQQ, " Sorry.\n" );
1791  cdEXIT(EXIT_FAILURE);
1792  }
1793 
1794  /* get next line image, this is tail of while loop that will
1795  * look for lgEOF or blank card */
1796  input_readarray(chCard,&lgEOF);
1797  }
1798 
1799  /*------------------------------------------------------------------- */
1800  /* fall through - hit lgEOF or blank line */
1801 
1802  /* >>chng 00 mar 31, removed dummy call to ParseQuantHeat, PvH */
1803 
1804  for( i=0; i < INPUT_LINE_LENGTH; i++ )
1805  {
1806  chCard[i] = ' ';
1807  }
1808  chCard[INPUT_LINE_LENGTH-1] = '\0';
1809 
1810  if( called.lgTalk )
1811  {
1812  fprintf( ioQQQ, "%23c*%81c*\n", ' ', ' ' );
1813  fprintf( ioQQQ, "%23c***********************************************************************************\n\n\n\n", ' ' );
1814  }
1815 
1816  /* this is an option to turn on trace printout on the nth
1817  * call from the optimizer - only allow trace if
1818  * this is the case and nOptimiz 1 below nTrOpt */
1819  if( optimize.lgTrOpt )
1820  {
1821  /* nTrOpt was set with "optimize trace" command,
1822  * is iteration to turn on trace */
1823  if( optimize.nTrOpt != optimize.nOptimiz )
1824  {
1825  trace.lgTrace = false;
1826  /* following overrides turning on trace elsewhere */
1827  trace.lgTrOvrd = false;
1828  }
1829  else
1830  {
1831  trace.lgTrace = true;
1832  called.lgTalk = true;
1833  trace.lgTrOvrd = true;
1834  fprintf( ioQQQ, " READR turns on trace from optimize option.\n" );
1835  }
1836  }
1837 
1838  /* set density from DLAW command, must be done here since it may depend on later input commands */
1839  if( strcmp(dense.chDenseLaw,"DLW1") == 0 )
1840  {
1842 
1843  if( dense.gas_phase[ipHYDROGEN] <= 0. )
1844  {
1845  fprintf( ioQQQ, " PROBLEM DISASTER Hydrogen density set by DLAW must be > 0.\n" );
1846  cdEXIT(EXIT_FAILURE);
1847  }
1848  }
1849  else if( strcmp(dense.chDenseLaw,"DLW2") == 0 )
1850  {
1852 
1853  if( dense.gas_phase[ipHYDROGEN] <= 0. )
1854  {
1855  fprintf( ioQQQ, " PROBLEM DISASTER Hydrogen density set by DLAW must be > 0.\n" );
1856  cdEXIT(EXIT_FAILURE);
1857  }
1858  }
1859 
1860  /* start checks on parameters set properly - this begins with same line saying start .. */
1861 
1862  /* lgStop_not_enough_info says that not enough info for model, so stop
1863  * set true in following tests if anything missing */
1864  lgStop_not_enough_info = false;
1865  lgStop = false;
1866 
1867  /* check whether hydrogen density has been set - this value was set to 0 in zero */
1868  if( dense.gas_phase[ipHYDROGEN] <= 0. )
1869  {
1870  fprintf( ioQQQ, " PROBLEM DISASTER Hydrogen density MUST be specified.\n" );
1871  lgStop_not_enough_info = true;
1872  lgStop = true;
1873  /* need to set something since used below - will abort
1874  * since lgStop is set */
1875  dense.gas_phase[ipHYDROGEN] = 1.;
1876  }
1877 
1879 
1880  /* mass flux for wind model - used for mass conservation */
1882 
1883  /* set converge criteria - limit number of iterations and zones */
1885  {
1887  for( j=0; j < iterations.iter_malloc; j++ )
1888  {
1890  }
1891  }
1892 
1893  /* NSAVE is number of lines saved, =0 if no commands entered */
1894  if( input.nSave < 0 )
1895  {
1896  fprintf( ioQQQ, " PROBLEM DISASTER No commands were entered - whats up?\n" );
1897  cdEXIT(EXIT_FAILURE);
1898  }
1899 
1900  /* iterate to convergence and wind models are mutually exclusive */
1901  if( wind.windv0 > 0. && conv.lgAutoIt )
1902  {
1903  if( called.lgTalk )
1904  {
1905  fprintf( ioQQQ, " NOTE PROBLEM Due to the nature of the Sobolev approximation, it makes no sense to converge a windy model.\n" );
1906  fprintf( ioQQQ, " NOTE Iterate to convergence is turned off\n\n\n" );
1907  }
1908  conv.lgAutoIt = false;
1909  iterations.itermx = 0;
1910  }
1911 
1912  /* iterate to convergence and case b do not make sense together */
1913  /* WJH 22 May 2004: unless we are doing i-front dynamics (negative wind.windv0) */
1914  if( opac.lgCaseB && conv.lgAutoIt && wind.windv0 >= 0. )
1915  {
1916  if( called.lgTalk )
1917  {
1918  fprintf( ioQQQ, " NOTE Case B is an artificial test, it makes no sense to converge this model.\n" );
1919  fprintf( ioQQQ, " NOTE Iterate to convergence is turned off.\n\n\n" );
1920  }
1921  conv.lgAutoIt = false;
1922  iterations.itermx = 0;
1923  }
1924 
1925  /* specifying a density power and constant pressure makes no sense */
1926  if( dense.DensityPower!=0. && strcmp( dense.chDenseLaw, "CPRE" )==0 )
1927  {
1928  if( called.lgTalk )
1929  {
1930  fprintf( ioQQQ, " NOTE Specifying both a density power law and constant pressure is impossible.\n" );
1931  }
1932  lgStop = true;
1933  }
1934 
1935  if( !rfield.lgIonizReevaluate && strcmp( dense.chDenseLaw, "CDEN" )!=0 )
1936  {
1937  if( called.lgTalk )
1938  {
1939  fprintf( ioQQQ, " NOTE NO REEVALUATE IONIZATION can only be used with constant density.\n" );
1940  fprintf( ioQQQ, " NOTE Resetting to reevaluate ionization.\n\n" );
1941  }
1942  rfield.lgIonizReevaluate = true;
1943  }
1944 
1945  /* check if the combination of stopping column density and H density are physically plausible */
1946  thickness = min( MIN3( StopCalc.tauend, StopCalc.colpls, StopCalc.colnut ),
1948  if( thickness < COLUMN_INIT )
1949  {
1950  /* a stop column density was specified - check on physical thickness this corresponds to */
1951  thickness /= (dense.gas_phase[ipHYDROGEN]*geometry.FillFac);
1952  /* >> chng 06 dec 21, don't complain if outer radius set small with `stop thickness' command. */
1953  if( thickness > 1e25 && radius.router[0] > 1e25 )
1954  {
1955  fprintf( ioQQQ,
1956  "NOTE The specified column density and hydrogen density correspond to a thickness of %.2e cm.\n",
1957  thickness);
1958  fprintf( ioQQQ,
1959  "NOTE This seems large to me.\n");
1960  fprintf(ioQQQ,"NOTE a very large radius may cause overflow.\n\n");
1961  }
1962  }
1963 
1965  {
1966  /* warn if constant grain temperature but gas-grain thermal effects
1967  * are still included */
1968  fprintf( ioQQQ,
1969  "NOTE The grain temperatures are set to a constant value with the "
1970  "CONSTANT GRAIN TEMPERATURE command, but "
1971  "energy exchange \n");
1972  fprintf( ioQQQ,
1973  "NOTE is still included. The grain-gas heating-cooling will be incorrect. "
1974  "Consider turning off gas-grain collisional energy\n");
1975  fprintf( ioQQQ,
1976  "NOTE exchange with the NO GRAIN GAS COLLISIONAL ENERGY EXCHANGE command.\n\n\n");
1977  }
1978 
1979  /* >>chng 04 nov 27, test on these two */
1981  {
1982  if( called.lgTalk )
1983  {
1984  fprintf( ioQQQ, " NOTE NO LINE TRANSER set but fine opacities still computed.\n" );
1985  fprintf( ioQQQ, " NOTE Turning off fine opacities.\n\n" );
1986  }
1987  rfield.lgOpacityFine = false;
1988  }
1989 
1990  /* >>chng 04 nov 27, test of these two */
1992  {
1993  if( called.lgTalk )
1994  {
1995  fprintf( ioQQQ, " NOTE Large H2 molecule turned on but line transfer and fine opacities are not.\n" );
1996  fprintf( ioQQQ, " NOTE Turning on line transfer and fine opacities.\n\n" );
1997  }
1998  rfield.lgOpacityFine = true;
1999  rfield.lgDoLineTrans = true;
2000  }
2001 
2003  {
2004  /* one of the input continua needs to have H-ionizing radiation
2005  * blocked with extinguish command, but it was not done */
2006  fprintf( ioQQQ, "\n NOTE\n"
2007  " NOTE One of the incident continuum is a form used when no H-ionizing radiation is produced.\n" );
2008  fprintf( ioQQQ, " NOTE You must also include the EXTINGUISH command to make sure this is done.\n" );
2009  fprintf( ioQQQ, " NOTE The EXTINGUISH command was not included.\n" );
2010  fprintf( ioQQQ, " NOTE YOU MAY BE MAKING A BIG MISTAKE!!\n NOTE\n\n\n\n" );
2011  }
2012 
2013  /* if stop temp set below default then we are going into cold and possibly molecular
2014  * gas - check some parameters in this case */
2016  /* thermal.ConstTemp def is zero, set pos when constant temperature is set */
2018  {
2019  /* print warning if temperature set below default but cosmic rays not turned on */
2020  /* >>chng 06 mar 02, do not print if molecules are off */
2021  if( (hextra.cryden == 0.) && !co.lgNoCOMole && !hmi.lgNoH2Mole )
2022  {
2023  fprintf( ioQQQ, "\n NOTE\n"
2024  " NOTE The simulation is going into neutral gas but cosmic rays are not included.\n" );
2025  fprintf( ioQQQ, " NOTE Ion-molecule chemistry will not occur without a source of ionization.\n" );
2026  fprintf( ioQQQ, " NOTE The chemistry network may collapse deep in molecular regions.\n" );
2027  fprintf( ioQQQ, " NOTE Consider adding galactic background cosmic rays with the COSMIC RAYS BACKGROUND command.\n" );
2028  fprintf( ioQQQ, " NOTE You may be making a BIG mistake.\n NOTE\n\n\n\n" );
2029  }
2030  }
2031 
2032  /* dense.gas_phase[ipHYDROGEN] is linear hydrogen density (cm-3) */
2033  /* test for hydrogen density properly set has already been done above */
2034  if( called.lgTalk && dense.gas_phase[ipHYDROGEN] < 1e-4 )
2035  {
2036  fprintf( ioQQQ, " NOTE Is the entered value of the hydrogen density (%.2e) reasonable?\n",
2038  fprintf( ioQQQ, " NOTE It seems pretty low to me.\n\n\n" );
2039  }
2040  else if( called.lgTalk && dense.gas_phase[ipHYDROGEN] > 1e15 )
2041  {
2042  fprintf( ioQQQ, " NOTE Is this value of the hydrogen density reasonable?\n" );
2043  fprintf( ioQQQ, " NOTE It seems pretty high to me.\n\n\n" );
2044  }
2045 
2046  /* is the model going to crash because of extreme density? */
2047  if( called.lgTalk && !lgStop && !lgStop_not_enough_info )
2048  {
2049  if( dense.gas_phase[ipHYDROGEN] < 1e-6 || dense.gas_phase[ipHYDROGEN] > 1e19 )
2050  {
2051  fprintf( ioQQQ, " NOTE Simulation may crash because of extreme "
2052  "density. The value was %.2e\n\n" ,
2054  }
2055  }
2056 
2057  if( rfield.nspec == 0 )
2058  {
2059  fprintf( ioQQQ, " PROBLEM DISASTER No incident radiation field was specified - "
2060  "at least put in the CMB.\n" );
2061  lgStop = true;
2062  lgStop_not_enough_info = true;
2063  }
2064 
2065  if( nqh == 0 )
2066  {
2067  fprintf( ioQQQ, " PROBLEM DISASTER Luminosity of continuum MUST be specified.\n" );
2068  lgStop = true;
2069  lgStop_not_enough_info = true;
2070  }
2071 
2072  /* Testing on 0 is safe since this it is initialized that way
2073  * only print comment if continuum has been specified,
2074  * if continuum not given then we are aborting anyway */
2075  if( radius.Radius == 0. && rfield.nspec> 0)
2076  {
2077  fprintf( ioQQQ, " PROBLEM DISASTER Starting radius MUST be specified.\n" );
2078  lgStop = true;
2079  lgStop_not_enough_info = true;
2080  }
2081 
2082  if( rfield.nspec != nqh )
2083  {
2084  fprintf( ioQQQ, " PROBLEM DISASTER There were not the same number of continuum shapes and luminosities entered.\n" );
2085  lgStop = true;
2086  }
2087 
2088  /* we only want to do this test on the first call to the command
2089  * parser - it will be called many more times but with no grid command
2090  * during the actual grid calculation */
2091  static bool lgFirstPass = true;
2092 
2093  /* the NO VARY option sets this flag, and can be used to turn off
2094  * the grid command, as well as the optimizer */
2095  if( optimize.lgNoVary && grid.lgGrid )
2096  {
2097  /* ignore grids */
2098  grid.lgGrid = false;
2099  optimize.nparm = 0;
2100  grid.nGridCommands = 0;
2101  }
2102 
2103  if( lgFirstPass && grid.lgGrid && (optimize.nparm!=grid.nGridCommands) )
2104  {
2105  /* number of grid vary options do match */
2106  fprintf( ioQQQ, " PROBLEM DISASTER The GRID command was entered "
2107  "but there were %li GRID commands and %li commands with a VARY option.\n" ,
2109  fprintf( ioQQQ, " There must be the same number of GRIDs and VARY.\n" );
2110  lgStop = true;
2111  }
2112  lgFirstPass = false;
2113 
2114  if( lgStop_not_enough_info )
2115  {
2116  fprintf( ioQQQ, " PROBLEM DISASTER I do not have enough information to do the simulation, I cannot go on.\n" );
2117  fprintf( ioQQQ, "\n\n Sorry.\n\n\n" );
2118  cdEXIT(EXIT_FAILURE);
2119  }
2120 
2121  if( lgStop )
2122  {
2123  cdEXIT(EXIT_FAILURE);
2124  }
2125 
2126  /* end checks on parameters set properly - this begins with same line saying start .. */
2127  return;
2128 }

Generated for cloudy by doxygen 1.8.1.2