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