00001
00002
00003
00004
00005
00006
00007 #include "cddefines.h"
00008 #include "physconst.h"
00009 #include "dense.h"
00010 #include "called.h"
00011 #include "thermal.h"
00012 #include "gammas.h"
00013 #include "colden.h"
00014 #include "thermal.h"
00015 #include "secondaries.h"
00016 #include "h2.h"
00017 #include "mole.h"
00018 #include "radius.h"
00019 #include "doppvel.h"
00020 #include "rfield.h"
00021 #include "ionbal.h"
00022 #include "rt.h"
00023 #include "opacity.h"
00024 #include "iso.h"
00025 #include "conv.h"
00026 #include "phycon.h"
00027 #include "hmi.h"
00028 #include "cosmology.h"
00029
00030
00031 #if 0
00032 #if !defined(NDEBUG)
00033
00034
00035
00036
00037 #define AUDIT(a) { \
00038 double total, mtotal; \
00039 for( i=0;i<N_H_MOLEC;i++) { \
00040 total = 0.; \
00041 for( j=0;j<N_H_MOLEC;j++) { \
00042 total += c[i][j]*nprot[j]; \
00043 } \
00044 if( fabs(total) > 1e-6*fabs(c[i][i]*nprot[i])) { \
00045 fprintf(ioQQQ,"PROBLEM Subtotal1 %c %.2e\n",a,fabs(total)/fabs(c[i][i]*nprot[i])); \
00046 fprintf(ioQQQ,"Species %li Total %g Diag %g\n",i,total,c[i][i]*nprot[i]); \
00047 } \
00048 } \
00049 total = mtotal = 0.;for( j=0;j<N_H_MOLEC;j++) { total += bvec[j]*nprot[j]; mtotal += fabs(bvec[j]*nprot[j]); }\
00050 if( fabs(total) > 1e-30 && fabs(total) > 1e-10*rtot) { \
00051 fprintf(ioQQQ,"PROBLEM Subtotal2 %c %.2e\n",a,fabs(total)/mtotal); \
00052 fprintf(ioQQQ,"RHS Total %g cf %g\n",total,mtotal); \
00053 } else if( a == '.' && fabs(total) > 1e-7*mtotal) { \
00054 fprintf(ioQQQ,"PROBLEM zone %li Hmole RHS conservation error %.2e of %.2e\n",nzone,total,mtotal); \
00055 fprintf(ioQQQ,"(may be due to high rate equilibrium reactions)\n"); \
00056 } \
00057 }
00058 #else
00059 #define AUDIT
00060 #endif
00061
00062 #endif
00063
00064 void hmole( void )
00065 {
00066 int nFixup, i;
00067 double error;
00068 realnum oxy_error=0.;
00069 static realnum abund0=-BIGFLOAT , abund1=-BIGFLOAT;
00070 realnum save1=dense.xIonDense[ipOXYGEN][1],
00071 save0=dense.xIonDense[ipOXYGEN][0];
00072
00073 DEBUG_ENTRY( "hmole()" );
00074
00075
00076 nFixup = 0;
00077 error = 1.;
00078
00079
00080
00081 if( conv.nTotalIoniz==0 && iteration==0 )
00082 {
00083 realnum pdr_mole_h2[N_H_MOLEC] = {1.00E+00f,
00084 3.18E-05f,
00085 1.95E-11f,
00086 4.00E-08f,
00087 1.08E-14f,
00088 1.08E-20f,
00089 3.85E-07f,
00090 8.04E-14f};
00091
00092
00093 ASSERT( dense.xIonDense[ipHYDROGEN][0]>SMALLFLOAT );
00094
00095 ASSERT( ipMH==0&&
00096 ipMHp == 1&&
00097 ipMHm == 2&&
00098 ipMH2g == 3&&
00099 ipMH2p == 4&&
00100 ipMH3p == 5&&
00101 ipMH2s == 6&&
00102 ipMHeHp== 7&&
00103 N_H_MOLEC==8 );
00104
00105 for( i=0; i<N_H_MOLEC; ++i )
00106 {
00107 hmi.Hmolec[i] = dense.xIonDense[ipHYDROGEN][0]*pdr_mole_h2[i];
00108 }
00109 }
00110
00111
00112 hmole_reactions();
00113
00114
00115 # define BIGERROR 1e-4
00116
00117
00118 # define LIM_LOOP 20
00119
00120 for(i=0; i < LIM_LOOP && ((error > BIGERROR||nFixup || oxy_error>conv.EdenErrorAllowed));i++)
00121 {
00122 {
00123
00124
00125 enum {DEBUG_LOC=false};
00126
00127 if( DEBUG_LOC && (nzone>140) )
00128 {
00129 fprintf(ioQQQ,"DEBUG hmole in \t%.2f\t%.5e\t%.5e",
00130 fnzone,
00131 phycon.te,
00132 dense.eden);
00133 for( i=0; i<N_H_MOLEC; i++ )
00134 fprintf(ioQQQ,"\t%.2e", hmi.Hmolec[i] );
00135 fprintf(ioQQQ,"\n" );
00136 }
00137 }
00138
00139
00140 nFixup = 0;
00141
00142 if( !conv.lgSearch )
00143 {
00144 IonOxyge();
00145 }
00146
00147 save1 = dense.xIonDense[ipOXYGEN][1];
00148 save0 = dense.xIonDense[ipOXYGEN][0];
00149
00150
00151
00152
00153 if( nzone )
00154 {
00155 # define OLD 0.75f
00156 abund0 = abund0*OLD+dense.xIonDense[ipOXYGEN][0]*(1.f-OLD);
00157 abund1 = abund1*OLD+dense.xIonDense[ipOXYGEN][1]*(1.f-OLD);
00158 }
00159 else
00160 {
00161 abund0 = dense.xIonDense[ipOXYGEN][0];
00162 abund1 = dense.xIonDense[ipOXYGEN][1];
00163 }
00164 dense.xIonDense[ipOXYGEN][0] = abund0;
00165
00166 dense.xIonDense[ipOXYGEN][1] = abund1;
00167
00168 oxy_error = (realnum)fabs(save1-abund1)/SDIV(dense.gas_phase[ipOXYGEN]);
00169
00170
00171
00172 hmole_step(&nFixup, &error);
00173 dense.xIonDense[ipOXYGEN][0] = save0;
00174 dense.xIonDense[ipOXYGEN][1] = save1;
00175 {
00176
00177
00178 enum {DEBUG_LOC=false};
00179
00180 if( DEBUG_LOC )
00181 {
00182 fprintf(ioQQQ,"DEBUG hmole out\t%i\t%.2f\t%.5e\t%.5e",
00183 i,
00184 fnzone,
00185 phycon.te,
00186 dense.eden);
00187 fprintf(ioQQQ,
00188 "\terror\t%.4e\tH0\t%.4e\tH+\t%.4e\tsink\t%.4e\t%.4e\tsour\t%.4e\t%.4e\tion\t%.4e\trec\t%.4e",
00189 error,
00190 dense.xIonDense[ipHYDROGEN][0],
00191 dense.xIonDense[ipHYDROGEN][1],
00192 mole.sink[ipHYDROGEN][0],
00193 mole.sink[ipHYDROGEN][1],
00194 mole.source[ipHYDROGEN][0] ,
00195 mole.source[ipHYDROGEN][1] ,
00196 ionbal.RateIonizTot(ipHYDROGEN,0),
00197 ionbal.RateRecomTot[ipHYDROGEN][0]);
00198
00199
00200
00201 fprintf(ioQQQ,"\n" );
00202 }
00203 }
00204 }
00205
00206 if( (i == LIM_LOOP && error > BIGERROR) || nFixup != 0 )
00207 {
00208
00209
00210 conv.lgConvPops = false;
00211
00212 if( !conv.lgSearch && called.lgTalk )
00213 {
00214 fprintf(ioQQQ," PROBLEM hmole, zone %li: %d iters, %d bad; final error: %g lgSearch %i\n",
00215 nzone,
00216 i,
00217 nFixup,
00218 error,
00219 conv.lgSearch);
00220 }
00221 ConvFail( "pops" , "Hmole");
00222 }
00223
00224
00225 if( strcmp( dense.chDenseLaw, "CDEN" )==0 && !cosmology.lgDo )
00226
00227 ASSERT( fabs( dense.gas_phase[ipHYDROGEN] - dense.den0 )/
00228 dense.gas_phase[ipHYDROGEN]<1e-4 );
00229
00230
00231
00232
00233 dense.xMolecules[ipHYDROGEN] = 0.;
00234 for(i=0;i<N_H_MOLEC;i++)
00235 {
00236 dense.xMolecules[ipHYDROGEN] += hmi.Hmolec[i]*hmi.nProton[i];
00237 }
00238
00239 dense.xMolecules[ipHYDROGEN] -= (hmi.Hmolec[ipMH]+hmi.Hmolec[ipMHp]);
00240
00241
00242 for( i=0; i < mole.num_comole_calc; i++ )
00243 {
00244 dense.xMolecules[ipHYDROGEN] += COmole[i]->hevmol*COmole[i]->nElem[ipHYDROGEN];
00245 }
00246
00247
00248 return;
00249 }
00250
00251
00252 double hmirat(double te)
00253 {
00254 double hmirat_v;
00255
00256 DEBUG_ENTRY( "hmirat()" );
00257
00258
00259 if( te < 31.62 )
00260 {
00261 hmirat_v = 8.934e-18*phycon.sqrte*phycon.te003*phycon.te001*
00262 phycon.te001;
00263 }
00264 else if( te < 90. )
00265 {
00266 hmirat_v = 5.159e-18*phycon.sqrte*phycon.te10*phycon.te03*
00267 phycon.te03*phycon.te003*phycon.te001;
00268 }
00269 else if( te < 1200. )
00270 {
00271 hmirat_v = 2.042e-18*te/phycon.te10/phycon.te03;
00272 }
00273 else if( te < 3800. )
00274 {
00275 hmirat_v = 8.861e-18*phycon.te70/phycon.te03/phycon.te01*
00276 phycon.te003;
00277 }
00278 else if( te <= 1.4e4 )
00279 {
00280
00281 hmirat_v = 8.204e-17*phycon.sqrte/phycon.te10/phycon.te01*
00282 phycon.te003;
00283 }
00284 else
00285 {
00286
00287 hmirat_v = 5.44e-16*phycon.te20/phycon.te01;
00288 }
00289 return( hmirat_v );
00290 }
00291
00292
00293
00294 void hmole_init(void)
00295 {
00296
00297 DEBUG_ENTRY( "hmole_init()" );
00298
00299
00300 strcpy( hmi.chLab[ipMH], "H0 " );
00301 strcpy( hmi.chLab[ipMHp], "H+ " );
00302 strcpy( hmi.chLab[ipMHm], "H- " );
00303 strcpy( hmi.chLab[ipMH2g], "H2g " );
00304 strcpy( hmi.chLab[ipMH2p], "H2+ " );
00305 strcpy( hmi.chLab[ipMH3p], "H3+ " );
00306 strcpy( hmi.chLab[ipMH2s], "H2* " );
00307 strcpy( hmi.chLab[ipMHeHp], "HeH+" );
00308
00309
00310 hmi.rheph2hpheh = 0.;
00311 return;
00312
00313 }
00314
00315
00316 void hmole_reactions( void )
00317 {
00318 static double teused=-1;
00319 double exph2,
00320 exph2s,
00321 exphp,
00322 ex3hp;
00323 long i;
00324 double h2esc,
00325 th2,
00326 cr_H2s ,
00327 cr_H2dis,
00328 cr_H2dis_ELWERT_H2g,
00329 cr_H2dis_ELWERT_H2s;
00330
00331 DEBUG_ENTRY( "hmole_reactions()" );
00332
00333
00334
00335 if( ! fp_equal( phycon.te, teused ) )
00336 {
00337 teused = phycon.te;
00338
00339
00340
00341
00342 hmi.exphmi = sexp(8.745e3/phycon.te);
00343 if( hmi.exphmi > 0. )
00344 {
00345
00346 hmi.rel_pop_LTE_Hmin = SAHA/(phycon.te32*hmi.exphmi)*(1./(2.*2.));
00347 }
00348 else
00349 {
00350 hmi.rel_pop_LTE_Hmin = 0.;
00351 }
00352
00353
00354
00355
00356 if( h2.lgH2ON && hmi.lgBigH2_evaluated && hmi.lgH2_Chemistry_BigH2 )
00357 {
00358
00359 hmi.rel_pop_LTE_H2g = hmi.H2g_LTE_bigH2;
00360 hmi.rel_pop_LTE_H2s = hmi.H2s_LTE_bigH2;
00361 # if 0
00362 exph2 = sexp((5.195e4-hmi.H2_BigH2_H2g_av*T1CM)/phycon.te);
00363
00364 if( exph2 > 0. )
00365 hmi.rel_pop_LTE_H2g = SAHA/SDIV(phycon.te32*exph2)*(1./(2.*2.))*3.634e-5;
00366 else
00367 hmi.rel_pop_LTE_H2g = 0.;
00368
00369
00370
00371
00372
00373 exph2s = sexp((5.195e4-hmi.H2_BigH2_H2s_av * T1CM)/phycon.te);
00374
00375 if( exph2s > 0. )
00376 hmi.rel_pop_LTE_H2s = SAHA/SDIV(phycon.te32*exph2s)*(1./(2.*2.))*3.634e-5;
00377 else
00378 hmi.rel_pop_LTE_H2s = 0.;
00379
00380
00381
00382
00383 # endif
00384 }
00385 else
00386 {
00387
00388
00389 exph2 = sexp((5.195e4)/phycon.te);
00390
00391
00392 if( exph2 > 0. )
00393 {
00394
00395 hmi.rel_pop_LTE_H2g = SAHA/(phycon.te32*exph2)*(1./(2.*2.))*3.634e-5;
00396 }
00397 else
00398 {
00399 hmi.rel_pop_LTE_H2g = 0.;
00400 }
00401
00402
00403
00404
00405
00406 exph2s = sexp(2.178e4/phycon.te);
00407
00408 if( exph2s > 0. )
00409 {
00410
00411 hmi.rel_pop_LTE_H2s = SAHA/(phycon.te32*exph2s)*(1./(2.*2.))*3.634e-5;
00412 }
00413 else
00414 {
00415 hmi.rel_pop_LTE_H2s = 0.;
00416 }
00417 }
00418 {
00419
00420
00421
00422
00423 enum {DEBUG_LOC=false};
00424
00425 if( DEBUG_LOC && nzone>187&& iteration > 1)
00426 {
00427
00428 fprintf(ioQQQ,"ph2lteee\t%.2e\t%.1e\t%.1e\n",
00429 hmi.rel_pop_LTE_H2g,
00430 sexp(2.178e4/phycon.te),
00431 phycon.te);
00432 }
00433 }
00434
00435
00436
00437 exphp = sexp(3.072e4/phycon.te);
00438 if( exphp > 0. )
00439 {
00440
00441
00442 hmi.rel_pop_LTE_H2p = SAHA/(phycon.te32*exphp)*(4./(2.*1.))*3.634e-5;
00443 }
00444 else
00445 {
00446 hmi.rel_pop_LTE_H2p = 0.;
00447 }
00448
00449
00450
00451 ex3hp = sexp(1.882e4/phycon.te);
00452 if( ex3hp > 0. )
00453 {
00454
00455
00456 hmi.rel_pop_LTE_H3p = SAHA/(phycon.te32*ex3hp)*(4./(2.*1.))*3.634e-5;
00457 }
00458 else
00459 {
00460 hmi.rel_pop_LTE_H3p = 0.;
00461 }
00462 }
00463
00464
00465
00466 hmi.hminus_rad_attach = hmirat(phycon.te);
00467
00468 hmi.hmicol = hmi.hminus_rad_attach*EN1RYD*phycon.te*1.15e-5;
00469
00470
00471
00472
00473
00474 hmi.hminus_rad_attach *= dense.eden;
00475 hmi.hmicol *= dense.eden*hmi.Hmolec[ipMH];
00476
00477
00478
00479
00480
00481
00484
00485
00486
00487 t_phoHeat photoHeat;
00488
00489 hmi.HMinus_photo_rate = GammaBn( hmi.iphmin-1 , iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0] , opac.iphmop ,
00490 0.055502 , &hmi.HMinus_induc_rec_rate , &hmi.HMinus_induc_rec_cooling, &photoHeat );
00491
00492
00493 hmi.HMinus_photo_heat = photoHeat.HeatNet;
00494
00495 {
00496
00497
00498 enum {DEBUG_LOC=false};
00499
00500 if( DEBUG_LOC)
00501 {
00502 fprintf(ioQQQ,"hminphoto\t%li\t%li\t%.2e\n", nzone, conv.nPres2Ioniz , hmi.HMinus_photo_rate );
00503 }
00504 }
00505
00506
00507 hmi.HMinus_induc_rec_rate *= hmi.rel_pop_LTE_Hmin*dense.eden;
00508
00509
00511 hmi.HMinus_induc_rec_cooling *= hmi.rel_pop_LTE_Hmin*dense.eden*hmi.Hmolec[ipMH];
00512
00513 {
00514
00515
00516 enum {DEBUG_LOC=false};
00517
00518 if( DEBUG_LOC && nzone>400)
00519 {
00520 fprintf(ioQQQ,"hmoledebugg %.2f ",fnzone);
00521 GammaPrt(
00522 hmi.iphmin-1 , iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0] , opac.iphmop ,
00523
00524 ioQQQ,
00525
00526 hmi.HMinus_photo_rate,
00527
00528 hmi.HMinus_photo_rate*0.05);
00529 }
00530 }
00531
00532
00533
00534
00535
00536
00537
00538
00539
00540
00541
00542 hmi.UV_Cont_rel2_Habing_TH85_depth = 0.;
00543
00544
00545
00546
00547
00548 for( i=rfield.ipG0_TH85_lo; i < rfield.ipG0_TH85_hi; ++i )
00549 {
00550 hmi.UV_Cont_rel2_Habing_TH85_depth += (rfield.flux[0][i-1] + rfield.ConInterOut[i-1]+
00551 rfield.outlin[0][i-1]+ rfield.outlin_noplot[i-1])*rfield.anu[i-1];
00552 }
00553
00554
00555
00556
00557
00558 hmi.UV_Cont_rel2_Habing_TH85_depth = (realnum)(hmi.UV_Cont_rel2_Habing_TH85_depth*EN1RYD/1.6e-3);
00559
00560 if( nzone<=1 )
00561 {
00562 hmi.UV_Cont_rel2_Habing_TH85_face = hmi.UV_Cont_rel2_Habing_TH85_depth;
00563 }
00564
00565
00566
00567 hmi.UV_Cont_rel2_Habing_spec_depth = 0.;
00568 for( i=rfield.ipG0_spec_lo; i < rfield.ipG0_spec_hi; ++i )
00569 {
00570 hmi.UV_Cont_rel2_Habing_spec_depth += (rfield.flux[0][i-1] + rfield.ConInterOut[i-1]+
00571 rfield.outlin[0][i-1]+ rfield.outlin_noplot[i-1])*rfield.anu[i-1];
00572 }
00573 hmi.UV_Cont_rel2_Habing_spec_depth = (realnum)(hmi.UV_Cont_rel2_Habing_spec_depth*EN1RYD/1.6e-3);
00574
00575
00576
00577 if( hmi.UV_Cont_rel2_Draine_DB96_face ==0 )
00578 {
00579
00580 for( i=rfield.ipG0_DB96_lo; i < rfield.ipG0_DB96_hi; ++i )
00581 {
00582 hmi.UV_Cont_rel2_Draine_DB96_face += (rfield.flux[0][i-1] + rfield.ConInterOut[i-1]+
00583 rfield.outlin[0][i-1]+ rfield.outlin_noplot[i-1]);
00584 }
00585
00586
00587
00588
00589
00590
00591
00592
00593 hmi.UV_Cont_rel2_Draine_DB96_face = hmi.UV_Cont_rel2_Draine_DB96_face/(1.232e7f * 1.71f);
00594 }
00595
00596
00597
00598
00599 hmi.H2Opacity = (realnum)(1.2e-14*(1e5/GetDopplerWidth(2.f*dense.AtomicWeight[ipHYDROGEN])));
00600
00601 th2 = (colden.colden[ipCOL_H2g]+ colden.colden[ipCOL_H2s])*hmi.H2Opacity;
00602
00603
00604 h2esc = esc_PRD_1side(th2,1e-4);
00605
00606
00607
00608
00609
00610
00611
00612 hmi.H2_Solomon_dissoc_rate_TH85_H2g = 3.4e-11 * hmi.UV_Cont_rel2_Habing_TH85_depth * h2esc;
00613 hmi.H2_Solomon_dissoc_rate_TH85_H2s = 3.4e-11 * hmi.UV_Cont_rel2_Habing_TH85_depth * h2esc;
00614 hmi.H2_H2g_to_H2s_rate_TH85 = hmi.H2_Solomon_dissoc_rate_TH85_H2g*9.;
00615
00616
00617 hmi.H2_Solomon_dissoc_rate_BHT90_H2g = 3.4e-11 * hmi.UV_Cont_rel2_Habing_TH85_depth * h2esc;
00618 hmi.H2_Solomon_dissoc_rate_BHT90_H2s = 3.4e-11 * hmi.UV_Cont_rel2_Habing_TH85_depth * h2esc;
00619 hmi.H2_H2g_to_H2s_rate_BHT90 = hmi.H2_Solomon_dissoc_rate_BHT90_H2g*9.;
00620
00621 {
00622
00623
00624
00625
00626 double x = (colden.colden[ipCOL_H2g]+colden.colden[ipCOL_H2s]) / 5e14;
00627 double sqrtx = sqrt(1.+x);
00628
00629 double b5 = GetDopplerWidth(2.f*dense.AtomicWeight[ipHYDROGEN])/1e5;
00630
00631 double fshield = 0.965/POW2(1.+x/b5) + 0.035/sqrtx *
00632 sexp(8.5e-4*sqrtx);
00633
00634
00635
00636
00637
00638
00639 hmi.UV_Cont_rel2_Draine_DB96_depth = hmi.UV_Cont_rel2_Draine_DB96_face *
00640 (realnum)(sexp( opac.TauAbsFace[rfield.ip1000A-1] )/radius.r1r0sq);
00641
00642
00643
00644
00645 if( co.lgUMISTrates )
00646 {
00647
00648 hmi.H2_Solomon_dissoc_rate_BD96_H2g = 4.6e-11 * hmi.UV_Cont_rel2_Draine_DB96_depth * fshield;
00649 hmi.H2_Solomon_dissoc_rate_BD96_H2s = 4.6e-11 * hmi.UV_Cont_rel2_Draine_DB96_depth * fshield;
00650 }
00651 else
00652 {
00653
00654 hmi.H2_Solomon_dissoc_rate_BD96_H2g = 5.18e-11* (hmi.UV_Cont_rel2_Habing_TH85_face/1.66f)
00655 *sexp(3.02*rfield.extin_mag_V_point)* fshield;
00656 hmi.H2_Solomon_dissoc_rate_BD96_H2s = 5.18e-11* (hmi.UV_Cont_rel2_Habing_TH85_face/1.66f)
00657 *sexp(3.02*rfield.extin_mag_V_point)* fshield;
00658 }
00659
00660
00661
00662
00663
00664 hmi.H2_H2g_to_H2s_rate_BD96 = 5.67* hmi.H2_Solomon_dissoc_rate_BD96_H2g;
00665 }
00666
00667
00668 if( hmi.UV_Cont_rel2_Draine_DB96_face > SMALLFLOAT )
00669 {
00670
00671
00672
00673
00674
00675 double b5 = GetDopplerWidth(2.f*dense.AtomicWeight[ipHYDROGEN])/1e5;
00676
00677
00678
00679
00680
00681
00682
00683 double x_H2g, x_H2s,
00684 fshield_H2g, fshield_H2s,
00685 f_H2s;
00686 static double a_H2g, a_H2s,
00687 e1_H2g, e1_H2s,
00688 e2_H2g,
00689 b_H2g,
00690 sl_H2g, sl_H2s,
00691 k_f_H2s,
00692 k_H2g_to_H2s,
00693 log_G0_face = -1;
00694
00695
00696
00697
00698
00699
00700
00701
00702 {
00703 if(hmi.UV_Cont_rel2_Draine_DB96_face <= 1.)
00704 {
00705 log_G0_face = 0.;
00706 }
00707 else if(hmi.UV_Cont_rel2_Draine_DB96_face >= 1e7)
00708 {
00709 log_G0_face = 7.;
00710 }
00711 else
00712 {
00713 log_G0_face = log10(hmi.UV_Cont_rel2_Draine_DB96_face);
00714 }
00715
00716
00717
00718
00719 a_H2g = 0.06 * log_G0_face + 1.32;
00720 a_H2s = 0.375 * log_G0_face + 2.125;
00721
00722 e1_H2g = -0.05 * log_G0_face + 2.25;
00723 e1_H2s = -0.125 * log_G0_face + 2.625;
00724
00725 e2_H2g = -0.005 * log_G0_face + 0.625;
00726
00727 b_H2g = -4.0e-3 * log_G0_face + 3.2e-2;
00728
00729
00730 sl_H2g = 4.0e14;
00731 sl_H2s = 9.0e15;
00732
00733
00734 k_f_H2s = MAX2(0.1,2.375 * log_G0_face - 1.875 );
00735
00736
00737 k_H2g_to_H2s = MAX2(1.,-1.75 * log_G0_face + 11.25);
00738
00739
00740
00741
00742 }
00743
00744
00745 f_H2s = k_f_H2s * pow((double)hmi.UV_Cont_rel2_Draine_DB96_depth, 0.2 );
00746
00747
00748 x_H2g = (colden.colden[ipCOL_H2g]) / sl_H2g;
00749 x_H2s = (colden.colden[ipCOL_H2s]) / sl_H2s;
00750
00751
00752 fshield_H2g = 0.965/pow(1.+x_H2g/b5,e1_H2g) + b_H2g/pow(1.+x_H2g/b5,e2_H2g);
00753 fshield_H2s = 0.965/pow(1.+x_H2s/b5,e1_H2s);
00754
00755
00756 hmi.H2_Solomon_dissoc_rate_ELWERT_H2g = a_H2g * 4.6e-11 * fshield_H2g * hmi.UV_Cont_rel2_Draine_DB96_depth;
00757 hmi.H2_Solomon_dissoc_rate_ELWERT_H2s = a_H2s * 4.6e-11 * fshield_H2s * (hmi.UV_Cont_rel2_Draine_DB96_depth + f_H2s);
00758
00759
00760 hmi.H2_H2g_to_H2s_rate_ELWERT = k_H2g_to_H2s * hmi.H2_Solomon_dissoc_rate_ELWERT_H2g;
00761
00762
00763
00764
00765 hmi.H2_photodissoc_ELWERT_H2s = hmi.UV_Cont_rel2_Draine_DB96_depth*1e-11;
00766 hmi.H2_photodissoc_ELWERT_H2g = hmi.H2_photodissoc_ELWERT_H2s * 1.0e-10;
00767 }
00768 else
00769 {
00770 hmi.H2_Solomon_dissoc_rate_ELWERT_H2g = 0.;
00771 hmi.H2_Solomon_dissoc_rate_ELWERT_H2s = 0.;
00772 hmi.H2_photodissoc_ELWERT_H2s = 0.;
00773 hmi.H2_photodissoc_ELWERT_H2g = 0.;
00774 }
00775
00776
00777 hmi.H2_photodissoc_TH85 = hmi.UV_Cont_rel2_Habing_TH85_depth*1e-11;
00778
00779
00780 hmi.H2_photodissoc_BHT90 = hmi.UV_Cont_rel2_Habing_TH85_depth*1e-11;
00781
00782
00783
00784
00785
00786
00787 cr_H2s = secondaries.x12tot*0.9 / 2. * hmi.lgLeidenCRHack;
00788
00789
00790 cr_H2dis = secondaries.x12tot*0.1 / 2. * hmi.lgLeidenCRHack;
00791
00792
00793
00794
00795
00796 cr_H2dis_ELWERT_H2g = secondaries.x12tot*5e-8 * hmi.lgLeidenCRHack;
00797 cr_H2dis_ELWERT_H2s = secondaries.x12tot*4e-2 * hmi.lgLeidenCRHack;
00798
00799
00800
00801
00802
00803
00804
00805
00806
00807
00808
00809 if( h2.lgH2ON && hmi.lgBigH2_evaluated && hmi.lgH2_Chemistry_BigH2 )
00810 {
00811
00812
00813
00814
00815
00816
00817 hmi.H2_Solomon_dissoc_rate_used_H2g = hmi.H2_Solomon_dissoc_rate_BigH2_H2g;
00818 if( hmi.H2_Solomon_dissoc_rate_used_H2g <= 0. )
00819 {
00820
00821 hmi.H2_Solomon_dissoc_rate_used_H2g = hmi.H2_Solomon_dissoc_rate_TH85_H2g;
00822 }
00823
00824 hmi.H2_Solomon_dissoc_rate_used_H2s = hmi.H2_Solomon_dissoc_rate_BigH2_H2s;
00825 if( hmi.H2_Solomon_dissoc_rate_used_H2s <= 0. )
00826 hmi.H2_Solomon_dissoc_rate_used_H2s = hmi.H2_Solomon_dissoc_rate_TH85_H2g;
00827
00828
00829 hmi.H2_H2g_to_H2s_rate_used = hmi.H2_H2g_to_H2s_rate_BigH2;
00830 if( hmi.H2_H2g_to_H2s_rate_used <= 0. )
00831 hmi.H2_H2g_to_H2s_rate_used = hmi.H2_H2g_to_H2s_rate_TH85;
00832
00833
00834
00835
00836
00837 hmi.H2_photodissoc_used_H2s = hmi.H2_photodissoc_BigH2_H2s;
00838 if( hmi.H2_photodissoc_used_H2s <= 0. )
00839 hmi.H2_photodissoc_used_H2s = hmi.H2_photodissoc_TH85;
00840
00841
00842 hmi.H2_photodissoc_used_H2g = hmi.H2_photodissoc_BigH2_H2g;
00843 if( hmi.H2_photodissoc_used_H2g <= 0. )
00844 hmi.H2_photodissoc_used_H2g = hmi.H2_photodissoc_TH85*1.0e-10f;
00845 }
00846 else if( hmi.chH2_small_model_type == 'T' )
00847 {
00848
00849
00850 hmi.H2_Solomon_dissoc_rate_used_H2g = hmi.H2_Solomon_dissoc_rate_TH85_H2g + cr_H2dis;
00851
00852 hmi.H2_Solomon_dissoc_rate_used_H2s = hmi.H2_Solomon_dissoc_rate_TH85_H2s + cr_H2dis;
00853 hmi.H2_H2g_to_H2s_rate_used = hmi.H2_H2g_to_H2s_rate_TH85 + cr_H2s;
00854
00855
00856
00857 hmi.H2_photodissoc_used_H2s = hmi.H2_photodissoc_TH85;
00858
00859
00860 hmi.H2_photodissoc_used_H2g = hmi.H2_photodissoc_TH85*1.0e-10f;
00861 }
00862
00863 else if( hmi.chH2_small_model_type == 'H' )
00864 {
00865
00866
00867
00868 hmi.H2_Solomon_dissoc_rate_used_H2g = hmi.H2_Solomon_dissoc_rate_BHT90_H2g + cr_H2dis;
00869
00870 hmi.H2_Solomon_dissoc_rate_used_H2s = hmi.H2_Solomon_dissoc_rate_BHT90_H2s + cr_H2dis;
00871 hmi.H2_H2g_to_H2s_rate_used = hmi.H2_H2g_to_H2s_rate_BHT90 + cr_H2s;
00872
00873
00874
00875 hmi.H2_photodissoc_used_H2s = hmi.H2_photodissoc_BHT90;
00876
00877
00878 hmi.H2_photodissoc_used_H2g = hmi.H2_photodissoc_BHT90*1.0e-10f;
00879 }
00880
00881 else if( hmi.chH2_small_model_type == 'B' )
00882 {
00883
00884
00885 hmi.H2_Solomon_dissoc_rate_used_H2g = hmi.H2_Solomon_dissoc_rate_BD96_H2g + cr_H2dis;
00886
00887 hmi.H2_Solomon_dissoc_rate_used_H2s = hmi.H2_Solomon_dissoc_rate_BD96_H2s + cr_H2dis;
00888
00889 hmi.H2_H2g_to_H2s_rate_used = hmi.H2_H2g_to_H2s_rate_BD96 + cr_H2s;
00890
00891
00892
00893
00894 hmi.H2_photodissoc_used_H2s = hmi.H2_photodissoc_TH85;
00895
00896
00897 hmi.H2_photodissoc_used_H2g = hmi.H2_photodissoc_TH85*1.0e-10f;
00898 }
00899 else if( hmi.chH2_small_model_type == 'E' )
00900 {
00901
00902
00903 hmi.H2_Solomon_dissoc_rate_used_H2g = hmi.H2_Solomon_dissoc_rate_ELWERT_H2g + cr_H2dis_ELWERT_H2g;
00904 hmi.H2_Solomon_dissoc_rate_used_H2s = hmi.H2_Solomon_dissoc_rate_ELWERT_H2s + cr_H2dis_ELWERT_H2s;
00905 hmi.H2_H2g_to_H2s_rate_used = hmi.H2_H2g_to_H2s_rate_ELWERT + cr_H2s;
00906
00907
00908
00909
00910 hmi.H2_photodissoc_used_H2s = hmi.H2_photodissoc_ELWERT_H2s;
00911 hmi.H2_photodissoc_used_H2g = hmi.H2_photodissoc_ELWERT_H2g;
00912 }
00913 else
00914 TotalInsanity();
00915
00916 {
00917
00918 enum {DEBUG_LOC=false};
00919
00920 if( DEBUG_LOC && h2.lgH2ON )
00921 {
00922 fprintf(ioQQQ," Solomon H2 dest rates: TH85 %.2e BD96 %.2e Big %.2e excit rates: TH85 %.2e Big %.2e\n",
00923 hmi.H2_Solomon_dissoc_rate_TH85_H2g,
00924 hmi.H2_Solomon_dissoc_rate_BD96_H2g,
00925 hmi.H2_Solomon_dissoc_rate_BigH2_H2g ,
00926 hmi.H2_H2g_to_H2s_rate_TH85 ,
00927 hmi.H2_H2g_to_H2s_rate_BigH2);
00928 }
00929 }
00930 return;
00931 }
00932