00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023 #define PRT_POPS false
00024
00025 #define LIM_H2_POP_LOOP 100
00026
00027
00028 #define H2_DISS_ALLISON_DALGARNO 6e-19f
00029 #include "cddefines.h"
00030 #include "cddrive.h"
00031 #include "physconst.h"
00032 #include "taulines.h"
00033 #include "atoms.h"
00034 #include "conv.h"
00035 #include "secondaries.h"
00036 #include "geometry.h"
00037 #include "iterations.h"
00038 #include "trace.h"
00039 #include "hmi.h"
00040 #include "hextra.h"
00041 #include "rt.h"
00042 #include "radius.h"
00043 #include "ipoint.h"
00044 #include "phycon.h"
00045 #include "thermal.h"
00046 #include "dense.h"
00047 #include "rfield.h"
00048 #include "lines_service.h"
00049 #include "mole.h"
00050 #include "h2.h"
00051 #include "h2_priv.h"
00052
00053
00054
00055
00056
00057
00058
00059 static long int loop_h2_pops;
00060
00061 float H2_te_hminus[nTE_HMINUS] = {10.,30.,100.,300.,1000.,3000.,10000.};
00062
00063
00064
00065 static float collider_density[N_X_COLLIDER];
00066 static float collider_density_total;
00067
00068
00069
00070
00071 int H2_nRot_add_ortho_para[N_H2_ELEC] = {0 , 1 , 1 , 0, 1, 1 , 0};
00072
00073
00074
00075
00076
00077
00078
00079 double H2_DissocEnergies[N_H2_ELEC] =
00080 { 36118.11, 118375.6, 118375.6, 118375.6, 118375.6,133608.6,133608.6 };
00081
00082
00083
00084 double exp_disoc;
00085
00086
00087
00088
00089 static void H2_X_coll_rate_evaluate( void )
00090 {
00091 long int nColl,
00092 ipHi ,
00093 ipLo,
00094 ip,
00095 iVibHi,
00096 iRotHi,
00097 iVibLo,
00098 iRotLo;
00099 float colldown;
00100 double exph2_big,
00101 exph2_big_0,
00102 rel_pop_LTE_H2_big;
00103
00104 DEBUG_ENTRY( "H2_X_coll_rate_evaluate()" );
00105
00106
00107
00108
00109 collider_density[0] = dense.xIonDense[ipHYDROGEN][0];
00110
00111 collider_density[2] = (float)h2.ortho_density;
00112
00113 collider_density[3] = (float)h2.para_density;
00114
00115 collider_density[4] = dense.xIonDense[ipHYDROGEN][1];
00116
00117
00118 collider_density[4] += hmi.Hmolec[ipMH3p];
00119
00120
00121
00122
00123
00124 ASSERT( !(mole.lgH2_He_Meudon * mole.lgH2_He_Stancil) );
00125 collider_density[1] = dense.xIonDense[ipHELIUM][0] * mole.lgH2_He_Meudon;
00126
00127
00128 collider_density[5] = dense.xIonDense[ipHELIUM][0] * mole.lgH2_He_Stancil;
00129
00130
00131 collider_density_total = 0.;
00132 #if 0
00133 for( nColl=0; nColl<N_X_COLLIDER; ++nColl )
00134 {
00135 collider_density_total += collider_density[nColl];
00136 }
00137
00138
00139 collider_density_total += (float)dense.eden;
00140 #endif
00141
00142 collider_density_total = collider_density[0]+ collider_density[1] + collider_density[4] + (float)dense.eden;
00143
00144 if(mole.nH2_TRACE >= mole.nH2_trace_full )
00145 {
00146 fprintf(ioQQQ," Collider densities are:");
00147 for( nColl=0; nColl<N_X_COLLIDER; ++nColl )
00148 {
00149 fprintf(ioQQQ,"\t%.3e", collider_density[nColl]);
00150 }
00151 fprintf(ioQQQ,"\n");
00152 }
00153
00154 H2_X_source[0] = 0.;
00155 H2_X_sink[0] = 0.;
00156 for( ipHi=1; ipHi<nLevels_per_elec[0]; ++ipHi )
00157 {
00158 H2_X_source[ipHi] = 0.;
00159 H2_X_sink[ipHi] = 0.;
00160 for( ipLo=0; ipLo<ipHi; ++ipLo )
00161 {
00162 H2_X_coll_rate[ipHi][ipLo] = 0.;
00163 }
00164 }
00165
00166 exp_disoc = sexp(H2_DissocEnergies[0]/phycon.te_wn);
00167 exph2_big_0 = exp_disoc/SDIV(H2_Boltzmann[0][0][0]);
00168
00169 rel_pop_LTE_H2_big =SAHA/SDIV((phycon.te32*exph2_big_0))*(H2_stat[0][0][0]/(2.*2.))*3.634e-5;
00170
00171
00172
00173
00174
00175
00176
00177
00178 H2_X_source[0] += H2_X_formation[0][0];
00179
00180
00181 H2_X_sink[0] += (float)(H2_X_Hmin_back[0][0]*hmi.rel_pop_LTE_Hmin/
00182 SDIV(rel_pop_LTE_H2_big)*dense.eden);
00183
00184
00185 H2_X_sink[0] += collider_density_total *
00186 H2_coll_dissoc_rate_coef[0][0] * mole.lgColl_deexec_Calc;
00187
00188 H2_X_sink[0] += (collider_density[2]+ collider_density[3])*
00189 H2_coll_dissoc_rate_coef_H2[0][0] * mole.lgColl_deexec_Calc;
00190
00191
00192 H2_X_sink[0] += secondaries.csupra[ipHYDROGEN][0]*2.02f * mole.lgColl_deexec_Calc;
00193
00194
00195 H2_X_sink[0] += (float)(3.9e-21 * hextra.cryden_ov_background * co.lgUMISTrates);
00196
00197
00198 H2_X_sink[0] += (float)(2.2e-19 * hextra.cryden_ov_background * co.lgUMISTrates);
00199
00200
00201 H2_X_sink[0] += (float)(secondaries.csupra[ipHYDROGEN][0]*0.0478);
00202
00203
00204
00205
00206 H2_X_sink[0] += 3.f*secondaries.x12tot;
00207
00208
00209 H2_X_sink[0] += (float)hmi.H2_photoionize_rate;
00210
00211
00212 H2_X_sink[0] += rfield.flux_accum[H2_ipPhoto[0][0]-1]*H2_DISS_ALLISON_DALGARNO;
00213
00214
00215
00216
00217
00218 H2_X_sink[0] += (float)(hmi.rh2h2p*dense.xIonDense[ipHYDROGEN][1]);
00219
00220
00221 ipHi = nLevels_per_elec[0];
00222 while( (--ipHi) > 0 )
00223 {
00224
00225 ip = H2_ipX_ener_sort[ipHi];
00226 iVibHi = ipVib_H2_energy_sort[ip];
00227 iRotHi = ipRot_H2_energy_sort[ip];
00228
00229
00230 exph2_big = exp_disoc/SDIV(H2_Boltzmann[0][iVibHi][iRotHi]);
00231
00232
00233 rel_pop_LTE_H2_big = SAHA/SDIV((phycon.te32*exph2_big))*(H2_stat[0][iVibHi][iRotHi]/(2.*2.))*3.634e-5;
00234
00235
00236
00237 H2_X_source[ipHi] += H2_X_formation[iVibHi][iRotHi];
00238
00239 H2_X_sink[ipHi] += (float)(H2_X_Hmin_back[iVibHi][iRotHi]*hmi.rel_pop_LTE_Hmin/
00240 SDIV(rel_pop_LTE_H2_big)*dense.eden);
00241
00242
00243 H2_X_sink[ipHi] += collider_density_total *
00244 H2_coll_dissoc_rate_coef[iVibHi][iRotHi] * mole.lgColl_deexec_Calc;
00245
00246
00247 H2_X_sink[ipHi] += (collider_density[2]+ collider_density[3])*
00248 H2_coll_dissoc_rate_coef_H2[iVibHi][iRotHi] * mole.lgColl_deexec_Calc;
00249
00250
00251 H2_X_sink[ipHi] += secondaries.csupra[ipHYDROGEN][0]*2.02f * mole.lgColl_deexec_Calc;
00252
00253
00254 H2_X_sink[ipHi] += (float)(3.9e-21 * hextra.cryden_ov_background * co.lgUMISTrates);
00255
00256
00257 H2_X_sink[ipHi] += (float)(2.2e-19 * hextra.cryden_ov_background * co.lgUMISTrates);
00258
00259
00260 H2_X_sink[ipHi] += (float)(secondaries.csupra[ipHYDROGEN][0]*0.0478);
00261
00262
00263
00264
00265 H2_X_sink[ipHi] += 3.f*secondaries.x12tot;
00266
00267
00268 H2_X_sink[ipHi] += (float)hmi.H2_photoionize_rate;
00269
00270
00271 H2_X_sink[ipHi] += rfield.flux_accum[H2_ipPhoto[iVibHi][iRotHi]-1]*H2_DISS_ALLISON_DALGARNO;
00272
00273
00274 for( ipLo=0; ipLo<ipHi; ++ipLo )
00275 {
00276
00277 ip = H2_ipX_ener_sort[ipLo];
00278 iVibLo = ipVib_H2_energy_sort[ip];
00279 iRotLo = ipRot_H2_energy_sort[ip];
00280
00281
00282 colldown = 0.;
00283 for( nColl=0; nColl<N_X_COLLIDER; ++nColl )
00284 {
00285
00286 float colldown1 =
00287 H2_CollRate[nColl][iVibHi][iRotHi][iVibLo][iRotLo]*collider_density[nColl] *
00288 mole.lgColl_deexec_Calc;
00289 ASSERT( colldown1>=0. );
00290 colldown += colldown1;
00291 }
00292
00293 H2_X_coll_rate[ipHi][ipLo] += colldown;
00294
00295 }
00296 }
00297
00298 DEBUG_EXIT( "H2_X_coll_rate_evaluate()" );
00299 return;
00300 }
00301
00302
00303 double H2_itrzn( void )
00304 {
00305 if( h2.lgH2ON && nH2_zone>0 )
00306 {
00307 return( (double)nH2_pops / (double)nH2_zone );
00308 }
00309 else
00310 {
00311 return 0.;
00312 }
00313 }
00314
00315
00316 void H2_ContPoint( void )
00317 {
00318 long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo;
00319
00320 if( !h2.lgH2ON )
00321 return;
00322
00323 DEBUG_ENTRY( "H2_ContPoint()" );
00324
00325
00326 for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi )
00327 {
00328 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
00329 {
00330 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
00331 {
00332
00333
00334
00335
00336 long int lim_elec_lo = 0;
00337 for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
00338 {
00339
00340
00341 long int nv = h2.nVib_hi[iElecLo];
00342 if( iElecLo==iElecHi )
00343 nv = iVibHi;
00344 for( iVibLo=0; iVibLo<=nv; ++iVibLo )
00345 {
00346 long nr = h2.nRot_hi[iElecLo][iVibLo];
00347 if( iElecLo==iElecHi && iVibHi==iVibLo )
00348 nr = iRotHi-1;
00349
00350 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
00351 {
00352
00353
00354
00355 if( lgH2_line_exists[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] )
00356
00357 {
00358 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].ipCont =
00359 (int)ipLineEnergy(
00360 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyWN * WAVNRYD ,
00361 "H2 " , 0 );
00362 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].ipFine =
00363 ipFineCont(
00364 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyWN * WAVNRYD );
00365 }
00366 }
00367 }
00368 }
00369 }
00370 }
00371 }
00372 DEBUG_EXIT( "H2_ContPoint()" );
00373 return;
00374 }
00375
00376
00377
00378 double H2_Accel(void)
00379 {
00380 long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo;
00381 double h2_drive;
00382
00383
00384 if( !h2.lgH2ON )
00385 return 0.;
00386
00387 DEBUG_ENTRY( "H2_Accel()" );
00388
00389
00390
00391
00392 h2_drive = 0.;
00393
00394 for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi )
00395 {
00396 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
00397 {
00398 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
00399 {
00400
00401
00402
00403
00404 long int lim_elec_lo = 0;
00405 for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
00406 {
00407
00408
00409 long int nv = h2.nVib_hi[iElecLo];
00410 if( iElecLo==iElecHi )
00411 nv = iVibHi;
00412 for( iVibLo=0; iVibLo<=nv; ++iVibLo )
00413 {
00414 long nr = h2.nRot_hi[iElecLo][iVibLo];
00415 if( iElecLo==iElecHi && iVibHi==iVibLo )
00416 nr = iRotHi-1;
00417
00418 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
00419 {
00420
00421
00422 if( lgH2_line_exists[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] )
00423
00424 {
00425 h2_drive += H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].pump*
00426 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyErg*
00427 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].PopOpc;
00428 }
00429 }
00430 }
00431 }
00432 }
00433 }
00434 }
00435 DEBUG_EXIT( "H2_Accel()" );
00436 return h2_drive;
00437 }
00438
00439
00440
00441 double H2_RadPress(void)
00442 {
00443 long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo;
00444 double press;
00445
00446
00447 float smallfloat=SMALLFLOAT*10.f;
00448
00449
00450 if( !h2.lgH2ON )
00451 return 0.;
00452
00453 DEBUG_ENTRY( "H2_RadPress()" );
00454
00455 press = 0.;
00456
00457 for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi )
00458 {
00459 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
00460 {
00461 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
00462 {
00463
00464
00465
00466
00467 long int lim_elec_lo = 0;
00468 for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
00469 {
00470
00471
00472 long int nv = h2.nVib_hi[iElecLo];
00473 if( iElecLo==iElecHi )
00474 nv = iVibHi;
00475 for( iVibLo=0; iVibLo<=nv; ++iVibLo )
00476 {
00477 long nr = h2.nRot_hi[iElecLo][iVibLo];
00478 if( iElecLo==iElecHi && iVibHi==iVibLo )
00479 nr = iRotHi-1;
00480
00481 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
00482 {
00483
00484
00485
00486 if( lgH2_line_exists[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] )
00487 {
00488 if( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].PopHi > smallfloat &&
00489 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].PopOpc > smallfloat )
00490 {
00491 double RadPres1 = 5.551e-2*(
00492 powi(H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyWN/1.e6,4))*
00493 (H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].PopHi/
00494 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].gHi)/
00495 (H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].PopLo/
00496 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].gLo)*
00497 RT_LineWidth(&H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo]);
00498
00499 press += RadPres1;
00500 }
00501 }
00502 }
00503 }
00504 }
00505 }
00506 }
00507 }
00508
00509 if(mole.nH2_TRACE >= mole.nH2_trace_full)
00510 fprintf(ioQQQ,
00511 " H2_RadPress returns, radiation pressure is %.2e\n",
00512 press );
00513
00514 DEBUG_EXIT( "H2_RadPress()" );
00515
00516 return press;
00517 }
00518
00519
00520
00521 double H2_InterEnergy(void)
00522 {
00523 long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo;
00524 double energy;
00525
00526
00527 if( !h2.lgH2ON )
00528 return 0.;
00529
00530 DEBUG_ENTRY( "H2_InterEnergy()" );
00531
00532 energy = 0.;
00533
00534 for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi )
00535 {
00536 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
00537 {
00538 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
00539 {
00540
00541
00542
00543
00544 long int lim_elec_lo = 0;
00545 for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
00546 {
00547
00548
00549 long int nv = h2.nVib_hi[iElecLo];
00550 if( iElecLo==iElecHi )
00551 nv = iVibHi;
00552 for( iVibLo=0; iVibLo<=nv; ++iVibLo )
00553 {
00554 long nr = h2.nRot_hi[iElecLo][iVibLo];
00555 if( iElecLo==iElecHi && iVibHi==iVibLo )
00556 nr = iRotHi-1;
00557
00558 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
00559 {
00560
00561
00562
00563 if( lgH2_line_exists[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] )
00564 {
00565 energy +=
00566 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].PopHi*
00567 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyErg;
00568 }
00569 }
00570 }
00571 }
00572 }
00573 }
00574 }
00575
00576 DEBUG_EXIT( "H2_InterEnergy()" );
00577
00578 return energy;
00579 }
00580
00581
00582 void H2_RT_diffuse(void)
00583 {
00584 long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo;
00585
00586 if( !h2.lgH2ON || !h2.nCallH2_this_zone )
00587 return;
00588
00589 DEBUG_ENTRY( "H2_RT_diffuse()" );
00590
00591
00592
00593 for( iElecHi=0; iElecHi<1; ++iElecHi )
00594 {
00595 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
00596 {
00597 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
00598 {
00599
00600
00601
00602
00603 long int lim_elec_lo = 0;
00604 for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
00605 {
00606
00607
00608 long int nv = h2.nVib_hi[iElecLo];
00609 if( iElecLo==iElecHi )
00610 nv = iVibHi;
00611 for( iVibLo=0; iVibLo<=nv; ++iVibLo )
00612 {
00613 long nr = h2.nRot_hi[iElecLo][iVibLo];
00614 if( iElecLo==iElecHi && iVibHi==iVibLo )
00615 nr = iRotHi-1;
00616
00617 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
00618 {
00619
00620
00621
00622
00623 if( lgH2_line_exists[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] )
00624 {
00625 outline( &H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] );
00626 }
00627 }
00628 }
00629 }
00630 }
00631 }
00632 }
00633
00634 DEBUG_EXIT( "H2_RT_diffuse()" );
00635
00636 return;
00637 }
00638
00639
00640 void H2_RTMake( bool lgDoEsc , bool lgUpdateFineOpac )
00641 {
00642 long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo;
00643
00644 if( !h2.lgH2ON )
00645 return;
00646
00647 DEBUG_ENTRY( "H2_RTMake()" );
00648
00649
00650
00651 for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi )
00652 {
00653 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
00654 {
00655 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
00656 {
00657
00658
00659
00660
00661 long int lim_elec_lo = 0;
00662 for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
00663 {
00664
00665
00666 long int nv = h2.nVib_hi[iElecLo];
00667 if( iElecLo==iElecHi )
00668 nv = iVibHi;
00669 for( iVibLo=0; iVibLo<=nv; ++iVibLo )
00670 {
00671 long nr = h2.nRot_hi[iElecLo][iVibLo];
00672 if( iElecLo==iElecHi && iVibHi==iVibLo )
00673 nr = iRotHi-1;
00674
00675 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
00676 {
00677
00678
00679
00680 if( lgH2_line_exists[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] )
00681 {
00682
00683
00684
00685
00686 RT_line_one( &H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] , lgDoEsc , lgUpdateFineOpac,false);
00687 }
00688 }
00689 }
00690 }
00691 }
00692 }
00693 }
00694
00695 {
00696
00697 enum {DEBUG_LOC=false};
00698
00699 if( DEBUG_LOC )
00700 {
00701 double sumpop = 0.;
00702 double sumpopA = 0.;
00703 for( iElecHi=1; iElecHi<mole.n_h2_elec_states; ++iElecHi )
00704 {
00705 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
00706 {
00707 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
00708 {
00709
00710
00711
00712
00713 iElecLo = 0;
00714
00715
00716 for( iVibLo=0; iVibLo<=h2.nVib_hi[iElecLo]; ++iVibLo )
00717 {
00718 long nr = h2.nRot_hi[iElecLo][iVibLo];
00719 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
00720 {
00721
00722
00723
00724 if( lgH2_line_exists[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] )
00725 {
00726 sumpop += H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].PopLo;
00727 sumpopA += H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].PopLo*
00728 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Aul;
00729 }
00730 }
00731 }
00732 }
00733 }
00734 }
00735 fprintf(ioQQQ,"sumpop = %.3e sumpopA= %.3e A=%.3e\n", sumpop, sumpopA,
00736 sumpopA/SDIV(sumpop) );
00737 }
00738 }
00739 DEBUG_EXIT( "H2_RTMake()" );
00740 return;
00741 }
00742
00743
00744
00745 void H2_RT_tau_inc(void)
00746 {
00747 long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo;
00748
00749
00750
00751 if( !h2.lgH2ON )
00752 return;
00753
00754 DEBUG_ENTRY( "H2_RT_tau_inc()" );
00755
00756
00757
00758
00759 if( nzone > 0 && nCallH2_this_iteration>2 )
00760 {
00761 h2.renorm_max = MAX2( H2_renorm_chemistry , h2.renorm_max );
00762 h2.renorm_min = MIN2( H2_renorm_chemistry , h2.renorm_min );
00763 }
00764
00765
00766 for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi )
00767 {
00768 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
00769 {
00770 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
00771 {
00772
00773
00774
00775
00776 long int lim_elec_lo = 0;
00777 for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
00778 {
00779
00780
00781 long int nv = h2.nVib_hi[iElecLo];
00782 if( iElecLo==iElecHi )
00783 nv = iVibHi;
00784 for( iVibLo=0; iVibLo<=nv; ++iVibLo )
00785 {
00786 long nr = h2.nRot_hi[iElecLo][iVibLo];
00787 if( iElecLo==iElecHi && iVibHi==iVibLo )
00788 nr = iRotHi-1;
00789
00790 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
00791 {
00792
00793
00794
00795
00796 if( lgH2_line_exists[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] )
00797 {
00798
00799
00800 RT_line_one_tauinc( &H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] ,
00801 -9 , iRotHi , iVibLo , iRotLo);
00802
00803 }
00804 }
00805 }
00806 }
00807 }
00808 }
00809 }
00810
00811
00812 DEBUG_EXIT( "H2_RT_tau_inc()" );
00813
00814 return;
00815 }
00816
00817
00818
00819 void H2_LineZero( void )
00820 {
00821 long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo;
00822
00823 if( !h2.lgH2ON )
00824 return;
00825
00826 DEBUG_ENTRY( "H2_LineZero()" );
00827
00828
00829 for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi )
00830 {
00831 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
00832 {
00833 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
00834 {
00835
00836
00837
00838
00839 long int lim_elec_lo = 0;
00840 for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
00841 {
00842
00843
00844 long int nv = h2.nVib_hi[iElecLo];
00845 if( iElecLo==iElecHi )
00846 nv = iVibHi;
00847 for( iVibLo=0; iVibLo<=nv; ++iVibLo )
00848 {
00849 long nr = h2.nRot_hi[iElecLo][iVibLo];
00850 if( iElecLo==iElecHi && iVibHi==iVibLo )
00851 nr = iRotHi-1;
00852
00853 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
00854 {
00855
00856
00857
00858
00859 if( lgH2_line_exists[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] )
00860 {
00861
00862 EmLineZero( &H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] );
00863 }
00864 }
00865 }
00866 }
00867 }
00868 }
00869 }
00870 DEBUG_EXIT( "H2_LineZero()" );
00871 return;
00872 }
00873
00874
00875 void H2_RT_tau_reset( void )
00876 {
00877 long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo;
00878
00879 if( !h2.lgH2ON )
00880 return;
00881
00882 DEBUG_ENTRY( "H2_RT_tau_reset()" );
00883
00884
00885 for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi )
00886 {
00887 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
00888 {
00889 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
00890 {
00891
00892
00893
00894
00895 long int lim_elec_lo = 0;
00896 for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
00897 {
00898
00899
00900 long int nv = h2.nVib_hi[iElecLo];
00901 if( iElecLo==iElecHi )
00902 nv = iVibHi;
00903 for( iVibLo=0; iVibLo<=nv; ++iVibLo )
00904 {
00905 long nr = h2.nRot_hi[iElecLo][iVibLo];
00906 if( iElecLo==iElecHi && iVibHi==iVibLo )
00907 nr = iRotHi-1;
00908
00909 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
00910 {
00911
00912
00913
00914 if( lgH2_line_exists[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] )
00915 {
00916
00917 RT_line_one_tau_reset( &H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] , 0.5);
00918 }
00919 }
00920 }
00921 }
00922 }
00923 }
00924 }
00925 DEBUG_EXIT( "H2_RT_tau_reset()" );
00926 return;
00927 }
00928
00929
00930 static double frac_matrix;
00931
00932
00933 static void H2_Level_low_matrix(
00934
00935 float abundance )
00936 {
00937
00938
00939 static double **AulEscp ,
00940 **col_str ,
00941 **AulDest,
00942
00943 **AulPump,
00944 **CollRate_levn,
00945 *pops,
00946 *create,
00947 *destroy,
00948 *depart,
00949
00950 *stat_levn ,
00951
00952 *excit;
00953 bool lgDoAs;
00954 static long int levelAsEval=-1;
00955 long int ip;
00956 static bool lgFirst=true;
00957 long int i,
00958 j,
00959 ilo ,
00960 ihi,
00961 iElec,
00962 iElecHi,
00963 iVib,
00964 iRot,
00965 iVibHi,
00966 iRotHi;
00967 int nNegPop;
00968 bool lgDeBug,
00969 lgZeroPop;
00970 double rot_cooling , dCoolDT;
00971 static long int ndimMalloced = 0;
00972 double rateout , ratein;
00973
00974 DEBUG_ENTRY( "H2_Level_low_matrix()" );
00975
00976
00977 if( nXLevelsMatrix <= 1 )
00978 {
00979 DEBUG_EXIT( "H2_Level_low_matrix()" );
00980 return;
00981 }
00982
00983 if( lgFirst )
00984 {
00985
00986 if( nXLevelsMatrix > nLevels_per_elec[0] )
00987 {
00988
00989 fprintf( ioQQQ,
00990 " The total number of levels used in the matrix solver must be <= %li, the number of levels within X.\n Sorry.\n",
00991 nLevels_per_elec[0]);
00992 puts( "[Stop in ParseAtomH2]" );
00993 cdEXIT(EXIT_FAILURE);
00994 }
00995
00996 lgFirst = false;
00997
00998
00999
01000 ndimMalloced = nLevels_per_elec[0];
01001
01002 excit = (double *)MALLOC( sizeof(double)*(size_t)(ndimMalloced) );
01003 stat_levn = (double *)MALLOC( sizeof(double)*(size_t)(ndimMalloced) );
01004 pops = (double *)MALLOC( sizeof(double)*(size_t)(ndimMalloced) );
01005 create = (double *)MALLOC( sizeof(double)*(size_t)(ndimMalloced) );
01006 destroy = (double *)MALLOC( sizeof(double)*(size_t)(ndimMalloced) );
01007 depart = (double *)MALLOC( sizeof(double)*(size_t)(ndimMalloced) );
01008
01009 AulPump = ((double **)MALLOC((size_t)(ndimMalloced)*sizeof(double *)));
01010 CollRate_levn = ((double **)MALLOC((size_t)(ndimMalloced)*sizeof(double *)));
01011 AulDest = ((double **)MALLOC((size_t)(ndimMalloced)*sizeof(double *)));
01012 AulEscp = ((double **)MALLOC((size_t)(ndimMalloced)*sizeof(double *)));
01013 col_str = ((double **)MALLOC((size_t)(ndimMalloced)*sizeof(double *)));
01014 for( i=0; i<(ndimMalloced); ++i )
01015 {
01016 AulPump[i] = ((double *)MALLOC((size_t)(ndimMalloced)*sizeof(double )));
01017 CollRate_levn[i] = ((double *)MALLOC((size_t)(ndimMalloced)*sizeof(double )));
01018 AulDest[i] = ((double *)MALLOC((size_t)(ndimMalloced)*sizeof(double )));
01019 AulEscp[i] = ((double *)MALLOC((size_t)(ndimMalloced)*sizeof(double )));
01020 col_str[i] = ((double *)MALLOC((size_t)(ndimMalloced)*sizeof(double )));
01021 }
01022
01023 for( j=0; j < ndimMalloced; j++ )
01024 {
01025 stat_levn[j]=0;
01026 excit[j] =0;
01027 }
01028
01029
01030 for( j=0; j < ndimMalloced; j++ )
01031 {
01032
01033 ip = H2_ipX_ener_sort[j];
01034 iVib = ipVib_H2_energy_sort[ip];
01035 iRot = ipRot_H2_energy_sort[ip];
01036
01037
01038 stat_levn[j] = H2_stat[0][iVib][iRot];
01039
01040 excit[j] = energy_wn[0][iVib][iRot]*T1CM;
01041 }
01042
01043 for( j=0; j < ndimMalloced-1; j++ )
01044 {
01045
01046 ASSERT( excit[j+1] > excit[j] );
01047 }
01048 }
01049
01050
01051
01052
01053 if( nXLevelsMatrix > ndimMalloced )
01054 {
01055 fprintf(ioQQQ," H2_Level_low_matrix has been called with the number of rotor levels greater than space allocated.\n");
01056 puts( "[Stop in H2_Level_low_matrix]" );
01057 cdEXIT(EXIT_FAILURE);
01058 }
01059
01060
01061 for( i=0; i < nXLevelsMatrix; i++ )
01062 {
01063 pops[i] = 0.;
01064 depart[i] = 0;
01065 for( j=0; j < nXLevelsMatrix; j++ )
01066 {
01067 col_str[j][i] = 0.;
01068 }
01069 }
01070
01071
01072 if( nzone!=nzoneAsEval || iteration!=iterationAsEval || nXLevelsMatrix!=levelAsEval)
01073 {
01074 lgDoAs = true;
01075 nzoneAsEval = nzone;
01076 iterationAsEval = iteration;
01077 levelAsEval = nXLevelsMatrix;
01078 ASSERT( levelAsEval <= ndimMalloced );
01079 }
01080 else
01081 {
01082 lgDoAs = false;
01083 }
01084
01085
01086 if( lgDoAs )
01087 {
01088 for( i=0; i < nXLevelsMatrix; i++ )
01089 {
01090 pops[i] = 0.;
01091 depart[i] = 0;
01092 for( j=0; j < nXLevelsMatrix; j++ )
01093 {
01094 AulEscp[j][i] = 0.;
01095 AulDest[j][i] = 0.;
01096 AulPump[j][i] = 0.;
01097 CollRate_levn[j][i] = 0.;
01098 }
01099 }
01100 }
01101
01102
01103
01104 iElec = 0;
01105 for( ilo=0; ilo < nXLevelsMatrix; ilo++ )
01106 {
01107 ip = H2_ipX_ener_sort[ilo];
01108 iRot = ipRot_H2_energy_sort[ip];
01109 iVib = ipVib_H2_energy_sort[ip];
01110
01111
01112
01113
01114
01115 destroy[ilo] = H2_X_sink[ilo];
01116
01117
01118 create[ilo] = H2_X_source[ilo];
01119
01120
01121
01122 if( lgDoAs )
01123 {
01124 for( ihi=ilo+1; ihi<nXLevelsMatrix; ++ihi )
01125 {
01126 ip = H2_ipX_ener_sort[ihi];
01127 iRotHi = ipRot_H2_energy_sort[ip];
01128 iVibHi = ipVib_H2_energy_sort[ip];
01129 ASSERT( H2_Xenergies[ip] <= H2_Xenergies[H2_ipX_ener_sort[nXLevelsMatrix-1]] );
01130
01131 if( (abs(iRotHi-iRot)==2 || (iRotHi-iRot)==0 ) && (iVib<=iVibHi) )
01132 {
01133
01134 if( lgH2_line_exists[0][iVibHi][iRotHi][0][iVib][iRot] )
01135 {
01136
01137
01138
01139 AulEscp[ihi][ilo] = H2Lines[0][iVibHi][iRotHi][0][iVib][iRot].Aul*(
01140 H2Lines[0][iVibHi][iRotHi][0][iVib][iRot].Pesc +
01141 H2Lines[0][iVibHi][iRotHi][0][iVib][iRot].Pdest +
01142 H2Lines[0][iVibHi][iRotHi][0][iVib][iRot].Pelec_esc);
01143 AulDest[ilo][ihi] = 0.;
01144 AulPump[ilo][ihi] = H2Lines[0][iVibHi][iRotHi][0][iVib][iRot].pump;
01145 }
01146 }
01147 }
01148 }
01149
01150 iElecHi = 0;
01151 iElec = 0;
01152 rateout = 0.;
01153 ratein = 0.;
01154
01155
01156 for( ihi=nXLevelsMatrix; ihi<nLevels_per_elec[0]; ++ihi )
01157 {
01158 ip = H2_ipX_ener_sort[ihi];
01159 iRotHi = ipRot_H2_energy_sort[ip];
01160 iVibHi = ipVib_H2_energy_sort[ip];
01161 if( (abs(iRotHi-iRot)==2 || (iRotHi-iRot)==0 ) && (iVib<=iVibHi) )
01162 {
01163 if( lgH2_line_exists[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot] )
01164 {
01165
01166
01167
01168 ratein +=
01169 H2_populations[iElecHi][iVibHi][iRotHi] *
01170 (H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Aul*
01171 (H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Pesc +
01172 H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Pelec_esc +
01173 H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Pdest)+H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].pump *
01174 H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].gLo/
01175 H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].gHi);
01176
01177 rateout +=
01178 H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].pump;
01179 }
01180 }
01181 }
01182
01183
01184 create[ilo] += ratein;
01185
01186
01187 destroy[ilo] += rateout;
01188
01189
01190
01191 create[ilo] += H2_X_rate_from_elec_excited[iVib][iRot];
01192
01193
01194 destroy[ilo] += H2_X_rate_to_elec_excited[iVib][iRot];
01195
01196 }
01197
01198
01199 if( mole.nH2_TRACE >= mole.nH2_trace_matrix )
01200 lgDeBug = true;
01201 else
01202 lgDeBug = false;
01203
01204
01205
01206
01207
01208 for( ilo=0; ilo < nXLevelsMatrix; ilo++ )
01209 {
01210 ip = H2_ipX_ener_sort[ilo];
01211 iRot = ipRot_H2_energy_sort[ip];
01212 iVib = ipVib_H2_energy_sort[ip];
01213 if( lgDoAs )
01214 {
01215 if(lgDeBug)fprintf(ioQQQ,"DEBUG H2_Level_low_matrix, ilo=%li",ilo);
01216 for( ihi=ilo+1; ihi < nXLevelsMatrix; ihi++ )
01217 {
01218 ip = H2_ipX_ener_sort[ihi];
01219 iRotHi = ipRot_H2_energy_sort[ip];
01220 iVibHi = ipVib_H2_energy_sort[ip];
01221
01222 CollRate_levn[ihi][ilo] = H2_X_coll_rate[ihi][ilo];
01223
01224
01225 if(lgDeBug)fprintf(ioQQQ,"\t%.1e",CollRate_levn[ihi][ilo]);
01226
01227
01228 CollRate_levn[ilo][ihi] = CollRate_levn[ihi][ilo]*
01229 H2_Boltzmann[0][iVibHi][iRotHi]/SDIV(H2_Boltzmann[0][iVib][iRot])*
01230 H2_stat[0][iVibHi][iRotHi] /
01231 H2_stat[0][iVib][iRot];
01232 }
01233 }
01234
01235 if(lgDeBug)fprintf(ioQQQ,"\n");
01236
01237
01238
01239 iElecHi = 0;
01240
01241 for( ihi=nXLevelsMatrix; ihi<nLevels_per_elec[0]; ++ihi )
01242 {
01243 ip = H2_ipX_ener_sort[ihi];
01244 iRotHi = ipRot_H2_energy_sort[ip];
01245 iVibHi = ipVib_H2_energy_sort[ip];
01246 rateout = 0;
01247
01248
01249
01250 ratein = H2_X_coll_rate[ihi][ilo];
01251 if(lgDeBug)fprintf(ioQQQ,"\t%.1e",ratein);
01252
01253
01254 rateout = ratein *
01255 H2_Boltzmann[0][iVibHi][iRotHi]/SDIV(H2_Boltzmann[0][iVib][iRot])*
01256 H2_stat[0][iVibHi][iRotHi]/H2_stat[0][iVib][iRot];
01257
01258
01259 create[ilo] += ratein*H2_populations[iElecHi][iVibHi][iRotHi];
01260 destroy[ilo] += rateout;
01261 }
01262 }
01263
01264
01265
01266
01267 if( lgDoAs )
01268 {
01269 for( ihi=2; ihi < nXLevelsMatrix; ihi++ )
01270 {
01271
01272 ip = H2_ipX_ener_sort[ihi];
01273 iVibHi = ipVib_H2_energy_sort[ip];
01274 iRotHi = ipRot_H2_energy_sort[ip];
01275
01276
01277
01278
01279
01280 CollRate_levn[ihi][H2_lgOrtho[0][iVibHi][iRotHi]] += hmi.rate_grain_h2_op_conserve;
01281 }
01282
01283
01284
01285 CollRate_levn[1][0] +=
01286 (float)(hmi.rate_grain_h2_J1_to_J0);
01287 }
01288
01289
01290 for( ihi=nXLevelsMatrix; ihi<nLevels_per_elec[0]; ++ihi )
01291 {
01292 ip = H2_ipX_ener_sort[ihi];
01293 iVibHi = ipVib_H2_energy_sort[ip];
01294 iRotHi = ipRot_H2_energy_sort[ip];
01295
01296
01297
01298 create[H2_lgOrtho[0][iVibHi][iRotHi]] += H2_populations[0][iVibHi][iRotHi]*hmi.rate_grain_h2_op_conserve;
01299 }
01300
01301 {
01302
01303 enum {DEBUG_LOC=false};
01304
01305 if( DEBUG_LOC || lgDeBug)
01306 {
01307 fprintf(ioQQQ,"DEBUG H2 matexcit");
01308 for(ilo=0; ilo<nXLevelsMatrix; ++ilo )
01309 {
01310 fprintf(ioQQQ,"\t%li",ilo );
01311 }
01312 fprintf(ioQQQ,"\n");
01313 for(ihi=0; ihi<nXLevelsMatrix;++ihi)
01314 {
01315 fprintf(ioQQQ,"\t%.2e",excit[ihi] );
01316 }
01317 fprintf(ioQQQ,"\n");
01318 for(ihi=0; ihi<nXLevelsMatrix;++ihi)
01319 {
01320 fprintf(ioQQQ,"\t%.2e",stat_levn[ihi] );
01321 }
01322 fprintf(ioQQQ,"\n");
01323
01324 fprintf(ioQQQ,"AulEscp[n][]\\[][n] = Aul*Pesc\n");
01325 for(ilo=0; ilo<nXLevelsMatrix; ++ilo )
01326 {
01327 fprintf(ioQQQ,"\t%li",ilo );
01328 }
01329 fprintf(ioQQQ,"\n");
01330 for(ihi=0; ihi<nXLevelsMatrix;++ihi)
01331 {
01332 fprintf(ioQQQ,"%li", ihi);
01333 for(ilo=0; ilo<nXLevelsMatrix; ++ilo )
01334 {
01335 fprintf(ioQQQ,"\t%.2e",AulEscp[ilo][ihi] );
01336 }
01337 fprintf(ioQQQ,"\n");
01338 }
01339
01340 fprintf(ioQQQ,"AulPump [n][]\\[][n]\n");
01341 for(ilo=0; ilo<nXLevelsMatrix; ++ilo )
01342 {
01343 fprintf(ioQQQ,"\t%li",ilo );
01344 }
01345 fprintf(ioQQQ,"\n");
01346 for(ihi=0; ihi<nXLevelsMatrix;++ihi)
01347 {
01348 fprintf(ioQQQ,"%li", ihi);
01349 for(ilo=0; ilo<nXLevelsMatrix; ++ilo )
01350 {
01351 fprintf(ioQQQ,"\t%.2e",AulPump[ihi][ilo] );
01352 }
01353 fprintf(ioQQQ,"\n");
01354 }
01355
01356 fprintf(ioQQQ,"CollRate_levn [n][]\\[][n]\n");
01357 for(ilo=0; ilo<nXLevelsMatrix; ++ilo )
01358 {
01359 fprintf(ioQQQ,"\t%li",ilo );
01360 }
01361 fprintf(ioQQQ,"\n");
01362 for(ihi=0; ihi<nXLevelsMatrix;++ihi)
01363 {
01364 fprintf(ioQQQ,"%li", ihi);
01365 for(ilo=0; ilo<nXLevelsMatrix; ++ilo )
01366 {
01367 fprintf(ioQQQ,"\t%.2e",CollRate_levn[ihi][ilo] );
01368 }
01369 fprintf(ioQQQ,"\n");
01370 }
01371 fprintf(ioQQQ,"SOURCE");
01372 for(ihi=0; ihi<nXLevelsMatrix;++ihi)
01373 {
01374 fprintf(ioQQQ,"\t%.2e",create[ihi]);
01375 }
01376 fprintf(ioQQQ,"\nSINK");
01377 for(ihi=0; ihi<nXLevelsMatrix;++ihi)
01378 {
01379 fprintf(ioQQQ,"\t%.2e",destroy[ihi]);
01380 }
01381 fprintf(ioQQQ,"\n");
01382 }
01383 }
01384
01385 atom_levelN(
01386
01387 nXLevelsMatrix,
01388 abundance,
01389 stat_levn,
01390 excit,
01391 'K',
01392 pops,
01393 depart,
01394
01395 &AulEscp,
01396
01397 &col_str,
01398 &AulDest,
01399 &AulPump,
01400 &CollRate_levn,
01401 create,
01402 destroy,
01403
01404 true,
01405
01406 &rot_cooling,
01407 &dCoolDT,
01408 " H2 ",
01409
01410 &nNegPop,
01411 &lgZeroPop,
01412 lgDeBug );
01413
01414 for( i=0; i< nXLevelsMatrix; ++i )
01415 {
01416 ip = H2_ipX_ener_sort[i];
01417 iRot = ipRot_H2_energy_sort[ip];
01418 iVib = ipVib_H2_energy_sort[ip];
01419
01420
01421
01422
01423 H2_populations[0][iVib][iRot] = pops[i];
01424 }
01425
01426 if( 0 && mole.nH2_TRACE >= mole.nH2_trace_full)
01427 {
01428
01429
01430 fprintf(ioQQQ,"\n DEBUG H2_Level_lowJ hmi.H2_total: %.3e matrix rel pops\n",hmi.H2_total);
01431 fprintf(ioQQQ,"v\tJ\tpop\n");
01432 for( i=0; i<nXLevelsMatrix; ++i )
01433 {
01434 ip = H2_ipX_ener_sort[i];
01435 iRot = ipRot_H2_energy_sort[ip];
01436 iVib = ipVib_H2_energy_sort[ip];
01437 fprintf(ioQQQ,"%3li\t%3li\t%.3e\t%.3e\t%.3e\n",
01438 iVib , iRot , H2_populations[0][iVib][iRot]/hmi.H2_total , create[i] , destroy[i]);
01439 }
01440 }
01441
01442
01443 if( nNegPop > 0 )
01444 {
01445 fprintf(ioQQQ," H2_Level_low_matrix called atom_levelN which returned negative H2_populations.\n");
01446 ConvFail( "pops" , "H2" );
01447 }
01448
01449 DEBUG_EXIT( "H2_Level_low_matrix()" );
01450
01451 return;
01452 }
01453
01454
01455 static void H2_collid_rates( void )
01456 {
01457 long int numb_coll_trans = 0;
01458 double excit;
01459 double t = phycon.te/1000. + 1.;
01460 double t2 = POW2(t);
01461 long int iElecHi , iElecLo , ipHi , iVibHi , iRotHi ,
01462 ipLo , iVibLo , iRotLo , nColl;
01463
01464 DEBUG_ENTRY( "H2_collid_rates()" );
01465
01466 iElecHi = 0;
01467 iElecLo = 0;
01468 if(mole.nH2_TRACE >= mole.nH2_trace_full)
01469 fprintf(ioQQQ,"H2 set collision rates\n");
01470
01471 # define PRT_COLL false
01472
01473
01474
01475
01476
01477 H2_coll_dissoc_rate_coef[0][0] = 0.;
01478 H2_coll_dissoc_rate_coef_H2[0][0] = 0.;
01479 for( ipHi=1; ipHi<nLevels_per_elec[0]; ++ipHi )
01480 {
01481 double energy;
01482
01483
01484 long int ip = H2_ipX_ener_sort[ipHi];
01485 iVibHi = ipVib_H2_energy_sort[ip];
01486 iRotHi = ipRot_H2_energy_sort[ip];
01487
01488
01489
01490 energy = H2_DissocEnergies[0] - energy_wn[0][iVibHi][iRotHi];
01491 ASSERT( energy > 0. );
01492
01493 H2_coll_dissoc_rate_coef[iVibHi][iRotHi] =
01494 1e-14f * (float)sexp(energy/phycon.te_wn) * mole.lgColl_dissoc_coll;
01495
01496
01497 H2_coll_dissoc_rate_coef_H2[iVibHi][iRotHi] =
01498 1e-11f * (float)sexp(energy/phycon.te_wn) * mole.lgColl_dissoc_coll;
01499
01500
01501
01502
01503 for( ipLo=0; ipLo<ipHi; ++ipLo )
01504 {
01505
01506 double gbarcoll[N_X_COLLIDER][3] =
01507 {
01508 {-9.9265 , -0.1048 , 0.456 },
01509 {-8.281 , -0.1303 , 0.4931 },
01510 {-10.0357, -0.0243 , 0.67 },
01511 {-8.6213 , -0.1004 , 0.5291 },
01512 {-9.2719 , -0.0001 , 1.0391 } ,
01513
01514 {-8.281 , -0.1303 , 0.4931 }};
01515 double fitted;
01516
01517 ip = H2_ipX_ener_sort[ipLo];
01518 iVibLo = ipVib_H2_energy_sort[ip];
01519 iRotLo = ipRot_H2_energy_sort[ip];
01520
01521
01522
01523
01524
01525
01526 ASSERT( energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo] > 0.);
01527
01528
01529
01530
01531
01532
01533
01534 ++numb_coll_trans;
01535
01536
01537
01538
01539 for( nColl=0; nColl<N_X_COLLIDER-2; ++nColl )
01540 {
01541
01542 if( CollRateFit[nColl][iVibHi][iRotHi][iVibLo][iRotLo][0]!= 0 )
01543 {
01544 double r = CollRateFit[nColl][iVibHi][iRotHi][iVibLo][iRotLo][0] +
01545 CollRateFit[nColl][iVibHi][iRotHi][iVibLo][iRotLo][1]/t +
01546 CollRateFit[nColl][iVibHi][iRotHi][iVibLo][iRotLo][2]/t2;
01547 H2_CollRate[nColl][iVibHi][iRotHi][iVibLo][iRotLo] =
01548 (float)pow(10.,r)*mole.lgColl_deexec_Calc;
01549 if( PRT_COLL )
01550 fprintf(ioQQQ,"col fit\t%li\t%li\t%li\t%li\t%li\t%.4e\t%.4e\n",
01551 nColl,
01552 iVibHi,iRotHi,iVibLo,iRotLo,
01553 energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo],
01554 H2_CollRate[nColl][iVibHi][iRotHi][iVibLo][iRotLo] );
01555 }
01556
01557
01558
01559 else if( mole.lgColl_gbar &&
01560 (H2_lgOrtho[0][iVibHi][iRotHi]-H2_lgOrtho[0][iVibLo][iRotLo]==0) )
01561 {
01562
01563
01564 double ediff = energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo];
01565
01566
01567 ediff = MAX2(100., ediff );
01568 H2_CollRate[nColl][iVibHi][iRotHi][iVibLo][iRotLo] =
01569 (float)pow(10. ,
01570 gbarcoll[nColl][0] + gbarcoll[nColl][1] *
01571 pow(ediff,gbarcoll[nColl][2]) )*mole.lgColl_deexec_Calc;
01572
01573 if( PRT_COLL )
01574 fprintf(ioQQQ,"col gbr\t%li\t%li\t%li\t%li\t%li\t%.4e\t%.4e\n",
01575 nColl+10,
01576 iVibHi,iRotHi,iVibLo,iRotLo,
01577 energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo],
01578 H2_CollRate[nColl][iVibHi][iRotHi][iVibLo][iRotLo] );
01579 }
01580 }
01581
01582
01583
01584
01585
01586 if( CollRateFit[N_X_COLLIDER-2][iVibHi][iRotHi][iVibLo][iRotLo][1] != 0 )
01587 {
01588 H2_CollRate[N_X_COLLIDER-2][iVibHi][iRotHi][iVibLo][iRotLo] =
01589 CollRateFit[N_X_COLLIDER-2][iVibHi][iRotHi][iVibLo][iRotLo][0] * 1e-10f *
01590
01591 (float)sexp( CollRateFit[N_X_COLLIDER-2][iVibHi][iRotHi][iVibLo][iRotLo][1]/1000./phycon.te_eV)*mole.lgColl_deexec_Calc;
01592 if( PRT_COLL )
01593 fprintf(ioQQQ,"col fit\t%i\t%li\t%li\t%li\t%li\t%.4e\t%.4e\n",
01594 N_X_COLLIDER-2,
01595 iVibHi,iRotHi,iVibLo,iRotLo,
01596 energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo],
01597 H2_CollRate[N_X_COLLIDER-2][iVibHi][iRotHi][iVibLo][iRotLo] );
01598 }
01599
01600
01601
01602 else if( mole.lgColl_gbar )
01603 {
01604
01605
01606 double ediff = energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo];
01607 ediff = MAX2(100., ediff );
01608 H2_CollRate[N_X_COLLIDER-2][iVibHi][iRotHi][iVibLo][iRotLo] =
01609 (float)pow(10. ,
01610 gbarcoll[N_X_COLLIDER-2][0] + gbarcoll[N_X_COLLIDER-2][1] *
01611 pow(ediff ,gbarcoll[N_X_COLLIDER-2][2]) )*mole.lgColl_deexec_Calc;
01612
01613 if( PRT_COLL )
01614 fprintf(ioQQQ,"col gbr\t%i\t%li\t%li\t%li\t%li\t%.4e\t%.4e\n",
01615 N_X_COLLIDER-2+10,
01616 iVibHi,iRotHi,iVibLo,iRotLo,
01617 energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo],
01618 H2_CollRate[N_X_COLLIDER-2][iVibHi][iRotHi][iVibLo][iRotLo] );
01619 }
01620
01621
01622
01623
01624
01625
01626
01627
01628 nColl = 5;
01629 if( (fitted=H2_He_coll(ipHi , ipLo , phycon.te ))>0. )
01630 {
01631
01632
01633
01634
01635 H2_CollRate[nColl][iVibHi][iRotHi][iVibLo][iRotLo] = (float)fitted*mole.lgColl_deexec_Calc;
01636 if( PRT_COLL )
01637 fprintf(ioQQQ,"col fit\t%li\t%li\t%.4e\t%li\t%li\t%li\t%li\t%li\t%.4e\t%.4e\n",
01638 ipLo,ipHi,phycon.te,nColl,
01639 iVibHi,iRotHi,iVibLo,iRotLo,
01640 energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo],
01641 H2_CollRate[nColl][iVibHi][iRotHi][iVibLo][iRotLo] );
01642 }
01643
01644
01645
01646
01647
01648
01649 else if( mole.lgColl_gbar &&
01650 (H2_lgOrtho[0][iVibHi][iRotHi]-H2_lgOrtho[0][iVibLo][iRotLo]==0) )
01651 {
01652
01653
01654 double ediff = energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo];
01655
01656
01657
01658 ediff = MAX2(100., ediff );
01659 H2_CollRate[nColl][iVibHi][iRotHi][iVibLo][iRotLo] =
01660 (float)pow(10. ,
01661 gbarcoll[nColl][0] + gbarcoll[nColl][1] *
01662 pow(ediff,gbarcoll[nColl][2]) )*mole.lgColl_deexec_Calc;
01663
01664 if( PRT_COLL )
01665 fprintf(ioQQQ,"col gbr\t%li\t%li\t%li\t%li\t%li\t%.4e\t%.4e\n",
01666 nColl+10,
01667 iVibHi,iRotHi,iVibLo,iRotLo,
01668 energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo],
01669 H2_CollRate[nColl][iVibHi][iRotHi][iVibLo][iRotLo] );
01670 }
01671
01672 if( !mole.lgH2_ortho_para_coll_on &&
01673 (H2_lgOrtho[0][iVibHi][iRotHi]-H2_lgOrtho[0][iVibLo][iRotLo]) )
01674 {
01675 H2_CollRate[nColl][iVibHi][iRotHi][iVibLo][iRotLo] = 0.;
01676 }
01677
01678 {
01679
01680 enum {DEBUG_LOC=false};
01681
01682 if( DEBUG_LOC )
01683 {
01684 fprintf(ioQQQ,"bugcoll\tiVibHi\t%li\tiRotHi\t%li\tiVibLo\t%li\tiRotLo\t%li\tcoll\t%.2e\n",
01685 iVibHi,iRotHi,iVibLo,iRotLo,
01686 H2_CollRate[nColl][iVibHi][iRotHi][iVibLo][iRotLo] );
01687 }
01688 }
01689 }
01690 }
01691
01692
01693
01694
01695
01696
01697 nColl = 0;
01698 iElecHi = 0;
01699 iElecLo = 0;
01700 iVibHi = 0;
01701 iVibLo = 0;
01702
01703
01704
01705
01706 {
01707 double excit1;
01708 double te_used = MAX2( 100. , phycon.te );
01709
01710 iRotLo = 0;
01711 iRotHi = 1;
01712 excit1 = sexp( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyK/te_used);
01713 excit = sexp( -(POW2(5.30-460./te_used)-21.2) )*1e-13;
01714
01715 H2_CollRate[0][iVibHi][iRotHi][iVibLo][iRotLo] = (float)(
01716 excit*H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].gLo/
01717 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].gHi /
01718
01719 SDIV(excit1) )*mole.lgColl_deexec_Calc *
01720
01721 mole.lgH2_ortho_para_coll_on;
01722
01723 if( PRT_COLL )
01724 fprintf(ioQQQ,"col o-p\t%li\t%li\t%li\t%li\t%li\t%.4e\t%.4e\n",
01725 nColl,
01726 iVibHi,iRotHi,iVibLo,iRotLo,
01727 energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo],
01728 H2_CollRate[nColl][iVibHi][iRotHi][iVibLo][iRotLo] );
01729
01730
01731 iRotLo = 0;
01732 iRotHi = 3;
01733 excit = sexp( -(POW2(6.36-373./te_used)-34.5) )*1e-13;
01734 excit1 = sexp( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyK/te_used);
01735 H2_CollRate[0][iVibHi][iRotHi][iVibLo][iRotLo] = (float)(
01736 excit*H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].gLo/
01737 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].gHi /
01738 SDIV(excit1) )*mole.lgColl_deexec_Calc *
01739
01740 mole.lgH2_ortho_para_coll_on;
01741
01742 if( PRT_COLL )
01743 fprintf(ioQQQ,"col o-p\t%li\t%li\t%li\t%li\t%li\t%.4e\t%.4e\n",
01744 nColl,
01745 iVibHi,iRotHi,iVibLo,iRotLo,
01746 energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo],
01747 H2_CollRate[nColl][iVibHi][iRotHi][iVibLo][iRotLo] );
01748
01749
01750 iRotLo = 1;
01751 iRotHi = 2;
01752 excit = sexp( -(POW2(5.35-454./te_used)-23.1 ) )*1e-13;
01753 excit1 = sexp( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyK/te_used);
01754 H2_CollRate[0][iVibHi][iRotHi][iVibLo][iRotLo] = (float)(
01755 excit*H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].gLo/
01756 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].gHi /
01757 SDIV(excit1) )*mole.lgColl_deexec_Calc *
01758
01759 mole.lgH2_ortho_para_coll_on;
01760
01761
01762
01763
01764
01765
01766
01767 if( phycon.te < 100. )
01768 {
01769
01770
01771 H2_CollRate[0][0][1][0][0] = (float)(H2_CollRate[0][0][1][0][0]*exp(-(3900-170.5)*(1./phycon.te - 1./100.)));
01772 H2_CollRate[0][0][3][0][0] = (float)(H2_CollRate[0][0][3][0][0]*exp(-(3900-1015.1)*(1./phycon.te - 1./100.)));
01773 H2_CollRate[0][0][2][0][1] = (float)(H2_CollRate[0][0][2][0][1]*exp(-(3900-339.3)*(1./phycon.te - 1./100.)));
01774 }
01775
01776 if( PRT_COLL )
01777 fprintf(ioQQQ,"col o-p\t%li\t%li\t%li\t%li\t%li\t%.4e\t%.4e\n",
01778 nColl,
01779 iVibHi,iRotHi,iVibLo,iRotLo,
01780 energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo],
01781 H2_CollRate[nColl][iVibHi][iRotHi][iVibLo][iRotLo] );
01782 }
01783
01784 if( mole.nH2_TRACE >= mole.nH2_trace_full )
01785 fprintf(ioQQQ,
01786 " collision rates updated for new temp, number of trans is %li\n",
01787 numb_coll_trans);
01788
01789 DEBUG_EXIT( "H2_collid_rates()" );
01790
01791 return;
01792 }
01793
01794
01795
01796 void H2_LevelPops( void )
01797 {
01798 static double TeUsedColl=-1.f;
01799 double H2_renorm_conserve=0.,
01800 H2_renorm_conserve_init=0. ,
01801 sumold,
01802 H2_BigH2_H2s,
01803 H2_BigH2_H2g;
01804 double old_solomon_rate=-1.;
01805 long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo;
01806 long int i;
01807 long int n_pop_oscil = 0;
01808 int kase=0;
01809 bool lgConv_h2_soln,
01810 lgPopsConv_total,
01811 lgPopsConv_relative,
01812 lgHeatConv,
01813 lgSolomonConv,
01814 lgOrthoParaRatioConv;
01815 double quant_old=-1.,
01816 quant_new=-1.;
01817
01818 bool lgH2_pops_oscil=false,
01819 lgH2_pops_ever_oscil=false;
01820 long int nEner,
01821 ipHi, ipLo;
01822 long int iElec , iVib , iRot,ip;
01823 double sum_pops_matrix;
01824 float collup;
01825
01826 static double ortho_para_old=0. , ortho_para_older=0. , ortho_para_current=0.;
01827 float frac_new_oscil=1.f;
01828
01829
01830 double PopChgMax_relative=0. , PopChgMaxOld_relative=0., PopChgMax_total=0., PopChgMaxOld_total=0.;
01831 long int iRotMaxChng_relative , iVibMaxChng_relative,
01832 iRotMaxChng_total , iVibMaxChng_total,
01833 nXLevelsMatrix_save;
01834 double popold_relative , popnew_relative , popold_total , popnew_total;
01835
01836 char chReason[100];
01837
01838 double flux_accum_photodissoc_BigH2_H2g, flux_accum_photodissoc_BigH2_H2s;
01839 long int ip_H2_level;
01840
01841
01842 double converge_pops_relative=0.1 ,
01843 converge_pops_total=1e-3,
01844 converge_ortho_para=1e-2;
01845
01846
01847
01848 if( !h2.lgH2ON || lgAbort )
01849 return;
01850
01851 DEBUG_ENTRY( "H2_LevelPops()" );
01852
01853 if(mole.nH2_TRACE >= mole.nH2_trace_full )
01854 {
01855 fprintf(ioQQQ,
01856 "\n***************H2_LevelPops call %li this iteration, zone is %.2f, H2/H:%.e Te:%e ne:%e\n",
01857 nCallH2_this_iteration,
01858 fnzone,
01859 hmi.H2_total/dense.gas_phase[ipHYDROGEN],
01860 phycon.te,
01861 dense.eden
01862 );
01863 }
01864 else if( mole.nH2_TRACE >= mole.nH2_trace_final )
01865 {
01866 static long int nzone_prt=-1;
01867 if( nzone!=nzone_prt )
01868 {
01869 nzone_prt = nzone;
01870 fprintf(ioQQQ,"DEBUG zone %li H2/H:%.3e Te:%.3e *ne:%.3e n(H2):%.3e\n",
01871 nzone,
01872 hmi.H2_total / dense.gas_phase[ipHYDROGEN],
01873 phycon.te,
01874 dense.eden,
01875 hmi.H2_total );
01876 }
01877 }
01878
01879
01880
01881 mole_H2_LTE();
01882
01883
01884
01885
01886 if( (!hmi.lgBigH2_evaluated && hmi.H2_total/dense.gas_phase[ipHYDROGEN] < mole.H2_to_H_limit )
01887 || hmi.H2_total < 1e-20 )
01888 {
01889
01890 if( mole.nH2_TRACE >= mole.nH2_trace_full )
01891 fprintf(ioQQQ,
01892 " H2_LevelPops pops too small, not computing, set to lte and return, H2/H is %.2e and mole.H2_to_H_limit is %.2e.",
01893 hmi.H2_total/dense.gas_phase[ipHYDROGEN] ,
01894 mole.H2_to_H_limit);
01895 H2_zero_pops_too_low();
01896
01897 DEBUG_EXIT( "H2_LevelPops()" );
01898 return;
01899 }
01900
01901
01902
01903
01904
01905
01906
01907 hmi.lgBigH2_evaluated = true;
01908
01909
01910
01911
01912
01913 nXLevelsMatrix_save = nXLevelsMatrix;
01914 if( conv.lgSearch )
01915 {
01916 nXLevelsMatrix = nLevels_per_elec[0];
01917 }
01918
01919
01920
01921
01922
01923
01924
01925
01926 if( fabs(1. - TeUsedColl / phycon.te ) > 0.005 )
01927 {
01928 H2_collid_rates();
01929 TeUsedColl = phycon.te;
01930 }
01931
01932
01933 if( PRT_COLL )
01934 exit(98);
01935
01936
01937
01938
01939 if( nCallH2_this_iteration==0 || mole.lgH2_LTE )
01940 {
01941
01942 if(mole.nH2_TRACE >= mole.nH2_trace_full )
01943 fprintf(ioQQQ,"H2 1st call - using lte level pops\n");
01944
01945 for( iElec=0; iElec<mole.n_h2_elec_states; ++iElec )
01946 {
01947 pops_per_elec[iElec] = 0.;
01948 for( iVib=0; iVib<=h2.nVib_hi[iElec]; ++iVib )
01949 {
01950 pops_per_vib[iElec][iVib] = 0.;
01951 for( iRot=h2.Jlowest[iElec]; iRot<=h2.nRot_hi[iElec][iVib]; ++iRot )
01952 {
01953
01954
01955 H2_old_populations[iElec][iVib][iRot] =
01956 (float)H2_populations_LTE[iElec][iVib][iRot]*hmi.H2_total;
01957 H2_populations[iElec][iVib][iRot] = H2_old_populations[iElec][iVib][iRot];
01958 }
01959 }
01960 }
01961
01962 h2.ortho_density = 0.75*hmi.H2_total;
01963 h2.para_density = 0.25*hmi.H2_total;
01964 ortho_para_current = h2.ortho_density / SDIV( h2.para_density );
01965 ortho_para_older = ortho_para_current;
01966 ortho_para_old = ortho_para_current;
01967
01968 frac_matrix = 1.;
01969 }
01970
01971
01972 iElec = 0;
01973 pops_per_elec[0] = 0.;
01974 for( iVib=0; iVib<=h2.nVib_hi[iElec]; ++iVib )
01975 {
01976 pops_per_vib[0][iVib] = 0.;
01977 for( iRot=h2.Jlowest[iElec]; iRot<=h2.nRot_hi[iElec][iVib]; ++iRot )
01978 {
01979 pops_per_elec[0] += H2_populations[iElec][iVib][iRot];
01980 pops_per_vib[0][iVib] += H2_populations[iElec][iVib][iRot];
01981 }
01982 }
01983 ASSERT( pops_per_elec[0]>SMALLFLOAT );
01984
01985
01986
01987
01988 H2_renorm_chemistry = hmi.H2_total/ SDIV(pops_per_elec[0]);
01989
01990
01991
01992 iElec = 0;
01993 hmi.H2_chem_BigH2_H2g = 0.;
01994 hmi.H2_chem_BigH2_H2s = 0.;
01995 for( iVib=0; iVib<=h2.nVib_hi[iElec]; ++iVib )
01996 {
01997 for( iRot=h2.Jlowest[iElec]; iRot<=h2.nRot_hi[iElec][iVib]; ++iRot )
01998 {
01999 if( energy_wn[0][iVib][iRot] > ENERGY_H2_STAR )
02000 {
02001 hmi.H2_chem_BigH2_H2s += H2_populations[iElec][iVib][iRot];
02002
02003 }
02004 else
02005 {
02006 hmi.H2_chem_BigH2_H2g += H2_populations[iElec][iVib][iRot];
02007
02008 }
02009 }
02010 }
02011
02012 hmi.H2_chem_BigH2_H2g = hmi.Hmolec[ipMH2g]/SDIV(hmi.H2_chem_BigH2_H2g);
02013 hmi.H2_chem_BigH2_H2s = hmi.Hmolec[ipMH2s]/SDIV(hmi.H2_chem_BigH2_H2s);
02014
02015
02016 if(mole.nH2_TRACE >= mole.nH2_trace_full)
02017 fprintf(ioQQQ,
02018 "H2 H2_renorm_chemistry is %.4e, hmi.H2_total is %.4e pops_per_elec[0] is %.4e\n",
02019 H2_renorm_chemistry ,
02020 hmi.H2_total,
02021 pops_per_elec[0]);
02022
02023
02024 iElec = 0;
02025 for( iVib=0; iVib<=h2.nVib_hi[iElec]; ++iVib )
02026 {
02027 for( iRot=h2.Jlowest[iElec]; iRot<=h2.nRot_hi[iElec][iVib]; ++iRot )
02028 {
02029 H2_populations[iElec][iVib][iRot] *= H2_renorm_chemistry;
02030 H2_old_populations[iElec][iVib][iRot] = H2_populations[iElec][iVib][iRot];
02031 }
02032 }
02033 ASSERT( fabs(h2.ortho_density+h2.para_density-hmi.H2_total)/hmi.H2_total < 0.001 );
02034
02035 if(mole.nH2_TRACE >= mole.nH2_trace_full )
02036 fprintf(ioQQQ,
02037 " H2 entry, old pops sumed to %.3e, renorm to htwo den of %.3e\n",
02038 pops_per_elec[0],
02039 hmi.H2_total);
02040
02041
02042 if( conv.lgSearch )
02043 {
02044 converge_pops_relative *= 2.;
02045 converge_pops_total *= 3.;
02046 converge_ortho_para *= 3.;
02047 }
02048
02049
02050 mole_H2_form();
02051
02052
02053 H2_X_coll_rate_evaluate();
02054
02055
02056
02057 lgConv_h2_soln = false;
02058
02059 loop_h2_pops = 0;
02060 {
02061 static long int nzoneEval=-1;
02062 if( nzone != nzoneEval )
02063 {
02064 nzoneEval = nzone;
02065
02066 ++nH2_zone;
02067 }
02068 }
02069
02070
02071
02072
02073
02074
02075
02076 while( loop_h2_pops < LIM_H2_POP_LOOP-n_pop_oscil && !lgConv_h2_soln && !mole.lgH2_LTE )
02077 {
02078 double rate_in , rate_out;
02079 static double old_HeatH2Dexc_BigH2=0., HeatChangeOld=0. , HeatChange=0.;
02080
02081
02082 ++loop_h2_pops;
02083
02084 ++nH2_pops;
02085
02086
02087
02088
02089 hmi.H2_H2g_to_H2s_rate_BigH2 = 0.;
02090
02091
02092
02093
02094
02095 rate_out = 0.;
02096
02097
02098
02099 iElec = 0;
02100 for( iVib=0; iVib<=h2.nVib_hi[iElec]; ++iVib )
02101 {
02102 for( iRot=h2.Jlowest[iElec]; iRot<=h2.nRot_hi[iElec][iVib]; ++iRot )
02103 {
02104
02105 H2_X_rate_from_elec_excited[iVib][iRot] = 0.;
02106
02107 H2_X_rate_to_elec_excited[iVib][iRot] = 0.;
02108 }
02109 }
02110
02111 iElecHi = -INT32_MAX;
02112
02113
02114 for( iElecHi=1; iElecHi<mole.n_h2_elec_states; ++iElecHi )
02115 {
02116
02117
02118 iVib = -INT32_MAX;
02119 iRot = -INT32_MAX;
02120 iVibLo = -INT32_MAX;
02121 iRotLo = -INT32_MAX;
02122 iVibHi = -INT32_MAX;
02123 iRotHi = -INT32_MAX;
02124
02125
02126 pops_per_elec[iElecHi] = 0.;
02127
02128 if(mole.nH2_TRACE >= mole.nH2_trace_full)
02129 fprintf(ioQQQ," Pop(e=%li):",iElecHi);
02130
02131
02132 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
02133 {
02134 pops_per_vib[iElecHi][iVibHi] = 0.;
02135
02136 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
02137 {
02138
02139
02140
02141 rate_in = 0.;
02142
02143
02144
02145 rate_out = H2_dissprob[iElecHi][iVibHi][iRotHi];
02146
02147
02148 iElecLo=0;
02149 for( iVibLo=0; iVibLo<=h2.nVib_hi[iElecLo]; ++iVibLo )
02150 {
02151 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=h2.nRot_hi[iElecLo][iVibLo]; ++iRotLo )
02152 {
02153
02154 if( lgH2_line_exists[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] )
02155 {
02156
02157
02158
02159 double rate_one =
02160 (H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].pump
02161
02162
02163
02164
02165
02166 +hmi.lgLeidenCRHack*secondaries.x12tot/2. );
02169
02170 ASSERT( H2_lgOrtho[iElecHi][iVibHi][iRotHi]-H2_lgOrtho[iElecLo][iVibLo][iRotLo]==0 );
02171
02172
02173 rate_in += H2_old_populations[iElecLo][iVibLo][iRotLo]*rate_one;
02174
02175
02176
02177
02178
02179
02180 H2_X_rate_to_elec_excited[iVibLo][iRotLo] += rate_one;
02181
02182
02183
02184
02185
02186
02187 rate_one = H2Lines[iElecHi][iVibHi][iRotHi][0][iVibLo][iRotLo].Aul*
02188
02189 (H2Lines[iElecHi][iVibHi][iRotHi][0][iVibLo][iRotLo].Pesc +
02190 H2Lines[iElecHi][iVibHi][iRotHi][0][iVibLo][iRotLo].Pelec_esc +
02191 H2Lines[iElecHi][iVibHi][iRotHi][0][iVibLo][iRotLo].Pdest) +
02192
02193 H2Lines[iElecHi][iVibHi][iRotHi][0][iVibLo][iRotLo].pump *
02194 H2Lines[iElecHi][iVibHi][iRotHi][0][iVibLo][iRotLo].gLo/
02195 H2Lines[iElecHi][iVibHi][iRotHi][0][iVibLo][iRotLo].gHi;
02196
02197 # if 0
02198
02199 if( energy_wn[0][iVibLo][iRotLo] > ENERGY_H2_STAR )
02200 {
02201
02202
02203 H2s_pump += H2_old_populations[0][iVibLo][iRotLo]*(H2Lines[iElecHi][iVibHi][iRotHi][0][iVibLo][iRotLo].pump+hmi.lgLeidenCRHack*secondaries.x12tot/2.f);
02204
02205 ratio_H2s_pump = H2s_pump/SDIV(rate_tot);
02206
02207 hmi.H2_H2g_to_H2s_rate_BigH2 += (H2Lines[iElecHi][iVibHi][iRotHi][0][iVibLo][iRotLo].Aul*
02208
02209 (H2Lines[iElecHi][iVibHi][iRotHi][0][iVibLo][iRotLo].Pesc +
02210 H2Lines[iElecHi][iVibHi][iRotHi][0][iVibLo][iRotLo].Pelec_esc +
02211 H2Lines[iElecHi][iVibHi][iRotHi][0][iVibLo][iRotLo].Pdest) +
02212
02213 H2Lines[iElecHi][iVibHi][iRotHi][0][iVibLo][iRotLo].pump *
02214 H2Lines[iElecHi][iVibHi][iRotHi][0][iVibLo][iRotLo].gLo/
02215 H2Lines[iElecHi][iVibHi][iRotHi][0][iVibLo][iRotLo].gHi)*H2_old_populations[iElecHi][iVibHi][iRotHi]*(1.-ratio_H2s_pump)/SDIV(hmi.H2_total);
02216
02217
02218
02219
02220 hmi.H2_H2s_to_H2s_rate_BigH2 += (H2Lines[iElecHi][iVibHi][iRotHi][0][iVibLo][iRotLo].Aul*
02221
02222 (H2Lines[iElecHi][iVibHi][iRotHi][0][iVibLo][iRotLo].Pesc +
02223 H2Lines[iElecHi][iVibHi][iRotHi][0][iVibLo][iRotLo].Pelec_esc +
02224 H2Lines[iElecHi][iVibHi][iRotHi][0][iVibLo][iRotLo].Pdest) +
02225
02226 H2Lines[iElecHi][iVibHi][iRotHi][0][iVibLo][iRotLo].pump *
02227 H2Lines[iElecHi][iVibHi][iRotHi][0][iVibLo][iRotLo].gLo/
02228 H2Lines[iElecHi][iVibHi][iRotHi][0][iVibLo][iRotLo].gHi)*H2_old_populations[iElecHi][iVibHi][iRotHi]*ratio_H2s_pump/SDIV(hmi.H2_total);
02229
02230 }
02231 # endif
02232
02233
02234
02235 rate_out += rate_one;
02236
02237 ASSERT( (rate_in >=0.) && rate_out >= 0. );
02238 }
02239 }
02240 }
02241
02242
02243
02244
02245
02246 H2_rad_rate_out[iElecHi][iVibHi][iRotHi] = rate_out;
02247 H2_populations[iElecHi][iVibHi][iRotHi] = rate_in / SDIV( rate_out);
02248 if( H2_old_populations[iElecHi][iVibHi][iRotHi]==0. )
02249 H2_old_populations[iElecHi][iVibHi][iRotHi] = H2_populations[iElecHi][iVibHi][iRotHi];
02250
02251
02252
02253 iElecLo=0;
02254 for( iVibLo=0; iVibLo<=h2.nVib_hi[iElecLo]; ++iVibLo )
02255 {
02256 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=h2.nRot_hi[iElecLo][iVibLo]; ++iRotLo )
02257 {
02258 if( lgH2_line_exists[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] )
02259 {
02260 double rate_one =
02261
02262
02263 H2_populations[iElecHi][iVibHi][iRotHi]*
02264 (H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Aul*
02265
02266 (H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Pesc +
02267 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Pelec_esc +
02268 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Pdest) +
02269
02270 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].pump *
02271 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].gLo/
02272 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].gHi);
02273
02274
02275 H2_X_rate_from_elec_excited[iVibLo][iRotLo] += rate_one;
02276 }
02277 }
02278 }
02279
02280 ASSERT( H2_populations[iElecHi][iVibHi][iRotHi] >= 0. &&
02281 H2_populations[iElecHi][iVibHi][iRotHi] <= hmi.H2_total );
02282
02283
02284 pops_per_vib[iElecHi][iVibHi] += H2_populations[iElecHi][iVibHi][iRotHi];
02285
02286
02287 }
02288
02289 if(mole.nH2_TRACE >= mole.nH2_trace_full)
02290 fprintf(ioQQQ,"\t%.2e",pops_per_vib[iElecHi][iVibHi]/hmi.H2_total);
02291
02292
02293 pops_per_elec[iElecHi] += pops_per_vib[iElecHi][iVibHi];
02294 }
02295
02296
02297 if(mole.nH2_TRACE >= mole.nH2_trace_full)
02298 fprintf(ioQQQ,"\n");
02299 }
02300
02301
02302
02303
02304
02305
02306 PopChgMaxOld_relative = PopChgMax_relative;
02307 PopChgMaxOld_total = PopChgMax_total;
02308 PopChgMax_relative = 0.;
02309 PopChgMax_total = 0.;
02310 iElec = 0;
02311 iElecHi = 0;
02312 iRotMaxChng_relative =-1;
02313 iVibMaxChng_relative = -1;
02314 iRotMaxChng_total =-1;
02315 iVibMaxChng_total = -1;
02316 popold_relative = 0.;
02317 popnew_relative = 0.;
02318 popold_total = 0.;
02319 popnew_total = 0.;
02320
02321
02322 for( nEner=0; nEner<nLevels_per_elec[0]; ++nEner )
02323 {
02324
02325 ip = H2_ipX_ener_sort[nEner];
02326 iVib = ipVib_H2_energy_sort[ip];
02327 iRot = ipRot_H2_energy_sort[ip];
02328
02329
02330 H2_col_rate_in[iVib][iRot] = 0.;
02331 H2_col_rate_out[iVib][iRot] = 0.;
02332 for( ipLo=0; ipLo<nEner; ++ipLo )
02333 {
02334 ip = H2_ipX_ener_sort[ipLo];
02335 iVibLo = ipVib_H2_energy_sort[ip];
02336 iRotLo = ipRot_H2_energy_sort[ip];
02337
02338
02339 H2_col_rate_out[iVib][iRot] += H2_X_coll_rate[nEner][ipLo];
02340
02341
02342
02343 collup = (float)(H2_old_populations[0][iVibLo][iRotLo] * H2_X_coll_rate[nEner][ipLo] *
02344 H2_stat[0][iVib][iRot] / H2_stat[0][iVibLo][iRotLo] *
02345 H2_Boltzmann[0][iVib][iRot] /
02346 SDIV( H2_Boltzmann[0][iVibLo][iRotLo] ) );
02347
02348 H2_col_rate_in[iVib][iRot] += collup;
02349 }
02350
02351 for( ipHi=nEner+1; ipHi<nLevels_per_elec[0]; ++ipHi )
02352 {
02353 double colldn;
02354 ip = H2_ipX_ener_sort[ipHi];
02355 iVibHi = ipVib_H2_energy_sort[ip];
02356 iRotHi = ipRot_H2_energy_sort[ip];
02357
02358
02359 H2_col_rate_out[iVib][iRot] += H2_X_coll_rate[ipHi][nEner] *
02360 H2_stat[0][iVibHi][iRotHi] / H2_stat[0][iVib][iRot] *
02361 (float)(H2_Boltzmann[0][iVibHi][iRotHi] /
02362 SDIV( H2_Boltzmann[0][iVib][iRot] ) );
02363
02364
02365 colldn = H2_old_populations[0][iVibHi][iRotHi] * H2_X_coll_rate[ipHi][nEner];
02366
02367
02368
02369
02370
02371 H2_col_rate_in[iVib][iRot] += colldn;
02372 }
02373 }
02374
02375
02376
02377
02378
02379
02380
02381
02382 nEner = nLevels_per_elec[0];
02383 while( (--nEner) >= nXLevelsMatrix )
02384 {
02385
02386
02387
02388 ip = H2_ipX_ener_sort[nEner];
02389 iVib = ipVib_H2_energy_sort[ip];
02390 iRot = ipRot_H2_energy_sort[ip];
02391
02392 if( nEner+1 < nLevels_per_elec[0] )
02393 ASSERT( H2_Xenergies[H2_ipX_ener_sort[nEner]] < H2_Xenergies[H2_ipX_ener_sort[nEner+1]] );
02394
02395
02396
02397
02398 if( nEner >1 )
02399 {
02400 H2_col_rate_out[iVib][iRot] +=
02401
02402
02403 (float)(hmi.rate_grain_h2_op_conserve);
02404
02405
02406
02407 H2_col_rate_in[0][H2_lgOrtho[0][iVib][iRot]] +=
02408
02409
02410
02411 (float)(hmi.rate_grain_h2_op_conserve*H2_old_populations[0][iVib][iRot]);
02412 }
02413 else if( nEner == 1 )
02414 {
02415
02416
02417 H2_col_rate_out[0][1] +=
02418
02419
02420
02421 (float)(hmi.rate_grain_h2_J1_to_J0);
02422
02423 H2_col_rate_in[0][0] +=
02424
02425
02426
02427 (float)(hmi.rate_grain_h2_J1_to_J0 *H2_old_populations[0][0][1]);
02428 }
02429
02430
02431 H2_rad_rate_in[iVib][iRot] = 0.;
02432 H2_rad_rate_out[0][iVib][iRot] = 0.;
02433
02434
02435
02436
02437
02438 H2_rad_rate_in[iVib][iRot] += H2_X_rate_from_elec_excited[iVib][iRot];
02439
02440
02441 H2_rad_rate_out[0][iVib][iRot] += H2_X_rate_to_elec_excited[iVib][iRot];
02442
02443
02444 iElecHi = 0;
02445 for( ipHi = nEner+1; ipHi<nLevels_per_elec[0]; ++ipHi )
02446 {
02447 ip = H2_ipX_ener_sort[ipHi];
02448 iVibHi = ipVib_H2_energy_sort[ip];
02449 iRotHi = ipRot_H2_energy_sort[ip];
02450
02451
02452
02453
02454
02455 if( (abs((iRotHi-iRot) )==2 || iRotHi==iRot )&& (iVib<=iVibHi) )
02456 {
02457 double rateone;
02458 rateone =
02459 H2_old_populations[iElecHi][iVibHi][iRotHi]*
02460 H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Aul*
02461 (H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Pesc +
02462 H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Pelec_esc +
02463 H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Pdest);
02464 ASSERT( rateone >=0 );
02465
02466
02467 H2_rad_rate_in[iVib][iRot] += rateone;
02468 }
02469 }
02470
02471
02472
02473
02474 iElecLo = 0;
02475 for( ipLo = 0; ipLo<nEner; ++ipLo )
02476 {
02477 ip = H2_ipX_ener_sort[ipLo];
02478 iVibLo = ipVib_H2_energy_sort[ip];
02479 iRotLo = ipRot_H2_energy_sort[ip];
02480
02481
02482
02483 if( ((abs((iRotLo-iRot)) == 2)||(abs((iRotLo-iRot)) == 0)) &&
02484 (iVibLo<=iVib) )
02485 {
02486 H2_rad_rate_out[0][iVib][iRot] +=
02487 H2Lines[iElec][iVib][iRot][iElecLo][iVibLo][iRotLo].Aul*
02488 (H2Lines[iElec][iVib][iRot][iElecLo][iVibLo][iRotLo].Pesc +
02489 H2Lines[iElec][iVib][iRot][iElecLo][iVibLo][iRotLo].Pelec_esc +
02490 H2Lines[iElec][iVib][iRot][iElecLo][iVibLo][iRotLo].Pdest);
02491 }
02492 }
02493
02494
02495
02496
02497 H2_populations[iElec][iVib][iRot] =
02498 (H2_col_rate_in[iVib][iRot]+ H2_rad_rate_in[iVib][iRot]+H2_X_source[nEner]) /
02499 SDIV(H2_col_rate_out[iVib][iRot]+H2_rad_rate_out[0][iVib][iRot]+H2_X_sink[nEner]);
02500
02501 ASSERT( H2_populations[iElec][iVib][iRot] >= 0. );
02502 }
02503
02504
02505
02506
02507 if( nXLevelsMatrix )
02508 {
02509 H2_Level_low_matrix(
02510
02511
02512 hmi.H2_total * (float)frac_matrix );
02513
02514
02515
02516
02517
02518
02519
02520 }
02521 iElecHi = 0;
02522 if(mole.nH2_TRACE >= mole.nH2_trace_full)
02523 {
02524 fprintf(ioQQQ," Rel pop(e=%li)" ,iElecHi);
02525 }
02526
02527
02528
02529 pops_per_elec[0] = 0.;
02530 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
02531 {
02532 double sumv;
02533 sumv = 0.;
02534 pops_per_vib[0][iVibHi] = 0.;
02535
02536 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
02537 {
02538 pops_per_elec[0] += H2_populations[iElecHi][iVibHi][iRotHi];
02539 sumv += H2_populations[iElecHi][iVibHi][iRotHi];
02540 pops_per_vib[0][iVibHi] += H2_populations[iElecHi][iVibHi][iRotHi];
02541 }
02542
02543 if(mole.nH2_TRACE >= mole.nH2_trace_full)
02544 fprintf(ioQQQ,"\t%.2e",sumv/hmi.H2_total);
02545 }
02546 ASSERT( pops_per_elec[0] > SMALLFLOAT );
02547
02548 if(mole.nH2_TRACE >= mole.nH2_trace_full)
02549 {
02550 fprintf(ioQQQ,"\n");
02551
02552 fprintf(ioQQQ," Rel pop(0,J)");
02553 iElecHi = 0;
02554 iVibHi = 0;
02555 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
02556 {
02557 fprintf(ioQQQ,"\t%.2e",H2_populations[iElecHi][iVibHi][iRotHi]/hmi.H2_total);
02558 }
02559 fprintf(ioQQQ,"\n");
02560 }
02561
02562
02563
02564 iElec = 0;
02565 sum_pops_matrix = 0.;
02566 ip =0;
02567 for( i=0; i<nXLevelsMatrix; ++i )
02568 {
02569 ip = H2_ipX_ener_sort[i];
02570 iVib = ipVib_H2_energy_sort[ip];
02571 iRot = ipRot_H2_energy_sort[ip];
02572 sum_pops_matrix += H2_populations[iElec][iVib][iRot];
02573 }
02574
02575
02576
02577
02578 frac_matrix = sum_pops_matrix / SDIV(pops_per_elec[0]);
02579
02580
02581
02582
02583
02584 H2_renorm_conserve = hmi.H2_total/ SDIV(pops_per_elec[0]);
02585
02586
02587
02588
02589
02590 H2_sum_excit_elec_den = 0.;
02591 for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi )
02592 {
02593 pops_per_elec[iElecHi] *= H2_renorm_conserve;
02594 if( iElecHi > 0 )
02595 H2_sum_excit_elec_den += pops_per_elec[iElecHi];
02596
02597 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
02598 {
02599 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
02600 {
02601 H2_populations[iElecHi][iVibHi][iRotHi] *= H2_renorm_conserve;
02602
02603 }
02604 }
02605 }
02606
02607
02608
02609
02610
02611 sumold = 0.;
02612 for( iElec=0; iElec<mole.n_h2_elec_states; ++iElec )
02613 {
02614 for( iVib=0; iVib<=h2.nVib_hi[iElec]; ++iVib )
02615 {
02616 for( iRot=h2.Jlowest[iElec]; iRot<=h2.nRot_hi[iElec][iVib]; ++iRot )
02617 {
02618 double rel_change;
02619
02620
02621 if( fabs(H2_populations[iElec][iVib][iRot] -
02622 H2_old_populations[iElec][iVib][iRot])/
02623
02624
02625
02626 SDIV(H2_populations[iElec][iVib][iRot]) > fabs(PopChgMax_relative) &&
02627
02628
02629
02630
02631
02632
02633 H2_populations[iElec][iVib][iRot]/SDIV(hmi.H2_total)>1e-6 )
02634 {
02635 PopChgMax_relative =
02636 (H2_populations[iElec][iVib][iRot] -
02637 H2_old_populations[iElec][iVib][iRot])/
02638 SDIV(H2_populations[iElec][iVib][iRot]);
02639 iRotMaxChng_relative = iRot;
02640 iVibMaxChng_relative = iVib;
02641 popold_relative = H2_old_populations[iElec][iVib][iRot];
02642 popnew_relative = H2_populations[iElec][iVib][iRot];
02643 }
02644
02645
02646
02647
02648 if( fabs(H2_populations[iElec][iVib][iRot] -
02649 H2_old_populations[iElec][iVib][iRot])/
02650
02651 SDIV(hmi.H2_total) > fabs(PopChgMax_total) )
02652 {
02653 PopChgMax_total =
02654 (H2_populations[iElec][iVib][iRot] -
02655 H2_old_populations[iElec][iVib][iRot])/
02656 SDIV(H2_populations[iElec][iVib][iRot]);
02657 iRotMaxChng_total = iRot;
02658 iVibMaxChng_total = iVib;
02659 popold_total = H2_old_populations[iElec][iVib][iRot];
02660 popnew_total = H2_populations[iElec][iVib][iRot];
02661 }
02662
02663 kase = -1;
02664
02665
02666
02667
02668
02669
02670 rel_change = fabs( H2_old_populations[iElec][iVib][iRot] -
02671 H2_populations[iElec][iVib][iRot] )/
02672 SDIV( H2_populations[iElec][iVib][iRot] );
02673
02674
02675 if( rel_change > 3. &&
02676 H2_old_populations[iElec][iVib][iRot]*H2_populations[iElec][iVib][iRot]>0 )
02677 {
02678
02679 H2_old_populations[iElec][iVib][iRot] = pow( 10. ,
02680 log10(H2_old_populations[iElec][iVib][iRot])/2. +
02681 log10(H2_populations[iElec][iVib][iRot])/2. );
02682 kase = 2;
02683 }
02684
02685
02686 else if( rel_change> 0.1 )
02687 {
02688 float frac_old=0.25f;
02689
02690 H2_old_populations[iElec][iVib][iRot] =
02691 frac_old*H2_old_populations[iElec][iVib][iRot] +
02692 (1.f-frac_old)*H2_populations[iElec][iVib][iRot];
02693 kase = 3;
02694 }
02695 else
02696 {
02697
02698 H2_old_populations[iElec][iVib][iRot] =
02699 H2_populations[iElec][iVib][iRot];
02700 kase = 4;
02701 }
02702 sumold += H2_old_populations[iElec][iVib][iRot];
02703 }
02704 }
02705 }
02706
02707 H2_renorm_conserve_init = hmi.H2_total/sumold;
02708
02709
02710
02711
02712 for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi )
02713 {
02714 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
02715 {
02716 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
02717 {
02718 H2_old_populations[iElecHi][iVibHi][iRotHi] *= H2_renorm_conserve_init;
02719
02720 }
02721 }
02722 }
02723
02724 iElecHi = 0;
02725 h2.ortho_density = 0.;
02726 h2.para_density = 0.;
02727 H2_den_s = 0.;
02728 H2_den_g = 0.;;
02729
02730 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
02731 {
02732 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
02733 {
02734
02735 if( energy_wn[0][iVibHi][iRotHi]> ENERGY_H2_STAR )
02736 {
02737 H2_den_s += H2_populations[iElecHi][iVibHi][iRotHi];
02738 }
02739 else
02740 {
02741 H2_den_g += H2_populations[iElecHi][iVibHi][iRotHi];
02742 }
02743 if( H2_lgOrtho[iElecHi][iVibHi][iRotHi] )
02744 {
02745 h2.ortho_density += H2_populations[iElecHi][iVibHi][iRotHi];
02746 }
02747 else
02748 {
02749 h2.para_density += H2_populations[iElecHi][iVibHi][iRotHi];
02750 }
02751
02752 }
02753 }
02754
02755 ortho_para_older = ortho_para_old;
02756 ortho_para_old = ortho_para_current;
02757 ortho_para_current = h2.ortho_density / SDIV( h2.para_density );
02758
02759
02760
02761 old_solomon_rate = hmi.H2_Solomon_dissoc_rate_BigH2_H2g;
02762
02763
02764
02765 H2_Solomon_rate( );
02766
02767
02768
02769
02770 HeatChangeOld = HeatChange;
02771 HeatChange = old_HeatH2Dexc_BigH2 - hmi.HeatH2Dexc_BigH2;
02772 {
02773 static long int loop_h2_oscil=-1;
02774
02775
02776 if( loop_h2_pops>2 && (
02777 (HeatChangeOld*HeatChange<0. ) ||
02778 (PopChgMax_relative*PopChgMaxOld_relative<0. ) ) )
02779 {
02780 lgH2_pops_oscil = true;
02781 if( loop_h2_pops > 6 )
02782 {
02783 loop_h2_oscil = loop_h2_pops;
02784 lgH2_pops_ever_oscil = true;
02785
02786
02787 frac_new_oscil *= 0.8f;
02788 frac_new_oscil = MAX2( frac_new_oscil , 0.1f);
02789 ++n_pop_oscil;
02790 }
02791 }
02792 else
02793 {
02794 lgH2_pops_oscil = false;
02795
02796 if( loop_h2_pops - loop_h2_oscil > 4 )
02797 {
02798 frac_new_oscil = 1.f;
02799 lgH2_pops_ever_oscil = false;
02800 }
02801 }
02802 }
02803
02804
02805
02806 old_HeatH2Dexc_BigH2 = hmi.HeatH2Dexc_BigH2;
02807 if(fabs(hmi.HeatH2Dexc_BigH2)/thermal.ctot > conv.HeatCoolRelErrorAllowed/10. ||
02808 hmi.HeatH2Dexc_BigH2==0. )
02809 H2_Cooling("H2lup");
02810
02811
02812 lgConv_h2_soln = true;
02813 lgPopsConv_total = true;;
02814 lgPopsConv_relative = true;;
02815 lgHeatConv = true;
02816 lgSolomonConv = true;
02817 lgOrthoParaRatioConv = true;
02818
02819
02820
02821 if( fabs(PopChgMax_relative)>converge_pops_relative )
02822 {
02823
02824 lgConv_h2_soln = false;
02825 lgPopsConv_relative = false;
02826
02827
02828 quant_old = PopChgMaxOld_relative;
02829
02830 quant_new = PopChgMax_relative;
02831
02832 strcpy( chReason , "rel pops changed" );
02833 }
02834
02835
02836
02837 else if( fabs(PopChgMax_total)>converge_pops_total)
02838 {
02839 lgConv_h2_soln = false;
02840 lgPopsConv_total = false;
02841
02842
02843 quant_old = PopChgMaxOld_total;
02844
02845 quant_new = PopChgMax_total;
02846
02847 strcpy( chReason , "tot pops changed" );
02848 }
02849
02850
02851
02852
02853
02854 else if( fabs(ortho_para_current-ortho_para_old) / SDIV(ortho_para_current)> converge_ortho_para )
02855
02856
02857 {
02858 lgConv_h2_soln = false;
02859 lgOrthoParaRatioConv = false;
02860 quant_old = ortho_para_old;
02861 quant_new = ortho_para_current;
02862 strcpy( chReason , "ortho/para ratio changed" );
02863 }
02864
02865
02866 else if( !thermal.lgTSetOn &&
02867 fabs(hmi.HeatH2Dexc_BigH2-old_HeatH2Dexc_BigH2)/MAX2(thermal.ctot,thermal.htot) >
02868 conv.HeatCoolRelErrorAllowed/2.
02869
02870 )
02871 {
02872
02873
02874
02875 lgConv_h2_soln = false;
02876 lgHeatConv = false;
02877 quant_old = old_HeatH2Dexc_BigH2/MAX2(thermal.ctot,thermal.htot);
02878 quant_new = hmi.HeatH2Dexc_BigH2/MAX2(thermal.ctot,thermal.htot);
02879 strcpy( chReason , "heating changed" );
02880
02881
02882
02883 }
02884
02885
02886
02887
02888
02889 else if( rfield.lgInducProcess &&
02890
02891
02892 hmi.H2_frac_abund_set==0 &&
02893
02894
02895 fabs( hmi.H2_Solomon_dissoc_rate_BigH2_H2g - old_solomon_rate)/SDIV(hmi.H2_rate_destroy) >
02896 conv.EdenErrorAllowed/5.)
02897 {
02898 lgConv_h2_soln = false;
02899 lgSolomonConv = false;
02900 quant_old = old_solomon_rate;
02901 quant_new = hmi.H2_Solomon_dissoc_rate_BigH2_H2g;
02902 strcpy( chReason , "Solomon rate changed" );
02903 }
02904
02905
02906 if( !lgConv_h2_soln )
02907
02908
02909
02910
02911
02912
02913
02914
02915
02918 {
02919
02920
02921
02922 if( PRT_POPS || mole.nH2_TRACE >=mole.nH2_trace_iterations )
02923 {
02924
02925
02926
02927
02928 fprintf(ioQQQ," loop %3li no conv oscl?%c why:%s ",
02929 loop_h2_pops,
02930 TorF(lgH2_pops_ever_oscil),
02931 chReason );
02932 if( !lgPopsConv_relative )
02933 fprintf(ioQQQ," PopChgMax_relative:%.4e v:%li J:%li old:%.4e new:%.4e",
02934 PopChgMax_relative,
02935 iVibMaxChng_relative,
02936 iRotMaxChng_relative ,
02937 popold_relative ,
02938 popnew_relative );
02939 else if( !lgPopsConv_total )
02940 fprintf(ioQQQ," PopChgMax_total:%.4e v:%li J:%li old:%.4e new:%.4e",
02941 PopChgMax_total,
02942 iVibMaxChng_total,
02943 iRotMaxChng_total ,
02944 popold_total ,
02945 popnew_total );
02946 else if( !lgHeatConv )
02947 fprintf(ioQQQ," heat:%.4e old:%.4e new:%.4e",
02948 (hmi.HeatH2Dexc_BigH2-old_HeatH2Dexc_BigH2)/MAX2(thermal.ctot,thermal.htot),
02949 quant_old ,
02950 quant_new);
02951
02952 else if( !lgSolomonConv )
02953 fprintf(ioQQQ," d(sol rate)/tot dest\t%2e",(old_solomon_rate - hmi.H2_Solomon_dissoc_rate_BigH2_H2g)/SDIV(hmi.H2_rate_destroy));
02954 else if( !lgOrthoParaRatioConv )
02955 fprintf(ioQQQ," current, old, older ratios are %.4e %.4e %.4e",
02956 ortho_para_current , ortho_para_old, ortho_para_older );
02957 else
02958 TotalInsanity();
02959 fprintf(ioQQQ,"\n");
02960 }
02961 }
02962
02963
02964
02965
02966
02967
02968
02969
02970
02971
02972
02973
02974
02975 if( trace.nTrConvg >= 5 )
02976 {
02977 fprintf( ioQQQ,
02978 " H2 5lev %li Conv?%c",
02979 loop_h2_pops ,
02980 TorF(lgConv_h2_soln) );
02981
02982 if( fabs(PopChgMax_relative)>0.1 )
02983 fprintf(ioQQQ," pops, rel chng %.3e",PopChgMax_relative);
02984 else
02985 fprintf(ioQQQ," rel heat %.3e rel chng %.3e H2 heat/cool %.2e",
02986 hmi.HeatH2Dexc_BigH2/thermal.ctot ,
02987 fabs(hmi.HeatH2Dexc_BigH2-old_HeatH2Dexc_BigH2)/thermal.ctot ,
02988 hmi.HeatH2Dexc_BigH2/thermal.ctot);
02989
02990 fprintf( ioQQQ,
02991 " Oscil?%c Ever Oscil?%c",
02992 TorF(lgH2_pops_oscil) ,
02993 TorF(lgH2_pops_ever_oscil) );
02994 if( lgH2_pops_ever_oscil )
02995 fprintf(ioQQQ," frac_new_oscil %.4f",frac_new_oscil);
02996 fprintf(ioQQQ,"\n");
02997 }
02998
02999 if( mole.nH2_TRACE >= mole.nH2_trace_full )
03000 {
03001 fprintf(ioQQQ,
03002 "H2 loop\t%li\tkase pop chng\t%i\tchem renorm fac\t%.4e\tortho/para ratio:\t%.3e\tfrac of pop in matrix: %.3f\n",
03003 loop_h2_pops,
03004 kase,
03005 H2_renorm_chemistry,
03006 h2.ortho_density / h2.para_density ,
03007 frac_matrix);
03008
03009
03010 if( iVibMaxChng_relative>=0 && iRotMaxChng_relative>=0 )
03011 fprintf(ioQQQ,
03012 "end loop %li H2 max rel chng=%.3e from %.3e to %.3e at v=%li J=%li\n\n",
03013 loop_h2_pops,
03014 PopChgMax_relative ,
03015 H2_old_populations[0][iVibMaxChng_relative][iRotMaxChng_relative],
03016 H2_populations[0][iVibMaxChng_relative][iRotMaxChng_relative],
03017 iVibMaxChng_relative , iRotMaxChng_relative
03018 );
03019 }
03020 }
03021
03022
03023
03024 H2_gs_rates();
03025
03026
03027 if( !lgConv_h2_soln && !conv.lgSearch )
03028 {
03029 conv.lgConvPops = false;
03030 strcpy( conv.chConvIoniz, "H2 pop cnv" );
03031 fprintf(ioQQQ,
03032 " H2_LevelPops: H2_populations not converged in %li tries; due to %s, old, new are %.4e %.4e, iteration %li zone %.2f.\n",
03033 loop_h2_pops,
03034 chReason,
03035 quant_old,
03036 quant_new ,
03037 iteration ,
03038 fnzone );
03039 ConvFail("pops","H2");
03040 }
03041
03042
03043
03044 for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi )
03045 {
03046 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
03047 {
03048 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
03049 {
03050 long int lim_elec_lo = 0;
03051
03052
03053
03054
03055 for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
03056 {
03057
03058
03059 long int nv = h2.nVib_hi[iElecLo];
03060 if( iElecLo==iElecHi )
03061 nv = iVibHi;
03062 for( iVibLo=0; iVibLo<=nv; ++iVibLo )
03063 {
03064 long nr = h2.nRot_hi[iElecLo][iVibLo];
03065 if( iElecLo==iElecHi && iVibHi==iVibLo )
03066 nr = iRotHi-1;
03067
03068 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
03069 {
03070 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].PopHi =
03071 H2_populations[iElecHi][iVibHi][iRotHi];
03072 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].PopLo =
03073 H2_populations[iElecLo][iVibLo][iRotLo];
03074 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].PopOpc =
03075 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].PopLo -
03076 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].PopHi*
03077 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].gLo /
03078 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].gHi;
03079 ASSERT(H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].PopHi >= 0. &&
03080 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].PopLo >= 0.);
03081
03082
03083
03084
03085 if( lgH2_line_exists[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] )
03086 {
03087
03088 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].cool = 0.;
03089 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].heat = 0.;
03090
03091
03092 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].phots =
03093 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Aul *
03094 (H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Pesc +
03095 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Pelec_esc) *
03096 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].PopHi;
03097
03098
03099 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].xIntensity =
03100 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].phots *
03101 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyErg;
03102
03103 if( iElecHi==0 )
03104 {
03106
03107
03108 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].ColOvTot = 1.;
03109 }
03110 else
03111 {
03112
03114 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].ColOvTot = 0.;
03115 }
03116 }
03117 }
03118 }
03119 }
03120 }
03121 }
03122 }
03123
03124
03125
03126
03127 hmi.H2_photodissoc_BigH2_H2s = 0.;
03128 hmi.H2_photodissoc_BigH2_H2g = 0.;
03129 hmi.H2_tripletdissoc_H2s =0.;
03130 hmi.H2_tripletdissoc_H2g =0.;
03131 hmi.H2_BigH2_H2g_av = 0.;
03132 hmi.H2_BigH2_H2s_av = 0.;
03133
03134 hmi.Average_collH2s_dissoc = 0.;
03135 hmi.Average_collH2g_dissoc = 0.;
03136
03137 iElec = 0;
03138 H2_BigH2_H2s = 0.;
03139 H2_BigH2_H2g = 0.;
03140 hmi.H2g_BigH2 =0;
03141 hmi.H2s_BigH2 = 0;
03142 hmi.H2_total_BigH2 =0;
03143 hmi.H2g_LTE_bigH2 =0.;
03144 hmi.H2s_LTE_bigH2 = 0.;
03145
03146
03147
03148
03149
03150
03151
03152 {
03153 static long ip_cut_off = -1;
03154 if( ip_cut_off < 0 )
03155 {
03156
03157 ip_cut_off = ipoint( 1.14 );
03158 }
03159
03160
03161
03162 flux_accum_photodissoc_BigH2_H2s = 0;
03163 ip_H2_level = ipoint( 1.07896 - 2.5 / EVRYD);
03164 for( i= ip_H2_level; i < ip_cut_off; ++i )
03165 {
03166 flux_accum_photodissoc_BigH2_H2s += ( rfield.flux[i-1] + rfield.ConInterOut[i-1]+
03167 rfield.outlin[i-1]+ rfield.outlin_noplot[i-1] );
03168 }
03169
03170
03171 for( iVib=0; iVib<=h2.nVib_hi[iElec]; ++iVib )
03172 {
03173 for( iRot=h2.Jlowest[iElec]; iRot<=h2.nRot_hi[iElec][iVib]; ++iRot )
03174 {
03175
03176
03177
03178
03179
03180
03181 if( energy_wn[0][iVib][iRot] > ENERGY_H2_STAR )
03182 {
03183 double arg;
03184 hmi.H2_photodissoc_BigH2_H2s +=
03185 H2_populations[iElec][iVib][iRot] * flux_accum_photodissoc_BigH2_H2s;
03186
03187 hmi.H2_tripletdissoc_H2s += H2_populations[iElec][iVib][iRot] * 3.f*secondaries.x12tot;
03188
03189
03190 H2_BigH2_H2s += H2_populations[iElec][iVib][iRot];
03191
03192
03193 hmi.H2_BigH2_H2s_av += (H2_populations[iElec][iVib][iRot] * energy_wn[0][iVib][iRot]);
03194
03195
03196 hmi.Average_collH2s_dissoc += H2_populations[iElec][iVib][iRot] * H2_coll_dissoc_rate_coef_H2[iVib][iRot];
03197
03198
03199 arg = exp_disoc/SDIV(H2_Boltzmann[0][iVib][iRot]);
03200 if( arg > 0. )
03201 {
03202
03203 hmi.H2s_LTE_bigH2 += H2_populations[0][iVib][iRot]*SAHA/SDIV(phycon.te32*arg)*
03204 (H2_stat[0][iVib][iRot]/(2.*2.))*3.634e-5;
03205 }
03206 }
03207 else
03208 {
03209 double arg;
03210
03211 flux_accum_photodissoc_BigH2_H2g = 0;
03212
03213 ip_H2_level = ipoint( 1.07896 - energy_wn[0][iVib][iRot] * WAVNRYD);
03214
03215 for( i= ip_H2_level; i < ip_cut_off; ++i )
03216 {
03217 flux_accum_photodissoc_BigH2_H2g += ( rfield.flux[i-1] + rfield.ConInterOut[i-1]+
03218 rfield.outlin[i-1]+ rfield.outlin_noplot[i-1] );
03219 }
03220
03221 hmi.H2_photodissoc_BigH2_H2g +=
03222 H2_populations[iElec][iVib][iRot] * flux_accum_photodissoc_BigH2_H2g;
03223
03224
03225 hmi.H2_tripletdissoc_H2g +=
03226 H2_populations[iElec][iVib][iRot] * 3.f*secondaries.x12tot;
03227
03228
03229 H2_BigH2_H2g += H2_populations[iElec][iVib][iRot];
03230
03231
03232 hmi.H2_BigH2_H2g_av += (H2_populations[iElec][iVib][iRot] * energy_wn[0][iVib][iRot]);
03233
03234
03235 hmi.Average_collH2g_dissoc += H2_populations[iElec][iVib][iRot] * H2_coll_dissoc_rate_coef_H2[iVib][iRot];
03236
03237
03238 arg = exp_disoc/SDIV(H2_Boltzmann[0][iVib][iRot]);
03239 if( arg > 0. )
03240 {
03241 hmi.H2g_LTE_bigH2 += H2_populations[0][iVib][iRot]*(SAHA/SDIV(phycon.te32*arg)*
03242 (H2_stat[0][iVib][iRot]/(2.*2.))*3.634e-5);
03243 }
03244 }
03245 }
03246 }
03247 }
03248 hmi.H2g_BigH2 = (float)H2_BigH2_H2g;
03249 hmi.H2s_BigH2 = (float)H2_BigH2_H2s;
03250 hmi.H2_total_BigH2 =hmi.H2g_BigH2+hmi.H2s_BigH2;
03251
03252
03253 hmi.H2_BigH2_H2s_av = hmi.H2_BigH2_H2s_av / H2_BigH2_H2s;
03254
03255 hmi.H2_BigH2_H2g_av = hmi.H2_BigH2_H2g_av / H2_BigH2_H2g;
03256
03257
03258
03259
03260
03261
03262
03263 hmi.H2_photodissoc_BigH2_H2s = hmi.H2_photodissoc_BigH2_H2s / SDIV(H2_BigH2_H2s) * H2_DISS_ALLISON_DALGARNO;
03264 hmi.H2_photodissoc_BigH2_H2g = hmi.H2_photodissoc_BigH2_H2g / SDIV(H2_BigH2_H2g) * H2_DISS_ALLISON_DALGARNO;
03265 hmi.H2_tripletdissoc_H2g = hmi.H2_tripletdissoc_H2g/SDIV(H2_BigH2_H2g);
03266 hmi.H2_tripletdissoc_H2s = hmi.H2_tripletdissoc_H2s/SDIV(H2_BigH2_H2s);
03267 hmi.Average_collH2g_dissoc = hmi.Average_collH2g_dissoc /SDIV(H2_BigH2_H2g);
03268 hmi.Average_collH2s_dissoc = hmi.Average_collH2s_dissoc /SDIV(H2_BigH2_H2s);
03269 hmi.H2s_LTE_bigH2 = hmi.H2s_LTE_bigH2/SDIV(H2_BigH2_H2s);
03270 hmi.H2g_LTE_bigH2 = hmi.H2g_LTE_bigH2/SDIV(H2_BigH2_H2g);
03271
03272
03273
03274 {
03275 double sumpop1 = 0.;
03276 double sumpopA1 = 0.;
03277 double sumpopcollH2O = 0.;
03278 double sumpopcollH2p = 0.;
03279 double sumpopcollH = 0.;
03280 double sumpop2 = 0.;
03281 double sumpopcollH2O_excit = 0.;
03282 double sumpopcollH2p_excit = 0.;
03283 double sumpopcollH_excit = 0.;
03284 double sumpop3 = 0.;
03285
03286
03287 iElecLo = 0;
03288 for( iVibHi=0; iVibHi<=h2.nVib_hi[0]; ++iVibHi )
03289 {
03290 long nr1 = h2.nRot_hi[0][iVibHi];
03291 for( iRotHi=h2.Jlowest[0]; iRotHi<=nr1; ++iRotHi )
03292 {
03293 for( iVibLo=0; iVibLo<=h2.nVib_hi[0]; ++iVibLo )
03294 {
03295 long nr2 = h2.nRot_hi[0][iVibLo];
03296 for( iRotLo=h2.Jlowest[0]; iRotLo<=nr2; ++iRotLo )
03297 {
03298
03299 if( (energy_wn[0][iVibHi][iRotHi] > ENERGY_H2_STAR) && (energy_wn[0][iVibLo][iRotLo] < ENERGY_H2_STAR ))
03300 {
03301
03302
03303
03304 if((H2_lgOrtho[0][iVibHi][iRotHi]-H2_lgOrtho[0][iVibLo][iRotLo])==0)
03305 {
03306 sumpop2 += H2_populations[0][iVibHi][iRotHi];
03307 sumpop3 += H2_populations[0][iVibLo][iRotLo];
03308 sumpopcollH += H2_populations[0][iVibHi][iRotHi]*H2_CollRate[0][iVibHi][iRotHi][iVibLo][iRotLo];
03309 sumpopcollH2O += H2_populations[0][iVibHi][iRotHi]*H2_CollRate[2][iVibHi][iRotHi][iVibLo][iRotLo];
03310 sumpopcollH2p += H2_populations[0][iVibHi][iRotHi]*H2_CollRate[3][iVibHi][iRotHi][iVibLo][iRotLo];
03311
03312 sumpopcollH_excit += H2_populations[0][iVibLo][iRotLo]*H2_CollRate[0][iVibHi][iRotHi][iVibLo][iRotLo]*H2_stat[0][iVibHi][iRotHi] / H2_stat[0][iVibLo][iRotLo] *
03313 H2_Boltzmann[0][iVibHi][iRotHi] /SDIV( H2_Boltzmann[0][iVibLo][iRotLo] );
03314 sumpopcollH2O_excit += H2_populations[0][iVibLo][iRotLo]*H2_CollRate[2][iVibHi][iRotHi][iVibLo][iRotLo]*H2_stat[0][iVibHi][iRotHi] / H2_stat[0][iVibLo][iRotLo] *
03315 H2_Boltzmann[0][iVibHi][iRotHi] /SDIV( H2_Boltzmann[0][iVibLo][iRotLo] );
03316 sumpopcollH2p_excit += H2_populations[0][iVibLo][iRotLo]*H2_CollRate[3][iVibHi][iRotHi][iVibLo][iRotLo]*H2_stat[0][iVibHi][iRotHi] / H2_stat[0][iVibLo][iRotLo] *
03317 H2_Boltzmann[0][iVibHi][iRotHi] /SDIV( H2_Boltzmann[0][iVibLo][iRotLo] );
03318
03319
03320
03321 if( (abs((iRotHi-iRotLo) ))==2 || (iRotHi==iRotLo ) )
03322
03323 {
03324 sumpop1 += H2Lines[0][iVibHi][iRotHi][0][iVibLo][iRotLo].PopHi;
03325 sumpopA1 += H2Lines[0][iVibHi][iRotHi][0][iVibLo][iRotLo].PopHi*
03326 H2Lines[0][iVibHi][iRotHi][0][iVibLo][iRotLo].Aul;
03327 }
03328 }
03329 }
03330 }
03331 }
03332 }
03333 }
03334 hmi.Average_A = sumpopA1/SDIV(sumpop1);
03335 hmi.Average_collH = sumpopcollH/SDIV(sumpop2);
03336 hmi.Average_collH2 = (sumpopcollH2O+sumpopcollH2p)/SDIV(sumpop2);
03337 hmi.Average_collH_excit = sumpopcollH_excit/SDIV(sumpop3);
03338 hmi.Average_collH2_excit = (sumpopcollH2O_excit+sumpopcollH2p_excit)/SDIV(sumpop3);
03339
03340
03341
03342
03343
03344
03345
03346 }
03347 if( mole.nH2_TRACE >= mole.nH2_trace_full|| (trace.lgTrace && trace.lgTr_H2_Mole) )
03348 {
03349 fprintf(ioQQQ," H2_LevelPops exit2 Sol dissoc %.2e (TH85 %.2e)",
03350 hmi.H2_Solomon_dissoc_rate_BigH2_H2g +
03351 hmi.H2_Solomon_dissoc_rate_BigH2_H2s ,
03352 hmi.H2_Solomon_dissoc_rate_TH85_H2g);
03353
03354
03355
03356 fprintf(ioQQQ," H2g Sol %.2e H2s Sol %.2e",
03357 hmi.H2_Solomon_dissoc_rate_used_H2g ,
03358 hmi.H2_Solomon_dissoc_rate_BigH2_H2s );
03359
03360
03361 fprintf(ioQQQ," H2g->H2s %.2e (TH85 %.2e)",
03362 hmi.H2_H2g_to_H2s_rate_BigH2 ,
03363 hmi.H2_H2g_to_H2s_rate_TH85);
03364
03365
03366
03367 fprintf(ioQQQ," H2 con diss %.2e (TH85 %.2e)\n",
03368 hmi.H2_photodissoc_BigH2_H2s ,
03369 hmi.H2_photodissoc_TH85);
03370 }
03371 else if( mole.nH2_TRACE )
03372 {
03373 fprintf(ioQQQ," H2_LevelPops exit1 %8.2f loops:%3li H2/H:%.3e Sol dis old %.3e new %.3e",
03374 fnzone ,
03375 loop_h2_pops ,
03376 hmi.H2_total / dense.gas_phase[ipHYDROGEN],
03377 old_solomon_rate,
03378 hmi.H2_Solomon_dissoc_rate_BigH2_H2g );
03379 fprintf(ioQQQ,"\n");
03380 }
03381
03382
03383
03384
03385
03386
03387 ++nCallH2_this_iteration;
03388
03389
03390
03391 ++h2.nCallH2_this_zone;
03392
03393
03394
03395
03396 nXLevelsMatrix = nXLevelsMatrix_save;
03397
03398
03399
03400
03401 if( nCallH2_this_iteration && nzone != nzone_nlevel_set )
03402 {
03403
03404 # define FRAC 0.99999
03405
03406 double sum_pop = 0.;
03407 nEner = 0;
03408 iElec = 0;
03409 # define PRT false
03410 if( PRT ) fprintf(ioQQQ,"DEBUG pops ");
03411 while( nEner < nLevels_per_elec[0] && sum_pop/hmi.H2_total < FRAC )
03412 {
03413
03414
03415 ip = H2_ipX_ener_sort[nEner];
03416 iVib = ipVib_H2_energy_sort[ip];
03417 iRot = ipRot_H2_energy_sort[ip];
03418 sum_pop += H2_old_populations[iElec][iVib][iRot];
03419 if( PRT ) fprintf(ioQQQ,"\t%.3e ", H2_old_populations[iElec][iVib][iRot]);
03420 ++nEner;
03421 }
03422 if( PRT ) fprintf(ioQQQ,"\n ");
03423 nzone_nlevel_set = nzone;
03424
03425
03426
03427 # undef FRAC
03428 # undef PRT
03429 }
03430
03431 DEBUG_EXIT( "H2_LevelPops()" );
03432 return;
03433 }
03434
03435
03436
03437
03438 void H2_Cooling(const char *chRoutine)
03439 {
03440 long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo;
03441 double heatone ,
03442 rate_dn_heat,
03443 rate_up_cool;
03444 long int nColl,
03445 ipHi, ipLo;
03446 double Big1;
03447 long int ipVib_big_hi,ipVib_big_lo ,ipRot_big_hi ,ipRot_big_lo;
03448 static double old_HeatH2Dexc=0.;
03449
03450 # ifdef DEBUG_DIS_DEAT
03451 double heatbig;
03452 long int iElecBig , iVibBig , iRotBig;
03453 # endif
03454
03455 DEBUG_ENTRY( "H2_Cooling()" );
03456
03457
03458
03459
03460
03461 if( !h2.lgH2ON || !nCallH2_this_iteration )
03462 {
03463 hmi.HeatH2Dexc_BigH2 = 0.;
03464 hmi.HeatH2Dish_BigH2 = 0.;
03465 hmi.deriv_HeatH2Dexc_BigH2 = 0.;
03466 DEBUG_EXIT( "H2_Cooling()" );
03467 return;
03468 }
03469
03470 old_HeatH2Dexc = hmi.HeatH2Dexc_BigH2;
03471
03472 hmi.HeatH2Dish_BigH2 = 0.;
03473 heatone = 0.;
03474 # ifdef DEBUG_DIS_DEAT
03475 heatbig = 0.;
03476 iElecBig = -1;
03477 iVibBig = -1;
03478 iRotBig = -1;
03479 # endif
03480
03481 for( iElecHi=1; iElecHi<mole.n_h2_elec_states; ++iElecHi )
03482 {
03483 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
03484 {
03485 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
03486 {
03487
03488 heatone = H2_populations[iElecHi][iVibHi][iRotHi] *
03489 H2_dissprob[iElecHi][iVibHi][iRotHi] *
03490 H2_disske[iElecHi][iVibHi][iRotHi];
03491 hmi.HeatH2Dish_BigH2 += heatone;
03492 # ifdef DEBUG_DIS_DEAT
03493 if( heatone > heatbig )
03494 {
03495 heatbig = heatone;
03496 iElecBig = iElecHi;
03497 iVibBig = iVibHi;
03498 iRotBig = iRotHi;
03499 }
03500 # endif
03501 }
03502 }
03503 }
03504 # ifdef DEBUG_DIS_DEAT
03505 fprintf(ioQQQ,"DEBUG H2 dis heat\t%.2f\t%.3f\t%li\t\t%li\t%li\n",
03506 fnzone ,
03507 heatbig / SDIV( hmi.HeatH2Dish_BigH2 ) ,
03508 iElecBig ,
03509 iVibBig ,
03510 iRotBig );
03511 # endif
03512
03513
03514 hmi.HeatH2Dish_BigH2 *= EN1EV;
03515
03516
03517
03518
03519
03520 hmi.HeatH2Dexc_BigH2 = 0.;
03521
03522
03523
03524
03525
03526 iElecHi = 0;
03527 iElecLo = 0;
03528 Big1 = 0.;
03529 ipVib_big_hi = 0;
03530 ipVib_big_lo = 0;
03531 ipRot_big_hi = 0;
03532 ipRot_big_lo = 0;
03533
03534 hmi.deriv_HeatH2Dexc_BigH2 = 0.;
03535 for( ipHi=1; ipHi<nLevels_per_elec[iElecHi]; ++ipHi )
03536 {
03537 long int ip = H2_ipX_ener_sort[ipHi];
03538 iVibHi = ipVib_H2_energy_sort[ip];
03539 iRotHi = ipRot_H2_energy_sort[ip];
03540 if( iVibHi > VIB_COLLID )
03541 continue;
03542
03543 for( ipLo=0; ipLo<ipHi; ++ipLo )
03544 {
03545 double oneline;
03546 ip = H2_ipX_ener_sort[ipLo];
03547 iVibLo = ipVib_H2_energy_sort[ip];
03548 iRotLo = ipRot_H2_energy_sort[ip];
03549 if( iVibLo > VIB_COLLID)
03550 continue;
03551
03552 rate_dn_heat = 0.;
03553 rate_up_cool = 0.;
03554
03555
03556 for( nColl=0; nColl<N_X_COLLIDER; ++nColl )
03557 {
03558 rate_dn_heat +=
03559 H2_populations[iElecHi][iVibHi][iRotHi] *
03560 H2_CollRate[nColl][iVibHi][iRotHi][iVibLo][iRotLo]*
03561 collider_density[nColl];
03562
03563
03564 rate_up_cool +=
03565 H2_populations[iElecLo][iVibLo][iRotLo] *
03566
03567 H2_CollRate[nColl][iVibHi][iRotHi][iVibLo][iRotLo]*
03568 collider_density[nColl]*
03569
03570 H2_stat[iElecHi][iVibHi][iRotHi] / H2_stat[iElecLo][iVibLo][iRotLo] *
03571 H2_Boltzmann[iElecHi][iVibHi][iRotHi] /
03572 SDIV( H2_Boltzmann[iElecLo][iVibLo][iRotLo] );
03573 }
03574
03575
03576
03577 oneline = (rate_dn_heat - rate_up_cool)*
03578 (energy_wn[iElecHi][iVibHi][iRotHi] - energy_wn[iElecLo][iVibLo][iRotLo]) *
03579 ERG1CM;
03580 hmi.HeatH2Dexc_BigH2 += oneline;
03581
03582
03583
03584
03585 hmi.deriv_HeatH2Dexc_BigH2 += (float)(oneline *energy_wn[iElecHi][iVibHi][iRotHi]);
03586
03587
03588 # define H2COOLDEBUG false
03589 if(H2COOLDEBUG)
03590 {
03591 if( fabs(oneline) > fabs(Big1 ) )
03592 {
03593 Big1 = oneline;
03594 ipVib_big_hi = iVibHi;
03595 ipVib_big_lo = iVibLo;
03596 ipRot_big_hi = iRotHi;
03597 ipRot_big_lo = iRotLo;
03598 }
03599 }
03600
03601
03602 ASSERT(
03603 (rate_up_cool==0 && rate_dn_heat==0) ||
03604 (energy_wn[iElecHi][iVibHi][iRotHi] > energy_wn[iElecLo][iVibLo][iRotLo]) );
03605 }
03606 }
03607 if(H2COOLDEBUG)
03608 {
03609 static long int nCount=0;
03610 if(nCount>55000 )
03611 {
03612 fprintf(ioQQQ,"DEBUG H2_Cooling %s, Te %.3e H+/0 %.2e H2 %.3e Sol %.3e grn %.3e total coll %.2e, frac 1 line %.2e %li %li %li %li ",
03613 chRoutine,
03614 phycon.te,
03615 dense.xIonDense[ipHYDROGEN][1]/SDIV(dense.xIonDense[ipHYDROGEN][1]),
03616 hmi.H2_total,
03617 hmi.H2_Solomon_dissoc_rate_BigH2_H2g,
03618 hmi.rate_grain_h2_J1_to_J0 ,
03619 hmi.HeatH2Dexc_BigH2 ,
03620 Big1/hmi.HeatH2Dexc_BigH2 ,
03621 ipVib_big_hi , ipRot_big_hi , ipVib_big_lo , ipRot_big_lo );
03622
03623 for( iRotLo=0; iRotLo<7; ++iRotLo )
03624 {
03625 fprintf(ioQQQ,"\t%li\t%.2e" , iRotLo ,
03626 H2_populations[0][0][iRotLo]/hmi.H2_total );
03627 }
03628 fprintf(ioQQQ,"\t%li\n",nCount);
03629 }
03630 ++nCount;
03631 if( nCount > 55402 )
03632 {
03633 cdEXIT( EXIT_FAILURE );
03634 }
03635 else if( nCount > 55390 )
03636 {
03637 mole.nH2_TRACE = mole.nH2_trace_full;
03638 trace.npsbug = 1;
03639 trace.lgTrace = true;
03640 geometry.nprint = 1;
03641 iterations.IterPrnt[0] = 1;
03642 }
03643 }
03644
03645
03646 if( PRT_POPS )
03647 fprintf(ioQQQ,
03648 " DEBUG H2 heat fnzone\t%.2f\trenrom\t%.3e\tte\t%.4e\tdexc\t%.3e\theat/tot\t%.3e\n",
03649 fnzone ,
03650 H2_renorm_chemistry ,
03651 phycon.te ,
03652 hmi.HeatH2Dexc_BigH2,
03653 hmi.HeatH2Dexc_BigH2/thermal.ctot);
03654
03655
03656
03657 hmi.deriv_HeatH2Dexc_BigH2 /= (float)POW2(phycon.te_wn);
03658
03659 {
03660
03661 enum {DEBUG_LOC=false};
03662
03663 if( DEBUG_LOC && (fabs(hmi.HeatH2Dexc_BigH2) > SMALLFLOAT) )
03664 {
03665 int iVib = 0;
03666
03667
03668
03669
03670
03671
03672
03673
03674
03675 iElecHi = iElecLo = 0;
03676 iVibHi = iVibLo = 0;
03677 iRotHi = 3;
03678 iRotLo = 1;
03679 rate_dn_heat = rate_up_cool = 0.;
03680
03681 for( nColl=0; nColl<N_X_COLLIDER; ++nColl )
03682 {
03683 rate_dn_heat +=
03684 H2_populations[iElecHi][iVibHi][iRotHi] *
03685 H2_CollRate[nColl][iVibHi][iRotHi][iVibLo][iRotLo]*
03686 collider_density[nColl];
03687
03688
03689 rate_up_cool +=
03690 H2_populations[iElecLo][iVibLo][iRotLo] *
03691
03692 H2_CollRate[nColl][iVibHi][iRotHi][iVibLo][iRotLo]*
03693 collider_density[nColl]*
03694
03695 H2_stat[iElecHi][iVibHi][iRotHi] / H2_stat[iElecLo][iVibLo][iRotLo] *
03696 H2_Boltzmann[iElecHi][iVibHi][iRotHi] /
03697 SDIV( H2_Boltzmann[iElecLo][iVibLo][iRotLo] );
03698 }
03699
03700 fprintf(ioQQQ," H2_cooling pop31\t%.3e\tdn up 31\t%.3e\t%.3e\n",
03701 H2_populations[0][iVib][3]/H2_populations[0][iVib][1],
03702 rate_dn_heat,rate_up_cool
03703 );
03704 }
03705 }
03706 # undef H2COOLDEBUG
03707 {
03708
03709 enum {DEBUG_LOC=false};
03710
03711 if( DEBUG_LOC )
03712 {
03713 static long nzdone=-1 , nzincre;
03714 if( nzone!=nzdone )
03715 {
03716 nzdone = nzone;
03717 nzincre = -1;
03718 }
03719 ++nzincre;
03720 fprintf(ioQQQ," H2 nz\t%.2f\tnzinc\t%li\tTe\t%.4e\tH2\t%.3e\tcXH\t%.2e\tdcXH/dt%.2e\tDish\t%.2e \n",
03721 fnzone,
03722 nzincre,
03723 phycon.te,
03724 hmi.H2_total ,
03725 hmi.HeatH2Dexc_BigH2,
03726 hmi.deriv_HeatH2Dexc_BigH2 ,
03727 hmi.HeatH2Dish_BigH2);
03728
03729 }
03730 }
03731
03732
03733
03734
03735 if( 1 || nzone <1 || old_HeatH2Dexc==0. || nCallH2_this_iteration <2)
03736 {
03737 old_HeatH2Dexc = hmi.HeatH2Dexc_BigH2;
03738 }
03739 else
03740 {
03741 hmi.HeatH2Dexc_BigH2 = (hmi.HeatH2Dexc_BigH2+old_HeatH2Dexc)/2.f;
03742 old_HeatH2Dexc = hmi.HeatH2Dexc_BigH2;
03743 }
03744 if( mole.nH2_TRACE >= mole.nH2_trace_full )
03745 fprintf(ioQQQ,
03746 " H2_Cooling Ctot\t%.4e\t HeatH2Dish_BigH2 \t%.4e\t HeatH2Dexc_BigH2 \t%.4e\n" ,
03747 thermal.ctot ,
03748 hmi.HeatH2Dish_BigH2 ,
03749 hmi.HeatH2Dexc_BigH2 );
03750
03751
03752
03753
03754
03755
03756
03757
03758 if( conv.lgSearch )
03759 hmi.HeatH2Dexc_BigH2 = 0.;
03760
03761 DEBUG_EXIT( "H2_Cooling()" );
03762 return;
03763
03764 }
03765
03766
03767
03768
03769 double cdH2_colden( long iVib , long iRot )
03770 {
03771
03772
03773
03774
03775
03776
03777 if( iVib < 0 )
03778 {
03779 if( iRot==0 )
03780 {
03781
03782 return( h2.ortho_colden + h2.para_colden );
03783 }
03784 else if( iRot==1 )
03785 {
03786
03787 return h2.ortho_colden;
03788 }
03789 else if( iRot==2 )
03790 {
03791
03792 return h2.para_colden;
03793 }
03794 else
03795 {
03796 fprintf(ioQQQ," iRot must be 0 (total), 1 (ortho), or 2 (para), returning -1.\n");
03797 return -1.;
03798 }
03799 }
03800 else if( h2.lgH2ON )
03801 {
03802
03803
03804 int iElec = 0;
03805 if( iRot <0 || iVib >h2.nVib_hi[iElec] || iRot > h2.nRot_hi[iElec][iVib])
03806 {
03807 fprintf(ioQQQ," iVib and iRot must lie within X, returning -2.\n");
03808 fprintf(ioQQQ," iVib must be <= %li and iRot must be <= %li.\n",
03809 h2.nVib_hi[iElec],h2.nRot_hi[iElec][iVib]);
03810 return -2.;
03811 }
03812 else
03813 {
03814 return H2_X_colden[iVib][iRot];
03815 }
03816 }
03817
03818 else
03819 return -1;
03820 }
03821
03822
03823 void H2_Colden( const char *chLabel )
03824 {
03825 long int iVib , iRot;
03826
03827
03828 if( !h2.lgH2ON )
03829 return;
03830
03831 DEBUG_ENTRY( "H2_Colden()" );
03832
03833 if( strcmp(chLabel,"ZERO") == 0 )
03834 {
03835
03836 for( iVib = 0; iVib <= h2.nVib_hi[0]; ++iVib )
03837 {
03838 for( iRot=h2.Jlowest[0]; iRot<=h2.nRot_hi[0][iVib]; ++iRot )
03839 {
03840
03841 H2_X_colden[iVib][iRot] = 0.;
03842 H2_X_colden_LTE[iVib][iRot] = 0.;
03843 }
03844 }
03845 }
03846
03847 else if( strcmp(chLabel,"ADD ") == 0 )
03848 {
03849
03850 for( iVib = 0; iVib <= h2.nVib_hi[0]; ++iVib )
03851 {
03852 for( iRot=h2.Jlowest[0]; iRot<=h2.nRot_hi[0][iVib]; ++iRot )
03853 {
03854
03855 H2_X_colden[iVib][iRot] += (float)(H2_populations[0][iVib][iRot]*radius.drad_x_fillfac);
03856
03857
03858 H2_X_colden_LTE[iVib][iRot] += (float)(H2_populations_LTE[0][iVib][iRot]*
03859 hmi.H2_total*radius.drad_x_fillfac);
03860 }
03861 }
03862 }
03863
03864
03865 else if( strcmp(chLabel,"PRIN") != 0 )
03866 {
03867 fprintf( ioQQQ, " H2_Colden does not understand the label %s\n",
03868 chLabel );
03869 puts( "[Stop in H2_Colden]" );
03870 cdEXIT(EXIT_FAILURE);
03871 }
03872
03873 DEBUG_EXIT( "H2_Colden()" );
03874
03875 return;
03876 }
03877
03878
03879 double H2_DR(void)
03880 {
03881 return BIGFLOAT;
03882 }
03883
03884
03885 void H2_RT_OTS( void )
03886 {
03887
03888 long int iElecHi , iVibHi , iRotHi , iElecLo , iVibLo , iRotLo;
03889
03890
03891 if( !h2.lgH2ON || !h2.nCallH2_this_zone )
03892 return;
03893
03894 DEBUG_ENTRY( "H2_RT_OTS()" );
03895
03896
03897 for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi )
03898 {
03899 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
03900 {
03901 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
03902 {
03903 long int lim_elec_lo = 0;
03904 for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
03905 {
03906
03907
03908 long int nv = h2.nVib_hi[iElecLo];
03909 if( iElecLo==iElecHi )
03910 nv = iVibHi;
03911 for( iVibLo=0; iVibLo<=nv; ++iVibLo )
03912 {
03913 long nr = h2.nRot_hi[iElecLo][iVibLo];
03914 if( iElecLo==iElecHi && iVibHi==iVibLo )
03915 nr = iRotHi-1;
03916
03917 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
03918 {
03919
03920
03921 if( lgH2_line_exists[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] )
03922 {
03923 if( iElecHi==0 )
03924 {
03925
03926 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].ots =
03927 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Aul *
03928 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].PopHi *
03929 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Pdest;
03930
03931
03932 RT_OTS_AddLine(
03933 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].ots,
03934 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].ipCont );
03935 }
03936 }
03937 }
03938 }
03939 }
03940 }
03941 }
03942 }
03943
03944 DEBUG_EXIT( "H2_RT_OTS()" );
03945
03946 return;
03947 }
03948
03949