cloudy  trunk
parse_stop.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 /*ParseStop parse the stop command */
4 #include "cddefines.h"
5 #include "physconst.h"
6 #include "optimize.h"
7 #include "phycon.h"
8 #include "predcont.h"
9 #include "rfield.h"
10 #include "radius.h"
11 #include "geometry.h"
12 #include "iterations.h"
13 #include "stopcalc.h"
14 #include "input.h"
15 #include "parser.h"
16 
17 void ParseStop(Parser &p)
18 {
19  long int j;
20 
21  double effcol,
22  tread;
23 
24  DEBUG_ENTRY( "ParseStop()" );
25 
26  /* time option, for stopping time dependent calculations, used to stop
27  * iterations rather than zones. Only some stop commands have this option */
28  bool lgStopZone = true;
29  if( p.nMatch("TIME") )
30  lgStopZone = false;
31 
32  if( p.nMatch("TEMP") )
33  {
34  double a = p.FFmtRead();
35 
36  if( p.lgEOL() && !p.nMatch(" OFF") )
37  {
38  p.NoNumb("temperature");
39  }
40 
41  /* off option disables this stopping criterion */
42  if( p.lgEOL() && p.nMatch(" OFF") )
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  if( lgStopZone )
48  StopCalc.TempLoStopZone = -1.f;
49  else
50  StopCalc.TempLoStopIteration = -1.f;
51  }
52  else
53  {
54  /* lowest electron temperature allowed before stopping
55  * assumed to be the log of the temperature if <10
56  * optional keyword LINEAR forces linear */
57  if( a <= 10. && !p.nMatch("LINE") )
58  {
59  tread = pow(10.,a);
60  }
61  else
62  {
63  tread = a;
64  }
65 
66  /* tread is linear temperature*/
67  if( tread < phycon.TEMP_LIMIT_LOW )
68  {
69  fprintf( ioQQQ,
70  " Temperatures below %.2e K not allowed. Reset to lowest value."
71  " I am doing this myself.\n" ,
72  phycon.TEMP_LIMIT_LOW );
73  /* set slightly off extreme limit for safety */
74  tread = phycon.TEMP_LIMIT_LOW*1.01;
75  }
76  else if( tread > phycon.TEMP_LIMIT_HIGH )
77  {
78  fprintf( ioQQQ,
79  " Temperatures is above %.2e K not allowed. Reset to highest value."
80  " I am doing this myself.\n" ,
81  phycon.TEMP_LIMIT_HIGH);
82  /* set slightly off extreme limit for safety */
83  tread = phycon.TEMP_LIMIT_HIGH*0.99;
84  }
85 
86  if( p.nMatch("EXCE") )
87  {
88  /* option for this to be highest allowed temperature,
89  * stop temperate exceeds */
90  if( lgStopZone )
91  StopCalc.TempHiStopZone = (realnum)tread;
92  else
93  StopCalc.TempHiStopIteration = (realnum)tread;
94  }
95  else
96  {
97  /* this is ending temperature - we stop if kinetic temperature
98  * falls below this */
99  if( lgStopZone )
100  StopCalc.TempLoStopZone = (realnum)tread;
101  else
102  StopCalc.TempLoStopIteration = (realnum)tread;
103  }
104  }
105  }
106 
107  /* stop at 21cm line center optical depth */
108  else if( p.nMatch("OPTI") && p.nMatch("21CM") )
109  {
110  /* default is for number to be log of optical depth */
111  bool lgLOG = true;
112  if( p.nMatch("LINE") )
113  {
114  /* force it to be linear not log */
115  lgLOG = false;
116  }
117  j = (long int)p.FFmtRead();
118  if( j!=21 )
119  {
120  fprintf( ioQQQ, " First number on STOP 21CM OPTICAL DEPTH command must be 21\n" );
121  cdEXIT(EXIT_FAILURE);
122  }
123  /* now get the next number, which is the optical depth */
124  double a = (long int)p.FFmtRead();
125 
126  /* tau must be a log, and second number is energy where tau specified */
127  if( lgLOG )
128  {
129  StopCalc.tauend = (realnum)pow(10.,a);
130  }
131  else
132  {
133  StopCalc.tauend = (realnum)a;
134  }
135  /* this flag says that 21cm line optical depth is the stop quantity */
136  StopCalc.lgStop21cm = true;
137  }
138  /* stop optical depth at some energy */
139  else if( p.nMatch("OPTI") )
140  {
141  double a = p.FFmtRead();
142 
143  if( p.lgEOL() && !p.nMatch(" OFF") )
144  {
145  p.NoNumb("optical depth");
146  }
147 
148  /* default is for number to be log of optical depth */
149  bool lgLOG = true;
150  if( p.nMatch("LINE") )
151  {
152  /* force it to be linear not log */
153  lgLOG = false;
154  }
155 
156  if( a > 37. )
157  {
158  fprintf( ioQQQ, " optical depth too big\n" );
159  cdEXIT(EXIT_FAILURE);
160  }
161 
162  /* tau entered as a log unless liner specified */
163  if( lgLOG )
164  {
165  StopCalc.tauend = (realnum)pow(10.,a);
166  }
167  else
168  {
169  StopCalc.tauend = (realnum)a;
170  }
171 
172  /* energy where tau specified */
173  StopCalc.taunu = (realnum)p.FFmtRead();
174 
175  if( p.lgEOL() )
176  {
177  if( p.nMatch("LYMA") )
178  {
179  /* largest Lyman limit optical depth */
180  StopCalc.taunu = 1.;
181  }
182  else if( p.nMatch("BALM") )
183  {
184  /* stop at this Balmer continuum optical depth */
185  StopCalc.taunu = 0.25;
186  }
187  else
188  {
189  fprintf( ioQQQ, " There must be a second number, the energy in Ryd. Sorry.\n" );
190  cdEXIT(EXIT_FAILURE);
191  }
192  }
193 
194  else
195  {
196  /* if second number is negative then log of energy in rydbergs */
197  if( StopCalc.taunu < 0. )
198  {
199  StopCalc.taunu = (realnum)pow((realnum)10.f,StopCalc.taunu);
200  }
201 
202  /* check that energy is within bounds of code */
203  if( StopCalc.taunu < rfield.emm || StopCalc.taunu > rfield.egamry )
204  {
205  fprintf( ioQQQ, " The energy must be in the range %10.2e to %10.2e. It was %10.2e. Sorry.\n",
206  rfield.emm, rfield.egamry, StopCalc.taunu );
207  cdEXIT(EXIT_FAILURE);
208  }
209  }
210 
211  /* vary option */
212  if( optimize.lgVarOn )
213  {
214  strcpy( optimize.chVarFmt[optimize.nparm], "STOP OPTICAL DEPTH = %f LOG AT %f RYD" );
215  /* pointer to where to write */
216  optimize.nvfpnt[optimize.nparm] = input.nRead;
217  optimize.vparm[0][optimize.nparm] = (realnum)log10(StopCalc.tauend);
218  optimize.vparm[1][optimize.nparm] = (realnum)StopCalc.taunu;
219  optimize.vincr[optimize.nparm] = 0.5;
220  optimize.nvarxt[optimize.nparm] = 2;
221  ++optimize.nparm;
222  }
223  }
224 
225  /* stop optical depth at extinction in V filter */
226  else if( p.nMatch(" AV ") )
227  {
228  double a = p.FFmtRead();
229 
230  if( p.lgEOL() && !p.nMatch(" OFF") )
231  {
232  p.NoNumb("optical depth in V");
233  }
234  /* default is for number to be A_V, log if negative */
235  if( a<=0. )
236  {
237  a = pow(10.,a);
238  }
239  /* A_V can be for either point or extended source, difference is (1-g) multiplied by scat opacity
240  * if keyword point occurs then for point source, otherwise for extended source */
241  if( p.nMatch("EXTE" ) )
242  {
243  StopCalc.AV_extended = (realnum)a;
244  }
245  else
246  {
247  /* default is point, as measured in ism work */
248  StopCalc.AV_point = (realnum)a;
249  }
250  }
251 
252  /* stop when a fraction of molecules frozen out on grain surfaces is reached */
253  else if( p.nMatch("MOLE") && p.nMatch("DEPL") )
254  {
255  double a = p.FFmtRead();
256 
257  if( p.lgEOL() && !p.nMatch(" OFF") )
258  {
259  p.NoNumb("molecular depletion");
260  }
261  if( a <= 0. )
262  {
263  StopCalc.StopDepleteFrac = (realnum)pow(10.,a);
264  }
265  else
266  {
267  StopCalc.StopDepleteFrac = (realnum)a;
268  }
269  }
270 
271  /* stop when absolute value of flow velocity falls below this value */
272  else if( p.nMatch("VELO") )
273  {
274  double a = p.FFmtRead();
275 
276  if( p.lgEOL() && !p.nMatch(" OFF") )
277  {
278  p.NoNumb("flow velocity");
279  }
280  /* entered in km/s but stored as cm/s */
281  StopCalc.StopVelocity = (realnum)(a*1e5);
282  }
283 
284  /* stop at a given computed mass */
285  else if( p.nMatch("MASS") )
286  {
287  double a = p.FFmtRead();
288 
289  if( p.lgEOL() && !p.nMatch(" OFF") )
290  {
291  p.NoNumb("mass");
292  }
293  /* number of log of mass in gm if inner radius is specified,
294  * mass per unit area, gm cm-2 if not
295  * leave it as a log since dare not deal with linear mass */
296  StopCalc.xMass = (realnum)a;
297  /* NB 0 is sentinel for not set, if a is zero we must reset it */
298  if( StopCalc.xMass == 0 )
299  StopCalc.xMass = SMALLFLOAT;
300 
301  /* vary option */
302  if( optimize.lgVarOn )
303  {
304  strcpy( optimize.chVarFmt[optimize.nparm], "STOP MASS = %f LOG" );
305  /* pointer to where to write */
306  optimize.nvfpnt[optimize.nparm] = input.nRead;
307  optimize.vparm[0][optimize.nparm] = (realnum)StopCalc.xMass;
308  optimize.vincr[optimize.nparm] = 0.5;
309  optimize.nvarxt[optimize.nparm] = 1;
310  ++optimize.nparm;
311  }
312  }
313 
314  /* stop thickness command, also stop depth, this must come after stop
315  * optical depth, since do not want to trigger on depth in optical depth */
316  else if( p.nMatch("THIC") || p.nMatch("DEPT") || p.nMatch("RADI") )
317  {
318  double a = p.FFmtRead();
319 
320  if( p.lgEOL() && !p.nMatch(" OFF") )
321  {
322  p.NoNumb("distance");
323  }
324  const double convl = p.nMatch("PARS") ? log10( PARSEC ) : 0.;
325  bool lgStopRadius = p.nMatch("RADI") ? true : false ;
326  const char* what = lgStopRadius ? "radius" : "thickness";
327 
328  if( p.nMatch("LINE") )
329  {
330  if( a > 0. )
331  {
332  a = log10(a) + convl;
333  }
334  else
335  {
336  fprintf(ioQQQ,"The %s is negative and linear is set - this is impossible.\n", what);
337  cdEXIT(EXIT_FAILURE);
338  }
339  }
340  else
341  {
342  a += convl;
343  }
344  if( a > 37. )
345  {
346  fprintf( ioQQQ, "DISASTER %s too large\n", what );
347  cdEXIT(EXIT_FAILURE);
348  }
349  if( lgStopRadius )
350  radius.StopRadius[0] = pow(10.,a);
351  else
352  radius.StopThickness[0] = pow(10.,a);
353 
354  /* can stop at different thickness on each iteration */
355  for( j=1; j < iterations.iter_malloc; j++ )
356  {
357  a = p.FFmtRead();
358  if( p.lgEOL() )
359  {
360  if( lgStopRadius )
361  radius.StopRadius[j] = radius.StopRadius[j-1];
362  else
363  radius.StopThickness[j] = radius.StopThickness[j-1];
364  }
365  else
366  {
367  if( p.nMatch("LINE") )
368  {
369  if( a > 0. )
370  {
371  a = log10(a) + convl;
372  }
373  else
374  {
375  fprintf(ioQQQ,"The %s is negative and linear is set -"
376  " this is impossible.\n", what);
377  cdEXIT(EXIT_FAILURE);
378  }
379  }
380  else
381  {
382  a += convl;
383  }
384  if( a > 37. )
385  {
386  fprintf( ioQQQ, "DISASTER %s too large\n", what );
387  cdEXIT(EXIT_FAILURE);
388  }
389  if( lgStopRadius )
390  radius.StopRadius[j] = pow(10.,a);
391  else
392  radius.StopThickness[j] = pow(10.,a);
393  }
394  }
395 
396  /* vary option */
397  if( optimize.lgVarOn )
398  {
399  optimize.nvarxt[optimize.nparm] = 1;
400  /* pointer to where to write */
401  optimize.nvfpnt[optimize.nparm] = input.nRead;
402  if( lgStopRadius )
403  {
404  strcpy( optimize.chVarFmt[optimize.nparm], "STOP RADIUS %f LOG" );
405  optimize.vparm[0][optimize.nparm] = (realnum)log10(radius.StopRadius[0]);
406  }
407  else
408  {
409  strcpy( optimize.chVarFmt[optimize.nparm], "STOP THICKNESS %f LOG" );
410  optimize.vparm[0][optimize.nparm] = (realnum)log10(radius.StopThickness[0]);
411  }
412  optimize.vincr[optimize.nparm] = 0.5f;
413  ++optimize.nparm;
414  }
415  }
416 
417  /* stop at a particular zone, for each iteration */
418  else if( p.nMatch("ZONE") )
419  {
420  double a = p.FFmtRead();
421 
422  if( p.lgEOL() && !p.nMatch(" OFF") )
423  {
424  p.NoNumb("zone number");
425  }
426  /* stop after computing this zone */
427  /* >>chng 03 jun 06, do not let fall below 1, stop zone 0 has same effect
428  * as stop zone 1, bug caught by Joop Schaye */
429  geometry.nend[0] = (long)MAX2(1.,a);
430  geometry.lgZoneSet = true;
431 
432  /* this tells code that we intend to stop at this zone, so caution not generated*/
433  geometry.lgEndDflt = false;
434 
435  long int nZoneMax = geometry.nend[0];
436  for( j=1; j < iterations.iter_malloc; j++ )
437  {
438  geometry.nend[j] = (long)p.FFmtRead();
439  /* if eol on this iteration, set to previous. In most cases
440  * all will be equal to the first */
441  if( p.lgEOL() )
442  {
443  geometry.nend[j] = geometry.nend[j-1];
444  }
445  else
446  {
447  /* do not let fall below 1, stop zone 0 has same effect
448  * as stop zone 1, bug caught by Joop Schaye */
449  geometry.nend[j] = MAX2( 1 , geometry.nend[j] );
450  }
451  nZoneMax = max( nZoneMax , geometry.nend[j] );
452  }
453 
454  if( nZoneMax>2000 )
455  fprintf(ioQQQ,"CAUTION - it will take a lot of memory to save"
456  " results for %li zones. Is this many zones really necessary?\n",
457  nZoneMax );
458  }
459 
460  /* stop when a prescribed continuum flux is reached */
461  else if( p.nMatch("CONT") && p.nMatch("FLUX") )
462  {
463  /* first read the continuum energy and add this point to PredCont */
464  double energy = p.FFmtRead();
465  if( p.lgEOL() )
466  p.NoNumb("energy");
467  const char* unit = p.StandardEnergyUnit();
468  long ind = t_PredCont::Inst().add( energy, unit );
469  Energy E( energy, unit );
470 
471  double flux = p.FFmtRead();
472  if( p.lgEOL() )
473  p.NoNumb("flux");
474  if( flux <= 0. || p.nMatch( " LOG") )
475  flux = pow(10.,flux);
476  Flux F( E, flux, p.StandardFluxUnit() );
477 
478  StopCalc.ContIndex.push_back( ind );
479  StopCalc.ContNFnu.push_back( F );
480  }
481 
482  /* stop at this electron fraction, relative to hydrogen */
483  else if( p.nMatch("EFRA") )
484  {
485  double a = p.FFmtRead();
486 
487  if( p.lgEOL() && !p.nMatch(" OFF") )
488  {
489  p.NoNumb("electron fraction");
490  }
491  if( a <= 0. )
492  {
493  StopCalc.StopElecFrac = (realnum)pow(10.,a);
494  }
495  else
496  {
497  StopCalc.StopElecFrac = (realnum)a;
498  }
499  }
500 
501  /* stop at a hydrogen molecular fraction, relative to total hydrogen,
502  * this is 2H_2 / H_total*/
503  else if( p.nMatch("MFRA") )
504  {
505  double a = p.FFmtRead();
506 
507  if( p.lgEOL() && !p.nMatch(" OFF") )
508  {
509  p.NoNumb("hydrogen molecular fraction");
510  }
511  if( a <= 0. )
512  {
513  StopCalc.StopH2MoleFrac = (realnum)pow(10.,a);
514  }
515  else
516  {
517  StopCalc.StopH2MoleFrac = (realnum)a;
518  }
519  }
520 
521  /* stop at a ionized hydrogen fraction, relative to total hydrogen,
522  * this is H+ / H_total */
523  else if( p.nMatch("PFRA") )
524  {
525  double a = p.FFmtRead();
526 
527  if( p.lgEOL() && !p.nMatch(" OFF") )
528  {
529  p.NoNumb("ionized hydrogen fraction");
530  }
531  if( a <= 0. )
532  {
533  StopCalc.StopHPlusFrac = (realnum)pow(10.,a);
534  }
535  else
536  {
537  StopCalc.StopHPlusFrac = (realnum)a;
538  }
539  }
540 
541  /* stop at a particular column density */
542  else if( p.nMatch("COLU") )
543  {
544  double a = p.FFmtRead();
545 
546  if( p.lgEOL() && !p.nMatch(" OFF") )
547  {
548  p.NoNumb("column density");
549  }
550  char chLabel[CHARS_SPECIES];
551 
552  /* check for linear option, if present take log since a being
553  * log column density is default */
554  if( p.nMatch( "LINE" ) )
555  a = log10(a);
556 
557  if( !p.GetQuote( chLabel , false ) )
558  {
559  /* species was given in double quotes */
560  StopCalc.col_species = (realnum)pow(10.,a);
561  StopCalc.lgStopSpeciesColumn = true;
562  strncpy( StopCalc.chSpeciesColumn, chLabel, CHARS_SPECIES );
563  }
564  /* stop at an effective column density */
565  else if( p.nMatch("EFFE") )
566  {
567  /* actually stop at certain optical depth at 1keV */
568  effcol = pow(10.,a);
569  StopCalc.tauend = (realnum)(effcol*2.14e-22);
570  StopCalc.taunu = (realnum)(1000./EVRYD);
571  /* vary option */
572  if( optimize.lgVarOn )
573  {
574  optimize.nvarxt[optimize.nparm] = 1;
575  strcpy( optimize.chVarFmt[optimize.nparm], "STOP EFFECTIVE COLUMN DENSITY %f LOG" );
576  /* pointer to where to write */
577  optimize.nvfpnt[optimize.nparm] = input.nRead;
578  /* log of temp will be pointer */
579  optimize.vparm[0][optimize.nparm] = (realnum)log10(effcol);
580  optimize.vincr[optimize.nparm] = 0.5f;
581  ++optimize.nparm;
582  }
583  }
584 
585  else if( p.nMatch("IONI") )
586  {
587  /* this is ionized column */
588  if( a > 37. )
589  {
590  fprintf( ioQQQ, " column too big\n" );
591  cdEXIT(EXIT_FAILURE);
592  }
593 
594  StopCalc.colpls = (realnum)pow(10.,a);
595 
596  /* vary option */
597  if( optimize.lgVarOn )
598  {
599  optimize.nvarxt[optimize.nparm] = 1;
600  strcpy( optimize.chVarFmt[optimize.nparm], "STOP IONIZED COLUMN DENSITY %f LOG" );
601  /* pointer to where to write */
602  optimize.nvfpnt[optimize.nparm] = input.nRead;
603  /* log of temp will be pointer */
604  optimize.vparm[0][optimize.nparm] = (realnum)log10(StopCalc.colpls);
605  optimize.vincr[optimize.nparm] = 0.5f;
606  ++optimize.nparm;
607  }
608  }
609 
610  /* stop at a neutral column */
611  else if( p.nMatch("NEUT") )
612  {
613  StopCalc.colnut = (realnum)pow(10.,a);
614 
615  /* vary option */
616  if( optimize.lgVarOn )
617  {
618  optimize.nvarxt[optimize.nparm] = 1;
619  strcpy( optimize.chVarFmt[optimize.nparm], "STOP NEUTRAL COLUMN DENSITY %f LOG");
620  /* pointer to where to write */
621  optimize.nvfpnt[optimize.nparm] = input.nRead;
622  /* log of temp will be pointer */
623  optimize.vparm[0][optimize.nparm] = (realnum)log10(StopCalc.colnut);
624  optimize.vincr[optimize.nparm] = 0.5f;
625  ++optimize.nparm;
626  }
627  }
628 
629  /* >>chng 03 apr 15, add this option
630  * stop at a molecular hydrogen column density, input parameter
631  * is log of the H2 column density */
632  else if( p.nMatch(" H2 ") )
633  {
634  /* this command has a 2 in the H2 label - must not parse the two by
635  * accident. Get the first number off the line image, and confirm that
636  * it is a 2 */
637  j = (long int)a;
638  if( j != 2 )
639  {
640  fprintf( ioQQQ, " Something is wrong with the order of the numbers on this line.\n" );
641  fprintf( ioQQQ, " The first number I encounter should be the 2 in H2.\n Sorry.\n" );
642  cdEXIT(EXIT_FAILURE);
643  }
644  a = p.FFmtRead();
645  StopCalc.col_h2 = (realnum)pow(10.,a);
646 
647  /* vary option */
648  if( optimize.lgVarOn )
649  {
650  optimize.nvarxt[optimize.nparm] = 1;
651  strcpy( optimize.chVarFmt[optimize.nparm], "STOP H2 COLUMN DENSITY %f LOG");
652  /* pointer to where to write */
653  optimize.nvfpnt[optimize.nparm] = input.nRead;
654  /* log of temp will be pointer */
655  optimize.vparm[0][optimize.nparm] = (realnum)log10(StopCalc.col_h2);
656  optimize.vincr[optimize.nparm] = 0.5f;
657  ++optimize.nparm;
658  }
659  }
660 
661  else if( p.nMatch("ATOM") )
662  {
663  StopCalc.col_h2_nut = (realnum)pow(10.,a);
664  /* vary option */
665  if( optimize.lgVarOn )
666  {
667  optimize.nvarxt[optimize.nparm] = 1;
668  strcpy( optimize.chVarFmt[optimize.nparm], "STOP ATOMIC COLUMN DENSITY %f LOG");
669  /* pointer to where to write */
670  optimize.nvfpnt[optimize.nparm] = input.nRead;
671  /* log of temp will be pointer */
672  optimize.vparm[0][optimize.nparm] = (realnum)log10(StopCalc.col_h2_nut);
673  optimize.vincr[optimize.nparm] = 0.5f;
674  ++optimize.nparm;
675  }
676  }
677 
678  else if( p.nMatch("H/TS") )
679  {
680  /* >> 05 jan 09, add stop integrated n(H0) / Tspin */
681  StopCalc.col_H0_ov_Tspin = (realnum)pow(10.,a);
682  /* vary option */
683  if( optimize.lgVarOn )
684  {
685  optimize.nvarxt[optimize.nparm] = 1;
686  strcpy( optimize.chVarFmt[optimize.nparm], "STOP H/TSPIN COLUMN DENSITY %f LOG");
687  /* pointer to where to write */
688  optimize.nvfpnt[optimize.nparm] = input.nRead;
689  /* log of temp will be pointer */
690  optimize.vparm[0][optimize.nparm] = (realnum)log10(StopCalc.col_H0_ov_Tspin);
691  optimize.vincr[optimize.nparm] = 0.5f;
692  ++optimize.nparm;
693  }
694  }
695 
696  else if( p.nMatch(" CO ") )
697  {
698  /* chng << 03 Oct. 27--Nick Abel, add this option */
699  /* stop at a carbon monoxide column density */
700  StopCalc.col_monoxco = (realnum)pow(10.,a);
701  /* vary option */
702  if( optimize.lgVarOn )
703  {
704  optimize.nvarxt[optimize.nparm] = 1;
705  strcpy( optimize.chVarFmt[optimize.nparm], "STOP CO COLUMN DENSITY %f LOG");
706  /* pointer to where to write */
707  optimize.nvfpnt[optimize.nparm] = input.nRead;
708  /* log of temp will be pointer */
709  optimize.vparm[0][optimize.nparm] = (realnum)log10(StopCalc.col_monoxco);
710  optimize.vincr[optimize.nparm] = 0.5f;
711  ++optimize.nparm;
712  }
713  }
714 
715  /* fall through default is total hydrogen column density */
716  else
717  {
718  /* both HII and HI */
719  if( a > 37. )
720  {
721  fprintf( ioQQQ, " column too big\n" );
722  cdEXIT(EXIT_FAILURE);
723  }
724 
725  StopCalc.HColStop = (realnum)pow(10.,a);
726 
727  /* vary option */
728  if( optimize.lgVarOn )
729  {
730  optimize.nvarxt[optimize.nparm] = 1;
731  strcpy( optimize.chVarFmt[optimize.nparm], "STOP COLUMN DENSITY %f LOG" );
732  /* pointer to where to write */
733  optimize.nvfpnt[optimize.nparm] = input.nRead;
734  /* log of temp will be pointer */
735  optimize.vparm[0][optimize.nparm] = (realnum)log10(StopCalc.HColStop);
736  optimize.vincr[optimize.nparm] = 0.5f;
737  ++optimize.nparm;
738  }
739  }
740  }
741 
742  /* stop when electron density falls below this value, linear or log */
743  else if( p.nMatch("EDEN") )
744  {
745  double a = p.FFmtRead();
746 
747  if( p.lgEOL() && !p.nMatch(" OFF") )
748  {
749  p.NoNumb("electron density");
750  }
751  /* stop if electron density falls below this value
752  * LINEAR option */
753  if( p.nMatch("LINE") )
754  {
755  StopCalc.StopElecDensity = (realnum)a;
756  }
757  else
758  {
759  StopCalc.StopElecDensity = (realnum)pow(10.,a);
760  }
761  }
762 
763  /* stop at a particular line ratio - this must come last since many commands
764  * have linear option - don't want to trigger on that */
765  else if( p.nMatch("LINE") )
766  {
767  char chLabel[5];
768  /* first line wavelength, then intensity relative to Hbeta for stop
769  * if third number is entered, it is wl of line in denominator */
770 
771  /* get label for the line - must do this first so we clear the label string before
772  * trying to read the wavelength */
773  p.GetQuote( chLabel , true );
774 
775  /* copy first four char of label into caps, and null terminate*/
776  strncpy( StopCalc.chStopLabel1[StopCalc.nstpl], chLabel , 4 );
777  StopCalc.chStopLabel1[StopCalc.nstpl][4] = 0;
778 
779  // default intrinsic intensity, accept emergent
780  StopCalc.nEmergent[StopCalc.nstpl] = 0;
781  if( p.nMatch("EMER") )
782  StopCalc.nEmergent[StopCalc.nstpl] = 1;
783 
784  /* get line wavelength */
785  StopCalc.StopLineWl1[StopCalc.nstpl] = (realnum)p.getWaveOpt();
786 
787  /* get relative intensity */
788  StopCalc.stpint[StopCalc.nstpl] = (realnum)p.FFmtRead();
789  if( p.lgEOL() )
790  {
791  fprintf( ioQQQ, " There MUST be a relative intensity entered "
792  "for first line in STOP LINE command. Sorry.\n" );
793  cdEXIT(EXIT_FAILURE);
794  }
795 
796  /* check for second line - use Hbeta is not specified */
797 
798  /* get label for the line - must do this first so we clear the label string before
799  * trying to read the wavelength */
800 
801  j = p.GetQuote( chLabel , false );
802 
803  if( j != 0 )
804  {
805  /* no line label, normalization line will be H beta */
806  strncpy( StopCalc.chStopLabel2[StopCalc.nstpl], "TOTL" , 4 );
807  StopCalc.chStopLabel2[StopCalc.nstpl][4] = 0;
808  StopCalc.StopLineWl2[StopCalc.nstpl] = 4861.f;
809  }
810  else
811  {
812  /* copy first four char of label into caps, and null terminate*/
813  strncpy( StopCalc.chStopLabel2[StopCalc.nstpl], chLabel , 4 );
814  StopCalc.chStopLabel2[StopCalc.nstpl][4] = 0;
815 
816  /* wavelength of second line, may be absent and so zero -
817  * we will use Hbeta if not specified */
818  StopCalc.StopLineWl2[StopCalc.nstpl] = (realnum)p.getWaveOpt();
819 
820  }
821  /* increment number of stop lines commands entered */
822  StopCalc.nstpl = MIN2(StopCalc.nstpl+1,MXSTPL-1);
823  }
824 
825  else if( p.nMatch("NTOTALIO" ) )
826  {
827  double a = p.FFmtRead();
828 
829  if( p.lgEOL() && !p.nMatch(" OFF") )
830  {
831  p.NoNumb("number of calls to conv_base");
832  }
833  /* this counter is incremented at end of conv_base - total number of
834  * times conv base was called
835  * this is a debugging aid - code aborts */
836  StopCalc.nTotalIonizStop = (long)a;
837  }
838 
839  /* oops! no keyword that we could find */
840  else
841  {
842  fprintf( ioQQQ, " I did not recognize a keyword on this STOP line, line image follows;\n" );
843  p.PrintLine(ioQQQ);
844  fprintf( ioQQQ, "Sorry.\n");
845  cdEXIT(EXIT_FAILURE);
846  }
847  return;
848 }
t_StopCalc::col_H0_ov_Tspin
realnum col_H0_ov_Tspin
Definition: stopcalc.h:83
MXSTPL
const int MXSTPL
Definition: stopcalc.h:10
t_StopCalc::col_h2_nut
realnum col_h2_nut
Definition: stopcalc.h:80
Parser::nMatch
bool nMatch(const char *chKey) const
Definition: parser.h:135
StopCalc
t_StopCalc StopCalc
Definition: stopcalc.cpp:5
t_StopCalc::HColStop
realnum HColStop
Definition: stopcalc.h:69
t_geometry::lgZoneSet
bool lgZoneSet
Definition: geometry.h:87
t_optimize::vincr
realnum vincr[LIMPAR]
Definition: optimize.h:191
Parser::FFmtRead
double FFmtRead(void)
Definition: parser.cpp:353
t_StopCalc::chStopLabel1
char chStopLabel1[MXSTPL][5]
Definition: stopcalc.h:115
t_StopCalc::colpls
realnum colpls
Definition: stopcalc.h:70
Parser::StandardFluxUnit
string StandardFluxUnit(void) const
Definition: parser.cpp:178
t_optimize::nparm
long int nparm
Definition: optimize.h:201
t_StopCalc::StopVelocity
realnum StopVelocity
Definition: stopcalc.h:65
t_radius::StopRadius
double * StopRadius
Definition: radius.h:58
t_StopCalc::nEmergent
int nEmergent[MXSTPL]
Definition: stopcalc.h:117
Singleton< t_PredCont >::Inst
static t_PredCont & Inst()
Definition: cddefines.h:175
t_input::nRead
long int nRead
Definition: input.h:49
rfield
t_rfield rfield
Definition: rfield.cpp:8
ioQQQ
FILE * ioQQQ
Definition: cddefines.cpp:7
geometry.h
realnum
float realnum
Definition: cddefines.h:103
iterations
t_iterations iterations
Definition: iterations.cpp:5
rfield.h
ParseStop
void ParseStop(Parser &p)
Definition: parse_stop.cpp:17
t_StopCalc::taunu
realnum taunu
Definition: stopcalc.h:26
Parser::StandardEnergyUnit
const char * StandardEnergyUnit(void) const
Definition: parser.cpp:174
t_StopCalc::lgStop21cm
bool lgStop21cm
Definition: stopcalc.h:120
t_optimize::lgVarOn
bool lgVarOn
Definition: optimize.h:203
t_geometry::lgEndDflt
bool lgEndDflt
Definition: geometry.h:96
phycon
t_phycon phycon
Definition: phycon.cpp:6
t_StopCalc::colnut
realnum colnut
Definition: stopcalc.h:71
t_StopCalc::ContIndex
vector< long > ContIndex
Definition: stopcalc.h:123
t_phycon::TEMP_LIMIT_LOW
const double TEMP_LIMIT_LOW
Definition: phycon.h:111
Parser::GetQuote
int GetQuote(char *chLabel, bool lgABORT)
Definition: parser.h:209
t_StopCalc::AV_point
realnum AV_point
Definition: stopcalc.h:89
input
t_input input
Definition: input.cpp:12
t_StopCalc::lgStopSpeciesColumn
bool lgStopSpeciesColumn
Definition: stopcalc.h:132
predcont.h
MIN2
#define MIN2
Definition: cddefines.h:761
radius
t_radius radius
Definition: radius.cpp:5
t_StopCalc::ContNFnu
vector< Flux > ContNFnu
Definition: stopcalc.h:124
t_StopCalc::StopDepleteFrac
realnum StopDepleteFrac
Definition: stopcalc.h:61
CHARS_SPECIES
@ CHARS_SPECIES
Definition: cddefines.h:274
optimize
t_optimize optimize
Definition: optimize.cpp:5
t_StopCalc::col_h2
realnum col_h2
Definition: stopcalc.h:74
t_StopCalc::StopHPlusFrac
realnum StopHPlusFrac
Definition: stopcalc.h:56
t_optimize::vparm
realnum vparm[LIMEXT][LIMPAR]
Definition: optimize.h:188
t_optimize::chVarFmt
char chVarFmt[LIMPAR][FILENAME_PATH_LENGTH_2]
Definition: optimize.h:263
t_StopCalc::AV_extended
realnum AV_extended
Definition: stopcalc.h:89
EXIT_FAILURE
#define EXIT_FAILURE
Definition: cddefines.h:140
t_rfield::egamry
realnum egamry
Definition: rfield.h:52
Parser::NoNumb
NORETURN void NoNumb(const char *chDesc) const
Definition: parser.cpp:233
PARSEC
const UNUSED double PARSEC
Definition: physconst.h:138
Parser
Definition: parser.h:31
t_StopCalc::StopElecDensity
realnum StopElecDensity
Definition: stopcalc.h:92
t_StopCalc::StopLineWl1
realnum StopLineWl1[MXSTPL]
Definition: stopcalc.h:111
t_StopCalc::col_monoxco
realnum col_monoxco
Definition: stopcalc.h:86
cddefines.h
optimize.h
t_StopCalc::TempHiStopZone
realnum TempHiStopZone
Definition: stopcalc.h:36
t_StopCalc::chSpeciesColumn
char chSpeciesColumn[CHARS_SPECIES]
Definition: stopcalc.h:133
t_optimize::nvarxt
long int nvarxt[LIMPAR]
Definition: optimize.h:194
radius.h
MAX2
#define MAX2
Definition: cddefines.h:782
t_StopCalc::stpint
realnum stpint[MXSTPL]
Definition: stopcalc.h:103
t_StopCalc::tauend
realnum tauend
Definition: stopcalc.h:23
cdEXIT
#define cdEXIT(FAIL)
Definition: cddefines.h:434
Parser::lgEOL
bool lgEOL(void) const
Definition: parser.h:98
t_StopCalc::TempLoStopZone
realnum TempLoStopZone
Definition: stopcalc.h:42
Flux
Definition: flux.h:9
t_PredCont::add
long add(double energy, const char *unit="Ryd")
Definition: predcont.cpp:122
Energy
Definition: energy.h:7
t_StopCalc::xMass
realnum xMass
Definition: stopcalc.h:77
t_StopCalc::TempHiStopIteration
realnum TempHiStopIteration
Definition: stopcalc.h:38
parser.h
t_radius::StopThickness
double * StopThickness
Definition: radius.h:55
physconst.h
t_optimize::nvfpnt
long int nvfpnt[LIMPAR]
Definition: optimize.h:195
t_StopCalc::nTotalIonizStop
long int nTotalIonizStop
Definition: stopcalc.h:127
t_StopCalc::TempLoStopIteration
realnum TempLoStopIteration
Definition: stopcalc.h:45
phycon.h
geometry
t_geometry geometry
Definition: geometry.cpp:5
iterations.h
t_StopCalc::col_species
realnum col_species
Definition: stopcalc.h:134
t_phycon::TEMP_LIMIT_HIGH
const double TEMP_LIMIT_HIGH
Definition: phycon.h:113
t_rfield::emm
realnum emm
Definition: rfield.h:49
t_StopCalc::StopElecFrac
realnum StopElecFrac
Definition: stopcalc.h:48
t_StopCalc::chStopLabel2
char chStopLabel2[MXSTPL][5]
Definition: stopcalc.h:116
t_iterations::iter_malloc
long int iter_malloc
Definition: iterations.h:29
EVRYD
const UNUSED double EVRYD
Definition: physconst.h:189
t_StopCalc::nstpl
long int nstpl
Definition: stopcalc.h:109
Parser::PrintLine
int PrintLine(FILE *fp) const
Definition: parser.h:204
t_StopCalc::StopH2MoleFrac
realnum StopH2MoleFrac
Definition: stopcalc.h:52
stopcalc.h
input.h
t_StopCalc::StopLineWl2
realnum StopLineWl2[MXSTPL]
Definition: stopcalc.h:112
max
long max(int a, long b)
Definition: cddefines.h:775
DEBUG_ENTRY
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:684
Parser::getWaveOpt
double getWaveOpt()
Definition: parser.cpp:244
SMALLFLOAT
const realnum SMALLFLOAT
Definition: cpu.h:191
t_geometry::nend
long int * nend
Definition: geometry.h:80
Generated by   doxygen 1.8.17