00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #include "cddefines.h"
00013 #include "phycon.h"
00014 #include "mole.h"
00015 #include "hmi.h"
00016 #include "taulines.h"
00017 #include "h2.h"
00018 #include "h2_priv.h"
00019
00020
00021
00022 void H2_Solomon_rate( void )
00023 {
00024
00025 long int iElecLo , iElecHi , iVibLo , iVibHi , iRotLo , iRotHi,
00026 iplo , ip;
00027
00028 DEBUG_ENTRY( "H2_Solomon_rate()" );
00029
00030
00031 iElecLo = 0;
00032
00033
00034
00035
00036
00037
00038
00039 hmi.H2_Solomon_dissoc_rate_BigH2_H2g = 0.;
00040 hmi.H2_Solomon_dissoc_rate_BigH2_H2s = 0.;
00041
00042
00043 hmi.H2_Solomon_elec_decay_H2g = 0.;
00044 hmi.H2_Solomon_elec_decay_H2s = 0.;
00045
00046
00047
00048
00049
00050
00051 for( iElecHi=1; iElecHi<mole.n_h2_elec_states; ++iElecHi )
00052 {
00053 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
00054 {
00055 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
00056 {
00057
00058
00059
00060 for( iplo=0; iplo < nEner_H2_ground; ++iplo )
00061 {
00062 ip = H2_ipX_ener_sort[iplo];
00063 iRotLo = ipRot_H2_energy_sort[ip];
00064 iVibLo = ipVib_H2_energy_sort[ip];
00065 if( lgH2_line_exists[iElecHi][iVibHi][iRotHi][0][iVibLo][iRotLo] )
00066 {
00067
00068
00069 double rate_up_cont =
00070 H2_populations[iElecLo][iVibLo][iRotLo]*
00071 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].pump*
00072 H2_dissprob[iElecHi][iVibHi][iRotHi] / H2_rad_rate_out[iElecHi][iVibHi][iRotHi];
00073
00074
00075 hmi.H2_Solomon_dissoc_rate_BigH2_H2g += rate_up_cont;
00076
00077
00078 hmi.H2_Solomon_elec_decay_H2g +=
00079 H2_populations[iElecHi][iVibHi][iRotHi]*
00080 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Aul*(
00081 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Pesc+
00082 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Pdest+
00083 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Pelec_esc);
00084 }
00085 }
00086
00087 for( iplo=nEner_H2_ground; iplo < nLevels_per_elec[0]; ++iplo )
00088 {
00089 ip = H2_ipX_ener_sort[iplo];
00090 iRotLo = ipRot_H2_energy_sort[ip];
00091 iVibLo = ipVib_H2_energy_sort[ip];
00092 if( lgH2_line_exists[iElecHi][iVibHi][iRotHi][0][iVibLo][iRotLo] )
00093 {
00094
00095
00096 double rate_up_cont =
00097 H2_populations[iElecLo][iVibLo][iRotLo]*
00098 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].pump*
00099 H2_dissprob[iElecHi][iVibHi][iRotHi] / H2_rad_rate_out[iElecHi][iVibHi][iRotHi];
00100
00101
00102 hmi.H2_Solomon_dissoc_rate_BigH2_H2s += rate_up_cont;
00103
00104
00105 hmi.H2_Solomon_elec_decay_H2s +=
00106 H2_populations[iElecHi][iVibHi][iRotHi]*
00107 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Aul*(
00108 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Pesc+
00109 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Pdest+
00110 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Pelec_esc);
00111 }
00112 }
00113 }
00114 }
00115 }
00116
00117
00118
00119 if( hmi.H2_total > SMALLFLOAT )
00120 {
00121 hmi.H2_Solomon_elec_decay_H2g /= SDIV( H2_sum_excit_elec_den );
00122 hmi.H2_Solomon_elec_decay_H2s /= SDIV( H2_sum_excit_elec_den );
00123
00124
00125 hmi.H2_Solomon_dissoc_rate_BigH2_H2s = hmi.H2_Solomon_dissoc_rate_BigH2_H2s / SDIV(H2_den_s);
00126
00127
00128 hmi.H2_Solomon_dissoc_rate_BigH2_H2g = hmi.H2_Solomon_dissoc_rate_BigH2_H2g / SDIV(H2_den_g);
00129
00130 }
00131 else
00132 {
00133 hmi.H2_Solomon_dissoc_rate_BigH2_H2s = 0;
00134 hmi.H2_Solomon_dissoc_rate_BigH2_H2g = 0;
00135
00136 }
00137
00138
00139
00140
00141
00142
00143
00144 hmi.H2_H2g_to_H2s_rate_BigH2 = 0.;
00145
00146
00147 DEBUG_EXIT( "H2_Solomon_rate()" );
00148
00149 return;
00150 }
00151
00152
00153 void H2_gs_rates( void )
00154 {
00155 long int ipLoX , ip , iRotLoX , iVibLoX ,
00156 iElecHi , iVibHi , ipOther , iRotOther,
00157 iRotHi , iVibOther;
00158
00159 DEBUG_ENTRY( "H2_gs_rates()" );
00160
00161
00162 hmi.H2_H2g_to_H2s_rate_BigH2 = 0.;
00163
00164
00165 for( ipLoX=0; ipLoX < nEner_H2_ground; ++ipLoX )
00166 {
00167 ip = H2_ipX_ener_sort[ipLoX];
00168 iRotLoX = ipRot_H2_energy_sort[ip];
00169 iVibLoX = ipVib_H2_energy_sort[ip];
00170
00171
00172 for( iElecHi=1; iElecHi<mole.n_h2_elec_states; ++iElecHi )
00173 {
00174 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
00175 {
00176 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
00177 {
00178 if( lgH2_line_exists[iElecHi][iVibHi][iRotHi][0][iVibLoX][iRotLoX] )
00179 {
00180
00181
00182 double rate_up_cont =
00183 H2_populations[0][iVibLoX][iRotLoX]*
00184 H2Lines[iElecHi][iVibHi][iRotHi][0][iVibLoX][iRotLoX].pump;
00185
00186 double decay_star = H2_rad_rate_out[iElecHi][iVibHi][iRotHi] - H2_dissprob[iElecHi][iVibHi][iRotHi];
00187
00188
00189 for( ipOther=0; ipOther < nEner_H2_ground; ++ipOther )
00190 {
00191 ip = H2_ipX_ener_sort[ipOther];
00192 iRotOther = ipRot_H2_energy_sort[ip];
00193 iVibOther = ipVib_H2_energy_sort[ip];
00194 if( lgH2_line_exists[iElecHi][iVibHi][iRotHi][0][iVibOther][iRotOther] )
00195 {
00196 decay_star -=
00197 H2Lines[iElecHi][iVibHi][iRotHi][0][iVibOther][iRotOther].Aul*(
00198 H2Lines[iElecHi][iVibHi][iRotHi][0][iVibOther][iRotOther].Pesc+
00199 H2Lines[iElecHi][iVibHi][iRotHi][0][iVibOther][iRotOther].Pdest+
00200 H2Lines[iElecHi][iVibHi][iRotHi][0][iVibOther][iRotOther].Pelec_esc);
00201 }
00202 }
00203
00204
00205 decay_star = MAX2(0., decay_star)/H2_rad_rate_out[iElecHi][iVibHi][iRotHi];
00206 hmi.H2_H2g_to_H2s_rate_BigH2 += rate_up_cont*decay_star;
00207
00208 }
00209 }
00210 }
00211 }
00212 }
00213
00214
00215 hmi.H2_H2g_to_H2s_rate_BigH2 /= SDIV( H2_den_g );
00216
00217 DEBUG_EXIT( "H2_gs_rates()" );
00218
00219 return;
00220 }
00221
00222
00223
00224 void H2_zero_pops_too_low( void )
00225 {
00226
00227 long int iElec, iElecHi, iElecLo, iVib , iVibLo , iVibHi ,
00228 iRot , iRotHi , iRotLo;
00229
00230 DEBUG_ENTRY( "H2_zero_pops_too_low()" );
00231
00232
00233 for( iElec=0; iElec<mole.n_h2_elec_states; ++iElec )
00234 {
00235 for( iVib=0; iVib<=h2.nVib_hi[iElec]; ++iVib )
00236 {
00237 for( iRot=h2.Jlowest[iElec]; iRot<=h2.nRot_hi[iElec][iVib]; ++iRot )
00238 {
00239
00240
00241 H2_old_populations[iElec][iVib][iRot] =
00242 (float)H2_populations_LTE[iElec][iVib][iRot]*hmi.H2_total;
00243 H2_populations[iElec][iVib][iRot] = H2_old_populations[iElec][iVib][iRot];
00244 }
00245 }
00246 }
00247
00248
00249
00250 for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi )
00251 {
00252 pops_per_elec[iElecHi] = 0.;
00253 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
00254 {
00255 pops_per_vib[iElecHi][iVibHi] = 0.;
00256 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
00257 {
00258 long int lim_elec_lo = 0;
00259
00260
00261
00262
00263
00264 for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
00265 {
00266
00267
00268 long int nv = h2.nVib_hi[iElecLo];
00269 if( iElecLo==iElecHi )
00270 nv = iVibHi;
00271 for( iVibLo=0; iVibLo<=nv; ++iVibLo )
00272 {
00273 long nr = h2.nRot_hi[iElecLo][iVibLo];
00274 if( iElecLo==iElecHi && iVibHi==iVibLo )
00275 nr = iRotHi-1;
00276
00277 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
00278 {
00279 if( lgH2_line_exists[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] )
00280 {
00281
00282
00283
00284 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].PopLo =
00285 H2_populations[iElecLo][iVibLo][iRotLo];
00286
00287
00288 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].PopHi =
00289 H2_populations[iElecHi][iVibHi][iRotHi];
00290
00291
00292 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].PopOpc =
00293 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].PopLo -
00294 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].PopHi*
00295 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].gLo /
00296 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].gHi;
00297
00298
00299 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].cool = 0.;
00300 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].heat = 0.;
00301
00302
00303 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].xIntensity = 0.;
00304
00305 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].phots = 0.;
00306 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].ots = 0.;
00307 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].ColOvTot = 0.;
00308 }
00309 }
00310 }
00311 }
00312 }
00313 }
00314 }
00315 hmi.H2_photodissoc_BigH2_H2s = 0.;
00316 hmi.H2_photodissoc_BigH2_H2g = 0.;
00317 hmi.HeatH2Dish_BigH2 = 0.;
00318 hmi.HeatH2Dexc_BigH2 = 0.;
00319 hmi.deriv_HeatH2Dexc_BigH2 = 0.;
00320 hmi.H2_Solomon_dissoc_rate_BigH2_H2g = 0.;
00321 hmi.H2_Solomon_dissoc_rate_BigH2_H2s = 0.;
00322 hmi.H2_H2g_to_H2s_rate_BigH2 = 0.;
00323
00324 DEBUG_EXIT( "H2_zero_pops_too_low()" );
00325
00326 return;
00327 }
00328
00329
00330 void mole_H2_LTE( void )
00331 {
00332
00333 static double TeUsedBoltz = -1.;
00334 double part_fun;
00335 long int iElec , iVib , iRot;
00336
00337 DEBUG_ENTRY( "mole_H2_LTE()" );
00338
00339
00340
00341 if( phycon.te != TeUsedBoltz )
00342
00343 {
00344 part_fun = 0.;
00345 TeUsedBoltz = phycon.te;
00346
00347 for( iElec=0; iElec<mole.n_h2_elec_states; ++iElec )
00348 {
00349 for( iVib=0; iVib<=h2.nVib_hi[iElec]; ++iVib )
00350 {
00351 for( iRot=h2.Jlowest[iElec]; iRot<=h2.nRot_hi[iElec][iVib]; ++iRot )
00352 {
00353 H2_Boltzmann[iElec][iVib][iRot] =
00354
00355
00356 sexp( energy_wn[iElec][iVib][iRot] / phycon.te_wn );
00357
00358 part_fun += H2_Boltzmann[iElec][iVib][iRot] * H2_stat[iElec][iVib][iRot];
00359 ASSERT( part_fun > 0 );
00360 }
00361 }
00362 }
00363
00364 for( iElec=0; iElec<mole.n_h2_elec_states; ++iElec )
00365 {
00366 for( iVib=0; iVib<=h2.nVib_hi[iElec]; ++iVib )
00367 {
00368 for( iRot=h2.Jlowest[iElec]; iRot<=h2.nRot_hi[iElec][iVib]; ++iRot )
00369 {
00370
00371
00372 H2_populations_LTE[iElec][iVib][iRot] =
00373 H2_Boltzmann[iElec][iVib][iRot] *
00374 H2_stat[iElec][iVib][iRot] / part_fun;
00375
00376
00377
00378 }
00379 }
00380 }
00381 if( mole.nH2_TRACE >= mole.nH2_trace_full )
00382 fprintf(ioQQQ,
00383 "mole_H2_LTE set H2_Boltzmann factors, T=%.2f, partition function is %.2f\n",
00384 phycon.te,
00385 part_fun);
00386 }
00387
00388
00389 DEBUG_EXIT( "mole_H2_LTE()" );
00390
00391 return;
00392
00393 }
00394
00395
00396 void H2_Init(void)
00397 {
00398
00399
00400
00401
00402
00403 long int nVib_hi_init[N_H2_ELEC] = {14 , 37 , 13 , 13, 9, 2 , 2};
00404
00405
00406
00407 long int Jlowest_init[N_H2_ELEC] = {0 , 0 , 1 , 1 , 0 , 1 , 1 };
00408
00409
00410
00411 long int nRot_hi_init[N_H2_ELEC][50]=
00412
00413 { {31, 30, 28, 27, 25,
00414 23, 22, 20, 18, 16,
00415 14, 12, 10, 7, 3 } ,
00416
00417 {25,25,25,25,25,25,25,25, 25,25,
00418 25,25,25,25,25,25,25,25, 25,25,
00419 25,25,25,25,25,25,25,25, 23,21,
00420 19,17,15,15,11,9,7, 7},
00421
00422 { 25, 25, 25, 25, 24, 23, 21, 19, 17, 14, 12, 10, 6, 2 },
00423
00424 { 25, 25, 25, 25, 24, 23, 21, 19, 17, 15, 13, 10, 7, 2 },
00425
00426 {19,17, 14, 12, 9, 8, 7, 7, 4, 1 },
00427
00428 {13, 10, 5},
00429
00430 {25 , 25 ,25 } }
00431 ;
00432
00433
00434 long int iElec;
00435
00436 DEBUG_ENTRY( "H2_Init()" );
00437
00438
00439
00440
00441 mole.n_h2_elec_states = N_H2_ELEC;
00442 h2.nCallH2_this_zone = 0;
00443
00444
00445
00446
00447
00448
00449 for( iElec=0; iElec<N_H2_ELEC; ++iElec )
00450 {
00451 int iVib;
00452 h2.nVib_hi[iElec] = nVib_hi_init[iElec];
00453 h2.Jlowest[iElec] = Jlowest_init[iElec];
00454 for( iVib=0; iVib<=h2.nVib_hi[iElec]; ++iVib )
00455 {
00456 h2.nRot_hi[iElec][iVib] = nRot_hi_init[iElec][iVib];
00457
00458
00459
00460
00461 }
00462 }
00463
00464 DEBUG_EXIT( "H2_Init()" );
00465
00466 return;
00467 }
00468
00469
00470
00471 void H2_Reset( void )
00472 {
00473
00474 DEBUG_ENTRY( "H2_Reset()" );
00475
00476 if(mole.nH2_TRACE)
00477 fprintf(ioQQQ,
00478 "\n***************H2_Reset called, resetting nCallH2_this_iteration, zone %.2f iteration %li\n",
00479 fnzone,
00480 iteration );
00481
00482
00483
00484 nCallH2_this_iteration = 0;
00485
00486
00487
00488 h2.renorm_max = 1.;
00489 h2.renorm_min = 1.;
00490
00491
00492 nH2_pops = 0;
00493 nH2_zone = 0;
00494
00495 nzone_nlevel_set = 0;
00496
00497 nzoneAsEval = -1;
00498 iterationAsEval = -1;
00499
00500 DEBUG_EXIT( "H2_Reset()" );
00501
00502 return;
00503
00504 }
00505
00506
00507
00508
00509
00510 void H2_Zero( void )
00511 {
00512
00513 DEBUG_ENTRY( "H2_Zero()" );
00514
00515
00516
00517
00518
00519
00520
00521 mole.H2_to_H_limit = 1e-8;
00522
00523
00524
00525 mole.nH2_trace_final = 1;
00526
00527 mole.nH2_trace_iterations = 2;
00528
00529 mole.nH2_trace_full = 3;
00530
00531 mole.nH2_trace_matrix = 4;
00532
00533 h2.lgH2ON = false;
00534 mole.nH2_TRACE = false;
00535
00536 mole.lgH2_LTE = false;
00537
00538
00539
00540
00541 h2.nElecLevelOutput = 1;
00542
00543
00544 nH2_pops = 0;
00545 nH2_zone = 0;
00546
00547 nzone_nlevel_set = 0;
00548
00549
00550 mole.lgH2_NOISE = false;
00551 mole.lgH2_NOISECOSMIC = false;
00552
00553
00554
00555
00556 mole.lgColl_gbar = true;
00557
00558
00559
00560 mole.lgColl_deexec_Calc = true;
00561
00562 mole.lgColl_dissoc_coll = true;
00563
00564
00565
00566 mole.lgH2_grain_deexcitation = false;
00567
00568
00569 mole.lgH2_ortho_para_coll_on = true;
00570
00571
00572
00573
00574 mole.lgH2_He_Meudon = true;
00575 mole.lgH2_He_Stancil = false;
00576
00577
00578
00579 h2.renorm_max = 1.;
00580 h2.renorm_min = 1.;
00581
00582 nCallH2_this_iteration = 0;
00583 h2.ortho_density = 0.;
00584 h2.para_density = 0.;
00585
00586 hmi.H2_Solomon_dissoc_rate_BigH2_H2s = 0.;
00587 hmi.H2_Solomon_dissoc_rate_BigH2_H2g = 0.;
00588
00589 hmi.H2_H2g_to_H2s_rate_BigH2 = 0.;
00590 hmi.H2_photodissoc_BigH2_H2s = 0.;
00591 hmi.H2_photodissoc_BigH2_H2g = 0.;
00592
00593 hmi.HeatH2Dexc_BigH2 = 0.;
00594
00595
00596 hmi.lgBigH2_evaluated = false;
00597
00598 hmi.lgH2_Thermal_BigH2 = true;
00599 hmi.lgH2_Chemistry_BigH2 = true;
00600
00601 if( !lgH2_READ_DATA )
00602 {
00603
00604
00605
00606
00607 mole.n_h2_elec_states = N_H2_ELEC;
00608 }
00609
00610
00611
00612
00613
00614 nXLevelsMatrix = 70;
00615
00616
00617 nzone_nlevel_set = -1;
00618
00619 nzoneAsEval = -1;
00620 iterationAsEval = -1;
00621
00622 DEBUG_EXIT( "H2_Zero()" );
00623
00624 return;
00625 }
00626
00627 #define SizeOf_line 300
00628 #define Num_fit_par 8
00629 #define nu 302
00630 #define nl 302
00631 static double fit_par[nu][nl][Num_fit_par];
00632 static int defn[nu][nl];
00633
00634
00635
00636
00637
00638 long int H2_He_coll_init(char FILE_NAME_IN[] )
00639 {
00640
00641 int i, j;
00642 long int magic;
00643 double par[Num_fit_par];
00644 char line[SizeOf_line];
00645 int h2_i, h2_f, he_i, he_f;
00646 char quality, space;
00647
00648
00649 double error;
00650
00651 FILE *ifp;
00652
00653 # if 0
00654 FILE *ioQQQ;
00655 ioQQQ = fopen(FILE_NAME_OUT, "w");
00656 if( ioQQQ == NULL )
00657 {
00658 printf("Can't open %s\n", FILE_NAME_OUT);
00659 exit(1);
00660 }
00661 # endif
00662
00663 #if 0
00664
00665 ifp = fopen(FILE_NAME_IN, "r");
00666 if( ifp == NULL )
00667 {
00668 printf("Can't open %s\n", FILE_NAME_IN);
00669 exit(1);
00670 }
00671 while( fgets(line, (int)sizeof(line), ifp) != NULL )
00672 {
00673 Num_lines++;
00674 }
00675 printf("The number of lines are %i.\n", Num_lines);
00676
00677 fclose(ifp);
00678 #endif
00679
00680
00681 for( i=0; i<nu; i++ )
00682 {
00683 for( j=0; j<nl; j++ )
00684 {
00685 defn[i][j] = 0;
00686 }
00687 }
00688
00689
00690 ifp = fopen(FILE_NAME_IN, "r");
00691 if( ifp == NULL )
00692 {
00693 printf("Can't open %s\n", FILE_NAME_IN);
00694 exit(1);
00695 }
00696
00697
00698 fgets(line, (int)sizeof(line), ifp);
00699 sscanf(line, "%li", &magic);
00700
00701
00702 while( fgets(line, (int)sizeof(line), ifp) != NULL )
00703 {
00704
00705 if( line[0]!='#' )
00706 {
00707 sscanf(line, "%i%i%i%i%c%c%c%c%lf%lf%lf%lf%lf%lf%lf%lf%lf", &h2_i, &h2_f, &he_i, &he_f,&space, &space, &space, &quality, &error, &par[0], &par[1], &par[2], &par[3], &par[4], &par[5], &par[6], &par[7]);
00708
00709
00710 defn[h2_i][h2_f] = 1;
00711 for( i=0; i<Num_fit_par; i++ )
00712
00713 fit_par[h2_i][h2_f][i] = par[i];
00714 }
00715 }
00716 fclose(ifp);
00717
00718 # if 0
00719 for( i=0; i<nu; i++ )
00720 {
00721 for( j=0; j<nl; j++ )
00722 {
00723 if( defn[i][j] == 1 )
00724 {
00725 fprintf(ioQQQ, "%i %i %e %e %e %e %e %e %e %e\n", i, j, fit_par[i][j][0], fit_par[i][j][1], fit_par[i][j][2], fit_par[i][j][3], fit_par[i][j][4], fit_par[i][j][5], fit_par[i][j][6], fit_par[i][j][7]);
00726 }
00727 }
00728 }
00729
00730 fclose(ofp);
00731 # endif
00732 return magic;
00733 }
00734
00735
00736
00737
00738
00739 double H2_He_coll(int init, int final, double temp)
00740 {
00741
00742 double k, t, b2, c2, f2, t2;
00743
00744
00745 if( temp<2 || temp > 1e8 )
00746 {
00747 k = -1;
00748 }
00749 else if( init <= final )
00750 {
00751 k = -1;
00752 }
00753
00754 else if( init < 0 || init >302 || final < 0 || final > 302 )
00755 {
00756 k = -1;
00757 }
00758
00759 else if( defn[init][final] == 0 )
00760 {
00761 k = -1;
00762 }
00763
00764 else if( defn[init][final] == 1 )
00765 {
00766 double sum1 , sum2;
00767
00768 t2 = temp/1e3;
00769
00770 t = t2+1;
00771 b2 = fit_par[init][final][1]/(fit_par[init][final][3]*t2+1);
00772 c2 = fit_par[init][final][2]/(t*t);
00773 {
00774
00775 enum {DEBUG_LOC=false};
00776
00777 if( DEBUG_LOC )
00778 {
00779 fprintf(ioQQQ,"bug H2 He coll\t%i %i %.3e %.3e %.3e \n",
00780 init,final,
00781 fit_par[init][final][5],
00782 fit_par[init][final][6]*t2+1,
00783 fit_par[init][final][7]
00784
00785 );
00786 }
00787 }
00788
00789 sum1 = fit_par[init][final][7] * log10( fit_par[init][final][6]*t2+1. );
00790
00791 if( fabs(sum1)< 38. )
00792 {
00793
00794
00795
00796
00797 f2 = fit_par[init][final][5]/pow(10. , sum1 );
00798 }
00799 else
00800 f2= 0.;
00801 {
00802
00803 enum {DEBUG_LOC=false };
00804
00805 if( DEBUG_LOC )
00806 {
00807 fprintf(ioQQQ,"bug H2 He coll\t%i %i %.3e %.3e %.3e %.3e %.3e sum %.3e %.3e \n",
00808 init,final,
00809 fit_par[init][final][0],
00810 b2,
00811 c2,
00812 fit_par[init][final][4],
00813 f2 ,
00814 fit_par[init][final][0]+b2+c2,
00815 fit_par[init][final][4]+f2);
00816 }
00817 }
00818
00819
00820 sum1 = fit_par[init][final][0]+b2+c2;
00821 sum2 = fit_par[init][final][4]+f2;
00822 k = 0.;
00823 if( fabs(sum1) < 38. )
00824 k += pow(10., sum1 );
00825 if( fabs(sum2) < 38. )
00826 k += pow(10., sum2 );
00827 {
00828
00829 enum {DEBUG_LOC=false};
00830
00831 if( DEBUG_LOC )
00832 {
00833 fprintf(ioQQQ,"DEBUG H2 He rate %.3e \n",
00834 k);
00835 }
00836 }
00837 }
00838
00839 else
00840 {
00841 k = 99;
00842 }
00843
00844 return k;
00845 }
00846