cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
mole_h2_io.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 /*H2_ParsePunch parse the punch h2 command */
4 /*H2_PunchDo punch 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 PunchLineStuff */
8 /*H2_Punch_line_data punch 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_ParsePunch parse the punch 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 "punch.h"
23 #include "hmi.h"
24 #include "prt.h"
25 #include "secondaries.h"
26 #include "grainvar.h"
27 #include "phycon.h"
28 #include "rfield.h"
29 #include "hyperfine.h"
30 #include "thermal.h"
31 #include "lines.h"
32 #include "lines_service.h"
33 #include "dense.h"
34 #include "radius.h"
35 #include "colden.h"
36 #include "taulines.h"
37 #include "h2.h"
38 #include "h2_priv.h"
39 #include "cddrive.h"
40 #include "mole.h"
41 
42 /* this will say whether ortho or para,
43  * H2_lgOrtho is 0 or 1 depending on whether or not ortho,
44  * so chlgPara[H2_lgOrtho] gives P or O for printing */
45 static char chlgPara[2]={'P','O'};
46 
47 /* intensity, relative to normalization line, for faintest line to punch */
49 
50 /*H2_LinesAdd add in explicit lines from the large H2 molecule, called by lines_molecules */
51 void H2_LinesAdd(void)
52 {
53  /* these are the quantum designations of the lines we will output */
54  int iRotHi, iVibHi, iElecHi ,iRotLo, iVibLo, iElecLo;
55 
56  /* H2 not on, so space not allocated */
57  if( !h2.lgH2ON )
58  return;
59 
60  DEBUG_ENTRY( "H2_LinesAdd()" );
61 
62  /* >>chng 05 nov 04, make info copies of these lines up here
63  * these are among the strongest lines in the 2 micron window and some have nearly the same
64  * wavelength as far weaker lines that may come before them in the line stack. in that case
65  * cdLine would find the much weaker line with the same wavelength.
66  * put strong H2 lines first so that line search will find these, and not far weaker
67  * lines with nearly the same wavelength - these will be duplicated in the output but
68  * these are here for into (the 'i) so does no harm
69  *
70  * the array indices in the following structures give upper and lower electronic, vib, rot quantum
71  * indices.
72  * H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo]
73  * >>chng 05 dec 22, had hand entered wavelength in A as second parameter. This gave
74  * rounded off result when set line precision 5 was used. now uses same logic that
75  * PutLine will eventually use - simply enter same wl in Ang
76  * 1-0 S(4) - 18910 */
77  lindst( H2Lines[0][1][6][0][0][4].Emis->xIntensity , H2Lines[0][1][6][0][0][4].WLAng , "H2 " ,
78  H2Lines[0][1][6][0][0][4].ipCont,'i',false ,
79  "H2 line");
80  /* 1-0 S(3) - 19570 */
81  lindst( H2Lines[0][1][5][0][0][3].Emis->xIntensity , H2Lines[0][1][5][0][0][3].WLAng , "H2 " ,
82  H2Lines[0][1][5][0][0][3].ipCont,'i',false ,
83  "H2 line");
84  /* 1-0 S(2) - 20330 */
85  lindst( H2Lines[0][1][4][0][0][2].Emis->xIntensity , H2Lines[0][1][4][0][0][2].WLAng , "H2 " ,
86  H2Lines[0][1][4][0][0][2].ipCont,'i',false ,
87  "H2 line");
88  /* 1-0 S(1) - 21210 */
89  lindst( H2Lines[0][1][3][0][0][1].Emis->xIntensity , H2Lines[0][1][3][0][0][1].WLAng , "H2 " ,
90  H2Lines[0][1][3][0][0][1].ipCont,'i',false ,
91  "H2 line");
92  /* 1-0 S(0) - 22230 */
93  lindst( H2Lines[0][1][2][0][0][0].Emis->xIntensity , H2Lines[0][1][2][0][0][0].WLAng , "H2 " ,
94  H2Lines[0][1][2][0][0][0].ipCont,'i',false ,
95  "H2 line");
96  /* start Q branch - selection rule requires that J be non-zero, so no Q(0) */
97  /* 1-0 Q(2) - 24130 */
98  lindst( H2Lines[0][1][2][0][0][2].Emis->xIntensity , H2Lines[0][1][2][0][0][2].WLAng , "H2 " ,
99  H2Lines[0][1][2][0][0][2].ipCont,'i',false ,
100  "H2 line");
101  /* 1-0 Q(1) - 24060 */
102  lindst( H2Lines[0][1][1][0][0][1].Emis->xIntensity , H2Lines[0][1][1][0][0][1].WLAng , "H2 " ,
103  H2Lines[0][1][1][0][0][1].ipCont,'i',false ,
104  "H2 line");
105 
106  /* print all lines from lowest n levels within X */
107  /* loop over all possible lines and set H2_populations,
108  * and quantities that depend on them */
109  for( iElecHi=0; iElecHi<h2.nElecLevelOutput; ++iElecHi )
110  {
111  for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
112  {
113  for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
114  {
115  long int lim_elec_lo = 0;
116  /* now the lower levels */
117  /* NB - X is the only lower level considered here, since we are only
118  * concerned with excited electronic levels as a photodissociation process
119  * code exists to relax this assumption - simply change following to iElecHi */
120  for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
121  {
122  /* want to include all vib states in lower level if different electronic level,
123  * but only lower vib levels if same electronic level */
124  long int nv = h2.nVib_hi[iElecLo];
125  if( iElecLo==iElecHi )
126  nv = iVibHi;
127  for( iVibLo=0; iVibLo<=nv; ++iVibLo )
128  {
129  long nr = h2.nRot_hi[iElecLo][iVibLo];
130  if( iElecLo==iElecHi && iVibHi==iVibLo )
131  nr = iRotHi-1;
132 
133  for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
134  {
135  if( lgH2_line_exists[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] )
136  {
137  /* all ground vib state rotation lines - first is J to J-2 */
138  PutLine(&H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo],
139  "H2 lines");
140  if( LineSave.ipass == 0 )
141  {
142  H2_SaveLine[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] = 0.;
143  }
144  else if( LineSave.ipass == 1 )
145  {
146  H2_SaveLine[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] += (realnum)(
147  radius.dVeff*H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->xIntensity);
148  }
149  }
150  }
151  }
152  }
153  }
154  }
155  }
156  return;
157 }
158 
159 /*H2_ParsePunch parse the punch h2 command */
160 void H2_ParsePunch( char *chCard ,
161  char *chHeader)
162 {
163  long int i , nelem;
164  bool lgEOL;
165 
166  DEBUG_ENTRY( "H2_ParsePunch()" );
167 
168  i = 5;
169  nelem = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
170  if( nelem != 2 )
171  {
172  fprintf( ioQQQ, " The first number on this line must be the 2 in H2\n Sorry.\n" );
173  cdEXIT(EXIT_FAILURE);
174  }
175  /* this provides info on the large H2 molecule */
176  if( nMatch("COLU",chCard) )
177  {
178  /* punch column density */
179  strcpy( punch.chPunch[punch.npunch], "H2cl" );
180 
181  /* this is an option to scan off highest vib and rot states
182  * to punch pops - first is limit to vibration, then rotation
183  * if no number is entered then 0 is set and all levels punched */
184  /* now get vib limit */
185  punch.punarg[punch.npunch][0] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
186 
187  /* highest rotation */
188  punch.punarg[punch.npunch][1] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
189  /* this says whether to punch triplets or a matrix for output -
190  * default is triplets, so only check for matrix */
191  if( nMatch( "MATR" , chCard ) )
192  {
193  /* matrix */
194  punch.punarg[punch.npunch][2] = 1;
195  sprintf( chHeader, "#vib\trot\tcolumn density\n" );
196  }
197  else
198  {
199  /* triplets */
200  punch.punarg[punch.npunch][2] = -1;
201  sprintf( chHeader, "#vib\trot\tEner(K)\tcolden\tcolden/stat wght\tLTE colden\tLTE colden/stat wght\n" );
202  }
203  }
204  else if( nMatch("COOL",chCard) )
205  {
206  /* heating and cooling rates */
207  strcpy( punch.chPunch[punch.npunch], "H2co" );
208  sprintf( chHeader,
209  "#H2 depth\ttot cool\tTH Sol\tBig Sol\tTH pht dis\tpht dis\tTH Xcool\tXcool \n" );
210  }
211 
212  else if( nMatch("CREA",chCard) )
213  {
214  /* H2 creation rates */
215  strcpy( punch.chPunch[punch.npunch], "H2cr" );
216  sprintf( chHeader,
217  "#H2 depth\tH2_rate_create\tH2_rate_destroy\trate_h2_form_grains_used_total\tassoc_detach"
218  "\tbh2dis\tbh2h2p\tradasc\th3ph2p\th2phmh2h\tbh2h22hh2\th3phmh2hh\th3phm2h2\th32h2\teh3_h2h\th3ph2hp"
219  "\tH_CH_C_H2\tH_CHP_CP_H2\tH_CH2_CH_H2\tH_CH3P_CH2P_H2\tH_OH_O_H2\tHminus_HCOP_CO_H2\tHminus_H3OP_H2O_H2\tHminus_H3OP_OH_H2_H"
220  "\tHP_CH2_CHP_H2\tHP_SiH_SiP_H2\tH2P_CH_CHP_H2\tH2P_CH2_CH2P_H2\tH2P_CO_COP_H2\tH2P_H2O_H2OP_H2\tH2P_O2_O2P_H2"
221  "\tH2P_OH_OHP_H2\tH3P_C_CHP_H2\tH3P_CH_CH2P_H2\tH3P_CH2_CH3P_H2\tH3P_OH_H2OP_H2\tH3P_H2O_H3OP_H2\tH3P_CO_HCOP_H2"
222  "\tH3P_O_OHP_H2\tH3P_SiH_SiH2P_H2\tH3P_SiO_SiOHP_H2\tH_CH3_CH2_H2\tH_CH4P_CH3P_H2\tH_CH5P_CH4P_H2\tH2P_CH4_CH3P_H2"
223  "\tH2P_CH4_CH4P_H2\tH3P_CH3_CH4P_H2\tH3P_CH4_CH5P_H2\tHP_CH4_CH3P_H2\tHP_HNO_NOP_H2\tHP_HS_SP_H2\tH_HSP_SP_H2"
224  "\tH3P_NH_NH2P_H2\tH3P_NH2_NH3P_H2\tH3P_NH3_NH4P_H2\tH3P_CN_HCNP_H2\tH3P_NO_HNOP_H2\tH3P_S_HSP_H2\tH3P_CS_HCSP_H2"
225  "\tH3P_NO2_NOP_OH_H2\tH2P_NH_NHP_H2\tH2P_NH2_NH2P_H2\tH2P_NH3_NH3P_H2\tH2P_CN_CNP_H2\tH2P_HCN_HCNP_H2\tH2P_NO_NOP_H2"
226  "\tH3P_Cl_HClP_H2\tH3P_HCl_H2ClP_H2\tH2P_C2_C2P_H2\tHminus_NH4P_NH3_H2\tH3P_HCN_HCNHP_H2"
227  "\tdestruction/creation\tav Einstein A\n");
228  }
229  else if( nMatch("DEST",chCard) )
230  {
231  /* punch H2 destruction - output destruction rates */
232  strcpy( punch.chPunch[punch.npunch], "H2ds" );
233  sprintf( chHeader,
234  "#depth\ttot H2 rate create\ttot H2 rate destroy\ttot H- backwards\tSolomon H2g\tSolomon H2s\tphotodissoc H2s\tphotodissoc H2g"
235  "\te- dissoc\trh2h2p\th2hph3p\tH0 dissoc\tCR\trheph2hpheh\theph2heh2p\thehph2h3phe\th3petc\tH2Ph3p"
236  "\th2sh2g\th2h22hh2\th2sh2sh2g2h\th2sh2sh2s2h\tH2_CHP_CH2P_H\tH2_CH2P_CH3P_H\tH2_OHP_H2OP_H\tH2_H2OP_H3OP_H\tH2_COP_HCOP_H"
237  "\tH2_OP_OHP_H\tH2_SiOP_SiOHP_H\tH2_C_CH_H\tH2_CP_CHP_H\tH2_CH_CH2_H\tH2_OH_H2O_H\tH2_O_OH_H"
238  "\th2s_ch_ch2_h\th2s_o_oh_h\th2s_oh_h2o_h\th2s_c_ch_h\th2s_cp_chp_h\tH2_CH2_CH3_H\tH2_CH3_CH4_H"
239  "\tH2_CH4P_CH5P_H\tH2s_CH2_CH3_H\tH2s_CH3_CH4_H\th2s_op_ohp_h\tH2_N_NH_H\tH2_NH_NH2_H\tH2_NH2_NH3_H\tH2_CN_HCN_H\tH2_NP_NHP_H"
240  "\tH2_NHP_N_H3P\tH2_NHP_NH2P_H\tH2_NH2P_NH3P_H\tH2_NH3P_NH4P_H\tH2_CNP_HCNP_H\tH2_SP_HSP_H\tH2_CSP_HCSP_H"
241  "\tH2_ClP_HClP_H\tH2_HClP_H2ClP_H\tH2_HCNP_HCNHP_H"
242  "\tfrac H2g\tfrac H2s\n");
243  }
244 
245  else if( nMatch("HEAT",chCard) )
246  {
247  /* heating and cooling rates */
248  strcpy( punch.chPunch[punch.npunch], "H2he" );
249  sprintf( chHeader,
250  "#H2 depth\ttot Heat\tHeat(big)\tHeat(TH85)\tDissoc(Big)\tDissoc(TH85) \n" );
251  }
252 
253  else if( nMatch("LEVE",chCard) )
254  {
255  /* punch H2 level energies */
256  strcpy( punch.chPunch[punch.npunch], "H2le" );
257  sprintf( chHeader,
258  "#H2 v\tJ\tenergy(wn)\tstat wght\tSum As" );
259  char chHoldit[chN_X_COLLIDER+12];
260  for( int nColl=0; nColl<N_X_COLLIDER; ++nColl )
261  {
262  /* labels for all colliders */
263  sprintf(chHoldit,"\tCritDen %s",chH2ColliderLabels[nColl]);
264  strcat( chHeader , chHoldit );
265  }
266  strcat( chHeader , "\n" );
267  }
268 
269  else if( nMatch("LINE",chCard) )
270  {
271  /* punch H2 lines - all in X */
272  strcpy( punch.chPunch[punch.npunch], "H2ln" );
273  sprintf( chHeader,
274  "#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" );
275  /* first optional number changes the threshold of weakest line to print*/
276  /* fe2thresh is intensity relative to normalization line,
277  * normally Hbeta, and is set to zero in zero.c */
278 
279  /* threshold for faintest line to punch, default is 1e-4 of norm line */
281  if( lgEOL )
282  {
283  /* this will be default relative intensity for faintest line to punch */
284  thresh_punline_h2 = 1e-4f;
285  }
286 
287  /* it is a log if negative */
288  if( thresh_punline_h2 < 0. )
289  {
291  }
292 
293  /* lines from how many electronic states? default is one, just X, and is
294  * obtained with GROUND keyword. ALL will produce all lines from all levels.
295  * else, if a number is present, will be the number. if no number, no keyword,
296  * appear then just ground */
297  if( nMatch( "ELEC" , chCard ) )
298  {
299  if( nMatch(" ALL",chCard) )
300  {
301  /* all electronic levels - when done, will set upper limit, the
302  * number of electronic levels actually computed, don't know this yet,
303  * so signify with negative number */
304  h2.nElecLevelOutput = -1;
305  }
306  else if( nMatch("GROU",chCard) )
307  {
308  /* just the ground electronic state */
309  h2.nElecLevelOutput = 1;
310  }
311  else
312  {
313  h2.nElecLevelOutput = (int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
314  if( lgEOL )
315  {
316  /* just the ground electronic state */
317  h2.nElecLevelOutput = 1;
318  }
319  }
320  }
321  }
322 
323  else if( nMatch(" PDR",chCard) )
324  {
325  /* creation and destruction processes */
326  strcpy( punch.chPunch[punch.npunch], "H2pd" );
327  sprintf( chHeader, "#H2 creation, destruction. \n" );
328  }
329  else if( nMatch("POPU",chCard) )
330  {
331  /* punch H2_populations */
332  strcpy( punch.chPunch[punch.npunch], "H2po" );
333 
334  /* this is an option to scan off highest vib and rot states
335  * to punch pops - first is limit to vibration, then rotation
336  * if no number is entered then 0 is set and all levels punched */
337  /* now get vib lim */
338  punch.punarg[punch.npunch][0] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
339 
340  /* this is limit to rotation quantum index */
341  punch.punarg[punch.npunch][1] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
342 
343  if( nMatch( "ZONE" , chCard ) )
344  {
345  /* punch v=0 pops for each zone, all along one line */
346  punch.punarg[punch.npunch][2] = 0;
347  sprintf( chHeader, "#depth\torth\tpar\te=1 rel pop\te=2 rel pop\tv=0 rel pops\n" );
348  }
349  else
350  {
351  /* will not do zone output, only output at the end of the calculation
352  * now check whether to punch triplets or a matrix for output -
353  * default is triplets, so only check for matrix */
354  if( nMatch( "MATR" , chCard ) )
355  {
356  /* matrix */
357  punch.punarg[punch.npunch][2] = 1;
358  sprintf( chHeader, "#vib\trot\tpops\n" );
359  }
360  else
361  {
362  /* triplets */
363  punch.punarg[punch.npunch][2] = -1;
364  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" );
365  }
366  }
367  }
368 
369  else if( nMatch("RATE",chCard) )
370  {
371  /* punch h2 rates - creation and destruction rates */
372  strcpy( punch.chPunch[punch.npunch], "H2ra" );
373  sprintf( chHeader,
374  "#depth\tN(H2)\tN(H2)/u(H2)\tA_V(star)\tn(Eval)"
375  "\tH2/Htot\trenorm\tfrm grn\tfrmH-\tdstTH85\tBD96\tELWERT\tBigH2\telec->H2g\telec->H2s"
376  "\tG(TH85)\tG(DB96)\tCR\tEleclife\tShield(BD96)\tShield(H2)\tBigh2/G0(spc)\ttot dest"
377  "\tHeatH2Dish_TH85\tHeatH2Dexc_TH85\tHeatH2Dish_BigH2\tHeatH2Dexc_BigH2\thtot\n" );
378  }
379  else if( nMatch("SOLO",chCard) )
380  {
381  /* rate of Solomon process then fracs of exits from each v, J level */
382  strcpy( punch.chPunch[punch.npunch], "H2so" );
383  sprintf( chHeader,
384  "#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" );
385  }
386  else if( nMatch("SPEC",chCard) )
387  {
388  /* special punch command*/
389  strcpy( punch.chPunch[punch.npunch], "H2sp" );
390  sprintf( chHeader,
391  "#depth\tspecial\n" );
392  }
393  else if( nMatch("TEMP",chCard) )
394  {
395  /* various temperatures for neutral/molecular gas */
396  strcpy( punch.chPunch[punch.npunch], "H2te" );
397  sprintf( chHeader,
398  "#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");
399  }
400  else if( nMatch("THER",chCard) )
401  {
402  /* thermal heating cooling processes involving H2 */
403  strcpy( punch.chPunch[punch.npunch], "H2th" );
404  sprintf( chHeader,
405  "#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");
406  }
407  else
408  {
409  fprintf( ioQQQ,
410  " There must be a second key; they are RATE, LINE, COOL, COLUMN, _PDR, SOLOmon, TEMP, and POPUlations\n" );
411  cdEXIT(EXIT_FAILURE);
412  }
413  return;
414 }
415 
416 
417 /*H2_Prt_Zone print H2 info into zone results, called from prtzone for each printed zone */
418 void H2_Prt_Zone(void)
419 {
420  int iElecHi , iVibHi;
421 
422  /* no print if H2 not turned on, or not computed for these conditions */
423  if( !h2.lgH2ON || !h2.nCallH2_this_zone )
424  return;
425 
426  DEBUG_ENTRY( "H2_Prt_Zone()" );
427 
428  fprintf( ioQQQ, " H2 density ");
429  fprintf(ioQQQ,PrintEfmt("%9.2e", hmi.H2_total));
430 
431  fprintf( ioQQQ, " orth/par");
432  fprintf(ioQQQ,PrintEfmt("%9.2e", h2.ortho_density / SDIV( h2.para_density )));
433 
434  iElecHi = 0;
435  iVibHi = 0;
436  fprintf( ioQQQ, " v0 J=0,3");
437  fprintf(ioQQQ,PrintEfmt("%9.2e", H2_populations[iElecHi][iVibHi][0] / hmi.H2_total));
438  fprintf(ioQQQ,PrintEfmt("%9.2e", H2_populations[iElecHi][iVibHi][1] / hmi.H2_total));
439  fprintf(ioQQQ,PrintEfmt("%9.2e", H2_populations[iElecHi][iVibHi][2] / hmi.H2_total));
440  fprintf(ioQQQ,PrintEfmt("%9.2e", H2_populations[iElecHi][iVibHi][3] / hmi.H2_total));
441 
442  fprintf( ioQQQ, " TOTv=0,3");
443  fprintf(ioQQQ,PrintEfmt("%9.2e", pops_per_vib[iElecHi][0] / hmi.H2_total));
444  fprintf(ioQQQ,PrintEfmt("%9.2e", pops_per_vib[iElecHi][1] / hmi.H2_total));
445  fprintf(ioQQQ,PrintEfmt("%9.2e", pops_per_vib[iElecHi][2] / hmi.H2_total));
446  fprintf(ioQQQ,PrintEfmt("%9.2e", pops_per_vib[iElecHi][3] / hmi.H2_total));
447  fprintf( ioQQQ, "\n");
448  return;
449 }
450 
451 /*H2_Prt_column_density print H2 info into zone results, called from prtzone for each printed zone */
453  /* this is stream used for io, is stdout when called by final,
454  * is punch unit when punch output generated */
455  FILE *ioMEAN )
456 
457 {
458  int iVibHi;
459 
460  /* no print if H2 not turned on, or not computed for these conditions */
461  if( !h2.lgH2ON || !h2.nCallH2_this_zone )
462  return;
463 
464  DEBUG_ENTRY( "H2_Prt_column_density()" );
465 
466  fprintf( ioMEAN, " H2 total ");
467  fprintf(ioMEAN,"%7.3f", log10(SDIV(colden.colden[ipCOL_H2g]+colden.colden[ipCOL_H2s])));
468 
469  fprintf( ioMEAN, " H2 ortho ");
470  fprintf(ioMEAN,"%7.3f", log10(SDIV(h2.ortho_colden)));
471 
472  fprintf( ioMEAN, " para");
473  fprintf(ioMEAN,"%7.3f", log10(SDIV(h2.para_colden)));
474 
475  iVibHi = 0;
476  fprintf( ioMEAN, " v0 J=0,3");
477  fprintf(ioMEAN,"%7.3f", log10(SDIV(H2_X_colden[iVibHi][0])));
478  fprintf(ioMEAN,"%7.3f", log10(SDIV(H2_X_colden[iVibHi][1])));
479  fprintf(ioMEAN,"%7.3f", log10(SDIV(H2_X_colden[iVibHi][2])));
480  fprintf(ioMEAN,"%7.3f", log10(SDIV(H2_X_colden[iVibHi][3])));
481 
482 # if 0
483  fprintf( ioMEAN, " v=0,3");
484  fprintf(ioMEAN,PrintEfmt("%9.2e", pops_per_vib[iElecHi][0] / hmi.H2_total));
485  fprintf(ioMEAN,PrintEfmt("%9.2e", pops_per_vib[iElecHi][1] / hmi.H2_total));
486  fprintf(ioMEAN,PrintEfmt("%9.2e", pops_per_vib[iElecHi][2] / hmi.H2_total));
487  fprintf(ioMEAN,PrintEfmt("%9.2e", pops_per_vib[iElecHi][3] / hmi.H2_total));
488  fprintf( ioMEAN, "\n");
489 # endif
490  return;
491 }
492 
493 
494 /*H2_ReadTransprob read transition probabilities */
495 void H2_ReadTransprob( long int nelec )
496 {
497  const char* cdDATAFILE[N_H2_ELEC] =
498  {
499  "H2_transprob_X.dat",
500  "H2_transprob_B.dat",
501  "H2_transprob_C_plus.dat",
502  "H2_transprob_C_minus.dat",
503  "H2_transprob_B_primed.dat",
504  "H2_transprob_D_plus.dat",
505  "H2_transprob_D_minus.dat"
506  };
507  FILE *ioDATA;
508  char chLine[FILENAME_PATH_LENGTH_2];
509  long int i, n1, n2, n3;
510  long int iVibHi , iVibLo , iRotHi , iRotLo , iElecHi , iElecLo;
511  bool lgEOL;
512 
513  DEBUG_ENTRY( "H2_ReadTransprob()" );
514 
515  /* now open the data file */
516  ioDATA = open_data( cdDATAFILE[nelec], "r" );
517 
518  /* read the first line and check that magic number is ok */
519  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
520  {
521  fprintf( ioQQQ, " H2_ReadTransprob could not read first line of %s\n", cdDATAFILE[nelec]);
522  cdEXIT(EXIT_FAILURE);
523  }
524  i = 1;
525  /* level 1 magic number */
526  n1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
527  n2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
528  n3 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
529 
530  /* magic number
531  * the following is the set of numbers that appear at the start of level1.dat 01 08 10 */
532  if( ( n1 != 2 ) || ( n2 != 4 ) || ( n3 != 29 ) )
533  {
534  fprintf( ioQQQ,
535  " H2_ReadTransprob: the version of %s is not the current version.\n", cdDATAFILE[nelec] );
536  fprintf( ioQQQ,
537  " I expected to find the number 2 4 29 and got %li %li %li instead.\n" ,
538  n1 , n2 , n3 );
539  fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
540  cdEXIT(EXIT_FAILURE);
541  }
542 
543  /* read until not a comment */
544  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
545  BadRead();
546 
547  while( chLine[0]=='#' )
548  {
549  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
550  BadRead();
551  }
552  iVibHi = 1;
553  while( iVibHi >= 0 )
554  {
555  double Aul;
556  sscanf(chLine,"%li\t%li\t%li\t%li\t%li\t%li\t%le",
557  &iElecHi , &iVibHi ,&iRotHi , &iElecLo , &iVibLo , &iRotLo , &Aul );
558  ASSERT( iElecHi == nelec );
559  /* negative iVibHi says end of data */
560  if( iVibHi < 0 )
561  continue;
562 
563  ASSERT( iElecHi < N_H2_ELEC );
564  ASSERT( iElecLo < N_H2_ELEC );
566  ASSERT( iVibHi < 50 );
567  ASSERT( iVibLo < 50 );
568 
569  /* check that we actually included the levels in the model representation */
570  if( iVibHi <= h2.nVib_hi[iElecHi] &&
571  iVibLo <= h2.nVib_hi[iElecLo] &&
572  iRotHi <= h2.nRot_hi[iElecHi][iVibHi] &&
573  iRotLo <= h2.nRot_hi[iElecLo][iVibLo])
574  {
575  double ener = energy_wn[iElecHi][iVibHi][iRotHi] - energy_wn[iElecLo][iVibLo][iRotLo];
576 
577  /* only lines that have real Aul are added to stack. */
578  H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis = AddLine2Stack( true );
579 
580  H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->Aul = (realnum)Aul;
581  /* say that this line exists */
582  lgH2_line_exists[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] = true;
583  /* prints transitions with negative energies - should not happen */
584  if( ener <= 0. )
585  {
586  fprintf(ioQQQ,"negative energy H2 transition\t%li\t%li\t%li\t%li\t%.2e\t%.2e\n",
587  iVibHi,iVibLo,iRotHi,iRotLo,Aul,
588  H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyWN);
589  ShowMe();
590  }
591  }
592 # if 0
593  /* this prints all levels with As but without energies */
594  else
595  {
596  fprintf(ioQQQ,"no\t%li\t%li\t%li\t%li\t%.2e\n",
597  iVibHi,iVibLo,iRotHi,iRotLo,Aul);
598  }
599 # endif
600 
601  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
602  BadRead();
603  while( chLine[0]=='#' )
604  {
605  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
606  BadRead();
607  }
608  }
609  fclose( ioDATA );
610  return;
611 }
612 
613 #if 0
614 /*H2_Read_Cosmicray_distribution read distribution function for H2 population following cosmic ray collisional excitation */
615 void H2_Read_Cosmicray_distribution(void)
616 {
617  /*>>refer H2 cr excit Tine, S., Lepp, S., Gredel, R., & Dalgarno, A. 1997, ApJ, 481, 282 */
618  FILE *ioDATA;
619  char chLine[FILENAME_PATH_LENGTH_2];
620  long int i, n1, n2, n3, iVib , iRot;
621  long neut_frac;
622  bool lgEOL;
623 
624  DEBUG_ENTRY( "H2_Read_Cosmicray_distribution()" );
625 
626  /* now open the data file */
627  ioDATA = open_data( "H2_CosmicRay_collision.dat", "r" );
628 
629  /* read the first line and check that magic number is ok */
630  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
631  {
632  fprintf( ioQQQ, " H2_Read_Cosmicray_distribution could not read first line of %s\n", "H2_Cosmic_collision.dat");
633  cdEXIT(EXIT_FAILURE);
634  }
635 
636  i = 1;
637  /* level 1 magic number */
638  n1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
639  n2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
640  n3 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
641 
642  /* magic number
643  * the following is the set of numbers that appear at the start of H2_Cosmic_collision.dat 01 21 03 */
644  if( ( n1 != 1 ) || ( n2 != 21 ) || ( n3 != 3 ) )
645  {
646  fprintf( ioQQQ,
647  " H2_Read_Cosmicray_distribution: the version of %s is not the current version.\n", "H2_Cosmic_collision.dat" );
648  fprintf( ioQQQ,
649  " I expected to find the number 1 21 3 and got %li %li %li instead.\n" ,
650  n1 , n2 , n3 );
651  fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
652  cdEXIT(EXIT_FAILURE);
653  }
654 
655  /* read until not a comment */
656  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
657  BadRead();
658 
659  while( chLine[0]=='#' )
660  {
661  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
662  BadRead();
663  }
664 
665  iRot = 1;
666  iVib = 1;
667  neut_frac = 0;
668  while( iVib >= 0 )
669  {
670  long int j_minus_ji;
671  double a[10];
672 
673  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",
674  &iVib ,&j_minus_ji , &a[0],&a[1],&a[2],&a[3],&a[4],&a[5],&a[6],&a[7],&a[8],&a[9]
675  );
676  /* negative iVib says end of data */
677  if( iVib < 0 )
678  continue;
679 
680  /* cr_rate[CR_X][CR_VIB][CR_J][CR_EXIT];*/
681  /* check that we actually included the levels in the model representation */
682  ASSERT( iVib < CR_VIB );
683  ASSERT( j_minus_ji == -2 || j_minus_ji == +2 || j_minus_ji == 0 );
684  ASSERT( neut_frac < CR_X );
685 
686  /* now make i_minus_ji an array index */
687  j_minus_ji = 1 + j_minus_ji/2;
688  ASSERT( j_minus_ji>=0 && j_minus_ji<=2 );
689 
690  /* option to add Gaussian random mole */
691  for( iRot=0; iRot<CR_J; ++iRot )
692  {
693  cr_rate[neut_frac][iVib][iRot][j_minus_ji] = (realnum)a[iRot];
694  }
695  if( mole.lgH2_NOISECOSMIC )
696  {
697  realnum r;
699 
700  for( iRot=0; iRot<CR_J; ++iRot )
701  {
702  cr_rate[neut_frac][iVib][iRot][j_minus_ji] *= (realnum)pow(10.,(double)r);
703  }
704  }
705 
706  if( CR_PRINT )
707  {
708  fprintf(ioQQQ,"cr rate\t%li\t%li", iVib , j_minus_ji );
709  for( iRot=0; iRot<CR_J; ++iRot )
710  {
711  fprintf(ioQQQ,"\t%.3e", cr_rate[neut_frac][iVib][iRot][j_minus_ji] );
712  }
713  fprintf(ioQQQ,"\n" );
714  }
715 
716  /* now get next line */
717  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
718  BadRead();
719  while( chLine[0]=='#' )
720  {
721  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
722  BadRead();
723  }
724  }
725  fclose( ioDATA );
726 
727  return;
728 }
729 #endif
730 
731 /*H2_ReadEnergies read energies for all electronic levels */
732 void H2_ReadEnergies( long int nelec )
733 {
734  const char* cdDATAFILE[N_H2_ELEC] =
735  {
736  "H2_energy_X.dat",
737  "H2_energy_B.dat",
738  "H2_energy_C_plus.dat",
739  "H2_energy_C_minus.dat",
740  "H2_energy_B_primed.dat",
741  "H2_energy_D_plus.dat",
742  "H2_energy_D_minus.dat"
743  };
744  FILE *ioDATA;
745  char chLine[FILENAME_PATH_LENGTH_2];
746  long int i, n1, n2, n3, iVib , iRot;
747  bool lgEOL;
748 
749  DEBUG_ENTRY( "H2_ReadEnergies()" );
750 
751  /* now open the data file */
752  ioDATA = open_data( cdDATAFILE[nelec], "r" );
753 
754  /* read the first line and check that magic number is ok */
755  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
756  {
757  fprintf( ioQQQ, " H2_ReadEnergies could not read first line of %s\n", cdDATAFILE[nelec]);
758  cdEXIT(EXIT_FAILURE);
759  }
760  i = 1;
761  /* level 1 magic number */
762  n1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
763  n2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
764  n3 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
765 
766  /* magic number
767  * the following is the set of numbers that appear at the start of level1.dat 01 08 10 */
768  if( ( n1 != 2 ) || ( n2 != 4 ) || ( n3 != 29 ) )
769  {
770  fprintf( ioQQQ,
771  " H2_ReadEnergies: the version of %s is not the current version.\n", cdDATAFILE[nelec] );
772  fprintf( ioQQQ,
773  " I expected to find the number 2 4 29 and got %li %li %li instead.\n" ,
774  n1 , n2 , n3 );
775  fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
776  cdEXIT(EXIT_FAILURE);
777  }
778 
779  /* read until not a comment */
780  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
781  BadRead();
782 
783  while( chLine[0]=='#' )
784  {
785  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
786  BadRead();
787  }
788 
789  /* this will count the number of levels within each electronic state */
790  nLevels_per_elec[nelec] = 0;
791 
792  for( iVib=0; iVib<=h2.nVib_hi[nelec]; ++iVib )
793  {
794  for( iRot=h2.Jlowest[nelec]; iRot<=h2.nRot_hi[nelec][iVib]; ++iRot )
795  {
796  i = 1;
797  sscanf(chLine,"%li\t%li\t%le", &n1 , &n2 , &energy_wn[nelec][iVib][iRot] );
798  ASSERT( n1 == iVib );
799  ASSERT( n2 == iRot );
800 # if 0
801  /* in atomic units, or 1 Hartree, or two rydbergs */
802  if( nelec == 0 )
803  {
804  /* only do this for Phillip Stancil's file */
805  /* corrections are to get lowest rotation level to have energy of zero */
806  energy_wn[0][iVib][iRot] = -( energy_wn[0][iVib][iRot]- 3.6118114E+04 );
807  }
808 # endif
809  ASSERT( energy_wn[nelec][iVib][iRot]> 0. || (nelec==0 && iVib==0 && iRot==0 ) );
810  /* increment number of levels within this electronic state */
811  ++nLevels_per_elec[nelec];
812 
813  /* now start reading next line */
814  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
815  BadRead();
816  while( chLine[0]=='#' )
817  {
818  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
819  BadRead();
820  }
821  }
822  }
823  fclose( ioDATA );
824  return;
825 }
826 
827 /*H2_ReadDissprob read dissociation probabilities and kinetic energies for all electronic levels */
828 void H2_ReadDissprob( long int nelec )
829 {
830  const char* cdDATAFILE[N_H2_ELEC] =
831  {
832  "H2_dissprob_X.dat",/* this does not exist and nelec == 0 is not valid */
833  "H2_dissprob_B.dat",
834  "H2_dissprob_C_plus.dat",
835  "H2_dissprob_C_minus.dat",
836  "H2_dissprob_B_primed.dat",
837  "H2_dissprob_D_plus.dat",
838  "H2_dissprob_D_minus.dat"
839  };
840  FILE *ioDATA;
841  char chLine[FILENAME_PATH_LENGTH_2];
842  long int i, n1, n2, n3, iVib , iRot;
843  bool lgEOL;
844 
845  DEBUG_ENTRY( "H2_ReadDissprob()" );
846 
847  ASSERT( nelec > 0 );
848 
849  /* now open the data file */
850  ioDATA = open_data( cdDATAFILE[nelec], "r" );
851 
852  /* read the first line and check that magic number is ok */
853  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
854  {
855  fprintf( ioQQQ, " H2_ReadDissprob could not read first line of %s\n", cdDATAFILE[nelec]);
856  cdEXIT(EXIT_FAILURE);
857  }
858  i = 1;
859  /* level 1 magic number */
860  n1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
861  n2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
862  n3 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
863 
864  /* magic number
865  * the following is the set of numbers that appear at the start of level1.dat 01 08 10 */
866  if( ( n1 != 3 ) || ( n2 != 2 ) || ( n3 != 11 ) )
867  {
868  fprintf( ioQQQ,
869  " H2_ReadDissprob: the version of %s is not the current version.\n", cdDATAFILE[nelec] );
870  fprintf( ioQQQ,
871  " I expected to find the number 3 2 11 and got %li %li %li instead.\n" ,
872  n1 , n2 , n3 );
873  fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
874  cdEXIT(EXIT_FAILURE);
875  }
876 
877  /* read until not a comment */
878  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
879  BadRead();
880 
881  while( chLine[0]=='#' )
882  {
883  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
884  BadRead();
885  }
886 
887  for( iVib=0; iVib<=h2.nVib_hi[nelec]; ++iVib )
888  {
889  for( iRot=h2.Jlowest[nelec]; iRot<=h2.nRot_hi[nelec][iVib]; ++iRot )
890  {
891  double a, b;
892  i = 1;
893  sscanf(chLine,"%li\t%li\t%le\t%le",
894  &n1 , &n2 ,
895  /* dissociation probability */
896  &a ,
897  /* dissociation kinetic energy - eV not ergs */
898  &b);
899 
900  /* these have to agree if data file is valid */
901  ASSERT( n1 == iVib );
902  ASSERT( n2 == iRot );
903 
904  /* dissociation probability */
905  H2_dissprob[nelec][iVib][iRot] = (realnum)a;
906  /* dissociation kinetic energy - eV not ergs */
907  H2_disske[nelec][iVib][iRot] = (realnum)b;
908 
909  /* now get next line */
910  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
911  BadRead();
912  while( chLine[0]=='#' )
913  {
914  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
915  BadRead();
916  }
917  }
918  }
919  fclose( ioDATA );
920  return;
921 }
922 
923 
924 /*H2_Read_hminus_distribution read distribution function for H2 population following formation from H minus */
926 {
927  FILE *ioDATA;
928  char chLine[FILENAME_PATH_LENGTH_2];
929  long int i, n1, n2, n3, iVib , iRot;
930  bool lgEOL;
931  double sumrate[nTE_HMINUS];
932  /* set true for lots of printout */
933 # define H2HMINUS_PRT false
934 
935  DEBUG_ENTRY( "H2_Read_hminus_distribution()" );
936 
937  /* now open the data file */
938  ioDATA = open_data( "H2_hminus_deposit.dat", "r" );
939 
940  /* read the first line and check that magic number is ok */
941  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
942  {
943  fprintf( ioQQQ, " H2_Read_hminus_distribution could not read first line of %s\n", "H2_hminus_deposit.dat");
944  cdEXIT(EXIT_FAILURE);
945  }
946 
947  i = 1;
948  /* level 1 magic number */
949  n1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
950  n2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
951  n3 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
952 
953  /* magic number
954  * the following is the set of numbers that appear at the start of H2_hminus_deposit.dat 01 08 10 */
955  if( ( n1 != 2 ) || ( n2 != 10 ) || ( n3 != 17 ) )
956  {
957  fprintf( ioQQQ,
958  " H2_Read_hminus_distribution: the version of %s is not the current version.\n", "H2_hminus_deposit.dat" );
959  fprintf( ioQQQ,
960  " I expected to find the number 2 10 17 and got %li %li %li instead.\n" ,
961  n1 , n2 , n3 );
962  fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
963  cdEXIT(EXIT_FAILURE);
964  }
965 
966  /* read until not a comment */
967  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
968  BadRead();
969 
970  while( chLine[0]=='#' )
971  {
972  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
973  BadRead();
974  }
975 
976  /* convert temps to log */
977  for(i=0; i<nTE_HMINUS; ++i )
978  {
979  H2_te_hminus[i] = (realnum)log10(H2_te_hminus[i]);
980  sumrate[i] = 0.;
981  }
982 
983  iRot = 1;
984  iVib = 1;
985  while( iVib >= 0 )
986  {
987  /* set true to print rates */
988 
989  double a[nTE_HMINUS] , ener;
990  sscanf(chLine,"%li\t%li\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf",
991  &iVib ,&iRot , &ener, &a[0],&a[1],&a[2] , &a[3],&a[4],&a[5] ,&a[6]
992  );
993  /* negative iVib says end of data */
994  if( iVib < 0 )
995  continue;
996 
997  /* check that we actually included the levels in the model representation */
998  ASSERT( iVib <= h2.nVib_hi[0] &&
999  iRot <= h2.nRot_hi[0][iVib] );
1000 
1001  if( H2HMINUS_PRT )
1002  fprintf(ioQQQ,"hminusss\t%li\t%li", iVib , iRot );
1003  for( i=0; i<nTE_HMINUS; ++i )
1004  {
1005  H2_X_hminus_formation_distribution[i][iVib][iRot] = (realnum)pow(10.,-a[i]);
1006  sumrate[i] += H2_X_hminus_formation_distribution[i][iVib][iRot];
1007  if( H2HMINUS_PRT )
1008  fprintf(ioQQQ,"\t%.3e", H2_X_hminus_formation_distribution[i][iVib][iRot] );
1009  }
1010  if( H2HMINUS_PRT )
1011  fprintf(ioQQQ,"\n" );
1012 
1013  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
1014  BadRead();
1015  while( chLine[0]=='#' )
1016  {
1017  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
1018  BadRead();
1019  }
1020  }
1021  fclose( ioDATA );
1022 
1023  if( H2HMINUS_PRT )
1024  {
1025  /* print total rate */
1026  fprintf(ioQQQ," total H- formation rate ");
1027  /* convert temps to log */
1028  for(i=0; i<nTE_HMINUS; ++i )
1029  {
1030  fprintf(ioQQQ,"\t%.3e" , sumrate[i]);
1031  }
1032  fprintf(ioQQQ,"\n" );
1033  }
1034 
1035  /* convert to dimensionless factors that add to unity */
1036  for( iVib=0; iVib<=h2.nVib_hi[0]; ++iVib )
1037  {
1038  for( iRot=h2.Jlowest[0]; iRot<=h2.nRot_hi[0][iVib]; ++iRot )
1039  {
1040  for(i=0; i<nTE_HMINUS; ++i )
1041  {
1042  H2_X_hminus_formation_distribution[i][iVib][iRot] /= (realnum)sumrate[i];
1043  }
1044  }
1045  }
1046 
1047  if( H2HMINUS_PRT )
1048  {
1049  /* print total rate */
1050  fprintf(ioQQQ," H- distribution function ");
1051  for( iVib=0; iVib<=h2.nVib_hi[0]; ++iVib )
1052  {
1053  for( iRot=h2.Jlowest[0]; iRot<=h2.nRot_hi[0][iVib]; ++iRot )
1054  {
1055  fprintf(ioQQQ,"%li\t%li", iVib , iRot );
1056  for(i=0; i<nTE_HMINUS; ++i )
1057  {
1058  fprintf(ioQQQ,"\t%.3e", H2_X_hminus_formation_distribution[i][iVib][iRot] );
1059  }
1060  fprintf(ioQQQ,"\n" );
1061  }
1062  }
1063  }
1064  return;
1065 }
1066 
1067 /* ===================================================================== */
1068 /*H2_Punch_line_data punch line data for H2 molecule */
1070  /* io unit for punch */
1071  FILE* ioPUN ,
1072  /* punch all levels if true, only subset if false */
1073  bool lgDoAll )
1074 {
1075  long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo;
1076 
1077  if( !h2.lgH2ON )
1078  return;
1079 
1080  DEBUG_ENTRY( "H2_Punch_line_data()" );
1081 
1082  if( lgDoAll )
1083  {
1084  fprintf( ioQQQ,
1085  " H2_Punch_line_data ALL option not implemented in H2_Punch_line_data yet 1\n" );
1086  cdEXIT(EXIT_FAILURE);
1087  }
1088  else
1089  {
1090 
1091  /* punch line date, looping over all possible lines */
1092  for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi )
1093  {
1094  for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
1095  {
1096  for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
1097  {
1098  /* now the lower levels */
1099  /* NB - X is the only lower level considered here, since we are only
1100  * concerned with excited electronic levels as a photodissociation process
1101  * code exists to relax this assumption - simply change following to iElecHi */
1102  long int lim_elec_lo = 0;
1103  for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
1104  {
1105  /* want to include all vib states in lower level if different electronic level,
1106  * but only lower vib levels if same electronic level */
1107  long int nv = h2.nVib_hi[iElecLo];
1108  if( iElecLo==iElecHi )
1109  nv = iVibHi;
1110  for( iVibLo=0; iVibLo<=nv; ++iVibLo )
1111  {
1112  long nr = h2.nRot_hi[iElecLo][iVibLo];
1113  if( iElecLo==iElecHi && iVibHi==iVibLo )
1114  nr = iRotHi-1;
1115 
1116  for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
1117  {
1118  /* only punch if radiative transition exists */
1119  if( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].ipCont > 0 )
1120  {
1122  H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Coll.cs = 0.;
1123  /* print quantum indices */
1124  fprintf(ioPUN,"%2li %2li %2li %2li %2li %2li ",
1125  iElecHi,iVibHi,iRotHi,iElecLo,iVibLo,iRotLo );
1126  Punch1LineData( &H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] , ioPUN , false);
1127  }
1128  }
1129  }
1130  }
1131  }
1132  }
1133  }
1134  fprintf( ioPUN , "\n");
1135  }
1136  return;
1137 }
1138 
1139 /*H2_PunchLineStuff include H2 lines in punched optical depths, etc, called from PunchLineStuff */
1140 void H2_PunchLineStuff( FILE * io , realnum xLimit , long index)
1141 {
1142  long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo;
1143 
1144  if( !h2.lgH2ON )
1145  return;
1146 
1147  DEBUG_ENTRY( "H2_PunchLineStuff()" );
1148 
1149  /* loop over all possible lines */
1150  for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi )
1151  {
1152  for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
1153  {
1154  for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
1155  {
1156  /* now the lower levels */
1157  /* NB - X is the only lower level considered here, since we are only
1158  * concerned with excited electronic levels as a photodissociation process
1159  * code exists to relax this assumption - simply change following to iElecHi */
1160  long int lim_elec_lo = 0;
1161  for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
1162  {
1163  /* want to include all vib states in lower level if different electronic level,
1164  * but only lower vib levels if same electronic level */
1165  long int nv = h2.nVib_hi[iElecLo];
1166  if( iElecLo==iElecHi )
1167  nv = iVibHi;
1168  for( iVibLo=0; iVibLo<=nv; ++iVibLo )
1169  {
1170  long nr = h2.nRot_hi[iElecLo][iVibLo];
1171  if( iElecLo==iElecHi && iVibHi==iVibLo )
1172  nr = iRotHi-1;
1173 
1174  for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
1175  {
1176  /* only punch if radiative transition exists */
1177  if( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].ipCont > 0 )
1178  {
1179  pun1Line( &H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] , io , xLimit , index , 1.);
1180  }
1181  }
1182  }
1183  }
1184  }
1185  }
1186  }
1187 
1188  return;
1189 }
1190 
1191 
1192 /*H2_Prt_line_tau print line optical depths, called from premet in response to print line optical depths command*/
1194 {
1195  long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo;
1196 
1197  if( !h2.lgH2ON )
1198  return;
1199 
1200  DEBUG_ENTRY( "H2_Prt_line_tau()" );
1201 
1202  /* loop over all possible lines */
1203  for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi )
1204  {
1205  for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
1206  {
1207  for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
1208  {
1209  /* now the lower levels */
1210  /* NB - X is the only lower level considered here, since we are only
1211  * concerned with excited electronic levels as a photodissociation process
1212  * code exists to relax this assumption - simply change following to iElecHi */
1213  long int lim_elec_lo = 0;
1214  for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
1215  {
1216  /* want to include all vib states in lower level if different electronic level,
1217  * but only lower vib levels if same electronic level */
1218  long int nv = h2.nVib_hi[iElecLo];
1219  if( iElecLo==iElecHi )
1220  nv = iVibHi;
1221  for( iVibLo=0; iVibLo<=nv; ++iVibLo )
1222  {
1223  long nr = h2.nRot_hi[iElecLo][iVibLo];
1224  if( iElecLo==iElecHi && iVibHi==iVibLo )
1225  nr = iRotHi-1;
1226 
1227  for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
1228  {
1229  /* only print if radiative transition exists */
1230  if( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].ipCont > 0 )
1231  {
1232  prme(" c",&H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] );
1233  }
1234  }
1235  }
1236  }
1237  }
1238  }
1239  }
1240 
1241  return;
1242 }
1243 
1244 
1245 /*chMolBranch returns a char with the spectroscopic branch of a transition */
1246 STATIC char chMolBranch( long iRotHi , long int iRotLo )
1247 {
1248  /* these are the spectroscopic branches */
1249  char chBranch[5] = {'O','P','Q','R','S'};
1250  /* this is the index within the chBranch array */
1251  int ip = 2 + (iRotHi - iRotLo);
1252  if( ip<0 || ip>=5 )
1253  {
1254  fprintf(ioQQQ," chMolBranch called with insane iRotHi=%li iRotLo=%li ip=%i\n",
1255  iRotHi , iRotLo , ip );
1256  ip = 0;
1257  }
1258 
1259  return( chBranch[ip] );
1260 }
1261 
1262 /*H2_PunchDo punch some properties of the large H2 molecule */
1263 void H2_PunchDo( FILE* io , char chJOB[] , const char chTime[] , long int ipPun )
1264 {
1265  long int iVibHi , iElecHi , iRotHi , iVibLo , iElecLo , iRotLo,
1266  ip;
1267  long int iRot , iVib;
1268  long int LimVib , LimRot;
1269 
1270  DEBUG_ENTRY( "H2_PunchDo()" );
1271 
1272  /* which job are we supposed to do? This routine is active even when H2 is not turned on
1273  * so do not test on h2.lgH2ON initially */
1274 
1275  /* H2 populations computed in last zone -
1276  * give all of molecule in either matrix or triplet format */
1277  if( (strcmp( chJOB , "H2po" ) == 0) && (strcmp(chTime,"LAST") == 0) &&
1278  (punch.punarg[ipPun][2] != 0) )
1279  {
1280  /* >>chng 04 feb 19, do not punch if H2 not yet evaluated */
1281  if( h2.lgH2ON && hmi.lgBigH2_evaluated )
1282  {
1283  iVibHi= 0;
1284  iRotHi = 0;
1285  iElecHi=0;
1286  /* the limit to the number of vibration levels punched -
1287  * default is all, but first two numbers on punch h2 pops command
1288  * reset limit */
1289  /* this is limit to vibration */
1290  if( punch.punarg[ipPun][0] > 0 )
1291  {
1292  LimVib = (long)punch.punarg[ipPun][0];
1293  }
1294  else
1295  {
1296  LimVib = h2.nVib_hi[iElecHi];
1297  }
1298 
1299  /* first punch the current ortho, para, and total H2 density */
1300  fprintf(io,"%i\t%i\t%.3e\tortho\n",
1301  103 ,
1302  103 ,
1303  h2.ortho_density );
1304  fprintf(io,"%i\t%i\t%.3e\tpara\n",
1305  101 ,
1306  101 ,
1307  h2.para_density );
1308  fprintf(io,"%i\t%i\t%.3e\ttotal\n",
1309  0 ,
1310  0 ,
1311  hmi.H2_total );
1312 
1313  /* now punch the actual H2_populations, first part both matrix and triplets */
1314  for( iVibHi=0; iVibHi<=LimVib; ++iVibHi )
1315  {
1316  /* this is limit to rotation quantum index */
1317  if( punch.punarg[ipPun][1] > 0 )
1318  {
1319  LimRot = (long)punch.punarg[ipPun][1];
1320  }
1321  else
1322  {
1323  LimRot = h2.nRot_hi[iElecHi][iVibHi];
1324  }
1325  if( punch.punarg[ipPun][2] > 0 )
1326  {
1327  long int i;
1328  /* this option punch matrix */
1329  if( iVibHi == 0 )
1330  {
1331  fprintf(io,"vib\\rot");
1332  /* this is first vib, so make row of rot numbs */
1333  for( i=0; i<=LimRot; ++i )
1334  {
1335  fprintf(io,"\t%li",i);
1336  }
1337  fprintf(io,"\n");
1338  }
1339  fprintf(io,"%li",iVibHi );
1340  for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=LimRot; ++iRotHi )
1341  {
1342  fprintf(io,"\t%.3e",
1343  H2_populations[iElecHi][iVibHi][iRotHi]/hmi.H2_total );
1344  }
1345  fprintf(io,"\n" );
1346  }
1347  else if( punch.punarg[ipPun][2] < 0 )
1348  {
1349  /* this option punch triplets - the default */
1350  for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=LimRot; ++iRotHi )
1351  {
1352  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",
1353  /* upper vibration and rotation quantum numbers */
1354  iVibHi , iRotHi ,
1355  /* an 'O' or 'P' for ortho or para */
1356  chlgPara[H2_lgOrtho[iElecHi][iVibHi][iRotHi]],
1357  /* the level excitation energy in wavenumbers */
1358  energy_wn[iElecHi][iVibHi][iRotHi],
1359  /* actual population relative to total H2 */
1360  H2_populations[iElecHi][iVibHi][iRotHi]/hmi.H2_total ,
1361  /* old level H2_populations for comparison */
1362  H2_old_populations[iElecHi][iVibHi][iRotHi]/hmi.H2_total ,
1363  /* H2_populations per h2 and per statistical weight */
1364  H2_populations[iElecHi][iVibHi][iRotHi]/hmi.H2_total/H2_stat[iElecHi][iVibHi][iRotHi] ,
1365  /* LTE departure coefficient */
1366  /* >>chng 05 jan 26, missing factor of H2 abundance LTE is norm to unity, not tot abund */
1367  H2_populations[iElecHi][iVibHi][iRotHi]/SDIV(H2_populations_LTE[iElecHi][iVibHi][iRotHi]*hmi.H2_total ) ,
1368  /* fraction of exits that were collisional */
1369  H2_col_rate_out[iVibHi][iRotHi]/SDIV(H2_col_rate_out[iVibHi][iRotHi]+H2_rad_rate_out[0][iVibHi][iRotHi]) ,
1370  /* fraction of entries that were collisional */
1371  H2_col_rate_in[iVibHi][iRotHi]/SDIV(H2_col_rate_in[iVibHi][iRotHi]+H2_rad_rate_in[iVibHi][iRotHi]),
1372  /* collisions out */
1373  H2_col_rate_out[iVibHi][iRotHi],
1374  /* radiation out */
1375  H2_rad_rate_out[0][iVibHi][iRotHi] ,
1376  /* radiation out */
1377  H2_col_rate_in[iVibHi][iRotHi],
1378  /* radiation in */
1379  H2_rad_rate_in[iVibHi][iRotHi]
1380  );
1381  }
1382  }
1383  }
1384  }
1385  }
1386  /* punch H2 populations for each zone
1387  * H2_populations of v=0 for each zone */
1388  else if( (strcmp( chJOB , "H2po" ) == 0) && (strcmp(chTime,"LAST") != 0) &&
1389  (punch.punarg[ipPun][2] == 0) )
1390  {
1391  /* >>chng 04 feb 19, do not punch if h2 not yet evaluated */
1392  if( h2.lgH2ON && hmi.lgBigH2_evaluated )
1393  {
1394  fprintf(io,"%.5e\t%.3e\t%.3e", radius.depth_mid_zone ,
1396  /* rel pops of first two excited electronic states */
1397  fprintf(io,"\t%.3e\t%.3e",
1398  pops_per_elec[1] , pops_per_elec[2]);
1399  iElecHi = 0;
1400  iVibHi = 0;
1401  /* this is limit to vibration quantum index */
1402  if( punch.punarg[ipPun][0] > 0 )
1403  {
1404  LimVib = (long)punch.punarg[ipPun][1];
1405  }
1406  else
1407  {
1408  LimVib = h2.nRot_hi[iElecHi][iVibHi];
1409  }
1410  LimVib = MIN2( LimVib , h2.nVib_hi[iElecHi] );
1411  /* this is limit to rotation quantum index */
1412  if( punch.punarg[ipPun][1] > 0 )
1413  {
1414  LimRot = (long)punch.punarg[ipPun][1];
1415  }
1416  else
1417  {
1418  LimRot = h2.nRot_hi[iElecHi][iVibHi];
1419  }
1420  for( iVibHi = 0; iVibHi<=LimVib; ++iVibHi )
1421  {
1422  long int LimRotVib = MIN2( LimRot , h2.nRot_hi[iElecHi][iVibHi] );
1423  for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=LimRotVib; ++iRotHi )
1424  {
1425  fprintf(io,"\t%.3e",
1426  H2_populations[iElecHi][iVibHi][iRotHi]/hmi.H2_total );
1427  }
1428  fprintf(io,"\t");
1429  }
1430  fprintf(io,"\n");
1431  }
1432  }
1433 
1434  /* punch column densities */
1435  else if( (strcmp( chJOB , "H2cl" ) == 0) && (strcmp(chTime,"LAST") == 0) )
1436  {
1437  iVibHi= 0;
1438  iRotHi = 0;
1439  iElecHi=0;
1440  /* the limit to the number of vibration levels punched -
1441  * default is all, but first two numbers on punch h2 pops command
1442  * reset limit */
1443  /* this is limit to vibration */
1444  if( punch.punarg[ipPun][0] > 0 )
1445  {
1446  LimVib = (long)punch.punarg[ipPun][0];
1447  }
1448  else
1449  {
1450  LimVib = h2.nVib_hi[iElecHi];
1451  }
1452 
1453  /* first punch ortho and para H2_populations */
1454  fprintf(io,"%i\t%i\t%.3e\tortho\n",
1455  103 ,
1456  103 ,
1457  h2.ortho_colden );
1458  fprintf(io,"%i\t%i\t%.3e\tpara\n",
1459  101 ,
1460  101 ,
1461  h2.para_colden );
1462  /* total H2 column density */
1463  fprintf(io,"%i\t%i\t%.3e\ttotal\n",
1464  0 ,
1465  0 ,
1467 
1468  /* punch level column densities */
1469  for( iVibHi=0; iVibHi<=LimVib; ++iVibHi )
1470  {
1471  if( h2.lgH2ON )
1472  {
1473  /* this is limit to rotation quantum index */
1474  if( punch.punarg[ipPun][1] > 0 )
1475  {
1476  LimRot = (long)punch.punarg[ipPun][1];
1477  }
1478  else
1479  {
1480  LimRot = h2.nRot_hi[iElecHi][iVibHi];
1481  }
1482  if( punch.punarg[ipPun][2] > 0 )
1483  {
1484  long int i;
1485  /* punch matrix */
1486  if( iVibHi == 0 )
1487  {
1488  fprintf(io,"vib\\rot");
1489  /* this is first vib, so make row of rot numbs */
1490  for( i=0; i<=LimRot; ++i )
1491  {
1492  fprintf(io,"\t%li",i);
1493  }
1494  fprintf(io,"\n");
1495  }
1496  fprintf(io,"%li",iVibHi );
1497  for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=LimRot; ++iRotHi )
1498  {
1499  fprintf(io,"\t%.3e",
1500  H2_X_colden[iVibHi][iRotHi]/hmi.H2_total );
1501  }
1502  fprintf(io,"\n" );
1503  }
1504  else
1505  {
1506  /* punch triplets - the default */
1507  for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=LimRot; ++iRotHi )
1508  {
1509  fprintf(io,"%li\t%li\t%.1f\t%.3e\t%.3e\t%.3e\t%.3e\n",
1510  iVibHi ,
1511  iRotHi ,
1512  /* energy relative to 0,0, T1CM converts wavenumber to K */
1513  energy_wn[iElecHi][iVibHi][iRotHi]*T1CM,
1514  /* these are column densities for actual molecule */
1515  H2_X_colden[iVibHi][iRotHi] ,
1516  H2_X_colden[iVibHi][iRotHi]/H2_stat[iElecHi][iVibHi][iRotHi] ,
1517  /* these are same column densities but for LTE populations */
1518  H2_X_colden_LTE[iVibHi][iRotHi] ,
1519  H2_X_colden_LTE[iVibHi][iRotHi]/H2_stat[iElecHi][iVibHi][iRotHi]);
1520  }
1521  }
1522  }
1523  }
1524  }
1525  else if( (strcmp(chJOB , "H2pd" ) == 0) && (strcmp(chTime,"LAST") != 0) )
1526  {
1527  /* punch PDR
1528  * output some PDR information (densities, rates) for each zone */
1529  fprintf(io,"%.5e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
1530  /* depth in cm */
1532  /* the computed ortho and para densities */
1533  h2.ortho_density ,
1534  h2.para_density ,
1535  /* the Lyman Werner band dissociation, Tielens & Hollenbach */
1537  /* the Lyman Werner band dissociation, Bertoldi & Draine */
1539  /* the Lyman Werner band dissociation, big H2 mole */
1541  }
1542  else if( (strcmp(chJOB , "H2co" ) == 0) && (strcmp(chTime,"LAST") != 0) )
1543  {
1544  /* punch H2 cooling - do heating cooling for each zone old new H2 */
1545  fprintf(io,"%.5e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
1546  /* depth in cm */
1548  /* total cooling, equal to total heating */
1549  thermal.ctot ,
1550  /* H2 destruction by Solomon process, TH85 rate */
1552  /* H2 destruction by Solomon process, big H2 model rate */
1555  /* H2 photodissociation heating, eqn A9 of Tielens & Hollenbach 1985a */
1557  /* heating due to dissociation of electronic excited states */
1559  /* cooling (usually neg and so heating) due to collisions within X */
1562  );
1563 
1564  }
1565  else if( (strcmp(chJOB , "H2cr" ) == 0) && (strcmp(chTime,"LAST") != 0) )
1566  {
1567  /* PUNCH H2 CREATION - show H2 creation processes for each zone */
1568  /* >>chng 05 jul 15, TE, punch all H2 creation processes, unit cm-3 s-1 */
1569  fprintf(io,"%.5e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e",
1570  /* depth in cm */
1572  /* creation cm-3 s-1, destruction rate, s-1 */
1588  );
1589  fprintf(io,"\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e",
1590  /*chemical network*/
1591  /*light elements*/
1625  /* heavy elements */
1648  );
1649  fprintf(io,"\t%.3e\t%.3e\n",
1652  );
1653  }
1654  else if( (strcmp(chJOB , "H2ds" ) == 0) && (strcmp(chTime,"LAST") != 0) )
1655  {
1656  /* punch H2 destruction - show H2 destruction processes for each zone
1657  * >>chng 05 nov 17, TE, added the new reaction H2s + O+ -> OH+ + H
1658  * >>chng 05 oct 04, TE, remove eh2hhm(was double) and include dissociation by electrons, H2g/H2s + e -> 2H
1659  * >>chng 05 jul 15, TE, punch all H2 destruction rates, weighted by fractional abundance of H2g, H2s, unit s-1 */
1660  fprintf(io,"%.5e\t%.2e\t%.2e",
1661  /* depth in cm */
1663  /* total H2 creation rate, cm-3 s-1 */
1665  /* destruction rate, s-1 */
1666  hmi.H2_rate_destroy );
1667  /* this can be zero following a high-Te init temperature search
1668  * abort */
1669  if( hmi.H2_total > 0. )
1670  {
1671  fprintf(io,"\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e",
1672  /* H2 + e -> H- + H0 */
1674  /*photons*/
1679  /*electrons*/
1681  /*H+*/
1684  /*H0*/
1686  /*CR*/
1688  /*He+,HeH+*/
1692  /*H3+*/
1694  /*H2+*/
1696  /*H2s+H2g -> H2g + 2H*/
1698  /*H2g+H2g -> H2g + 2H*/
1700  /*H2s+H2s -> H2g + 2H*/
1702  /*H2s+H2s -> H2s + 2H*/
1704  /*chemical network*/
1705  /*light elements*/
1729  /*heavy elements*/
1745  );
1746  fprintf(io,"\t%.4e\t%.4e",
1747  /*H2g/Htot, H2s/Htot chemical network and from big molecule model*/
1750  );
1751  }
1752  fprintf(io,"\n");
1753  }
1754 
1755  else if( (strcmp(chJOB , "H2le" ) == 0) && (strcmp(chTime,"LAST") == 0) )
1756  {
1757  /* punch H2 levels */
1758  for( long int ipHi=0; ipHi < nLevels_per_elec[0]; ipHi++ )
1759  {
1760  long int ipVJ_Hi = H2_ipX_ener_sort[ipHi];
1761  iRotHi = ipRot_H2_energy_sort[ipVJ_Hi];
1762  iVibHi = ipVib_H2_energy_sort[ipVJ_Hi];
1763  double Asum , Csum[N_X_COLLIDER];
1764  long int nColl;
1765  Asum = 0;
1766  for( nColl=0; nColl<N_X_COLLIDER; ++nColl )
1767  Csum[nColl] = 0.;
1768  for( long int ipLo=0; ipLo<ipHi; ++ipLo )
1769  {
1770  /* all lower levels */
1771  long int ipVJ_Lo = H2_ipX_ener_sort[ipLo];
1772  iRotLo = ipRot_H2_energy_sort[ipVJ_Lo];
1773  iVibLo = ipVib_H2_energy_sort[ipVJ_Lo];
1774 
1775  /* radiative decays down */
1776  if( ( abs(iRotHi-iRotLo) == 2 || (iRotHi-iRotLo) == 0 ) && iVibLo <= iVibHi &&
1777  lgH2_line_exists[0][iVibHi][iRotHi][0][iVibLo][iRotLo] )
1778  {
1779  Asum += H2Lines[0][iVibHi][iRotHi][0][iVibLo][iRotLo].Emis->Aul*(
1780  H2Lines[0][iVibHi][iRotHi][0][iVibLo][iRotLo].Emis->Pesc +
1781  H2Lines[0][iVibHi][iRotHi][0][iVibLo][iRotLo].Emis->Pdest +
1782  H2Lines[0][iVibHi][iRotHi][0][iVibLo][iRotLo].Emis->Pelec_esc);
1783  }
1784  /* all collisions down */
1785  mr5ci H2cr = H2_CollRate.begin(iVibHi,iRotHi,iVibLo,iRotLo);
1786  for( nColl=0; nColl<N_X_COLLIDER; ++nColl )
1787  Csum[nColl] += H2cr[nColl];
1788  }
1789 
1790  /* punch H2 level energies */
1791  fprintf(io,"%li\t%li\t%.2f\t%li\t%.3e",
1792  iVibHi , iRotHi,
1793  energy_wn[0][iVibHi][iRotHi],
1794  (long)H2_stat[0][iVibHi][iRotHi],
1795  Asum );
1796  for( nColl=0; nColl<N_X_COLLIDER; ++nColl )
1797  /* sum over all lower levels */
1798  fprintf(io,"\t%.3e",Csum[nColl]);
1799  fprintf(io,"\n");
1800  }
1801  }
1802 
1803  else if( (strcmp(chJOB , "H2ra" ) == 0) && (strcmp(chTime,"LAST") != 0) )
1804  {
1805  /* punch h2 rates - some rates and lifetimes */
1806  double sumpop = 0. , sumlife = 0.;
1807 
1808  /* this block, find lifetime against photo excitation into excited electronic states */
1809  iElecLo = 0;
1810  iVibLo = 0;
1811  if( h2.lgH2ON && hmi.lgBigH2_evaluated )
1812  {
1813  for( iElecHi=1; iElecHi<mole.n_h2_elec_states; ++iElecHi )
1814  {
1815  for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
1816  {
1817  for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
1818  {
1819  for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=h2.nRot_hi[iElecLo][iVibLo]; ++iRotLo )
1820  {
1821  /* only do if radiative transition exists */
1822  if( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].ipCont > 0 )
1823  {
1824  sumlife +=
1825  H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->pump *
1826  H2_populations[iElecLo][iVibLo][iRotLo];
1827  sumpop +=
1828  H2_populations[iElecLo][iVibLo][iRotLo];
1829  }
1830  }
1831  }
1832  }
1833  }
1834  }
1835 
1836  /* continue output from punch h2 rates command */
1837  /* find photoexcitation rates from v=0 */
1838  /* PDR information for each zone */
1839  fprintf(io,
1840  "%.5e\t%.3e\t%.3e\t%.3e\t%li",
1841  /* depth in cm */
1843  /* the column density (cm^-2) in H2 */
1845  /* this is a special form of column density - should be proportional
1846  * to total shielding */
1848  /* visual extinction due to dust alone, of point source (star)*/
1850  /* number of large molecule evaluations in this zone */
1852  fprintf(io,
1853  "\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e",
1854  /* total H2 fraction */
1856  /* chemistry renorm factor */
1858  /* rate H2 forms on grains */
1860  /* rate H2 forms by H minus route */
1861  hmi.Hmolec[ipMHm]*1.35e-9,
1862  /* H2 destruction by Solomon process, TH85 rate */
1864  /* H2 destruction by Solomon process, Bertoldi & Draine rate */
1866  /* H2 destruction by Solomon process, Elwert et al. in preparation */
1868  /* H2 destruction by Solomon process, big H2 model rate */
1870  /* rate s-1 H2 electronic excit states decay into H2g */
1872  /* rate s-1 H2 electronic excit states decay into H2s */
1874  );
1875  fprintf(io,
1876  "\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e",
1877  /* The TH85 estimate of the radiation field relative to the Habing value */
1879  /* The DB96 estimate of the radiation field relative to the Habing value */
1881  /* cosmic ray ionization rate */
1882  secondaries.csupra[ipHYDROGEN][0]*0.93,
1883  sumlife/SDIV( sumpop ) ,
1888  fprintf(io,
1889  "\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
1894  thermal.htot);
1895  }
1896  /* punch h2 solomon */
1897  else if( (strcmp(chJOB , "H2so" ) == 0) && (strcmp(chTime,"LAST") != 0) )
1898  {
1899  /* remember as many as NSOL lines contributing to total Solomon process */
1900 # define NSOL 100
1901  double sum, one;
1902  long int jlosave[NSOL] , ivlosave[NSOL],
1903  iehisave[NSOL] ,jhisave[NSOL] , ivhisave[NSOL],
1904  nsave,
1905  ipOrdered[NSOL];
1906  int nFail;
1907  int i;
1908  realnum fsave[NSOL], wlsave[NSOL];
1909  /* Solomon process, and where it came from */
1910  fprintf(io,"%.5e\t%.3e",
1911  /* depth in cm */
1913  /* H2 destruction by Solomon process, big H2 model rate */
1916  sum = 0.;
1917  iElecLo = 0;
1918  /* find sum of all radiative exits from X into excited electronic states */
1919  if( h2.lgH2ON && hmi.lgBigH2_evaluated )
1920  {
1921  for( iElecHi=1; iElecHi<mole.n_h2_elec_states; ++iElecHi )
1922  {
1923  for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
1924  {
1925  for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
1926  {
1927  long int nv = h2.nVib_hi[iElecLo];
1928  for( iVibLo=0; iVibLo<=nv; ++iVibLo )
1929  {
1930  long nr = h2.nRot_hi[iElecLo][iVibLo];
1931  for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
1932  {
1933  /* only do if radiative transition exists */
1934  if( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].ipCont > 0 )
1935  {
1936  one = H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Lo->Pop *
1937  H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->pump;
1938  sum += one;
1939  }
1940  }
1941  }
1942  }
1943  }
1944  }
1945 
1946  /* make sure it is safe to div by sum */
1947  sum = SDIV( sum );
1948  nsave = 0;
1949  /* now loop back over X and print all those which contribute more than FRAC of the total */
1950 # define FRAC 0.01
1951  for( iElecHi=1; iElecHi<mole.n_h2_elec_states; ++iElecHi )
1952  {
1953  for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
1954  {
1955  for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
1956  {
1957  long int nv = h2.nVib_hi[iElecLo];
1958  for( iVibLo=0; iVibLo<=nv; ++iVibLo )
1959  {
1960  long nr = h2.nRot_hi[iElecLo][iVibLo];
1961  for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
1962  {
1963  /* only do if radiative transition exists */
1964  if( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].ipCont > 0 )
1965  {
1966  one = H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Lo->Pop *
1967  H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->pump;
1968  if( one/sum > FRAC && nsave<NSOL)
1969  {
1970  fsave[nsave] = (realnum)(one/sum);
1971  jlosave[nsave] = iRotLo;
1972  ivlosave[nsave] = iVibLo;
1973  jhisave[nsave] = iRotHi;
1974  ivhisave[nsave] = iVibHi;
1975  iehisave[nsave] = iElecHi;
1976  wlsave[nsave] = H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].WLAng;
1977  ++nsave;
1978  /*fprintf(io,"\t%li\t%li\t%li\t%li\t%li\t%.3f",
1979  iElecHi,iVibHi,iRotHi,iVibLo , iRotLo , one/sum );*/
1980  }
1981  }
1982  }
1983  }
1984  }
1985  }
1986  }/* iElecHi */
1987  /* now sort these into decreasing order */
1988 
1989  /* now sort by decreasing importance */
1990  /*spsort netlib routine to sort array returning sorted indices */
1991  spsort(
1992  /* input array to be sorted */
1993  fsave,
1994  /* number of values in x */
1995  nsave,
1996  /* permutation output array */
1997  ipOrdered,
1998  /* flag saying what to do - 1 sorts into increasing order, not changing
1999  * the original routine */
2000  -1,
2001  /* error condition, should be 0 */
2002  &nFail);
2003 
2004  /* print ratio of pumps to dissociations - this is 9:1 in TH85 */
2005  /*>>chng 05 jul 20, TE, punch average energy in H2s and renormalization factors for H2g and H2s */
2006  /* >>chng 05 sep 16, TE, chng denominator to do g and s with proper dissoc rates */
2007  fprintf(io,"\t%.3f\t%.3f\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e",
2008  /* this is sum of photons and CRs */
2011  /* this is sum of photons and CRs */
2017  );
2018  for( i=0; i<nsave; ++i )
2019  {
2020  ip = ipOrdered[i];
2021  /*lint -e644 not init */
2022  fprintf(io,"\t%li\t%li\t%li\t%li\t%li\t%.3f\t%.3f",
2023  iehisave[ip],ivhisave[ip],jhisave[ip],ivlosave[ip] , jlosave[ip] , fsave[ip] , wlsave[ip] );
2024  /*lint +e644 not init */
2025  }
2026  fprintf(io,"\n");
2027  }
2028  /*fprintf(io,"DEBUG tau\t%.3e\t%.3f\n",
2029  H2Lines[1][0][1][0][0][0].Emis->TauIn, H2Lines[1][0][1][0][0][0].WLAng); */
2030 # undef NSOL
2031  }
2032 
2033  else if( (strcmp(chJOB , "H2te" ) == 0) && (strcmp(chTime,"LAST") != 0) )
2034  {
2035  /* punch h2 temperatures */
2036  double pop_ratio10,pop_ratio20,pop_ratio30,pop_ratio31,pop_ratio40;
2037  double T10,T20,T30,T31,T40;
2038  /* subscript"sum" denotes integrated quantities */
2039  double T10_sum,T20_sum,T30_sum,T31_sum,T40_sum;
2040  double pop_ratio10_sum,pop_ratio20_sum,pop_ratio30_sum,pop_ratio31_sum,pop_ratio40_sum;
2041  if( h2.lgH2ON && h2.nCallH2_this_zone )
2042  {
2043  double energyK = T1CM*(energy_wn[0][0][1] - energy_wn[0][0][0]);
2044  /* the ratio of H2_populations of J=1 to 0 */
2045  pop_ratio10 = H2_populations[0][0][1]/SDIV(H2_populations[0][0][0]);
2046  pop_ratio10_sum = H2_X_colden[0][1]/SDIV(H2_X_colden[0][0]);
2047  /* the corresponding temperature */
2048  T10 = -170.5/log(SDIV(pop_ratio10) * H2_stat[0][0][0]/H2_stat[0][0][1]);
2049  T10_sum = -170.5/log(SDIV(pop_ratio10_sum) * H2_stat[0][0][0]/H2_stat[0][0][1]);
2050 
2051  energyK = T1CM*(energy_wn[0][0][2] - energy_wn[0][0][0]);
2052  pop_ratio20 = H2_populations[0][0][2]/SDIV(H2_populations[0][0][0]);
2053  T20 = -energyK/log(SDIV(pop_ratio20) * H2_stat[0][0][0]/H2_stat[0][0][2]);
2054 
2055  pop_ratio20_sum = H2_X_colden[0][2]/SDIV(H2_X_colden[0][0]);
2056  T20_sum = -energyK/log(SDIV(pop_ratio20_sum) * H2_stat[0][0][0]/H2_stat[0][0][2]);
2057 
2058  energyK = T1CM*(energy_wn[0][0][3] - energy_wn[0][0][0]);
2059  pop_ratio30 = H2_populations[0][0][3]/SDIV(H2_populations[0][0][0]);
2060  T30 = -energyK/log(SDIV(pop_ratio30) * H2_stat[0][0][0]/H2_stat[0][0][3]);
2061 
2062  pop_ratio30_sum = H2_X_colden[0][3]/SDIV(H2_X_colden[0][0]);
2063  T30_sum = -energyK/log(SDIV(pop_ratio30_sum) * H2_stat[0][0][0]/H2_stat[0][0][3]);
2064 
2065  energyK = T1CM*(energy_wn[0][0][3] - energy_wn[0][0][1]);
2066  pop_ratio31 = H2_populations[0][0][3]/SDIV(H2_populations[0][0][1]);
2067  T31 = -energyK/log(SDIV(pop_ratio31) * H2_stat[0][0][1]/H2_stat[0][0][3]);
2068 
2069  pop_ratio31_sum = H2_X_colden[0][3]/SDIV(H2_X_colden[0][1]);
2070  T31_sum = -energyK/log(SDIV(pop_ratio31_sum) * H2_stat[0][0][1]/H2_stat[0][0][3]);
2071 
2072  energyK = T1CM*(energy_wn[0][0][4] - energy_wn[0][0][0]);
2073  pop_ratio40 = H2_populations[0][0][4]/SDIV(H2_populations[0][0][0]);
2074  T40 = -energyK/log(SDIV(pop_ratio40) * H2_stat[0][0][0]/H2_stat[0][0][4]);
2075 
2076  pop_ratio40_sum = H2_X_colden[0][4]/SDIV(H2_X_colden[0][0]);
2077  T40_sum = -energyK/log(SDIV(pop_ratio40_sum) * H2_stat[0][0][0]/H2_stat[0][0][4]);
2078  }
2079  else
2080  {
2081  pop_ratio10 = 0.;
2082  pop_ratio10_sum = 0.;
2083  T10 = 0.;
2084  T20 = 0.;
2085  T30 = 0.;
2086  T31 = 0.;
2087  T40 = 0.;
2088  T10_sum = 0.;
2089  T20_sum = 0.;
2090  T30_sum = 0.;
2091  T31_sum = 0.;
2092  T40_sum = 0.;
2093  }
2094 
2095  /* various temperatures for neutral/molecular gas */
2096  fprintf( io,
2097  "%.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" ,
2098  /* depth in cm */
2100  /* total H2 fraction */
2102  /* ratio of H2_populations of 1 to 0 only */
2103  pop_ratio10 ,
2104  /* sum of all ortho and para */
2106  T10,T20,T30,T31,T40,
2107  phycon.te ,
2108  hyperfine.Tspin21cm,T10_sum,T20_sum,T30_sum,T31_sum,T40_sum );
2109  }
2110  else if( (strcmp(chJOB , "H2ln" ) == 0) && (strcmp(chTime,"LAST") == 0) )
2111  {
2112  /* punch H2 lines - output the full emission-line spectrum */
2113  double thresh;
2114  double renorm;
2115  /* first test, is H2 turned on? Second test, have lines arrays
2116  * been set up - nsum is negative if abort occurs before lines
2117  * are set up */
2118  if( h2.lgH2ON && LineSave.nsum > 0)
2119  {
2120  ASSERT( LineSave.ipNormWavL >= 0 );
2121  /* get the normalization line */
2124  else
2125  renorm = 1.;
2126 
2127  if( renorm > SMALLFLOAT )
2128  {
2129  /* this is threshold for faintest line, normally 0, set with
2130  * number on punch H2 command */
2131  thresh = thresh_punline_h2/(realnum)renorm;
2132  }
2133  else
2134  thresh = 0.f;
2135 
2136  /* punch H2 line intensities at end of iteration
2137  * h2.nElecLevelOutput is electronic level with 1 for ground, so this loop is < h2.nElecLevelOutput */
2138  for( iElecHi=0; iElecHi < h2.nElecLevelOutput; ++iElecHi )
2139  {
2140  for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
2141  {
2142  for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
2143  {
2144  /* now the lower levels */
2145  /* NB - X is the only lower level considered here, since we are only
2146  * concerned with excited electronic levels as a photodissociation process
2147  * code exists to relax this assumption - simply change following to iElecHi */
2148  long int lim_elec_lo = 0;
2149  for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
2150  {
2151  long int nv = h2.nVib_hi[iElecLo];
2152  if( iElecLo==iElecHi )
2153  nv = iVibHi;
2154  for( iVibLo=0; iVibLo<=nv; ++iVibLo )
2155  {
2156  long nr = h2.nRot_hi[iElecLo][iVibLo];
2157  if( iElecLo==iElecHi && iVibHi==iVibLo )
2158  nr = iRotHi-1;
2159  for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<nr; ++iRotLo )
2160  {
2161  if( H2_SaveLine[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] > thresh )
2162  {
2163  /* air wavelength in microns */
2164  /* WLAng contains correction for index of refraction of air */
2165  double wl = H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].WLAng/1e4;
2166  /*ASSERT( abs(iRotHi-iRotLo)<=2 );*/
2167 
2168  fprintf(io, "%li-%li %c(%li)",
2169  iVibHi ,
2170  iVibLo ,
2171  chMolBranch( iRotHi , iRotLo ) ,
2172  iRotLo );
2173  fprintf( io, "\t%ld\t%ld\t%ld\t%ld\t%ld\t%ld",
2174  iElecHi , iVibHi , iRotHi , iElecLo , iVibLo , iRotLo);
2175  /* WLAng contains correction for index of refraction of air */
2176  fprintf( io, "\t%.7f\t", wl );
2177  /*prt_wl print floating wavelength in Angstroms, in output format */
2178  prt_wl( io , H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].WLAng );
2179  fprintf( io, "\t%.3f\t%.3e",
2180  log10(MAX2(1e-37,H2_SaveLine[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo])) + radius.Conv2PrtInten,
2181  H2_SaveLine[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo]*renorm );
2182  /* excitation energy of upper level in K */
2183  fprintf( io, "\t%.3f", energy_wn[iElecHi][iVibHi][iRotHi]*T1CM );
2184  /* the product h nu * Aul */
2185  ASSERT( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].ipCont > 0 );
2186  fprintf( io, "\t%.3e", H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->Aul*
2187  H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyErg *
2188  H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Hi->g);
2189  fprintf( io, "\n");
2190  }
2191  }
2192  }
2193  }
2194  }
2195  }
2196  }
2197  }
2198  }
2199  else if( (strcmp(chJOB , "H2sp" ) == 0) )
2200  {
2201  iVib = 0;
2202  iRot = 0;
2203 # if 0
2204  /* punch h2 special */
2205  fprintf(io,"%.4e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
2207  H2_populations[0][iVib][iRot] ,
2208  radius.depth_mid_zone * H2_populations[0][iVib][iRot] ,
2209  H2_rad_rate_out[0][iVib][iRot] ,
2210  H2_rad_rate_in[iVib][iRot] ,
2211  H2_col_rate_out_ old[iVib][iRot] ,
2212  H2_col_rate_in_ old[iVib][iRot] );
2213 # endif
2214  fprintf(io,"%.4e\t%.2e\t%.2e\t%.2e\t%.2e\n",
2216  H2_populations[0][iVib][iRot] ,
2217  H2Lines[1][1][1][0][iVib][iRot].Emis->pump,
2218  H2Lines[1][1][1][0][iVib][iRot].Emis->TauIn,
2219  H2Lines[1][1][1][0][iVib][iRot].Emis->TauCon);
2220  }
2221  return;
2222 }
2223  /*cdH2_Line determines intensity and luminosity of and H2 line. The first
2224  * six arguments give the upper and lower quantum designation of the levels.
2225  * The function returns 1 if we found the line,
2226  * and false==0 if we did not find the line because ohoto-para transition
2227  * or upper level has lower energy than lower level */
2228 long int cdH2_Line(
2229  /* indices for the upper level */
2230  long int iElecHi,
2231  long int iVibHi ,
2232  long int iRotHi ,
2233  /* indices for lower level */
2234  long int iElecLo,
2235  long int iVibLo ,
2236  long int iRotLo ,
2237  /* linear intensity relative to normalization line*/
2238  double *relint,
2239  /* log of luminosity or intensity of line */
2240  double *absint )
2241 {
2242 
2243  DEBUG_ENTRY( "cdH2_Line()" );
2244 
2245  /* these will be return values if we can't find the line */
2246  *relint = 0.;
2247  *absint = 0.;
2248 
2249  /* for now both electronic levels must be zero */
2250  if( iElecHi!=0 || iElecLo!=0 )
2251  {
2252  return 0;
2253  }
2254 
2255  /* check that energy of first level is higher than energy of second level */
2256  if( energy_wn[iElecHi][iVibHi][iRotHi] < energy_wn[iElecLo][iVibHi][iRotHi] )
2257  {
2258  return 0;
2259  }
2260 
2261  /* check that ortho-para does not change */
2262  if( H2_lgOrtho[iElecHi][iVibHi][iRotHi] - H2_lgOrtho[iElecLo][iVibLo][iRotLo] != 0 )
2263  {
2264  return 0;
2265  }
2266 
2267  /* exit if lines does not exist */
2268  if( !lgH2_line_exists[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] )
2269  {
2270  return 0;
2271  }
2272 
2273  ASSERT( LineSave.ipNormWavL >= 0 );
2274  /* does the normalization line have a positive intensity*/
2276  {
2277  *relint = H2_SaveLine[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo]/
2279  }
2280  else
2281  {
2282  *relint = 0.;
2283  }
2284 
2285  /* return log of line intensity if it is positive */
2286  if( H2_SaveLine[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] > 0. )
2287  {
2288  *absint = log10(H2_SaveLine[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo]) +
2290  }
2291  else
2292  {
2293  /* line intensity is actually zero, return small number */
2294  *absint = -37.;
2295  }
2296  /* this indicates success */
2297  return 1;
2298 }

Generated for cloudy by doxygen 1.8.3.1