cloudy trunk
Loading...
Searching...
No Matches
cddrive.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/*cdDrive main routine to call cloudy under all circumstances) */
4/*cdReasonGeo write why the model stopped and type of geometry on io file */
5/*cdWarnings write all warnings entered into comment stack */
6/*cdEms obtain the local emissivity for a line, for the last computed zone */
7/*cdColm get the column density for a constituent */
8/*cdLine get the predicted line intensity, also index for line in stack */
9/*cdLine_ip get the predicted line intensity, using index for line in stack */
10/*cdCautions print out all cautions after calculation, on arbitrary io unit */
11/*cdTemp_last routine to query results and return temperature of last zone */
12/*cdDepth_depth get depth structure from previous iteration */
13/*cdTimescales returns thermal, recombination, and H2 formation timescales */
14/*cdSurprises print out all surprises on arbitrary unit number */
15/*cdNotes print stack of notes about current calculation */
16/*cdPressure_last routine to query results and return pressure of last zone */
17/*cdTalk tells the code whether to print results or be silent */
18/*cdOutput redirect output to arbitrary file */
19/*cdInput redirect input from arbitrary file */
20/*cdRead routine to read in command lines when cloudy used as subroutine */
21/*cdErrors produce summary of all warnings, cautions, etc, on arbitrary io unit */
22/*cdIonFrac get ionization fractions for a constituent */
23/*cdTemp get mean electron temperature for any element */
24/*cdCooling_last routine to query results and return cooling of last zone */
25/*cdHeating_last routine to query results and return heating of last zone */
26/*cdEDEN_last return electron density of last zone */
27/*cdNoExec call this routine to tell code not to actually execute */
28/*cdDate - puts date of code into string */
29/*cdVersion produces string that gives version number of the code */
30/*cdExecTime any routine can call this, find the time [s] since cdInit was called */
31/*cdPrintCommands( FILE *) prints all input commands into file */
32/*cdDrive main routine to call cloudy under all circumstances) */
33/*cdNwcns get the number of cautions and warnings, to tell if calculation is ok */
34/*debugLine provides a debugging hook into the main line array */
35/*cdEms_ip obtain the local emissivity for a line with known index */
36/*cdnZone gets number of zones */
37/*cdClosePunchFiles closes all the save files that have been used */
38/*cdB21cm - returns B as measured by 21 cm */
39/*cdPrtWL print line wavelengths in Angstroms in the standard format */
40#include "cddefines.h"
41#include "trace.h"
42#include "conv.h"
43#include "init.h"
44#include "lines.h"
45#include "pressure.h"
46#include "prt.h"
47#include "colden.h"
48#include "dense.h"
49#include "radius.h"
50#include "struc.h"
51#include "mole.h"
52#include "elementnames.h"
53#include "mean.h"
54#include "phycon.h"
55#include "called.h"
56#include "parse.h"
57#include "input.h"
58#include "taulines.h"
59#include "version.h"
60#include "thermal.h"
61#include "optimize.h"
62#include "grid.h"
63#include "timesc.h"
64#include "cloudy.h"
65#include "warnings.h"
66#include "lines_service.h"
67#include "cddrive.h"
68#include "iso.h"
69#include "save.h"
70
71/*************************************************************************
72 *
73 * cdDrive - main routine to call cloudy - returns 0 if all ok, 1 if problems
74 *
75 ************************************************************************/
76
78{
79 bool lgBAD;
80
81 DEBUG_ENTRY( "cdDrive()" );
82 /*********************************
83 * main routine to call cloudy *
84 *********************************/
85
86 /* this is set false when code loaded, set true when cdInit called,
87 * this is check that cdInit was called first */
88 if( !lgcdInitCalled )
89 {
90 printf(" cdInit was not called first - this must be the first call.\n");
92 }
93
94 if( trace.lgTrace )
95 {
96 fprintf( ioQQQ,
97 "cdDrive: lgOptimr=%1i lgVaryOn=%1i lgNoVary=%1i input.nSave:%li\n",
98 optimize.lgOptimr , optimize.lgVaryOn , optimize.lgNoVary, input.nSave );
99 }
100
101 /* should we call cloudy, or the optimization driver?
102 * possible to have VARY on line without OPTIMIZE being set
103 * lgNoVary set with "no optimize" command */
104 if( optimize.lgOptimr && optimize.lgVaryOn && !optimize.lgNoVary )
105 optimize.lgVaryOn = true;
106 else
107 optimize.lgVaryOn = false;
108
109 /* one time initialization of core load - returns if already called
110 * called here rather than in cdInit since at this point we know if
111 * single sim or grid */
112 InitCoreload();
113
114 if( optimize.lgVaryOn )
115 {
116 /* this branch called if optimizing or grid calculation */
117 if( trace.lgTrace )
118 fprintf( ioQQQ, "cdDrive: calling grid_do\n");
119 /* option to drive optimizer set if OPTIMIZE was in input stream */
120 lgBAD = grid_do();
121 }
122 else
123 {
124 if( trace.lgTrace )
125 fprintf( ioQQQ, "cdDrive: calling cloudy\n");
126
127 /* optimize did not occur, only compute one model, call cloudy */
128 lgBAD = cloudy();
129 }
130
131 /* reset flag saying that cdInit has not been called */
132 lgcdInitCalled = false;
133
134 if( lgAbort || lgBAD )
135 {
136 if( trace.lgTrace )
137 fprintf( ioQQQ, "cdDrive: returning failure during call. \n");
138 /* lgAbort set true if something wrong, so return lgBAD false. */
139 return 1;
140 }
141 else
142 {
143 /* everything is ok, return 0 */
144 return 0;
145 }
146}
147
148/*************************************************************************
149*
150* cdPrtWL write emission line wavelength
151*
152************************************************************************/
153
154/*cdPrtWL print line wavelengths in Angstroms in the standard format -
155 * a wrapper for prt_wl which must be kept parallel with sprt_wl
156 * both of those live in pdt.c */
157void cdPrtWL( FILE *io , realnum wl )
158{
159
160 DEBUG_ENTRY( "cdPrtWL()" );
161
162 prt_wl( io , wl );
163 return;
164}
165
166
167/*************************************************************************
168 *
169 * cdReasonGeo write why the model stopped and type of geometry on io file
170 *
171 ************************************************************************/
172
173
174/*cdReasonGeo write why the model stopped and type of geometry on io file */
175void cdReasonGeo(FILE * ioOUT)
176{
177
178 DEBUG_ENTRY( "cdReasonGeo()" );
179
180 /*this is the reason the calculation stopped*/
181 fprintf( ioOUT, "%s", warnings.chRgcln[0] );
182 fprintf( ioOUT , "\n" );
183 /* this is the geometry */
184 fprintf( ioOUT, "%s", warnings.chRgcln[1] );
185 fprintf( ioOUT , "\n" );
186 return;
187}
188
189
190/*************************************************************************
191 *
192 * cdWarnings write all warnings entered into comment stack
193 *
194 ************************************************************************/
195
196/*cdWarnings write all warnings entered into comment stack */
197
198void cdWarnings(FILE *ioPNT )
199{
200 long int i;
201
202 DEBUG_ENTRY( "cdWarnings()" );
203
204 if( warnings.nwarn > 0 )
205 {
206 for( i=0; i < warnings.nwarn; i++ )
207 {
208 /* these are all warnings that were entered in comment */
209 fprintf( ioPNT, "%s", warnings.chWarnln[i] );
210 fprintf( ioPNT, "\n" );
211 }
212 }
213
214 return;
215}
216
217
218/*************************************************************************
219 *
220 * cdCautions print out all cautions after calculation, on arbitrary io unit
221 *
222 ************************************************************************/
223
224/*cdCautions print out all cautions after calculation, on arbitrary io unit */
225
226void cdCautions(FILE * ioOUT)
227{
228 long int i;
229
230 DEBUG_ENTRY( "cdCautions()" );
231
232 if( warnings.ncaun > 0 )
233 {
234 for( i=0; i < warnings.ncaun; i++ )
235 {
236 fprintf( ioOUT, "%s", warnings.chCaunln[i] );
237 fprintf( ioOUT, "\n" );
238 }
239 }
240 return;
241}
242
243/*************************************************************************
244 *
245 * cdTimescales returns thermal, recombination, and H2 formation timescales
246 *
247 ************************************************************************/
248
250 /* the thermal cooling timescale */
251 double *TTherm ,
252 /* the hydrogen recombination timescale */
253 double *THRecom ,
254 /* the H2 formation timescale */
255 double *TH2 )
256{
257
258 DEBUG_ENTRY( "cdTimescales()" );
259
260 /* these were all evaluated in AgeCheck, which was called by PrtComment */
261
262 /* thermal or cooling timescale */
263 *TTherm = timesc.time_therm_long;
264
265 /* the hydrogen recombination timescale */
266 *THRecom = timesc.time_Hrecom_long;
267
268 /* longer of the the H2 formation and destruction timescales */
269 *TH2 = MAX2( timesc.time_H2_Dest_longest , timesc.time_H2_Form_longest );
270 return;
271}
272
273
274/*************************************************************************
275 *
276 * cdSurprises print out all surprises on arbitrary unit number
277 *
278 ************************************************************************/
279
280/*cdSurprises print out all surprises on arbitrary unit number */
281
282void cdSurprises(FILE * ioOUT)
283{
284 long int i;
285
286 DEBUG_ENTRY( "cdSurprises()" );
287
288 if( warnings.nbang > 0 )
289 {
290 for( i=0; i < warnings.nbang; i++ )
291 {
292 fprintf( ioOUT, "%s", warnings.chBangln[i] );
293 fprintf( ioOUT, "\n" );
294 }
295 }
296
297 return;
298}
299
300
301/*************************************************************************
302 *
303 * cdNotes print stack of notes about current calculation
304 *
305 ************************************************************************/
306
307/*cdNotes print stack of notes about current calculation */
308
309void cdNotes(FILE * ioOUT)
310{
311 long int i;
312
313 DEBUG_ENTRY( "cdNotes()" );
314
315 if( warnings.nnote > 0 )
316 {
317 for( i=0; i < warnings.nnote; i++ )
318 {
319 fprintf( ioOUT, "%s", warnings.chNoteln[i] );
320 fprintf( ioOUT, "\n" );
321 }
322 }
323 return;
324}
325
326/*************************************************************************
327 *
328 * cdCooling_last routine to query results and return cooling of last zone
329 *
330 ************************************************************************/
331
332/*cdCooling_last routine to query results and return cooling of last zone */
333double cdCooling_last() /* return cooling for last zone */
334{
335 return thermal.ctot;
336}
337
338/*************************************************************************
339 *
340 * cdVersion - puts version number of code into string
341 * incoming string must have sufficient length and will become null
342 * terminated string
343 *
344 ************************************************************************/
345
346void cdVersion(char chString[])
347{
348 strcpy( chString , t_version::Inst().chVersion );
349 return;
350}
351
352/*************************************************************************
353 *
354 * cdDate - puts date of code into string
355 * incoming string must have at least 8 char and will become null
356 * terminated string
357 *
358 ************************************************************************/
359
360/* cdDate - puts date of code into string */
361void cdDate(char chString[])
362{
363 strcpy( chString , t_version::Inst().chDate );
364 return;
365}
366
367
368/*************************************************************************
369 *
370 * cdHeating_last routine to query results and return heating of last zone
371 *
372 ************************************************************************/
373
374/*cdHeating_last routine to query results and return heating of last zone */
375
376double cdHeating_last() /* return heating for last zone */
377{
378 return thermal.htot;
379}
380
381
382/*************************************************************************
383 *
384 * cdEDEN_last return electron density of last zone
385 *
386 ************************************************************************/
387
388/*cdEDEN_last return electron density of last zone */
389
390double cdEDEN_last() /* return electron density for last zone */
391{
392 return dense.eden;
393}
394
395/*************************************************************************
396 *
397 * cdNoExec call this routine to tell code not to actually execute
398 *
399 ************************************************************************/
400
401/*cdNoExec call this routine to tell code not to actually execute */
402#include "noexec.h"
403
405{
406
407 DEBUG_ENTRY( "cdNoExec()" );
408
409 /* option to read in all input quantiles and NOT execute the actual model
410 * only check on input parameters - set by calling cdNoExec */
411 noexec.lgNoExec = true;
412
413 return;
414}
415
416
417/*************************************************************************
418 *
419 * cdSetExecTime routine to initialize variables keeping track of time at start of calculation
420 *
421 ************************************************************************/
422
423/* set to false initially, then to true when cdSetExecTime() is called to
424 * initialize the clock */
425static bool lgCalled=false;
426
427/* >>chng 06 dec 19, RP rm "|| defined(__HP_aCC)" to run on hp */
428#if defined(_MSC_VER)
429/* _MSC_VER branches assume getrusage isn't implemented by MS
430 * also is not implemented on our HP superdome */
431struct timeval {
432 long tv_sec; /* seconds */
433 long tv_usec; /* microseconds */
434};
435#else
436#include <sys/time.h>
437#include <sys/resource.h>
438#endif
439
440/* will be used to save initial time */
441static struct timeval before;
442
443/* cdClock stores time since arbitrary datum in clock_dat */
444STATIC void cdClock(struct timeval *clock_dat)
445{
446 DEBUG_ENTRY( "cdClock()" );
447
448/* >>chng 06 sep 2 rjrw: use long recurrence, fine grain UNIX clock *
449 * -- maintain system dependences in a single routine */
450#if defined(_MSC_VER) || defined(__HP_aCC)
451 clock_t clock_val;
452 /* >>chng 05 dec 21, from above to below due to negative exec times */
453 /*return (double)(clock() - before) / (double)CLOCKS_PER_SEC;*/
454 clock_val = clock();
455 clock_dat->tv_sec = clock_val/CLOCKS_PER_SEC;
456 clock_dat->tv_usec = 1000000*(clock_val-(clock_dat->tv_sec*CLOCKS_PER_SEC))/CLOCKS_PER_SEC;
457 /*>>chng 06 oct 05, this produced very large number, time typically 50% too long
458 clock_dat->tv_usec = 0;*/
459#else
460 struct rusage rusage;
461 if(getrusage(RUSAGE_SELF,&rusage) != 0)
462 {
463 fprintf( ioQQQ, "DISASTER cdClock called getrusage with invalid arguments.\n" );
464 fprintf( ioQQQ, "Sorry.\n" );
466 }
467 clock_dat->tv_sec = rusage.ru_utime.tv_sec;
468 clock_dat->tv_usec = rusage.ru_utime.tv_usec;
469#endif
470
471}
472/* cdSetExecTime is called by cdInit when everything is initialized,
473 * so that every time cdExecTime is called the elapsed time is returned */
475{
476 cdClock(&before);
477 lgCalled = true;
478}
479/* cdExecTime returns the elapsed time cpu time (sec) that has elapsed
480 * since cdInit called cdSetExecTime.*/
482{
483 DEBUG_ENTRY( "cdExecTime()" );
484
485 struct timeval clock_dat;
486 /* check that we were properly initialized */
487 if( lgCalled)
488 {
489 cdClock(&clock_dat);
490 /*fprintf(ioQQQ,"\n DEBUG sec %.2e usec %.2e\n",
491 (double)(clock_dat.tv_sec-before.tv_sec),
492 1e-6*(double)(clock_dat.tv_usec-before.tv_usec));*/
493 return (double)(clock_dat.tv_sec-before.tv_sec)+1e-6*(double)(clock_dat.tv_usec-before.tv_usec);
494 }
495 else
496 {
497 /* this is a big problem, we were asked for the elapsed time but
498 * the timer was not initialized by calling SetExecTime */
499 fprintf( ioQQQ, "DISASTER cdExecTime was called before SetExecTime, impossible.\n" );
500 fprintf( ioQQQ, "Sorry.\n" );
502 }
503}
504
505/*************************************************************************
506 *
507 * cdPrintCommands prints all input commands into file
508 *
509 ************************************************************************/
510
511/* cdPrintCommands( FILE *)
512 * prints all input commands into file */
513void cdPrintCommands( FILE * ioOUT )
514{
515 long int i;
516 fprintf( ioOUT, " Input commands follow:\n" );
517 fprintf( ioOUT, "c ======================\n" );
518
519 for( i=0; i <= input.nSave; i++ )
520 {
521 fprintf( ioOUT, "%s\n", input.chCardSav[i] );
522 }
523 fprintf( ioOUT, "c ======================\n" );
524}
525
526
527/*************************************************************************
528 *
529 * cdEms obtain the local emissivity for a line, for the last computed zone
530 *
531 ************************************************************************/
532// if called without last parameter, return intrinsic emission
533long int cdEmis(
534 /* return value will be index of line within stack,
535 * negative of number of lines in the stack if the line could not be found*/
536 /* 4 char null terminated string label */
537 char *chLabel,
538 /* line wavelength */
540 /* the vol emissivity of this line in last computed zone */
541 double *emiss )
542{
543 DEBUG_ENTRY( "cdEms()" );
544 long int i = cdEmis( chLabel , wavelength , emiss, false );
545 return i;
546}
547
548/* \todo 2 This routine, cdLine, cdEmis_ip, and cdLine_ip should be consolidated somewhat.
549 * in particular so that this routine has the same "closest line" reporting that cdLine has. */
550long int cdEmis(
551 /* return value will be index of line within stack,
552 * negative of number of lines in the stack if the line could not be found*/
553 /* 4 char null terminated string label */
554 char *chLabel,
555 /* line wavelength */
557 /* the vol emissivity of this line in last computed zone */
558 double *emiss ,
559 // false to return intrinsic emission, true for emergent
560 bool lgEmergent )
561{
562 /* use following to store local image of character strings */
563 char chCARD[INPUT_LINE_LENGTH];
564 char chCaps[5];
565 long int j;
566 realnum errorwave;
567
568 DEBUG_ENTRY( "cdEms()" );
569
570 /* returns the emissivity in the desired line
571 * used internally by the code, to do save lines structure */
572
573 // can return intrinsic or emergent emission, default is intrinsic
574 long int ipEmisType = 0;
575 if( lgEmergent )
576 ipEmisType = 1;
577
578 strcpy( chCARD, chLabel );
579
580 /* make sure chLabel is all caps */
581 caps(chCARD);/* convert to caps */
582
583 /* get the error associated with 4 significant figures */
584 errorwave = WavlenErrorGet( wavelength );
585
586 for( j=0; j < LineSave.nsum; j++ )
587 {
588 /* change chLabel to all caps to be like input chLineLabel */
589 cap4(chCaps, LineSv[j].chALab);
590
591 /* check wavelength and chLabel for a match */
592 /*if( fabs(LineSv[j].wavelength- wavelength)/MAX2(DELTA,wavelength)<errorwave &&
593 strcmp(chCaps,chCARD) == 0 ) */
594 if( fabs(LineSv[j].wavelength-wavelength) < errorwave && strcmp(chCaps,chCARD) == 0 )
595 {
596 /* match, so set emiss to emissivity in line */
597 *emiss = LineSv[j].emslin[ipEmisType];
598 /* and announce success by returning line index within stack */
599 return j;
600 }
601 }
602
603 /* we fell through without finding the line - return false */
604 return -LineSave.nsum;
605}
606
607/*************************************************************************
608 *
609 * cdEms_ip obtain the local emissivity for a line with known index
610 *
611 ************************************************************************/
612
614 /* index of the line in the stack */
615 long int ipLine,
616 /* the vol emissivity of this line in last computed zone */
617 double *emiss ,
618 // intrinsic or emergent
619 bool lgEmergent )
620{
621 DEBUG_ENTRY( "cdEmis_ip()" );
622
623 /* returns the emissivity in a line - implements save lines structure
624 * this version uses previously stored line index to save speed */
625 ASSERT( ipLine >= 0 && ipLine < LineSave.nsum );
626 *emiss = LineSv[ipLine].emslin[lgEmergent];
627 return;
628}
629
630/*************************************************************************
631 *
632 * cdColm get the column density for a constituent - 0 return if ok, 1 if problems
633 *
634 ************************************************************************/
635
637 /* return value is zero if all ok, 1 if errors happened */
638 /* 4-char + eol string that is first
639 * 4 char of element name as spelled by Cloudy, upper or lower case */
640 const char *chLabel,
641
642 /* integer stage of ionization, 1 for atom, 2 for A+, etc,
643 * 0 is special flag for CO, H2, OH, or excited state */
644 long int ion,
645
646 /* the column density derived by the code [cm-2] */
647 double *theocl )
648{
649 long int nelem;
650 /* use following to store local image of character strings */
651 char chLABEL_CAPS[20];
652
653 DEBUG_ENTRY( "cdColm()" );
654
655 /* check that chLabel[4] is null - supposed to be 4 char + end */
656 if( chLabel[4] != '\0' || strlen(chLabel) != 4 )
657 {
658 fprintf( ioQQQ, " cdColm called with insane chLabel (between quotes) \"%s\", must be exactly 4 characters long.\n",
659 chLabel );
660 return 1;
661 }
662
663 strcpy( chLABEL_CAPS, chLabel );
664 /* convert element label to all caps */
665 caps(chLABEL_CAPS);
666
667 /* zero ionization stage has special meaning. The quantities recognized are
668 * the molecules, "H2 ", "OH ", "CO ", etc
669 * "CII*" excited state C+ */
670 if( ion < 0 )
671 {
672 fprintf( ioQQQ, " cdColm called with insane ion, =%li\n",
673 ion );
674 return 1;
675 }
676
677 else if( ion == 0 )
678 {
679 /* this case molecular column density */
680 /* want the molecular hydrogen H2 column density */
681 if( strcmp( chLABEL_CAPS , "H2 " )==0 )
682 {
683 *theocl = colden.colden[ipCOL_H2g] + colden.colden[ipCOL_H2s];
684 }
685
686 /* H- column density */
687 else if( strcmp( chLABEL_CAPS , "H- " )==0 )
688 {
689 *theocl = colden.colden[ipCOL_HMIN];
690 }
691
692 /* H2+ column density ipCOL_H2p is 4 */
693 else if( strcmp( chLABEL_CAPS , "H2+ " )==0 )
694 {
695 *theocl = colden.colden[ipCOL_H2p];
696 }
697
698 /* H3+ column density */
699 else if( strcmp( chLABEL_CAPS , "H3+ " )==0 )
700 {
701 *theocl = colden.colden[ipCOL_H3p];
702 }
703
704 /* H2g - ground H2 column density */
705 else if( strcmp( chLABEL_CAPS , "H2G " )==0 )
706 {
707 *theocl = colden.colden[ipCOL_H2g];
708 }
709
710 /* H2* - excited H2 - column density */
711 else if( strcmp( chLABEL_CAPS , "H2* " )==0 )
712 {
713 *theocl = colden.colden[ipCOL_H2s];
714 }
715
716 /* HeH+ column density */
717 else if( strcmp( chLABEL_CAPS , "HEH+" )==0 )
718 {
719 *theocl = colden.colden[ipCOL_HeHp];
720 }
721
722 /* carbon monoxide column density */
723 else if( strcmp( chLABEL_CAPS , "CO " )==0 )
724 {
725 *theocl = findspecieslocal("CO")->column;
726 }
727
728 /* OH column density */
729 else if( strcmp( chLABEL_CAPS , "OH " )==0 )
730 {
731 *theocl = findspecieslocal("OH")->column;
732 }
733
734 /* H2O column density */
735 else if( strcmp( chLABEL_CAPS , "H2O " )==0 )
736 {
737 *theocl = findspecieslocal("H2O")->column;
738 }
739
740 /* O2 column density */
741 else if( strcmp( chLABEL_CAPS , "O2 " )==0 )
742 {
743 *theocl = findspecieslocal("O2")->column;
744 }
745
746 /* SiO column density */
747 else if( strcmp( chLABEL_CAPS , "SIO " )==0 )
748 {
749 *theocl = findspecieslocal("SiO")->column;
750 }
751
752 /* C2 column density */
753 else if( strcmp( chLABEL_CAPS , "C2 " )==0 )
754 {
755 *theocl = findspecieslocal("C2")->column;
756 }
757
758 /* C3 column density */
759 else if( strcmp( chLABEL_CAPS , "C3 " )==0 )
760 {
761 *theocl = findspecieslocal("C3")->column;
762 }
763
764 /* CN column density */
765 else if( strcmp( chLABEL_CAPS , "CN " )==0 )
766 {
767 *theocl = findspecieslocal("CN")->column;
768 }
769
770 /* CH column density */
771 else if( strcmp( chLABEL_CAPS , "CH " )==0 )
772 {
773 *theocl = findspecieslocal("CH")->column;
774 }
775
776 /* CH+ column density */
777 else if( strcmp( chLABEL_CAPS , "CH+ " )==0 )
778 {
779 *theocl = findspecieslocal("CH+")->column;
780 }
781
782 /* ===========================================================*/
783 /* end special case molecular column densities, start special cases
784 * excited state column densities */
785 /* CII^* column density, population of J=3/2 upper level of split ground term */
786 else if( strcmp( chLABEL_CAPS , "CII*" )==0 )
787 {
788 *theocl = colden.C2Colden[1];
789 }
790 else if( strcmp( chLABEL_CAPS , "C11*" )==0 )
791 {
792 *theocl = colden.C1Colden[0];
793 }
794 else if( strcmp( chLABEL_CAPS , "C12*" )==0 )
795 {
796 *theocl = colden.C1Colden[1];
797 }
798 else if( strcmp( chLABEL_CAPS , "C13*" )==0 )
799 {
800 *theocl = colden.C1Colden[2];
801 }
802 else if( strcmp( chLABEL_CAPS , "O11*" )==0 )
803 {
804 *theocl = colden.O1Colden[0];
805 }
806 else if( strcmp( chLABEL_CAPS , "O12*" )==0 )
807 {
808 *theocl = colden.O1Colden[1];
809 }
810 else if( strcmp( chLABEL_CAPS , "O13*" )==0 )
811 {
812 *theocl = colden.O1Colden[2];
813 }
814 /* CIII excited states, upper level of 1909 */
815 else if( strcmp( chLABEL_CAPS , "C30*" )==0 )
816 {
817 *theocl = colden.C3Colden[1];
818 }
819 else if( strcmp( chLABEL_CAPS , "C31*" )==0 )
820 {
821 *theocl = colden.C3Colden[2];
822 }
823 else if( strcmp( chLABEL_CAPS , "C32*" )==0 )
824 {
825 *theocl = colden.C3Colden[3];
826 }
827 else if( strcmp( chLABEL_CAPS , "SI2*" )==0 )
828 {
829 *theocl = colden.Si2Colden[1];
830 }
831 else if( strcmp( chLABEL_CAPS , "HE1*" )==0 )
832 {
833 *theocl = colden.He123S;
834 }
835 /* special option, "H2vJ" */
836 else if( strncmp(chLABEL_CAPS , "H2" , 2 ) == 0 )
837 {
838 long int iVib = chLABEL_CAPS[2] - '0';
839 long int iRot = chLABEL_CAPS[3] - '0';
840 if( iVib<0 || iRot < 0 )
841 {
842 fprintf( ioQQQ, " cdColm called with insane v,J for H2=\"%4.4s\" caps=\"%4.4s\"\n",
843 chLabel , chLABEL_CAPS );
844 return 1;
845 }
846 *theocl = cdH2_colden( iVib , iRot );
847 }
848
849 /* clueless as to what was meant - bail */
850 else
851 {
852 fprintf( ioQQQ, " cdColm called with unknown element chLabel, org=\"%4.4s \" 0 caps=\"%4.4s\" 0\n",
853 chLabel , chLABEL_CAPS );
854 return 1;
855 }
856 }
857 else
858 {
859 /* this case, ionization stage of some element */
860 /* find which element this is */
861 nelem = 0;
862 while( nelem < LIMELM &&
863 strncmp(chLABEL_CAPS,elementnames.chElementNameShort[nelem],4) != 0 )
864 {
865 ++nelem;
866 }
867
868 /* this is true if we have one of the first 30 elements in the label,
869 * nelem is on C scale */
870 if( nelem < LIMELM )
871 {
872
873 /* sanity check - does this ionization stage exist?
874 * max2 is to pick up H2 as H 3 */
875 if( ion > MAX2(3,nelem + 2) )
876 {
877 fprintf( ioQQQ,
878 " cdColm asked to return ionization stage %ld for element %s but this is not physical.\n",
879 ion, chLabel );
880 return 1;
881 }
882
883 /* the column density, ion is on physics scale, but means are on C scale */
884 *theocl = mean.xIonMean[0][nelem][ion-1][0];
885 /*>>chng 06 jan 23, div by factor of two
886 * special case of H2 when being tricked as H 3 - this stores 2H_2 so that
887 * the fraction of H in H0 and H+ is correct - need to remove this extra
888 * factor of two here */
889 if( nelem==ipHYDROGEN && ion==3 )
890 *theocl /= 2.;
891 }
892 else
893 {
894 fprintf( ioQQQ,
895 " cdColm did not understand this combination of ion %4ld and element %4.4s.\n",
896 ion, chLabel );
897 return 1;
898 }
899 }
900 return 0;
901}
902
903
904/*************************************************************************
905 *
906 *cdErrors produce summary of all warnings, cautions, etc, on arbitrary io unit
907 *
908 ************************************************************************/
909
910void cdErrors(FILE *ioOUT)
911{
912 long int nc,
913 nn,
914 npe,
915 ns,
916 nte,
917 nw ,
918 nIone,
919 nEdene;
920 bool lgAbort_loc;
921
922 DEBUG_ENTRY( "cdErrors()" );
923
924 /* first check for number of warnings, cautions, etc */
925 cdNwcns(&lgAbort_loc,&nw,&nc,&nn,&ns,&nte,&npe, &nIone, &nEdene );
926
927 /* only say something is one of these problems is nonzero */
928 if( nw || nc || nte || npe || nIone || nEdene || lgAbort_loc )
929 {
930 /* say the title of the model */
931 fprintf( ioOUT, "%75.75s\n", input.chTitle );
932
933 if( lgAbort_loc )
934 fprintf(ioOUT," Calculation ended with abort!\n");
935
936 /* print warnings on the io unit */
937 if( nw != 0 )
938 {
939 cdWarnings(ioOUT);
940 }
941
942 /* print cautions on the io unit */
943 if( nc != 0 )
944 {
945 cdCautions(ioOUT);
946 }
947
948 if( nte != 0 )
949 {
950 fprintf( ioOUT , "Te failures=%4ld\n", nte );
951 }
952
953 if( npe != 0 )
954 {
955 fprintf( ioOUT , "Pressure failures=%4ld\n", npe );
956 }
957
958 if( nIone != 0 )
959 {
960 fprintf( ioOUT , "Ionization failures=%4ld\n", nte );
961 }
962
963 if( nEdene != 0 )
964 {
965 fprintf( ioOUT , "Electron density failures=%4ld\n", npe );
966 }
967 }
968 return;
969}
970
971/*************************************************************************
972 *
973 *cdDepth_depth get depth structure from previous iteration
974 *
975 ************************************************************************/
976void cdDepth_depth( double cdDepth[] )
977{
978 long int nz;
979
980 DEBUG_ENTRY( "cdDepth_depth()" );
981
982 for( nz = 0; nz<nzone; ++nz )
983 {
984 cdDepth[nz] = struc.depth[nz];
985 }
986 return;
987}
988
989/*************************************************************************
990 *
991 *cdPressure_depth routine to query results and return pressure of last iteration
992 *
993 ************************************************************************/
994
995/*
996 * cdPressure_depth
997 * This returns the pressure and its constituents for the last iteration.
998 * space was allocated in the calling routine for the vectors -
999 * before calling this, cdnZone should have been called to get the number of
1000 * zones, then space allocated for the arrays */
1002 /* total pressure, all forms*/
1003 double TotalPressure[],
1004 /* gas pressure */
1005 double GasPressure[],
1006 /* line radiation pressure */
1007 double RadiationPressure[])
1008{
1009 long int nz;
1010
1011 DEBUG_ENTRY( "cdPressure_depth()" );
1012
1013 for( nz = 0; nz<nzone; ++nz )
1014 {
1015 TotalPressure[nz] = struc.pressure[nz];
1016 GasPressure[nz] = struc.GasPressure[nz];
1017 RadiationPressure[nz] = struc.pres_radiation_lines_curr[nz];
1018 }
1019 return;
1020}
1021
1022/*************************************************************************
1023 *
1024 *cdPressure_last routine to query results and return pressure of last zone
1025 *
1026 ************************************************************************/
1027
1029 double *PresTotal, /* total pressure, all forms, for the last computed zone*/
1030 double *PresGas, /* gas pressure */
1031 double *PresRad) /* line radiation pressure */
1032{
1033
1034 DEBUG_ENTRY( "cdPressure_last()" );
1035
1036 *PresGas = pressure.PresGasCurr;
1037 *PresRad = pressure.pres_radiation_lines_curr;
1038 *PresTotal = pressure.PresTotlCurr;
1039 return;
1040}
1041
1042/*************************************************************************
1043 *
1044 *cdnZone gets number of zones
1045 *
1046 ************************************************************************/
1047
1048/* returns number of zones */
1049long int cdnZone()
1050{
1051 return nzone;
1052}
1053
1054/*************************************************************************
1055 *
1056 *cdTemp_last routine to query results and return temperature of last zone
1057 *
1058 ************************************************************************/
1059
1060
1062{
1063 return phycon.te;
1064}
1065
1066/*************************************************************************
1067 *
1068 *cdIonFrac get ionization fractions for a constituent
1069 *
1070 ************************************************************************/
1071
1073 /* four char string, null terminated, giving the element name */
1074 const char *chLabel,
1075 /* IonStage is ionization stage, 1 for atom, up to N+1 where N is atomic number,
1076 * 0 says special case */
1077 long int IonStage,
1078 /* will be fractional ionization */
1079 double *fracin,
1080 /* how to weight the average, must be "VOLUME" or "RADIUS" */
1081 const char *chWeight ,
1082 /* if true then weighting also has electron density, if false then only volume or radius */
1083 bool lgDensity )
1084 /* return value is 0 if element was found, non-zero if failed */
1085{
1086 long int ip,
1087 ion, /* used as index within aaa vector*/
1088 nelem;
1089 realnum aaa[LIMELM + 1];
1090 /* use following to store local image of character strings */
1091 char chCARD[INPUT_LINE_LENGTH];
1092
1093 DEBUG_ENTRY( "cdIonFrac()" );
1094
1095 strcpy( chCARD, chWeight );
1096 /* make sure chWeight is all caps */
1097 caps(chCARD);/* convert to caps */
1098
1099 int dim;
1100 if( strcmp(chCARD,"RADIUS") == 0 )
1101 dim = 0;
1102 else if( strcmp(chCARD,"AREA") == 0 )
1103 dim = 1;
1104 else if( strcmp(chCARD,"VOLUME") == 0 )
1105 dim = 2;
1106 else
1107 {
1108 fprintf( ioQQQ, " cdIonFrac: chWeight=%6.6s makes no sense to me, valid options are RADIUS, AREA, and VOLUME\n",
1109 chWeight );
1110 *fracin = 0.;
1111 return 1;
1112 }
1113
1114 /* first ensure that chLabel is all caps */
1115 strcpy( chCARD, chLabel );
1116 /* make sure chLabel is all caps */
1117 caps(chCARD);/* convert to caps */
1118
1119 if( IonStage==0 )
1120 {
1121 /* special case */
1122 if( strcmp(chCARD,"H2 " ) == 0 )
1123 {
1124 /* this will be request for H2, treated as third stage of hydrogen */
1125 nelem = 0;
1126 IonStage = 3;
1127 }
1128 else
1129 {
1130 fprintf( ioQQQ, " cdIonFrac: ion stage of zero and element %s makes no sense to me\n",
1131 chCARD );
1132 *fracin = 0.;
1133 return 1;
1134 }
1135 }
1136
1137 else
1138 {
1139 /* find which element this is */
1140 nelem = 0;
1141 while( nelem < LIMELM &&
1142 strcmp(chCARD,elementnames.chElementNameShort[nelem]) != 0 )
1143 {
1144 ++nelem;
1145 }
1146 }
1147
1148 /* if element not recognized and above loop falls through, nelem is LIMELM */
1149 if( nelem >= LIMELM )
1150 {
1151 fprintf( ioQQQ, " cdIonFrac called with unknown element chLabel, =%4.4s\n",
1152 chLabel );
1153 return 1;
1154 }
1155
1156 /* sanity check - does this ionization stage exist?
1157 * IonStage is on spectroscopic scale and nelem is on C scale */
1158
1159 /* ion will be used as pointer within the aaa array that contains average values,
1160 * convert to C scale */
1161 ion = IonStage - 1;
1162
1163 if( (ion > nelem+1 || ion < 0 ) && !(nelem==ipHYDROGEN&&ion==2))
1164 {
1165 fprintf( ioQQQ, " cdIonFrac asked to return ionization stage %ld for element %4.4s but this is not physical.\n",
1166 IonStage, chLabel );
1167 *fracin = -1.;
1168 return 1;
1169 }
1170
1171 /* get average, aaa is filled in from 0 */
1172 /* aaa is dim'd LIMELM+1 so largest argument is LIMELM
1173 * 'i' means ionization, not temperature */
1174 /* last argument says whether to include electron density */
1175 /* MeanIon uses c scale for nelem */
1176 mean.MeanIon('i',nelem,dim,&ip,aaa,lgDensity);
1177 *fracin = pow((realnum)10.f,aaa[ion]);
1178
1179 /* we succeeded - say so */
1180 return 0;
1181}
1182
1183/*************************************************************************
1184 *
1185 * debugLine provides a debugging hook into the main line array
1186 *
1187 ************************************************************************/
1188
1189 /* debugLine provides a debugging hook into the main line array
1190 * loops over whole array and finds every line that matches length,
1191 * the wavelength, the argument to the function
1192 * put breakpoint inside if test
1193 * return value is number of matches, also prints all matches*/
1195{
1196 long j, kount;
1197 realnum errorwave;
1198
1199 kount = 0;
1200
1201 /* get the error associated with 4 significant figures */
1202 errorwave = WavlenErrorGet( wavelength );
1203
1204 for( j=0; j < LineSave.nsum; j++ )
1205 {
1206 /* check wavelength and chLabel for a match */
1207 /* if( fabs(LineSv[j].wavelength- wavelength)/MAX2(DELTA,wavelength) < errorwave ) */
1208 if( fabs(LineSv[j].wavelength-wavelength) < errorwave )
1209 {
1210 printf("%s\n", LineSv[j].chALab);
1211 ++kount;
1212 }
1213 }
1214 printf(" hits = %li\n", kount );
1215 return kount;
1216}
1217
1218/*************************************************************************
1219 *
1220 *cdLine get the predicted line intensity, also index for line in stack
1221 *
1222 ************************************************************************/
1223
1224// returns array index for line in array stack if we found the line,
1225// return negative of total number of lines as debugging aid if line not found
1226// return <0 if problems
1227// emergent or intrinsic not specified - use intrinsic
1228long int cdLine(
1229 const char *chLabel,
1230 /* wavelength of line in angstroms, not format printed by code */
1232 /* linear intensity relative to normalization line*/
1233 double *relint,
1234 /* log of luminosity or intensity of line */
1235 double *absint )
1236{
1237 DEBUG_ENTRY( "cdLine()" );
1238 long int i = cdLine( chLabel , wavelength , relint , absint, 0 );
1239 return i;
1240}
1241long int cdLine(
1242 const char *chLabel,
1243 /* wavelength of line in angstroms, not format printed by code */
1245 /* linear intensity relative to normalization line*/
1246 double *relint,
1247 /* log of luminosity or intensity of line */
1248 double *absint ,
1249 // 0 is intrinsic,
1250 // 1 emergent
1251 // 2 is intrinsic cumulative,
1252 // 3 emergent cumulative
1253 int LineType )
1254{
1255 char chCaps[5],
1256 chFind[5];
1257 long int ipobs,
1258 j, index_of_closest=LONG_MIN,
1259 index_of_closest_w_correct_label=-1;
1260 realnum errorwave, smallest_error=BIGFLOAT,
1261 smallest_error_w_correct_label=BIGFLOAT;
1262
1263 DEBUG_ENTRY( "cdLine()" );
1264
1265 if( LineType<0 || LineType>3 )
1266 {
1267 fprintf( ioQQQ, " cdLine called with insane nLineType - it must be between 0 and 3.\n");
1268 return 0;
1269 }
1270
1271 /* this is zero when cdLine called with cdNoExec called too
1272 * will be zero if code aborts during search for initial conditions */
1273 if( LineSave.nsum == 0 )
1274 {
1275 *relint = 0.;
1276 *absint = 0.;
1277 return 0;
1278 }
1279 ASSERT( LineSave.ipNormWavL >= 0);
1280 ASSERT( LineSave.nsum > 0);
1281
1282 /* check that chLabel[4] is null - supposed to be 4 char + end */
1283 if( chLabel[4] != '\0' || strlen(chLabel) != 4 )
1284 {
1285 fprintf( ioQQQ, " cdLine called with insane chLabel (between quotes) \"%s\", must be exactly 4 characters long.\n",
1286 chLabel );
1287 return 1;
1288 }
1289
1290 /* change chLabel to all caps */
1291 cap4(chFind, chLabel);
1292
1293 /* this replaces tabs with spaces. */
1294 /* \todo 2 keep this in, do it elsewhere, just warn and bail? */
1295 for( j=0; j<=3; j++ )
1296 {
1297 if( chFind[j] == '\t' )
1298 {
1299 chFind[j] = ' ';
1300 }
1301 }
1302
1303 /* get the error associated with 4 significant figures */
1304 errorwave = WavlenErrorGet( wavelength );
1305
1306 {
1307 enum{DEBUG_LOC=false};
1308 if( DEBUG_LOC && fabs(wavelength-1000.)<0.01 )
1309 {
1310 fprintf(ioQQQ,"cdDrive wl %.4e error %.3e\n",
1311 wavelength, errorwave );
1312 }
1313 }
1314
1315 /* now go through entire line stack, do not do 0, which is unity integration */
1316 for( j=1; j < LineSave.nsum; j++ )
1317 {
1318 /* find closest line to requested wavelength to
1319 * report when we don't get exact match */
1320 realnum current_error;
1321 current_error = (realnum)fabs(LineSv[j].wavelength-wavelength);
1322
1323 /* change chLabel to all caps to be like input chALab */
1324 cap4(chCaps, LineSv[j].chALab);
1325
1326 if( current_error < smallest_error )
1327 {
1328 index_of_closest = j;
1329 smallest_error = current_error;
1330 }
1331
1332 if( current_error < smallest_error_w_correct_label &&
1333 (strcmp(chCaps,chFind) == 0) )
1334 {
1335 index_of_closest_w_correct_label = j;
1336 smallest_error_w_correct_label = current_error;
1337 }
1338
1339 /* check wavelength and chLabel for a match */
1340 /* DELTA since often want wavelength of zero */
1341 if( current_error <= errorwave ||
1342 fp_equal( wavelength + errorwave, LineSv[j].wavelength ) ||
1343 fp_equal( wavelength - errorwave, LineSv[j].wavelength ) )
1344 {
1345 /* now see if labels agree */
1346 if( strcmp(chCaps,chFind) == 0 )
1347 {
1348 /* match, so set pointer */
1349 ipobs = j;
1350
1351 /* does the normalization line have a positive intensity*/
1352 if( LineSv[LineSave.ipNormWavL].SumLine[LineType] > 0. )
1353 {
1354 *relint = LineSv[ipobs].SumLine[LineType]/
1355 LineSv[LineSave.ipNormWavL].SumLine[LineType]*
1356 LineSave.ScaleNormLine;
1357 }
1358 else
1359 {
1360 *relint = 0.;
1361 }
1362
1363 /* return log of current line intensity if it is positive */
1364 if( LineSv[ipobs].SumLine[LineType] > 0. )
1365 {
1366 *absint = log10(LineSv[ipobs].SumLine[LineType]) +
1367 radius.Conv2PrtInten;
1368 }
1369 else
1370 {
1371 /* line intensity is actually zero, return small number */
1372 *absint = -37.;
1373 }
1374 /* we found the line, return pointer to its location */
1375 return ipobs;
1376 }
1377 }
1378 }
1379
1380 /* >>chng 05 dec 21, report closest line if we did not find exact match, note that
1381 * exact match returns above, where we will return negative number of lines in stack */
1382 fprintf( ioQQQ," PROBLEM cdLine did not find line with label (between quotes) \"%4s\" and wavelength ", chFind );
1384 if( index_of_closest >= 0 )
1385 {
1386 fprintf( ioQQQ,".\n The closest line (any label) was \"%4s\"\t",
1387 LineSv[index_of_closest].chALab );
1388 prt_wl(ioQQQ,LineSv[index_of_closest].wavelength );
1389 if( index_of_closest_w_correct_label >= 0 )
1390 {
1391 fprintf( ioQQQ,"\n The closest with correct label was \"%4s\"\t", chFind );
1392 prt_wl(ioQQQ,LineSv[index_of_closest_w_correct_label].wavelength );
1393 fprintf( ioQQQ,"\n" );
1394 }
1395 else
1396 fprintf( ioQQQ,"\n No line found with label \"%s\".\n", chFind );
1397 fprintf( ioQQQ,"\n" );
1398 }
1399 else
1400 {
1401 fprintf( ioQQQ,".\n PROBLEM No close line was found\n" );
1402 TotalInsanity();
1403 }
1404
1405 *absint = 0.;
1406 *relint = 0.;
1407
1408 /* if we fell down to here we did not find the line
1409 * return negative of total number of lines as debugging aid */
1410 return -LineSave.nsum;
1411}
1412
1413/*cdLine_ip get the predicted line intensity, using index for line in stack */
1414void cdLine_ip(long int ipLine,
1415 /* linear intensity relative to normalization line*/
1416 double *relint,
1417 /* log of luminosity or intensity of line */
1418 double *absint )
1419{
1420
1421 DEBUG_ENTRY( "cdLine_ip()" );
1422 cdLine_ip( ipLine , relint , absint , 0 );
1423}
1424void cdLine_ip(long int ipLine,
1425 /* linear intensity relative to normalization line*/
1426 double *relint,
1427 /* log of luminosity or intensity of line */
1428 double *absint ,
1429 // 0 is intrinsic,
1430 // 1 emergent
1431 // 2 is intrinsic cumulative,
1432 // 3 emergent cumulative
1433 int LineType )
1434{
1435
1436 DEBUG_ENTRY( "cdLine_ip()" );
1437
1438 if( LineType<0 || LineType>3 )
1439 {
1440 fprintf( ioQQQ, " cdLine_ip called with insane nLineType - it must be between 0 and 3.\n");
1441 *relint = 0.;
1442 *absint = 0.;
1443 return;
1444 }
1445
1446 /* this is zero when cdLine called with cdNoExec called too */
1447 if( LineSave.nsum == 0 )
1448 {
1449 *relint = 0.;
1450 *absint = 0.;
1451 return;
1452 }
1453 ASSERT( LineSave.ipNormWavL >= 0);
1454 ASSERT( LineSave.nsum > 0);
1455
1456 /* does the normalization line have a positive intensity*/
1457 if( LineSv[LineSave.ipNormWavL].SumLine[LineType] > 0. )
1458 *relint = LineSv[ipLine].SumLine[LineType]/
1459 LineSv[LineSave.ipNormWavL].SumLine[LineType]*
1460 LineSave.ScaleNormLine;
1461 else
1462 *relint = 0.;
1463
1464 /* return log of current line intensity if it is positive */
1465 if( LineSv[ipLine].SumLine[LineType] > 0. )
1466 *absint = log10(LineSv[ipLine].SumLine[LineType]) +
1467 radius.Conv2PrtInten;
1468 else
1469 /* line intensity is actually zero, return small number */
1470 *absint = -37.;
1471
1472 return;
1473}
1474
1475/*************************************************************************
1476 *
1477 *cdNwcns get the number of cautions and warnings, to tell if calculation is ok
1478 *
1479 ************************************************************************/
1480
1482 /* abort status, this better be false */
1483 bool *lgAbort_ret ,
1484 /* the number of warnings, cautions, notes, and surprises */
1485 long int *NumberWarnings,
1486 long int *NumberCautions,
1487 long int *NumberNotes,
1488 long int *NumberSurprises,
1489 /* the number of temperature convergence failures */
1490 long int *NumberTempFailures,
1491 /* the number of pressure convergence failures */
1492 long int *NumberPresFailures,
1493 /* the number of ionization convergence failures */
1494 long int *NumberIonFailures,
1495 /* the number of electron density convergence failures */
1496 long int *NumberNeFailures )
1497{
1498
1499 DEBUG_ENTRY( "cdNwcns()" );
1500
1501 /* this would be set true if code ended with abort - very very bad */
1502 *lgAbort_ret = lgAbort;
1503 /* this sub is called after comment, to find the number of various comments */
1504 *NumberWarnings = warnings.nwarn;
1505 *NumberCautions = warnings.ncaun;
1506 *NumberNotes = warnings.nnote;
1507 *NumberSurprises = warnings.nbang;
1508
1509 /* these are counters that were incremented during convergence failures */
1510 *NumberTempFailures = conv.nTeFail;
1511 *NumberPresFailures = conv.nPreFail;
1512 *NumberIonFailures = conv.nIonFail;
1513 *NumberNeFailures = conv.nNeFail;
1514 return;
1515}
1516
1517void cdOutput( const char* filename, const char *mode )
1518{
1519 DEBUG_ENTRY( "cdOutput()" );
1520
1521 if( ioQQQ != stdout && ioQQQ != NULL )
1522 fclose(ioQQQ);
1523 FILE *fp = stdout;
1524 if( *filename != '\0' )
1525 fp = open_data( filename, mode, AS_LOCAL_ONLY );
1526 save.chOutputFile = filename;
1527 ioQQQ = fp;
1528}
1529
1530void cdInput( const char* filename, const char *mode )
1531{
1532 DEBUG_ENTRY( "cdInput()" );
1533
1534 if( ioStdin != stdin && ioStdin != NULL )
1535 fclose(ioStdin);
1536 FILE *fp = stdin;
1537 if( *filename != '\0' )
1538 fp = open_data( filename, mode, AS_LOCAL_ONLY );
1539 ioStdin = fp;
1540}
1541
1542/*************************************************************************
1543 *
1544 * cdTalk tells the code whether to print results or be silent
1545 *
1546 ************************************************************************/
1547
1548void cdTalk(bool lgTOn)
1549{
1550
1551 DEBUG_ENTRY( "cdTalk()" );
1552
1553 /* MPI talking rules must be obeyed, otherwise loss of output may result */
1554 called.lgTalk = lgTOn && cpu.i().lgMPI_talk();
1555 /* has talk been forced off? */
1556 called.lgTalkForcedOff = !lgTOn;
1557 return;
1558}
1559
1560/*cdB21cm - returns B as measured by 21 cm */
1561double cdB21cm()
1562{
1563 double ret;
1564
1565 DEBUG_ENTRY( "cdB21cm()" );
1566
1567 // return average over radius
1568 if( mean.TempB_HarMean[0][1] > SMALLFLOAT )
1569 {
1570 ret = mean.TempB_HarMean[0][0]/mean.TempB_HarMean[0][1];
1571 }
1572 else
1573 {
1574 ret = 0.;
1575 }
1576 return ret;
1577}
1578
1579/*************************************************************************
1580 *
1581 * cdTemp get mean electron temperature for any element
1582 *
1583 ************************************************************************/
1584
1585/* This routine finds the mean electron temperature for any ionization stage
1586 * It returns 0 if it could find the species, 1 if it could not find the species.
1587 * The first argument is a null terminated 4 char string that gives the element
1588 * name as spelled by Cloudy.
1589 * The second argument is ion stage, 1 for atom, 2 for first ion, etc
1590 * This third argument will be returned as result,
1591 * Last parameter is either "RADIUS", "AREA", or "VOLUME" to give weighting
1592 *
1593 * if the ion stage is zero then the element label will have a special meaning.
1594 * The string "21CM" is will return the 21 cm temperature.
1595 * The string "H2 " will return the temperature weighted wrt the H2 density
1596 * There are several other options as listed below */
1599
1600/* return value is o if things are ok and element was found,
1601 * non-zero if element not found or there are problems */
1603 /* four char string, null terminated, giving the element name */
1604 const char *chLabel,
1605 /* IonStage is ionization stage, 1 for atom, up to Z+1 where Z is atomic number,
1606 * 0 means that chLabel is a special case */
1607 long int IonStage,
1608 /* will be temperature */
1609 double *TeMean,
1610 /* how to weight the average, must be "RADIUS", "AREA, or "VOLUME" */
1611 const char *chWeight )
1612{
1613 long int ip,
1614 ion, /* used as pointer within aaa vector*/
1615 nelem;
1616 realnum aaa[LIMELM + 1];
1617 /* use following to store local image of character strings */
1618 char chWGHT[INPUT_LINE_LENGTH] , chELEM[INPUT_LINE_LENGTH];
1619
1620 DEBUG_ENTRY( "cdTemp()" );
1621
1622 /* make sure strings are all caps */
1623 strcpy( chWGHT, chWeight );
1624 caps(chWGHT);
1625 strcpy( chELEM, chLabel );
1626 caps(chELEM);
1627
1628 /* now see which weighting */
1629 int dim;
1630 if( strcmp(chWGHT,"RADIUS") == 0 )
1631 dim = 0;
1632 else if( strcmp(chWGHT,"AREA") == 0 )
1633 dim = 1;
1634 else if( strcmp(chWGHT,"VOLUME") == 0 )
1635 dim = 2;
1636 else
1637 {
1638 fprintf( ioQQQ, " cdTemp: chWeight=%6.6s makes no sense to me, the options are RADIUS, AREA, and VOLUME.\n",
1639 chWeight );
1640 *TeMean = 0.;
1641 return 1;
1642 }
1643
1644 if( IonStage == 0 )
1645 {
1646 /* return atomic hydrogen weighted harmonic mean of gas kinetic temperature */
1647 if( strcmp(chELEM,"21CM") == 0 )
1648 {
1649 if( mean.TempHarMean[dim][1] > SMALLFLOAT )
1650 *TeMean = mean.TempHarMean[dim][0]/mean.TempHarMean[dim][1];
1651 else
1652 *TeMean = 0.;
1653 }
1654 /* return atomic hydrogen weighted harmonic mean 21 cm spin temperature,
1655 * using actual level populations with 1s of H0 */
1656 else if( strcmp(chELEM,"SPIN") == 0 )
1657 {
1658 *TeMean = mean.TempH_21cmSpinMean[dim][0] / SDIV(mean.TempH_21cmSpinMean[dim][1]);
1659 }
1660 /* return temperature deduced from ratio of 21 cm and Lya optical depths */
1661 else if( strcmp(chELEM,"OPTI") == 0 )
1662 {
1663 *TeMean =
1664 3.84e-7 * iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().TauCon() /
1665 SDIV( HFLines[0].Emis().TauCon() );
1666 }
1667 /* mean temp weighted by H2 density */
1668 else if( strcmp(chELEM,"H2 ") == 0 )
1669 {
1670 if( mean.TempIonMean[dim][ipHYDROGEN][2][1] > SMALLFLOAT )
1671 *TeMean = mean.TempIonMean[dim][ipHYDROGEN][2][0] /
1672 mean.TempIonMean[dim][ipHYDROGEN][2][1];
1673
1674 else
1675 *TeMean = 0.;
1676 }
1677 /* this is temperature weighted by eden */
1678 else if( strcmp(chELEM,"TENE") == 0 )
1679 {
1680 if( mean.TempEdenMean[dim][1] > SMALLFLOAT )
1681 *TeMean = mean.TempEdenMean[dim][0]/mean.TempEdenMean[dim][1];
1682 else
1683 *TeMean = 0.;
1684 }
1685 /* four spaces mean to return simple mean of Te */
1686 else if( strcmp(chELEM," ") == 0 )
1687 {
1688 if( mean.TempMean[dim][1] > SMALLFLOAT )
1689 *TeMean = mean.TempMean[dim][0]/mean.TempMean[dim][1];
1690 else
1691 *TeMean = 0.;
1692 }
1693 else
1694 {
1695 fprintf( ioQQQ, " cdTemp called with ion=0 and unknown quantity, =%4.4s\n",
1696 chLabel );
1697 return 1;
1698 }
1699
1700 /* say things are ok */
1701 return 0;
1702 }
1703
1704 /* find which element this is */
1705 nelem = 0;
1706 while( nelem < LIMELM &&
1707 strcmp(chELEM,elementnames.chElementNameShort[nelem]) != 0 )
1708 {
1709 ++nelem;
1710 }
1711
1712 /* if element not recognized and above loop falls through, nelem is LIMELM */
1713 if( nelem >= LIMELM )
1714 {
1715 fprintf( ioQQQ, " cdTemp called with unknown element chLabel, =%4.4s\n",
1716 chLabel );
1717 return 1;
1718 }
1719
1720 /* sanity check - does this ionization stage exist?
1721 * IonStage on spectroscopic scale, nelem on c */
1722
1723 /* ion will be used as pointer within the aaa array that contains average values,
1724 * done this way to prevent lint from false problem in access to aaa array */
1725 ion = IonStage - 1;
1726
1727 if( ion > nelem+1 || ion < 0 )
1728 {
1729 fprintf( ioQQQ, " cdTemp asked to return ionization stage %ld for element %4.4s but this is not physical.\n",
1730 IonStage, chLabel );
1731 return 1;
1732 }
1733
1734 /* get average, aaa is filled in from 0 */
1735 /* aaa is dim'd LIMELM+1 so largest arg is LIMELM */
1736 /* MeanIon uses C scale for nelem */
1737 mean.MeanIon('t', nelem,dim,&ip,aaa,false);
1738 *TeMean = pow((realnum)10.f,aaa[ion]);
1739 return 0;
1740}
1741
1742/*************************************************************************
1743 *
1744 * cdRead routine to read in command lines
1745 * called by user when cloudy used as subroutine
1746 * called by maincl when used as a routine
1747 *
1748 ************************************************************************/
1749
1750/* returns the number of lines available in command stack
1751 * this is limit to how many more commands can be read */
1753 /* the string containing the command */
1754 const char *chInputLine )
1755{
1756 char *chEOL , /* will be used to search for end of line symbols */
1757 chCARD[INPUT_LINE_LENGTH],
1758 chLocal[INPUT_LINE_LENGTH];
1759 bool lgComment;
1760
1761 DEBUG_ENTRY( "cdRead()" );
1762
1763 /* this is set false when code loaded, set true when cdInit called,
1764 * this is check that cdInit was called first */
1765 if( !lgcdInitCalled )
1766 {
1767 printf(" cdInit was not called first - this must be the first call.\n");
1769 }
1770
1771 /* totally ignore any card starting with a #, *, space, //, or %
1772 * but want to include special "c " type of comment
1773 * >>chng 06 sep 04 use routine to check for comments */
1774 if( ( lgInputComment( chInputLine ) ||
1775 /* these two are end-of-input-stream sentinels */
1776 chInputLine[0] == '\n' || chInputLine[0] == ' ' ) &&
1777 /* option to allow "C" lines through */
1778 ! ( chInputLine[0] == 'C' || chInputLine[0] == 'c' ) )
1779 {
1780 /* return value is number of lines that can still be stuffed in */
1781 return NKRD - input.nSave;
1782 }
1783
1784 /***************************************************************************
1785 * validate a location to store this line image, then store the version *
1786 * that has been truncated from special end of line characters *
1787 * stored image will have < INPUT_LINE_LENGTH valid characters *
1788 ***************************************************************************/
1789
1790 /* this will now point to the next free slot in the line image save array
1791 * this is where we will stuff this line image */
1792 ++input.nSave;
1793
1794 if( input.nSave >= NKRD )
1795 {
1796 /* too many input commands were entered - bail */
1797 fprintf( ioQQQ,
1798 " Too many line images entered to cdRead. The limit is %d\n",
1799 NKRD );
1800 fprintf( ioQQQ,
1801 " The limit to the number of allowed input lines is %d. This limit was exceeded. Sorry.\n",
1802 NKRD );
1803 fprintf( ioQQQ,
1804 " This limit is set by the variable NKRD which appears in input.h \n" );
1805 fprintf( ioQQQ,
1806 " Increase it everywhere it appears.\n" );
1808 }
1809
1810 strncpy( chLocal, chInputLine, INPUT_LINE_LENGTH );
1811 // strncpy will pad chLocal with null bytes if chInputLine is shorter than
1812 // INPUT_LINE_LENGTH characters, so this indicates an overlong input string
1813 if( chLocal[INPUT_LINE_LENGTH-1] != '\0' )
1814 {
1815 chLocal[INPUT_LINE_LENGTH-1] = '\0';
1816 fprintf(ioQQQ," PROBLEM cdRead, while parsing the following input line:\n %s\n",
1817 chInputLine);
1818 fprintf(ioQQQ," found that the line is longer than %i characters, the longest possible line.\n",
1820 fprintf(ioQQQ," Please make the command line no longer than this limit.\n");
1821 }
1822
1823 /* now kill any part of line image after special end of line character,
1824 * this stops info user wants ignored from entering after here */
1825 if( (chEOL = strchr_s(chLocal, '\n' ) ) != NULL )
1826 {
1827 *chEOL = '\0';
1828 }
1829 if( (chEOL = strchr_s(chLocal, '%' ) ) != NULL )
1830 {
1831 *chEOL = '\0';
1832 }
1833 /* >>chng 02 apr 10, add this char */
1834 if( (chEOL = strchr_s(chLocal, '#' ) ) != NULL )
1835 {
1836 *chEOL = '\0';
1837 }
1838 if( (chEOL = strchr_s(chLocal, ';' ) ) != NULL )
1839 {
1840 *chEOL = '\0';
1841 }
1842 if( (chEOL = strstr_s(chLocal, "//" ) ) != NULL )
1843 {
1844 *chEOL = '\0';
1845 }
1846
1847 /* now do it again, since we now want to make sure that there is a trailing space
1848 * if the line is shorter than 80 char, test on null is to keep lint happy */
1849 if( (chEOL = strchr_s( chLocal, '\0' )) == NULL )
1850 TotalInsanity();
1851
1852 /* pad with two spaces if short enough,
1853 * if not short enough for this to be done, then up to user to create correct input */
1854 if( chEOL-chLocal < INPUT_LINE_LENGTH-2 )
1855 {
1856 strcat( chLocal, " " );
1857 }
1858
1859 /* save string in master array for later use in readr */
1860 strcpy( input.chCardSav[input.nSave], chLocal );
1861
1862 /* copy string to chCARD, then convert this to caps to check for keywords*/
1863 strcpy( chCARD, chLocal );
1864 caps(chCARD);/* convert to caps */
1865
1866 // is this a comment? if so do not check for keywords
1867 lgComment = false;
1868 if( strncmp(chCARD , "C ", 2 )==0 )
1869 lgComment = true;
1870 bool lgTitle = ( strncmp(chCARD, "TITL", 4) == 0 );
1871
1872 /* check whether this is a trace command, turn on printout if so */
1873 if( strncmp(chCARD,"TRACE",5) == 0 )
1874 trace.lgTrace = true;
1875
1876 /* print input lines if trace specified */
1877 if( trace.lgTrace )
1878 fprintf( ioQQQ,"cdRead=%s=\n",input.chCardSav[input.nSave] );
1879
1880 // remove string between double quotes
1881 char chFilename[INPUT_LINE_LENGTH],
1882 chTemp[INPUT_LINE_LENGTH];
1883 // this has to have copy of line for GetQuote to function
1884 strcpy( chTemp , input.chCardSav[input.nSave] );
1885 GetQuote( chFilename , chCARD , chTemp , false );
1886
1887 /* now check whether VARY is specified */
1888 if( !lgComment && !lgTitle && nMatch("VARY",chCARD) )
1889 /* a optimize flag was on the line image */
1890 optimize.lgVaryOn = true;
1891
1892 /* now check whether line is "no vary" command */
1893 if( strncmp(chCARD,"NO VARY",7) == 0 )
1894 optimize.lgNoVary = true;
1895
1896 /* now check whether line is "grid" command */
1897 if( strncmp(chCARD,"GRID",4) == 0 )
1898 {
1899 grid.lgGrid = true;
1900 ++grid.nGridCommands;
1901 }
1902
1903 /* NO BUFFERING turn off buffered io for standard output,
1904 * used to get complete output when code crashes */
1905 if( strncmp(chCARD,"NO BUFF",7) == 0 )
1906 {
1907 /* if output has already been redirected (e.g. using cdOutput()) then
1908 * ignore this command, a warning will be printed in ParseDont() */
1909 /* stdout may be a preprocessor macro, so lets be really careful here */
1910 FILE *test = stdout;
1911 if( ioQQQ == test )
1912 {
1913 fprintf( ioQQQ, " cdRead found NO BUFFERING command, redirecting output to stderr now.\n" );
1914 /* make sure output is not lost */
1915 fflush( ioQQQ );
1916 /* stderr is always unbuffered */
1917 //ioQQQ = stderr;
1918 setvbuf( ioQQQ, NULL, _IONBF, 0);
1919 /* will be used to generate comment at end */
1920 input.lgSetNoBuffering = true;
1921 }
1922 else if( !save.chOutputFile.empty() )
1923 {
1924 fprintf( ioQQQ, " cdRead found NO BUFFERING command, reopening file %s now.\n",
1925 save.chOutputFile.c_str() );
1926 fclose( ioQQQ );
1927 // using AS_SILENT_TRY assures that open_data will not write to ioQQQ
1928 ioQQQ = open_data( save.chOutputFile.c_str(), "a", AS_SILENT_TRY );
1929 if( ioQQQ == NULL )
1930 {
1931 // ioQQQ is no longer valid, so switch to stderr and abort...
1932 ioQQQ = stderr;
1933 fprintf( ioQQQ, " cdRead failed to reopen %s, aborting!\n",
1934 save.chOutputFile.c_str() );
1936 }
1937 if( setvbuf( ioQQQ, NULL, _IONBF, 0 ) != 0 )
1938 fprintf( ioQQQ, " PROBLEM -- cdRead failed to set NO BUFFERING mode.\n" );
1939 else
1940 input.lgSetNoBuffering = true;
1941 }
1942 }
1943
1944 /* use grid command as substitute for optimize command */
1945 if( strncmp(chCARD,"OPTI",4) == 0 || strncmp(chCARD,"GRID",4) == 0 )
1946 /* optimize command read in */
1947 optimize.lgOptimr = true;
1948
1949 return NKRD - input.nSave;
1950}
1951
1952/* wrapper to close all save files */
1954{
1955
1956 DEBUG_ENTRY( "cdClosePunchFiles()" );
1957
1958 CloseSaveFiles( true );
1959 return;
1960}
t_called called
Definition called.cpp:5
long int nzone
Definition cddefines.cpp:14
FILE * ioQQQ
Definition cddefines.cpp:7
FILE * ioStdin
Definition cddefines.cpp:8
bool lgAbort
Definition cddefines.cpp:10
#define ASSERT(exp)
Definition cddefines.h:578
#define STATIC
Definition cddefines.h:97
long nMatch(const char *chKey, const char *chCard)
Definition service.cpp:451
const int LIMELM
Definition cddefines.h:258
const int INPUT_LINE_LENGTH
Definition cddefines.h:254
#define EXIT_FAILURE
Definition cddefines.h:140
void cap4(char *chCAP, const char *chLab)
Definition service.cpp:240
#define cdEXIT(FAIL)
Definition cddefines.h:434
const char * strstr_s(const char *haystack, const char *needle)
Definition cddefines.h:1429
int GetQuote(char *chLabel, char *chCard, char *chCardRaw, bool lgABORT)
Definition service.cpp:513
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
const char * strchr_s(const char *s, int c)
Definition cddefines.h:1439
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition cddefines.h:812
void caps(char *chCard)
Definition service.cpp:280
NORETURN void TotalInsanity(void)
Definition service.cpp:886
sys_float SDIV(sys_float x)
Definition cddefines.h:952
const int ipHYDROGEN
Definition cddefines.h:305
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
void cdDepth_depth(double cdDepth[])
Definition cddrive.cpp:976
double cdCooling_last()
Definition cddrive.cpp:333
int cdColm(const char *chLabel, long int ion, double *theocl)
Definition cddrive.cpp:636
void cdCautions(FILE *ioOUT)
Definition cddrive.cpp:226
void cdPrintCommands(FILE *ioOUT)
Definition cddrive.cpp:513
long int cdEmis(char *chLabel, realnum wavelength, double *emiss)
Definition cddrive.cpp:533
void cdLine_ip(long int ipLine, double *relint, double *absint)
Definition cddrive.cpp:1414
long int cdnZone()
Definition cddrive.cpp:1049
double cdEDEN_last()
Definition cddrive.cpp:390
void cdPressure_last(double *PresTotal, double *PresGas, double *PresRad)
Definition cddrive.cpp:1028
void cdSurprises(FILE *ioOUT)
Definition cddrive.cpp:282
void cdNotes(FILE *ioOUT)
Definition cddrive.cpp:309
static bool lgCalled
Definition cddrive.cpp:425
int cdRead(const char *chInputLine)
Definition cddrive.cpp:1752
int cdIonFrac(const char *chLabel, long int IonStage, double *fracin, const char *chWeight, bool lgDensity)
Definition cddrive.cpp:1072
void cdClosePunchFiles()
Definition cddrive.cpp:1953
void cdEmis_ip(long int ipLine, double *emiss, bool lgEmergent)
Definition cddrive.cpp:613
static struct timeval before
Definition cddrive.cpp:441
void cdPrtWL(FILE *io, realnum wl)
Definition cddrive.cpp:157
double cdB21cm()
Definition cddrive.cpp:1561
void cdInput(const char *filename, const char *mode)
Definition cddrive.cpp:1530
void cdTimescales(double *TTherm, double *THRecom, double *TH2)
Definition cddrive.cpp:249
int cdDrive()
Definition cddrive.cpp:77
void cdNwcns(bool *lgAbort_ret, long int *NumberWarnings, long int *NumberCautions, long int *NumberNotes, long int *NumberSurprises, long int *NumberTempFailures, long int *NumberPresFailures, long int *NumberIonFailures, long int *NumberNeFailures)
Definition cddrive.cpp:1481
void cdDate(char chString[])
Definition cddrive.cpp:361
void cdSetExecTime()
Definition cddrive.cpp:474
void cdNoExec()
Definition cddrive.cpp:404
void cdVersion(char chString[])
Definition cddrive.cpp:346
void cdWarnings(FILE *ioPNT)
Definition cddrive.cpp:198
STATIC void cdClock(struct timeval *clock_dat)
Definition cddrive.cpp:444
long int cdLine(const char *chLabel, realnum wavelength, double *relint, double *absint)
Definition cddrive.cpp:1228
void cdErrors(FILE *ioOUT)
Definition cddrive.cpp:910
void cdPressure_depth(double TotalPressure[], double GasPressure[], double RadiationPressure[])
Definition cddrive.cpp:1001
void cdTalk(bool lgTOn)
Definition cddrive.cpp:1548
int cdTemp(const char *chLabel, long int IonStage, double *TeMean, const char *chWeight)
Definition cddrive.cpp:1602
long debugLine(realnum wavelength)
Definition cddrive.cpp:1194
double cdTemp_last()
Definition cddrive.cpp:1061
void cdReasonGeo(FILE *ioOUT)
Definition cddrive.cpp:175
double cdExecTime()
Definition cddrive.cpp:481
void cdOutput(const char *filename, const char *mode)
Definition cddrive.cpp:1517
double cdHeating_last()
Definition cddrive.cpp:376
double cdH2_colden(long iVib, long iRot)
Definition mole_h2.cpp:2318
bool lgcdInitCalled
Definition cdinit.cpp:34
LinSv * LineSv
Definition cdinit.cpp:70
static t_version & Inst()
Definition cddefines.h:175
realnum column
Definition mole.h:359
bool cloudy()
Definition cloudy.cpp:38
t_colden colden
Definition colden.cpp:5
#define ipCOL_HeHp
Definition colden.h:24
#define ipCOL_H2s
Definition colden.h:18
#define ipCOL_H2g
Definition colden.h:16
#define ipCOL_HMIN
Definition colden.h:14
#define ipCOL_H2p
Definition colden.h:20
#define ipCOL_H3p
Definition colden.h:28
t_conv conv
Definition conv.cpp:5
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition cpu.cpp:625
static t_cpu cpu
Definition cpu.h:355
UNUSED const realnum BIGFLOAT
Definition cpu.h:189
@ AS_LOCAL_ONLY
Definition cpu.h:208
@ AS_SILENT_TRY
Definition cpu.h:208
const realnum SMALLFLOAT
Definition cpu.h:191
t_dense dense
Definition dense.cpp:24
t_elementnames elementnames
t_grid grid
Definition grid.cpp:5
bool grid_do(void)
Definition grid_do.cpp:19
void InitCoreload(void)
t_input input
Definition input.cpp:12
bool lgInputComment(const char *chLine)
Definition input.cpp:18
#define NKRD
Definition input.h:10
t_iso_sp iso_sp[NISO][LIMELM]
Definition iso.cpp:8
const int ipH1s
Definition iso.h:27
const int ipH2p
Definition iso.h:29
const int ipH_LIKE
Definition iso.h:62
t_LineSave LineSave
Definition lines.cpp:5
realnum WavlenErrorGet(realnum wavelength)
t_mean mean
Definition mean.cpp:17
molezone * findspecieslocal(const char buf[])
static realnum * wavelength
t_noexec noexec
Definition noexec.cpp:5
t_optimize optimize
Definition optimize.cpp:5
void CloseSaveFiles(bool lgFinal)
t_phycon phycon
Definition phycon.cpp:6
t_pressure pressure
Definition pressure.cpp:5
void prt_wl(FILE *ioOUT, realnum wl)
Definition prt.cpp:13
static long int * ipLine
t_radius radius
Definition radius.cpp:5
t_save save
Definition save.cpp:5
t_struc struc
Definition struc.cpp:6
TransitionList HFLines("HFLines", &AnonStates)
t_thermal thermal
Definition thermal.cpp:5
t_timesc timesc
Definition timesc.cpp:5
t_trace trace
Definition trace.cpp:5
t_warnings warnings
Definition warnings.cpp:11