cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
parse_stop.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 /*ParseStop parse the stop command */
4 #include "cddefines.h"
5 #include "physconst.h"
6 #include "optimize.h"
7 #include "phycon.h"
8 #include "rfield.h"
9 #include "radius.h"
10 #include "geometry.h"
11 #include "iterations.h"
12 #include "stopcalc.h"
13 #include "input.h"
14 #include "parse.h"
15 
16 void ParseStop(char *chCard)
17 {
18  bool lgEOL;
19 
20  long int i,
21  j;
22 
23  double a,
24  effcol,
25  tread;
26 
27  DEBUG_ENTRY( "ParseStop()" );
28 
29  /* first scan off a generic number, all subcommands use at least one */
30  i = 5;
31  a = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
32 
33  /* >>chng 04 nov 09, introduce off option */
34  if( lgEOL && !nMatch(" OFF",chCard) )
35  {
36  NoNumb(chCard);
37  }
38 
39  if( nMatch("TEMP",chCard) )
40  {
41  /* >>chng 04 nov 09, introduce off option to disable this stopping criterion */
42  if( lgEOL && nMatch(" OFF",chCard) )
43  {
44  /* this is special case for ending temperature - do not use -
45  * but will still stop if Te falls below TeLowest, the lowest
46  * possible temperature */
47  StopCalc.tend = -1.f;
48  }
49  else
50  {
51  /* lowest electron temperature allowed before stopping
52  * assumed to be the log of the temperature if <10
53  * optional keyword LINEAR forces linear */
54  if( a <= 10. && !nMatch("LINE",chCard) )
55  {
56  tread = pow(10.,a);
57  }
58  else
59  {
60  tread = a;
61  }
62 
63  /* tread is linear temperature*/
64  if( tread < phycon.TEMP_LIMIT_LOW )
65  {
66  fprintf( ioQQQ,
67  " Temperatures below %.2e K not allowed. Reset to lowest value."
68  " I am doing this myself.\n" ,
70  /* set slightly off extreme limit for safety */
71  tread = phycon.TEMP_LIMIT_LOW*1.01;
72  }
73  else if( tread > phycon.TEMP_LIMIT_HIGH )
74  {
75  fprintf( ioQQQ,
76  " Temperatures is above %.2e K not allowed. Reset to highest value."
77  " I am doing this myself.\n" ,
79  /* set slightly off extreme limit for safety */
80  tread = phycon.TEMP_LIMIT_HIGH*0.99;
81  }
82 
83  if( nMatch("EXCE",chCard) )
84  {
85  /* option for this to be highest allowed temperature,
86  * stop temperate exceeds */
87  StopCalc.T2High = (realnum)tread;
88  }
89  else
90  {
91  /* this is ending temperature - we stop if kinetic temperature
92  * falls below this */
93  StopCalc.tend = (realnum)tread;
94  }
95  }
96  }
97 
98  /* stop at 21cm line center optical depth */
99  else if( nMatch("OPTI",chCard) && nMatch("21CM",chCard) )
100  {
101  /* default is for number to be log of optical depth */
102  bool lgLOG = true;
103  if( nMatch("LINE",chCard) )
104  {
105  /* force it to be linear not log */
106  lgLOG = false;
107  }
108  i = 4;
109  j = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
110  if( j!=21 )
111  {
112  fprintf( ioQQQ, " First number on STOP 21CM OPTICAL DEPTH command must be 21\n" );
113  cdEXIT(EXIT_FAILURE);
114  }
115  /* now get the next number, which is the optical depth */
116  a = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
117 
118  /* tau must be a log, and second number is energy where tau specified */
119  if( lgLOG )
120  {
121  StopCalc.tauend = (realnum)pow(10.,a);
122  }
123  else
124  {
125  StopCalc.tauend = (realnum)a;
126  }
127  /* this flag says that 21cm line optical depth is the stop quantity */
128  StopCalc.lgStop21cm = true;
129  }
130  /* stop optical depth at some energy */
131  else if( nMatch("OPTI",chCard) )
132  {
133  /* default is for number to be log of optical depth */
134  bool lgLOG = true;
135  if( nMatch("LINE",chCard) )
136  {
137  /* force it to be linear not log */
138  lgLOG = false;
139  }
140 
141  if( a > 37. )
142  {
143  fprintf( ioQQQ, " optical depth too big\n" );
144  cdEXIT(EXIT_FAILURE);
145  }
146 
147  /* tau must be a log, and second number is energy where tau specified */
148  if( lgLOG )
149  {
150  StopCalc.tauend = (realnum)pow(10.,a);
151  }
152  else
153  {
154  StopCalc.tauend = (realnum)a;
155  }
156 
157  /* energy where tau specified */
158  StopCalc.taunu = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
159 
160  if( lgEOL )
161  {
162  if( nMatch("LYMA",chCard) )
163  {
164  /* largest Lyman limit optical depth */
165  StopCalc.taunu = 1.;
166  }
167  else if( nMatch("BALM",chCard) )
168  {
169  /* stop at this Balmer continuum optical depth */
170  StopCalc.taunu = 0.25;
171  }
172  else
173  {
174  fprintf( ioQQQ, " There must be a second number, the energy in Ryd. Sorry.\n" );
175  cdEXIT(EXIT_FAILURE);
176  }
177  }
178 
179  else
180  {
181  /* if second number is negative then log of energy in rydbergs */
182  if( StopCalc.taunu < 0. )
183  {
184  StopCalc.taunu = (realnum)pow((realnum)10.f,StopCalc.taunu);
185  }
186 
187  /* check that energy is within bounds of code */
189  {
190  fprintf( ioQQQ, " The energy must be in the range %10.2e to %10.2e. It was %10.2e. Sorry.\n",
192  cdEXIT(EXIT_FAILURE);
193  }
194  }
195  }
196 
197  /* stop optical depth at extinction in V filter */
198  else if( nMatch(" AV ",chCard) )
199  {
200  /* default is for number to be A_V, log if negative */
201  if( a<=0. )
202  {
203  a = pow(10.,a);
204  }
205  /* A_V can be for either point or extended source, difference is (1-g) multiplied by scat opacity
206  * if keyword point occurs then for point source, otherwise for extended source */
207  if( nMatch("EXTE" , chCard ) )
208  {
210  }
211  else
212  {
213  /* default is point, as measured in ism work */
215  }
216  }
217 
218  /* stop when a fraction of molecules frozen out on grain surfaces is reached */
219  else if( nMatch("MOLE",chCard) && nMatch("DEPL",chCard) )
220  {
221  if( a <= 0. )
222  {
223  StopCalc.StopDepleteFrac = (realnum)pow(10.,a);
224  }
225  else
226  {
228  }
229  }
230 
231  /* stop when absolute value of flow velocity falls below this value */
232  else if( nMatch("VELO",chCard) )
233  {
234  /* entered in km/s but stored as cm/s */
235  StopCalc.StopVelocity = (realnum)(a*1e5);
236  }
237 
238  /* stop at a given computed mass */
239  else if( nMatch("MASS",chCard) )
240  {
241  /* number of log of mass in gm if inner radius is specified,
242  * mass per unit area, gm cm-2 if not
243  * leave it as a log since dare not deal with linear mass */
244  StopCalc.xMass = (realnum)a;
245  /* NB 0 is sentinel for not set, if a is zero we must reset it */
246  if( StopCalc.xMass == 0 )
248  }
249 
250  /* stop thickness command, also stop depth, this must come after stop
251  * optical depth, since do not want to trigger on depth in optical depth */
252  else if( nMatch("THIC",chCard) || nMatch("DEPT",chCard) )
253  {
254  /* decide whether linear or log,
255  * this branch if key linear appears */
256  if( nMatch("LINE",chCard) )
257  {
258  radius.router[0] = a;
259  }
260  else
261  {
262  /* this is default branch, no key, so log */
263  if( a > 37. )
264  {
265  fprintf( ioQQQ, " thickness too large\n" );
266  cdEXIT(EXIT_FAILURE);
267  }
268  radius.router[0] = pow(10.,a);
269  }
270 
271  /* optional parsec unit */
272  if( nMatch("PARS",chCard) )
273  {
274  radius.router[0] *= PARSEC;
275  }
276 
277  /* can stop at different thickness on each iteration */
278  for( j=1; j < iterations.iter_malloc; j++ )
279  {
280  a = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
281  if( lgEOL )
282  {
283  radius.router[j] = radius.router[j-1];
284  }
285  else
286  {
287  if( nMatch("LINE",chCard) )
288  {
289  radius.router[j] = a;
290  }
291  else
292  {
293  if( a > 37. )
294  {
295  fprintf( ioQQQ, " thickness too large\n" );
296  cdEXIT(EXIT_FAILURE);
297  }
298  radius.router[j] = pow(10.,a);
299  }
300  if( nMatch("PARS",chCard) )
301  {
302  radius.router[j] *= PARSEC;
303  }
304  }
305  }
306 
307  /* vary option */
308  if( optimize.lgVarOn )
309  {
311  strcpy( optimize.chVarFmt[optimize.nparm], "STOP THICKNESS %f" );
312 
313  /* pointer to where to write */
315 
316  /* log of temp will be pointer */
317  optimize.vparm[0][optimize.nparm] = (realnum)log10(radius.router[0]);
318  optimize.vincr[optimize.nparm] = 0.5;
319 
320  ++optimize.nparm;
321  }
322  }
323 
324  /* stop at a particular zone, for each iteration */
325  else if( nMatch("ZONE",chCard) )
326  {
327  /* stop after computing this zone */
328  /* >>chng 03 jun 06, do not let fall below 1, stop zone 0 has same effect
329  * as stop zone 1, bug caught by Joop Schaye */
330  geometry.nend[0] = (long)MAX2(1.,a);
331  geometry.lgZoneSet = true;
332 
333  /* this tells code that we intend to stop at this zone, so caution not generated*/
334  geometry.lgEndDflt = false;
335 
336  for( j=1; j < iterations.iter_malloc; j++ )
337  {
338  geometry.nend[j] = (long)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
339  /* if eol on this iteration, set to previous. In most cases
340  * all will be equal to the first */
341  if( lgEOL )
342  {
343  geometry.nend[j] = geometry.nend[j-1];
344  }
345  else
346  {
347  /* >>chng 03 jun 06, do not let fall below 1, stop zone 0 has same effect
348  * as stop zone 1, bug caught by Joop Schaye */
349  geometry.nend[j] = MAX2( 1 , geometry.nend[j] );
350  }
351  }
352  }
353 
354  /* stop at this electron fraction, relative to hydrogen */
355  else if( nMatch("EFRA",chCard) )
356  {
357  if( a <= 0. )
358  {
359  StopCalc.StopElecFrac = (realnum)pow(10.,a);
360  }
361  else
362  {
364  }
365  }
366 
367  /* stop at a hydrogen molecular fraction, relative to total hydrogen,
368  * this is 2H_2 / H_total*/
369  else if( nMatch("MFRA",chCard) )
370  {
371  if( a <= 0. )
372  {
373  StopCalc.StopH2MoleFrac = (realnum)pow(10.,a);
374  }
375  else
376  {
378  }
379  }
380 
381  /* stop at a ionized hydrogen fraction, relative to total hydrogen,
382  * this is H+ / H_total */
383  else if( nMatch("PFRA",chCard) )
384  {
385  if( a <= 0. )
386  {
387  StopCalc.StopHPlusFrac = (realnum)pow(10.,a);
388  }
389  else
390  {
392  }
393  }
394 
395  /* stop at a particular column density */
396  else if( nMatch("COLU",chCard) )
397  {
398  /* check for linear option, if present take log since a being
399  * log column density is default */
400  if( nMatch( "LINE" , chCard ) )
401  a = log10(a);
402 
403  /* stop at an effective column density */
404  if( nMatch("EFFE",chCard) )
405  {
406  /* actually stop at certain optical depth at 1keV */
407  effcol = pow(10.,a);
408  StopCalc.tauend = (realnum)(effcol*2.14e-22);
409  StopCalc.taunu = 73.5;
410  /* vary option */
411  if( optimize.lgVarOn )
412  {
414  strcpy( optimize.chVarFmt[optimize.nparm], "STOP EFFECTIVE COLUMN DENSITY %f" );
415  /* pointer to where to write */
417  /* log of temp will be pointer */
418  optimize.vparm[0][optimize.nparm] = (realnum)log10(effcol);
419  optimize.vincr[optimize.nparm] = 0.5;
420  ++optimize.nparm;
421  }
422  }
423 
424  else if( nMatch("IONI",chCard) )
425  {
426  /* this is ionized column */
427  if( a > 37. )
428  {
429  fprintf( ioQQQ, " column too big\n" );
430  cdEXIT(EXIT_FAILURE);
431  }
432 
433  StopCalc.colpls = (realnum)pow(10.,a);
434 
435  /* vary option */
436  if( optimize.lgVarOn )
437  {
439  strcpy( optimize.chVarFmt[optimize.nparm], "STOP IONIZED COLUMN DENSITY %f" );
440  /* pointer to where to write */
442  /* log of temp will be pointer */
444  optimize.vincr[optimize.nparm] = 0.5;
445  ++optimize.nparm;
446  }
447  }
448 
449  /* stop at a neutral column */
450  else if( nMatch("NEUT",chCard) )
451  {
452  StopCalc.colnut = (realnum)pow(10.,a);
453 
454  /* vary option */
455  if( optimize.lgVarOn )
456  {
458  strcpy( optimize.chVarFmt[optimize.nparm], "STOP NEUTRAL COLUMN DENSITY %f");
459  /* pointer to where to write */
461  /* log of temp will be pointer */
463  optimize.vincr[optimize.nparm] = 0.5;
464  ++optimize.nparm;
465  }
466  }
467 
468  /* >>chng 03 apr 15, add this option
469  * stop at a molecular hydrogen column density, input parameter
470  * is log of the H2 column density */
471  else if( nMatch(" H2 ",chCard) )
472  {
473  /* this command has a 2 in the H2 label - must not parse the two by
474  * accident. Get the first number off the line image, and confirm that
475  * it is a 2 */
476  i = 5;
477  j = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
478  if( j != 2 )
479  {
480  fprintf( ioQQQ, " Something is wrong with the order of the numbers on this line.\n" );
481  fprintf( ioQQQ, " The first number I encounter should be the 2 in H2.\n Sorry.\n" );
482  cdEXIT(EXIT_FAILURE);
483  }
484  a = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
485  StopCalc.col_h2 = (realnum)pow(10.,a);
486 
487  /* vary option */
488  if( optimize.lgVarOn )
489  {
491  strcpy( optimize.chVarFmt[optimize.nparm], "STOP H2 COLUMN DENSITY %f");
492  /* pointer to where to write */
494  /* log of temp will be pointer */
496  optimize.vincr[optimize.nparm] = 0.5;
497  ++optimize.nparm;
498  }
499  }
500 
501  else if( nMatch("ATOM",chCard) )
502  {
503  StopCalc.col_h2_nut = (realnum)pow(10.,a);
504  /* vary option */
505  if( optimize.lgVarOn )
506  {
508  strcpy( optimize.chVarFmt[optimize.nparm], "STOP ATOMIC COLUMN DENSITY %f");
509  /* pointer to where to write */
511  /* log of temp will be pointer */
513  optimize.vincr[optimize.nparm] = 0.5;
514  ++optimize.nparm;
515  }
516  }
517 
518  else if( nMatch("H/TS",chCard) )
519  {
520  /* >> 05 jan 09, add stop integrated n(H0) / Tspin */
521  StopCalc.col_H0_ov_Tspin = (realnum)pow(10.,a);
522  /* vary option */
523  if( optimize.lgVarOn )
524  {
526  strcpy( optimize.chVarFmt[optimize.nparm], "STOP H/TSPIN COLUMN DENSITY %f");
527  /* pointer to where to write */
529  /* log of temp will be pointer */
531  optimize.vincr[optimize.nparm] = 0.5;
532  ++optimize.nparm;
533  }
534  }
535 
536  else if( nMatch(" CO ",chCard) )
537  {
538  /* chng << 03 Oct. 27--Nick Abel, add this option */
539  /* stop at a carbon monoxide column density */
540  StopCalc.col_monoxco = (realnum)pow(10.,a);
541  /* vary option */
542  if( optimize.lgVarOn )
543  {
545  strcpy( optimize.chVarFmt[optimize.nparm], "STOP CO COLUMN DENSITY %f");
546  /* pointer to where to write */
548  /* log of temp will be pointer */
550  optimize.vincr[optimize.nparm] = 0.5;
551  ++optimize.nparm;
552  }
553  }
554 
555  /* fall through default is total hydrogen column density */
556  else
557  {
558  /* both HII and HI */
559  if( a > 37. )
560  {
561  fprintf( ioQQQ, " column too big\n" );
562  cdEXIT(EXIT_FAILURE);
563  }
564 
565  StopCalc.HColStop = (realnum)pow(10.,a);
566 
567  /* vary option */
568  if( optimize.lgVarOn )
569  {
571  strcpy( optimize.chVarFmt[optimize.nparm], "STOP COLUMN DENSITY %f" );
572  /* pointer to where to write */
574  /* log of temp will be pointer */
576  optimize.vincr[optimize.nparm] = 0.5;
577  ++optimize.nparm;
578  }
579  }
580  }
581 
582  /* stop when electron density falls below this value, linear or log */
583  else if( nMatch("EDEN",chCard) )
584  {
585  /* stop if electron density falls below this value
586  * LINEAR option */
587  if( nMatch("LINE",chCard) )
588  {
590  }
591  else
592  {
593  StopCalc.StopElecDensity = (realnum)pow(10.,a);
594  }
595  }
596 
597  /* stop at a particular line ratio - this must come last since many commands
598  * have linear option - don't want to trigger on that */
599  else if( nMatch("LINE",chCard) )
600  {
601  char chLabel[5];
602  /* first line wavelength, then intensity relative to Hbeta for stop
603  * if third number is entered, it is wl of line in denominator */
604 
605  /* get label for the line - must do this first so we clear the label string before
606  * trying to read the wavelength */
607  GetQuote( chLabel , chCard , true );
608 
609  /* copy first four char of label into caps, and null terminate*/
610  strncpy( StopCalc.chStopLabel1[StopCalc.nstpl], chLabel , 4 );
612 
613  i = 5;
614  /* get line wavelength */
616  (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH, &lgEOL);
617 
618  /* line was entered, look for possible micron or cm label */
619  if( input.chCARDCAPS[i-1] == 'M' )
620  {
621  /* microns */
623  }
624  else if( input.chCARDCAPS[i-1] == 'C' )
625  {
626  /* microns */
628  }
629 
630  /* get relative intensity */
632  (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH, &lgEOL);
633  if( lgEOL )
634  {
635  fprintf( ioQQQ, " There MUST be a relative intensity entered "
636  "for first line in STOP LINE command. Sorry.\n" );
637  cdEXIT(EXIT_FAILURE);
638  }
639 
640  /* check for second line - use Hbeta is not specified */
641  j = i;
642  a = FFmtRead(chCard,&j,INPUT_LINE_LENGTH, &lgEOL);
643 
644  if( lgEOL )
645  {
646  /* hit the eol, normalization line will be H beta */
647  strncpy( StopCalc.chStopLabel2[StopCalc.nstpl], "TOTL" , 4 );
650  }
651  else
652  {
653  /* get label for the line - must do this first so we clear the label string before
654  * trying to read the wavelength */
655  GetQuote( chLabel , chCard , true );
656  /* copy first four char of label into caps, and null terminate*/
657  strncpy( StopCalc.chStopLabel2[StopCalc.nstpl], chLabel , 4 );
659 
660  /* wavelength of second line, may be absent and so zero -
661  * we will use Hbeta if not specified */
663  (realnum)FFmtRead(chCard,&i, INPUT_LINE_LENGTH,&lgEOL);
664 
665  /* line was entered, look for possible micron or cm label */
666  if( input.chCARDCAPS[i-1] == 'M' )
667  {
668  /* microns */
670  }
671  else if( input.chCARDCAPS[i-1] == 'C' )
672  {
673  /* microns */
675  }
676  }
677  /* increment number of stop lines commands entered */
679  }
680 
681  else if( nMatch("NTOT" , chCard ) && nMatch("ALIO" , chCard ) )
682  {
683  /* this counter is incremented at end of conv_base - total number of
684  * times conv base was called
685  * this is a debugging aid - code aborts */
686  StopCalc.nTotalIonizStop = (long)a;
687  }
688 
689  /* oops! no keyword that we could find */
690  else
691  {
692  fprintf( ioQQQ, " I did not recognize a keyword on this STOP line, line image follows;\n" );
693  fprintf( ioQQQ, " \"%s\"\n Sorry.\n" , chCard);
694  cdEXIT(EXIT_FAILURE);
695  }
696  return;
697 }

Generated for cloudy by doxygen 1.8.4