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

Generated for cloudy by doxygen 1.8.4