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