00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015 #include "cddefines.h"
00016 #include "phycon.h"
00017 #include "mole.h"
00018 #include "hmi.h"
00019 #include "taulines.h"
00020 #include "h2.h"
00021 #include "h2_priv.h"
00022
00023
00024
00025 void H2_Solomon_rate( void )
00026 {
00027
00028 long int iElecLo , iElecHi , iVibLo , iVibHi , iRotLo , iRotHi;
00029
00030 DEBUG_ENTRY( "H2_Solomon_rate()" );
00031
00032
00033 iElecLo = 0;
00034
00035
00036
00037
00038
00039
00040
00041 hmi.H2_Solomon_dissoc_rate_BigH2_H2g = 0.;
00042 hmi.H2_Solomon_dissoc_rate_BigH2_H2s = 0.;
00043
00044
00045 hmi.H2_Solomon_elec_decay_H2g = 0.;
00046 hmi.H2_Solomon_elec_decay_H2s = 0.;
00047
00048
00049
00050
00051
00052
00053 for( iElecHi=1; iElecHi<mole.n_h2_elec_states; ++iElecHi )
00054 {
00055 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
00056 {
00057 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
00058 {
00059 double factor = (double)H2_dissprob[iElecHi][iVibHi][iRotHi]/
00060 H2_rad_rate_out[iElecHi][iVibHi][iRotHi];
00061 double H2popHi = H2_populations[iElecHi][iVibHi][iRotHi];
00062
00063
00064
00065
00066
00067 iElecLo = 0;
00068
00069 for( iVibLo=0; iVibLo<=h2.nVib_hi[iElecLo]; ++iVibLo )
00070 {
00071 long nr = h2.nRot_hi[iElecLo][iVibLo];
00072
00073 mb6ci lgH2le = lgH2_line_exists.ptr(iElecHi,iVibHi,iRotHi,iElecLo,iVibLo,h2.Jlowest[iElecLo]);
00074 mt6ci H2L = H2Lines.ptr(iElecHi,iVibHi,iRotHi,iElecLo,iVibLo,h2.Jlowest[iElecLo]);
00075 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
00076 {
00077 if( *lgH2le++ )
00078 {
00079
00080
00081
00082 double rate_up_cont =
00083 H2_populations[iElecLo][iVibLo][iRotLo]*
00084 H2L->Emis->pump*factor;
00085
00086
00087 double elec_decay =
00088 H2popHi*
00089 H2L->Emis->Aul*(H2L->Emis->Pesc+H2L->Emis->Pdest+H2L->Emis->Pelec_esc);
00090
00091 if( energy_wn[0][iVibLo][iRotLo] > ENERGY_H2_STAR )
00092 {
00093
00094
00095 hmi.H2_Solomon_dissoc_rate_BigH2_H2s += rate_up_cont;
00096
00097 hmi.H2_Solomon_elec_decay_H2s += elec_decay;
00098 }
00099 else
00100 {
00101
00102
00103 hmi.H2_Solomon_dissoc_rate_BigH2_H2g += rate_up_cont;
00104
00105 hmi.H2_Solomon_elec_decay_H2g += elec_decay;
00106 }
00107 }
00108 ++H2L;
00109 }
00110 }
00111 }
00112 }
00113 }
00114
00115
00116
00117 if( hmi.H2_total > SMALLFLOAT )
00118 {
00119 hmi.H2_Solomon_elec_decay_H2g /= SDIV( H2_sum_excit_elec_den );
00120 hmi.H2_Solomon_elec_decay_H2s /= SDIV( H2_sum_excit_elec_den );
00121
00122
00123 hmi.H2_Solomon_dissoc_rate_BigH2_H2s = hmi.H2_Solomon_dissoc_rate_BigH2_H2s / SDIV(H2_den_s);
00124
00125
00126 hmi.H2_Solomon_dissoc_rate_BigH2_H2g = hmi.H2_Solomon_dissoc_rate_BigH2_H2g / SDIV(H2_den_g);
00127
00128 }
00129 else
00130 {
00131 hmi.H2_Solomon_dissoc_rate_BigH2_H2s = 0;
00132 hmi.H2_Solomon_dissoc_rate_BigH2_H2g = 0;
00133
00134 }
00135
00136
00137
00138
00139
00140
00141
00142 hmi.H2_H2g_to_H2s_rate_BigH2 = 0.;
00143 return;
00144 }
00145
00146
00147 void H2_gs_rates( void )
00148 {
00149 long int ipLoX , ip , iRotLoX , iVibLoX ,
00150 iElecHi , iVibHi , ipOther , iRotOther,
00151 iRotHi , iVibOther;
00152
00153 DEBUG_ENTRY( "H2_gs_rates()" );
00154
00155
00156 hmi.H2_H2g_to_H2s_rate_BigH2 = 0.;
00157
00158
00159 for( ipLoX=0; ipLoX < nEner_H2_ground; ++ipLoX )
00160 {
00161 ip = H2_ipX_ener_sort[ipLoX];
00162 iRotLoX = ipRot_H2_energy_sort[ip];
00163 iVibLoX = ipVib_H2_energy_sort[ip];
00164
00165
00166 for( iElecHi=1; iElecHi<mole.n_h2_elec_states; ++iElecHi )
00167 {
00168 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
00169 {
00170 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
00171 {
00172 if( lgH2_line_exists[iElecHi][iVibHi][iRotHi][0][iVibLoX][iRotLoX] )
00173 {
00174
00175
00176
00177 double rate_up_cont =
00178 H2_populations[0][iVibLoX][iRotLoX]*
00179 H2Lines[iElecHi][iVibHi][iRotHi][0][iVibLoX][iRotLoX].Emis->pump;
00180
00181 double decay_star = H2_rad_rate_out[iElecHi][iVibHi][iRotHi] - H2_dissprob[iElecHi][iVibHi][iRotHi];
00182
00183
00184
00185 for( ipOther=0; ipOther < nEner_H2_ground; ++ipOther )
00186 {
00187 ip = H2_ipX_ener_sort[ipOther];
00188 iRotOther = ipRot_H2_energy_sort[ip];
00189 iVibOther = ipVib_H2_energy_sort[ip];
00190 if( lgH2_line_exists[iElecHi][iVibHi][iRotHi][0][iVibOther][iRotOther] )
00191 {
00192 decay_star -=
00193 H2Lines[iElecHi][iVibHi][iRotHi][0][iVibOther][iRotOther].Emis->Aul*(
00194 H2Lines[iElecHi][iVibHi][iRotHi][0][iVibOther][iRotOther].Emis->Pesc+
00195 H2Lines[iElecHi][iVibHi][iRotHi][0][iVibOther][iRotOther].Emis->Pdest+
00196 H2Lines[iElecHi][iVibHi][iRotHi][0][iVibOther][iRotOther].Emis->Pelec_esc);
00197 }
00198 }
00199
00200
00201 decay_star = MAX2(0., decay_star)/SDIV(H2_rad_rate_out[iElecHi][iVibHi][iRotHi]);
00202 hmi.H2_H2g_to_H2s_rate_BigH2 += rate_up_cont*decay_star;
00203
00204 }
00205 }
00206 }
00207 }
00208 }
00209
00210
00211 hmi.H2_H2g_to_H2s_rate_BigH2 /= SDIV( H2_den_g );
00212 return;
00213 }
00214
00215
00216
00217 void H2_zero_pops_too_low( void )
00218 {
00219
00220 long int iElec, iElecHi, iElecLo, iVib , iVibLo , iVibHi ,
00221 iRot , iRotHi , iRotLo;
00222
00223 DEBUG_ENTRY( "H2_zero_pops_too_low()" );
00224
00225
00226 for( iElec=0; iElec<mole.n_h2_elec_states; ++iElec )
00227 {
00228 for( iVib=0; iVib<=h2.nVib_hi[iElec]; ++iVib )
00229 {
00230 for( iRot=h2.Jlowest[iElec]; iRot<=h2.nRot_hi[iElec][iVib]; ++iRot )
00231 {
00232
00233
00234 H2_old_populations[iElec][iVib][iRot] =
00235 (realnum)H2_populations_LTE[iElec][iVib][iRot]*hmi.H2_total;
00236 H2_populations[iElec][iVib][iRot] = H2_old_populations[iElec][iVib][iRot];
00237 }
00238 }
00239 }
00240
00241
00242
00243 for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi )
00244 {
00245 pops_per_elec[iElecHi] = 0.;
00246 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
00247 {
00248 pops_per_vib[iElecHi][iVibHi] = 0.;
00249 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
00250 {
00251 long int lim_elec_lo = 0;
00252
00253
00254 H2Lines[iElecHi][iVibHi][iRotHi][0][0][0].Hi->Pop = H2_populations[iElecHi][iVibHi][iRotHi];
00255
00256
00257
00258
00259
00260 for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
00261 {
00262
00263
00264 long int nv = h2.nVib_hi[iElecLo];
00265 if( iElecLo==iElecHi )
00266 nv = iVibHi;
00267 for( iVibLo=0; iVibLo<=nv; ++iVibLo )
00268 {
00269 long nr = h2.nRot_hi[iElecLo][iVibLo];
00270 if( iElecLo==iElecHi && iVibHi==iVibLo )
00271 nr = iRotHi-1;
00272
00273 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
00274 {
00275 if( lgH2_line_exists[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] )
00276 {
00277
00278 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->PopOpc =
00279 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Lo->Pop -
00280 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Hi->Pop*
00281 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Lo->g /
00282 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Hi->g;
00283
00284
00285 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Coll.cool = 0.;
00286 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Coll.heat = 0.;
00287
00288
00289 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->xIntensity = 0.;
00290
00291 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->phots = 0.;
00292 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->ots = 0.;
00293 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->ColOvTot = 0.;
00294 }
00295 }
00296 }
00297 }
00298 }
00299 }
00300 }
00301 hmi.H2_photodissoc_BigH2_H2s = 0.;
00302 hmi.H2_photodissoc_BigH2_H2g = 0.;
00303 hmi.HeatH2Dish_BigH2 = 0.;
00304 hmi.HeatH2Dexc_BigH2 = 0.;
00305 hmi.deriv_HeatH2Dexc_BigH2 = 0.;
00306 hmi.H2_Solomon_dissoc_rate_BigH2_H2g = 0.;
00307 hmi.H2_Solomon_dissoc_rate_BigH2_H2s = 0.;
00308 hmi.H2_H2g_to_H2s_rate_BigH2 = 0.;
00309 return;
00310 }
00311
00312
00313 void mole_H2_LTE( void )
00314 {
00315
00316 static double TeUsedBoltz = -1.;
00317 double part_fun;
00318 long int iElec , iVib , iRot;
00319
00320 DEBUG_ENTRY( "mole_H2_LTE()" );
00321
00322
00323 if( ! fp_equal( phycon.te, TeUsedBoltz ) )
00324 {
00325 part_fun = 0.;
00326 TeUsedBoltz = phycon.te;
00327
00328 for( iElec=0; iElec<mole.n_h2_elec_states; ++iElec )
00329 {
00330 for( iVib=0; iVib<=h2.nVib_hi[iElec]; ++iVib )
00331 {
00332 for( iRot=h2.Jlowest[iElec]; iRot<=h2.nRot_hi[iElec][iVib]; ++iRot )
00333 {
00334 H2_Boltzmann[iElec][iVib][iRot] =
00335
00336
00337 sexp( energy_wn[iElec][iVib][iRot] / phycon.te_wn );
00338
00339 part_fun += H2_Boltzmann[iElec][iVib][iRot] * H2_stat[iElec][iVib][iRot];
00340 ASSERT( part_fun > 0 );
00341 }
00342 }
00343 }
00344
00345 for( iElec=0; iElec<mole.n_h2_elec_states; ++iElec )
00346 {
00347 for( iVib=0; iVib<=h2.nVib_hi[iElec]; ++iVib )
00348 {
00349 for( iRot=h2.Jlowest[iElec]; iRot<=h2.nRot_hi[iElec][iVib]; ++iRot )
00350 {
00351
00352
00353 H2_populations_LTE[iElec][iVib][iRot] =
00354 H2_Boltzmann[iElec][iVib][iRot] *
00355 H2_stat[iElec][iVib][iRot] / part_fun;
00356
00357
00358
00359 }
00360 }
00361 }
00362 if( mole.nH2_TRACE >= mole.nH2_trace_full )
00363 fprintf(ioQQQ,
00364 "mole_H2_LTE set H2_Boltzmann factors, T=%.2f, partition function is %.2f\n",
00365 phycon.te,
00366 part_fun);
00367 }
00368
00369 return;
00370 }
00371
00372
00373 void H2_init_coreload( void )
00374 {
00375
00376
00377
00378
00379
00380 long int nVib_hi_init[N_H2_ELEC] = {14 , 37 , 13 , 13, 9, 2 , 2};
00381
00382
00383
00384 long int Jlowest_init[N_H2_ELEC] = {0 , 0 , 1 , 1 , 0 , 1 , 1 };
00385
00386
00387
00388 long int nRot_hi_init[N_H2_ELEC][50]=
00389
00390 { {31, 30, 28, 27, 25,
00391 23, 22, 20, 18, 16,
00392 14, 12, 10, 7, 3 } ,
00393
00394 {25,25,25,25,25,25,25,25, 25,25,
00395 25,25,25,25,25,25,25,25, 25,25,
00396 25,25,25,25,25,25,25,25, 23,21,
00397 19,17,15,15,11,9,7, 7},
00398
00399 { 25, 25, 25, 25, 24, 23, 21, 19, 17, 14, 12, 10, 6, 2 },
00400
00401 { 25, 25, 25, 25, 24, 23, 21, 19, 17, 15, 13, 10, 7, 2 },
00402
00403 {19,17, 14, 12, 9, 8, 7, 7, 4, 1 },
00404
00405 {13, 10, 5},
00406
00407 {25 , 25 ,25 } }
00408 ;
00409
00410
00411
00412 long int iElec;
00413
00414 DEBUG_ENTRY( "H2_init_coreload()" );
00415
00416
00417
00418
00419
00420
00421 for( iElec=0; iElec<N_H2_ELEC; ++iElec )
00422 {
00423 int iVib;
00424 h2.nVib_hi[iElec] = nVib_hi_init[iElec];
00425 h2.Jlowest[iElec] = Jlowest_init[iElec];
00426 for( iVib=0; iVib<=h2.nVib_hi[iElec]; ++iVib )
00427 {
00428 h2.nRot_hi[iElec][iVib] = nRot_hi_init[iElec][iVib];
00429
00430
00431
00432
00433 }
00434 }
00435 strcpy( chH2ColliderLabels[0] , "H0" );
00436 strcpy( chH2ColliderLabels[1] , "He" );
00437 strcpy( chH2ColliderLabels[2] , "H2 o" );
00438 strcpy( chH2ColliderLabels[3] , "H2 p" );
00439 strcpy( chH2ColliderLabels[4] , "H+" );
00440
00441 return;
00442 }
00443
00444
00445 void H2_Init(void)
00446 {
00447
00448 DEBUG_ENTRY( "H2_Init()" );
00449
00450
00451
00452
00453 mole.n_h2_elec_states = N_H2_ELEC;
00454 h2.nCallH2_this_zone = 0;
00455 return;
00456 }
00457
00458
00459
00460 void H2_Reset( void )
00461 {
00462
00463 DEBUG_ENTRY( "H2_Reset()" );
00464
00465 if(mole.nH2_TRACE)
00466 fprintf(ioQQQ,
00467 "\n***************H2_Reset called, resetting nCallH2_this_iteration, zone %.2f iteration %li\n",
00468 fnzone,
00469 iteration );
00470
00471
00472
00473 nCallH2_this_iteration = 0;
00474
00475
00476
00477 h2.renorm_max = 1.;
00478 h2.renorm_min = 1.;
00479
00480
00481 nH2_pops = 0;
00482 nH2_zone = 0;
00483
00484 nzone_nlevel_set = 0;
00485
00486 nzoneAsEval = -1;
00487 iterationAsEval = -1;
00488
00489
00490 H2_SaveLine.zero();
00491
00492 return;
00493 }
00494
00495
00496
00497
00498
00499 void H2_Zero( void )
00500 {
00501
00502 DEBUG_ENTRY( "H2_Zero()" );
00503
00504
00505
00506
00507
00508
00509
00510 mole.H2_to_H_limit = 1e-8;
00511
00512 h2.lgH2ON = false;
00513
00514 mole.lgH2_LTE = false;
00515
00516
00517 nH2_pops = 0;
00518 nH2_zone = 0;
00519
00520 nzone_nlevel_set = 0;
00521
00522
00523
00524 h2.renorm_max = 1.;
00525 h2.renorm_min = 1.;
00526
00527 nCallH2_this_iteration = 0;
00528 h2.ortho_density = 0.;
00529 h2.para_density = 0.;
00530
00531 hmi.H2_Solomon_dissoc_rate_BigH2_H2s = 0.;
00532 hmi.H2_Solomon_dissoc_rate_BigH2_H2g = 0.;
00533
00534 hmi.H2_H2g_to_H2s_rate_BigH2 = 0.;
00535 hmi.H2_photodissoc_BigH2_H2s = 0.;
00536 hmi.H2_photodissoc_BigH2_H2g = 0.;
00537
00538 hmi.HeatH2Dexc_BigH2 = 0.;
00539
00540
00541 hmi.lgBigH2_evaluated = false;
00542
00543 hmi.lgH2_Thermal_BigH2 = true;
00544 hmi.lgH2_Chemistry_BigH2 = true;
00545
00546 if( !lgH2_READ_DATA )
00547 {
00548
00549
00550
00551
00552 mole.n_h2_elec_states = N_H2_ELEC;
00553 }
00554
00555
00556
00557
00558
00559 nXLevelsMatrix = 70;
00560
00561
00562 nzone_nlevel_set = -1;
00563
00564 nzoneAsEval = -1;
00565 iterationAsEval = -1;
00566 return;
00567 }