cloudy trunk
Loading...
Searching...
No Matches
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
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" );
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" );
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" );
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 );
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);
338 }
339 }
340 else
341 {
342 a += convl;
343 }
344 if( a > 37. )
345 {
346 fprintf( ioQQQ, "DISASTER %s too large\n", what );
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);
378 }
379 }
380 else
381 {
382 a += convl;
383 }
384 if( a > 37. )
385 {
386 fprintf( ioQQQ, "DISASTER %s too large\n", what );
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" );
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" );
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" );
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" );
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");
846 }
847 return;
848}
FILE * ioQQQ
Definition cddefines.cpp:7
@ CHARS_SPECIES
Definition cddefines.h:274
#define MIN2
Definition cddefines.h:761
#define EXIT_FAILURE
Definition cddefines.h:140
#define cdEXIT(FAIL)
Definition cddefines.h:434
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
long max(int a, long b)
Definition cddefines.h:775
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
Definition energy.h:8
Definition flux.h:10
double FFmtRead(void)
Definition parser.cpp:353
bool nMatch(const char *chKey) const
Definition parser.h:135
double getWaveOpt()
Definition parser.cpp:244
const char * StandardEnergyUnit(void) const
Definition parser.cpp:174
bool lgEOL(void) const
Definition parser.h:98
string StandardFluxUnit(void) const
Definition parser.cpp:178
NORETURN void NoNumb(const char *chDesc) const
Definition parser.cpp:233
int GetQuote(char *chLabel, bool lgABORT)
Definition parser.h:209
int PrintLine(FILE *fp) const
Definition parser.h:204
static t_PredCont & Inst()
Definition cddefines.h:175
long add(double energy, const char *unit="Ryd")
Definition predcont.cpp:122
const realnum SMALLFLOAT
Definition cpu.h:191
t_geometry geometry
Definition geometry.cpp:5
t_input input
Definition input.cpp:12
t_iterations iterations
Definition iterations.cpp:5
t_optimize optimize
Definition optimize.cpp:5
void ParseStop(Parser &p)
t_phycon phycon
Definition phycon.cpp:6
UNUSED const double EVRYD
Definition physconst.h:189
UNUSED const double PARSEC
Definition physconst.h:138
t_radius radius
Definition radius.cpp:5
t_rfield rfield
Definition rfield.cpp:8
t_StopCalc StopCalc
Definition stopcalc.cpp:5
const int MXSTPL
Definition stopcalc.h:10