00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #include "cddefines.h"
00022 #include "physconst.h"
00023 #include "punch.h"
00024 #include "hmi.h"
00025 #include "path.h"
00026 #include "prt.h"
00027 #include "secondaries.h"
00028 #include "grainvar.h"
00029 #include "phycon.h"
00030 #include "rfield.h"
00031 #include "hyperfine.h"
00032 #include "thermal.h"
00033 #include "lines.h"
00034 #include "lines_service.h"
00035 #include "dense.h"
00036 #include "radius.h"
00037 #include "colden.h"
00038 #include "taulines.h"
00039 #include "h2.h"
00040 #include "h2_priv.h"
00041 #include "cddrive.h"
00042 #include "mole.h"
00043
00044
00045
00046
00047
00048 static char chlgPara[2]={'P','O'};
00049
00050
00051 static float thresh_punline_h2;
00052
00053
00054 void H2_LinesAdd(void)
00055 {
00056
00057 int iRotHi, iVibHi, iElecHi ,iRotLo, iVibLo, iElecLo;
00058
00059
00060 if( !h2.lgH2ON )
00061 return;
00062
00063 DEBUG_ENTRY( "H2_LinesAdd()" );
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080 lindst( H2Lines[0][1][6][0][0][4].xIntensity , H2Lines[0][1][6][0][0][4].WLAng , "H2 " ,
00081 H2Lines[0][1][6][0][0][4].ipCont,'i',false );
00082
00083 lindst( H2Lines[0][1][5][0][0][3].xIntensity , H2Lines[0][1][5][0][0][3].WLAng , "H2 " ,
00084 H2Lines[0][1][5][0][0][3].ipCont,'i',false );
00085
00086 lindst( H2Lines[0][1][4][0][0][2].xIntensity , H2Lines[0][1][4][0][0][2].WLAng , "H2 " ,
00087 H2Lines[0][1][4][0][0][2].ipCont,'i',false );
00088
00089 lindst( H2Lines[0][1][3][0][0][1].xIntensity , H2Lines[0][1][3][0][0][1].WLAng , "H2 " ,
00090 H2Lines[0][1][3][0][0][1].ipCont,'i',false );
00091
00092 lindst( H2Lines[0][1][2][0][0][0].xIntensity , H2Lines[0][1][2][0][0][0].WLAng , "H2 " ,
00093 H2Lines[0][1][2][0][0][0].ipCont,'i',false );
00094
00095
00096 lindst( H2Lines[0][1][2][0][0][2].xIntensity , H2Lines[0][1][2][0][0][2].WLAng , "H2 " ,
00097 H2Lines[0][1][2][0][0][2].ipCont,'i',false );
00098
00099 lindst( H2Lines[0][1][1][0][0][1].xIntensity , H2Lines[0][1][1][0][0][1].WLAng , "H2 " ,
00100 H2Lines[0][1][1][0][0][1].ipCont,'i',false );
00101
00102
00103
00104
00105 for( iElecHi=0; iElecHi<h2.nElecLevelOutput; ++iElecHi )
00106 {
00107 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
00108 {
00109 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
00110 {
00111 long int lim_elec_lo = 0;
00112
00113
00114
00115
00116 for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
00117 {
00118
00119
00120 long int nv = h2.nVib_hi[iElecLo];
00121 if( iElecLo==iElecHi )
00122 nv = iVibHi;
00123 for( iVibLo=0; iVibLo<=nv; ++iVibLo )
00124 {
00125 long nr = h2.nRot_hi[iElecLo][iVibLo];
00126 if( iElecLo==iElecHi && iVibHi==iVibLo )
00127 nr = iRotHi-1;
00128
00129 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
00130 {
00131 if( lgH2_line_exists[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] )
00132 {
00133
00134 PutLine(&H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo]);
00135 if( LineSave.ipass == 0 )
00136 {
00137 H2_SaveLine[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] = 0.;
00138 }
00139 else if( LineSave.ipass == 1 )
00140 {
00141 H2_SaveLine[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] += (float)(
00142 radius.dVeff*H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].xIntensity);
00143 }
00144 }
00145 }
00146 }
00147 }
00148 }
00149 }
00150 }
00151 DEBUG_EXIT( "H2_LinesAdd()" );
00152 return;
00153 }
00154
00155
00156 void H2_ParsePunch( char *chCard )
00157 {
00158 long int i , nelem;
00159 bool lgEOL;
00160
00161 DEBUG_ENTRY( "H2_ParsePunch()" );
00162
00163 i = 5;
00164 nelem = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
00165 if( nelem != 2 )
00166 {
00167 fprintf( ioQQQ, " The first number on this line must be the 2 in H2\n Sorry.\n" );
00168 puts( "[Stop in ParsePunch]" );
00169 cdEXIT(EXIT_FAILURE);
00170 }
00171
00172 if( nMatch("COLU",chCard) )
00173 {
00174
00175 strcpy( punch.chPunch[punch.npunch], "H2cl" );
00176
00177
00178
00179
00180
00181 punch.punarg[punch.npunch][0] = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
00182
00183
00184 punch.punarg[punch.npunch][1] = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
00185
00186
00187 if( nMatch( "MATR" , chCard ) )
00188 {
00189
00190 punch.punarg[punch.npunch][2] = 1;
00191 fprintf( punch.ipPnunit[punch.npunch], "#vib\trot\tcolumn density\n" );
00192 }
00193 else
00194 {
00195
00196 punch.punarg[punch.npunch][2] = -1;
00197 fprintf( punch.ipPnunit[punch.npunch], "#vib\trot\tEner(K)\tcolden\tcolden/stat wght\tLTE colden\tLTE colden/stat wght\n" );
00198 }
00199 }
00200 else if( nMatch("COOL",chCard) )
00201 {
00202
00203 strcpy( punch.chPunch[punch.npunch], "H2co" );
00204 fprintf( punch.ipPnunit[punch.npunch],
00205 "#H2 depth\ttot cool\tTH Sol\tBig Sol\tTH pht dis\tpht dis\tTH Xcool\tXcool \n" );
00206 }
00207
00208 else if( nMatch("CREA",chCard) )
00209 {
00210
00211 strcpy( punch.chPunch[punch.npunch], "H2cr" );
00212 fprintf( punch.ipPnunit[punch.npunch],
00213 "#H2 depth\tH2_rate_create\tH2_rate_destroy\trate_h2_form_grains_used_total\tassoc_detach");
00214 fprintf( punch.ipPnunit[punch.npunch],
00215 "\tbh2dis\tbh2h2p\tradasc\th3ph2p\th2phmh2h\tbh2h22hh2\th3phmh2hh\th3phm2h2\th32h2\teh3_h2h\th3ph2hp");
00216 fprintf( punch.ipPnunit[punch.npunch],
00217 "\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");
00218 fprintf( punch.ipPnunit[punch.npunch],
00219 "\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");
00220 fprintf( punch.ipPnunit[punch.npunch],
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 fprintf( punch.ipPnunit[punch.npunch],
00223 "\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");
00224 fprintf( punch.ipPnunit[punch.npunch],
00225 "\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");
00226 fprintf( punch.ipPnunit[punch.npunch],
00227 "\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");
00228 fprintf( punch.ipPnunit[punch.npunch],
00229 "\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");
00230 fprintf( punch.ipPnunit[punch.npunch],
00231 "\tH3P_Cl_HClP_H2\tH3P_HCl_H2ClP_H2\tH2P_C2_C2P_H2\tHminus_NH4P_NH3_H2\tH3P_HCN_HCNHP_H2");
00232 fprintf( punch.ipPnunit[punch.npunch],
00233 "\tdestruction/creation\tav Einstein A\n");
00234 }
00235 else if( nMatch("DEST",chCard) )
00236 {
00237
00238 strcpy( punch.chPunch[punch.npunch], "H2ds" );
00239 fprintf( punch.ipPnunit[punch.npunch],
00240 "#depth\ttot H2 rate create\ttot H2 rate destroy\ttot H- backwards\tSolomon H2g\tSolomon H2s\tphotodissoc H2s\tphotodissoc H2g");
00241 fprintf( punch.ipPnunit[punch.npunch],
00242 "\te- dissoc\trh2h2p\th2hph3p\tH0 dissoc\tCR\trheph2hpheh\theph2heh2p\thehph2h3phe\th3petc\tH2Ph3p");
00243 fprintf( punch.ipPnunit[punch.npunch],
00244 "\th2sh2g\th2h22hh2\th2sh2sh2g2h\th2sh2sh2s2h\tH2_CHP_CH2P_H\tH2_CH2P_CH3P_H\tH2_OHP_H2OP_H\tH2_H2OP_H3OP_H\tH2_COP_HCOP_H");
00245 fprintf( punch.ipPnunit[punch.npunch],
00246 "\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");
00247 fprintf( punch.ipPnunit[punch.npunch],
00248 "\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");
00249 fprintf( punch.ipPnunit[punch.npunch],
00250 "\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");
00251 fprintf( punch.ipPnunit[punch.npunch],
00252 "\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");
00253 fprintf( punch.ipPnunit[punch.npunch],
00254 "\tH2_ClP_HClP_H\tH2_HClP_H2ClP_H\tH2_HCNP_HCNHP_H");
00255 fprintf( punch.ipPnunit[punch.npunch],
00256 "\tfrac H2g\tfrac H2s\n");
00257 }
00258
00259 else if( nMatch("HEAT",chCard) )
00260 {
00261
00262 strcpy( punch.chPunch[punch.npunch], "H2he" );
00263 fprintf( punch.ipPnunit[punch.npunch],
00264 "#H2 depth\ttot Heat\tHeat(big)\tHeat(TH85)\tDissoc(Big)\tDissoc(TH85) \n" );
00265 }
00266
00267 else if( nMatch("LEVE",chCard) )
00268 {
00269
00270 strcpy( punch.chPunch[punch.npunch], "H2le" );
00271 fprintf( punch.ipPnunit[punch.npunch],
00272 "#H2 energy(wn)\tstat wght\tv\tJ\n" );
00273 }
00274
00275 else if( nMatch("LINE",chCard) )
00276 {
00277
00278 strcpy( punch.chPunch[punch.npunch], "H2ln" );
00279 fprintf( punch.ipPnunit[punch.npunch],
00280 "#H2 line\tEhi\tVhi\tJhi\tElo\tVlo\tJlo\twl(mic)\twl(lab)\tlog flux\tI/Inorm\tExcit(hi, K)\tg_u h\\nu * Aul\n" );
00281
00282
00283
00284
00285
00286 thresh_punline_h2 = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
00287 if( lgEOL )
00288 {
00289
00290 thresh_punline_h2 = 1e-4f;
00291 }
00292
00293
00294 if( thresh_punline_h2 < 0. )
00295 {
00296 thresh_punline_h2 = (float)pow(10.f,thresh_punline_h2);
00297 }
00298
00299
00300
00301
00302
00303 if( nMatch( "ELEC" , chCard ) )
00304 {
00305 if( nMatch(" ALL",chCard) )
00306 {
00307
00308
00309
00310 h2.nElecLevelOutput = -1;
00311 }
00312 else if( nMatch("GROU",chCard) )
00313 {
00314
00315 h2.nElecLevelOutput = 1;
00316 }
00317 else
00318 {
00319 h2.nElecLevelOutput = (int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
00320 if( lgEOL )
00321 {
00322
00323 h2.nElecLevelOutput = 1;
00324 }
00325 }
00326 }
00327 }
00328
00329 else if( nMatch(" PDR",chCard) )
00330 {
00331
00332 strcpy( punch.chPunch[punch.npunch], "H2pd" );
00333 fprintf( punch.ipPnunit[punch.npunch], "#H2 creation, destruction. \n" );
00334 }
00335 else if( nMatch("POPU",chCard) )
00336 {
00337
00338 strcpy( punch.chPunch[punch.npunch], "H2po" );
00339
00340
00341
00342
00343
00344 punch.punarg[punch.npunch][0] = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
00345
00346
00347 punch.punarg[punch.npunch][1] = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
00348
00349 if( nMatch( "ZONE" , chCard ) )
00350 {
00351
00352 punch.punarg[punch.npunch][2] = 0;
00353 fprintf( punch.ipPnunit[punch.npunch], "#depth\torth\tpar\te=1 rel pop\te=2 rel pop\tv=0 rel pops\n" );
00354 }
00355 else
00356 {
00357
00358
00359
00360 if( nMatch( "MATR" , chCard ) )
00361 {
00362
00363 punch.punarg[punch.npunch][2] = 1;
00364 fprintf( punch.ipPnunit[punch.npunch], "#vib\trot\tpops\n" );
00365 }
00366 else
00367 {
00368
00369 punch.punarg[punch.npunch][2] = -1;
00370 fprintf( punch.ipPnunit[punch.npunch], "#vib\trot\ts\tenergy(wn)\tpops/H2\told/H2\tpops/g\tdep coef\tFin(Col)\tFout(col)\tRCout\tRRout\tRCin\tRRin\n" );
00371 }
00372 }
00373 }
00374
00375 else if( nMatch("RATE",chCard) )
00376 {
00377
00378 strcpy( punch.chPunch[punch.npunch], "H2ra" );
00379 fprintf( punch.ipPnunit[punch.npunch],
00380 "#depth\tN(H2)\tN(H2)/u(H2)\tA_V(star)\tn(Eval)"
00381 "\tH2/Htot\trenorm\tfrm grn\tfrmH-\tdstTH85\tBD96\tELWERT\tBigH2\telec->H2g\telec->H2s"
00382 "\tG(TH85)\tG(DB96)\tCR\tEleclife\tShield(BD96)\tShield(H2)\tBigh2/G0(spc)\ttot dest"
00383 "\tHeatH2Dish_TH85\tHeatH2Dexc_TH85\tHeatH2Dish_BigH2\tHeatH2Dexc_BigH2\thtot\n" );
00384 }
00385 else if( nMatch("SOLO",chCard) )
00386 {
00387
00388 strcpy( punch.chPunch[punch.npunch], "H2so" );
00389 fprintf( punch.ipPnunit[punch.npunch],
00390 "#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" );
00391 }
00392 else if( nMatch("SPEC",chCard) )
00393 {
00394
00395 strcpy( punch.chPunch[punch.npunch], "H2sp" );
00396 fprintf( punch.ipPnunit[punch.npunch],
00397 "#depth\tspecial\n" );
00398 }
00399 else if( nMatch("TEMP",chCard) )
00400 {
00401
00402 strcpy( punch.chPunch[punch.npunch], "H2te" );
00403 fprintf( punch.ipPnunit[punch.npunch],
00404 "#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");
00405 }
00406 else if( nMatch("THER",chCard) )
00407 {
00408
00409 strcpy( punch.chPunch[punch.npunch], "H2th" );
00410 fprintf( punch.ipPnunit[punch.npunch],
00411 "#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");
00412 }
00413 else
00414 {
00415 fprintf( ioQQQ,
00416 " There must be a second key; they are RATE, LINE, COOL, COLUMN, _PDR, SOLOmon, TEMP, and POPUlations\n" );
00417 puts( "[Stop in ParsePunch]" );
00418 cdEXIT(EXIT_FAILURE);
00419 }
00420 DEBUG_EXIT( "H2_ParsePunch()" );
00421 return;
00422 }
00423
00424
00425
00426 void H2_Prt_Zone(void)
00427 {
00428 int iElecHi , iVibHi;
00429
00430
00431 if( !h2.lgH2ON || !h2.nCallH2_this_zone )
00432 return;
00433
00434 DEBUG_ENTRY( "H2_Prt_Zone()" );
00435
00436 fprintf( ioQQQ, " H2 density ");
00437 fprintf(ioQQQ,PrintEfmt("%9.2e", hmi.H2_total));
00438
00439 fprintf( ioQQQ, " orth/par");
00440 fprintf(ioQQQ,PrintEfmt("%9.2e", h2.ortho_density / SDIV( h2.para_density )));
00441
00442 iElecHi = 0;
00443 iVibHi = 0;
00444 fprintf( ioQQQ, " v0 J=0,3");
00445 fprintf(ioQQQ,PrintEfmt("%9.2e", H2_populations[iElecHi][iVibHi][0] / hmi.H2_total));
00446 fprintf(ioQQQ,PrintEfmt("%9.2e", H2_populations[iElecHi][iVibHi][1] / hmi.H2_total));
00447 fprintf(ioQQQ,PrintEfmt("%9.2e", H2_populations[iElecHi][iVibHi][2] / hmi.H2_total));
00448 fprintf(ioQQQ,PrintEfmt("%9.2e", H2_populations[iElecHi][iVibHi][3] / hmi.H2_total));
00449
00450 fprintf( ioQQQ, " TOTv=0,3");
00451 fprintf(ioQQQ,PrintEfmt("%9.2e", pops_per_vib[iElecHi][0] / hmi.H2_total));
00452 fprintf(ioQQQ,PrintEfmt("%9.2e", pops_per_vib[iElecHi][1] / hmi.H2_total));
00453 fprintf(ioQQQ,PrintEfmt("%9.2e", pops_per_vib[iElecHi][2] / hmi.H2_total));
00454 fprintf(ioQQQ,PrintEfmt("%9.2e", pops_per_vib[iElecHi][3] / hmi.H2_total));
00455 fprintf( ioQQQ, "\n");
00456
00457 DEBUG_EXIT( "H2_Prt_Zone()" );
00458 return;
00459 }
00460
00461
00462 void H2_Prt_column_density(
00463
00464
00465 FILE *ioMEAN )
00466
00467 {
00468 int iVibHi;
00469
00470
00471 if( !h2.lgH2ON || !h2.nCallH2_this_zone )
00472 return;
00473
00474 DEBUG_ENTRY( "H2_Prt_column_density()" );
00475
00476 fprintf( ioMEAN, " H2 total ");
00477 fprintf(ioMEAN,"%7.3f", log10(SDIV(colden.colden[ipCOL_H2g]+colden.colden[ipCOL_H2s])));
00478
00479 fprintf( ioMEAN, " H2 ortho ");
00480 fprintf(ioMEAN,"%7.3f", log10(SDIV(h2.ortho_colden)));
00481
00482 fprintf( ioMEAN, " para");
00483 fprintf(ioMEAN,"%7.3f", log10(SDIV(h2.para_colden)));
00484
00485 iVibHi = 0;
00486 fprintf( ioMEAN, " v0 J=0,3");
00487 fprintf(ioMEAN,"%7.3f", log10(SDIV(H2_X_colden[iVibHi][0])));
00488 fprintf(ioMEAN,"%7.3f", log10(SDIV(H2_X_colden[iVibHi][1])));
00489 fprintf(ioMEAN,"%7.3f", log10(SDIV(H2_X_colden[iVibHi][2])));
00490 fprintf(ioMEAN,"%7.3f", log10(SDIV(H2_X_colden[iVibHi][3])));
00491
00492 # if 0
00493 fprintf( ioMEAN, " v=0,3");
00494 fprintf(ioMEAN,PrintEfmt("%9.2e", pops_per_vib[iElecHi][0] / hmi.H2_total));
00495 fprintf(ioMEAN,PrintEfmt("%9.2e", pops_per_vib[iElecHi][1] / hmi.H2_total));
00496 fprintf(ioMEAN,PrintEfmt("%9.2e", pops_per_vib[iElecHi][2] / hmi.H2_total));
00497 fprintf(ioMEAN,PrintEfmt("%9.2e", pops_per_vib[iElecHi][3] / hmi.H2_total));
00498 fprintf( ioMEAN, "\n");
00499 # endif
00500
00501 DEBUG_EXIT( "H2_Prt_column_density()" );
00502 return;
00503 }
00504
00505
00506
00507 void H2_ReadCollRates( long int nColl )
00508 {
00509
00510 char cdDATAFILE[N_X_COLLIDER][FILENAME_PATH_LENGTH_2] =
00511 {
00512 "H2_coll_H.dat" ,
00513 "H2_coll_He.dat" ,
00514 "H2_coll_H2ortho.dat" ,
00515 "H2_coll_H2para.dat",
00516 "H2_coll_Hp.dat" ,
00517 "H2_coll_He_Stancil.dat"
00518 };
00519 FILE *ioDATA;
00520 char chLine[FILENAME_PATH_LENGTH_2] ,
00521 chFilename[FILENAME_PATH_LENGTH_2];
00522 long int i, n1, n2, n3;
00523 long int iVibHi , iVibLo , iRotHi , iRotLo;
00524 bool lgEOL;
00525
00526 DEBUG_ENTRY( "H2_ReadCollRates()" );
00527
00528
00529
00530 if( lgDataPathSet == true )
00531 {
00532
00533 strcpy( chFilename , chDataPath );
00534 strcat( chFilename , cdDATAFILE[nColl] );
00535 }
00536 else
00537 {
00538
00539 strcpy( chFilename , cdDATAFILE[nColl] );
00540 }
00541
00542 if( nColl==5 )
00543 {
00544
00545
00546
00547
00548
00549
00550 if( H2_He_coll_init( chFilename ) != 51018 )
00551 {
00552 fprintf(ioQQQ,"the H2 - He collision data file H2_coll_He_Stancil.dat does not have the correct magic number.\n");
00553 fprintf(ioQQQ,"I expected %i\n", 51018 );
00554
00555 puts( "[Stop in H2_ReadCollRates]" );
00556 cdEXIT(EXIT_FAILURE);
00557 }
00558
00559 DEBUG_EXIT( "H2_ReadCollRates()" );
00560 return;
00561 }
00562
00563
00564 if( ( ioDATA = fopen( chFilename , "r" ) ) == NULL )
00565 {
00566 fprintf( ioQQQ, " H2_ReadCollRates could not open %s\n", cdDATAFILE[nColl] );
00567 if( lgDataPathSet == true )
00568 fprintf( ioQQQ, " even tried path\n" );
00569
00570 if( lgDataPathSet == true )
00571 {
00572 fprintf( ioQQQ, " H2_ReadCollRates could not open %s\n", cdDATAFILE[nColl]);
00573 fprintf( ioQQQ, " path is ==%s==\n",chDataPath );
00574 fprintf( ioQQQ, " final path is ==%s==\n",chFilename );
00575 }
00576
00577 path_not_set();
00578
00579 puts( "[Stop in H2_ReadCollRates]" );
00580 cdEXIT(EXIT_FAILURE);
00581 }
00582
00583
00584 if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00585 {
00586 fprintf( ioQQQ, " H2_ReadCollRates could not read first line of %s\n", cdDATAFILE[nColl]);
00587 puts( "[Stop in H2_ReadCollRates]" );
00588 cdEXIT(EXIT_FAILURE);
00589 }
00590 i = 1;
00591
00592 n1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00593 n2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00594 n3 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00595
00596
00597
00598 if( ( n1 != 2 ) || ( n2 != 4 ) || ( n3 != 29 ) )
00599 {
00600 fprintf( ioQQQ,
00601 " H2_ReadCollRates: the version of %s is not the current version.\n", cdDATAFILE[nColl] );
00602 fprintf( ioQQQ,
00603 " I expected to find the number 2 4 29 and got %li %li %li instead.\n" ,
00604 n1 , n2 , n3 );
00605 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
00606 puts( "[Stop in H2_ReadCollRates]" );
00607 cdEXIT(EXIT_FAILURE);
00608 }
00609
00610
00611 if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00612 BadRead();
00613
00614 while( chLine[0]=='#' )
00615 {
00616 if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00617 BadRead();
00618 }
00619 iRotHi = 1;
00620 while( iRotHi >= 0 )
00621 {
00622 double a[3];
00623 sscanf(chLine,"%li\t%li\t%li\t%li\t%le\t%le\t%le",
00624 &iVibHi ,&iRotHi , &iVibLo , &iRotLo , &a[0],&a[1],&a[2] );
00625
00626 if( iRotHi < 0 )
00627 continue;
00628
00629
00630 ASSERT( iVibHi <= VIB_COLLID &&
00631 iVibLo <= VIB_COLLID &&
00632 iRotHi <= h2.nRot_hi[0][iVibHi] &&
00633 iRotLo <= h2.nRot_hi[0][iVibLo]);
00634
00635
00636 if( !( (iVibHi == iVibLo) && (iRotHi == iRotLo )) )
00637 {
00638
00639 ASSERT( (energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo] ) > 0. );
00640 for( i=0; i<3; ++i )
00641 {
00642 CollRateFit[nColl][iVibHi][iRotHi][iVibLo][iRotLo][i] = (float)a[i];
00643 }
00644
00645
00646
00647
00648 }
00649
00650 if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00651 BadRead();
00652 while( chLine[0]=='#' )
00653 {
00654 if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00655 BadRead();
00656 }
00657 }
00658 fclose( ioDATA );
00659
00660 DEBUG_EXIT( "H2_ReadCollRates()" );
00661
00662 return;
00663 }
00664
00665
00666
00667 void H2_ReadTransprob( long int nelec )
00668 {
00669 char cdDATAFILE[N_H2_ELEC][FILENAME_PATH_LENGTH_2] =
00670 {
00671 "H2_transprob_X.dat" ,
00672 "H2_transprob_B.dat" ,
00673 "H2_transprob_C_plus.dat" ,
00674 "H2_transprob_C_minus.dat" ,
00675 "H2_transprob_B_primed.dat" ,
00676 "H2_transprob_D_plus.dat" ,
00677 "H2_transprob_D_minus.dat"
00678 };
00679 FILE *ioDATA;
00680 char chLine[FILENAME_PATH_LENGTH_2] ,
00681 chFilename[FILENAME_PATH_LENGTH_2];
00682 long int i, n1, n2, n3;
00683 long int iVibHi , iVibLo , iRotHi , iRotLo , iElecHi , iElecLo;
00684 bool lgEOL;
00685
00686 DEBUG_ENTRY( "H2_ReadTransprob()" );
00687
00688
00689
00690 if( lgDataPathSet == true )
00691 {
00692
00693 strcpy( chFilename , chDataPath );
00694 strcat( chFilename , cdDATAFILE[nelec] );
00695 }
00696 else
00697 {
00698
00699 strcpy( chFilename , cdDATAFILE[nelec] );
00700 }
00701
00702
00703 if( ( ioDATA = fopen( chFilename , "r" ) ) == NULL )
00704 {
00705 fprintf( ioQQQ, " H2_ReadTransprob could not open %s\n", cdDATAFILE[nelec] );
00706 if( lgDataPathSet == true )
00707 fprintf( ioQQQ, " even tried path\n" );
00708
00709 if( lgDataPathSet == true )
00710 {
00711 fprintf( ioQQQ, " H2_ReadTransprob could not open %s\n", cdDATAFILE[nelec]);
00712 fprintf( ioQQQ, " path is ==%s==\n",chDataPath );
00713 fprintf( ioQQQ, " final path is ==%s==\n",chFilename );
00714 }
00715
00716 puts( "[Stop in H2_ReadTransprob]" );
00717 cdEXIT(EXIT_FAILURE);
00718 }
00719
00720
00721 if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00722 {
00723 fprintf( ioQQQ, " H2_ReadTransprob could not read first line of %s\n", cdDATAFILE[nelec]);
00724 puts( "[Stop in H2_ReadTransprob]" );
00725 cdEXIT(EXIT_FAILURE);
00726 }
00727 i = 1;
00728
00729 n1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00730 n2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00731 n3 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00732
00733
00734
00735 if( ( n1 != 2 ) || ( n2 != 4 ) || ( n3 != 29 ) )
00736 {
00737 fprintf( ioQQQ,
00738 " H2_ReadTransprob: the version of %s is not the current version.\n", cdDATAFILE[nelec] );
00739 fprintf( ioQQQ,
00740 " I expected to find the number 2 4 29 and got %li %li %li instead.\n" ,
00741 n1 , n2 , n3 );
00742 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
00743 puts( "[Stop in H2_ReadTransprob]" );
00744 cdEXIT(EXIT_FAILURE);
00745 }
00746
00747
00748 if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00749 BadRead();
00750
00751 while( chLine[0]=='#' )
00752 {
00753 if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00754 BadRead();
00755 }
00756 iVibHi = 1;
00757 while( iVibHi >= 0 )
00758 {
00759 double Aul;
00760 sscanf(chLine,"%li\t%li\t%li\t%li\t%li\t%li\t%le",
00761 &iElecHi , &iVibHi ,&iRotHi , &iElecLo , &iVibLo , &iRotLo , &Aul );
00762 ASSERT( iElecHi == nelec );
00763
00764 if( iVibHi < 0 )
00765 continue;
00766
00767
00768 if( iVibHi <= h2.nVib_hi[iElecHi] &&
00769 iVibLo <= h2.nVib_hi[iElecLo] &&
00770 iRotHi <= h2.nRot_hi[iElecHi][iVibHi] &&
00771 iRotLo <= h2.nRot_hi[iElecLo][iVibLo])
00772 {
00773 double ener = energy_wn[iElecHi][iVibHi][iRotHi] - energy_wn[iElecLo][iVibLo][iRotLo];
00774
00775
00776 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Aul = (float)Aul;
00777
00778 lgH2_line_exists[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] = true;
00779
00780 if( ener <= 0. )
00781 {
00782 fprintf(ioQQQ,"negative energy H2 transition\t%li\t%li\t%li\t%li\t%.2e\t%.2e\n",
00783 iVibHi,iVibLo,iRotHi,iRotLo,Aul,
00784 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyWN);
00785 ShowMe();
00786 }
00787 }
00788 # if 0
00789
00790 else
00791 {
00792 fprintf(ioQQQ,"no\t%li\t%li\t%li\t%li\t%.2e\n",
00793 iVibHi,iVibLo,iRotHi,iRotLo,Aul);
00794 }
00795 # endif
00796
00797 if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00798 BadRead();
00799 while( chLine[0]=='#' )
00800 {
00801 if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00802 BadRead();
00803 }
00804 }
00805 fclose( ioDATA );
00806 DEBUG_EXIT( "H2_ReadTransprob()" );
00807 return;
00808 }
00809
00810 #if 0
00811
00812 void H2_Read_Cosmicray_distribution(void)
00813 {
00814
00815 FILE *ioDATA;
00816 char chLine[FILENAME_PATH_LENGTH_2] ,
00817 chFilename[FILENAME_PATH_LENGTH_2];
00818 long int i, n1, n2, n3, iVib , iRot;
00819 long neut_frac;
00820 bool lgEOL;
00821
00822 DEBUG_ENTRY( "H2_Read_Cosmicray_distribution()" );
00823
00824
00825
00826 if( lgDataPathSet == true )
00827 {
00828
00829 strcpy( chFilename , chDataPath );
00830 strcat( chFilename , "H2_CosmicRay_collision.dat" );
00831 }
00832 else
00833 {
00834
00835 strcpy( chFilename , "H2_CosmicRay_collision.dat" );
00836 }
00837
00838
00839 if( ( ioDATA = fopen( chFilename , "r" ) ) == NULL )
00840 {
00841 fprintf( ioQQQ, " H2_Read_CosmicRay_distribution could not open %s\n", chFilename );
00842 if( lgDataPathSet == true )
00843 fprintf( ioQQQ, " even tried path\n" );
00844
00845 if( lgDataPathSet == true )
00846 {
00847 fprintf( ioQQQ, " H2_Read_CosmicRay_distribution could not open %s\n", chFilename );
00848 fprintf( ioQQQ, " path is ==%s==\n",chDataPath );
00849 fprintf( ioQQQ, " final path is ==%s==\n",chFilename );
00850 }
00851
00852 puts( "[Stop in H2_Read_Cosmicray_distribution]" );
00853 cdEXIT(EXIT_FAILURE);
00854 }
00855
00856
00857 if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00858 {
00859 fprintf( ioQQQ, " H2_Read_Cosmicray_distribution could not read first line of %s\n", "H2_Cosmic_collision.dat");
00860 puts( "[Stop in H2_Read_Cosmicray_distribution]" );
00861 cdEXIT(EXIT_FAILURE);
00862 }
00863
00864 i = 1;
00865
00866 n1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00867 n2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00868 n3 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00869
00870
00871
00872 if( ( n1 != 1 ) || ( n2 != 21 ) || ( n3 != 3 ) )
00873 {
00874 fprintf( ioQQQ,
00875 " H2_Read_Cosmicray_distribution: the version of %s is not the current version.\n", "H2_Cosmic_collision.dat" );
00876 fprintf( ioQQQ,
00877 " I expected to find the number 1 21 3 and got %li %li %li instead.\n" ,
00878 n1 , n2 , n3 );
00879 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
00880 puts( "[Stop in H2_Read_Cosmicray_distribution]" );
00881 cdEXIT(EXIT_FAILURE);
00882 }
00883
00884
00885 if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00886 BadRead();
00887
00888 while( chLine[0]=='#' )
00889 {
00890 if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00891 BadRead();
00892 }
00893
00894 iRot = 1;
00895 iVib = 1;
00896 neut_frac = 0;
00897 while( iVib >= 0 )
00898 {
00899 long int j_minus_ji;
00900 double a[10];
00901
00902 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",
00903 &iVib ,&j_minus_ji , &a[0],&a[1],&a[2],&a[3],&a[4],&a[5],&a[6],&a[7],&a[8],&a[9]
00904 );
00905
00906 if( iVib < 0 )
00907 continue;
00908
00909
00910
00911 ASSERT( iVib < CR_VIB );
00912 ASSERT( j_minus_ji == -2 || j_minus_ji == +2 || j_minus_ji == 0 );
00913 ASSERT( neut_frac < CR_X );
00914
00915
00916 j_minus_ji = 1 + j_minus_ji/2;
00917 ASSERT( j_minus_ji>=0 && j_minus_ji<=2 );
00918
00919
00920 for( iRot=0; iRot<CR_J; ++iRot )
00921 {
00922 cr_rate[neut_frac][iVib][iRot][j_minus_ji] = (float)a[iRot];
00923 }
00924 if( mole.lgH2_NOISECOSMIC )
00925 {
00926 float r;
00927 r = (float)RandGauss( mole.xMeanNoise , mole.xSTDNoise );
00928
00929 for( iRot=0; iRot<CR_J; ++iRot )
00930 {
00931 cr_rate[neut_frac][iVib][iRot][j_minus_ji] *= (float)pow(10.,(double)r);
00932 }
00933 }
00934
00935 if( CR_PRINT )
00936 {
00937 fprintf(ioQQQ,"cr rate\t%li\t%li", iVib , j_minus_ji );
00938 for( iRot=0; iRot<CR_J; ++iRot )
00939 {
00940 fprintf(ioQQQ,"\t%.3e", cr_rate[neut_frac][iVib][iRot][j_minus_ji] );
00941 }
00942 fprintf(ioQQQ,"\n" );
00943 }
00944
00945
00946 if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00947 BadRead();
00948 while( chLine[0]=='#' )
00949 {
00950 if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00951 BadRead();
00952 }
00953 }
00954 fclose( ioDATA );
00955
00956 DEBUG_EXIT( "H2_Read_Cosmicray_distribution()" );
00957
00958 return;
00959 }
00960 #endif
00961
00962
00963 void H2_ReadEnergies( long int nelec )
00964 {
00965 char cdDATAFILE[N_H2_ELEC][FILENAME_PATH_LENGTH_2] =
00966 {
00967 "H2_energy_X.dat" ,
00968 "H2_energy_B.dat" ,
00969 "H2_energy_C_plus.dat" ,
00970 "H2_energy_C_minus.dat" ,
00971 "H2_energy_B_primed.dat" ,
00972 "H2_energy_D_plus.dat" ,
00973 "H2_energy_D_minus.dat"
00974 };
00975 FILE *ioDATA;
00976 char chLine[FILENAME_PATH_LENGTH_2] ,
00977 chFilename[FILENAME_PATH_LENGTH_2];
00978 long int i, n1, n2, n3, iVib , iRot;
00979 bool lgEOL;
00980
00981 DEBUG_ENTRY( "H2_ReadEnergies()" );
00982
00983
00984
00985 if( lgDataPathSet == true )
00986 {
00987
00988 strcpy( chFilename , chDataPath );
00989 strcat( chFilename , cdDATAFILE[nelec] );
00990 }
00991 else
00992 {
00993
00994 strcpy( chFilename , cdDATAFILE[nelec] );
00995 }
00996
00997
00998 if( ( ioDATA = fopen( chFilename , "r" ) ) == NULL )
00999 {
01000 fprintf( ioQQQ, " H2_ReadEnergies could not open %s\n", cdDATAFILE[nelec] );
01001 if( lgDataPathSet == true )
01002 fprintf( ioQQQ, " even tried path\n" );
01003
01004 if( lgDataPathSet == true )
01005 {
01006 fprintf( ioQQQ, " H2_ReadEnergies could not open %s\n", cdDATAFILE[nelec]);
01007 fprintf( ioQQQ, " path is ==%s==\n",chDataPath );
01008 fprintf( ioQQQ, " final path is ==%s==\n",chFilename );
01009 }
01010
01011 puts( "[Stop in H2_ReadEnergies]" );
01012 cdEXIT(EXIT_FAILURE);
01013 }
01014
01015
01016 if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
01017 {
01018 fprintf( ioQQQ, " H2_ReadEnergies could not read first line of %s\n", cdDATAFILE[nelec]);
01019 puts( "[Stop in H2_ReadEnergies]" );
01020 cdEXIT(EXIT_FAILURE);
01021 }
01022 i = 1;
01023
01024 n1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
01025 n2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
01026 n3 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
01027
01028
01029
01030 if( ( n1 != 2 ) || ( n2 != 4 ) || ( n3 != 29 ) )
01031 {
01032 fprintf( ioQQQ,
01033 " H2_ReadEnergies: the version of %s is not the current version.\n", cdDATAFILE[nelec] );
01034 fprintf( ioQQQ,
01035 " I expected to find the number 2 4 29 and got %li %li %li instead.\n" ,
01036 n1 , n2 , n3 );
01037 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
01038 puts( "[Stop in H2_ReadEnergies]" );
01039 cdEXIT(EXIT_FAILURE);
01040 }
01041
01042
01043 if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
01044 BadRead();
01045
01046 while( chLine[0]=='#' )
01047 {
01048 if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
01049 BadRead();
01050 }
01051
01052
01053 nLevels_per_elec[nelec] = 0;
01054
01055 for( iVib=0; iVib<=h2.nVib_hi[nelec]; ++iVib )
01056 {
01057 for( iRot=h2.Jlowest[nelec]; iRot<=h2.nRot_hi[nelec][iVib]; ++iRot )
01058 {
01059 i = 1;
01060 sscanf(chLine,"%li\t%li\t%le", &n1 , &n2 , &energy_wn[nelec][iVib][iRot] );
01061 ASSERT( n1 == iVib );
01062 ASSERT( n2 == iRot );
01063 # if 0
01064
01065 if( nelec == 0 )
01066 {
01067
01068
01069 energy_wn[0][iVib][iRot] = -( energy_wn[0][iVib][iRot]- 3.6118114E+04 );
01070 }
01071 # endif
01072 ASSERT( energy_wn[nelec][iVib][iRot]> 0. || (nelec==0 && iVib==0 && iRot==0 ) );
01073
01074 ++nLevels_per_elec[nelec];
01075
01076
01077 if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
01078 BadRead();
01079 while( chLine[0]=='#' )
01080 {
01081 if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
01082 BadRead();
01083 }
01084 }
01085 }
01086 fclose( ioDATA );
01087 DEBUG_EXIT( "H2_ReadEnergies()" );
01088 return;
01089 }
01090
01091
01092 void H2_ReadDissprob( long int nelec )
01093 {
01094 char cdDATAFILE[N_H2_ELEC][FILENAME_PATH_LENGTH_2] =
01095 {
01096 "H2_dissprob_X.dat" ,
01097 "H2_dissprob_B.dat" ,
01098 "H2_dissprob_C_plus.dat" ,
01099 "H2_dissprob_C_minus.dat" ,
01100 "H2_dissprob_B_primed.dat" ,
01101 "H2_dissprob_D_plus.dat" ,
01102 "H2_dissprob_D_minus.dat"
01103 };
01104 FILE *ioDATA;
01105 char chLine[FILENAME_PATH_LENGTH_2] ,
01106 chFilename[FILENAME_PATH_LENGTH_2];
01107 long int i, n1, n2, n3, iVib , iRot;
01108 bool lgEOL;
01109
01110 DEBUG_ENTRY( "H2_ReadDissprob()" );
01111
01112 ASSERT( nelec > 0 );
01113
01114
01115
01116 if( lgDataPathSet == true )
01117 {
01118
01119 strcpy( chFilename , chDataPath );
01120 strcat( chFilename , cdDATAFILE[nelec] );
01121 }
01122 else
01123 {
01124
01125 strcpy( chFilename , cdDATAFILE[nelec] );
01126 }
01127
01128
01129 if( ( ioDATA = fopen( chFilename , "r" ) ) == NULL )
01130 {
01131 fprintf( ioQQQ, " H2_ReadDissprob could not open %s\n", cdDATAFILE[nelec] );
01132 if( lgDataPathSet == true )
01133 fprintf( ioQQQ, " even tried path\n" );
01134
01135 if( lgDataPathSet == true )
01136 {
01137 fprintf( ioQQQ, " H2_ReadDissprob could not open %s\n", cdDATAFILE[nelec]);
01138 fprintf( ioQQQ, " path is ==%s==\n",chDataPath );
01139 fprintf( ioQQQ, " final path is ==%s==\n",chFilename );
01140 }
01141
01142 puts( "[Stop in H2_ReadDissprob]" );
01143 cdEXIT(EXIT_FAILURE);
01144 }
01145
01146
01147 if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
01148 {
01149 fprintf( ioQQQ, " H2_ReadDissprob could not read first line of %s\n", cdDATAFILE[nelec]);
01150 puts( "[Stop in H2_ReadDissprob]" );
01151 cdEXIT(EXIT_FAILURE);
01152 }
01153 i = 1;
01154
01155 n1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
01156 n2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
01157 n3 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
01158
01159
01160
01161 if( ( n1 != 3 ) || ( n2 != 2 ) || ( n3 != 11 ) )
01162 {
01163 fprintf( ioQQQ,
01164 " H2_ReadDissprob: the version of %s is not the current version.\n", cdDATAFILE[nelec] );
01165 fprintf( ioQQQ,
01166 " I expected to find the number 3 2 11 and got %li %li %li instead.\n" ,
01167 n1 , n2 , n3 );
01168 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
01169 puts( "[Stop in H2_ReadDissprob]" );
01170 cdEXIT(EXIT_FAILURE);
01171 }
01172
01173
01174 if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
01175 BadRead();
01176
01177 while( chLine[0]=='#' )
01178 {
01179 if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
01180 BadRead();
01181 }
01182
01183 for( iVib=0; iVib<=h2.nVib_hi[nelec]; ++iVib )
01184 {
01185 for( iRot=h2.Jlowest[nelec]; iRot<=h2.nRot_hi[nelec][iVib]; ++iRot )
01186 {
01187 double a, b;
01188 i = 1;
01189 sscanf(chLine,"%li\t%li\t%le\t%le",
01190 &n1 , &n2 ,
01191
01192 &a ,
01193
01194 &b);
01195
01196
01197 ASSERT( n1 == iVib );
01198 ASSERT( n2 == iRot );
01199
01200
01201 H2_dissprob[nelec][iVib][iRot] = (float)a;
01202
01203 H2_disske[nelec][iVib][iRot] = (float)b;
01204
01205
01206 if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
01207 BadRead();
01208 while( chLine[0]=='#' )
01209 {
01210 if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
01211 BadRead();
01212 }
01213 }
01214 }
01215 fclose( ioDATA );
01216 }
01217
01218
01219
01220 void H2_Read_hminus_distribution(void)
01221 {
01222 FILE *ioDATA;
01223 char chLine[FILENAME_PATH_LENGTH_2] ,
01224 chFilename[FILENAME_PATH_LENGTH_2];
01225 long int i, n1, n2, n3, iVib , iRot;
01226 bool lgEOL;
01227 double sumrate[nTE_HMINUS];
01228
01229 # define H2HMINUS_PRT false
01230
01231
01232
01233 if( lgDataPathSet == true )
01234 {
01235
01236 strcpy( chFilename , chDataPath );
01237 strcat( chFilename , "H2_hminus_deposit.dat" );
01238 }
01239 else
01240 {
01241
01242 strcpy( chFilename , "H2_hminus_deposit.dat" );
01243 }
01244
01245
01246 if( ( ioDATA = fopen( chFilename , "r" ) ) == NULL )
01247 {
01248 fprintf( ioQQQ, " H2_Read_hminus_distribution could not open %s\n", "H2_hminus_deposit.dat" );
01249 if( lgDataPathSet == true )
01250 fprintf( ioQQQ, " even tried path\n" );
01251
01252 if( lgDataPathSet == true )
01253 {
01254 fprintf( ioQQQ, " H2_Read_hminus_distribution could not open %s\n", "H2_hminus_deposit.dat");
01255 fprintf( ioQQQ, " path is ==%s==\n",chDataPath );
01256 fprintf( ioQQQ, " final path is ==%s==\n",chFilename );
01257 }
01258
01259 puts( "[Stop in H2_Read_hminus_distribution]" );
01260 cdEXIT(EXIT_FAILURE);
01261 }
01262
01263
01264 if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
01265 {
01266 fprintf( ioQQQ, " H2_Read_hminus_distribution could not read first line of %s\n", "H2_hminus_deposit.dat");
01267 puts( "[Stop in H2_Read_hminus_distribution]" );
01268 cdEXIT(EXIT_FAILURE);
01269 }
01270
01271 i = 1;
01272
01273 n1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
01274 n2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
01275 n3 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
01276
01277
01278
01279 if( ( n1 != 2 ) || ( n2 != 10 ) || ( n3 != 17 ) )
01280 {
01281 fprintf( ioQQQ,
01282 " H2_Read_hminus_distribution: the version of %s is not the current version.\n", "H2_hminus_deposit.dat" );
01283 fprintf( ioQQQ,
01284 " I expected to find the number 2 10 17 and got %li %li %li instead.\n" ,
01285 n1 , n2 , n3 );
01286 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
01287 puts( "[Stop in H2_Read_hminus_distribution]" );
01288 cdEXIT(EXIT_FAILURE);
01289 }
01290
01291
01292 if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
01293 BadRead();
01294
01295 while( chLine[0]=='#' )
01296 {
01297 if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
01298 BadRead();
01299 }
01300
01301
01302 for(i=0; i<nTE_HMINUS; ++i )
01303 {
01304 H2_te_hminus[i] = (float)log10(H2_te_hminus[i]);
01305 sumrate[i] = 0.;
01306 }
01307
01308 iRot = 1;
01309 iVib = 1;
01310 while( iVib >= 0 )
01311 {
01312
01313
01314 double a[nTE_HMINUS] , ener;
01315 sscanf(chLine,"%li\t%li\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf",
01316 &iVib ,&iRot , &ener, &a[0],&a[1],&a[2] , &a[3],&a[4],&a[5] ,&a[6]
01317 );
01318
01319 if( iVib < 0 )
01320 continue;
01321
01322
01323 ASSERT( iVib <= h2.nVib_hi[0] &&
01324 iRot <= h2.nRot_hi[0][iVib] );
01325
01326 if( H2HMINUS_PRT )
01327 fprintf(ioQQQ,"hminusss\t%li\t%li", iVib , iRot );
01328 for( i=0; i<nTE_HMINUS; ++i )
01329 {
01330 H2_X_hminus_formation_distribution[i][iVib][iRot] = (float)pow(10.,-a[i]);
01331 sumrate[i] += H2_X_hminus_formation_distribution[i][iVib][iRot];
01332 if( H2HMINUS_PRT )
01333 fprintf(ioQQQ,"\t%.3e", H2_X_hminus_formation_distribution[i][iVib][iRot] );
01334 }
01335 if( H2HMINUS_PRT )
01336 fprintf(ioQQQ,"\n" );
01337
01338 if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
01339 BadRead();
01340 while( chLine[0]=='#' )
01341 {
01342 if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
01343 BadRead();
01344 }
01345 }
01346 fclose( ioDATA );
01347
01348 if( H2HMINUS_PRT )
01349 {
01350
01351 fprintf(ioQQQ," total H- formation rate ");
01352
01353 for(i=0; i<nTE_HMINUS; ++i )
01354 {
01355 fprintf(ioQQQ,"\t%.3e" , sumrate[i]);
01356 }
01357 fprintf(ioQQQ,"\n" );
01358 }
01359
01360
01361 for( iVib=0; iVib<=h2.nVib_hi[0]; ++iVib )
01362 {
01363 for( iRot=h2.Jlowest[0]; iRot<=h2.nRot_hi[0][iVib]; ++iRot )
01364 {
01365 for(i=0; i<nTE_HMINUS; ++i )
01366 {
01367 H2_X_hminus_formation_distribution[i][iVib][iRot] /= (float)sumrate[i];
01368 }
01369 }
01370 }
01371
01372 if( H2HMINUS_PRT )
01373 {
01374
01375 fprintf(ioQQQ," H- distribution function ");
01376 for( iVib=0; iVib<=h2.nVib_hi[0]; ++iVib )
01377 {
01378 for( iRot=h2.Jlowest[0]; iRot<=h2.nRot_hi[0][iVib]; ++iRot )
01379 {
01380 fprintf(ioQQQ,"%li\t%li", iVib , iRot );
01381 for(i=0; i<nTE_HMINUS; ++i )
01382 {
01383 fprintf(ioQQQ,"\t%.3e", H2_X_hminus_formation_distribution[i][iVib][iRot] );
01384 }
01385 fprintf(ioQQQ,"\n" );
01386 }
01387 }
01388 }
01389
01390 DEBUG_EXIT( "H2_ReadDissprob()" );
01391 return;
01392 }
01393
01394
01395
01396 void H2_Punch_line_data(
01397
01398 FILE* ioPUN ,
01399
01400 bool lgDoAll )
01401 {
01402 long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo;
01403
01404 if( !h2.lgH2ON )
01405 return;
01406
01407 DEBUG_ENTRY( "H2_Punch_line_data()" );
01408
01409 if( lgDoAll )
01410 {
01411 fprintf( ioQQQ,
01412 " H2_Punch_line_data ALL option not implemented in H2_Punch_line_data yet 1\n" );
01413 puts( "[Stop in H2_Punch_line_data]" );
01414 cdEXIT(EXIT_FAILURE);
01415 }
01416 else
01417 {
01418
01419
01420 for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi )
01421 {
01422 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
01423 {
01424 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
01425 {
01426
01427
01428
01429
01430 long int lim_elec_lo = 0;
01431 for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
01432 {
01433
01434
01435 long int nv = h2.nVib_hi[iElecLo];
01436 if( iElecLo==iElecHi )
01437 nv = iVibHi;
01438 for( iVibLo=0; iVibLo<=nv; ++iVibLo )
01439 {
01440 long nr = h2.nRot_hi[iElecLo][iVibLo];
01441 if( iElecLo==iElecHi && iVibHi==iVibLo )
01442 nr = iRotHi-1;
01443
01444 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
01445 {
01446
01447 if( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Aul > 0. )
01448 {
01450 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].cs = 0.;
01451
01452 fprintf(ioPUN,"%2li %2li %2li %2li %2li %2li ",
01453 iElecHi,iVibHi,iRotHi,iElecLo,iVibLo,iRotLo );
01454 Punch1LineData( &H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] , ioPUN , false);
01455 }
01456 }
01457 }
01458 }
01459 }
01460 }
01461 }
01462 fprintf( ioPUN , "\n");
01463 }
01464
01465 DEBUG_EXIT( "H2_Punch_line_data()" );
01466 return;
01467 }
01468
01469
01470 void H2_PunchLineStuff( FILE * io , float xLimit , long index)
01471 {
01472 long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo;
01473
01474 if( !h2.lgH2ON )
01475 return;
01476
01477 DEBUG_ENTRY( "H2_PunchLineStuff()" );
01478
01479
01480 for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi )
01481 {
01482 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
01483 {
01484 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
01485 {
01486
01487
01488
01489
01490 long int lim_elec_lo = 0;
01491 for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
01492 {
01493
01494
01495 long int nv = h2.nVib_hi[iElecLo];
01496 if( iElecLo==iElecHi )
01497 nv = iVibHi;
01498 for( iVibLo=0; iVibLo<=nv; ++iVibLo )
01499 {
01500 long nr = h2.nRot_hi[iElecLo][iVibLo];
01501 if( iElecLo==iElecHi && iVibHi==iVibLo )
01502 nr = iRotHi-1;
01503
01504 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
01505 {
01506
01507
01508 if( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Aul > 0. )
01509 {
01510 pun1Line( &H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] , io , xLimit , index , 1.);
01511 }
01512 }
01513 }
01514 }
01515 }
01516 }
01517 }
01518
01519 DEBUG_EXIT( "H2_PunchLineStuff()" );
01520
01521 return;
01522 }
01523
01524
01525
01526 void H2_Prt_line_tau(void)
01527 {
01528 long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo;
01529
01530 if( !h2.lgH2ON )
01531 return;
01532
01533 DEBUG_ENTRY( "H2_Prt_line_tau()" );
01534
01535
01536 for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi )
01537 {
01538 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
01539 {
01540 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
01541 {
01542
01543
01544
01545
01546 long int lim_elec_lo = 0;
01547 for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
01548 {
01549
01550
01551 long int nv = h2.nVib_hi[iElecLo];
01552 if( iElecLo==iElecHi )
01553 nv = iVibHi;
01554 for( iVibLo=0; iVibLo<=nv; ++iVibLo )
01555 {
01556 long nr = h2.nRot_hi[iElecLo][iVibLo];
01557 if( iElecLo==iElecHi && iVibHi==iVibLo )
01558 nr = iRotHi-1;
01559
01560 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
01561 {
01562
01563 if( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Aul > 0. )
01564 {
01565 prme(" c",&H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] );
01566 }
01567 }
01568 }
01569 }
01570 }
01571 }
01572 }
01573
01574 DEBUG_EXIT( "H2_Prt_line_tau()" );
01575
01576 return;
01577 }
01578
01579
01580
01581 static char chMolBranch( long iRotHi , long int iRotLo )
01582 {
01583
01584 char chBranch[5] = {'O','P','Q','R','S'};
01585
01586 int ip = 2 + (iRotHi - iRotLo);
01587 if( ip<0 || ip>=5 )
01588 {
01589 fprintf(ioQQQ," chMolBranch called with insane iRotHi=%li iRotLo=%li ip=%i\n",
01590 iRotHi , iRotLo , ip );
01591 ip = 0;
01592 }
01593
01594 return( chBranch[ip] );
01595 }
01596
01597
01598 void H2_PunchDo( FILE* io , char chJOB[] , const char chTime[] , long int ipPun )
01599 {
01600 long int iVibHi , iElecHi , iRotHi , iVibLo , iElecLo , iRotLo,
01601 ip;
01602 long int ilo , iRot , iVib;
01603 long int LimVib , LimRot;
01604
01605 DEBUG_ENTRY( "H2_PunchDo()" );
01606
01607
01608
01609
01610
01611
01612 if( (strcmp( chJOB , "H2po" ) == 0) && (strcmp(chTime,"LAST") == 0) &&
01613 (punch.punarg[ipPun][2] != 0) )
01614 {
01615
01616 if( h2.lgH2ON && hmi.lgBigH2_evaluated )
01617 {
01618 iVibHi= 0;
01619 iRotHi = 0;
01620 iElecHi=0;
01621
01622
01623
01624
01625 if( punch.punarg[ipPun][0] > 0 )
01626 {
01627 LimVib = (long)punch.punarg[ipPun][0];
01628 }
01629 else
01630 {
01631 LimVib = h2.nVib_hi[iElecHi];
01632 }
01633
01634
01635 fprintf(io,"%i\t%i\t%.3e\tortho\n",
01636 103 ,
01637 103 ,
01638 h2.ortho_density );
01639 fprintf(io,"%i\t%i\t%.3e\tpara\n",
01640 101 ,
01641 101 ,
01642 h2.para_density );
01643 fprintf(io,"%i\t%i\t%.3e\ttotal\n",
01644 0 ,
01645 0 ,
01646 hmi.H2_total );
01647
01648
01649 for( iVibHi=0; iVibHi<=LimVib; ++iVibHi )
01650 {
01651
01652 if( punch.punarg[ipPun][1] > 0 )
01653 {
01654 LimRot = (long)punch.punarg[ipPun][1];
01655 }
01656 else
01657 {
01658 LimRot = h2.nRot_hi[iElecHi][iVibHi];
01659 }
01660 if( punch.punarg[ipPun][2] > 0 )
01661 {
01662 long int i;
01663
01664 if( iVibHi == 0 )
01665 {
01666 fprintf(io,"vib\\rot");
01667
01668 for( i=0; i<=LimRot; ++i )
01669 {
01670 fprintf(io,"\t%li",i);
01671 }
01672 fprintf(io,"\n");
01673 }
01674 fprintf(io,"%li",iVibHi );
01675 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=LimRot; ++iRotHi )
01676 {
01677 fprintf(io,"\t%.3e",
01678 H2_populations[iElecHi][iVibHi][iRotHi]/hmi.H2_total );
01679 }
01680 fprintf(io,"\n" );
01681 }
01682 else if( punch.punarg[ipPun][2] < 0 )
01683 {
01684
01685 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=LimRot; ++iRotHi )
01686 {
01687 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",
01688
01689 iVibHi , iRotHi ,
01690
01691 chlgPara[H2_lgOrtho[iElecHi][iVibHi][iRotHi]],
01692
01693 energy_wn[iElecHi][iVibHi][iRotHi],
01694
01695 H2_populations[iElecHi][iVibHi][iRotHi]/hmi.H2_total ,
01696
01697 H2_old_populations[iElecHi][iVibHi][iRotHi]/hmi.H2_total ,
01698
01699 H2_populations[iElecHi][iVibHi][iRotHi]/hmi.H2_total/H2_stat[iElecHi][iVibHi][iRotHi] ,
01700
01701
01702 H2_populations[iElecHi][iVibHi][iRotHi]/SDIV(H2_populations_LTE[iElecHi][iVibHi][iRotHi]*hmi.H2_total ) ,
01703
01704 H2_col_rate_out[iVibHi][iRotHi]/SDIV(H2_col_rate_out[iVibHi][iRotHi]+H2_rad_rate_out[0][iVibHi][iRotHi]) ,
01705
01706 H2_col_rate_in[iVibHi][iRotHi]/SDIV(H2_col_rate_in[iVibHi][iRotHi]+H2_rad_rate_in[iVibHi][iRotHi]),
01707
01708 H2_col_rate_out[iVibHi][iRotHi],
01709
01710 H2_rad_rate_out[0][iVibHi][iRotHi] ,
01711
01712 H2_col_rate_in[iVibHi][iRotHi],
01713
01714 H2_rad_rate_in[iVibHi][iRotHi]
01715 );
01716 }
01717 }
01718 }
01719 }
01720 }
01721
01722
01723 else if( (strcmp( chJOB , "H2po" ) == 0) && (strcmp(chTime,"LAST") != 0) &&
01724 (punch.punarg[ipPun][2] == 0) )
01725 {
01726
01727 if( h2.lgH2ON && hmi.lgBigH2_evaluated )
01728 {
01729 fprintf(io,"%.5e\t%.3e\t%.3e", radius.depth_mid_zone ,
01730 h2.ortho_density , h2.para_density);
01731
01732 fprintf(io,"\t%.3e\t%.3e",
01733 pops_per_elec[1] , pops_per_elec[2]);
01734 iElecHi = 0;
01735 iVibHi = 0;
01736
01737 if( punch.punarg[ipPun][0] > 0 )
01738 {
01739 LimVib = (long)punch.punarg[ipPun][1];
01740 }
01741 else
01742 {
01743 LimVib = h2.nRot_hi[iElecHi][iVibHi];
01744 }
01745 LimVib = MIN2( LimVib , h2.nVib_hi[iElecHi] );
01746
01747 if( punch.punarg[ipPun][1] > 0 )
01748 {
01749 LimRot = (long)punch.punarg[ipPun][1];
01750 }
01751 else
01752 {
01753 LimRot = h2.nRot_hi[iElecHi][iVibHi];
01754 }
01755 for( iVibHi = 0; iVibHi<=LimVib; ++iVibHi )
01756 {
01757 long int LimRotVib = MIN2( LimRot , h2.nRot_hi[iElecHi][iVibHi] );
01758 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=LimRotVib; ++iRotHi )
01759 {
01760 fprintf(io,"\t%.3e",
01761 H2_populations[iElecHi][iVibHi][iRotHi]/hmi.H2_total );
01762 }
01763 fprintf(io,"\t");
01764 }
01765 fprintf(io,"\n");
01766 }
01767 }
01768
01769
01770 else if( (strcmp( chJOB , "H2cl" ) == 0) && (strcmp(chTime,"LAST") == 0) )
01771 {
01772 iVibHi= 0;
01773 iRotHi = 0;
01774 iElecHi=0;
01775
01776
01777
01778
01779 if( punch.punarg[ipPun][0] > 0 )
01780 {
01781 LimVib = (long)punch.punarg[ipPun][0];
01782 }
01783 else
01784 {
01785 LimVib = h2.nVib_hi[iElecHi];
01786 }
01787
01788
01789 fprintf(io,"%i\t%i\t%.3e\tortho\n",
01790 103 ,
01791 103 ,
01792 h2.ortho_colden );
01793 fprintf(io,"%i\t%i\t%.3e\tpara\n",
01794 101 ,
01795 101 ,
01796 h2.para_colden );
01797
01798 fprintf(io,"%i\t%i\t%.3e\ttotal\n",
01799 0 ,
01800 0 ,
01801 colden.colden[ipCOL_H2g]+colden.colden[ipCOL_H2s]);
01802
01803
01804 for( iVibHi=0; iVibHi<=LimVib; ++iVibHi )
01805 {
01806 if( h2.lgH2ON )
01807 {
01808
01809 if( punch.punarg[ipPun][1] > 0 )
01810 {
01811 LimRot = (long)punch.punarg[ipPun][1];
01812 }
01813 else
01814 {
01815 LimRot = h2.nRot_hi[iElecHi][iVibHi];
01816 }
01817 if( punch.punarg[ipPun][2] > 0 )
01818 {
01819 long int i;
01820
01821 if( iVibHi == 0 )
01822 {
01823 fprintf(io,"vib\\rot");
01824
01825 for( i=0; i<=LimRot; ++i )
01826 {
01827 fprintf(io,"\t%li",i);
01828 }
01829 fprintf(io,"\n");
01830 }
01831 fprintf(io,"%li",iVibHi );
01832 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=LimRot; ++iRotHi )
01833 {
01834 fprintf(io,"\t%.3e",
01835 H2_X_colden[iVibHi][iRotHi]/hmi.H2_total );
01836 }
01837 fprintf(io,"\n" );
01838 }
01839 else
01840 {
01841
01842 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=LimRot; ++iRotHi )
01843 {
01844 fprintf(io,"%li\t%li\t%.1f\t%.3e\t%.3e\t%.3e\t%.3e\n",
01845 iVibHi ,
01846 iRotHi ,
01847
01848 energy_wn[iElecHi][iVibHi][iRotHi]*T1CM,
01849
01850 H2_X_colden[iVibHi][iRotHi] ,
01851 H2_X_colden[iVibHi][iRotHi]/H2_stat[iElecHi][iVibHi][iRotHi] ,
01852
01853 H2_X_colden_LTE[iVibHi][iRotHi] ,
01854 H2_X_colden_LTE[iVibHi][iRotHi]/H2_stat[iElecHi][iVibHi][iRotHi]);
01855 }
01856 }
01857 }
01858 }
01859 }
01860 else if( (strcmp(chJOB , "H2pd" ) == 0) && (strcmp(chTime,"LAST") != 0) )
01861 {
01862
01863
01864 fprintf(io,"%.5e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
01865
01866 radius.depth_mid_zone ,
01867
01868 h2.ortho_density ,
01869 h2.para_density ,
01870
01871 hmi.H2_Solomon_dissoc_rate_TH85_H2g ,
01872
01873 hmi.H2_Solomon_dissoc_rate_BD96_H2g,
01874
01875 hmi.H2_Solomon_dissoc_rate_BigH2_H2g);
01876 }
01877 else if( (strcmp(chJOB , "H2co" ) == 0) && (strcmp(chTime,"LAST") != 0) )
01878 {
01879
01880 fprintf(io,"%.5e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
01881
01882 radius.depth_mid_zone ,
01883
01884 thermal.ctot ,
01885
01886 hmi.H2_Solomon_dissoc_rate_TH85_H2g,
01887
01888 hmi.H2_Solomon_dissoc_rate_BigH2_H2g +
01889 hmi.H2_Solomon_dissoc_rate_BigH2_H2s,
01890
01891 hmi.HeatH2Dish_TH85,
01892
01893 hmi.HeatH2Dish_BigH2 ,
01894
01895 hmi.HeatH2Dexc_TH85,
01896 hmi.HeatH2Dexc_BigH2
01897 );
01898
01899 }
01900 else if( (strcmp(chJOB , "H2cr" ) == 0) && (strcmp(chTime,"LAST") != 0) )
01901 {
01902
01903
01904 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",
01905
01906 radius.depth_mid_zone ,
01907
01908 hmi.H2_rate_create,
01909 hmi.H2_rate_destroy,
01910 gv.rate_h2_form_grains_used_total * dense.xIonDense[ipHYDROGEN][0] / hmi.H2_rate_create,
01911 hmi.assoc_detach * hmi.Hmolec[ipMH]*hmi.Hmolec[ipMHm] / hmi.H2_rate_create,
01912 hmi.bh2dis * dense.xIonDense[ipHYDROGEN][0] / hmi.H2_rate_create,
01913 hmi.bh2h2p * dense.xIonDense[ipHYDROGEN][0] * hmi.Hmolec[ipMH2p] / hmi.H2_rate_create,
01914 hmi.radasc * dense.xIonDense[ipHYDROGEN][0] / hmi.H2_rate_create,
01915 hmi.h3ph2p * dense.xIonDense[ipHYDROGEN][0] * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create,
01916 hmi.h2phmh2h * hmi.Hmolec[ipMH2p] * hmi.Hmolec[ipMHm] / hmi.H2_rate_create,
01917 hmi.bh2h22hh2 * 2 * dense.xIonDense[ipHYDROGEN][0] * hmi.Hmolec[ipMH2g] / hmi.H2_rate_create,
01918 hmi.h3phmh2hh * hmi.Hmolec[ipMH3p] * hmi.Hmolec[ipMHm] / hmi.H2_rate_create,
01919 hmi.h3phm2h2 * hmi.Hmolec[ipMH3p] * hmi.Hmolec[ipMHm] / hmi.H2_rate_create,
01920 hmi.h32h2 * hmi.Hmolec[ipMH2p] * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create,
01921 hmi.eh3_h2h * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create,
01922 hmi.h3ph2hp * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create
01923 );
01924 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",
01925
01926
01927 co.H_CH_C_H2 * dense.xIonDense[ipHYDROGEN][0] / hmi.H2_rate_create,
01928 co.H_CHP_CP_H2 * dense.xIonDense[ipHYDROGEN][0] / hmi.H2_rate_create,
01929 co.H_CH2_CH_H2 * dense.xIonDense[ipHYDROGEN][0] / hmi.H2_rate_create,
01930 co.H_CH3P_CH2P_H2 * dense.xIonDense[ipHYDROGEN][0] / hmi.H2_rate_create,
01931 co.H_OH_O_H2 * dense.xIonDense[ipHYDROGEN][0] / hmi.H2_rate_create,
01932 co.Hminus_HCOP_CO_H2 * hmi.Hmolec[ipMHm] / hmi.H2_rate_create,
01933 co.Hminus_H3OP_H2O_H2 * hmi.Hmolec[ipMHm] / hmi.H2_rate_create,
01934 co.Hminus_H3OP_OH_H2_H * hmi.Hmolec[ipMHm] / hmi.H2_rate_create,
01935 co.HP_CH2_CHP_H2 * hmi.Hmolec[ipMHp] / hmi.H2_rate_create,
01936 co.HP_SiH_SiP_H2* hmi.Hmolec[ipMHp] / hmi.H2_rate_create,
01937 co.H2P_CH_CHP_H2 * hmi.Hmolec[ipMH2p] / hmi.H2_rate_create,
01938 co.H2P_CH2_CH2P_H2 * hmi.Hmolec[ipMH2p] / hmi.H2_rate_create,
01939 co.H2P_CO_COP_H2 * hmi.Hmolec[ipMH2p] / hmi.H2_rate_create,
01940 co.H2P_H2O_H2OP_H2 * hmi.Hmolec[ipMH2p] / hmi.H2_rate_create,
01941 co.H2P_O2_O2P_H2 * hmi.Hmolec[ipMH2p] / hmi.H2_rate_create,
01942 co.H2P_OH_OHP_H2 * hmi.Hmolec[ipMH2p] / hmi.H2_rate_create,
01943 co.H3P_C_CHP_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create,
01944 co.H3P_CH_CH2P_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create,
01945 co.H3P_CH2_CH3P_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create,
01946 co.H3P_OH_H2OP_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create,
01947 co.H3P_H2O_H3OP_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create,
01948 co.H3P_CO_HCOP_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create,
01949 co.H3P_O_OHP_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create,
01950 co.H3P_SiH_SiH2P_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create,
01951 co.H3P_SiO_SiOHP_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create,
01952 co.H_CH3_CH2_H2 * dense.xIonDense[ipHYDROGEN][0] / hmi.H2_rate_create,
01953 co.H_CH4P_CH3P_H2 * dense.xIonDense[ipHYDROGEN][0] / hmi.H2_rate_create,
01954 co.H_CH5P_CH4P_H2 * dense.xIonDense[ipHYDROGEN][0] / hmi.H2_rate_create,
01955 co.H2P_CH4_CH3P_H2 * hmi.Hmolec[ipMH2p] / hmi.H2_rate_create,
01956 co.H2P_CH4_CH4P_H2 * hmi.Hmolec[ipMH2p] / hmi.H2_rate_create,
01957 co.H3P_CH3_CH4P_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create,
01958 co.H3P_CH4_CH5P_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create,
01959 co.HP_CH4_CH3P_H2 * hmi.Hmolec[ipMHp] / hmi.H2_rate_create,
01960
01961 co.HP_HNO_NOP_H2 * hmi.Hmolec[ipMHp] / hmi.H2_rate_create,
01962 co.HP_HS_SP_H2 * hmi.Hmolec[ipMHp] / hmi.H2_rate_create,
01963 co.H_HSP_SP_H2 * dense.xIonDense[ipHYDROGEN][0] / hmi.H2_rate_create,
01964 co.H3P_NH_NH2P_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create,
01965 co.H3P_NH2_NH3P_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create,
01966 co.H3P_NH3_NH4P_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create,
01967 co.H3P_CN_HCNP_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create,
01968 co.H3P_NO_HNOP_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create,
01969 co.H3P_S_HSP_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create,
01970 co.H3P_CS_HCSP_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create,
01971 co.H3P_NO2_NOP_OH_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create,
01972 co.H2P_NH_NHP_H2 * hmi.Hmolec[ipMH2p] / hmi.H2_rate_create,
01973 co.H2P_NH2_NH2P_H2 * hmi.Hmolec[ipMH2p] / hmi.H2_rate_create,
01974 co.H2P_NH3_NH3P_H2 * hmi.Hmolec[ipMH2p] / hmi.H2_rate_create,
01975 co.H2P_CN_CNP_H2 * hmi.Hmolec[ipMH2p] / hmi.H2_rate_create,
01976 co.H2P_HCN_HCNP_H2 * hmi.Hmolec[ipMH2p] / hmi.H2_rate_create,
01977 co.H2P_NO_NOP_H2 * hmi.Hmolec[ipMH2p] / hmi.H2_rate_create,
01978 co.H3P_Cl_HClP_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create,
01979 co.H3P_HCl_H2ClP_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create,
01980 co.H2P_C2_C2P_H2 * hmi.Hmolec[ipMH2p] / hmi.H2_rate_create,
01981 co.Hminus_NH4P_NH3_H2 * hmi.Hmolec[ipMHm] / hmi.H2_rate_create,
01982 co.H3P_HCN_HCNHP_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create
01983 );
01984 fprintf(io,"\t%.3e\t%.3e\n",
01985 hmi.H2_rate_destroy * hmi.H2_total / hmi.H2_rate_create,
01986 hmi.h2s_sp_decay
01987 );
01988 }
01989 else if( (strcmp(chJOB , "H2ds" ) == 0) && (strcmp(chTime,"LAST") != 0) )
01990 {
01991
01992
01993
01994
01995 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\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",
01996
01997 radius.depth_mid_zone ,
01998
01999 hmi.H2_rate_create,
02000
02001 hmi.H2_rate_destroy ,
02002
02003 (hmi.assoc_detach_backwards_grnd * hmi.Hmolec[ipMH2g] + hmi.assoc_detach_backwards_exct * hmi.Hmolec[ipMH2s]) / hmi.H2_rate_destroy / hmi.H2_total,
02004
02005 hmi.H2_Solomon_dissoc_rate_used_H2g / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
02006 hmi.H2_Solomon_dissoc_rate_used_H2s / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2s] / hmi.H2_total,
02007 hmi.H2_photodissoc_used_H2s / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2s] / hmi.H2_total,
02008 hmi.H2_photodissoc_used_H2g / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
02009
02010 (hmi.h2ge2h * hmi.Hmolec[ipMH2g] + hmi.h2se2h * hmi.Hmolec[ipMH2s]) / hmi.H2_rate_destroy / hmi.H2_total,
02011
02012 hmi.rh2h2p*dense.xIonDense[ipHYDROGEN][1] / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
02013 hmi.h2hph3p*dense.xIonDense[ipHYDROGEN][1] / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
02014
02015 (hmi.rh2dis * hmi.Hmolec[ipMH2g] + hmi.h2sh * hmi.Hmolec[ipMH2s]) * dense.xIonDense[ipHYDROGEN][0] / hmi.H2_rate_destroy / hmi.H2_total,
02016
02017 (hmi.CR_reac_H2g * hmi.Hmolec[ipMH2g] + hmi.CR_reac_H2s * hmi.Hmolec[ipMH2s]) / hmi.H2_rate_destroy / hmi.H2_total,
02018
02019 hmi.rheph2hpheh*dense.xIonDense[ipHELIUM][1] / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
02020 hmi.heph2heh2p*dense.xIonDense[ipHELIUM][1] / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
02021 hmi.hehph2h3phe*hmi.Hmolec[ipMHeHp] / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
02022
02023 hmi.h3petc*hmi.Hmolec[ipMH3p] / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
02024
02025 hmi.h2ph3p*hmi.Hmolec[ipMH2p] / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
02026
02027 hmi.h2sh2g * hmi.Hmolec[ipMH2g] / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2s] / hmi.H2_total,
02028
02029 hmi.h2h22hh2*2*hmi.Hmolec[ipMH2g] / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
02030
02031 hmi.h2sh2sh2g2h*2*hmi.Hmolec[ipMH2s] / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2s] / hmi.H2_total,
02032
02033 hmi.h2sh2sh2s2h*2*hmi.Hmolec[ipMH2s] / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2s] / hmi.H2_total,
02034
02035
02036 co.H2_CHP_CH2P_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
02037 co.H2_CH2P_CH3P_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
02038 co.H2_OHP_H2OP_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
02039 co.H2_H2OP_H3OP_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
02040 co.H2_COP_HCOP_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
02041 co.H2_OP_OHP_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
02042 co.H2_SiOP_SiOHP_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
02043 co.H2_C_CH_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
02044 co.H2_CP_CHP_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
02045 co.H2_CH_CH2_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
02046 co.H2_OH_H2O_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
02047 co.H2_O_OH_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
02048 co.H2s_CH_CH2_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2s] / hmi.H2_total,
02049 co.H2s_O_OH_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2s] / hmi.H2_total,
02050 co.H2s_OH_H2O_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2s] / hmi.H2_total,
02051 co.H2s_C_CH_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2s] / hmi.H2_total,
02052 co.H2s_CP_CHP_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2s] / hmi.H2_total,
02053 co.H2_CH2_CH3_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
02054 co.H2_CH3_CH4_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
02055 co.H2_CH4P_CH5P_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
02056 co.H2s_CH2_CH3_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2s] / hmi.H2_total,
02057 co.H2s_CH3_CH4_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2s] / hmi.H2_total,
02058 co.H2s_OP_OHP_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2s] / hmi.H2_total,
02059
02060 co.H2_N_NH_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
02061 co.H2_NH_NH2_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
02062 co.H2_NH2_NH3_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
02063 co.H2_CN_HCN_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
02064 co.H2_NP_NHP_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
02065 co.H2_NHP_N_H3P / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
02066 co.H2_NHP_NH2P_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
02067 co.H2_NH2P_NH3P_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
02068 co.H2_NH3P_NH4P_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
02069 co.H2_CNP_HCNP_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
02070 co.H2_SP_HSP_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
02071 co.H2_CSP_HCSP_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
02072 co.H2_ClP_HClP_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
02073 co.H2_HClP_H2ClP_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
02074 co.H2_HCNP_HCNHP_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total
02075 );
02076 fprintf(io,"\t%.4e\t%.4e\n",
02077
02078 hmi.Hmolec[ipMH2g] / hmi.H2_total,
02079 hmi.Hmolec[ipMH2s] / hmi.H2_total
02080 );
02081
02082 }
02083 else if( (strcmp(chJOB , "H2le" ) == 0) && (strcmp(chTime,"LAST") != 0) )
02084 {
02085 for( ilo=0; ilo < nLevels_per_elec[0]; ilo++ )
02086 {
02087 ip = H2_ipX_ener_sort[ilo];
02088 iRot = ipRot_H2_energy_sort[ip];
02089 iVib = ipVib_H2_energy_sort[ip];
02090
02091
02092 fprintf(io,"%.2f\t%li\t%li\t%li\n",
02093 energy_wn[0][iVib][iRot],
02094 (long)H2_stat[0][iVib][iRot],
02095 iVib , iRot );
02096 }
02097
02098 }
02099 else if( (strcmp(chJOB , "H2ra" ) == 0) && (strcmp(chTime,"LAST") != 0) )
02100 {
02101
02102 double sumpop = 0. , sumlife = 0.;
02103
02104
02105 iElecLo = 0;
02106 iVibLo = 0;
02107 if( h2.lgH2ON && hmi.lgBigH2_evaluated )
02108 {
02109 for( iElecHi=1; iElecHi<mole.n_h2_elec_states; ++iElecHi )
02110 {
02111 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
02112 {
02113 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
02114 {
02115 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=h2.nRot_hi[iElecLo][iVibLo]; ++iRotLo )
02116 {
02117
02118 if( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Aul > 0. )
02119 {
02120 sumlife +=
02121 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].pump *
02122 H2_populations[iElecLo][iVibLo][iRotLo];
02123 sumpop +=
02124 H2_populations[iElecLo][iVibLo][iRotLo];
02125 }
02126 }
02127 }
02128 }
02129 }
02130 }
02131
02132
02133
02134
02135 fprintf(io,
02136 "%.5e\t%.3e\t%.3e\t%.3e\t%li",
02137
02138 radius.depth_mid_zone ,
02139
02140 colden.colden[ipCOL_H2g]+colden.colden[ipCOL_H2s],
02141
02142
02143 colden.coldenH2_ov_vel ,
02144
02145 rfield.extin_mag_V_point,
02146
02147 h2.nCallH2_this_zone );
02148 fprintf(io,
02149 "\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e",
02150
02151 hmi.H2_total/dense.gas_phase[ipHYDROGEN] ,
02152
02153 H2_renorm_chemistry,
02154
02155 gv.rate_h2_form_grains_used_total ,
02156
02157 hmi.Hmolec[ipMHm]*1.35e-9,
02158
02159 hmi.H2_Solomon_dissoc_rate_TH85_H2g + hmi.H2_Solomon_dissoc_rate_TH85_H2s,
02160
02161 hmi.H2_Solomon_dissoc_rate_BD96_H2g + hmi.H2_Solomon_dissoc_rate_BD96_H2s,
02162
02163 hmi.H2_Solomon_dissoc_rate_ELWERT_H2g + hmi.H2_Solomon_dissoc_rate_ELWERT_H2g,
02164
02165 hmi.H2_Solomon_dissoc_rate_BigH2_H2g + hmi.H2_Solomon_dissoc_rate_BigH2_H2s,
02166
02167 hmi.H2_Solomon_elec_decay_H2g ,
02168
02169 hmi.H2_Solomon_elec_decay_H2s
02170 );
02171 fprintf(io,
02172 "\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e",
02173
02174 hmi.UV_Cont_rel2_Habing_TH85_depth,
02175
02176 hmi.UV_Cont_rel2_Draine_DB96_depth,
02177
02178 secondaries.csupra[ipHYDROGEN][0]*0.93,
02179 sumlife/SDIV( sumpop ) ,
02180 hmi.H2_Solomon_dissoc_rate_BD96_H2g/SDIV(hmi.UV_Cont_rel2_Habing_TH85_depth) ,
02181 hmi.H2_Solomon_dissoc_rate_BigH2_H2g/SDIV(hmi.UV_Cont_rel2_Habing_TH85_depth),
02182 hmi.H2_Solomon_dissoc_rate_BigH2_H2g/SDIV(hmi.UV_Cont_rel2_Habing_spec_depth),
02183 hmi.H2_rate_destroy);
02184 fprintf(io,
02185 "\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
02186 hmi.HeatH2Dish_TH85,
02187 hmi.HeatH2Dexc_TH85,
02188 hmi.HeatH2Dish_BigH2,
02189 hmi.HeatH2Dexc_BigH2,
02190 thermal.htot);
02191 }
02192
02193 else if( (strcmp(chJOB , "H2so" ) == 0) && (strcmp(chTime,"LAST") != 0) )
02194 {
02195
02196 # define NSOL 100
02197 double sum, one;
02198 long int jlosave[NSOL] , ivlosave[NSOL],
02199 iehisave[NSOL] ,jhisave[NSOL] , ivhisave[NSOL],
02200 nsave,
02201 ipOrdered[NSOL];
02202 int nFail;
02203 int i;
02204 float fsave[NSOL], wlsave[NSOL];
02205
02206 fprintf(io,"%.5e\t%.3e",
02207
02208 radius.depth_mid_zone ,
02209
02210 hmi.H2_Solomon_dissoc_rate_BigH2_H2g +
02211 hmi.H2_Solomon_dissoc_rate_BigH2_H2s);
02212 sum = 0.;
02213 iElecLo = 0;
02214
02215 if( h2.lgH2ON && hmi.lgBigH2_evaluated )
02216 {
02217 for( iElecHi=1; iElecHi<mole.n_h2_elec_states; ++iElecHi )
02218 {
02219 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
02220 {
02221 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
02222 {
02223 long int nv = h2.nVib_hi[iElecLo];
02224 for( iVibLo=0; iVibLo<=nv; ++iVibLo )
02225 {
02226 long nr = h2.nRot_hi[iElecLo][iVibLo];
02227 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
02228 {
02229 if( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Aul > 0. )
02230 {
02231 one = H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].PopLo *
02232 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].pump;
02233 sum += one;
02234 }
02235 }
02236 }
02237 }
02238 }
02239 }
02240
02241
02242 sum = SDIV( sum );
02243 nsave = 0;
02244
02245 # define FRAC 0.01
02246 for( iElecHi=1; iElecHi<mole.n_h2_elec_states; ++iElecHi )
02247 {
02248 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
02249 {
02250 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
02251 {
02252 long int nv = h2.nVib_hi[iElecLo];
02253 for( iVibLo=0; iVibLo<=nv; ++iVibLo )
02254 {
02255 long nr = h2.nRot_hi[iElecLo][iVibLo];
02256 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
02257 {
02258 if( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Aul > 0. )
02259 {
02260 one = H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].PopLo *
02261 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].pump;
02262 if( one/sum > FRAC && nsave<NSOL)
02263 {
02264 fsave[nsave] = (float)(one/sum);
02265 jlosave[nsave] = iRotLo;
02266 ivlosave[nsave] = iVibLo;
02267 jhisave[nsave] = iRotHi;
02268 ivhisave[nsave] = iVibHi;
02269 iehisave[nsave] = iElecHi;
02270 wlsave[nsave] = H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].WLAng;
02271 ++nsave;
02272
02273
02274 }
02275 }
02276 }
02277 }
02278 }
02279 }
02280 }
02281
02282
02283
02284
02285 spsort(
02286
02287 fsave,
02288
02289 nsave,
02290
02291 ipOrdered,
02292
02293
02294 -1,
02295
02296 &nFail);
02297
02298
02299
02300
02301 fprintf(io,"\t%.3f\t%.3f\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e",
02302
02303 (sum + secondaries.csupra[ipHYDROGEN][0]*2.02f)/SDIV((hmi.H2_Solomon_dissoc_rate_BigH2_H2g * hmi.Hmolec[ipMH2g] +
02304 hmi.H2_Solomon_dissoc_rate_BigH2_H2s * hmi.Hmolec[ipMH2s]) ),
02305
02306 (sum + secondaries.csupra[ipHYDROGEN][0]*2.02f) /SDIV((hmi.H2_Solomon_dissoc_rate_BigH2_H2g * hmi.H2g_BigH2 +
02307 hmi.H2_Solomon_dissoc_rate_BigH2_H2s * hmi.H2s_BigH2) ),
02308 hmi.H2_BigH2_H2g_av, hmi.H2_BigH2_H2s_av,
02309 hmi.H2_chem_BigH2_H2g, hmi.H2_chem_BigH2_H2s,
02310 hmi.H2g_BigH2/SDIV(hmi.H2_total_BigH2), hmi.H2s_BigH2/SDIV(hmi.H2_total_BigH2)
02311 );
02312 for( i=0; i<nsave; ++i )
02313 {
02314 ip = ipOrdered[i];
02315
02316 fprintf(io,"\t%li\t%li\t%li\t%li\t%li\t%.3f\t%.3f",
02317 iehisave[ip],ivhisave[ip],jhisave[ip],ivlosave[ip] , jlosave[ip] , fsave[ip] , wlsave[ip] );
02318
02319 }
02320 fprintf(io,"\n");
02321 }
02322
02323
02324 # undef NSOL
02325 }
02326
02327 else if( (strcmp(chJOB , "H2te" ) == 0) && (strcmp(chTime,"LAST") != 0) )
02328 {
02329
02330 double pop_ratio10,pop_ratio20,pop_ratio30,pop_ratio31,pop_ratio40;
02331 double T10,T20,T30,T31,T40;
02332
02333 double T10_sum,T20_sum,T30_sum,T31_sum,T40_sum;
02334 double pop_ratio10_sum,pop_ratio20_sum,pop_ratio30_sum,pop_ratio31_sum,pop_ratio40_sum;
02335 if( h2.lgH2ON && h2.nCallH2_this_zone )
02336 {
02337 double energyK = T1CM*(energy_wn[0][0][1] - energy_wn[0][0][0]);
02338
02339 pop_ratio10 = H2_populations[0][0][1]/SDIV(H2_populations[0][0][0]);
02340 pop_ratio10_sum = H2_X_colden[0][1]/SDIV(H2_X_colden[0][0]);
02341
02342 T10 = -170.5/log(SDIV(pop_ratio10) * H2_stat[0][0][0]/H2_stat[0][0][1]);
02343 T10_sum = -170.5/log(SDIV(pop_ratio10_sum) * H2_stat[0][0][0]/H2_stat[0][0][1]);
02344
02345 energyK = T1CM*(energy_wn[0][0][2] - energy_wn[0][0][0]);
02346 pop_ratio20 = H2_populations[0][0][2]/SDIV(H2_populations[0][0][0]);
02347 T20 = -energyK/log(SDIV(pop_ratio20) * H2_stat[0][0][0]/H2_stat[0][0][2]);
02348
02349 pop_ratio20_sum = H2_X_colden[0][2]/SDIV(H2_X_colden[0][0]);
02350 T20_sum = -energyK/log(SDIV(pop_ratio20_sum) * H2_stat[0][0][0]/H2_stat[0][0][2]);
02351
02352 energyK = T1CM*(energy_wn[0][0][3] - energy_wn[0][0][0]);
02353 pop_ratio30 = H2_populations[0][0][3]/SDIV(H2_populations[0][0][0]);
02354 T30 = -energyK/log(SDIV(pop_ratio30) * H2_stat[0][0][0]/H2_stat[0][0][3]);
02355
02356 pop_ratio30_sum = H2_X_colden[0][3]/SDIV(H2_X_colden[0][0]);
02357 T30_sum = -energyK/log(SDIV(pop_ratio30_sum) * H2_stat[0][0][0]/H2_stat[0][0][3]);
02358
02359 energyK = T1CM*(energy_wn[0][0][3] - energy_wn[0][0][1]);
02360 pop_ratio31 = H2_populations[0][0][3]/SDIV(H2_populations[0][0][1]);
02361 T31 = -energyK/log(SDIV(pop_ratio31) * H2_stat[0][0][1]/H2_stat[0][0][3]);
02362
02363 pop_ratio31_sum = H2_X_colden[0][3]/SDIV(H2_X_colden[0][1]);
02364 T31_sum = -energyK/log(SDIV(pop_ratio31_sum) * H2_stat[0][0][1]/H2_stat[0][0][3]);
02365
02366 energyK = T1CM*(energy_wn[0][0][4] - energy_wn[0][0][0]);
02367 pop_ratio40 = H2_populations[0][0][4]/SDIV(H2_populations[0][0][0]);
02368 T40 = -energyK/log(SDIV(pop_ratio40) * H2_stat[0][0][0]/H2_stat[0][0][4]);
02369
02370 pop_ratio40_sum = H2_X_colden[0][4]/SDIV(H2_X_colden[0][0]);
02371 T40_sum = -energyK/log(SDIV(pop_ratio40_sum) * H2_stat[0][0][0]/H2_stat[0][0][4]);
02372 }
02373 else
02374 {
02375 pop_ratio10 = 0.;
02376 pop_ratio10_sum = 0.;
02377 T10 = 0.;
02378 T20 = 0.;
02379 T30 = 0.;
02380 T31 = 0.;
02381 T40 = 0.;
02382 T10_sum = 0.;
02383 T20_sum = 0.;
02384 T30_sum = 0.;
02385 T31_sum = 0.;
02386 T40_sum = 0.;
02387 }
02388
02389
02390 fprintf( io,
02391 "%.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" ,
02392
02393 radius.depth_mid_zone ,
02394
02395 hmi.H2_total/dense.gas_phase[ipHYDROGEN] ,
02396
02397 pop_ratio10 ,
02398
02399 h2.ortho_density / SDIV(h2.para_density),
02400 T10,T20,T30,T31,T40,
02401 phycon.te ,
02402 hyperfine.Tspin21cm,T10_sum,T20_sum,T30_sum,T31_sum,T40_sum );
02403 }
02404 else if( (strcmp(chJOB , "H2ln" ) == 0) && (strcmp(chTime,"LAST") == 0) )
02405 {
02406
02407 double thresh;
02408 double renorm;
02409 if( h2.lgH2ON )
02410 {
02411
02412 if( LineSv[LineSave.ipNormWavL].sumlin[LineSave.lgLineEmergent] > SMALLFLOAT )
02413 {
02414 renorm = LineSave.ScaleNormLine/LineSv[LineSave.ipNormWavL].sumlin[LineSave.lgLineEmergent];
02415 }
02416 else
02417 {
02418 renorm = 1.;
02419 }
02420 if( renorm > SMALLFLOAT )
02421 {
02422
02423
02424 thresh = thresh_punline_h2/(float)renorm;
02425 }
02426 else
02427 {
02428 thresh = 0.f;
02429 }
02430
02431
02432
02433 for( iElecHi=0; iElecHi < h2.nElecLevelOutput; ++iElecHi )
02434 {
02435 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
02436 {
02437 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
02438 {
02439
02440
02441
02442
02443 long int lim_elec_lo = 0;
02444 for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
02445 {
02446 long int nv = h2.nVib_hi[iElecLo];
02447 if( iElecLo==iElecHi )
02448 nv = iVibHi;
02449 for( iVibLo=0; iVibLo<=nv; ++iVibLo )
02450 {
02451 long nr = h2.nRot_hi[iElecLo][iVibLo];
02452 if( iElecLo==iElecHi && iVibHi==iVibLo )
02453 nr = iRotHi-1;
02454 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<nr; ++iRotLo )
02455 {
02456 if( H2_SaveLine[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] > thresh )
02457 {
02458
02459
02460 double wl = H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].WLAng/1e4;
02461
02462
02463 fprintf(io, "%li-%li %c(%li)",
02464 iVibHi ,
02465 iVibLo ,
02466 chMolBranch( iRotHi , iRotLo ) ,
02467 iRotLo );
02468 fprintf( io, "\t%ld\t%ld\t%ld\t%ld\t%ld\t%ld",
02469 iElecHi , iVibHi , iRotHi , iElecLo , iVibLo , iRotLo);
02470
02471 fprintf( io, "\t%.7f\t", wl );
02472
02473 prt_wl( io , H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].WLAng );
02474 fprintf( io, "\t%.3f\t%.3e",
02475 log10(MAX2(1e-37,H2_SaveLine[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo])) + radius.Conv2PrtInten,
02476 H2_SaveLine[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo]*renorm );
02477
02478 fprintf( io, "\t%.3f", energy_wn[iElecHi][iVibHi][iRotHi]*T1CM );
02479
02480 fprintf( io, "\t%.3e", H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Aul*
02481 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyErg *
02482 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].gHi);
02483 fprintf( io, "\n");
02484 }
02485 }
02486 }
02487 }
02488 }
02489 }
02490 }
02491 }
02492 }
02493 else if( (strcmp(chJOB , "H2sp" ) == 0) )
02494 {
02495 iVib = 0;
02496 iRot = 0;
02497 # if 0
02498
02499 fprintf(io,"%.4e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
02500 radius.depth_mid_zone ,
02501 H2_populations[0][iVib][iRot] ,
02502 radius.depth_mid_zone * H2_populations[0][iVib][iRot] ,
02503 H2_rad_rate_out[0][iVib][iRot] ,
02504 H2_rad_rate_in[iVib][iRot] ,
02505 H2_col_rate_out_ old[iVib][iRot] ,
02506 H2_col_rate_in_ old[iVib][iRot] );
02507 # endif
02508 fprintf(io,"%.4e\t%.2e\t%.2e\t%.2e\t%.2e\n",
02509 radius.depth_mid_zone ,
02510 H2_populations[0][iVib][iRot] ,
02511 H2Lines[1][1][1][0][iVib][iRot].pump,
02512 H2Lines[1][1][1][0][iVib][iRot].TauIn,
02513 H2Lines[1][1][1][0][iVib][iRot].TauCon);;
02514 }
02515
02516 DEBUG_EXIT( "H2_PunchDo()" );
02517 return;
02518 }
02519
02520
02521
02522
02523
02524 long int cdH2_Line(
02525
02526 long int iElecHi,
02527 long int iVibHi ,
02528 long int iRotHi ,
02529
02530 long int iElecLo,
02531 long int iVibLo ,
02532 long int iRotLo ,
02533
02534 double *relint,
02535
02536 double *absint )
02537 {
02538
02539 DEBUG_ENTRY( "cdH2_Line()" );
02540
02541
02542 *relint = 0.;
02543 *absint = 0.;
02544
02545
02546 if( iElecHi!=0 || iElecLo!=0 )
02547 {
02548 DEBUG_EXIT( "cdH2_Line()" );
02549 return 0;
02550 }
02551
02552
02553 if( energy_wn[iElecHi][iVibHi][iRotHi] < energy_wn[iElecLo][iVibHi][iRotHi] )
02554 {
02555 DEBUG_EXIT( "cdH2_Line()" );
02556 return 0;
02557 }
02558
02559
02560 if( H2_lgOrtho[iElecHi][iVibHi][iRotHi] - H2_lgOrtho[iElecLo][iVibLo][iRotLo] != 0 )
02561 {
02562 DEBUG_EXIT( "cdH2_Line()" );
02563 return 0;
02564 }
02565
02566
02567 if( !lgH2_line_exists[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] )
02568 {
02569 DEBUG_EXIT( "cdH2_Line()" );
02570 return 0;
02571 }
02572
02573
02574 if( LineSv[LineSave.ipNormWavL].sumlin[LineSave.lgLineEmergent] > 0. )
02575 {
02576 *relint = H2_SaveLine[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo]/
02577 LineSv[LineSave.ipNormWavL].sumlin[LineSave.lgLineEmergent] * LineSave.ScaleNormLine;
02578 }
02579 else
02580 {
02581 *relint = 0.;
02582 }
02583
02584
02585 if( H2_SaveLine[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] > 0. )
02586 {
02587 *absint = log10(H2_SaveLine[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo]) +
02588 radius.Conv2PrtInten;
02589 }
02590 else
02591 {
02592
02593 *absint = -37.;
02594 }
02595
02596 DEBUG_EXIT( "cdH2_Line()" );
02597
02598
02599 return 1;
02600 }
02601
02602
02603