cloudy  trunk
mole_h2_io.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 /*H2_ParseSave parse the save h2 command */
4 /*H2_PunchDo save some properties of the large H2 molecule */
5 /*chMolBranch returns a char with the spectroscopic branch of a transition */
6 /*H2_Prt_line_tau print line optical depths, called from premet in response to print line optical depths command*/
7 /*H2_PunchLineStuff include H2 lines in punched optical depths, etc, called from SaveLineStuff */
8 /*H2_Punch_line_data save line data for H2 molecule */
9 /*H2_Read_hminus_distribution read distribution function for H2 population following formation from H minus */
10 /*H2_ReadDissprob read dissociation probabilities and kinetic energies for all electronic levels */
11 /*H2_ReadEnergies read energies for all electronic levels */
12 /*H2_ReadTransprob read transition probabilities */
13 /*H2_Prt_Zone print H2 info into zone results, called from prtzone for each printed zone */
14 /*H2_ParseSave parse the save h2 command */
15 /*H2_Prt_column_density print H2 info into zone results, called from prtzone for each printed zone */
16 /*H2_LinesAdd add in explicit lines from the large H2 molecule, called by lines_molecules */
17  /*cdH2_Line returns 1 if we found the line,
18  * or false==0 if we did not find the line because ohoto-para transition
19  * or upper level has lower energy than lower level */
20 #include "cddefines.h"
21 #include "physconst.h"
22 #include "save.h"
23 #include "hmi.h"
24 #include "prt.h"
25 #include "secondaries.h"
26 #include "grainvar.h"
27 #include "input.h"
28 #include "phycon.h"
29 #include "rfield.h"
30 #include "hyperfine.h"
31 #include "thermal.h"
32 #include "lines.h"
33 #include "lines_service.h"
34 #include "mole.h"
35 #include "dense.h"
36 #include "radius.h"
37 #include "colden.h"
38 #include "taulines.h"
39 #include "h2.h"
40 #include "h2_priv.h"
41 #include "cddrive.h"
42 #include "doppvel.h"
43 #include "doppvel.h"
44 #include "parser.h"
45 
46 static realnum thresh_punline_h2;
47 
48 /*H2_LinesAdd add in explicit lines from the large H2 molecule, called by lines_molecules */
49 void diatomics::H2_LinesAdd(void)
50 {
51  /* H2 not on, so space not allocated */
52  if( !lgEnabled )
53  return;
54 
55  DEBUG_ENTRY( "H2_LinesAdd()" );
56 
57  if( strcmp( "H2 ", label.c_str() ) == 0 )
58  {
59  /* >>chng 05 nov 04, make info copies of these lines up here
60  * these are among the strongest lines in the 2 micron window and some have nearly the same
61  * wavelength as far weaker lines that may come before them in the line stack. in that case
62  * cdLine would find the much weaker line with the same wavelength.
63  * put strong H2 lines first so that line search will find these, and not far weaker
64  * lines with nearly the same wavelength - these will be duplicated in the output but
65  * these are here for into (the 'i') so does no harm
66  */
67 
68  /* >>chng 05 dec 22, had hand entered wavelength in A as second parameter. This gave
69  * rounded off result when set line precision 5 was used. now uses same logic that
70  * PutLine will eventually use - simply enter same wl in Ang */
71  /* 1-0 S(4) - 18910 */
72  lindst( trans[ ipTransitionSort[ ipEnergySort[0][1][6] ][ ipEnergySort[0][0][4] ] ], "H2 ", 'i', false, "H2 line");
73  /* 1-0 S(3) - 19570 */
74  lindst( trans[ ipTransitionSort[ ipEnergySort[0][1][5] ][ ipEnergySort[0][0][3] ] ], "H2 ", 'i', false, "H2 line");
75  /* 1-0 S(2) - 20330 */
76  lindst( trans[ ipTransitionSort[ ipEnergySort[0][1][4] ][ ipEnergySort[0][0][2] ] ], "H2 ", 'i', false, "H2 line");
77  /* 1-0 S1) - 21210 */
78  lindst( trans[ ipTransitionSort[ ipEnergySort[0][1][3] ][ ipEnergySort[0][0][1] ] ], "H2 ", 'i', false, "H2 line");
79  /* 1-0 S(0) - 22230 */
80  lindst( trans[ ipTransitionSort[ ipEnergySort[0][1][2] ][ ipEnergySort[0][0][0] ] ], "H2 ", 'i', false, "H2 line");
81  /* start Q branch - selection rule requires that J be non-zero, so no Q(0) */
82  /* 1-0 Q(2) - 24130 */
83  lindst( trans[ ipTransitionSort[ ipEnergySort[0][1][2] ][ ipEnergySort[0][0][2] ] ], "H2 ", 'i', false, "H2 line");
84  /* 1-0 Q(1) - 24060 */
85  lindst( trans[ ipTransitionSort[ ipEnergySort[0][1][1] ][ ipEnergySort[0][0][1] ] ], "H2 ", 'i', false, "H2 line");
86  }
87 
88 
89  for( TransitionList::iterator tr = trans.begin(); tr != rad_end; ++tr )
90  {
91  qList::iterator Hi = ( (*tr).Hi() );
92  if( (*Hi).n() >= nElecLevelOutput ) continue;
93  qList::iterator Lo = ( (*tr).Lo() );
94  /* all ground vib state rotation lines - first is J to J-2 */
95  PutLine( *tr, "diatoms lines", label.c_str() );
96  if( LineSave.ipass == 0 )
97  {
98  H2_SaveLine[(*Hi).n()][(*Hi).v()][(*Hi).J()][(*Lo).n()][(*Lo).v()][(*Lo).J()] = 0.;
99  }
100  else if( LineSave.ipass == 1 )
101  {
102  H2_SaveLine[(*Hi).n()][(*Hi).v()][(*Hi).J()][(*Lo).n()][(*Lo).v()][(*Lo).J()] +=
103  (realnum)( radius.dVeffAper * (*tr).Emis().xIntensity() );
104  }
105  }
106 
107  return;
108 }
109 
110 /*H2_ParseSave parse the save h2 command */
111 void diatomics::H2_ParseSave( Parser &p ,
112  char *chHeader)
113 {
114  DEBUG_ENTRY( "H2_ParseSave()" );
115 
116  save.whichDiatomToPrint[save.nsave] = &(*this);
117 
118  /* this provides info on the large H2 molecule */
119  if( p.nMatch("COLU") )
120  {
121  /* save column density */
122  strcpy( save.chSave[save.nsave], "H2cl" );
123 
124  /* this is an option to scan off highest vib and rot states
125  * to save pops - first is limit to vibration, then rotation
126  * if no number is entered then 0 is set and all levels punched */
127  /* now get vib limit */
128  save.punarg[save.nsave][0] = (realnum)p.getNumberDefault(
129  "H2 vibration state",0.0);
130 
131  /* highest rotation */
132  save.punarg[save.nsave][1] = (realnum)p.getNumberDefault(
133  "H2 rotation state",0.0);
134  /* this says whether to save triplets or a matrix for output -
135  * default is triplets, so only check for matrix */
136  if( p.nMatch( "MATR" ) )
137  {
138  /* matrix */
139  save.punarg[save.nsave][2] = 1;
140  sprintf( chHeader, "#vib\trot\tcolumn density\n" );
141  }
142  else
143  {
144  /* triplets */
145  save.punarg[save.nsave][2] = -1;
146  sprintf( chHeader, "#vib\trot\tEner(K)\tcolden\tcolden/stat wght\tLTE colden\tLTE colden/stat wght\n" );
147  }
148  }
149  else if( p.nMatch("COOL") )
150  {
151  /* heating and cooling rates */
152  strcpy( save.chSave[save.nsave], "H2co" );
153  sprintf( chHeader,
154  "#H2 depth\ttot cool\tTH Sol\tBig Sol\tTH pht dis\tpht dis\tTH Xcool\tXcool \n" );
155  }
156 
157  else if( p.nMatch("CREA") )
158  {
159  /* H2 creation rates */
160  fprintf( ioQQQ, " This command has been superseded by the \"creation\" option of the \"save chemistry rates\" command.\n" );
161  fprintf( ioQQQ, " Sorry.\n" );
162  cdEXIT(EXIT_FAILURE);
163  }
164  else if( p.nMatch("DEST") )
165  {
166  /* save H2 destruction - output destruction rates */
167  fprintf( ioQQQ, " This command has been superseded by the \"destruction\" option of the \"save chemistry rates\" command.\n" );
168  fprintf( ioQQQ, " Sorry.\n" );
169  cdEXIT(EXIT_FAILURE);
170  }
171 
172  else if( p.nMatch("HEAT") )
173  {
174  /* heating and cooling rates */
175  strcpy( save.chSave[save.nsave], "H2he" );
176  sprintf( chHeader,
177  "#H2 depth\ttot Heat\tHeat(big)\tHeat(TH85)\tDissoc(Big)\tDissoc(TH85) \n" );
178  }
179 
180  else if( p.nMatch("LEVE") )
181  {
182  /* save H2 level energies */
183  strcpy( save.chSave[save.nsave], "H2le" );
184  sprintf( chHeader,
185  "#H2 v\tJ\tenergy(wn)\tstat wght\tSum As" );
186  char chHoldit[chN_X_COLLIDER+12];
187  for( int nColl=0; nColl<N_X_COLLIDER; ++nColl )
188  {
189  /* labels for all colliders */
190  sprintf(chHoldit,"\tCritDen %s",chH2ColliderLabels[nColl]);
191  strcat( chHeader , chHoldit );
192  }
193  strcat( chHeader , "\n" );
194  }
195 
196  else if( p.nMatch("LINE") )
197  {
198  /* save H2 lines - all in X */
199  strcpy( save.chSave[save.nsave], "H2ln" );
200  sprintf( chHeader,
201  "#H2 line\tEhi\tVhi\tJhi\tElo\tVlo\tJlo\twl(mic)\twl(lab)\tlog L or I\tI/Inorm\tExcit(hi, K)\tg_u h nu * Aul\n" );
202  /* first optional number changes the threshold of weakest line to print*/
203  /* fe2thresh is intensity relative to normalization line,
204  * normally Hbeta, and is set to zero in zero.c */
205 
206  /* threshold for faintest line to save, default is 1e-4 of norm line */
207  thresh_punline_h2 = (realnum)p.getNumberDefaultNegImplLog(
208  "faintest line to save",1e-4);
209 
210  /* lines from how many electronic states? default is one, just X, and is
211  * obtained with GROUND keyword. ALL will produce all lines from all levels.
212  * else, if a number is present, will be the number. if no number, no keyword,
213  * appear then just ground */
214  if( p.nMatch( "ELEC" ) )
215  {
216  if( p.nMatch(" ALL") )
217  {
218  /* all electronic levels - when done, will set upper limit, the
219  * number of electronic levels actually computed, don't know this yet,
220  * so signify with negative number */
221  nElecLevelOutput = -1;
222  }
223  else if( p.nMatch("GROU") )
224  {
225  /* just the ground electronic state */
226  nElecLevelOutput = 1;
227  }
228  else
229  {
230  nElecLevelOutput = (int)p.getNumberDefault(
231  "electronic levels for output",1.0);
232  }
233  }
234  }
235 
236  else if( p.nMatch(" PDR") )
237  {
238  /* creation and destruction processes */
239  strcpy( save.chSave[save.nsave], "H2pd" );
240  sprintf( chHeader, "#H2 creation, destruction. \n" );
241  }
242  else if( p.nMatch("POPU") )
243  {
244  /* save populations */
245  strcpy( save.chSave[save.nsave], "H2po" );
246 
247  /* this is an option to scan off highest vib and rot states
248  * to save pops - first is limit to vibration, then rotation
249  * if no number is entered then 0 is set and all levels punched */
250  /* now get vib lim */
251  save.punarg[save.nsave][0] = (realnum)p.getNumberDefault(
252  "highest H2 save vibration state",0.0);
253 
254  /* this is limit to rotation quantum index */
255  save.punarg[save.nsave][1] = (realnum)p.getNumberDefault(
256  "highest H2 save rotation state",0.0);
257 
258  if( p.nMatch( "ZONE" ) )
259  {
260  /* save v=0 pops for each zone, all along one line */
261  save.punarg[save.nsave][2] = 0;
262  sprintf( chHeader, "#depth\torth\tpar\te=1 rel pop\te=2 rel pop\tv,J rel pops\n" );
263  }
264  else
265  {
266  /* will not do zone output, only output at the end of the calculation
267  * now check whether to save triplets or a matrix for output -
268  * default is triplets, so only check for matrix */
269  if( p.nMatch( "MATR" ) )
270  {
271  /* matrix */
272  save.punarg[save.nsave][2] = 1;
273  sprintf( chHeader, "#vib\trot\tpops\n" );
274  }
275  else
276  {
277  /* triplets */
278  save.punarg[save.nsave][2] = -1;
279  sprintf( chHeader, "#vib\trot\ts\tenergy(wn)\tpops/H2\told/H2\tpops/g/H2\tdep coef\tFin(Col)\tFout(col)\tRCout\tRRout\tRCin\tRRin\n" );
280  }
281  }
282  }
283 
284  else if( p.nMatch("RATE") )
285  {
286  /* save h2 rates - creation and destruction rates */
287  strcpy( save.chSave[save.nsave], "H2ra" );
288  sprintf( chHeader,
289  "#depth\tN(H2)\tN(H2)/u(H2)\tA_V(star)\tn(Eval)"
290  "\tH2/Htot\trenorm\tfrm grn\tfrmH-\tdstTH85\tBD96\tELWERT\tBigH2\telec->H2g\telec->H2s"
291  "\tG(TH85)\tG(DB96)\tCR\tEleclife\tShield(BD96)\tShield(H2)\tBigh2/G0(spc)\ttot dest"
292  "\tHeatH2Dish_TH85\tHeatH2Dexc_TH85\tHeatDish_BigH2\tHeatDexc_BigH2\thtot\n" );
293  }
294  else if( p.nMatch("SOLO") )
295  {
296  /* rate of Solomon process then fracs of exits from each v, J level */
297  strcpy( save.chSave[save.nsave], "H2so" );
298  sprintf( chHeader,
299  "#depth\tSol tot\tpump/dissoc\tpump/dissoc BigH2\tavH2g\tavH2s\tH2g chem/big H2\tH2s chem/big H2\tfrac H2g BigH2\tfrac H2s BigH2\teHi\tvHi\tJHi\tvLo\tJLo\tfrac\twl(A)\n" );
300  }
301  else if( p.nMatch("SPEC") )
302  {
303  /* special save command*/
304  strcpy( save.chSave[save.nsave], "H2sp" );
305  sprintf( chHeader,
306  "#depth\tspecial\n" );
307  }
308  else if( p.nMatch("TEMP") )
309  {
310  /* various temperatures for neutral/molecular gas */
311  strcpy( save.chSave[save.nsave], "H2te" );
312  sprintf( chHeader,
313  "#depth\tH2/H\tn(1/0)\tn(ortho/para)\tT(1/0)\tT(2/0)\tT(3/0)\tT(3/1)\tT(4/0)\tT(kin)\tT(21cm)\tT_sum(1/0)\tT_sum(2/0)\tT_sum(3/0)\tT_sum(3/1)\tT_sum(4/0) \n");
314  }
315  else if( p.nMatch("THER") )
316  {
317  /* thermal heating cooling processes involving H2 */
318  strcpy( save.chSave[save.nsave], "H2th" );
319  sprintf( chHeader,
320  "#depth\tH2/H\tn(1/0)\tn(ortho/para)\tT(1/0)\tT(2/0)\tT(3/0)\tT(3/1)\tT(4/0)\tT(kin)\tT(21cm)\tT_sum(1/0)\tT_sum(2/0)\tT_sum(3/0)\tT_sum(3/1)\tT_sum(4/0) \n");
321  }
322  else
323  {
324  fprintf( ioQQQ,
325  " There must be a second key; they are RATE, LINE, COOL, COLUMN, _PDR, SOLOmon, TEMP, and POPUlations\n" );
326  cdEXIT(EXIT_FAILURE);
327  }
328  return;
329 }
330 
331 
332 /*H2_Prt_Zone print H2 info into zone results, called from prtzone for each printed zone */
333 void diatomics::H2_Prt_Zone(void)
334 {
335  /* no print if H2 not turned on, or not computed for these conditions */
336  if( !lgEnabled || !nCall_this_zone )
337  return;
338 
339  DEBUG_ENTRY( "H2_Prt_Zone()" );
340 
341  fprintf( ioQQQ, " %s density ", label.c_str() );
342  fprintf(ioQQQ,PrintEfmt("%9.2e", (*dense_total)));
343 
344  fprintf( ioQQQ, " orth/par");
345  fprintf(ioQQQ,PrintEfmt("%9.2e", ortho_density / SDIV( para_density )));
346 
347  fprintf( ioQQQ, " v0 J=0,3");
348  fprintf(ioQQQ,PrintEfmt("%9.2e", states[ ipEnergySort[0][0][0] ].Pop() / (*dense_total)));
349  fprintf(ioQQQ,PrintEfmt("%9.2e", states[ ipEnergySort[0][0][1] ].Pop() / (*dense_total)));
350  fprintf(ioQQQ,PrintEfmt("%9.2e", states[ ipEnergySort[0][0][2] ].Pop() / (*dense_total)));
351  fprintf(ioQQQ,PrintEfmt("%9.2e", states[ ipEnergySort[0][0][3] ].Pop() / (*dense_total)));
352 
353  fprintf( ioQQQ, " TOTv=0,3");
354  fprintf(ioQQQ,PrintEfmt("%9.2e", pops_per_vib[0][0] / (*dense_total)));
355  fprintf(ioQQQ,PrintEfmt("%9.2e", pops_per_vib[0][1] / (*dense_total)));
356  fprintf(ioQQQ,PrintEfmt("%9.2e", pops_per_vib[0][2] / (*dense_total)));
357  fprintf(ioQQQ,PrintEfmt("%9.2e", pops_per_vib[0][3] / (*dense_total)));
358  fprintf( ioQQQ, "\n");
359  return;
360 }
361 
362 void diatomics::H2_PrtDepartCoef(void)
363 {
364  /* no print if H2 not turned on, or not computed for these conditions */
365  if( !lgEnabled || !nCall_this_zone )
366  return;
367 
368  DEBUG_ENTRY( "H2_PrtDepartCoef()" );
369 
370  // print departure coefficients
371  fprintf( ioQQQ, " %s departure coefficients\n", label.c_str() );
372  for( long iElec=0; iElec<n_elec_states; ++iElec )
373  {
374  fprintf( ioQQQ, "%li electronic\n", iElec );
375  for( long iVib=0; iVib<=nVib_hi[iElec]; ++iVib )
376  {
377  for( long iRot=0; iRot<Jlowest[iElec]; ++iRot )
378  fprintf( ioQQQ, " -----" );
379  for( long iRot=Jlowest[iElec]; iRot<=nRot_hi[iElec][iVib]; ++iRot )
380  {
381  long i = ipEnergySort[iElec][iVib][iRot];
382  fprintf( ioQQQ, " %5.3f", depart[i] );
383  }
384  fprintf( ioQQQ, "\n" );
385  }
386  fprintf( ioQQQ, "\n" );
387  if( iElec==0 )
388  break;
389  }
390 
391  return;
392 }
393 
394 /*H2_Prt_column_density print H2 info into zone results, called from prtzone for each printed zone */
395 void diatomics::H2_Prt_column_density(
396  /* this is stream used for io, is stdout when called by final,
397  * is save unit when save output generated */
398  FILE *ioMEAN )
399 
400 {
401  int iVibHi;
402 
403  /* no print if H2 not turned on, or not computed for these conditions */
404  if( !lgEnabled || !nCall_this_zone )
405  return;
406 
407  DEBUG_ENTRY( "H2_Prt_column_density()" );
408 
409  fprintf( ioMEAN, " H2 total ");
410  fprintf(ioMEAN,"%7.3f", log10(SDIV(ortho_colden + para_colden)));
411 
412  fprintf( ioMEAN, " H2 ortho ");
413  fprintf(ioMEAN,"%7.3f", log10(SDIV(ortho_colden)));
414 
415  fprintf( ioMEAN, " para");
416  fprintf(ioMEAN,"%7.3f", log10(SDIV(para_colden)));
417 
418  iVibHi = 0;
419  fprintf( ioMEAN, " v0 J=0,3");
420  fprintf(ioMEAN,"%7.3f", log10(SDIV(H2_X_colden[iVibHi][0])));
421  fprintf(ioMEAN,"%7.3f", log10(SDIV(H2_X_colden[iVibHi][1])));
422  fprintf(ioMEAN,"%7.3f", log10(SDIV(H2_X_colden[iVibHi][2])));
423  fprintf(ioMEAN,"%7.3f", log10(SDIV(H2_X_colden[iVibHi][3])));
424 
425  return;
426 }
427 
428 
429 /*H2_ReadTransprob read transition probabilities */
430 void diatomics::H2_ReadTransprob( long int nelec , TransitionList &trns)
431 {
432  const char* cdDATAFILE[N_ELEC] =
433  {
434  "transprob_X.dat",
435  "transprob_B.dat",
436  "transprob_C_plus.dat",
437  "transprob_C_minus.dat",
438  "transprob_B_primed.dat",
439  "transprob_D_plus.dat",
440  "transprob_D_minus.dat"
441  };
442  FILE *ioDATA;
443  char chLine[FILENAME_PATH_LENGTH_2];
444  long int i, n1, n2, n3;
445  long int iVibHi , iVibLo , iRotHi , iRotLo , iElecHi , iElecLo;
446  bool lgEOL;
447 
448  DEBUG_ENTRY( "H2_ReadTransprob()" );
449 
450  /* now open the data file */
451  char chPath[FILENAME_PATH_LENGTH_2];
452  strcpy( chPath, path.c_str() );
453  strcat( chPath, input.chDelimiter );
454  strcat( chPath, cdDATAFILE[nelec] );
455  ioDATA = open_data( chPath , "r" );
456 
457  /* read the first line and check that magic number is ok */
458  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
459  {
460  fprintf( ioQQQ, " H2_ReadTransprob could not read first line of %s\n", cdDATAFILE[nelec]);
461  cdEXIT(EXIT_FAILURE);
462  }
463  i = 1;
464  /* magic number */
465  n1 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
466  n2 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
467  n3 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
468 
469  /* magic number
470  * the following is the set of numbers that appear at the start of level1.dat 01 08 10 */
471  if( ( n1 != 2 ) || ( n2 != 4 ) || ( n3 != 29 ) )
472  {
473  fprintf( ioQQQ,
474  " H2_ReadTransprob: the version of %s is not the current version.\n", cdDATAFILE[nelec] );
475  fprintf( ioQQQ,
476  " I expected to find the number 2 4 29 and got %li %li %li instead.\n" ,
477  n1 , n2 , n3 );
478  fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
479  cdEXIT(EXIT_FAILURE);
480  }
481 
482  long nlines = 0;
483  while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
484  {
485  /* skip comment */
486  if( chLine[0]=='#' )
487  continue;
488  if( chLine[0]=='\n' || chLine[0]=='\0' || chLine[0]==' ' )
489  break;
490 
491  double Aul;
492  int n = sscanf(chLine,"%li\t%li\t%li\t%li\t%li\t%li\t%le",
493  &iElecHi , &iVibHi ,&iRotHi , &iElecLo , &iVibLo , &iRotLo , &Aul );
494  ASSERT( n == 7 );
495  ASSERT( iElecHi == nelec );
496  ASSERT( iElecHi < N_ELEC );
497  ASSERT( iElecLo < N_ELEC );
498 
499  /* check that we actually included the levels in the model representation */
500  if( iVibHi <= nVib_hi[iElecHi] &&
501  iVibLo <= nVib_hi[iElecLo] &&
502  iRotHi <= nRot_hi[iElecHi][iVibHi] &&
503  iRotLo <= nRot_hi[iElecLo][iVibLo])
504  {
505  long ipHi = ipEnergySort[iElecHi][iVibHi][iRotHi];
506  long ipLo = ipEnergySort[iElecLo][iVibLo][iRotLo];
507  double ener = states[ipHi].energy().WN() - states[ipLo].energy().WN();
508  long lineIndex = ipTransitionSort[ipHi][ipLo];
509 
510  /* only lines that have real Aul are added to stack. */
511  trns[lineIndex].AddLine2Stack();
512  trns[lineIndex].Emis().Aul() = (realnum)Aul;
513 
514  /* say that this line exists */
515  lgH2_radiative[ipHi][ipLo] = true;
516  ++nlines;
517 
518  /* prints transitions with negative energies - should not happen */
519  if( ener <= 0. )
520  {
521  fprintf(ioQQQ,"negative energy H2 transition\t%li\t%li\t%li\t%li\t%.2e\t%.2e\n",
522  iVibHi,iVibLo,iRotHi,iRotLo,Aul,ener);
523  ShowMe();
524  cdEXIT(EXIT_FAILURE);
525  }
526  }
527  }
528 
529  fclose( ioDATA );
530  return;
531 }
532 
533 #if 0
534 /*H2_Read_Cosmicray_distribution read distribution function for H2 population following cosmic ray collisional excitation */
535 void H2_Read_Cosmicray_distribution(void)
536 {
537  //CR_PRINT (false), CR_X (1), CR_VIB(15), CR_J(10), CR_EXIT(3)
538 
539  /*>>refer H2 cr excit Tine, S., Lepp, S., Gredel, R., & Dalgarno, A. 1997, ApJ, 481, 282 */
540  FILE *ioDATA;
541  char chLine[FILENAME_PATH_LENGTH_2];
542  long int i, n1, n2, n3, iVib , iRot;
543  long neut_frac;
544  bool lgEOL;
545 
546  DEBUG_ENTRY( "H2_Read_Cosmicray_distribution()" );
547 
548  /* now open the data file */
549  char chPath[FILENAME_PATH_LENGTH_2];
550  strcpy( chPath, path.c_str() );
551  strcat( chPath, input.chDelimiter );
552  strcat( chPath, "H2_CosmicRay_collision.dat" );
553  ioDATA = open_data( chPath, "r" );
554 
555  /* read the first line and check that magic number is ok */
556  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
557  {
558  fprintf( ioQQQ, " H2_Read_Cosmicray_distribution could not read first line of %s\n", "H2_Cosmic_collision.dat");
559  cdEXIT(EXIT_FAILURE);
560  }
561 
562  i = 1;
563  /* magic number */
564  n1 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
565  n2 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
566  n3 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
567 
568  /* magic number
569  * the following is the set of numbers that appear at the start of H2_Cosmic_collision.dat 01 21 03 */
570  if( ( n1 != 1 ) || ( n2 != 21 ) || ( n3 != 3 ) )
571  {
572  fprintf( ioQQQ,
573  " H2_Read_Cosmicray_distribution: the version of %s is not the current version.\n", "H2_Cosmic_collision.dat" );
574  fprintf( ioQQQ,
575  " I expected to find the number 1 21 3 and got %li %li %li instead.\n" ,
576  n1 , n2 , n3 );
577  fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
578  cdEXIT(EXIT_FAILURE);
579  }
580 
581  /* read until not a comment */
582  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
583  BadRead();
584 
585  while( chLine[0]=='#' )
586  {
587  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
588  BadRead();
589  }
590 
591  iRot = 1;
592  iVib = 1;
593  neut_frac = 0;
594  while( iVib >= 0 )
595  {
596  long int j_minus_ji;
597  double a[10];
598 
599  sscanf(chLine,"%li\t%li\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf",
600  &iVib ,&j_minus_ji , &a[0],&a[1],&a[2],&a[3],&a[4],&a[5],&a[6],&a[7],&a[8],&a[9]
601  );
602  /* negative iVib says end of data */
603  if( iVib < 0 )
604  continue;
605 
606  /* cr_rate[CR_X][CR_VIB][CR_J][CR_EXIT];*/
607  /* check that we actually included the levels in the model representation */
608  ASSERT( iVib < CR_VIB );
609  ASSERT( j_minus_ji == -2 || j_minus_ji == +2 || j_minus_ji == 0 );
610  ASSERT( neut_frac < CR_X );
611 
612  /* now make i_minus_ji an array index */
613  j_minus_ji = 1 + j_minus_ji/2;
614  ASSERT( j_minus_ji>=0 && j_minus_ji<=2 );
615 
616  /* option to add Gaussian random mole */
617  for( iRot=0; iRot<CR_J; ++iRot )
618  {
619  cr_rate[neut_frac][iVib][iRot][j_minus_ji] = (realnum)a[iRot];
620  }
621  if( lgH2_NOISECOSMIC )
622  {
623  realnum r;
624  r = (realnum)RandGauss( xMeanNoise , xSTDNoise );
625 
626  for( iRot=0; iRot<CR_J; ++iRot )
627  {
628  cr_rate[neut_frac][iVib][iRot][j_minus_ji] *= (realnum)pow(10.,(double)r);
629  }
630  }
631 
632  if( CR_PRINT )
633  {
634  fprintf(ioQQQ,"cr rate\t%li\t%li", iVib , j_minus_ji );
635  for( iRot=0; iRot<CR_J; ++iRot )
636  {
637  fprintf(ioQQQ,"\t%.3e", cr_rate[neut_frac][iVib][iRot][j_minus_ji] );
638  }
639  fprintf(ioQQQ,"\n" );
640  }
641 
642  /* now get next line */
643  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
644  BadRead();
645  while( chLine[0]=='#' )
646  {
647  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
648  BadRead();
649  }
650  }
651  fclose( ioDATA );
652 
653  return;
654 }
655 #endif
656 
657 class level_tmp
658 {
659 public:
660  bool operator<( const level_tmp& second ) const
661  {
662  if( eWN < second.eWN )
663  return true;
664  else
665  return false;
666  }
667  long n, v, J;
668  double eWN;
669 };
670 
671 /*H2_ReadEnergies read energies for all electronic levels */
672 void diatomics::H2_ReadEnergies( )
673 {
674  DEBUG_ENTRY( "H2_ReadEnergies()" );
675 
676  vector<int> n, v, J;
677  vector<double>eWN;
678 
679  for( long nelec=0; nelec<n_elec_states; ++nelec )
680  {
681  /* get energies out of files */
682  H2_ReadEnergies(nelec,n,v,J,eWN);
683  }
684 
685  vector<level_tmp> levels;
686  levels.resize( n.size() );
687  ASSERT( levels.size() > 0 );
688  for( unsigned i = 0; i < n.size(); ++i )
689  {
690  levels[i].n = n[i];
691  levels[i].v = v[i];
692  levels[i].J = J[i];
693  levels[i].eWN = eWN[i];
694  }
695 
696  // now get levels in energy order (comparison done by operator< of level_tmp class)
697  sort( levels.begin(), levels.end() );
698 
699  // populate states
700  for( vector<level_tmp>::iterator lev = levels.begin(); lev != levels.end(); ++lev )
701  {
702  states.resize( states.size() + 1 );
703  long i = states.size() - 1;
704  states[i].n() = lev->n;
705  states[i].v() = lev->v;
706  states[i].J() = lev->J;
707  states[i].energy().set( lev->eWN, "cm^-1" );
708  /* NB this must be kept parallel with nelem and ionstag in transition struc,
709  * since that struc expects to find the abundances here - abund set in hmole.c */
710  states[i].nelem() = -1;
711  /* this does not mean anything for a molecule */
712  states[i].IonStg() = -1;
713  strcpy( states[i].chLabel(), label.c_str() );
714  }
715 
716  ASSERT( states.size() > 0 );
717  ASSERT( states.size() == levels.size() );
718 
719  for( long nelec=0; nelec<n_elec_states; ++nelec )
720  {
721  ASSERT( nLevels_per_elec[nelec] > 0 );
722  ASSERT( nVib_hi[nelec] > 0 );
723  ASSERT( nVib_hi[nelec] > Jlowest[nelec] );
724 
725  nRot_hi[nelec].resize( nVib_hi[nelec]+1 );
726  nRot_hi[nelec] = 0;
727  }
728 
729  for( qList::const_iterator st = states.begin(); st != states.end(); ++st )
730  {
731  nRot_hi[ (*st).n() ][ (*st).v() ] = MAX2( nRot_hi[ (*st).n() ][ (*st).v() ], (*st).J() );
732  }
733 
734  return;
735 }
736 
737 void diatomics::H2_ReadEnergies( long int nelec, vector<int>& n, vector<int>& v, vector<int>&J, vector<double>& eWN )
738 {
739  DEBUG_ENTRY("diatomics::H2_ReadEnergies()");
740  const char* cdDATAFILE[N_ELEC] =
741  {
742  "energy_X.dat",
743  "energy_B.dat",
744  "energy_C_plus.dat",
745  "energy_C_minus.dat",
746  "energy_B_primed.dat",
747  "energy_D_plus.dat",
748  "energy_D_minus.dat"
749  };
750  /* now open the data file */
751  char chPath[FILENAME_PATH_LENGTH_2];
752  strcpy( chPath, path.c_str() );
753  strcat( chPath, input.chDelimiter );
754  strcat( chPath, cdDATAFILE[nelec] );
755  FILE *ioDATA = open_data( chPath, "r" );
756 
757  char chLine[FILENAME_PATH_LENGTH_2];
758  long int i, n1, n2, n3;
759  bool lgEOL;
760 
761  /* read the first line and check that magic number is ok */
762  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
763  {
764  fprintf( ioQQQ, " H2_ReadEnergies could not read first line of %s\n", cdDATAFILE[nelec]);
765  cdEXIT(EXIT_FAILURE);
766  }
767  i = 1;
768  /* magic number */
769  n1 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
770  n2 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
771  n3 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
772 
773  /* magic number
774  * the following is the set of numbers that appear at the start of level1.dat 01 08 10 */
775  if( ( n1 != 2 ) || ( n2 != 4 ) || ( n3 != 29 ) )
776  {
777  fprintf( ioQQQ,
778  " H2_ReadEnergies: the version of %s is not the current version.\n", cdDATAFILE[nelec] );
779  fprintf( ioQQQ,
780  " I expected to find the number 2 4 29 and got %li %li %li instead.\n" ,
781  n1 , n2 , n3 );
782  fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
783  cdEXIT(EXIT_FAILURE);
784  }
785 
786  /* this will count the number of levels within each electronic state */
787  nLevels_per_elec[nelec] = 0;
788  nVib_hi[nelec] = 0;
789  Jlowest[nelec] = LONG_MAX;
790 
791  while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
792  {
793  /* skip comment */
794  if( chLine[0]=='#' )
795  continue;
796  if( chLine[0]=='\n' || chLine[0]=='\0' || chLine[0]==' ' )
797  break;
798  long iVib, iRot;
799  double energyWN;
800  int nReads = sscanf(chLine,"%li\t%li\t%le", &iVib, &iRot, &energyWN );
801  ASSERT( nReads == 3 );
802  ASSERT( iVib >= 0 );
803  ASSERT( iRot >= 0 );
804  ASSERT( energyWN > 0. || (nelec==0 && iVib==0 && iRot==0 ) );
805 
806  n.push_back( nelec );
807  v.push_back( iVib );
808  J.push_back( iRot );
809  eWN.push_back( energyWN );
810 
811  // update limits
812  nVib_hi[nelec] = MAX2( nVib_hi[nelec], iVib );
813  Jlowest[nelec] = MIN2( Jlowest[nelec], iRot );
814  /* increment number of levels within this electronic state */
815  ++nLevels_per_elec[nelec];
816  }
817 
818  ASSERT( n.size() > 0 );
819  ASSERT( nLevels_per_elec[nelec] > 0 );
820  ASSERT( nVib_hi[nelec] > 0 );
821  ASSERT( nVib_hi[nelec] > Jlowest[nelec] );
822 
823  fclose( ioDATA );
824 
825  return;
826 }
827 
828 /*H2_ReadDissocEnergies read energies for all electronic levels */
829 void diatomics::H2_ReadDissocEnergies( void )
830 {
831  const char* cdDATAFILE = "energy_dissoc.dat";
832  FILE *ioDATA;
833  char chLine[FILENAME_PATH_LENGTH_2];
834  long int i, n1, n2, n3;
835  bool lgEOL;
836 
837  DEBUG_ENTRY( "H2_ReadDissocEnergies()" );
838 
839  /* now open the data file */
840  char chPath[FILENAME_PATH_LENGTH_2];
841  strcpy( chPath, path.c_str() );
842  strcat( chPath, input.chDelimiter );
843  strcat( chPath, cdDATAFILE );
844  ioDATA = open_data( chPath, "r" );
845 
846  /* read the first line and check that magic number is ok */
847  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
848  {
849  fprintf( ioQQQ, " H2_ReadDissocEnergies could not read first line of %s\n", cdDATAFILE );
850  cdEXIT(EXIT_FAILURE);
851  }
852  i = 1;
853  /* magic number */
854  n1 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
855  n2 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
856  n3 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
857 
858  /* magic number
859  * the following is the set of numbers that appear at the start of level1.dat 01 08 10 */
860  if( ( n1 != 2 ) || ( n2 != 4 ) || ( n3 != 29 ) )
861  {
862  fprintf( ioQQQ,
863  " H2_ReadDissocEnergies: the version of %s is not the current version.\n", cdDATAFILE );
864  fprintf( ioQQQ,
865  " I expected to find the number 2 4 29 and got %li %li %li instead.\n" ,
866  n1 , n2 , n3 );
867  fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
868  cdEXIT(EXIT_FAILURE);
869  }
870 
871  while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
872  {
873  /* skip comment */
874  if( chLine[0]=='#' )
875  continue;
876  if( chLine[0]=='\n' || chLine[0]=='\0' || chLine[0]==' ' )
877  break;
878  long iElec;
879  double energyWN;
880  int n = sscanf(chLine,"%li\t%le", &iElec, &energyWN );
881  ASSERT( n == 2 );
882  ASSERT( iElec >= 0 );
883  ASSERT( iElec < N_ELEC );
884  ASSERT( energyWN > 0. );
885  H2_DissocEnergies[iElec] = energyWN;
886  }
887  fclose( ioDATA );
888 
889  return;
890 }
891 
892 /*H2_ReadDissprob read dissociation probabilities and kinetic energies for all electronic levels */
893 void diatomics::H2_ReadDissprob( long int nelec )
894 {
895  const char* cdDATAFILE[N_ELEC] =
896  {
897  "dissprob_X.dat",/* this does not exist and nelec == 0 is not valid */
898  "dissprob_B.dat",
899  "dissprob_C_plus.dat",
900  "dissprob_C_minus.dat",
901  "dissprob_B_primed.dat",
902  "dissprob_D_plus.dat",
903  "dissprob_D_minus.dat"
904  };
905  char chLine[FILENAME_PATH_LENGTH_2];
906  bool lgEOL;
907 
908  DEBUG_ENTRY( "H2_ReadDissprob()" );
909 
910  ASSERT( nelec > 0 );
911 
912  /* now open the data file */
913  char chPath[FILENAME_PATH_LENGTH_2];
914  strcpy( chPath, path.c_str() );
915  strcat( chPath, input.chDelimiter );
916  strcat( chPath, cdDATAFILE[nelec] );
917  FILE *ioDATA = open_data( chPath, "r" );
918 
919  /* read the first line and check that magic number is ok */
920  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
921  {
922  fprintf( ioQQQ, " H2_ReadDissprob could not read first line of %s\n", cdDATAFILE[nelec]);
923  cdEXIT(EXIT_FAILURE);
924  }
925  long i = 1;
926  /* magic number */
927  long n1 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
928  long n2 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
929  long n3 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
930 
931  /* magic number
932  * the following is the set of numbers that appear at the start of level1.dat 01 08 10 */
933  if( ( n1 != 3 ) || ( n2 != 2 ) || ( n3 != 11 ) )
934  {
935  fprintf( ioQQQ,
936  " H2_ReadDissprob: the version of %s is not the current version.\n", cdDATAFILE[nelec] );
937  fprintf( ioQQQ,
938  " I expected to find the number 3 2 11 and got %li %li %li instead.\n" ,
939  n1 , n2 , n3 );
940  fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
941  cdEXIT(EXIT_FAILURE);
942  }
943 
944  while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
945  {
946  /* skip comment */
947  if( chLine[0]=='#' )
948  continue;
949  if( chLine[0]=='\n' || chLine[0]=='\0' || chLine[0]==' ' )
950  break;
951 
952  long iVib, iRot;
953  double a, b;
954  i = 1;
955  sscanf(chLine,"%li\t%li\t%le\t%le",
956  &iVib, &iRot,
957  /* dissociation probability */
958  &a ,
959  /* dissociation kinetic energy - eV not ergs */
960  &b);
961 
962  /* these have to agree if data file is valid */
963  //ASSERT( iVib >= 0 );
964  //ASSERT( iVib <= nVib_hi[nelec] );
965  //ASSERT( iRot >= Jlowest[nelec] );
966  //ASSERT( iRot <= nRot_hi[nelec][iVib] );
967  if( ( iVib < 0 ) ||
968  ( iVib > nVib_hi[nelec] ) ||
969  ( iRot < Jlowest[nelec] ) ||
970  ( iRot > nRot_hi[nelec][iVib] ) )
971  continue;
972 
973  /* dissociation probability */
974  H2_dissprob[nelec][iVib][iRot] = (realnum)a;
975  /* dissociation kinetic energy - eV not ergs */
976  H2_disske[nelec][iVib][iRot] = (realnum)b;
977  }
978  fclose( ioDATA );
979  return;
980 }
981 
982 
983 /*H2_Read_hminus_distribution read distribution function for H2 population following formation from H minus */
984 void diatomics::H2_Read_hminus_distribution(void)
985 {
986  FILE *ioDATA;
987  char chLine[FILENAME_PATH_LENGTH_2];
988  long int i, n1, n2, n3, iVib , iRot;
989  bool lgEOL;
990  double sumrate[nTE_HMINUS] = {0};
991  /* set true for lots of printout */
992  const bool lgH2HMINUS_PRT = false;
993 
994  DEBUG_ENTRY( "H2_Read_hminus_distribution()" );
995 
996  /* now open the data file */
997  char chPath[FILENAME_PATH_LENGTH_2];
998  strcpy( chPath, path.c_str() );
999  strcat( chPath, input.chDelimiter );
1000  strcat( chPath, "hminus_deposit.dat" );
1001  ioDATA = open_data( chPath, "r" );
1002 
1003  /* read the first line and check that magic number is ok */
1004  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
1005  {
1006  fprintf( ioQQQ, " H2_Read_hminus_distribution could not read first line of %s\n", chPath );
1007  cdEXIT(EXIT_FAILURE);
1008  }
1009 
1010  i = 1;
1011  /* magic number */
1012  n1 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
1013  n2 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
1014  n3 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
1015 
1016  /* magic number
1017  * the following is the set of numbers that appear at the start of H2_hminus_deposit.dat 01 08 10 */
1018  if( ( n1 != 2 ) || ( n2 != 10 ) || ( n3 != 17 ) )
1019  {
1020  fprintf( ioQQQ,
1021  " H2_Read_hminus_distribution: the version of %s is not the current version.\n", chPath );
1022  fprintf( ioQQQ,
1023  " I expected to find the number 2 10 17 and got %li %li %li instead.\n" ,
1024  n1 , n2 , n3 );
1025  fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
1026  cdEXIT(EXIT_FAILURE);
1027  }
1028 
1029  /* read until not a comment */
1030  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
1031  BadRead();
1032 
1033  while( chLine[0]=='#' )
1034  {
1035  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
1036  BadRead();
1037  }
1038 
1039  iRot = 1;
1040  iVib = 1;
1041  while( iVib >= 0 )
1042  {
1043  /* set true to print rates */
1044 
1045  double a[nTE_HMINUS] , ener;
1046  sscanf(chLine,"%li\t%li\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf",
1047  &iVib ,&iRot , &ener, &a[0],&a[1],&a[2] , &a[3],&a[4],&a[5] ,&a[6]
1048  );
1049  /* negative iVib says end of data */
1050  if( iVib < 0 )
1051  continue;
1052 
1053  /* check that we actually included the levels in the model representation */
1054  ASSERT( iVib <= nVib_hi[0] &&
1055  iRot <= nRot_hi[0][iVib] );
1056 
1057  if( lgH2HMINUS_PRT )
1058  fprintf(ioQQQ,"hminusss\t%li\t%li", iVib , iRot );
1059  for( i=0; i<nTE_HMINUS; ++i )
1060  {
1061  H2_X_hminus_formation_distribution[i][iVib][iRot] = (realnum)pow(10.,-a[i]);
1062  sumrate[i] += H2_X_hminus_formation_distribution[i][iVib][iRot];
1063  if( lgH2HMINUS_PRT )
1064  fprintf(ioQQQ,"\t%.3e", H2_X_hminus_formation_distribution[i][iVib][iRot] );
1065  }
1066  if( lgH2HMINUS_PRT )
1067  fprintf(ioQQQ,"\n" );
1068 
1069  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
1070  BadRead();
1071  while( chLine[0]=='#' )
1072  {
1073  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
1074  BadRead();
1075  }
1076  }
1077  fclose( ioDATA );
1078 
1079  if( lgH2HMINUS_PRT )
1080  {
1081  /* print total rate */
1082  fprintf(ioQQQ," total H- formation rate ");
1083  /* convert temps to log */
1084  for(i=0; i<nTE_HMINUS; ++i )
1085  {
1086  fprintf(ioQQQ,"\t%.3e" , sumrate[i]);
1087  }
1088  fprintf(ioQQQ,"\n" );
1089  }
1090 
1091  /* convert to dimensionless factors that add to unity */
1092  for( iVib=0; iVib<=nVib_hi[0]; ++iVib )
1093  {
1094  for( iRot=Jlowest[0]; iRot<=nRot_hi[0][iVib]; ++iRot )
1095  {
1096  for(i=0; i<nTE_HMINUS; ++i )
1097  {
1098  H2_X_hminus_formation_distribution[i][iVib][iRot] /= (realnum)sumrate[i];
1099  }
1100  }
1101  }
1102 
1103  if( lgH2HMINUS_PRT )
1104  {
1105  /* print total rate */
1106  fprintf(ioQQQ," H- distribution function ");
1107  for( iVib=0; iVib<=nVib_hi[0]; ++iVib )
1108  {
1109  for( iRot=Jlowest[0]; iRot<=nRot_hi[0][iVib]; ++iRot )
1110  {
1111  fprintf(ioQQQ,"%li\t%li", iVib , iRot );
1112  for(i=0; i<nTE_HMINUS; ++i )
1113  {
1114  fprintf(ioQQQ,"\t%.3e", H2_X_hminus_formation_distribution[i][iVib][iRot] );
1115  }
1116  fprintf(ioQQQ,"\n" );
1117  }
1118  }
1119  }
1120  return;
1121 }
1122 
1123 /* ===================================================================== */
1124 /*H2_Punch_line_data save line data for H2 molecule */
1125 void diatomics::H2_Punch_line_data(
1126  /* io unit for save */
1127  FILE* ioPUN ,
1128  /* save all levels if true, only subset if false */
1129  bool lgDoAll )
1130 {
1131  if( !lgEnabled )
1132  return;
1133 
1134  DEBUG_ENTRY( "H2_Punch_line_data()" );
1135 
1136  if( lgDoAll )
1137  {
1138  fprintf( ioQQQ,
1139  " H2_Punch_line_data ALL option not implemented in H2_Punch_line_data yet 1\n" );
1140  cdEXIT(EXIT_FAILURE);
1141  }
1142  else
1143  {
1144  bool lgPrint = false;
1145  fprintf( ioPUN, "#Eu\tVu\tJu\tEl\tVl\tJl\tWL\tgl\tgu\tgf\tA\tCS\tn(crt)\n" );
1146  /* save line date, looping over all possible lines */
1147  for( TransitionList::iterator tr = trans.begin(); tr != trans.end(); ++tr )
1148  {
1149  if( (*tr).ipCont() <= 0 )
1150  continue;
1151  (*tr).Coll().col_str() = 0.;
1152  qList::iterator Hi = ( (*tr).Hi() );
1153  qList::iterator Lo = ( (*tr).Lo() );
1154  /* print quantum indices */
1155  fprintf(ioPUN,"%2li\t%2li\t%2li\t%2li\t%2li\t%2li\t",
1156  (*Hi).n(), (*Hi).v(), (*Hi).J(),
1157  (*Lo).n(), (*Lo).v(), (*Lo).J() );
1158  Save1LineData( *tr, ioPUN, false , lgPrint);
1159  }
1160 
1161  fprintf( ioPUN , "\n");
1162  }
1163  return;
1164 }
1165 
1166 /*H2_PunchLineStuff include H2 lines in punched optical depths, etc, called from SaveLineStuff */
1167 void diatomics::H2_PunchLineStuff( FILE * io , realnum xLimit , long index)
1168 {
1169  if( !lgEnabled )
1170  return;
1171 
1172  DEBUG_ENTRY( "H2_PunchLineStuff()" );
1173 
1174  /* loop over all possible lines */
1175  for( TransitionList::iterator tr = trans.begin(); tr != trans.end(); ++tr )
1176  {
1177  if( (*tr).ipCont() <= 0 )
1178  continue;
1179  Save1Line( *tr, io, xLimit, index, GetDopplerWidth(2.f*dense.AtomicWeight[ipHYDROGEN]));
1180  }
1181 
1182  return;
1183 }
1184 
1185 
1186 /*H2_Prt_line_tau print line optical depths, called from premet in response to print line optical depths command*/
1187 void diatomics::H2_Prt_line_tau(void)
1188 {
1189  if( !lgEnabled )
1190  return;
1191 
1192  DEBUG_ENTRY( "H2_Prt_line_tau()" );
1193 
1194  /* loop over all possible lines */
1195  for( TransitionList::iterator tr = trans.begin(); tr != trans.end(); ++tr )
1196  {
1197  if( (*tr).ipCont() <= 0 )
1198  continue;
1199  prme( false, *tr );
1200  }
1201 
1202  return;
1203 }
1204 
1205 
1206 /*chMolBranch returns a char with the spectroscopic branch of a transition */
1207 STATIC char chMolBranch( long iRotHi , long int iRotLo )
1208 {
1209  /* these are the spectroscopic branches */
1210  char chBranch[5] = {'O','P','Q','R','S'};
1211  /* this is the index within the chBranch array */
1212  int ip = 2 + (iRotHi - iRotLo);
1213  if( ip<0 || ip>=5 )
1214  {
1215  fprintf(ioQQQ," chMolBranch called with insane iRotHi=%li iRotLo=%li ip=%i\n",
1216  iRotHi , iRotLo , ip );
1217  ip = 0;
1218  }
1219 
1220  return( chBranch[ip] );
1221 }
1222 
1223 /*H2_PunchDo save some properties of the large H2 molecule */
1224 void diatomics::H2_PunchDo( FILE* io , char chJOB[] , const char chTime[] , long int ipPun )
1225 {
1226  DEBUG_ENTRY( "H2_PunchDo()" );
1227 
1228  /* which job are we supposed to do? This routine is active even when H2 is not turned on
1229  * so do not test on lgEnabled initially */
1230 
1231  /* H2 populations computed in last zone -
1232  * give all of molecule in either matrix or triplet format */
1233  if( (strcmp( chJOB , "H2po" ) == 0) && (strcmp(chTime,"LAST") == 0) &&
1234  (save.punarg[ipPun][2] != 0) )
1235  {
1236  /* >>chng 04 feb 19, do not save if H2 not yet evaluated */
1237  if( lgEnabled && lgEvaluated )
1238  {
1239  long iVibHi= 0;
1240  long iRotHi = 0;
1241  long iElecHi=0;
1242  long LimVib, LimRot;
1243  /* the limit to the number of vibration levels punched -
1244  * default is all, but first two numbers on save h2 pops command
1245  * reset limit */
1246  /* this is limit to vibration */
1247  if( save.punarg[ipPun][0] > 0 )
1248  {
1249  LimVib = (long)save.punarg[ipPun][0];
1250  }
1251  else
1252  {
1253  LimVib = nVib_hi[iElecHi];
1254  }
1255 
1256  /* first save the current ortho, para, and total H2 density */
1257  fprintf(io,"%i\t%i\t%.3e\tortho\n",
1258  103 ,
1259  103 ,
1260  ortho_density );
1261  fprintf(io,"%i\t%i\t%.3e\tpara\n",
1262  101 ,
1263  101 ,
1264  para_density );
1265  fprintf(io,"%i\t%i\t%.3e\ttotal\n",
1266  0 ,
1267  0 ,
1268  (*dense_total) );
1269 
1270  /* now save the actual populations, first part both matrix and triplets */
1271  for( iVibHi=0; iVibHi<=LimVib; ++iVibHi )
1272  {
1273  /* this is limit to rotation quantum index */
1274  if( save.punarg[ipPun][1] > 0 )
1275  {
1276  LimRot = (long)MIN2(
1277  save.punarg[ipPun][1] , (realnum)nRot_hi[iElecHi][iVibHi]);
1278  }
1279  else
1280  {
1281  LimRot = nRot_hi[iElecHi][iVibHi];
1282  }
1283  if( save.punarg[ipPun][2] > 0 )
1284  {
1285  long int i;
1286  /* this option save matrix */
1287  if( iVibHi == 0 )
1288  {
1289  fprintf(io,"vib\\rot");
1290  /* this is first vib, so make row of rot numbs */
1291  for( i=0; i<=LimRot; ++i )
1292  {
1293  fprintf(io,"\t%li",i);
1294  }
1295  fprintf(io,"\n");
1296  }
1297  fprintf(io,"%li",iVibHi );
1298  for( iRotHi=Jlowest[iElecHi]; iRotHi<=LimRot; ++iRotHi )
1299  {
1300  fprintf(io,"\t%.3e",
1301  states[ ipEnergySort[iElecHi][iVibHi][iRotHi] ].Pop()/(*dense_total) );
1302  }
1303  fprintf(io,"\n" );
1304  }
1305  else if( save.punarg[ipPun][2] < 0 )
1306  {
1307  /* this option save triplets - the default */
1308  for( iRotHi=Jlowest[iElecHi]; iRotHi<=LimRot; ++iRotHi )
1309  {
1310  /* this will say whether ortho or para,
1311  * H2_lgOrtho is 0 or 1 depending on whether or not ortho,
1312  * so chlgPara[H2_lgOrtho] gives P or O for printing */
1313  const char chlgPara[2]={'P','O'};
1314  const long ipHi = ipEnergySort[iElecHi][iVibHi][iRotHi];
1315 
1316  /* intensity, relative to normalization line, for faintest line to save */
1317  fprintf(io,"%li\t%li\t%c\t%.1f\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
1318  /* upper vibration and rotation quantum numbers */
1319  iVibHi , iRotHi ,
1320  /* an 'O' or 'P' for ortho or para */
1321  chlgPara[H2_lgOrtho[iElecHi][iVibHi][iRotHi]],
1322  /* the level excitation energy in wavenumbers */
1323  states[ipHi].energy().WN(),
1324  /* actual population relative to total H2 */
1325  states[ipHi].Pop()/(*dense_total) ,
1326  /* old level populations for comparison */
1327  H2_old_populations[iElecHi][iVibHi][iRotHi]/(*dense_total) ,
1328  /* populations per h2 and per statistical weight */
1329  states[ipHi].Pop()/(*dense_total)/H2_stat[iElecHi][iVibHi][iRotHi] ,
1330  /* LTE departure coefficient */
1331  /* >>chng 05 jan 26, missing factor of H2 abundance LTE is norm to unity, not tot abund */
1332  states[ipHi].Pop()/SDIV(H2_populations_LTE[iElecHi][iVibHi][iRotHi]*(*dense_total) ) ,
1333  /* fraction of exits that were collisional */
1334  H2_col_rate_out[iVibHi][iRotHi]/SDIV(H2_col_rate_out[iVibHi][iRotHi]+H2_rad_rate_out[0][iVibHi][iRotHi]) ,
1335  /* fraction of entries that were collisional */
1336  H2_col_rate_in[iVibHi][iRotHi]/SDIV(H2_col_rate_in[iVibHi][iRotHi]+H2_rad_rate_in[iVibHi][iRotHi]),
1337  /* collisions out */
1338  H2_col_rate_out[iVibHi][iRotHi],
1339  /* radiation out */
1340  H2_rad_rate_out[0][iVibHi][iRotHi] ,
1341  /* radiation out */
1342  H2_col_rate_in[iVibHi][iRotHi],
1343  /* radiation in */
1344  H2_rad_rate_in[iVibHi][iRotHi]
1345  );
1346  }
1347  }
1348  }
1349  }
1350  }
1351  /* save H2 populations for each zone
1352  * populations of v=0 for each zone */
1353  else if( (strcmp( chJOB , "H2po" ) == 0) && (strcmp(chTime,"LAST") != 0) &&
1354  (save.punarg[ipPun][2] == 0) )
1355  {
1356  /* >>chng 04 feb 19, do not save if h2 not yet evaluated */
1357  if( lgEnabled && lgEvaluated )
1358  {
1359  fprintf(io,"%.5e\t%.3e\t%.3e", radius.depth_mid_zone ,
1360  ortho_density , para_density);
1361  /* rel pops of first two excited electronic states */
1362  fprintf(io,"\t%.3e\t%.3e",
1363  pops_per_elec[1] , pops_per_elec[2]);
1364  long iElecHi = 0;
1365  long iVibHi = 0;
1366  long LimVib, LimRot;
1367  /* this is limit to vibration quantum index */
1368  if( save.punarg[ipPun][0] > 0 )
1369  {
1370  LimVib = (long)save.punarg[ipPun][1];
1371  }
1372  else
1373  {
1374  LimVib = nRot_hi[iElecHi][iVibHi];
1375  }
1376  LimVib = MIN2( LimVib , nVib_hi[iElecHi] );
1377  /* this is limit to rotation quantum index */
1378  if( save.punarg[ipPun][1] > 0 )
1379  {
1380  LimRot = (long)save.punarg[ipPun][1];
1381  }
1382  else
1383  {
1384  LimRot = nRot_hi[iElecHi][iVibHi];
1385  }
1386  for( iVibHi = 0; iVibHi<=LimVib; ++iVibHi )
1387  {
1388  fprintf(io,"\tv=%li",iVibHi);
1389  long int LimRotVib = MIN2( LimRot , nRot_hi[iElecHi][iVibHi] );
1390  for( long iRotHi=Jlowest[iElecHi]; iRotHi<=LimRotVib; ++iRotHi )
1391  {
1392  fprintf(io,"\t%.3e",
1393  states[ ipEnergySort[iElecHi][iVibHi][iRotHi] ].Pop()/(*dense_total) );
1394  }
1395  }
1396  fprintf(io,"\n");
1397  }
1398  }
1399 
1400  /* save column densities */
1401  else if( (strcmp( chJOB , "H2cl" ) == 0) && (strcmp(chTime,"LAST") == 0) )
1402  {
1403  long iVibHi= 0;
1404  long iRotHi = 0;
1405  long iElecHi=0;
1406  long LimVib, LimRot;
1407  /* the limit to the number of vibration levels punched -
1408  * default is all, but first two numbers on save h2 pops command
1409  * reset limit */
1410  /* this is limit to vibration */
1411  if( save.punarg[ipPun][0] > 0 )
1412  {
1413  LimVib = (long)save.punarg[ipPun][0];
1414  }
1415  else
1416  {
1417  LimVib = nVib_hi[iElecHi];
1418  }
1419 
1420  /* first save ortho and para populations */
1421  fprintf(io,"%i\t%i\t%.3e\tortho\n",
1422  103 ,
1423  103 ,
1424  ortho_colden );
1425  fprintf(io,"%i\t%i\t%.3e\tpara\n",
1426  101 ,
1427  101 ,
1428  para_colden );
1429  /* total H2 column density */
1430  fprintf(io,"%i\t%i\t%.3e\ttotal\n",
1431  0 ,
1432  0 ,
1433  colden.colden[ipCOL_H2g]+colden.colden[ipCOL_H2s]);
1434 
1435  /* save level column densities */
1436  for( iVibHi=0; iVibHi<=LimVib; ++iVibHi )
1437  {
1438  if( lgEnabled )
1439  {
1440  /* this is limit to rotation quantum index */
1441  if( save.punarg[ipPun][1] > 0 )
1442  {
1443  LimRot = (long)save.punarg[ipPun][1];
1444  }
1445  else
1446  {
1447  LimRot = nRot_hi[iElecHi][iVibHi];
1448  }
1449  if( save.punarg[ipPun][2] > 0 )
1450  {
1451  long int i;
1452  /* save matrix */
1453  if( iVibHi == 0 )
1454  {
1455  fprintf(io,"vib\\rot");
1456  /* this is first vib, so make row of rot numbs */
1457  for( i=0; i<=LimRot; ++i )
1458  {
1459  fprintf(io,"\t%li",i);
1460  }
1461  fprintf(io,"\n");
1462  }
1463  fprintf(io,"%li",iVibHi );
1464  for( iRotHi=Jlowest[iElecHi]; iRotHi<=LimRot; ++iRotHi )
1465  {
1466  fprintf(io,"\t%.3e",
1467  H2_X_colden[iVibHi][iRotHi]/(*dense_total) );
1468  }
1469  fprintf(io,"\n" );
1470  }
1471  else
1472  {
1473  /* save triplets - the default */
1474  for( iRotHi=Jlowest[iElecHi]; iRotHi<=LimRot; ++iRotHi )
1475  {
1476  fprintf(io,"%li\t%li\t%.1f\t%.3e\t%.3e\t%.3e\t%.3e\n",
1477  iVibHi ,
1478  iRotHi ,
1479  /* energy relative to 0,0, T1CM converts wavenumber to K */
1480  states[ ipEnergySort[iElecHi][iVibHi][iRotHi] ].energy().K(),
1481  /* these are column densities for actual molecule */
1482  H2_X_colden[iVibHi][iRotHi] ,
1483  H2_X_colden[iVibHi][iRotHi]/H2_stat[iElecHi][iVibHi][iRotHi] ,
1484  /* these are same column densities but for LTE populations */
1485  H2_X_colden_LTE[iVibHi][iRotHi] ,
1486  H2_X_colden_LTE[iVibHi][iRotHi]/H2_stat[iElecHi][iVibHi][iRotHi]);
1487  }
1488  }
1489  }
1490  }
1491  }
1492  else if( (strcmp(chJOB , "H2pd" ) == 0) && (strcmp(chTime,"LAST") != 0) )
1493  {
1494  /* save PDR
1495  * output some PDR information (densities, rates) for each zone */
1496  fprintf(io,"%.5e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
1497  /* depth in cm */
1498  radius.depth_mid_zone ,
1499  /* the computed ortho and para densities */
1500  ortho_density ,
1501  para_density ,
1502  /* the Lyman Werner band dissociation, Tielens & Hollenbach */
1503  hmi.H2_Solomon_dissoc_rate_TH85_H2g ,
1504  /* the Lyman Werner band dissociation, Bertoldi & Draine */
1505  hmi.H2_Solomon_dissoc_rate_BD96_H2g,
1506  /* the Lyman Werner band dissociation, big H2 mole */
1507  Solomon_dissoc_rate_g);
1508  }
1509  else if( (strcmp(chJOB , "H2co" ) == 0) && (strcmp(chTime,"LAST") != 0) )
1510  {
1511  /* save H2 cooling - do heating cooling for each zone old new H2 */
1512  fprintf(io,"%.5e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
1513  /* depth in cm */
1514  radius.depth_mid_zone ,
1515  /* total cooling, equal to total heating */
1516  thermal.ctot ,
1517  /* H2 destruction by Solomon process, TH85 rate */
1518  hmi.H2_Solomon_dissoc_rate_TH85_H2g,
1519  /* H2 destruction by Solomon process, big H2 model rate */
1520  Solomon_dissoc_rate_g +
1521  Solomon_dissoc_rate_s,
1522  /* H2 photodissociation heating, eqn A9 of Tielens & Hollenbach 1985a */
1523  hmi.HeatH2Dish_TH85,
1524  /* heating due to dissociation of electronic excited states */
1525  HeatDiss ,
1526  /* cooling (usually neg and so heating) due to collisions within X */
1527  hmi.HeatH2Dexc_TH85,
1528  HeatDexc
1529  );
1530 
1531  }
1532 
1533  else if( (strcmp(chJOB , "H2le" ) == 0) && (strcmp(chTime,"LAST") == 0) )
1534  {
1535  /* save H2 levels */
1536  for( long int ipHi=0; ipHi < nLevels_per_elec[0]; ipHi++ )
1537  {
1538  long iRotHi = ipRot_H2_energy_sort[ipHi];
1539  long iVibHi = ipVib_H2_energy_sort[ipHi];
1540  long int nColl;
1541  double Asum , Csum[N_X_COLLIDER];
1542  Asum = 0;
1543  for( nColl=0; nColl<N_X_COLLIDER; ++nColl )
1544  Csum[nColl] = 0.;
1545  for( long int ipLo=0; ipLo<ipHi; ++ipLo )
1546  {
1547  /* all lower levels */
1548  long iRotLo = ipRot_H2_energy_sort[ipLo];
1549  long iVibLo = ipVib_H2_energy_sort[ipLo];
1550  EmissionProxy em = trans[ ipTransitionSort[ipHi][ipLo] ].Emis();
1551 
1552  /* radiative decays down */
1553  if( ( abs(iRotHi-iRotLo) == 2 || (iRotHi-iRotLo) == 0 ) && iVibLo <= iVibHi &&
1554  lgH2_radiative[ipHi][ipLo] )
1555  {
1556  Asum += em.Aul() * ( em.Pesc() + em.Pdest() + em.Pelec_esc() );
1557  }
1558  /* all collisions down */
1559  mr3ci H2cr = CollRateCoeff.begin(ipHi,ipLo);
1560  for( nColl=0; nColl<N_X_COLLIDER; ++nColl )
1561  Csum[nColl] += H2cr[nColl];
1562  }
1563 
1564  /* save H2 level energies */
1565  fprintf(io,"%li\t%li\t%.2f\t%li\t%.3e",
1566  iVibHi , iRotHi,
1567  states[ipHi].energy().WN(),
1568  (long)states[ipHi].g(),
1569  Asum );
1570  for( nColl=0; nColl<N_X_COLLIDER; ++nColl )
1571  /* sum over all lower levels */
1572  fprintf(io,"\t%.3e",Csum[nColl]);
1573  fprintf(io,"\n");
1574  }
1575  }
1576 
1577  else if( (strcmp(chJOB , "H2ra" ) == 0) && (strcmp(chTime,"LAST") != 0) )
1578  {
1579  /* save h2 rates - some rates and lifetimes */
1580  double sumpop = 0. , sumlife = 0.;
1581 
1582  /* this block, find lifetime against photo excitation into excited electronic states */
1583  if( lgEnabled && lgEvaluated )
1584  {
1585  /* only do if radiative transition exists */
1586  for( TransitionList::iterator tr = trans.begin(); tr != rad_end; ++tr )
1587  {
1588  qList::iterator Lo = ( (*tr).Lo() );
1589  if( (*Lo).n() > 0 || (*Lo).v() > 0 )
1590  continue;
1591  sumlife += (*tr).Emis().pump() * (*(*tr).Lo()).Pop();
1592  sumpop += (*(*tr).Lo()).Pop();
1593  }
1594  }
1595 
1596  /* continue output from save h2 rates command */
1597  /* find photoexcitation rates from v=0 */
1598  /* PDR information for each zone */
1599  fprintf(io,
1600  "%.5e\t%.3e\t%.3e\t%.3e\t%li",
1601  /* depth in cm */
1602  radius.depth_mid_zone ,
1603  /* the column density (cm^-2) in H2 */
1604  colden.colden[ipCOL_H2g]+colden.colden[ipCOL_H2s],
1605  /* this is a special form of column density - should be proportional
1606  * to total shielding */
1607  colden.coldenH2_ov_vel ,
1608  /* visual extinction due to dust alone, of point source (star)*/
1609  rfield.extin_mag_V_point,
1610  /* number of large molecule evaluations in this zone */
1611  nCall_this_zone );
1612  fprintf(io,
1613  "\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e",
1614  /* total H2 fraction */
1615  (*dense_total)/dense.gas_phase[ipHYDROGEN] ,
1616  /* chemistry renorm factor */
1617  H2_renorm_chemistry,
1618  /* rate H2 forms on grains */
1619  gv.rate_h2_form_grains_used_total ,
1620  /* rate H2 forms by H minus route */
1621  findspecieslocal("H-")->den*1.35e-9,
1622  /* H2 destruction by Solomon process, TH85 rate */
1623  hmi.H2_Solomon_dissoc_rate_TH85_H2g + hmi.H2_Solomon_dissoc_rate_TH85_H2s,
1624  /* H2 destruction by Solomon process, Bertoldi & Draine rate */
1625  hmi.H2_Solomon_dissoc_rate_BD96_H2g + hmi.H2_Solomon_dissoc_rate_BD96_H2s,
1626  /* H2 destruction by Solomon process, Elwert et al. in preparation */
1627  hmi.H2_Solomon_dissoc_rate_ELWERT_H2g + hmi.H2_Solomon_dissoc_rate_ELWERT_H2g,
1628  /* H2 destruction by Solomon process, big H2 model rate */
1629  Solomon_dissoc_rate_g + Solomon_dissoc_rate_s,
1630  /* rate s-1 H2 electronic excit states decay into H2g */
1631  Solomon_elec_decay_g ,
1632  /* rate s-1 H2 electronic excit states decay into H2s */
1633  Solomon_elec_decay_s
1634  );
1635  fprintf(io,
1636  "\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e",
1637  /* The TH85 estimate of the radiation field relative to the Habing value */
1638  hmi.UV_Cont_rel2_Habing_TH85_depth,
1639  /* The DB96 estimate of the radiation field relative to the Habing value */
1640  hmi.UV_Cont_rel2_Draine_DB96_depth,
1641  /* cosmic ray ionization rate */
1642  secondaries.csupra[ipHYDROGEN][0]*0.93,
1643  sumlife/SDIV( sumpop ) ,
1644  hmi.H2_Solomon_dissoc_rate_BD96_H2g/SDIV(hmi.UV_Cont_rel2_Habing_TH85_depth) ,
1645  Solomon_dissoc_rate_g/SDIV(hmi.UV_Cont_rel2_Habing_TH85_depth),
1646  Solomon_dissoc_rate_g/SDIV(hmi.UV_Cont_rel2_Habing_spec_depth),
1647  hmi.H2_rate_destroy);
1648  fprintf(io,
1649  "\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
1650  hmi.HeatH2Dish_TH85,
1651  hmi.HeatH2Dexc_TH85,
1652  HeatDiss,
1653  HeatDexc,
1654  thermal.htot);
1655  }
1656  /* save h2 solomon */
1657  else if( (strcmp(chJOB , "H2so" ) == 0) && (strcmp(chTime,"LAST") != 0) )
1658  {
1659  /* remember as many as nSOL lines contributing to total Solomon process */
1660  const int nSOL = 100;
1661  double sum, one;
1662  long int jlosave[nSOL] , ivlosave[nSOL],
1663  iehisave[nSOL] ,jhisave[nSOL] , ivhisave[nSOL],
1664  nsave,
1665  ipOrdered[nSOL];
1666  int nFail;
1667  realnum fsave[nSOL], wlsave[nSOL];
1668  /* Solomon process, and where it came from */
1669  fprintf(io,"%.5e\t%.3e",
1670  /* depth in cm */
1671  radius.depth_mid_zone ,
1672  /* H2 destruction by Solomon process, big H2 model rate */
1673  Solomon_dissoc_rate_g +
1674  Solomon_dissoc_rate_s);
1675  sum = 0.;
1676  /* find sum of all radiative exits from X into excited electronic states */
1677  if( lgEnabled && lgEvaluated )
1678  {
1679  for( TransitionList::iterator tr = trans.begin(); tr != rad_end; ++tr )
1680  {
1681  qList::iterator Lo = ( (*tr).Lo() );
1682  if( (*Lo).n() > 0 )
1683  continue;
1684  sum += (*(*tr).Lo()).Pop() * (*tr).Emis().pump();
1685  }
1686 
1687  /* make sure it is safe to div by sum */
1688  sum = SDIV( sum );
1689  nsave = 0;
1690  /* now loop back over X and print all those which contribute more than frac of the total */
1691  const double frac = 0.01;
1692  /* only do if radiative transition exists */
1693  for( TransitionList::iterator tr = trans.begin(); tr != rad_end; ++tr )
1694  {
1695  qList::iterator Lo = ( (*tr).Lo() );
1696  if( (*Lo).n() > 0 )
1697  continue;
1698  one = (*(*tr).Lo()).Pop() * (*tr).Emis().pump();
1699  if( one/sum > frac && nsave<nSOL)
1700  {
1701  qList::iterator Hi = ( (*tr).Hi() );
1702  fsave[nsave] = (realnum)(one/sum);
1703  jlosave[nsave] = (*Lo).J();
1704  ivlosave[nsave] = (*Lo).v();
1705  jhisave[nsave] = (*Hi).J();
1706  ivhisave[nsave] = (*Hi).v();
1707  iehisave[nsave] = (*Hi).n();
1708  wlsave[nsave] = (*tr).WLAng();
1709  ++nsave;
1710  }
1711  }
1712 
1713  /* now sort by decreasing importance */
1714  /*spsort netlib routine to sort array returning sorted indices */
1715  spsort(
1716  /* input array to be sorted */
1717  fsave,
1718  /* number of values in x */
1719  nsave,
1720  /* permutation output array */
1721  ipOrdered,
1722  /* flag saying what to do - 1 sorts into increasing order, not changing
1723  * the original routine */
1724  -1,
1725  /* error condition, should be 0 */
1726  &nFail);
1727 
1728  /* print ratio of pumps to dissociations - this is 9:1 in TH85 */
1729  /*>>chng 05 jul 20, TE, save average energy in H2s and renormalization factors for H2g and H2s */
1730  /* >>chng 05 sep 16, TE, chng denominator to do g and s with proper dissoc rates */
1731  fprintf(io,"\t%.3f\t%.3f\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e",
1732  /* this is sum of photons and CRs */
1733  (sum + secondaries.csupra[ipHYDROGEN][0]*2.02f)/SDIV((Solomon_dissoc_rate_g * findspecieslocal("H2")->den +
1734  Solomon_dissoc_rate_s * findspecieslocal("H2*")->den) ),
1735  /* this is sum of photons and CRs */
1736  (sum + secondaries.csupra[ipHYDROGEN][0]*2.02f) /SDIV((Solomon_dissoc_rate_g * H2_den_g+
1737  Solomon_dissoc_rate_s * H2_den_s) ),
1738  average_energy_g, average_energy_s,
1739  findspecieslocal("H2")->den/SDIV(H2_den_g), findspecieslocal("H2*")->den/SDIV(H2_den_s),
1740  H2_den_g/SDIV((*dense_total)), H2_den_s/SDIV((*dense_total))
1741  );
1742  for( long i=0; i<nsave; ++i )
1743  {
1744  long ip = ipOrdered[i];
1745  /*lint -e644 not init */
1746  fprintf(io,"\t%li\t%li\t%li\t%li\t%li\t%.3f\t%.3f",
1747  iehisave[ip],ivhisave[ip],jhisave[ip],ivlosave[ip] , jlosave[ip] , fsave[ip] , wlsave[ip] );
1748  /*lint +e644 not init */
1749  }
1750  fprintf(io,"\n");
1751  }
1752  }
1753 
1754  else if( (strcmp(chJOB , "H2te" ) == 0) && (strcmp(chTime,"LAST") != 0) )
1755  {
1756  /* save h2 temperatures */
1757  double pop_ratio10,pop_ratio20,pop_ratio30,pop_ratio31,pop_ratio40;
1758  double T10,T20,T30,T31,T40;
1759  /* subscript"sum" denotes integrated quantities */
1760  double T10_sum,T20_sum,T30_sum,T31_sum,T40_sum;
1761  double pop_ratio10_sum,pop_ratio20_sum,pop_ratio30_sum,pop_ratio31_sum,pop_ratio40_sum;
1762  if( lgEnabled && nCall_this_zone )
1763  {
1764  double pop0 = states[0].Pop();
1765  double pop1 = states[1].Pop();
1766  double pop2 = states[2].Pop();
1767  double pop3 = states[3].Pop();
1768  double pop4 = states[4].Pop();
1769 
1770  double energyK = states[1].energy().K() - states[0].energy().K();
1771  /* the ratio of populations of J=1 to 0 */
1772  pop_ratio10 = pop1/SDIV(pop0);
1773  pop_ratio10_sum = H2_X_colden[0][1]/SDIV(H2_X_colden[0][0]);
1774  /* the corresponding temperature */
1775  T10 = -170.5/log(SDIV(pop_ratio10) * H2_stat[0][0][0]/H2_stat[0][0][1]);
1776  T10_sum = -170.5/log(SDIV(pop_ratio10_sum) * H2_stat[0][0][0]/H2_stat[0][0][1]);
1777 
1778  energyK = states[2].energy().K() - states[0].energy().K();
1779  pop_ratio20 = pop2/SDIV(pop0);
1780  T20 = -energyK/log(SDIV(pop_ratio20) * H2_stat[0][0][0]/H2_stat[0][0][2]);
1781 
1782  pop_ratio20_sum = H2_X_colden[0][2]/SDIV(H2_X_colden[0][0]);
1783  T20_sum = -energyK/log(SDIV(pop_ratio20_sum) * H2_stat[0][0][0]/H2_stat[0][0][2]);
1784 
1785  energyK = states[3].energy().K() - states[0].energy().K();
1786  pop_ratio30 = pop3/SDIV(pop0);
1787  T30 = -energyK/log(SDIV(pop_ratio30) * H2_stat[0][0][0]/H2_stat[0][0][3]);
1788 
1789  pop_ratio30_sum = H2_X_colden[0][3]/SDIV(H2_X_colden[0][0]);
1790  T30_sum = -energyK/log(SDIV(pop_ratio30_sum) * H2_stat[0][0][0]/H2_stat[0][0][3]);
1791 
1792  energyK = states[3].energy().K() - states[1].energy().K();
1793  pop_ratio31 = pop3/SDIV(pop1);
1794  T31 = -energyK/log(SDIV(pop_ratio31) * H2_stat[0][0][1]/H2_stat[0][0][3]);
1795 
1796  pop_ratio31_sum = H2_X_colden[0][3]/SDIV(H2_X_colden[0][1]);
1797  T31_sum = -energyK/log(SDIV(pop_ratio31_sum) * H2_stat[0][0][1]/H2_stat[0][0][3]);
1798 
1799  energyK = states[4].energy().K() - states[0].energy().K();
1800  pop_ratio40 = pop4/SDIV(pop0);
1801  T40 = -energyK/log(SDIV(pop_ratio40) * H2_stat[0][0][0]/H2_stat[0][0][4]);
1802 
1803  pop_ratio40_sum = H2_X_colden[0][4]/SDIV(H2_X_colden[0][0]);
1804  T40_sum = -energyK/log(SDIV(pop_ratio40_sum) * H2_stat[0][0][0]/H2_stat[0][0][4]);
1805  }
1806  else
1807  {
1808  pop_ratio10 = 0.;
1809  pop_ratio10_sum = 0.;
1810  T10 = 0.;
1811  T20 = 0.;
1812  T30 = 0.;
1813  T31 = 0.;
1814  T40 = 0.;
1815  T10_sum = 0.;
1816  T20_sum = 0.;
1817  T30_sum = 0.;
1818  T31_sum = 0.;
1819  T40_sum = 0.;
1820  }
1821 
1822  /* various temperatures for neutral/molecular gas */
1823  fprintf( io,
1824  "%.5e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n" ,
1825  /* depth in cm */
1826  radius.depth_mid_zone ,
1827  /* total H2 fraction */
1828  (*dense_total)/dense.gas_phase[ipHYDROGEN] ,
1829  /* ratio of populations of 1 to 0 only */
1830  pop_ratio10 ,
1831  /* sum of all ortho and para */
1832  ortho_density / SDIV(para_density),
1833  T10,T20,T30,T31,T40,
1834  phycon.te ,
1835  hyperfine.Tspin21cm,T10_sum,T20_sum,T30_sum,T31_sum,T40_sum );
1836  }
1837  else if( (strcmp(chJOB , "H2ln" ) == 0) && (strcmp(chTime,"LAST") == 0) )
1838  {
1839  /* save H2 lines - output the full emission-line spectrum */
1840  double thresh;
1841  double renorm;
1842  /* first test, is H2 turned on? Second test, have lines arrays
1843  * been set up - nsum is negative if abort occurs before lines
1844  * are set up */
1845  if( lgEnabled && LineSave.nsum > 0)
1846  {
1847  ASSERT( LineSave.ipNormWavL >= 0 );
1848  /* get the normalization line */
1849  if( LineSv[LineSave.ipNormWavL].SumLine[0] > SMALLFLOAT )
1850  renorm = LineSave.ScaleNormLine/
1851  LineSv[LineSave.ipNormWavL].SumLine[0];
1852  else
1853  renorm = 1.;
1854 
1855  if( renorm > SMALLFLOAT )
1856  {
1857  /* this is threshold for faintest line, normally 0, set with
1858  * number on save H2 command */
1859  thresh = thresh_punline_h2/(realnum)renorm;
1860  }
1861  else
1862  thresh = 0.f;
1863 
1864  /* save H2 line intensities at end of iteration
1865  * nElecLevelOutput is electronic level with 1 for ground, so this loop is < nElecLevelOutput */
1866  for( TransitionList::iterator tr = trans.begin(); tr != rad_end; ++tr )
1867  {
1868  qList::iterator Hi = ( (*tr).Hi() );
1869  qList::iterator Lo = ( (*tr).Lo() );
1870  long iElecHi = (*Hi).n();
1871  long iVibHi = (*Hi).v();
1872  long iRotHi = (*Hi).J();
1873  long iElecLo = (*Lo).n();
1874  long iVibLo = (*Lo).v();
1875  long iRotLo = (*Lo).J();
1876  if( iElecHi >= nElecLevelOutput )
1877  continue;
1878  if( H2_SaveLine[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] > thresh )
1879  {
1880  /* air wavelength in microns */
1881  /* WLAng contains correction for index of refraction of air */
1882  double wl = (*tr).WLAng()/1e4;
1883  fprintf(io, "%li-%li %c(%li)", iVibHi, iVibLo, chMolBranch( iRotHi, iRotLo ), iRotLo );
1884  fprintf( io, "\t%ld\t%ld\t%ld\t%ld\t%ld\t%ld", iElecHi , iVibHi , iRotHi , iElecLo , iVibLo , iRotLo);
1885  /* WLAng contains correction for index of refraction of air */
1886  fprintf( io, "\t%.7f\t", wl );
1887  /*prt_wl print floating wavelength in Angstroms, in output format */
1888  prt_wl( io , (*tr).WLAng() );
1889  /* the log of the line intensity or luminosity */
1890  fprintf( io, "\t%.3f\t%.3e",
1891  log10(MAX2(1e-37,H2_SaveLine[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo])) + radius.Conv2PrtInten,
1892  H2_SaveLine[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo]*renorm );
1893  /* excitation energy of upper level in K */
1894  fprintf( io, "\t%.3f", (*Hi).energy().K() );
1895  /* the product g_hi h nu * Aul */
1896  fprintf( io, "\t%.3e", (*tr).Emis().Aul() * (*tr).EnergyErg() * (*(*tr).Hi()).g() );
1897  fprintf( io, "\n");
1898  }
1899  }
1900  }
1901  }
1902  else if( (strcmp(chJOB , "H2sp" ) == 0) )
1903  {
1904  fprintf( io, "PUT SOMETHING HERE!\n" );
1905  }
1906  return;
1907 }
1908  /*cdH2_Line determines intensity and luminosity of and H2 line. The first
1909  * six arguments give the upper and lower quantum designation of the levels.
1910  * The function returns 1 if we found the line,
1911  * and false==0 if we did not find the line because ohoto-para transition
1912  * or upper level has lower energy than lower level */
1913 long int cdH2_Line(
1914  /* indices for the upper level */
1915  long int iElecHi,
1916  long int iVibHi ,
1917  long int iRotHi ,
1918  /* indices for lower level */
1919  long int iElecLo,
1920  long int iVibLo ,
1921  long int iRotLo ,
1922  /* linear intensity relative to normalization line*/
1923  double *relint,
1924  /* log of luminosity or intensity of line */
1925  double *absint )
1926 {
1927  DEBUG_ENTRY( "cdH2_Line()" );
1928  return h2.getLine( iElecHi, iVibHi, iRotHi, iElecLo, iVibLo, iRotLo, relint, absint );
1929 };
1930 
1931 long int diatomics::getLine( long iElecHi, long iVibHi, long iRotHi, long iElecLo, long iVibLo, long iRotLo, double *relint, double *absint )
1932 {
1933 
1934  DEBUG_ENTRY( "diatomics::getline()" );
1935 
1936  /* these will be return values if we can't find the line */
1937  *relint = 0.;
1938  *absint = 0.;
1939 
1940  /* for now both electronic levels must be zero */
1941  if( iElecHi!=0 || iElecLo!=0 )
1942  {
1943  return 0;
1944  }
1945 
1946  long ipHi = ipEnergySort[iElecHi][iVibHi][iRotHi];
1947  long ipLo = ipEnergySort[iElecLo][iVibLo][iRotLo];
1948 
1949  /* check that energy of first level is higher than energy of second level */
1950  if( states[ipHi].energy().WN() < states[ipLo].energy().WN() )
1951  {
1952  return 0;
1953  }
1954 
1955  /* check that ortho-para does not change */
1956  if( H2_lgOrtho[iElecHi][iVibHi][iRotHi] - H2_lgOrtho[iElecLo][iVibLo][iRotLo] != 0 )
1957  {
1958  return 0;
1959  }
1960 
1961  /* exit if lines does not exist */
1962  if( !lgH2_radiative[ipHi][ipLo] )
1963  {
1964  return 0;
1965  }
1966 
1967  ASSERT( LineSave.ipNormWavL >= 0 );
1968  /* does the normalization line have a positive intensity*/
1969  if( LineSv[LineSave.ipNormWavL].SumLine[0] > 0. )
1970  {
1971  *relint = H2_SaveLine[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo]/
1972  LineSv[LineSave.ipNormWavL].SumLine[0] * LineSave.ScaleNormLine;
1973  }
1974  else
1975  {
1976  *relint = 0.;
1977  }
1978 
1979  /* return log of line intensity if it is positive */
1980  if( H2_SaveLine[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] > 0. )
1981  {
1982  *absint = log10(H2_SaveLine[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo]) +
1983  radius.Conv2PrtInten;
1984  }
1985  else
1986  {
1987  /* line intensity is actually zero, return small number */
1988  *absint = -37.;
1989  }
1990  /* this indicates success */
1991  return 1;
1992 }
1993 
1994 void diatomics::set_numLevelsMatrix( long numLevels )
1995 {
1996  if( !lgREAD_DATA )
1997  nXLevelsMatrix = numLevels;
1998 }
1999 
diatomics::depart
vector< double > depart
Definition: h2_priv.h:702
colden.h
thermal.h
diatomics::H2_LinesAdd
void H2_LinesAdd(void)
Definition: mole_h2_io.cpp:49
t_hyperfine::Tspin21cm
double Tspin21cm
Definition: hyperfine.h:47
diatomics::nVib_hi
long int nVib_hi[N_ELEC]
Definition: h2_priv.h:611
findspecieslocal
molezone * findspecieslocal(const char buf[])
Definition: mole_species.cpp:833
ipCOL_H2s
#define ipCOL_H2s
Definition: colden.h:18
Parser::nMatch
bool nMatch(const char *chKey) const
Definition: parser.h:135
prme
void prme(const bool lgReset, const TransitionProxy &t)
Definition: prt_met.cpp:97
diatomics::H2_disske
multi_arr< realnum, 3 > H2_disske
Definition: h2_priv.h:633
open_data
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition: cpu.cpp:625
lines.h
h2
diatomics h2("h2", 4100., &hmi.H2_total, Yan_H2_CS)
diatomics::average_energy_g
double average_energy_g
Definition: h2_priv.h:286
Save1LineData
void Save1LineData(const TransitionProxy &t, FILE *io, bool lgCS_2, bool &lgPrint)
Definition: save_linedata.cpp:278
diatomics::H2_PunchDo
void H2_PunchDo(FILE *io, char chJOB[], const char chTime[], long int ipPun)
Definition: mole_h2_io.cpp:1224
diatomics::H2_SaveLine
multi_arr< realnum, 6 > H2_SaveLine
Definition: h2_priv.h:710
dense
t_dense dense
Definition: dense.cpp:24
secondaries.h
BadRead
NORETURN void BadRead(void)
Definition: service.cpp:901
diatomics::pops_per_vib
multi_arr< double, 2 > pops_per_vib
Definition: h2_priv.h:600
diatomics::lgEvaluated
bool lgEvaluated
Definition: h2_priv.h:310
rfield
t_rfield rfield
Definition: rfield.cpp:8
ioQQQ
FILE * ioQQQ
Definition: cddefines.cpp:7
level_tmp
Definition: mole_h2_io.cpp:657
diatomics::lgEnabled
bool lgEnabled
Definition: h2_priv.h:345
FFmtRead
double FFmtRead(const char *chCard, long int *ipnt, long int last, bool *lgEOL)
Definition: service.cpp:381
t_radius::Conv2PrtInten
double Conv2PrtInten
Definition: radius.h:147
realnum
float realnum
Definition: cddefines.h:103
t_radius::dVeffAper
double dVeffAper
Definition: radius.h:87
diatomics::label
string label
Definition: h2_priv.h:571
rfield.h
STATIC
#define STATIC
Definition: cddefines.h:97
mole.h
diatomics::states
qList states
Definition: h2_priv.h:565
diatomics::H2_PunchLineStuff
void H2_PunchLineStuff(FILE *io, realnum xLimit, long index)
Definition: mole_h2_io.cpp:1167
RandGauss
double RandGauss(double xMean, double s)
Definition: service.cpp:1643
spsort
void spsort(realnum x[], long int n, long int iperm[], int kflag, int *ier)
Definition: service.cpp:1100
diatomics::H2_ReadTransprob
void H2_ReadTransprob(long int nelec, TransitionList &trans)
Definition: mole_h2_io.cpp:430
t_hmi::H2_Solomon_dissoc_rate_BD96_H2g
double H2_Solomon_dissoc_rate_BD96_H2g
Definition: hmi.h:95
level_tmp::n
long n
Definition: mole_h2_io.cpp:667
level_tmp::J
long J
Definition: mole_h2_io.cpp:667
phycon
t_phycon phycon
Definition: phycon.cpp:6
t_hmi::UV_Cont_rel2_Habing_spec_depth
realnum UV_Cont_rel2_Habing_spec_depth
Definition: hmi.h:66
diatomics::H2_col_rate_in
multi_arr< double, 2 > H2_col_rate_in
Definition: h2_priv.h:651
t_LineSave::ipNormWavL
long int ipNormWavL
Definition: lines.h:81
TransitionList::begin
iterator begin(void)
Definition: transition.h:305
diatomics::trans
TransitionList trans
Definition: h2_priv.h:566
diatomics::H2_Read_hminus_distribution
void H2_Read_hminus_distribution(void)
Definition: mole_h2_io.cpp:984
ProxyIterator
Definition: proxy_iterator.h:58
SDIV
sys_float SDIV(sys_float x)
Definition: cddefines.h:952
diatomics::lgREAD_DATA
bool lgREAD_DATA
Definition: h2_priv.h:252
diatomics::H2_Prt_Zone
void H2_Prt_Zone(void)
Definition: mole_h2_io.cpp:333
diatomics::H2_rad_rate_in
multi_arr< double, 2 > H2_rad_rate_in
Definition: h2_priv.h:653
t_dense::gas_phase
realnum gas_phase[LIMELM]
Definition: dense.h:71
lines_service.h
EmissionProxy::Pdest
realnum & Pdest() const
Definition: emission.h:543
t_LineSave::nsum
long int nsum
Definition: lines.h:62
diatomics::pops_per_elec
double pops_per_elec[N_ELEC]
Definition: h2_priv.h:620
diatomics::average_energy_s
double average_energy_s
Definition: h2_priv.h:287
diatomics::ortho_density
double ortho_density
Definition: h2_priv.h:319
t_hmi::H2_Solomon_dissoc_rate_ELWERT_H2g
double H2_Solomon_dissoc_rate_ELWERT_H2g
Definition: hmi.h:96
input
t_input input
Definition: input.cpp:12
ASSERT
#define ASSERT(exp)
Definition: cddefines.h:578
diatomics::H2_X_colden
multi_arr< realnum, 2 > H2_X_colden
Definition: h2_priv.h:660
t_rfield::extin_mag_V_point
double extin_mag_V_point
Definition: rfield.h:277
qList::end
iterator end()
Definition: quantumstate.h:345
diatomics::rad_end
TransitionList::iterator rad_end
Definition: h2_priv.h:567
t_save::whichDiatomToPrint
diatomics * whichDiatomToPrint[LIMPUN]
Definition: save.h:226
ipHYDROGEN
const int ipHYDROGEN
Definition: cddefines.h:305
t_thermal::ctot
double ctot
Definition: thermal.h:112
mr3ci
multi_arr< realnum, 3 >::const_iterator mr3ci
Definition: container_classes.h:1825
MIN2
#define MIN2
Definition: cddefines.h:761
lindst
void lindst(double xInten, realnum wavelength, const char *chLab, long int ipnt, char chInfo, bool lgOutToo, const char *chComment)
Definition: lines_service.cpp:468
hyperfine
t_hyperfine hyperfine
Definition: hyperfine.cpp:5
diatomics::H2_ReadEnergies
void H2_ReadEnergies()
Definition: mole_h2_io.cpp:672
cddrive.h
diatomics::nLevels_per_elec
long int nLevels_per_elec[N_ELEC]
Definition: h2_priv.h:618
LineSave
t_LineSave LineSave
Definition: lines.cpp:5
diatomics::getLine
long int getLine(long iElecHi, long iVibHi, long iRotHi, long iElecLo, long iVibLo, long iRotLo, double *relint, double *absint)
Definition: mole_h2_io.cpp:1931
radius
t_radius radius
Definition: radius.cpp:5
diatomics::H2_populations_LTE
multi_arr< double, 3 > H2_populations_LTE
Definition: h2_priv.h:639
diatomics::path
string path
Definition: h2_priv.h:573
diatomics::H2_renorm_chemistry
double H2_renorm_chemistry
Definition: h2_priv.h:603
TransitionList::end
iterator end(void)
Definition: transition.h:309
diatomics::H2_den_g
double H2_den_g
Definition: h2_priv.h:682
t_hmi::UV_Cont_rel2_Habing_TH85_depth
realnum UV_Cont_rel2_Habing_TH85_depth
Definition: hmi.h:64
level_tmp::eWN
double eWN
Definition: mole_h2_io.cpp:668
EXIT_FAILURE
#define EXIT_FAILURE
Definition: cddefines.h:140
diatomics::nCall_this_zone
long int nCall_this_zone
Definition: h2_priv.h:341
diatomics::dense_total
const double *const dense_total
Definition: h2_priv.h:589
Parser
Definition: parser.h:31
dense.h
t_input::chDelimiter
char chDelimiter[3]
Definition: input.h:42
diatomics::ipTransitionSort
multi_arr< long int, 2 > ipTransitionSort
Definition: h2_priv.h:691
diatomics::H2_ReadDissprob
void H2_ReadDissprob(long int nelec)
Definition: mole_h2_io.cpp:893
diatomics::nRot_hi
valarray< long > nRot_hi[N_ELEC]
Definition: h2_priv.h:613
cddefines.h
qList::resize
void resize(size_t i)
Definition: quantumstate.h:83
diatomics::H2_X_colden_LTE
multi_arr< realnum, 2 > H2_X_colden_LTE
Definition: h2_priv.h:662
diatomics::Solomon_elec_decay_s
double Solomon_elec_decay_s
Definition: h2_priv.h:269
diatomics::Solomon_dissoc_rate_g
double Solomon_dissoc_rate_g
Definition: h2_priv.h:264
diatomics::n_elec_states
long int n_elec_states
Definition: h2_priv.h:409
t_LineSave::ScaleNormLine
double ScaleNormLine
Definition: lines.h:94
t_LineSave::ipass
long int ipass
Definition: lines.h:75
thermal
t_thermal thermal
Definition: thermal.cpp:5
diatomics::H2_ReadDissocEnergies
void H2_ReadDissocEnergies(void)
Definition: mole_h2_io.cpp:829
diatomics::H2_PrtDepartCoef
void H2_PrtDepartCoef(void)
Definition: mole_h2_io.cpp:362
thresh_punline_h2
static realnum thresh_punline_h2
Definition: mole_h2_io.cpp:46
t_hmi::HeatH2Dish_TH85
double HeatH2Dish_TH85
Definition: hmi.h:130
diatomics::para_density
double para_density
Definition: h2_priv.h:321
EmissionProxy
Definition: emission.h:14
diatomics::H2_col_rate_out
multi_arr< double, 2 > H2_col_rate_out
Definition: h2_priv.h:652
diatomics::Solomon_elec_decay_g
double Solomon_elec_decay_g
Definition: h2_priv.h:268
t_tag_LineSv::SumLine
double SumLine[4]
Definition: lines.h:125
diatomics::HeatDexc
double HeatDexc
Definition: h2_priv.h:290
hyperfine.h
level_tmp::operator<
bool operator<(const level_tmp &second) const
Definition: mole_h2_io.cpp:660
Parser::getNumberDefaultNegImplLog
double getNumberDefaultNegImplLog(const char *chDesc, double fdef)
Definition: parser.cpp:336
radius.h
diatomics::H2_lgOrtho
multi_arr< bool, 3 > H2_lgOrtho
Definition: h2_priv.h:643
colden
t_colden colden
Definition: colden.cpp:5
diatomics::ortho_colden
double ortho_colden
Definition: h2_priv.h:328
diatomics::chH2ColliderLabels
char chH2ColliderLabels[N_X_COLLIDER][chN_X_COLLIDER]
Definition: h2_priv.h:591
hmi.h
diatomics::nElecLevelOutput
int nElecLevelOutput
Definition: h2_priv.h:349
t_colden::colden
realnum colden[NCOLD]
Definition: colden.h:38
diatomics::H2_dissprob
multi_arr< realnum, 3 > H2_dissprob
Definition: h2_priv.h:632
diatomics::H2_DissocEnergies
double H2_DissocEnergies[N_ELEC]
Definition: h2_priv.h:609
MAX2
#define MAX2
Definition: cddefines.h:782
Save1Line
void Save1Line(const TransitionProxy &t, FILE *io, realnum xLimit, long index, realnum DopplerWidth)
Definition: save_do.cpp:4347
multi_arr::begin
iterator begin(size_type i1)
Definition: container_classes.h:1519
level_tmp::v
long v
Definition: mole_h2_io.cpp:667
diatomics::CollRateCoeff
multi_arr< realnum, 3 > CollRateCoeff
Definition: h2_priv.h:621
cdEXIT
#define cdEXIT(FAIL)
Definition: cddefines.h:434
t_save::nsave
long int nsave
Definition: save.h:222
N_ELEC
const int N_ELEC
Definition: h2_priv.h:27
save.h
diatomics::H2_Punch_line_data
void H2_Punch_line_data(FILE *ioPUN, bool lgDoAll)
Definition: mole_h2_io.cpp:1125
chN_X_COLLIDER
const int chN_X_COLLIDER
Definition: h2_priv.h:15
t_save::chSave
char chSave[LIMPUN][5]
Definition: save.h:225
diatomics::H2_Prt_line_tau
void H2_Prt_line_tau(void)
Definition: mole_h2_io.cpp:1187
prt.h
qList::size
size_t size() const
Definition: quantumstate.h:116
EmissionProxy::Pelec_esc
realnum & Pelec_esc() const
Definition: emission.h:533
doppvel.h
t_thermal::htot
double htot
Definition: thermal.h:149
t_hmi::UV_Cont_rel2_Draine_DB96_depth
realnum UV_Cont_rel2_Draine_DB96_depth
Definition: hmi.h:74
FILENAME_PATH_LENGTH_2
const int FILENAME_PATH_LENGTH_2
Definition: cddefines.h:249
grainvar.h
GrainVar::rate_h2_form_grains_used_total
double rate_h2_form_grains_used_total
Definition: grainvar.h:574
t_secondaries::csupra
realnum ** csupra
Definition: secondaries.h:21
PutLine
void PutLine(const TransitionProxy &t, const char *chComment, const char *chLabelTemp)
Definition: transition.cpp:449
h2_priv.h
t_dense::AtomicWeight
realnum AtomicWeight[LIMELM]
Definition: dense.h:75
diatomics::para_colden
double para_colden
Definition: h2_priv.h:329
TransitionList
Definition: transition.h:274
diatomics::H2_Prt_column_density
void H2_Prt_column_density(FILE *ioMEAN)
Definition: mole_h2_io.cpp:395
parser.h
Parser::getNumberDefault
double getNumberDefault(const char *chDesc, double fdef)
Definition: parser.cpp:282
EmissionProxy::Pesc
realnum & Pesc() const
Definition: emission.h:523
diatomics::ipEnergySort
multi_arr< long int, 3 > ipEnergySort
Definition: h2_priv.h:690
nTE_HMINUS
const int nTE_HMINUS
Definition: h2_priv.h:24
hmi
t_hmi hmi
Definition: hmi.cpp:5
t_hmi::HeatH2Dexc_TH85
double HeatH2Dexc_TH85
Definition: hmi.h:138
physconst.h
secondaries
t_secondaries secondaries
Definition: secondaries.cpp:5
GetDopplerWidth
realnum GetDopplerWidth(realnum massAMU)
Definition: temp_change.cpp:499
ipCOL_H2g
#define ipCOL_H2g
Definition: colden.h:16
diatomics::HeatDiss
double HeatDiss
Definition: h2_priv.h:289
gv
GrainVar gv
Definition: grainvar.cpp:5
cdH2_Line
long int cdH2_Line(long int iElecHi, long int iVibHi, long int iRotHi, long int iElecLo, long int iVibLo, long int iRotLo, double *relint, double *absint)
Definition: mole_h2_io.cpp:1913
diatomics::ipRot_H2_energy_sort
valarray< long > ipRot_H2_energy_sort
Definition: h2_priv.h:689
PrintEfmt
#define PrintEfmt(F, V)
Definition: cddefines.h:1472
prt_wl
void prt_wl(FILE *ioOUT, realnum wl)
Definition: prt.cpp:13
read_whole_line
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
Definition: service.cpp:70
taulines.h
diatomics::Solomon_dissoc_rate_s
double Solomon_dissoc_rate_s
Definition: h2_priv.h:265
t_hmi::H2_Solomon_dissoc_rate_TH85_H2g
double H2_Solomon_dissoc_rate_TH85_H2g
Definition: hmi.h:93
t_colden::coldenH2_ov_vel
realnum coldenH2_ov_vel
Definition: colden.h:43
t_radius::depth_mid_zone
double depth_mid_zone
Definition: radius.h:41
phycon.h
ShowMe
void ShowMe(void)
Definition: service.cpp:181
t_save::punarg
realnum punarg[LIMPUN][3]
Definition: save.h:254
chMolBranch
STATIC char chMolBranch(long iRotHi, long int iRotLo)
Definition: mole_h2_io.cpp:1207
N_X_COLLIDER
const int N_X_COLLIDER
Definition: h2_priv.h:13
diatomics::H2_stat
multi_arr< realnum, 3 > H2_stat
Definition: h2_priv.h:641
diatomics::H2_ParseSave
void H2_ParseSave(Parser &p, char *chHeader)
Definition: mole_h2_io.cpp:111
diatomics::ipVib_H2_energy_sort
valarray< long > ipVib_H2_energy_sort
Definition: h2_priv.h:687
diatomics::H2_old_populations
multi_arr< double, 3 > H2_old_populations
Definition: h2_priv.h:637
h2.h
LineSv
LinSv * LineSv
Definition: cdinit.cpp:70
diatomics::lgH2_radiative
multi_arr< bool, 2 > lgH2_radiative
Definition: h2_priv.h:714
diatomics::nXLevelsMatrix
long int nXLevelsMatrix
Definition: h2_priv.h:695
EmissionProxy::Aul
realnum & Aul() const
Definition: emission.h:613
t_phycon::te
double te
Definition: phycon.h:11
diatomics::Jlowest
long int Jlowest[N_ELEC]
Definition: h2_priv.h:616
input.h
t_hmi::H2_rate_destroy
double H2_rate_destroy
Definition: hmi.h:21
diatomics::H2_rad_rate_out
multi_arr< double, 3 > H2_rad_rate_out
Definition: h2_priv.h:634
diatomics::H2_X_hminus_formation_distribution
multi_arr< realnum, 3 > H2_X_hminus_formation_distribution
Definition: h2_priv.h:685
diatomics::set_numLevelsMatrix
void set_numLevelsMatrix(long numLevels)
Definition: mole_h2_io.cpp:1994
DEBUG_ENTRY
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:684
qList::begin
iterator begin()
Definition: quantumstate.h:337
save
t_save save
Definition: save.cpp:5
t_hmi::H2_Solomon_dissoc_rate_BD96_H2s
double H2_Solomon_dissoc_rate_BD96_H2s
Definition: hmi.h:101
SMALLFLOAT
const realnum SMALLFLOAT
Definition: cpu.h:191
diatomics::H2_den_s
double H2_den_s
Definition: h2_priv.h:682
t_hmi::H2_Solomon_dissoc_rate_TH85_H2s
double H2_Solomon_dissoc_rate_TH85_H2s
Definition: hmi.h:99
TransitionList::Emis
EmissionList & Emis()
Definition: transition.h:329
g
static double * g
Definition: species2.cpp:28
Generated by   doxygen 1.8.17