cloudy trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /*H2_ParsePunch parse the punch h2 command */ 00004 /*H2_PunchDo punch some properties of the large H2 molecule */ 00005 /*chMolBranch returns a char with the spectroscopic branch of a transition */ 00006 /*H2_Prt_line_tau print line optical depths, called from premet in response to print line optical depths command*/ 00007 /*H2_PunchLineStuff include H2 lines in punched optical depths, etc, called from PunchLineStuff */ 00008 /*H2_Punch_line_data punch line data for H2 molecule */ 00009 /*H2_Read_hminus_distribution read distribution function for H2 population following formation from H minus */ 00010 /*H2_ReadDissprob read dissociation probabilities and kinetic energies for all electronic levels */ 00011 /*H2_ReadEnergies read energies for all electronic levels */ 00012 /*H2_ReadTransprob read transition probabilities */ 00013 /*H2_Prt_Zone print H2 info into zone results, called from prtzone for each printed zone */ 00014 /*H2_ParsePunch parse the punch h2 command */ 00015 /*H2_Prt_column_density print H2 info into zone results, called from prtzone for each printed zone */ 00016 /*H2_LinesAdd add in explicit lines from the large H2 molecule, called by lines_molecules */ 00017 /*cdH2_Line returns 1 if we found the line, 00018 * or false==0 if we did not find the line because ohoto-para transition 00019 * or upper level has lower energy than lower level */ 00020 #include "cddefines.h" 00021 #include "physconst.h" 00022 #include "punch.h" 00023 #include "hmi.h" 00024 #include "prt.h" 00025 #include "secondaries.h" 00026 #include "grainvar.h" 00027 #include "phycon.h" 00028 #include "rfield.h" 00029 #include "hyperfine.h" 00030 #include "thermal.h" 00031 #include "lines.h" 00032 #include "lines_service.h" 00033 #include "dense.h" 00034 #include "radius.h" 00035 #include "colden.h" 00036 #include "taulines.h" 00037 #include "h2.h" 00038 #include "h2_priv.h" 00039 #include "cddrive.h" 00040 #include "mole.h" 00041 00042 /* this will say whether ortho or para, 00043 * H2_lgOrtho is 0 or 1 depending on whether or not ortho, 00044 * so chlgPara[H2_lgOrtho] gives P or O for printing */ 00045 static char chlgPara[2]={'P','O'}; 00046 00047 /* intensity, relative to normalization line, for faintest line to punch */ 00048 static realnum thresh_punline_h2; 00049 00050 /*H2_LinesAdd add in explicit lines from the large H2 molecule, called by lines_molecules */ 00051 void H2_LinesAdd(void) 00052 { 00053 /* these are the quantum designations of the lines we will output */ 00054 int iRotHi, iVibHi, iElecHi ,iRotLo, iVibLo, iElecLo; 00055 00056 /* H2 not on, so space not allocated */ 00057 if( !h2.lgH2ON ) 00058 return; 00059 00060 DEBUG_ENTRY( "H2_LinesAdd()" ); 00061 00062 /* >>chng 05 nov 04, make info copies of these lines up here 00063 * these are among the strongest lines in the 2 micron window and some have nearly the same 00064 * wavelength as far weaker lines that may come before them in the line stack. in that case 00065 * cdLine would find the much weaker line with the same wavelength. 00066 * put strong H2 lines first so that line search will find these, and not far weaker 00067 * lines with nearly the same wavelength - these will be duplicated in the output but 00068 * these are here for into (the 'i) so does no harm 00069 * 00070 * the array indices in the following structures give upper and lower electronic, vib, rot quantum 00071 * indices. 00072 * H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] 00073 * >>chng 05 dec 22, had hand entered wavelength in A as second parameter. This gave 00074 * rounded off result when set line precision 5 was used. now uses same logic that 00075 * PutLine will eventually use - simply enter same wl in Ang 00076 * 1-0 S(4) - 18910 */ 00077 lindst( H2Lines[0][1][6][0][0][4].Emis->xIntensity , H2Lines[0][1][6][0][0][4].WLAng , "H2 " , 00078 H2Lines[0][1][6][0][0][4].ipCont,'i',false , 00079 "H2 line"); 00080 /* 1-0 S(3) - 19570 */ 00081 lindst( H2Lines[0][1][5][0][0][3].Emis->xIntensity , H2Lines[0][1][5][0][0][3].WLAng , "H2 " , 00082 H2Lines[0][1][5][0][0][3].ipCont,'i',false , 00083 "H2 line"); 00084 /* 1-0 S(2) - 20330 */ 00085 lindst( H2Lines[0][1][4][0][0][2].Emis->xIntensity , H2Lines[0][1][4][0][0][2].WLAng , "H2 " , 00086 H2Lines[0][1][4][0][0][2].ipCont,'i',false , 00087 "H2 line"); 00088 /* 1-0 S(1) - 21210 */ 00089 lindst( H2Lines[0][1][3][0][0][1].Emis->xIntensity , H2Lines[0][1][3][0][0][1].WLAng , "H2 " , 00090 H2Lines[0][1][3][0][0][1].ipCont,'i',false , 00091 "H2 line"); 00092 /* 1-0 S(0) - 22230 */ 00093 lindst( H2Lines[0][1][2][0][0][0].Emis->xIntensity , H2Lines[0][1][2][0][0][0].WLAng , "H2 " , 00094 H2Lines[0][1][2][0][0][0].ipCont,'i',false , 00095 "H2 line"); 00096 /* start Q branch - selection rule requires that J be non-zero, so no Q(0) */ 00097 /* 1-0 Q(2) - 24130 */ 00098 lindst( H2Lines[0][1][2][0][0][2].Emis->xIntensity , H2Lines[0][1][2][0][0][2].WLAng , "H2 " , 00099 H2Lines[0][1][2][0][0][2].ipCont,'i',false , 00100 "H2 line"); 00101 /* 1-0 Q(1) - 24060 */ 00102 lindst( H2Lines[0][1][1][0][0][1].Emis->xIntensity , H2Lines[0][1][1][0][0][1].WLAng , "H2 " , 00103 H2Lines[0][1][1][0][0][1].ipCont,'i',false , 00104 "H2 line"); 00105 00106 /* print all lines from lowest n levels within X */ 00107 /* loop over all possible lines and set H2_populations, 00108 * and quantities that depend on them */ 00109 for( iElecHi=0; iElecHi<h2.nElecLevelOutput; ++iElecHi ) 00110 { 00111 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi ) 00112 { 00113 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi ) 00114 { 00115 long int lim_elec_lo = 0; 00116 /* now the lower levels */ 00117 /* NB - X is the only lower level considered here, since we are only 00118 * concerned with excited electronic levels as a photodissociation process 00119 * code exists to relax this assumption - simply change following to iElecHi */ 00120 for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo ) 00121 { 00122 /* want to include all vib states in lower level if different electronic level, 00123 * but only lower vib levels if same electronic level */ 00124 long int nv = h2.nVib_hi[iElecLo]; 00125 if( iElecLo==iElecHi ) 00126 nv = iVibHi; 00127 for( iVibLo=0; iVibLo<=nv; ++iVibLo ) 00128 { 00129 long nr = h2.nRot_hi[iElecLo][iVibLo]; 00130 if( iElecLo==iElecHi && iVibHi==iVibLo ) 00131 nr = iRotHi-1; 00132 00133 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo ) 00134 { 00135 if( lgH2_line_exists[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] ) 00136 { 00137 /* all ground vib state rotation lines - first is J to J-2 */ 00138 PutLine(&H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo], 00139 "H2 lines"); 00140 if( LineSave.ipass == 0 ) 00141 { 00142 H2_SaveLine[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] = 0.; 00143 } 00144 else if( LineSave.ipass == 1 ) 00145 { 00146 H2_SaveLine[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] += (realnum)( 00147 radius.dVeff*H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->xIntensity); 00148 } 00149 } 00150 } 00151 } 00152 } 00153 } 00154 } 00155 } 00156 return; 00157 } 00158 00159 /*H2_ParsePunch parse the punch h2 command */ 00160 void H2_ParsePunch( char *chCard , 00161 char *chHeader) 00162 { 00163 long int i , nelem; 00164 bool lgEOL; 00165 00166 DEBUG_ENTRY( "H2_ParsePunch()" ); 00167 00168 i = 5; 00169 nelem = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00170 if( nelem != 2 ) 00171 { 00172 fprintf( ioQQQ, " The first number on this line must be the 2 in H2\n Sorry.\n" ); 00173 cdEXIT(EXIT_FAILURE); 00174 } 00175 /* this provides info on the large H2 molecule */ 00176 if( nMatch("COLU",chCard) ) 00177 { 00178 /* punch column density */ 00179 strcpy( punch.chPunch[punch.npunch], "H2cl" ); 00180 00181 /* this is an option to scan off highest vib and rot states 00182 * to punch pops - first is limit to vibration, then rotation 00183 * if no number is entered then 0 is set and all levels punched */ 00184 /* now get vib limit */ 00185 punch.punarg[punch.npunch][0] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00186 00187 /* highest rotation */ 00188 punch.punarg[punch.npunch][1] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00189 /* this says whether to punch triplets or a matrix for output - 00190 * default is triplets, so only check for matrix */ 00191 if( nMatch( "MATR" , chCard ) ) 00192 { 00193 /* matrix */ 00194 punch.punarg[punch.npunch][2] = 1; 00195 sprintf( chHeader, "#vib\trot\tcolumn density\n" ); 00196 } 00197 else 00198 { 00199 /* triplets */ 00200 punch.punarg[punch.npunch][2] = -1; 00201 sprintf( chHeader, "#vib\trot\tEner(K)\tcolden\tcolden/stat wght\tLTE colden\tLTE colden/stat wght\n" ); 00202 } 00203 } 00204 else if( nMatch("COOL",chCard) ) 00205 { 00206 /* heating and cooling rates */ 00207 strcpy( punch.chPunch[punch.npunch], "H2co" ); 00208 sprintf( chHeader, 00209 "#H2 depth\ttot cool\tTH Sol\tBig Sol\tTH pht dis\tpht dis\tTH Xcool\tXcool \n" ); 00210 } 00211 00212 else if( nMatch("CREA",chCard) ) 00213 { 00214 /* H2 creation rates */ 00215 strcpy( punch.chPunch[punch.npunch], "H2cr" ); 00216 sprintf( chHeader, 00217 "#H2 depth\tH2_rate_create\tH2_rate_destroy\trate_h2_form_grains_used_total\tassoc_detach" 00218 "\tbh2dis\tbh2h2p\tradasc\th3ph2p\th2phmh2h\tbh2h22hh2\th3phmh2hh\th3phm2h2\th32h2\teh3_h2h\th3ph2hp" 00219 "\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" 00220 "\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" 00221 "\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" 00222 "\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" 00223 "\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" 00224 "\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" 00225 "\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" 00226 "\tH3P_Cl_HClP_H2\tH3P_HCl_H2ClP_H2\tH2P_C2_C2P_H2\tHminus_NH4P_NH3_H2\tH3P_HCN_HCNHP_H2" 00227 "\tdestruction/creation\tav Einstein A\n"); 00228 } 00229 else if( nMatch("DEST",chCard) ) 00230 { 00231 /* punch H2 destruction - output destruction rates */ 00232 strcpy( punch.chPunch[punch.npunch], "H2ds" ); 00233 sprintf( chHeader, 00234 "#depth\ttot H2 rate create\ttot H2 rate destroy\ttot H- backwards\tSolomon H2g\tSolomon H2s\tphotodissoc H2s\tphotodissoc H2g" 00235 "\te- dissoc\trh2h2p\th2hph3p\tH0 dissoc\tCR\trheph2hpheh\theph2heh2p\thehph2h3phe\th3petc\tH2Ph3p" 00236 "\th2sh2g\th2h22hh2\th2sh2sh2g2h\th2sh2sh2s2h\tH2_CHP_CH2P_H\tH2_CH2P_CH3P_H\tH2_OHP_H2OP_H\tH2_H2OP_H3OP_H\tH2_COP_HCOP_H" 00237 "\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" 00238 "\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" 00239 "\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" 00240 "\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" 00241 "\tH2_ClP_HClP_H\tH2_HClP_H2ClP_H\tH2_HCNP_HCNHP_H" 00242 "\tfrac H2g\tfrac H2s\n"); 00243 } 00244 00245 else if( nMatch("HEAT",chCard) ) 00246 { 00247 /* heating and cooling rates */ 00248 strcpy( punch.chPunch[punch.npunch], "H2he" ); 00249 sprintf( chHeader, 00250 "#H2 depth\ttot Heat\tHeat(big)\tHeat(TH85)\tDissoc(Big)\tDissoc(TH85) \n" ); 00251 } 00252 00253 else if( nMatch("LEVE",chCard) ) 00254 { 00255 /* punch H2 level energies */ 00256 strcpy( punch.chPunch[punch.npunch], "H2le" ); 00257 sprintf( chHeader, 00258 "#H2 v\tJ\tenergy(wn)\tstat wght\tSum As" ); 00259 char chHoldit[chN_X_COLLIDER+12]; 00260 for( int nColl=0; nColl<N_X_COLLIDER; ++nColl ) 00261 { 00262 /* labels for all colliders */ 00263 sprintf(chHoldit,"\tCritDen %s",chH2ColliderLabels[nColl]); 00264 strcat( chHeader , chHoldit ); 00265 } 00266 strcat( chHeader , "\n" ); 00267 } 00268 00269 else if( nMatch("LINE",chCard) ) 00270 { 00271 /* punch H2 lines - all in X */ 00272 strcpy( punch.chPunch[punch.npunch], "H2ln" ); 00273 sprintf( chHeader, 00274 "#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" ); 00275 /* first optional number changes the threshold of weakest line to print*/ 00276 /* fe2thresh is intensity relative to normalization line, 00277 * normally Hbeta, and is set to zero in zero.c */ 00278 00279 /* threshold for faintest line to punch, default is 1e-4 of norm line */ 00280 thresh_punline_h2 = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00281 if( lgEOL ) 00282 { 00283 /* this will be default relative intensity for faintest line to punch */ 00284 thresh_punline_h2 = 1e-4f; 00285 } 00286 00287 /* it is a log if negative */ 00288 if( thresh_punline_h2 < 0. ) 00289 { 00290 thresh_punline_h2 = (realnum)pow((realnum)10.f,thresh_punline_h2); 00291 } 00292 00293 /* lines from how many electronic states? default is one, just X, and is 00294 * obtained with GROUND keyword. ALL will produce all lines from all levels. 00295 * else, if a number is present, will be the number. if no number, no keyword, 00296 * appear then just ground */ 00297 if( nMatch( "ELEC" , chCard ) ) 00298 { 00299 if( nMatch(" ALL",chCard) ) 00300 { 00301 /* all electronic levels - when done, will set upper limit, the 00302 * number of electronic levels actually computed, don't know this yet, 00303 * so signify with negative number */ 00304 h2.nElecLevelOutput = -1; 00305 } 00306 else if( nMatch("GROU",chCard) ) 00307 { 00308 /* just the ground electronic state */ 00309 h2.nElecLevelOutput = 1; 00310 } 00311 else 00312 { 00313 h2.nElecLevelOutput = (int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00314 if( lgEOL ) 00315 { 00316 /* just the ground electronic state */ 00317 h2.nElecLevelOutput = 1; 00318 } 00319 } 00320 } 00321 } 00322 00323 else if( nMatch(" PDR",chCard) ) 00324 { 00325 /* creation and destruction processes */ 00326 strcpy( punch.chPunch[punch.npunch], "H2pd" ); 00327 sprintf( chHeader, "#H2 creation, destruction. \n" ); 00328 } 00329 else if( nMatch("POPU",chCard) ) 00330 { 00331 /* punch H2_populations */ 00332 strcpy( punch.chPunch[punch.npunch], "H2po" ); 00333 00334 /* this is an option to scan off highest vib and rot states 00335 * to punch pops - first is limit to vibration, then rotation 00336 * if no number is entered then 0 is set and all levels punched */ 00337 /* now get vib lim */ 00338 punch.punarg[punch.npunch][0] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00339 00340 /* this is limit to rotation quantum index */ 00341 punch.punarg[punch.npunch][1] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00342 00343 if( nMatch( "ZONE" , chCard ) ) 00344 { 00345 /* punch v=0 pops for each zone, all along one line */ 00346 punch.punarg[punch.npunch][2] = 0; 00347 sprintf( chHeader, "#depth\torth\tpar\te=1 rel pop\te=2 rel pop\tv=0 rel pops\n" ); 00348 } 00349 else 00350 { 00351 /* will not do zone output, only output at the end of the calculation 00352 * now check whether to punch triplets or a matrix for output - 00353 * default is triplets, so only check for matrix */ 00354 if( nMatch( "MATR" , chCard ) ) 00355 { 00356 /* matrix */ 00357 punch.punarg[punch.npunch][2] = 1; 00358 sprintf( chHeader, "#vib\trot\tpops\n" ); 00359 } 00360 else 00361 { 00362 /* triplets */ 00363 punch.punarg[punch.npunch][2] = -1; 00364 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" ); 00365 } 00366 } 00367 } 00368 00369 else if( nMatch("RATE",chCard) ) 00370 { 00371 /* punch h2 rates - creation and destruction rates */ 00372 strcpy( punch.chPunch[punch.npunch], "H2ra" ); 00373 sprintf( chHeader, 00374 "#depth\tN(H2)\tN(H2)/u(H2)\tA_V(star)\tn(Eval)" 00375 "\tH2/Htot\trenorm\tfrm grn\tfrmH-\tdstTH85\tBD96\tELWERT\tBigH2\telec->H2g\telec->H2s" 00376 "\tG(TH85)\tG(DB96)\tCR\tEleclife\tShield(BD96)\tShield(H2)\tBigh2/G0(spc)\ttot dest" 00377 "\tHeatH2Dish_TH85\tHeatH2Dexc_TH85\tHeatH2Dish_BigH2\tHeatH2Dexc_BigH2\thtot\n" ); 00378 } 00379 else if( nMatch("SOLO",chCard) ) 00380 { 00381 /* rate of Solomon process then fracs of exits from each v, J level */ 00382 strcpy( punch.chPunch[punch.npunch], "H2so" ); 00383 sprintf( chHeader, 00384 "#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" ); 00385 } 00386 else if( nMatch("SPEC",chCard) ) 00387 { 00388 /* special punch command*/ 00389 strcpy( punch.chPunch[punch.npunch], "H2sp" ); 00390 sprintf( chHeader, 00391 "#depth\tspecial\n" ); 00392 } 00393 else if( nMatch("TEMP",chCard) ) 00394 { 00395 /* various temperatures for neutral/molecular gas */ 00396 strcpy( punch.chPunch[punch.npunch], "H2te" ); 00397 sprintf( chHeader, 00398 "#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"); 00399 } 00400 else if( nMatch("THER",chCard) ) 00401 { 00402 /* thermal heating cooling processes involving H2 */ 00403 strcpy( punch.chPunch[punch.npunch], "H2th" ); 00404 sprintf( chHeader, 00405 "#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"); 00406 } 00407 else 00408 { 00409 fprintf( ioQQQ, 00410 " There must be a second key; they are RATE, LINE, COOL, COLUMN, _PDR, SOLOmon, TEMP, and POPUlations\n" ); 00411 cdEXIT(EXIT_FAILURE); 00412 } 00413 return; 00414 } 00415 00416 00417 /*H2_Prt_Zone print H2 info into zone results, called from prtzone for each printed zone */ 00418 void H2_Prt_Zone(void) 00419 { 00420 int iElecHi , iVibHi; 00421 00422 /* no print if H2 not turned on, or not computed for these conditions */ 00423 if( !h2.lgH2ON || !h2.nCallH2_this_zone ) 00424 return; 00425 00426 DEBUG_ENTRY( "H2_Prt_Zone()" ); 00427 00428 fprintf( ioQQQ, " H2 density "); 00429 fprintf(ioQQQ,PrintEfmt("%9.2e", hmi.H2_total)); 00430 00431 fprintf( ioQQQ, " orth/par"); 00432 fprintf(ioQQQ,PrintEfmt("%9.2e", h2.ortho_density / SDIV( h2.para_density ))); 00433 00434 iElecHi = 0; 00435 iVibHi = 0; 00436 fprintf( ioQQQ, " v0 J=0,3"); 00437 fprintf(ioQQQ,PrintEfmt("%9.2e", H2_populations[iElecHi][iVibHi][0] / hmi.H2_total)); 00438 fprintf(ioQQQ,PrintEfmt("%9.2e", H2_populations[iElecHi][iVibHi][1] / hmi.H2_total)); 00439 fprintf(ioQQQ,PrintEfmt("%9.2e", H2_populations[iElecHi][iVibHi][2] / hmi.H2_total)); 00440 fprintf(ioQQQ,PrintEfmt("%9.2e", H2_populations[iElecHi][iVibHi][3] / hmi.H2_total)); 00441 00442 fprintf( ioQQQ, " TOTv=0,3"); 00443 fprintf(ioQQQ,PrintEfmt("%9.2e", pops_per_vib[iElecHi][0] / hmi.H2_total)); 00444 fprintf(ioQQQ,PrintEfmt("%9.2e", pops_per_vib[iElecHi][1] / hmi.H2_total)); 00445 fprintf(ioQQQ,PrintEfmt("%9.2e", pops_per_vib[iElecHi][2] / hmi.H2_total)); 00446 fprintf(ioQQQ,PrintEfmt("%9.2e", pops_per_vib[iElecHi][3] / hmi.H2_total)); 00447 fprintf( ioQQQ, "\n"); 00448 return; 00449 } 00450 00451 /*H2_Prt_column_density print H2 info into zone results, called from prtzone for each printed zone */ 00452 void H2_Prt_column_density( 00453 /* this is stream used for io, is stdout when called by final, 00454 * is punch unit when punch output generated */ 00455 FILE *ioMEAN ) 00456 00457 { 00458 int iVibHi; 00459 00460 /* no print if H2 not turned on, or not computed for these conditions */ 00461 if( !h2.lgH2ON || !h2.nCallH2_this_zone ) 00462 return; 00463 00464 DEBUG_ENTRY( "H2_Prt_column_density()" ); 00465 00466 fprintf( ioMEAN, " H2 total "); 00467 fprintf(ioMEAN,"%7.3f", log10(SDIV(colden.colden[ipCOL_H2g]+colden.colden[ipCOL_H2s]))); 00468 00469 fprintf( ioMEAN, " H2 ortho "); 00470 fprintf(ioMEAN,"%7.3f", log10(SDIV(h2.ortho_colden))); 00471 00472 fprintf( ioMEAN, " para"); 00473 fprintf(ioMEAN,"%7.3f", log10(SDIV(h2.para_colden))); 00474 00475 iVibHi = 0; 00476 fprintf( ioMEAN, " v0 J=0,3"); 00477 fprintf(ioMEAN,"%7.3f", log10(SDIV(H2_X_colden[iVibHi][0]))); 00478 fprintf(ioMEAN,"%7.3f", log10(SDIV(H2_X_colden[iVibHi][1]))); 00479 fprintf(ioMEAN,"%7.3f", log10(SDIV(H2_X_colden[iVibHi][2]))); 00480 fprintf(ioMEAN,"%7.3f", log10(SDIV(H2_X_colden[iVibHi][3]))); 00481 00482 # if 0 00483 fprintf( ioMEAN, " v=0,3"); 00484 fprintf(ioMEAN,PrintEfmt("%9.2e", pops_per_vib[iElecHi][0] / hmi.H2_total)); 00485 fprintf(ioMEAN,PrintEfmt("%9.2e", pops_per_vib[iElecHi][1] / hmi.H2_total)); 00486 fprintf(ioMEAN,PrintEfmt("%9.2e", pops_per_vib[iElecHi][2] / hmi.H2_total)); 00487 fprintf(ioMEAN,PrintEfmt("%9.2e", pops_per_vib[iElecHi][3] / hmi.H2_total)); 00488 fprintf( ioMEAN, "\n"); 00489 # endif 00490 return; 00491 } 00492 00493 00494 /*H2_ReadTransprob read transition probabilities */ 00495 void H2_ReadTransprob( long int nelec ) 00496 { 00497 const char* cdDATAFILE[N_H2_ELEC] = 00498 { 00499 "H2_transprob_X.dat", 00500 "H2_transprob_B.dat", 00501 "H2_transprob_C_plus.dat", 00502 "H2_transprob_C_minus.dat", 00503 "H2_transprob_B_primed.dat", 00504 "H2_transprob_D_plus.dat", 00505 "H2_transprob_D_minus.dat" 00506 }; 00507 FILE *ioDATA; 00508 char chLine[FILENAME_PATH_LENGTH_2]; 00509 long int i, n1, n2, n3; 00510 long int iVibHi , iVibLo , iRotHi , iRotLo , iElecHi , iElecLo; 00511 bool lgEOL; 00512 00513 DEBUG_ENTRY( "H2_ReadTransprob()" ); 00514 00515 /* now open the data file */ 00516 ioDATA = open_data( cdDATAFILE[nelec], "r" ); 00517 00518 /* read the first line and check that magic number is ok */ 00519 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL ) 00520 { 00521 fprintf( ioQQQ, " H2_ReadTransprob could not read first line of %s\n", cdDATAFILE[nelec]); 00522 cdEXIT(EXIT_FAILURE); 00523 } 00524 i = 1; 00525 /* level 1 magic number */ 00526 n1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL); 00527 n2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL); 00528 n3 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL); 00529 00530 /* magic number 00531 * the following is the set of numbers that appear at the start of level1.dat 01 08 10 */ 00532 if( ( n1 != 2 ) || ( n2 != 4 ) || ( n3 != 29 ) ) 00533 { 00534 fprintf( ioQQQ, 00535 " H2_ReadTransprob: the version of %s is not the current version.\n", cdDATAFILE[nelec] ); 00536 fprintf( ioQQQ, 00537 " I expected to find the number 2 4 29 and got %li %li %li instead.\n" , 00538 n1 , n2 , n3 ); 00539 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine ); 00540 cdEXIT(EXIT_FAILURE); 00541 } 00542 00543 /* read until not a comment */ 00544 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL ) 00545 BadRead(); 00546 00547 while( chLine[0]=='#' ) 00548 { 00549 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL ) 00550 BadRead(); 00551 } 00552 iVibHi = 1; 00553 while( iVibHi >= 0 ) 00554 { 00555 double Aul; 00556 sscanf(chLine,"%li\t%li\t%li\t%li\t%li\t%li\t%le", 00557 &iElecHi , &iVibHi ,&iRotHi , &iElecLo , &iVibLo , &iRotLo , &Aul ); 00558 ASSERT( iElecHi == nelec ); 00559 /* negative iVibHi says end of data */ 00560 if( iVibHi < 0 ) 00561 continue; 00562 00563 ASSERT( iElecHi < N_H2_ELEC ); 00564 ASSERT( iElecLo < N_H2_ELEC ); 00566 ASSERT( iVibHi < 50 ); 00567 ASSERT( iVibLo < 50 ); 00568 00569 /* check that we actually included the levels in the model representation */ 00570 if( iVibHi <= h2.nVib_hi[iElecHi] && 00571 iVibLo <= h2.nVib_hi[iElecLo] && 00572 iRotHi <= h2.nRot_hi[iElecHi][iVibHi] && 00573 iRotLo <= h2.nRot_hi[iElecLo][iVibLo]) 00574 { 00575 double ener = energy_wn[iElecHi][iVibHi][iRotHi] - energy_wn[iElecLo][iVibLo][iRotLo]; 00576 00577 /* only lines that have real Aul are added to stack. */ 00578 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis = AddLine2Stack( true ); 00579 00580 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->Aul = (realnum)Aul; 00581 /* say that this line exists */ 00582 lgH2_line_exists[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] = true; 00583 /* prints transitions with negative energies - should not happen */ 00584 if( ener <= 0. ) 00585 { 00586 fprintf(ioQQQ,"negative energy H2 transition\t%li\t%li\t%li\t%li\t%.2e\t%.2e\n", 00587 iVibHi,iVibLo,iRotHi,iRotLo,Aul, 00588 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyWN); 00589 ShowMe(); 00590 } 00591 } 00592 # if 0 00593 /* this prints all levels with As but without energies */ 00594 else 00595 { 00596 fprintf(ioQQQ,"no\t%li\t%li\t%li\t%li\t%.2e\n", 00597 iVibHi,iVibLo,iRotHi,iRotLo,Aul); 00598 } 00599 # endif 00600 00601 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL ) 00602 BadRead(); 00603 while( chLine[0]=='#' ) 00604 { 00605 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL ) 00606 BadRead(); 00607 } 00608 } 00609 fclose( ioDATA ); 00610 return; 00611 } 00612 00613 #if 0 00614 /*H2_Read_Cosmicray_distribution read distribution function for H2 population following cosmic ray collisional excitation */ 00615 void H2_Read_Cosmicray_distribution(void) 00616 { 00617 /*>>refer H2 cr excit Tine, S., Lepp, S., Gredel, R., & Dalgarno, A. 1997, ApJ, 481, 282 */ 00618 FILE *ioDATA; 00619 char chLine[FILENAME_PATH_LENGTH_2]; 00620 long int i, n1, n2, n3, iVib , iRot; 00621 long neut_frac; 00622 bool lgEOL; 00623 00624 DEBUG_ENTRY( "H2_Read_Cosmicray_distribution()" ); 00625 00626 /* now open the data file */ 00627 ioDATA = open_data( "H2_CosmicRay_collision.dat", "r" ); 00628 00629 /* read the first line and check that magic number is ok */ 00630 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL ) 00631 { 00632 fprintf( ioQQQ, " H2_Read_Cosmicray_distribution could not read first line of %s\n", "H2_Cosmic_collision.dat"); 00633 cdEXIT(EXIT_FAILURE); 00634 } 00635 00636 i = 1; 00637 /* level 1 magic number */ 00638 n1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL); 00639 n2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL); 00640 n3 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL); 00641 00642 /* magic number 00643 * the following is the set of numbers that appear at the start of H2_Cosmic_collision.dat 01 21 03 */ 00644 if( ( n1 != 1 ) || ( n2 != 21 ) || ( n3 != 3 ) ) 00645 { 00646 fprintf( ioQQQ, 00647 " H2_Read_Cosmicray_distribution: the version of %s is not the current version.\n", "H2_Cosmic_collision.dat" ); 00648 fprintf( ioQQQ, 00649 " I expected to find the number 1 21 3 and got %li %li %li instead.\n" , 00650 n1 , n2 , n3 ); 00651 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine ); 00652 cdEXIT(EXIT_FAILURE); 00653 } 00654 00655 /* read until not a comment */ 00656 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL ) 00657 BadRead(); 00658 00659 while( chLine[0]=='#' ) 00660 { 00661 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL ) 00662 BadRead(); 00663 } 00664 00665 iRot = 1; 00666 iVib = 1; 00667 neut_frac = 0; 00668 while( iVib >= 0 ) 00669 { 00670 long int j_minus_ji; 00671 double a[10]; 00672 00673 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", 00674 &iVib ,&j_minus_ji , &a[0],&a[1],&a[2],&a[3],&a[4],&a[5],&a[6],&a[7],&a[8],&a[9] 00675 ); 00676 /* negative iVib says end of data */ 00677 if( iVib < 0 ) 00678 continue; 00679 00680 /* cr_rate[CR_X][CR_VIB][CR_J][CR_EXIT];*/ 00681 /* check that we actually included the levels in the model representation */ 00682 ASSERT( iVib < CR_VIB ); 00683 ASSERT( j_minus_ji == -2 || j_minus_ji == +2 || j_minus_ji == 0 ); 00684 ASSERT( neut_frac < CR_X ); 00685 00686 /* now make i_minus_ji an array index */ 00687 j_minus_ji = 1 + j_minus_ji/2; 00688 ASSERT( j_minus_ji>=0 && j_minus_ji<=2 ); 00689 00690 /* option to add Gaussian random mole */ 00691 for( iRot=0; iRot<CR_J; ++iRot ) 00692 { 00693 cr_rate[neut_frac][iVib][iRot][j_minus_ji] = (realnum)a[iRot]; 00694 } 00695 if( mole.lgH2_NOISECOSMIC ) 00696 { 00697 realnum r; 00698 r = (realnum)RandGauss( mole.xMeanNoise , mole.xSTDNoise ); 00699 00700 for( iRot=0; iRot<CR_J; ++iRot ) 00701 { 00702 cr_rate[neut_frac][iVib][iRot][j_minus_ji] *= (realnum)pow(10.,(double)r); 00703 } 00704 } 00705 00706 if( CR_PRINT ) 00707 { 00708 fprintf(ioQQQ,"cr rate\t%li\t%li", iVib , j_minus_ji ); 00709 for( iRot=0; iRot<CR_J; ++iRot ) 00710 { 00711 fprintf(ioQQQ,"\t%.3e", cr_rate[neut_frac][iVib][iRot][j_minus_ji] ); 00712 } 00713 fprintf(ioQQQ,"\n" ); 00714 } 00715 00716 /* now get next line */ 00717 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL ) 00718 BadRead(); 00719 while( chLine[0]=='#' ) 00720 { 00721 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL ) 00722 BadRead(); 00723 } 00724 } 00725 fclose( ioDATA ); 00726 00727 return; 00728 } 00729 #endif 00730 00731 /*H2_ReadEnergies read energies for all electronic levels */ 00732 void H2_ReadEnergies( long int nelec ) 00733 { 00734 const char* cdDATAFILE[N_H2_ELEC] = 00735 { 00736 "H2_energy_X.dat", 00737 "H2_energy_B.dat", 00738 "H2_energy_C_plus.dat", 00739 "H2_energy_C_minus.dat", 00740 "H2_energy_B_primed.dat", 00741 "H2_energy_D_plus.dat", 00742 "H2_energy_D_minus.dat" 00743 }; 00744 FILE *ioDATA; 00745 char chLine[FILENAME_PATH_LENGTH_2]; 00746 long int i, n1, n2, n3, iVib , iRot; 00747 bool lgEOL; 00748 00749 DEBUG_ENTRY( "H2_ReadEnergies()" ); 00750 00751 /* now open the data file */ 00752 ioDATA = open_data( cdDATAFILE[nelec], "r" ); 00753 00754 /* read the first line and check that magic number is ok */ 00755 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL ) 00756 { 00757 fprintf( ioQQQ, " H2_ReadEnergies could not read first line of %s\n", cdDATAFILE[nelec]); 00758 cdEXIT(EXIT_FAILURE); 00759 } 00760 i = 1; 00761 /* level 1 magic number */ 00762 n1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL); 00763 n2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL); 00764 n3 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL); 00765 00766 /* magic number 00767 * the following is the set of numbers that appear at the start of level1.dat 01 08 10 */ 00768 if( ( n1 != 2 ) || ( n2 != 4 ) || ( n3 != 29 ) ) 00769 { 00770 fprintf( ioQQQ, 00771 " H2_ReadEnergies: the version of %s is not the current version.\n", cdDATAFILE[nelec] ); 00772 fprintf( ioQQQ, 00773 " I expected to find the number 2 4 29 and got %li %li %li instead.\n" , 00774 n1 , n2 , n3 ); 00775 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine ); 00776 cdEXIT(EXIT_FAILURE); 00777 } 00778 00779 /* read until not a comment */ 00780 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL ) 00781 BadRead(); 00782 00783 while( chLine[0]=='#' ) 00784 { 00785 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL ) 00786 BadRead(); 00787 } 00788 00789 /* this will count the number of levels within each electronic state */ 00790 nLevels_per_elec[nelec] = 0; 00791 00792 for( iVib=0; iVib<=h2.nVib_hi[nelec]; ++iVib ) 00793 { 00794 for( iRot=h2.Jlowest[nelec]; iRot<=h2.nRot_hi[nelec][iVib]; ++iRot ) 00795 { 00796 i = 1; 00797 sscanf(chLine,"%li\t%li\t%le", &n1 , &n2 , &energy_wn[nelec][iVib][iRot] ); 00798 ASSERT( n1 == iVib ); 00799 ASSERT( n2 == iRot ); 00800 # if 0 00801 /* in atomic units, or 1 Hartree, or two rydbergs */ 00802 if( nelec == 0 ) 00803 { 00804 /* only do this for Phillip Stancil's file */ 00805 /* corrections are to get lowest rotation level to have energy of zero */ 00806 energy_wn[0][iVib][iRot] = -( energy_wn[0][iVib][iRot]- 3.6118114E+04 ); 00807 } 00808 # endif 00809 ASSERT( energy_wn[nelec][iVib][iRot]> 0. || (nelec==0 && iVib==0 && iRot==0 ) ); 00810 /* increment number of levels within this electronic state */ 00811 ++nLevels_per_elec[nelec]; 00812 00813 /* now start reading next line */ 00814 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL ) 00815 BadRead(); 00816 while( chLine[0]=='#' ) 00817 { 00818 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL ) 00819 BadRead(); 00820 } 00821 } 00822 } 00823 fclose( ioDATA ); 00824 return; 00825 } 00826 00827 /*H2_ReadDissprob read dissociation probabilities and kinetic energies for all electronic levels */ 00828 void H2_ReadDissprob( long int nelec ) 00829 { 00830 const char* cdDATAFILE[N_H2_ELEC] = 00831 { 00832 "H2_dissprob_X.dat",/* this does not exist and nelec == 0 is not valid */ 00833 "H2_dissprob_B.dat", 00834 "H2_dissprob_C_plus.dat", 00835 "H2_dissprob_C_minus.dat", 00836 "H2_dissprob_B_primed.dat", 00837 "H2_dissprob_D_plus.dat", 00838 "H2_dissprob_D_minus.dat" 00839 }; 00840 FILE *ioDATA; 00841 char chLine[FILENAME_PATH_LENGTH_2]; 00842 long int i, n1, n2, n3, iVib , iRot; 00843 bool lgEOL; 00844 00845 DEBUG_ENTRY( "H2_ReadDissprob()" ); 00846 00847 ASSERT( nelec > 0 ); 00848 00849 /* now open the data file */ 00850 ioDATA = open_data( cdDATAFILE[nelec], "r" ); 00851 00852 /* read the first line and check that magic number is ok */ 00853 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL ) 00854 { 00855 fprintf( ioQQQ, " H2_ReadDissprob could not read first line of %s\n", cdDATAFILE[nelec]); 00856 cdEXIT(EXIT_FAILURE); 00857 } 00858 i = 1; 00859 /* level 1 magic number */ 00860 n1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL); 00861 n2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL); 00862 n3 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL); 00863 00864 /* magic number 00865 * the following is the set of numbers that appear at the start of level1.dat 01 08 10 */ 00866 if( ( n1 != 3 ) || ( n2 != 2 ) || ( n3 != 11 ) ) 00867 { 00868 fprintf( ioQQQ, 00869 " H2_ReadDissprob: the version of %s is not the current version.\n", cdDATAFILE[nelec] ); 00870 fprintf( ioQQQ, 00871 " I expected to find the number 3 2 11 and got %li %li %li instead.\n" , 00872 n1 , n2 , n3 ); 00873 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine ); 00874 cdEXIT(EXIT_FAILURE); 00875 } 00876 00877 /* read until not a comment */ 00878 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL ) 00879 BadRead(); 00880 00881 while( chLine[0]=='#' ) 00882 { 00883 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL ) 00884 BadRead(); 00885 } 00886 00887 for( iVib=0; iVib<=h2.nVib_hi[nelec]; ++iVib ) 00888 { 00889 for( iRot=h2.Jlowest[nelec]; iRot<=h2.nRot_hi[nelec][iVib]; ++iRot ) 00890 { 00891 double a, b; 00892 i = 1; 00893 sscanf(chLine,"%li\t%li\t%le\t%le", 00894 &n1 , &n2 , 00895 /* dissociation probability */ 00896 &a , 00897 /* dissociation kinetic energy - eV not ergs */ 00898 &b); 00899 00900 /* these have to agree if data file is valid */ 00901 ASSERT( n1 == iVib ); 00902 ASSERT( n2 == iRot ); 00903 00904 /* dissociation probability */ 00905 H2_dissprob[nelec][iVib][iRot] = (realnum)a; 00906 /* dissociation kinetic energy - eV not ergs */ 00907 H2_disske[nelec][iVib][iRot] = (realnum)b; 00908 00909 /* now get next line */ 00910 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL ) 00911 BadRead(); 00912 while( chLine[0]=='#' ) 00913 { 00914 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL ) 00915 BadRead(); 00916 } 00917 } 00918 } 00919 fclose( ioDATA ); 00920 return; 00921 } 00922 00923 00924 /*H2_Read_hminus_distribution read distribution function for H2 population following formation from H minus */ 00925 void H2_Read_hminus_distribution(void) 00926 { 00927 FILE *ioDATA; 00928 char chLine[FILENAME_PATH_LENGTH_2]; 00929 long int i, n1, n2, n3, iVib , iRot; 00930 bool lgEOL; 00931 double sumrate[nTE_HMINUS]; 00932 /* set true for lots of printout */ 00933 # define H2HMINUS_PRT false 00934 00935 DEBUG_ENTRY( "H2_Read_hminus_distribution()" ); 00936 00937 /* now open the data file */ 00938 ioDATA = open_data( "H2_hminus_deposit.dat", "r" ); 00939 00940 /* read the first line and check that magic number is ok */ 00941 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL ) 00942 { 00943 fprintf( ioQQQ, " H2_Read_hminus_distribution could not read first line of %s\n", "H2_hminus_deposit.dat"); 00944 cdEXIT(EXIT_FAILURE); 00945 } 00946 00947 i = 1; 00948 /* level 1 magic number */ 00949 n1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL); 00950 n2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL); 00951 n3 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL); 00952 00953 /* magic number 00954 * the following is the set of numbers that appear at the start of H2_hminus_deposit.dat 01 08 10 */ 00955 if( ( n1 != 2 ) || ( n2 != 10 ) || ( n3 != 17 ) ) 00956 { 00957 fprintf( ioQQQ, 00958 " H2_Read_hminus_distribution: the version of %s is not the current version.\n", "H2_hminus_deposit.dat" ); 00959 fprintf( ioQQQ, 00960 " I expected to find the number 2 10 17 and got %li %li %li instead.\n" , 00961 n1 , n2 , n3 ); 00962 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine ); 00963 cdEXIT(EXIT_FAILURE); 00964 } 00965 00966 /* read until not a comment */ 00967 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL ) 00968 BadRead(); 00969 00970 while( chLine[0]=='#' ) 00971 { 00972 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL ) 00973 BadRead(); 00974 } 00975 00976 /* convert temps to log */ 00977 for(i=0; i<nTE_HMINUS; ++i ) 00978 { 00979 H2_te_hminus[i] = (realnum)log10(H2_te_hminus[i]); 00980 sumrate[i] = 0.; 00981 } 00982 00983 iRot = 1; 00984 iVib = 1; 00985 while( iVib >= 0 ) 00986 { 00987 /* set true to print rates */ 00988 00989 double a[nTE_HMINUS] , ener; 00990 sscanf(chLine,"%li\t%li\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf", 00991 &iVib ,&iRot , &ener, &a[0],&a[1],&a[2] , &a[3],&a[4],&a[5] ,&a[6] 00992 ); 00993 /* negative iVib says end of data */ 00994 if( iVib < 0 ) 00995 continue; 00996 00997 /* check that we actually included the levels in the model representation */ 00998 ASSERT( iVib <= h2.nVib_hi[0] && 00999 iRot <= h2.nRot_hi[0][iVib] ); 01000 01001 if( H2HMINUS_PRT ) 01002 fprintf(ioQQQ,"hminusss\t%li\t%li", iVib , iRot ); 01003 for( i=0; i<nTE_HMINUS; ++i ) 01004 { 01005 H2_X_hminus_formation_distribution[i][iVib][iRot] = (realnum)pow(10.,-a[i]); 01006 sumrate[i] += H2_X_hminus_formation_distribution[i][iVib][iRot]; 01007 if( H2HMINUS_PRT ) 01008 fprintf(ioQQQ,"\t%.3e", H2_X_hminus_formation_distribution[i][iVib][iRot] ); 01009 } 01010 if( H2HMINUS_PRT ) 01011 fprintf(ioQQQ,"\n" ); 01012 01013 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL ) 01014 BadRead(); 01015 while( chLine[0]=='#' ) 01016 { 01017 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL ) 01018 BadRead(); 01019 } 01020 } 01021 fclose( ioDATA ); 01022 01023 if( H2HMINUS_PRT ) 01024 { 01025 /* print total rate */ 01026 fprintf(ioQQQ," total H- formation rate "); 01027 /* convert temps to log */ 01028 for(i=0; i<nTE_HMINUS; ++i ) 01029 { 01030 fprintf(ioQQQ,"\t%.3e" , sumrate[i]); 01031 } 01032 fprintf(ioQQQ,"\n" ); 01033 } 01034 01035 /* convert to dimensionless factors that add to unity */ 01036 for( iVib=0; iVib<=h2.nVib_hi[0]; ++iVib ) 01037 { 01038 for( iRot=h2.Jlowest[0]; iRot<=h2.nRot_hi[0][iVib]; ++iRot ) 01039 { 01040 for(i=0; i<nTE_HMINUS; ++i ) 01041 { 01042 H2_X_hminus_formation_distribution[i][iVib][iRot] /= (realnum)sumrate[i]; 01043 } 01044 } 01045 } 01046 01047 if( H2HMINUS_PRT ) 01048 { 01049 /* print total rate */ 01050 fprintf(ioQQQ," H- distribution function "); 01051 for( iVib=0; iVib<=h2.nVib_hi[0]; ++iVib ) 01052 { 01053 for( iRot=h2.Jlowest[0]; iRot<=h2.nRot_hi[0][iVib]; ++iRot ) 01054 { 01055 fprintf(ioQQQ,"%li\t%li", iVib , iRot ); 01056 for(i=0; i<nTE_HMINUS; ++i ) 01057 { 01058 fprintf(ioQQQ,"\t%.3e", H2_X_hminus_formation_distribution[i][iVib][iRot] ); 01059 } 01060 fprintf(ioQQQ,"\n" ); 01061 } 01062 } 01063 } 01064 return; 01065 } 01066 01067 /* ===================================================================== */ 01068 /*H2_Punch_line_data punch line data for H2 molecule */ 01069 void H2_Punch_line_data( 01070 /* io unit for punch */ 01071 FILE* ioPUN , 01072 /* punch all levels if true, only subset if false */ 01073 bool lgDoAll ) 01074 { 01075 long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo; 01076 01077 if( !h2.lgH2ON ) 01078 return; 01079 01080 DEBUG_ENTRY( "H2_Punch_line_data()" ); 01081 01082 if( lgDoAll ) 01083 { 01084 fprintf( ioQQQ, 01085 " H2_Punch_line_data ALL option not implemented in H2_Punch_line_data yet 1\n" ); 01086 cdEXIT(EXIT_FAILURE); 01087 } 01088 else 01089 { 01090 01091 /* punch line date, looping over all possible lines */ 01092 for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi ) 01093 { 01094 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi ) 01095 { 01096 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi ) 01097 { 01098 /* now the lower levels */ 01099 /* NB - X is the only lower level considered here, since we are only 01100 * concerned with excited electronic levels as a photodissociation process 01101 * code exists to relax this assumption - simply change following to iElecHi */ 01102 long int lim_elec_lo = 0; 01103 for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo ) 01104 { 01105 /* want to include all vib states in lower level if different electronic level, 01106 * but only lower vib levels if same electronic level */ 01107 long int nv = h2.nVib_hi[iElecLo]; 01108 if( iElecLo==iElecHi ) 01109 nv = iVibHi; 01110 for( iVibLo=0; iVibLo<=nv; ++iVibLo ) 01111 { 01112 long nr = h2.nRot_hi[iElecLo][iVibLo]; 01113 if( iElecLo==iElecHi && iVibHi==iVibLo ) 01114 nr = iRotHi-1; 01115 01116 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo ) 01117 { 01118 /* only punch if radiative transition exists */ 01119 if( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].ipCont > 0 ) 01120 { 01122 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Coll.cs = 0.; 01123 /* print quantum indices */ 01124 fprintf(ioPUN,"%2li %2li %2li %2li %2li %2li ", 01125 iElecHi,iVibHi,iRotHi,iElecLo,iVibLo,iRotLo ); 01126 Punch1LineData( &H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] , ioPUN , false); 01127 } 01128 } 01129 } 01130 } 01131 } 01132 } 01133 } 01134 fprintf( ioPUN , "\n"); 01135 } 01136 return; 01137 } 01138 01139 /*H2_PunchLineStuff include H2 lines in punched optical depths, etc, called from PunchLineStuff */ 01140 void H2_PunchLineStuff( FILE * io , realnum xLimit , long index) 01141 { 01142 long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo; 01143 01144 if( !h2.lgH2ON ) 01145 return; 01146 01147 DEBUG_ENTRY( "H2_PunchLineStuff()" ); 01148 01149 /* loop over all possible lines */ 01150 for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi ) 01151 { 01152 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi ) 01153 { 01154 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi ) 01155 { 01156 /* now the lower levels */ 01157 /* NB - X is the only lower level considered here, since we are only 01158 * concerned with excited electronic levels as a photodissociation process 01159 * code exists to relax this assumption - simply change following to iElecHi */ 01160 long int lim_elec_lo = 0; 01161 for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo ) 01162 { 01163 /* want to include all vib states in lower level if different electronic level, 01164 * but only lower vib levels if same electronic level */ 01165 long int nv = h2.nVib_hi[iElecLo]; 01166 if( iElecLo==iElecHi ) 01167 nv = iVibHi; 01168 for( iVibLo=0; iVibLo<=nv; ++iVibLo ) 01169 { 01170 long nr = h2.nRot_hi[iElecLo][iVibLo]; 01171 if( iElecLo==iElecHi && iVibHi==iVibLo ) 01172 nr = iRotHi-1; 01173 01174 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo ) 01175 { 01176 /* only punch if radiative transition exists */ 01177 if( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].ipCont > 0 ) 01178 { 01179 pun1Line( &H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] , io , xLimit , index , 1.); 01180 } 01181 } 01182 } 01183 } 01184 } 01185 } 01186 } 01187 01188 return; 01189 } 01190 01191 01192 /*H2_Prt_line_tau print line optical depths, called from premet in response to print line optical depths command*/ 01193 void H2_Prt_line_tau(void) 01194 { 01195 long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo; 01196 01197 if( !h2.lgH2ON ) 01198 return; 01199 01200 DEBUG_ENTRY( "H2_Prt_line_tau()" ); 01201 01202 /* loop over all possible lines */ 01203 for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi ) 01204 { 01205 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi ) 01206 { 01207 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi ) 01208 { 01209 /* now the lower levels */ 01210 /* NB - X is the only lower level considered here, since we are only 01211 * concerned with excited electronic levels as a photodissociation process 01212 * code exists to relax this assumption - simply change following to iElecHi */ 01213 long int lim_elec_lo = 0; 01214 for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo ) 01215 { 01216 /* want to include all vib states in lower level if different electronic level, 01217 * but only lower vib levels if same electronic level */ 01218 long int nv = h2.nVib_hi[iElecLo]; 01219 if( iElecLo==iElecHi ) 01220 nv = iVibHi; 01221 for( iVibLo=0; iVibLo<=nv; ++iVibLo ) 01222 { 01223 long nr = h2.nRot_hi[iElecLo][iVibLo]; 01224 if( iElecLo==iElecHi && iVibHi==iVibLo ) 01225 nr = iRotHi-1; 01226 01227 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo ) 01228 { 01229 /* only print if radiative transition exists */ 01230 if( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].ipCont > 0 ) 01231 { 01232 prme(" c",&H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] ); 01233 } 01234 } 01235 } 01236 } 01237 } 01238 } 01239 } 01240 01241 return; 01242 } 01243 01244 01245 /*chMolBranch returns a char with the spectroscopic branch of a transition */ 01246 STATIC char chMolBranch( long iRotHi , long int iRotLo ) 01247 { 01248 /* these are the spectroscopic branches */ 01249 char chBranch[5] = {'O','P','Q','R','S'}; 01250 /* this is the index within the chBranch array */ 01251 int ip = 2 + (iRotHi - iRotLo); 01252 if( ip<0 || ip>=5 ) 01253 { 01254 fprintf(ioQQQ," chMolBranch called with insane iRotHi=%li iRotLo=%li ip=%i\n", 01255 iRotHi , iRotLo , ip ); 01256 ip = 0; 01257 } 01258 01259 return( chBranch[ip] ); 01260 } 01261 01262 /*H2_PunchDo punch some properties of the large H2 molecule */ 01263 void H2_PunchDo( FILE* io , char chJOB[] , const char chTime[] , long int ipPun ) 01264 { 01265 long int iVibHi , iElecHi , iRotHi , iVibLo , iElecLo , iRotLo, 01266 ip; 01267 long int iRot , iVib; 01268 long int LimVib , LimRot; 01269 01270 DEBUG_ENTRY( "H2_PunchDo()" ); 01271 01272 /* which job are we supposed to do? This routine is active even when H2 is not turned on 01273 * so do not test on h2.lgH2ON initially */ 01274 01275 /* H2 populations computed in last zone - 01276 * give all of molecule in either matrix or triplet format */ 01277 if( (strcmp( chJOB , "H2po" ) == 0) && (strcmp(chTime,"LAST") == 0) && 01278 (punch.punarg[ipPun][2] != 0) ) 01279 { 01280 /* >>chng 04 feb 19, do not punch if H2 not yet evaluated */ 01281 if( h2.lgH2ON && hmi.lgBigH2_evaluated ) 01282 { 01283 iVibHi= 0; 01284 iRotHi = 0; 01285 iElecHi=0; 01286 /* the limit to the number of vibration levels punched - 01287 * default is all, but first two numbers on punch h2 pops command 01288 * reset limit */ 01289 /* this is limit to vibration */ 01290 if( punch.punarg[ipPun][0] > 0 ) 01291 { 01292 LimVib = (long)punch.punarg[ipPun][0]; 01293 } 01294 else 01295 { 01296 LimVib = h2.nVib_hi[iElecHi]; 01297 } 01298 01299 /* first punch the current ortho, para, and total H2 density */ 01300 fprintf(io,"%i\t%i\t%.3e\tortho\n", 01301 103 , 01302 103 , 01303 h2.ortho_density ); 01304 fprintf(io,"%i\t%i\t%.3e\tpara\n", 01305 101 , 01306 101 , 01307 h2.para_density ); 01308 fprintf(io,"%i\t%i\t%.3e\ttotal\n", 01309 0 , 01310 0 , 01311 hmi.H2_total ); 01312 01313 /* now punch the actual H2_populations, first part both matrix and triplets */ 01314 for( iVibHi=0; iVibHi<=LimVib; ++iVibHi ) 01315 { 01316 /* this is limit to rotation quantum index */ 01317 if( punch.punarg[ipPun][1] > 0 ) 01318 { 01319 LimRot = (long)punch.punarg[ipPun][1]; 01320 } 01321 else 01322 { 01323 LimRot = h2.nRot_hi[iElecHi][iVibHi]; 01324 } 01325 if( punch.punarg[ipPun][2] > 0 ) 01326 { 01327 long int i; 01328 /* this option punch matrix */ 01329 if( iVibHi == 0 ) 01330 { 01331 fprintf(io,"vib\\rot"); 01332 /* this is first vib, so make row of rot numbs */ 01333 for( i=0; i<=LimRot; ++i ) 01334 { 01335 fprintf(io,"\t%li",i); 01336 } 01337 fprintf(io,"\n"); 01338 } 01339 fprintf(io,"%li",iVibHi ); 01340 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=LimRot; ++iRotHi ) 01341 { 01342 fprintf(io,"\t%.3e", 01343 H2_populations[iElecHi][iVibHi][iRotHi]/hmi.H2_total ); 01344 } 01345 fprintf(io,"\n" ); 01346 } 01347 else if( punch.punarg[ipPun][2] < 0 ) 01348 { 01349 /* this option punch triplets - the default */ 01350 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=LimRot; ++iRotHi ) 01351 { 01352 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", 01353 /* upper vibration and rotation quantum numbers */ 01354 iVibHi , iRotHi , 01355 /* an 'O' or 'P' for ortho or para */ 01356 chlgPara[H2_lgOrtho[iElecHi][iVibHi][iRotHi]], 01357 /* the level excitation energy in wavenumbers */ 01358 energy_wn[iElecHi][iVibHi][iRotHi], 01359 /* actual population relative to total H2 */ 01360 H2_populations[iElecHi][iVibHi][iRotHi]/hmi.H2_total , 01361 /* old level H2_populations for comparison */ 01362 H2_old_populations[iElecHi][iVibHi][iRotHi]/hmi.H2_total , 01363 /* H2_populations per h2 and per statistical weight */ 01364 H2_populations[iElecHi][iVibHi][iRotHi]/hmi.H2_total/H2_stat[iElecHi][iVibHi][iRotHi] , 01365 /* LTE departure coefficient */ 01366 /* >>chng 05 jan 26, missing factor of H2 abundance LTE is norm to unity, not tot abund */ 01367 H2_populations[iElecHi][iVibHi][iRotHi]/SDIV(H2_populations_LTE[iElecHi][iVibHi][iRotHi]*hmi.H2_total ) , 01368 /* fraction of exits that were collisional */ 01369 H2_col_rate_out[iVibHi][iRotHi]/SDIV(H2_col_rate_out[iVibHi][iRotHi]+H2_rad_rate_out[0][iVibHi][iRotHi]) , 01370 /* fraction of entries that were collisional */ 01371 H2_col_rate_in[iVibHi][iRotHi]/SDIV(H2_col_rate_in[iVibHi][iRotHi]+H2_rad_rate_in[iVibHi][iRotHi]), 01372 /* collisions out */ 01373 H2_col_rate_out[iVibHi][iRotHi], 01374 /* radiation out */ 01375 H2_rad_rate_out[0][iVibHi][iRotHi] , 01376 /* radiation out */ 01377 H2_col_rate_in[iVibHi][iRotHi], 01378 /* radiation in */ 01379 H2_rad_rate_in[iVibHi][iRotHi] 01380 ); 01381 } 01382 } 01383 } 01384 } 01385 } 01386 /* punch H2 populations for each zone 01387 * H2_populations of v=0 for each zone */ 01388 else if( (strcmp( chJOB , "H2po" ) == 0) && (strcmp(chTime,"LAST") != 0) && 01389 (punch.punarg[ipPun][2] == 0) ) 01390 { 01391 /* >>chng 04 feb 19, do not punch if h2 not yet evaluated */ 01392 if( h2.lgH2ON && hmi.lgBigH2_evaluated ) 01393 { 01394 fprintf(io,"%.5e\t%.3e\t%.3e", radius.depth_mid_zone , 01395 h2.ortho_density , h2.para_density); 01396 /* rel pops of first two excited electronic states */ 01397 fprintf(io,"\t%.3e\t%.3e", 01398 pops_per_elec[1] , pops_per_elec[2]); 01399 iElecHi = 0; 01400 iVibHi = 0; 01401 /* this is limit to vibration quantum index */ 01402 if( punch.punarg[ipPun][0] > 0 ) 01403 { 01404 LimVib = (long)punch.punarg[ipPun][1]; 01405 } 01406 else 01407 { 01408 LimVib = h2.nRot_hi[iElecHi][iVibHi]; 01409 } 01410 LimVib = MIN2( LimVib , h2.nVib_hi[iElecHi] ); 01411 /* this is limit to rotation quantum index */ 01412 if( punch.punarg[ipPun][1] > 0 ) 01413 { 01414 LimRot = (long)punch.punarg[ipPun][1]; 01415 } 01416 else 01417 { 01418 LimRot = h2.nRot_hi[iElecHi][iVibHi]; 01419 } 01420 for( iVibHi = 0; iVibHi<=LimVib; ++iVibHi ) 01421 { 01422 long int LimRotVib = MIN2( LimRot , h2.nRot_hi[iElecHi][iVibHi] ); 01423 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=LimRotVib; ++iRotHi ) 01424 { 01425 fprintf(io,"\t%.3e", 01426 H2_populations[iElecHi][iVibHi][iRotHi]/hmi.H2_total ); 01427 } 01428 fprintf(io,"\t"); 01429 } 01430 fprintf(io,"\n"); 01431 } 01432 } 01433 01434 /* punch column densities */ 01435 else if( (strcmp( chJOB , "H2cl" ) == 0) && (strcmp(chTime,"LAST") == 0) ) 01436 { 01437 iVibHi= 0; 01438 iRotHi = 0; 01439 iElecHi=0; 01440 /* the limit to the number of vibration levels punched - 01441 * default is all, but first two numbers on punch h2 pops command 01442 * reset limit */ 01443 /* this is limit to vibration */ 01444 if( punch.punarg[ipPun][0] > 0 ) 01445 { 01446 LimVib = (long)punch.punarg[ipPun][0]; 01447 } 01448 else 01449 { 01450 LimVib = h2.nVib_hi[iElecHi]; 01451 } 01452 01453 /* first punch ortho and para H2_populations */ 01454 fprintf(io,"%i\t%i\t%.3e\tortho\n", 01455 103 , 01456 103 , 01457 h2.ortho_colden ); 01458 fprintf(io,"%i\t%i\t%.3e\tpara\n", 01459 101 , 01460 101 , 01461 h2.para_colden ); 01462 /* total H2 column density */ 01463 fprintf(io,"%i\t%i\t%.3e\ttotal\n", 01464 0 , 01465 0 , 01466 colden.colden[ipCOL_H2g]+colden.colden[ipCOL_H2s]); 01467 01468 /* punch level column densities */ 01469 for( iVibHi=0; iVibHi<=LimVib; ++iVibHi ) 01470 { 01471 if( h2.lgH2ON ) 01472 { 01473 /* this is limit to rotation quantum index */ 01474 if( punch.punarg[ipPun][1] > 0 ) 01475 { 01476 LimRot = (long)punch.punarg[ipPun][1]; 01477 } 01478 else 01479 { 01480 LimRot = h2.nRot_hi[iElecHi][iVibHi]; 01481 } 01482 if( punch.punarg[ipPun][2] > 0 ) 01483 { 01484 long int i; 01485 /* punch matrix */ 01486 if( iVibHi == 0 ) 01487 { 01488 fprintf(io,"vib\\rot"); 01489 /* this is first vib, so make row of rot numbs */ 01490 for( i=0; i<=LimRot; ++i ) 01491 { 01492 fprintf(io,"\t%li",i); 01493 } 01494 fprintf(io,"\n"); 01495 } 01496 fprintf(io,"%li",iVibHi ); 01497 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=LimRot; ++iRotHi ) 01498 { 01499 fprintf(io,"\t%.3e", 01500 H2_X_colden[iVibHi][iRotHi]/hmi.H2_total ); 01501 } 01502 fprintf(io,"\n" ); 01503 } 01504 else 01505 { 01506 /* punch triplets - the default */ 01507 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=LimRot; ++iRotHi ) 01508 { 01509 fprintf(io,"%li\t%li\t%.1f\t%.3e\t%.3e\t%.3e\t%.3e\n", 01510 iVibHi , 01511 iRotHi , 01512 /* energy relative to 0,0, T1CM converts wavenumber to K */ 01513 energy_wn[iElecHi][iVibHi][iRotHi]*T1CM, 01514 /* these are column densities for actual molecule */ 01515 H2_X_colden[iVibHi][iRotHi] , 01516 H2_X_colden[iVibHi][iRotHi]/H2_stat[iElecHi][iVibHi][iRotHi] , 01517 /* these are same column densities but for LTE populations */ 01518 H2_X_colden_LTE[iVibHi][iRotHi] , 01519 H2_X_colden_LTE[iVibHi][iRotHi]/H2_stat[iElecHi][iVibHi][iRotHi]); 01520 } 01521 } 01522 } 01523 } 01524 } 01525 else if( (strcmp(chJOB , "H2pd" ) == 0) && (strcmp(chTime,"LAST") != 0) ) 01526 { 01527 /* punch PDR 01528 * output some PDR information (densities, rates) for each zone */ 01529 fprintf(io,"%.5e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n", 01530 /* depth in cm */ 01531 radius.depth_mid_zone , 01532 /* the computed ortho and para densities */ 01533 h2.ortho_density , 01534 h2.para_density , 01535 /* the Lyman Werner band dissociation, Tielens & Hollenbach */ 01536 hmi.H2_Solomon_dissoc_rate_TH85_H2g , 01537 /* the Lyman Werner band dissociation, Bertoldi & Draine */ 01538 hmi.H2_Solomon_dissoc_rate_BD96_H2g, 01539 /* the Lyman Werner band dissociation, big H2 mole */ 01540 hmi.H2_Solomon_dissoc_rate_BigH2_H2g); 01541 } 01542 else if( (strcmp(chJOB , "H2co" ) == 0) && (strcmp(chTime,"LAST") != 0) ) 01543 { 01544 /* punch H2 cooling - do heating cooling for each zone old new H2 */ 01545 fprintf(io,"%.5e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n", 01546 /* depth in cm */ 01547 radius.depth_mid_zone , 01548 /* total cooling, equal to total heating */ 01549 thermal.ctot , 01550 /* H2 destruction by Solomon process, TH85 rate */ 01551 hmi.H2_Solomon_dissoc_rate_TH85_H2g, 01552 /* H2 destruction by Solomon process, big H2 model rate */ 01553 hmi.H2_Solomon_dissoc_rate_BigH2_H2g + 01554 hmi.H2_Solomon_dissoc_rate_BigH2_H2s, 01555 /* H2 photodissociation heating, eqn A9 of Tielens & Hollenbach 1985a */ 01556 hmi.HeatH2Dish_TH85, 01557 /* heating due to dissociation of electronic excited states */ 01558 hmi.HeatH2Dish_BigH2 , 01559 /* cooling (usually neg and so heating) due to collisions within X */ 01560 hmi.HeatH2Dexc_TH85, 01561 hmi.HeatH2Dexc_BigH2 01562 ); 01563 01564 } 01565 else if( (strcmp(chJOB , "H2cr" ) == 0) && (strcmp(chTime,"LAST") != 0) ) 01566 { 01567 /* PUNCH H2 CREATION - show H2 creation processes for each zone */ 01568 /* >>chng 05 jul 15, TE, punch all H2 creation processes, unit cm-3 s-1 */ 01569 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", 01570 /* depth in cm */ 01571 radius.depth_mid_zone , 01572 /* creation cm-3 s-1, destruction rate, s-1 */ 01573 hmi.H2_rate_create, 01574 hmi.H2_rate_destroy, 01575 gv.rate_h2_form_grains_used_total * dense.xIonDense[ipHYDROGEN][0] / hmi.H2_rate_create, 01576 hmi.assoc_detach * hmi.Hmolec[ipMH]*hmi.Hmolec[ipMHm] / hmi.H2_rate_create, 01577 hmi.bh2dis * dense.xIonDense[ipHYDROGEN][0] / hmi.H2_rate_create, 01578 hmi.bh2h2p * dense.xIonDense[ipHYDROGEN][0] * hmi.Hmolec[ipMH2p] / hmi.H2_rate_create, 01579 hmi.radasc * dense.xIonDense[ipHYDROGEN][0] / hmi.H2_rate_create, 01580 hmi.h3ph2p * dense.xIonDense[ipHYDROGEN][0] * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create, 01581 hmi.h2phmh2h * hmi.Hmolec[ipMH2p] * hmi.Hmolec[ipMHm] / hmi.H2_rate_create, 01582 hmi.bh2h22hh2 * 2 * dense.xIonDense[ipHYDROGEN][0] * hmi.Hmolec[ipMH2g] / hmi.H2_rate_create, 01583 hmi.h3phmh2hh * hmi.Hmolec[ipMH3p] * hmi.Hmolec[ipMHm] / hmi.H2_rate_create, 01584 hmi.h3phm2h2 * hmi.Hmolec[ipMH3p] * hmi.Hmolec[ipMHm] / hmi.H2_rate_create, 01585 hmi.h32h2 * hmi.Hmolec[ipMH2p] * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create, 01586 hmi.eh3_h2h * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create, 01587 hmi.h3ph2hp * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create 01588 ); 01589 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", 01590 /*chemical network*/ 01591 /*light elements*/ 01592 co.H_CH_C_H2 * dense.xIonDense[ipHYDROGEN][0] / hmi.H2_rate_create, 01593 co.H_CHP_CP_H2 * dense.xIonDense[ipHYDROGEN][0] / hmi.H2_rate_create, 01594 co.H_CH2_CH_H2 * dense.xIonDense[ipHYDROGEN][0] / hmi.H2_rate_create, 01595 co.H_CH3P_CH2P_H2 * dense.xIonDense[ipHYDROGEN][0] / hmi.H2_rate_create, 01596 co.H_OH_O_H2 * dense.xIonDense[ipHYDROGEN][0] / hmi.H2_rate_create, 01597 co.Hminus_HCOP_CO_H2 * hmi.Hmolec[ipMHm] / hmi.H2_rate_create, 01598 co.Hminus_H3OP_H2O_H2 * hmi.Hmolec[ipMHm] / hmi.H2_rate_create, 01599 co.Hminus_H3OP_OH_H2_H * hmi.Hmolec[ipMHm] / hmi.H2_rate_create, 01600 co.HP_CH2_CHP_H2 * hmi.Hmolec[ipMHp] / hmi.H2_rate_create, 01601 co.HP_SiH_SiP_H2* hmi.Hmolec[ipMHp] / hmi.H2_rate_create, 01602 co.H2P_CH_CHP_H2 * hmi.Hmolec[ipMH2p] / hmi.H2_rate_create, 01603 co.H2P_CH2_CH2P_H2 * hmi.Hmolec[ipMH2p] / hmi.H2_rate_create, 01604 co.H2P_CO_COP_H2 * hmi.Hmolec[ipMH2p] / hmi.H2_rate_create, 01605 co.H2P_H2O_H2OP_H2 * hmi.Hmolec[ipMH2p] / hmi.H2_rate_create, 01606 co.H2P_O2_O2P_H2 * hmi.Hmolec[ipMH2p] / hmi.H2_rate_create, 01607 co.H2P_OH_OHP_H2 * hmi.Hmolec[ipMH2p] / hmi.H2_rate_create, 01608 co.H3P_C_CHP_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create, 01609 co.H3P_CH_CH2P_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create, 01610 co.H3P_CH2_CH3P_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create, 01611 co.H3P_OH_H2OP_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create, 01612 co.H3P_H2O_H3OP_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create, 01613 co.H3P_CO_HCOP_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create, 01614 co.H3P_O_OHP_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create, 01615 co.H3P_SiH_SiH2P_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create, 01616 co.H3P_SiO_SiOHP_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create, 01617 co.H_CH3_CH2_H2 * dense.xIonDense[ipHYDROGEN][0] / hmi.H2_rate_create, 01618 co.H_CH4P_CH3P_H2 * dense.xIonDense[ipHYDROGEN][0] / hmi.H2_rate_create, 01619 co.H_CH5P_CH4P_H2 * dense.xIonDense[ipHYDROGEN][0] / hmi.H2_rate_create, 01620 co.H2P_CH4_CH3P_H2 * hmi.Hmolec[ipMH2p] / hmi.H2_rate_create, 01621 co.H2P_CH4_CH4P_H2 * hmi.Hmolec[ipMH2p] / hmi.H2_rate_create, 01622 co.H3P_CH3_CH4P_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create, 01623 co.H3P_CH4_CH5P_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create, 01624 co.HP_CH4_CH3P_H2 * hmi.Hmolec[ipMHp] / hmi.H2_rate_create, 01625 /* heavy elements */ 01626 co.HP_HNO_NOP_H2 * hmi.Hmolec[ipMHp] / hmi.H2_rate_create, 01627 co.HP_HS_SP_H2 * hmi.Hmolec[ipMHp] / hmi.H2_rate_create, 01628 co.H_HSP_SP_H2 * dense.xIonDense[ipHYDROGEN][0] / hmi.H2_rate_create, 01629 co.H3P_NH_NH2P_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create, 01630 co.H3P_NH2_NH3P_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create, 01631 co.H3P_NH3_NH4P_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create, 01632 co.H3P_CN_HCNP_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create, 01633 co.H3P_NO_HNOP_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create, 01634 co.H3P_S_HSP_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create, 01635 co.H3P_CS_HCSP_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create, 01636 co.H3P_NO2_NOP_OH_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create, 01637 co.H2P_NH_NHP_H2 * hmi.Hmolec[ipMH2p] / hmi.H2_rate_create, 01638 co.H2P_NH2_NH2P_H2 * hmi.Hmolec[ipMH2p] / hmi.H2_rate_create, 01639 co.H2P_NH3_NH3P_H2 * hmi.Hmolec[ipMH2p] / hmi.H2_rate_create, 01640 co.H2P_CN_CNP_H2 * hmi.Hmolec[ipMH2p] / hmi.H2_rate_create, 01641 co.H2P_HCN_HCNP_H2 * hmi.Hmolec[ipMH2p] / hmi.H2_rate_create, 01642 co.H2P_NO_NOP_H2 * hmi.Hmolec[ipMH2p] / hmi.H2_rate_create, 01643 co.H3P_Cl_HClP_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create, 01644 co.H3P_HCl_H2ClP_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create, 01645 co.H2P_C2_C2P_H2 * hmi.Hmolec[ipMH2p] / hmi.H2_rate_create, 01646 co.Hminus_NH4P_NH3_H2 * hmi.Hmolec[ipMHm] / hmi.H2_rate_create, 01647 co.H3P_HCN_HCNHP_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create 01648 ); 01649 fprintf(io,"\t%.3e\t%.3e\n", 01650 hmi.H2_rate_destroy * hmi.H2_total / hmi.H2_rate_create, 01651 hmi.h2s_sp_decay 01652 ); 01653 } 01654 else if( (strcmp(chJOB , "H2ds" ) == 0) && (strcmp(chTime,"LAST") != 0) ) 01655 { 01656 /* punch H2 destruction - show H2 destruction processes for each zone 01657 * >>chng 05 nov 17, TE, added the new reaction H2s + O+ -> OH+ + H 01658 * >>chng 05 oct 04, TE, remove eh2hhm(was double) and include dissociation by electrons, H2g/H2s + e -> 2H 01659 * >>chng 05 jul 15, TE, punch all H2 destruction rates, weighted by fractional abundance of H2g, H2s, unit s-1 */ 01660 fprintf(io,"%.5e\t%.2e\t%.2e", 01661 /* depth in cm */ 01662 radius.depth_mid_zone , 01663 /* total H2 creation rate, cm-3 s-1 */ 01664 hmi.H2_rate_create, 01665 /* destruction rate, s-1 */ 01666 hmi.H2_rate_destroy ); 01667 /* this can be zero following a high-Te init temperature search 01668 * abort */ 01669 if( hmi.H2_total > 0. ) 01670 { 01671 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", 01672 /* H2 + e -> H- + H0 */ 01673 (hmi.assoc_detach_backwards_grnd * hmi.Hmolec[ipMH2g] + hmi.assoc_detach_backwards_exct * hmi.Hmolec[ipMH2s]) / hmi.H2_rate_destroy / hmi.H2_total, 01674 /*photons*/ 01675 hmi.H2_Solomon_dissoc_rate_used_H2g / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total, 01676 hmi.H2_Solomon_dissoc_rate_used_H2s / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2s] / hmi.H2_total, 01677 hmi.H2_photodissoc_used_H2s / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2s] / hmi.H2_total, 01678 hmi.H2_photodissoc_used_H2g / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total, 01679 /*electrons*/ 01680 (hmi.h2ge2h * hmi.Hmolec[ipMH2g] + hmi.h2se2h * hmi.Hmolec[ipMH2s]) / hmi.H2_rate_destroy / hmi.H2_total, 01681 /*H+*/ 01682 hmi.rh2h2p*dense.xIonDense[ipHYDROGEN][1] / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total, 01683 hmi.h2hph3p*dense.xIonDense[ipHYDROGEN][1] / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total, 01684 /*H0*/ 01685 (hmi.rh2dis * hmi.Hmolec[ipMH2g] + hmi.h2sh * hmi.Hmolec[ipMH2s]) * dense.xIonDense[ipHYDROGEN][0] / hmi.H2_rate_destroy / hmi.H2_total, 01686 /*CR*/ 01687 (hmi.CR_reac_H2g * hmi.Hmolec[ipMH2g] + hmi.CR_reac_H2s * hmi.Hmolec[ipMH2s]) / hmi.H2_rate_destroy / hmi.H2_total, 01688 /*He+,HeH+*/ 01689 hmi.rheph2hpheh*dense.xIonDense[ipHELIUM][1] / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total, 01690 hmi.heph2heh2p*dense.xIonDense[ipHELIUM][1] / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total, 01691 hmi.hehph2h3phe*hmi.Hmolec[ipMHeHp] / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total, 01692 /*H3+*/ 01693 hmi.h3petc*hmi.Hmolec[ipMH3p] / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total, 01694 /*H2+*/ 01695 hmi.h2ph3p*hmi.Hmolec[ipMH2p] / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total, 01696 /*H2s+H2g -> H2g + 2H*/ 01697 hmi.h2sh2g * hmi.Hmolec[ipMH2g] / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2s] / hmi.H2_total, 01698 /*H2g+H2g -> H2g + 2H*/ 01699 hmi.h2h22hh2*2*hmi.Hmolec[ipMH2g] / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total, 01700 /*H2s+H2s -> H2g + 2H*/ 01701 hmi.h2sh2sh2g2h*2*hmi.Hmolec[ipMH2s] / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2s] / hmi.H2_total, 01702 /*H2s+H2s -> H2s + 2H*/ 01703 hmi.h2sh2sh2s2h*2*hmi.Hmolec[ipMH2s] / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2s] / hmi.H2_total, 01704 /*chemical network*/ 01705 /*light elements*/ 01706 co.H2_CHP_CH2P_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total, 01707 co.H2_CH2P_CH3P_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total, 01708 co.H2_OHP_H2OP_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total, 01709 co.H2_H2OP_H3OP_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total, 01710 co.H2_COP_HCOP_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total, 01711 co.H2_OP_OHP_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total, 01712 co.H2_SiOP_SiOHP_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total, 01713 co.H2_C_CH_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total, 01714 co.H2_CP_CHP_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total, 01715 co.H2_CH_CH2_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total, 01716 co.H2_OH_H2O_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total, 01717 co.H2_O_OH_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total, 01718 co.H2s_CH_CH2_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2s] / hmi.H2_total, 01719 co.H2s_O_OH_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2s] / hmi.H2_total, 01720 co.H2s_OH_H2O_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2s] / hmi.H2_total, 01721 co.H2s_C_CH_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2s] / hmi.H2_total, 01722 co.H2s_CP_CHP_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2s] / hmi.H2_total, 01723 co.H2_CH2_CH3_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total, 01724 co.H2_CH3_CH4_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total, 01725 co.H2_CH4P_CH5P_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total, 01726 co.H2s_CH2_CH3_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2s] / hmi.H2_total, 01727 co.H2s_CH3_CH4_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2s] / hmi.H2_total, 01728 co.H2s_OP_OHP_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2s] / hmi.H2_total, 01729 /*heavy elements*/ 01730 co.H2_N_NH_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total, 01731 co.H2_NH_NH2_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total, 01732 co.H2_NH2_NH3_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total, 01733 co.H2_CN_HCN_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total, 01734 co.H2_NP_NHP_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total, 01735 co.H2_NHP_N_H3P / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total, 01736 co.H2_NHP_NH2P_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total, 01737 co.H2_NH2P_NH3P_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total, 01738 co.H2_NH3P_NH4P_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total, 01739 co.H2_CNP_HCNP_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total, 01740 co.H2_SP_HSP_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total, 01741 co.H2_CSP_HCSP_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total, 01742 co.H2_ClP_HClP_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total, 01743 co.H2_HClP_H2ClP_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total, 01744 co.H2_HCNP_HCNHP_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total 01745 ); 01746 fprintf(io,"\t%.4e\t%.4e", 01747 /*H2g/Htot, H2s/Htot chemical network and from big molecule model*/ 01748 hmi.Hmolec[ipMH2g] / hmi.H2_total, 01749 hmi.Hmolec[ipMH2s] / hmi.H2_total 01750 ); 01751 } 01752 fprintf(io,"\n"); 01753 } 01754 01755 else if( (strcmp(chJOB , "H2le" ) == 0) && (strcmp(chTime,"LAST") == 0) ) 01756 { 01757 /* punch H2 levels */ 01758 for( long int ipHi=0; ipHi < nLevels_per_elec[0]; ipHi++ ) 01759 { 01760 long int ipVJ_Hi = H2_ipX_ener_sort[ipHi]; 01761 iRotHi = ipRot_H2_energy_sort[ipVJ_Hi]; 01762 iVibHi = ipVib_H2_energy_sort[ipVJ_Hi]; 01763 double Asum , Csum[N_X_COLLIDER]; 01764 long int nColl; 01765 Asum = 0; 01766 for( nColl=0; nColl<N_X_COLLIDER; ++nColl ) 01767 Csum[nColl] = 0.; 01768 for( long int ipLo=0; ipLo<ipHi; ++ipLo ) 01769 { 01770 /* all lower levels */ 01771 long int ipVJ_Lo = H2_ipX_ener_sort[ipLo]; 01772 iRotLo = ipRot_H2_energy_sort[ipVJ_Lo]; 01773 iVibLo = ipVib_H2_energy_sort[ipVJ_Lo]; 01774 01775 /* radiative decays down */ 01776 if( ( abs(iRotHi-iRotLo) == 2 || (iRotHi-iRotLo) == 0 ) && iVibLo <= iVibHi && 01777 lgH2_line_exists[0][iVibHi][iRotHi][0][iVibLo][iRotLo] ) 01778 { 01779 Asum += H2Lines[0][iVibHi][iRotHi][0][iVibLo][iRotLo].Emis->Aul*( 01780 H2Lines[0][iVibHi][iRotHi][0][iVibLo][iRotLo].Emis->Pesc + 01781 H2Lines[0][iVibHi][iRotHi][0][iVibLo][iRotLo].Emis->Pdest + 01782 H2Lines[0][iVibHi][iRotHi][0][iVibLo][iRotLo].Emis->Pelec_esc); 01783 } 01784 /* all collisions down */ 01785 mr5ci H2cr = H2_CollRate.begin(iVibHi,iRotHi,iVibLo,iRotLo); 01786 for( nColl=0; nColl<N_X_COLLIDER; ++nColl ) 01787 Csum[nColl] += H2cr[nColl]; 01788 } 01789 01790 /* punch H2 level energies */ 01791 fprintf(io,"%li\t%li\t%.2f\t%li\t%.3e", 01792 iVibHi , iRotHi, 01793 energy_wn[0][iVibHi][iRotHi], 01794 (long)H2_stat[0][iVibHi][iRotHi], 01795 Asum ); 01796 for( nColl=0; nColl<N_X_COLLIDER; ++nColl ) 01797 /* sum over all lower levels */ 01798 fprintf(io,"\t%.3e",Csum[nColl]); 01799 fprintf(io,"\n"); 01800 } 01801 } 01802 01803 else if( (strcmp(chJOB , "H2ra" ) == 0) && (strcmp(chTime,"LAST") != 0) ) 01804 { 01805 /* punch h2 rates - some rates and lifetimes */ 01806 double sumpop = 0. , sumlife = 0.; 01807 01808 /* this block, find lifetime against photo excitation into excited electronic states */ 01809 iElecLo = 0; 01810 iVibLo = 0; 01811 if( h2.lgH2ON && hmi.lgBigH2_evaluated ) 01812 { 01813 for( iElecHi=1; iElecHi<mole.n_h2_elec_states; ++iElecHi ) 01814 { 01815 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi ) 01816 { 01817 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi ) 01818 { 01819 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=h2.nRot_hi[iElecLo][iVibLo]; ++iRotLo ) 01820 { 01821 /* only do if radiative transition exists */ 01822 if( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].ipCont > 0 ) 01823 { 01824 sumlife += 01825 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->pump * 01826 H2_populations[iElecLo][iVibLo][iRotLo]; 01827 sumpop += 01828 H2_populations[iElecLo][iVibLo][iRotLo]; 01829 } 01830 } 01831 } 01832 } 01833 } 01834 } 01835 01836 /* continue output from punch h2 rates command */ 01837 /* find photoexcitation rates from v=0 */ 01838 /* PDR information for each zone */ 01839 fprintf(io, 01840 "%.5e\t%.3e\t%.3e\t%.3e\t%li", 01841 /* depth in cm */ 01842 radius.depth_mid_zone , 01843 /* the column density (cm^-2) in H2 */ 01844 colden.colden[ipCOL_H2g]+colden.colden[ipCOL_H2s], 01845 /* this is a special form of column density - should be proportional 01846 * to total shielding */ 01847 colden.coldenH2_ov_vel , 01848 /* visual extinction due to dust alone, of point source (star)*/ 01849 rfield.extin_mag_V_point, 01850 /* number of large molecule evaluations in this zone */ 01851 h2.nCallH2_this_zone ); 01852 fprintf(io, 01853 "\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e", 01854 /* total H2 fraction */ 01855 hmi.H2_total/dense.gas_phase[ipHYDROGEN] , 01856 /* chemistry renorm factor */ 01857 H2_renorm_chemistry, 01858 /* rate H2 forms on grains */ 01859 gv.rate_h2_form_grains_used_total , 01860 /* rate H2 forms by H minus route */ 01861 hmi.Hmolec[ipMHm]*1.35e-9, 01862 /* H2 destruction by Solomon process, TH85 rate */ 01863 hmi.H2_Solomon_dissoc_rate_TH85_H2g + hmi.H2_Solomon_dissoc_rate_TH85_H2s, 01864 /* H2 destruction by Solomon process, Bertoldi & Draine rate */ 01865 hmi.H2_Solomon_dissoc_rate_BD96_H2g + hmi.H2_Solomon_dissoc_rate_BD96_H2s, 01866 /* H2 destruction by Solomon process, Elwert et al. in preparation */ 01867 hmi.H2_Solomon_dissoc_rate_ELWERT_H2g + hmi.H2_Solomon_dissoc_rate_ELWERT_H2g, 01868 /* H2 destruction by Solomon process, big H2 model rate */ 01869 hmi.H2_Solomon_dissoc_rate_BigH2_H2g + hmi.H2_Solomon_dissoc_rate_BigH2_H2s, 01870 /* rate s-1 H2 electronic excit states decay into H2g */ 01871 hmi.H2_Solomon_elec_decay_H2g , 01872 /* rate s-1 H2 electronic excit states decay into H2s */ 01873 hmi.H2_Solomon_elec_decay_H2s 01874 ); 01875 fprintf(io, 01876 "\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e", 01877 /* The TH85 estimate of the radiation field relative to the Habing value */ 01878 hmi.UV_Cont_rel2_Habing_TH85_depth, 01879 /* The DB96 estimate of the radiation field relative to the Habing value */ 01880 hmi.UV_Cont_rel2_Draine_DB96_depth, 01881 /* cosmic ray ionization rate */ 01882 secondaries.csupra[ipHYDROGEN][0]*0.93, 01883 sumlife/SDIV( sumpop ) , 01884 hmi.H2_Solomon_dissoc_rate_BD96_H2g/SDIV(hmi.UV_Cont_rel2_Habing_TH85_depth) , 01885 hmi.H2_Solomon_dissoc_rate_BigH2_H2g/SDIV(hmi.UV_Cont_rel2_Habing_TH85_depth), 01886 hmi.H2_Solomon_dissoc_rate_BigH2_H2g/SDIV(hmi.UV_Cont_rel2_Habing_spec_depth), 01887 hmi.H2_rate_destroy); 01888 fprintf(io, 01889 "\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n", 01890 hmi.HeatH2Dish_TH85, 01891 hmi.HeatH2Dexc_TH85, 01892 hmi.HeatH2Dish_BigH2, 01893 hmi.HeatH2Dexc_BigH2, 01894 thermal.htot); 01895 } 01896 /* punch h2 solomon */ 01897 else if( (strcmp(chJOB , "H2so" ) == 0) && (strcmp(chTime,"LAST") != 0) ) 01898 { 01899 /* remember as many as NSOL lines contributing to total Solomon process */ 01900 # define NSOL 100 01901 double sum, one; 01902 long int jlosave[NSOL] , ivlosave[NSOL], 01903 iehisave[NSOL] ,jhisave[NSOL] , ivhisave[NSOL], 01904 nsave, 01905 ipOrdered[NSOL]; 01906 int nFail; 01907 int i; 01908 realnum fsave[NSOL], wlsave[NSOL]; 01909 /* Solomon process, and where it came from */ 01910 fprintf(io,"%.5e\t%.3e", 01911 /* depth in cm */ 01912 radius.depth_mid_zone , 01913 /* H2 destruction by Solomon process, big H2 model rate */ 01914 hmi.H2_Solomon_dissoc_rate_BigH2_H2g + 01915 hmi.H2_Solomon_dissoc_rate_BigH2_H2s); 01916 sum = 0.; 01917 iElecLo = 0; 01918 /* find sum of all radiative exits from X into excited electronic states */ 01919 if( h2.lgH2ON && hmi.lgBigH2_evaluated ) 01920 { 01921 for( iElecHi=1; iElecHi<mole.n_h2_elec_states; ++iElecHi ) 01922 { 01923 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi ) 01924 { 01925 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi ) 01926 { 01927 long int nv = h2.nVib_hi[iElecLo]; 01928 for( iVibLo=0; iVibLo<=nv; ++iVibLo ) 01929 { 01930 long nr = h2.nRot_hi[iElecLo][iVibLo]; 01931 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo ) 01932 { 01933 /* only do if radiative transition exists */ 01934 if( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].ipCont > 0 ) 01935 { 01936 one = H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Lo->Pop * 01937 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->pump; 01938 sum += one; 01939 } 01940 } 01941 } 01942 } 01943 } 01944 } 01945 01946 /* make sure it is safe to div by sum */ 01947 sum = SDIV( sum ); 01948 nsave = 0; 01949 /* now loop back over X and print all those which contribute more than FRAC of the total */ 01950 # define FRAC 0.01 01951 for( iElecHi=1; iElecHi<mole.n_h2_elec_states; ++iElecHi ) 01952 { 01953 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi ) 01954 { 01955 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi ) 01956 { 01957 long int nv = h2.nVib_hi[iElecLo]; 01958 for( iVibLo=0; iVibLo<=nv; ++iVibLo ) 01959 { 01960 long nr = h2.nRot_hi[iElecLo][iVibLo]; 01961 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo ) 01962 { 01963 /* only do if radiative transition exists */ 01964 if( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].ipCont > 0 ) 01965 { 01966 one = H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Lo->Pop * 01967 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->pump; 01968 if( one/sum > FRAC && nsave<NSOL) 01969 { 01970 fsave[nsave] = (realnum)(one/sum); 01971 jlosave[nsave] = iRotLo; 01972 ivlosave[nsave] = iVibLo; 01973 jhisave[nsave] = iRotHi; 01974 ivhisave[nsave] = iVibHi; 01975 iehisave[nsave] = iElecHi; 01976 wlsave[nsave] = H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].WLAng; 01977 ++nsave; 01978 /*fprintf(io,"\t%li\t%li\t%li\t%li\t%li\t%.3f", 01979 iElecHi,iVibHi,iRotHi,iVibLo , iRotLo , one/sum );*/ 01980 } 01981 } 01982 } 01983 } 01984 } 01985 } 01986 }/* iElecHi */ 01987 /* now sort these into decreasing order */ 01988 01989 /* now sort by decreasing importance */ 01990 /*spsort netlib routine to sort array returning sorted indices */ 01991 spsort( 01992 /* input array to be sorted */ 01993 fsave, 01994 /* number of values in x */ 01995 nsave, 01996 /* permutation output array */ 01997 ipOrdered, 01998 /* flag saying what to do - 1 sorts into increasing order, not changing 01999 * the original routine */ 02000 -1, 02001 /* error condition, should be 0 */ 02002 &nFail); 02003 02004 /* print ratio of pumps to dissociations - this is 9:1 in TH85 */ 02005 /*>>chng 05 jul 20, TE, punch average energy in H2s and renormalization factors for H2g and H2s */ 02006 /* >>chng 05 sep 16, TE, chng denominator to do g and s with proper dissoc rates */ 02007 fprintf(io,"\t%.3f\t%.3f\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e", 02008 /* this is sum of photons and CRs */ 02009 (sum + secondaries.csupra[ipHYDROGEN][0]*2.02f)/SDIV((hmi.H2_Solomon_dissoc_rate_BigH2_H2g * hmi.Hmolec[ipMH2g] + 02010 hmi.H2_Solomon_dissoc_rate_BigH2_H2s * hmi.Hmolec[ipMH2s]) ), 02011 /* this is sum of photons and CRs */ 02012 (sum + secondaries.csupra[ipHYDROGEN][0]*2.02f) /SDIV((hmi.H2_Solomon_dissoc_rate_BigH2_H2g * hmi.H2g_BigH2 + 02013 hmi.H2_Solomon_dissoc_rate_BigH2_H2s * hmi.H2s_BigH2) ), 02014 hmi.H2_BigH2_H2g_av, hmi.H2_BigH2_H2s_av, 02015 hmi.H2_chem_BigH2_H2g, hmi.H2_chem_BigH2_H2s, 02016 hmi.H2g_BigH2/SDIV(hmi.H2_total_BigH2), hmi.H2s_BigH2/SDIV(hmi.H2_total_BigH2) 02017 ); 02018 for( i=0; i<nsave; ++i ) 02019 { 02020 ip = ipOrdered[i]; 02021 /*lint -e644 not init */ 02022 fprintf(io,"\t%li\t%li\t%li\t%li\t%li\t%.3f\t%.3f", 02023 iehisave[ip],ivhisave[ip],jhisave[ip],ivlosave[ip] , jlosave[ip] , fsave[ip] , wlsave[ip] ); 02024 /*lint +e644 not init */ 02025 } 02026 fprintf(io,"\n"); 02027 } 02028 /*fprintf(io,"DEBUG tau\t%.3e\t%.3f\n", 02029 H2Lines[1][0][1][0][0][0].Emis->TauIn, H2Lines[1][0][1][0][0][0].WLAng); */ 02030 # undef NSOL 02031 } 02032 02033 else if( (strcmp(chJOB , "H2te" ) == 0) && (strcmp(chTime,"LAST") != 0) ) 02034 { 02035 /* punch h2 temperatures */ 02036 double pop_ratio10,pop_ratio20,pop_ratio30,pop_ratio31,pop_ratio40; 02037 double T10,T20,T30,T31,T40; 02038 /* subscript"sum" denotes integrated quantities */ 02039 double T10_sum,T20_sum,T30_sum,T31_sum,T40_sum; 02040 double pop_ratio10_sum,pop_ratio20_sum,pop_ratio30_sum,pop_ratio31_sum,pop_ratio40_sum; 02041 if( h2.lgH2ON && h2.nCallH2_this_zone ) 02042 { 02043 double energyK = T1CM*(energy_wn[0][0][1] - energy_wn[0][0][0]); 02044 /* the ratio of H2_populations of J=1 to 0 */ 02045 pop_ratio10 = H2_populations[0][0][1]/SDIV(H2_populations[0][0][0]); 02046 pop_ratio10_sum = H2_X_colden[0][1]/SDIV(H2_X_colden[0][0]); 02047 /* the corresponding temperature */ 02048 T10 = -170.5/log(SDIV(pop_ratio10) * H2_stat[0][0][0]/H2_stat[0][0][1]); 02049 T10_sum = -170.5/log(SDIV(pop_ratio10_sum) * H2_stat[0][0][0]/H2_stat[0][0][1]); 02050 02051 energyK = T1CM*(energy_wn[0][0][2] - energy_wn[0][0][0]); 02052 pop_ratio20 = H2_populations[0][0][2]/SDIV(H2_populations[0][0][0]); 02053 T20 = -energyK/log(SDIV(pop_ratio20) * H2_stat[0][0][0]/H2_stat[0][0][2]); 02054 02055 pop_ratio20_sum = H2_X_colden[0][2]/SDIV(H2_X_colden[0][0]); 02056 T20_sum = -energyK/log(SDIV(pop_ratio20_sum) * H2_stat[0][0][0]/H2_stat[0][0][2]); 02057 02058 energyK = T1CM*(energy_wn[0][0][3] - energy_wn[0][0][0]); 02059 pop_ratio30 = H2_populations[0][0][3]/SDIV(H2_populations[0][0][0]); 02060 T30 = -energyK/log(SDIV(pop_ratio30) * H2_stat[0][0][0]/H2_stat[0][0][3]); 02061 02062 pop_ratio30_sum = H2_X_colden[0][3]/SDIV(H2_X_colden[0][0]); 02063 T30_sum = -energyK/log(SDIV(pop_ratio30_sum) * H2_stat[0][0][0]/H2_stat[0][0][3]); 02064 02065 energyK = T1CM*(energy_wn[0][0][3] - energy_wn[0][0][1]); 02066 pop_ratio31 = H2_populations[0][0][3]/SDIV(H2_populations[0][0][1]); 02067 T31 = -energyK/log(SDIV(pop_ratio31) * H2_stat[0][0][1]/H2_stat[0][0][3]); 02068 02069 pop_ratio31_sum = H2_X_colden[0][3]/SDIV(H2_X_colden[0][1]); 02070 T31_sum = -energyK/log(SDIV(pop_ratio31_sum) * H2_stat[0][0][1]/H2_stat[0][0][3]); 02071 02072 energyK = T1CM*(energy_wn[0][0][4] - energy_wn[0][0][0]); 02073 pop_ratio40 = H2_populations[0][0][4]/SDIV(H2_populations[0][0][0]); 02074 T40 = -energyK/log(SDIV(pop_ratio40) * H2_stat[0][0][0]/H2_stat[0][0][4]); 02075 02076 pop_ratio40_sum = H2_X_colden[0][4]/SDIV(H2_X_colden[0][0]); 02077 T40_sum = -energyK/log(SDIV(pop_ratio40_sum) * H2_stat[0][0][0]/H2_stat[0][0][4]); 02078 } 02079 else 02080 { 02081 pop_ratio10 = 0.; 02082 pop_ratio10_sum = 0.; 02083 T10 = 0.; 02084 T20 = 0.; 02085 T30 = 0.; 02086 T31 = 0.; 02087 T40 = 0.; 02088 T10_sum = 0.; 02089 T20_sum = 0.; 02090 T30_sum = 0.; 02091 T31_sum = 0.; 02092 T40_sum = 0.; 02093 } 02094 02095 /* various temperatures for neutral/molecular gas */ 02096 fprintf( io, 02097 "%.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" , 02098 /* depth in cm */ 02099 radius.depth_mid_zone , 02100 /* total H2 fraction */ 02101 hmi.H2_total/dense.gas_phase[ipHYDROGEN] , 02102 /* ratio of H2_populations of 1 to 0 only */ 02103 pop_ratio10 , 02104 /* sum of all ortho and para */ 02105 h2.ortho_density / SDIV(h2.para_density), 02106 T10,T20,T30,T31,T40, 02107 phycon.te , 02108 hyperfine.Tspin21cm,T10_sum,T20_sum,T30_sum,T31_sum,T40_sum ); 02109 } 02110 else if( (strcmp(chJOB , "H2ln" ) == 0) && (strcmp(chTime,"LAST") == 0) ) 02111 { 02112 /* punch H2 lines - output the full emission-line spectrum */ 02113 double thresh; 02114 double renorm; 02115 /* first test, is H2 turned on? Second test, have lines arrays 02116 * been set up - nsum is negative if abort occurs before lines 02117 * are set up */ 02118 if( h2.lgH2ON && LineSave.nsum > 0) 02119 { 02120 ASSERT( LineSave.ipNormWavL >= 0 ); 02121 /* get the normalization line */ 02122 if( LineSv[LineSave.ipNormWavL].sumlin[LineSave.lgLineEmergent] > SMALLFLOAT ) 02123 renorm = LineSave.ScaleNormLine/LineSv[LineSave.ipNormWavL].sumlin[LineSave.lgLineEmergent]; 02124 else 02125 renorm = 1.; 02126 02127 if( renorm > SMALLFLOAT ) 02128 { 02129 /* this is threshold for faintest line, normally 0, set with 02130 * number on punch H2 command */ 02131 thresh = thresh_punline_h2/(realnum)renorm; 02132 } 02133 else 02134 thresh = 0.f; 02135 02136 /* punch H2 line intensities at end of iteration 02137 * h2.nElecLevelOutput is electronic level with 1 for ground, so this loop is < h2.nElecLevelOutput */ 02138 for( iElecHi=0; iElecHi < h2.nElecLevelOutput; ++iElecHi ) 02139 { 02140 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi ) 02141 { 02142 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi ) 02143 { 02144 /* now the lower levels */ 02145 /* NB - X is the only lower level considered here, since we are only 02146 * concerned with excited electronic levels as a photodissociation process 02147 * code exists to relax this assumption - simply change following to iElecHi */ 02148 long int lim_elec_lo = 0; 02149 for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo ) 02150 { 02151 long int nv = h2.nVib_hi[iElecLo]; 02152 if( iElecLo==iElecHi ) 02153 nv = iVibHi; 02154 for( iVibLo=0; iVibLo<=nv; ++iVibLo ) 02155 { 02156 long nr = h2.nRot_hi[iElecLo][iVibLo]; 02157 if( iElecLo==iElecHi && iVibHi==iVibLo ) 02158 nr = iRotHi-1; 02159 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<nr; ++iRotLo ) 02160 { 02161 if( H2_SaveLine[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] > thresh ) 02162 { 02163 /* air wavelength in microns */ 02164 /* WLAng contains correction for index of refraction of air */ 02165 double wl = H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].WLAng/1e4; 02166 /*ASSERT( abs(iRotHi-iRotLo)<=2 );*/ 02167 02168 fprintf(io, "%li-%li %c(%li)", 02169 iVibHi , 02170 iVibLo , 02171 chMolBranch( iRotHi , iRotLo ) , 02172 iRotLo ); 02173 fprintf( io, "\t%ld\t%ld\t%ld\t%ld\t%ld\t%ld", 02174 iElecHi , iVibHi , iRotHi , iElecLo , iVibLo , iRotLo); 02175 /* WLAng contains correction for index of refraction of air */ 02176 fprintf( io, "\t%.7f\t", wl ); 02177 /*prt_wl print floating wavelength in Angstroms, in output format */ 02178 prt_wl( io , H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].WLAng ); 02179 fprintf( io, "\t%.3f\t%.3e", 02180 log10(MAX2(1e-37,H2_SaveLine[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo])) + radius.Conv2PrtInten, 02181 H2_SaveLine[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo]*renorm ); 02182 /* excitation energy of upper level in K */ 02183 fprintf( io, "\t%.3f", energy_wn[iElecHi][iVibHi][iRotHi]*T1CM ); 02184 /* the product h nu * Aul */ 02185 ASSERT( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].ipCont > 0 ); 02186 fprintf( io, "\t%.3e", H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->Aul* 02187 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyErg * 02188 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Hi->g); 02189 fprintf( io, "\n"); 02190 } 02191 } 02192 } 02193 } 02194 } 02195 } 02196 } 02197 } 02198 } 02199 else if( (strcmp(chJOB , "H2sp" ) == 0) ) 02200 { 02201 iVib = 0; 02202 iRot = 0; 02203 # if 0 02204 /* punch h2 special */ 02205 fprintf(io,"%.4e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n", 02206 radius.depth_mid_zone , 02207 H2_populations[0][iVib][iRot] , 02208 radius.depth_mid_zone * H2_populations[0][iVib][iRot] , 02209 H2_rad_rate_out[0][iVib][iRot] , 02210 H2_rad_rate_in[iVib][iRot] , 02211 H2_col_rate_out_ old[iVib][iRot] , 02212 H2_col_rate_in_ old[iVib][iRot] ); 02213 # endif 02214 fprintf(io,"%.4e\t%.2e\t%.2e\t%.2e\t%.2e\n", 02215 radius.depth_mid_zone , 02216 H2_populations[0][iVib][iRot] , 02217 H2Lines[1][1][1][0][iVib][iRot].Emis->pump, 02218 H2Lines[1][1][1][0][iVib][iRot].Emis->TauIn, 02219 H2Lines[1][1][1][0][iVib][iRot].Emis->TauCon); 02220 } 02221 return; 02222 } 02223 /*cdH2_Line determines intensity and luminosity of and H2 line. The first 02224 * six arguments give the upper and lower quantum designation of the levels. 02225 * The function returns 1 if we found the line, 02226 * and false==0 if we did not find the line because ohoto-para transition 02227 * or upper level has lower energy than lower level */ 02228 long int cdH2_Line( 02229 /* indices for the upper level */ 02230 long int iElecHi, 02231 long int iVibHi , 02232 long int iRotHi , 02233 /* indices for lower level */ 02234 long int iElecLo, 02235 long int iVibLo , 02236 long int iRotLo , 02237 /* linear intensity relative to normalization line*/ 02238 double *relint, 02239 /* log of luminosity or intensity of line */ 02240 double *absint ) 02241 { 02242 02243 DEBUG_ENTRY( "cdH2_Line()" ); 02244 02245 /* these will be return values if we can't find the line */ 02246 *relint = 0.; 02247 *absint = 0.; 02248 02249 /* for now both electronic levels must be zero */ 02250 if( iElecHi!=0 || iElecLo!=0 ) 02251 { 02252 return 0; 02253 } 02254 02255 /* check that energy of first level is higher than energy of second level */ 02256 if( energy_wn[iElecHi][iVibHi][iRotHi] < energy_wn[iElecLo][iVibHi][iRotHi] ) 02257 { 02258 return 0; 02259 } 02260 02261 /* check that ortho-para does not change */ 02262 if( H2_lgOrtho[iElecHi][iVibHi][iRotHi] - H2_lgOrtho[iElecLo][iVibLo][iRotLo] != 0 ) 02263 { 02264 return 0; 02265 } 02266 02267 /* exit if lines does not exist */ 02268 if( !lgH2_line_exists[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] ) 02269 { 02270 return 0; 02271 } 02272 02273 ASSERT( LineSave.ipNormWavL >= 0 ); 02274 /* does the normalization line have a positive intensity*/ 02275 if( LineSv[LineSave.ipNormWavL].sumlin[LineSave.lgLineEmergent] > 0. ) 02276 { 02277 *relint = H2_SaveLine[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo]/ 02278 LineSv[LineSave.ipNormWavL].sumlin[LineSave.lgLineEmergent] * LineSave.ScaleNormLine; 02279 } 02280 else 02281 { 02282 *relint = 0.; 02283 } 02284 02285 /* return log of line intensity if it is positive */ 02286 if( H2_SaveLine[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] > 0. ) 02287 { 02288 *absint = log10(H2_SaveLine[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo]) + 02289 radius.Conv2PrtInten; 02290 } 02291 else 02292 { 02293 /* line intensity is actually zero, return small number */ 02294 *absint = -37.; 02295 } 02296 /* this indicates success */ 02297 return 1; 02298 }