00001 #include "cdstd.h"
00002 #include <stdlib.h>
00003 #include "cddefines.h"
00004 #include "colden.h"
00005 #include "physconst.h"
00006 #include "mole.h"
00007 #include "version.h"
00008 #include "mole_priv.h"
00009 #include "hmi.h"
00010 #include "conv.h"
00011 #include "dense.h"
00012 #include "thermal.h"
00013 #include "gammas.h"
00014 #include "grainvar.h"
00015 #include "ionbal.h"
00016 #include "opacity.h"
00017 #include "path.h"
00018 #include "radius.h"
00019 #include "thermal.h"
00020 #include "atmdat.h"
00021 #include "taulines.h"
00022 #include "trace.h"
00023 #include "deuterium.h"
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047 enum {UDFA = 0};
00048
00049 STATIC void newreact(const char label[], const char fun[], double a, double b, double c);
00050 STATIC long parse_reaction( count_ptr<mole_reaction>& rate, const char label[] );
00051 STATIC string canonicalize_reaction_label( const char label[] );
00052 STATIC void canonicalize_reaction( count_ptr<mole_reaction>& rate );
00053 STATIC void register_reaction_vectors( count_ptr<mole_reaction> rate );
00054
00055 STATIC void mole_check_reverse_reactions(void);
00056 STATIC double mole_get_equilibrium_constant( const char buf[] );
00057 STATIC double mole_partition_function( const molecule* const sp);
00058 STATIC void mole_generate_isotopologue_reactions( string atom_old, string atom_new );
00059
00060
00061 STATIC void plot_sparsity(void);
00062
00063
00064 STATIC bool lgReactBalance(const count_ptr<mole_reaction> &rate);
00065
00066
00067 STATIC void read_data(const char file[], void (*parse)(char *s));
00068 STATIC void parse_base(char *s);
00069 STATIC void parse_udfa(char *s);
00070 STATIC void compare_udfa(const count_ptr<mole_reaction> &rate);
00071
00072
00073
00074
00075 static realnum albedo = 0.5;
00076
00077
00078 namespace {
00079
00080 typedef count_ptr<mole_reaction> ratefunc;
00081 template<class T> void newfunc()
00082 {
00083 ratefunc fun = count_ptr<mole_reaction>(new T);
00084 ASSERT( mole_priv::functab.find(fun->name()) == mole_priv::functab.end() );
00085 mole_priv::functab[fun->name()] = fun;
00086 }
00087 ratefunc findfunc(const char name[]);
00088
00089
00090 class mole_reaction_hmrate;
00091 double hmrate(const mole_reaction *rate);
00092 class mole_reaction_th85rate;
00093 double th85rate(const mole_reaction *rate);
00094 class mole_reaction_crnurate_noalbedo;
00095 double crnurate_noalbedo(const mole_reaction *rate);
00096
00097
00098 class mole_reaction_grn_react;
00099 double grn_react(const mole_reaction *rate);
00100 class mole_reaction_h2_collh_excit;
00101 double h2_collh_excit(const mole_reaction *rate);
00102 class mole_reaction_h2_collh2_excit;
00103 double h2_collh2_excit(const mole_reaction *rate);
00104 class mole_reaction_h2_collh_deexc;
00105 double h2_collh_deexc(const mole_reaction *rate);
00106 class mole_reaction_h2_collh2_deexc;
00107 double h2_collh2_deexc(const mole_reaction *rate);
00108 class mole_reaction_rh2g_dis_h;
00109 double rh2g_dis_h(const mole_reaction *rate);
00110 class mole_reaction_rh2s_dis_h;
00111 double rh2s_dis_h(const mole_reaction *rate);
00112 class mole_reaction_rh2g_dis_h2;
00113 double rh2g_dis_h2(const mole_reaction *rate);
00114 class mole_reaction_rh2s_dis_h2;
00115 double rh2s_dis_h2(const mole_reaction *rate);
00116 class mole_reaction_rh2s_dis_h2_nodeexcit;
00117 double rh2s_dis_h2_nodeexcit(const mole_reaction *rate);
00118 class mole_reaction_bh2g_dis_h;
00119 double bh2g_dis_h(const mole_reaction *rate);
00120 class mole_reaction_bh2s_dis_h;
00121 double bh2s_dis_h(const mole_reaction *rate);
00122 class mole_reaction_bh2g_dis_h2;
00123 double bh2g_dis_h2(const mole_reaction *rate);
00124 class mole_reaction_hneut;
00125 double hneut(const mole_reaction *rate);
00126 class mole_reaction_cionhm;
00127 double cionhm(const mole_reaction *rate);
00128 }
00129
00130
00131
00132
00133
00134
00135
00136 #include "phycon.h"
00137 #include "doppvel.h"
00138
00139 namespace
00140 {
00141 class mole_reaction_null : public mole_reaction
00142 {
00143 typedef mole_reaction_null T;
00144 public:
00145 virtual T* Create() const {return new T;}
00146 virtual const char* name() {return "null";}
00147 double rk() const
00148 {
00149 ASSERT( false );
00150 return 0.0;
00151 }
00152 };
00153 }
00154
00155
00156
00157
00158
00159
00160
00161
00162 namespace {
00163
00164
00165
00166 STATIC double noneq_offset(const mole_reaction *rate)
00167 {
00168
00169 int nreact, n;
00170 bool lgFact;
00171
00172 DEBUG_ENTRY("noneq_offset()");
00173
00174 lgFact = false;
00175 if(mole_global.lgNonEquilChem)
00176 {
00177 if(mole_global.lgNeutrals)
00178 {
00179 lgFact = true;
00180 }
00181 else
00182 {
00183 nreact = rate->nreactants;
00184 for(n=0;n<nreact;n++)
00185 {
00186 if (rate->reactants[n]->charge != 0)
00187 {
00188 lgFact = true;
00189 break;
00190 }
00191 }
00192 }
00193 }
00194
00195 if( lgFact )
00196 return 0.333f*POW2(DoppVel.TurbVel)/BOLTZMANN*rate->reduced_mass;
00197 else
00198 return 0.;
00199 }
00200 double hmrate(const mole_reaction *rate)
00201 {
00202 double te;
00203
00204 DEBUG_ENTRY("hmrate()");
00205
00206 te = phycon.te+noneq_offset(rate);
00207
00208 if( rate->c < 0. )
00209 ASSERT( -rate->c/te < 10. );
00210
00211 return pow(te/300.,rate->b)*exp(-rate->c/te);
00212 }
00213
00214 class mole_reaction_hmrate_exo : public mole_reaction
00215 {
00216 typedef mole_reaction_hmrate_exo T;
00217 public:
00218 virtual T* Create() const {return new T;}
00219
00220 virtual const char* name() {return "hmrate_exo";}
00221
00222 double rk() const
00223 {
00224 double te;
00225
00226 DEBUG_ENTRY("hmrate_exo()");
00227
00228 te = phycon.te+noneq_offset(this);
00229
00230 if( this->c < 0. )
00231 te = MIN2( te, -10. * this->c );
00232
00233 return pow(te/300.,this->b)*exp(-te/this->c);
00234 }
00235
00236 };
00237
00238 class mole_reaction_hmrate : public mole_reaction
00239 {
00240 typedef mole_reaction_hmrate T;
00241 public:
00242 virtual T* Create() const {return new T;}
00243 virtual const char* name() {return "hmrate";}
00244 double rk() const
00245 {
00246 return hmrate(this);
00247 }
00248 };
00249
00250 class mole_reaction_constrate : public mole_reaction
00251 {
00252 typedef mole_reaction_constrate T;
00253 public:
00254 virtual T* Create() const {return new T;}
00255 virtual const char* name() {return "constrate";}
00256 double rk() const
00257 {
00258 return 1.;
00259 }
00260 };
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274 }
00275 #include "rfield.h"
00276 namespace {
00277 double th85rate(const mole_reaction *rate)
00278 {
00279 double rk;
00280
00281 DEBUG_ENTRY("th85rate()");
00282
00283 if (!mole_global.lgLeidenHack || rate->c == 0.0)
00284 {
00285 rk = hmi.UV_Cont_rel2_Habing_TH85_depth/1.66;
00286 }
00287 else
00288 {
00289 rk = hmi.UV_Cont_rel2_Habing_TH85_face/1.66*exp(-(rate->c*rfield.extin_mag_V_point));
00290 }
00291
00292 return rk;
00293 }
00294 class mole_reaction_th85rate : public mole_reaction
00295 {
00296 typedef mole_reaction_th85rate T;
00297 public:
00298 virtual T* Create() const {return new T;}
00299 virtual const char* name() {return "th85rate";}
00300 double rk() const
00301 {
00302 return th85rate(this);
00303 }
00304 };
00305 }
00306
00307 #include "secondaries.h"
00308 namespace {
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322 double crnurate_noalbedo(const mole_reaction *)
00323 {
00324 return 2.17*secondaries.csupra[ipHYDROGEN][0];
00325 }
00326 class mole_reaction_crnurate_noalbedo : public mole_reaction
00327 {
00328 typedef mole_reaction_crnurate_noalbedo T;
00329 public:
00330 virtual T* Create() const {return new T;}
00331
00332 virtual const char* name() {return "crnurate_noalbedo";}
00333 double rk() const
00334 {
00335 return crnurate_noalbedo(this);
00336 }
00337 };
00338
00339 class mole_reaction_crnurate : public mole_reaction
00340 {
00341 typedef mole_reaction_crnurate T;
00342 public:
00343 virtual T* Create() const {return new T;}
00344 virtual const char* name() {return "crnurate";}
00345 double rk() const
00346 {
00347 return crnurate_noalbedo(this)/(1.0-albedo);
00348 }
00349 };
00350 }
00351 #include "hextra.h"
00352 namespace {
00353
00354 class mole_reaction_cryden_ov_bg : public mole_reaction
00355 {
00356 typedef mole_reaction_cryden_ov_bg T;
00357 public:
00358 virtual T* Create() const {return new T;}
00359 virtual const char* name() {return "cryden_ov_bg";}
00360 double rk() const
00361 {
00362 return hextra.cryden_ov_background;
00363 }
00364 };
00365
00366 class mole_reaction_co_lnu_c_o_lnu : public mole_reaction
00367 {
00368 typedef mole_reaction_co_lnu_c_o_lnu T;
00369 public:
00370 virtual T* Create() const {return new T;}
00371 virtual const char* name() {return "co_lnu_c_o_lnu";}
00372 double rk() const
00373 {
00374 double val = 0;
00375 int ns, ion;
00376
00377
00378
00379
00380
00381 DEBUG_ENTRY("co_lnu_c_o_lnu()");
00382
00383 for( ns=0; ns<2; ++ns )
00384 {
00385 ion = 0;
00386 val += ionbal.PhotoRate_Shell[ipCARBON][ion][ns][0];
00387 val += ionbal.PhotoRate_Shell[ipOXYGEN][ion][ns][0];
00388 }
00389
00390 return val;
00391 }
00392 };
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427 class mole_reaction_vib_evap : public mole_reaction
00428 {
00429 typedef mole_reaction_vib_evap T;
00430 public:
00431 virtual T* Create() const {return new T;}
00432 virtual const char* name() {return "vib_evap";}
00433 double rk() const
00434 {
00435 double binding_energy, exponent, vib_freq;
00436
00437 DEBUG_ENTRY("vib_evap()");
00438
00439 exponent = 0.0;
00440
00441 binding_energy = this->b;
00442 double bin_total=0.0;
00443 for( size_t nd=0; nd < gv.bin.size() ; nd++ )
00444 {
00445 double bin_area = gv.bin[nd]->IntArea*gv.bin[nd]->cnv_H_pCM3;
00446 exponent += exp(-binding_energy/gv.bin[nd]->tedust)*bin_area;
00447 bin_total += bin_area;
00448 }
00449 exponent /= bin_total;
00450 const double surface_density_of_sites = 1.5e15;
00451
00452
00453 vib_freq = sqrt(2*surface_density_of_sites*(0.3*BOLTZMANN)*binding_energy/(PI*PI*this->reactants[0]->mole_mass));
00454
00455
00456
00457
00458
00459
00460
00461
00463
00464
00465
00466
00467 return vib_freq*exponent +sexp( 555.89/phycon.sqrte - 5.55 );
00468 }
00469 };
00470
00471 class mole_reaction_grn_abs : public mole_reaction
00472 {
00473 typedef mole_reaction_grn_abs T;
00474 public:
00475 virtual T* Create() const {return new T;}
00476 virtual const char* name() {return "grn_abs";}
00477 double rk() const
00478 {
00479 double mass;
00480
00481 DEBUG_ENTRY("grn_abs()");
00482
00483
00484 if (this->reactants[0]->n_nuclei() != 0)
00485 mass = this->reactants[0]->mole_mass;
00486 else
00487 mass = this->reactants[1]->mole_mass;
00488
00489 return sqrt(8.*BOLTZMANN*phycon.te/(PI*mass));
00490 }
00491 };
00492
00493 class mole_reaction_grn_react : public mole_reaction
00494 {
00495 typedef mole_reaction_grn_react T;
00496 public:
00497 virtual T* Create() const {return new T;}
00498 virtual const char* name() {return "grn_react";}
00499 double rk() const
00500 {
00501 return grn_react(this);
00502 }
00503 };
00504 double grn_react(const mole_reaction *rate)
00505 {
00506
00507
00508
00509
00510
00511
00512
00513
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523 DEBUG_ENTRY("grn_react()");
00524
00525 double quant_barrier = 1e-8;
00526 double surface_density_of_sites = 1.5e15;
00527 fixit();
00528 ASSERT( rate->nreactants == 2 );
00529 double E_i = rate->reactants[0]->form_enthalpy;
00530 double E_j = rate->reactants[1]->form_enthalpy;
00531 double activ_barrier = rate->c;
00532
00533
00534 fixit();
00535 double vib_freq_i = sqrt(2*surface_density_of_sites*(0.3*BOLTZMANN)*E_i/(PI*PI*rate->reactants[0]->mole_mass));
00536 double vib_freq_j = sqrt(2*surface_density_of_sites*(0.3*BOLTZMANN)*E_j/(PI*PI*rate->reactants[1]->mole_mass));
00537
00538 double Exp_i = 0.0;
00539 double Exp_j = 0.0;
00540 double dust_density = 0.0;
00541
00542
00543 fixit();
00544
00545 for( size_t nd=0; nd < gv.bin.size(); nd++ )
00546 {
00547 double bin_density = gv.bin[nd]->IntArea*gv.bin[nd]->cnv_H_pCM3;
00548 Exp_i += exp(-E_i/gv.bin[nd]->tedust)*bin_density;
00549 Exp_j += exp(-E_j/gv.bin[nd]->tedust)*bin_density;
00550 dust_density += bin_density/(4*1e-10);
00551 }
00552
00553 ASSERT(fp_equal((realnum)dust_density,(realnum)(mole.grain_area/1e-10)));
00554
00555 double total_sites = 4.*mole.grain_area*surface_density_of_sites;
00556
00557
00558 double diff_i = vib_freq_i*Exp_i/total_sites;
00559 double diff_j = vib_freq_j*Exp_j/total_sites;
00560
00561
00562
00563 double scale = exp(-2*(quant_barrier/(HBAReV*EN1EV))*sqrt(2.0*rate->reduced_mass*0.3*BOLTZMANN*activ_barrier));
00564
00565
00566 return scale*(diff_i + diff_j)/SDIV(dust_density);
00567 }
00568
00569 class mole_reaction_grn_photo : public mole_reaction
00570 {
00571 typedef mole_reaction_grn_photo T;
00572 public:
00573 virtual T* Create() const {return new T;}
00574 virtual const char* name() {return "grn_photo";}
00575 double rk() const
00576 {
00577
00578
00579
00580
00581
00582 DEBUG_ENTRY("grn_photo()");
00583
00584
00585
00586
00587
00588
00589 fixit();
00590
00591
00592
00593
00594 return 2e-15 * hmi.UV_Cont_rel2_Draine_DB96_depth *(1.232e7f * 1.71f);
00595 }
00596 };
00597 }
00598 #include "rt.h"
00599 namespace {
00600 class mole_reaction_th85rate_co : public mole_reaction
00601 {
00602 typedef mole_reaction_th85rate_co T;
00603 public:
00604 virtual T* Create() const {return new T;}
00605
00606 virtual const char* name() {return "th85rate_co";}
00607
00608 double rk() const
00609 {
00610 double esc_co, column;
00611
00612
00613
00614
00615
00616
00617
00618
00619
00620
00621
00622
00623
00624
00625
00626
00627
00628
00629
00630 if (this->reactants[0]->n_nuclei() != 0)
00631 column = mole.species[ this->reactants[0]->index ].column;
00632 else
00633 column = mole.species[ this->reactants[1]->index ].column;
00634
00635 esc_co = 4.4e-15 * column /
00636
00637 (GetDopplerWidth(dense.AtomicWeight[5]+dense.AtomicWeight[7])/1e5) /
00638
00639 (1. + phycon.sqrte*0.6019);
00640 return esca0k2(esc_co)*th85rate(this);
00641 }
00642 };
00643
00644 class mole_reaction_oh_c2h2_co_ch3 : public mole_reaction
00645 {
00646 typedef mole_reaction_oh_c2h2_co_ch3 T;
00647 public:
00648 virtual T* Create() const {return new T;}
00649 virtual const char* name() {return "oh_c2h2_co_ch3";}
00650 double rk() const
00651 {
00652
00653
00654 if(phycon.te > 500)
00655 {
00656 return hmrate(this);
00657 }
00658 else
00659 {
00660 return 6.3E-18;
00661 }
00662 }
00663 };
00664
00665
00666 class mole_reaction_h_hnc_hcn_h : public mole_reaction
00667 {
00668 typedef mole_reaction_h_hnc_hcn_h T;
00669 public:
00670 virtual T* Create() const {return new T;}
00671
00672 virtual const char* name() {return "h_hnc_hcn_h";}
00673 double rk() const
00674 {
00675 if(phycon.te > 100)
00676 {
00677 return hmrate(this);
00678 }
00679 else
00680 {
00681 return 1e-15;
00682 }
00683 }
00684 };
00685 }
00686 #include "iso.h"
00687 #include "h2.h"
00688 double frac_H2star_hminus(void)
00689 {
00690
00691 if( h2.lgEnabled && h2.lgEvaluated && hmi.lgH2_Chemistry_BigH2 )
00692 {
00693 return hmi.H2star_forms_hminus /
00694 SDIV(hmi.H2star_forms_hminus+hmi.H2_forms_hminus);
00695
00696
00697
00698
00699 }
00700 else
00701 {
00702
00703
00704
00705
00706
00707
00708 return 1. - 4.938e-6;
00709 }
00710 }
00711
00712 namespace {
00713 double frac_H2star_grains(void)
00714 {
00715
00716 if( h2.lgEnabled && h2.lgEvaluated && hmi.lgH2_Chemistry_BigH2 )
00717 {
00718 return hmi.H2star_forms_grains /
00719 SDIV(hmi.H2star_forms_grains+hmi.H2_forms_grains);
00720
00721
00722
00723
00724 }
00725 else
00726 {
00727
00728
00729
00730
00731
00732
00733 return 0.9416;
00734 }
00735 }
00736
00737 class mole_reaction_h2ph3p : public mole_reaction
00738 {
00739 typedef mole_reaction_h2ph3p T;
00740 public:
00741 virtual T* Create() const {return new T;}
00742 virtual const char* name() {return "h2ph3p";}
00743 double rk() const
00744 {
00745
00749 if (!mole_global.lgLeidenHack)
00750 return 1.40e-9*(1. - sexp(9940./phycon.te));
00751 else
00752 return 2.08e-9;
00753 }
00754 };
00755 }
00756
00757 #include "opacity.h"
00758 #include "gammas.h"
00759 #include "atmdat.h"
00760
00761 namespace {
00762 class mole_reaction_hopexch : public mole_reaction
00763 {
00764 typedef mole_reaction_hopexch T;
00765 public:
00766 virtual T* Create() const {return new T;}
00767 virtual const char* name() {return "hopexch";}
00768 double rk() const
00769 {
00770 return atmdat.HCharExcRecTo[ipOXYGEN][0];
00771 }
00772 };
00773
00774 class mole_reaction_hpoexch : public mole_reaction
00775 {
00776 typedef mole_reaction_hpoexch T;
00777 public:
00778 virtual T* Create() const {return new T;}
00779 virtual const char* name() {return "hpoexch";}
00780 double rk() const
00781 {
00782 return atmdat.HCharExcIonOf[ipOXYGEN][0];
00783 }
00784 };
00785
00786 class mole_reaction_hmattach : public mole_reaction
00787 {
00788 typedef mole_reaction_hmattach T;
00789 public:
00790 virtual T* Create() const {return new T;}
00791 virtual const char* name() {return "hmattach";}
00792 double rk() const
00793 {
00794 return hmirat(phycon.te) + (hmi.HMinus_induc_rec_rate)/SDIV(dense.eden);
00795 }
00796 };
00797
00798 class mole_reaction_hmihph2p : public mole_reaction
00799 {
00800 typedef mole_reaction_hmihph2p T;
00801 public:
00802 virtual T* Create() const {return new T;}
00803 virtual const char* name() {return "hmihph2p";}
00804 double rk() const
00805 {
00806
00807
00808
00809
00810
00811 if(phycon.te <= 7891.)
00812 {
00813
00814 return 6.9e-9 / (phycon.te30 * phycon.te05);
00815 }
00816 else
00817 {
00818
00819
00820 return 9.6e-7 / phycon.te90;
00821 }
00822 }
00823 };
00824
00825 class mole_reaction_hmphoto : public mole_reaction
00826 {
00827 typedef mole_reaction_hmphoto T;
00828 public:
00829 virtual T* Create() const {return new T;}
00830 virtual const char* name() {return "hmphoto";}
00831 double rk() const
00832 {
00833 return hmi.HMinus_photo_rate;
00834 }
00835 };
00836
00837 double cionhm(const mole_reaction *)
00838 {
00839
00840 if( phycon.te < 3074. )
00841 {
00842 return 1.46e-32*(powi(phycon.te,6))*phycon.sqrte*hmi.exphmi;
00843 }
00844 else if( phycon.te >= 3074. && phycon.te < 30000. )
00845 {
00846
00847
00848
00849
00850 return 5.9e-19*phycon.tesqrd*phycon.sqrte*phycon.te05;
00851 }
00852 else
00853 {
00854 return 1.54e-7;
00855 }
00856
00857 }
00858 class mole_reaction_cionhm : public mole_reaction
00859 {
00860 typedef mole_reaction_cionhm T;
00861 public:
00862 virtual T* Create() const {return new T;}
00863 virtual const char* name() {return "cionhm";}
00864 double rk() const
00865 {
00866 return cionhm(this);
00867 }
00868 };
00869
00870 double assoc_detach(void)
00871 {
00872
00873
00874
00875
00876
00877
00878
00879
00880
00881
00882 double y , x;
00883 x = MAX2(10., phycon.te );
00884 x = MIN2(1e4, x );
00885 y=545969508.1323510+x*71239.23653059864;
00886 return 1./y;
00887 }
00888
00889 class mole_reaction_c3bod : public mole_reaction
00890 {
00891 typedef mole_reaction_c3bod T;
00892 public:
00893 virtual T* Create() const {return new T;}
00894 virtual const char* name() {return "c3bod";}
00895 double rk() const
00896 {
00897 return cionhm(this)*hmi.rel_pop_LTE_Hmin;
00898 }
00899 };
00900
00901 class mole_reaction_asdfg : public mole_reaction
00902 {
00903 typedef mole_reaction_asdfg T;
00904 public:
00905 virtual T* Create() const {return new T;}
00906 virtual const char* name() {return "asdfg";}
00907 double rk() const
00908 {
00909 return assoc_detach()*(1-frac_H2star_hminus());
00910 }
00911 };
00912
00913 class mole_reaction_asdfs : public mole_reaction
00914 {
00915 typedef mole_reaction_asdfs T;
00916 public:
00917 virtual T* Create() const {return new T;}
00918 virtual const char* name() {return "asdfs";}
00919 double rk() const
00920 {
00921 return assoc_detach()*frac_H2star_hminus();
00922 }
00923 };
00924
00925 class mole_reaction_asdbg : public mole_reaction
00926 {
00927 typedef mole_reaction_asdbg T;
00928 public:
00929 virtual T* Create() const {return new T;}
00930 virtual const char* name() {return "asdbg";}
00931 double rk() const
00932 {
00933 if( hmi.rel_pop_LTE_H2g > 0. )
00934 {
00935 return assoc_detach()*hmi.rel_pop_LTE_Hmin/hmi.rel_pop_LTE_H2g*
00936 (1.-frac_H2star_hminus());
00937 }
00938 else
00939 {
00940 return 0.;
00941 }
00942 }
00943 };
00944
00945 class mole_reaction_asdbs : public mole_reaction
00946 {
00947 typedef mole_reaction_asdbs T;
00948 public:
00949 virtual T* Create() const {return new T;}
00950 virtual const char* name() {return "asdbs";}
00951 double rk() const
00952 {
00953
00954
00955
00956
00957
00958
00959
00960 if( hmi.rel_pop_LTE_H2s > 0. )
00961 {
00962 double ratio = mole_get_equilibrium_constant("H-,H=>H2,e-");
00963 return assoc_detach() * ratio *
00964 frac_H2star_hminus();
00965 }
00966 else
00967 {
00968 return 0.;
00969 }
00970 }
00971 };
00972
00973 class mole_reaction_bhneut : public mole_reaction
00974 {
00975 typedef mole_reaction_bhneut T;
00976 public:
00977 virtual T* Create() const {return new T;}
00978 virtual const char* name() {return "bhneut";}
00979 double rk() const
00980 {
00981 if( phycon.te > 1000. && dense.xIonDense[ipHYDROGEN][0] > 0.0 )
00982 {
00983 double ratio = mole_get_equilibrium_constant("H-,H+=>H,H");
00984
00985
00986
00988
00989 return (hneut(this)*ratio*
00990 (iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH3s].Pop()+
00991 iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH3p].Pop()+
00992 iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH3d].Pop()) /
00993 SDIV(dense.xIonDense[ipHYDROGEN][0]));
00994 }
00995 else
00996 {
00997 return 0.;
00998 }
00999 }
01000 };
01001
01002 double hneut(const mole_reaction *)
01003 {
01004 if( phycon.te < 14125. )
01005 {
01006
01007
01008 return 1.4e-7*pow(phycon.te/300,-0.487)*exp(phycon.te/29300);
01009 }
01010 else
01011 {
01012 return 3.4738192887404660e-008;
01013 }
01014 }
01015 class mole_reaction_hneut : public mole_reaction
01016 {
01017 typedef mole_reaction_hneut T;
01018 public:
01019 virtual T* Create() const {return new T;}
01020 virtual const char* name() {return "hneut";}
01021
01022 double rk() const
01023 {
01024 return hneut(this);
01025 }
01026 };
01027
01028 class mole_reaction_h2_spon_diss : public mole_reaction
01029 {
01030 typedef mole_reaction_h2_spon_diss T;
01031 public:
01032 virtual T* Create() const {return new T;}
01033 virtual const char* name() {return "h2_spon_diss";}
01034 double rk() const
01035 {
01036 return h2.spon_diss_tot;
01037 }
01038 };
01039
01040 double grnh2tot(const mole_reaction *)
01041 {
01042
01043
01044
01045 fixit();
01046 if( mole.grain_area*dense.xIonDense[ipHYDROGEN][0]>0 )
01047 return gv.rate_h2_form_grains_used_total/(mole.grain_area*dense.xIonDense[ipHYDROGEN][0]);
01048 else
01049 return 0.;
01050 }
01051 class mole_reaction_grnh2tot : public mole_reaction
01052 {
01053 public:
01054 virtual const char* name() {return "grnh2tot";}
01055 double rk() const
01056 {
01057 return grnh2tot(this);
01058 }
01059 };
01060
01061 class mole_reaction_grnh2 : public mole_reaction
01062 {
01063 typedef mole_reaction_grnh2 T;
01064 public:
01065 virtual T* Create() const {return new T;}
01066 virtual const char* name() {return "grnh2";}
01067 double rk() const
01068 {
01069 return grnh2tot(this)*(1.-frac_H2star_grains());
01070 }
01071 };
01072
01073 class mole_reaction_grnh2s : public mole_reaction
01074 {
01075 typedef mole_reaction_grnh2s T;
01076 public:
01077 virtual T* Create() const {return new T;}
01078 virtual const char* name() {return "grnh2s";}
01079 double rk() const
01080 {
01081 return grnh2tot(this)*frac_H2star_grains();
01082 }
01083 };
01084
01085 class mole_reaction_radasc : public mole_reaction
01086 {
01087 typedef mole_reaction_radasc T;
01088 public:
01089 virtual T* Create() const {return new T;}
01090 virtual const char* name() {return "radasc";}
01091 double rk() const
01092 {
01093
01094
01095
01096 if( dense.xIonDense[ipHYDROGEN][0] > 0. )
01097 {
01098 return hmrate(this) *
01099 (iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop() / dense.xIonDense[ipHYDROGEN][0] ) *
01100 (iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH2p].Pop() + iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH2s].Pop())/
01101 dense.xIonDense[ipHYDROGEN][0];
01102 }
01103 else
01104 return 0.;
01105 }
01106 };
01107
01108 class mole_reaction_assoc_ion : public mole_reaction
01109 {
01110 typedef mole_reaction_assoc_ion T;
01111 public:
01112 virtual T* Create() const {return new T;}
01113 virtual const char* name() {return "assoc_ion";}
01114 double rk() const
01115 {
01116
01117
01118 if( dense.xIonDense[ipHYDROGEN][0] > 0. )
01119 {
01120 return hmrate(this) *
01121 (iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop() / dense.xIonDense[ipHYDROGEN][0] ) *
01122 (iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH2p].Pop() + iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH2s].Pop())/
01123 dense.xIonDense[ipHYDROGEN][0];
01124 }
01125 else
01126 return 0.;
01127 }
01128 };
01129
01130
01131 double rh2g_dis_h(const mole_reaction *)
01132 {
01133 if( h2.lgEnabled && h2.lgEvaluated && hmi.lgH2_Chemistry_BigH2 )
01134 {
01135 return h2.Average_collH_dissoc_g;
01136 }
01137 else
01138 {
01139 double corr = MIN2(6.,14.44-phycon.alogte*3.08);
01140
01141 if(corr > 0.)
01142 corr = pow(10.,corr*findspecieslocal("H")->den/(findspecieslocal("H")->den+1.6e4));
01143 else
01144 corr = 1.;
01145
01146 return 1.55e-8/phycon.sqrte*sexp(65107./phycon.te)* corr;
01147 }
01148 }
01149 class mole_reaction_rh2g_dis_h : public mole_reaction
01150 {
01151 typedef mole_reaction_rh2g_dis_h T;
01152 public:
01153 virtual T* Create() const {return new T;}
01154 virtual const char* name() {return "rh2g_dis_h";}
01155 double rk() const
01156 {
01157 return rh2g_dis_h(this);
01158 }
01159 };
01160
01161 double rh2s_dis_h(const mole_reaction *rate)
01162 {
01163 if( h2.lgEnabled && h2.lgEvaluated && hmi.lgH2_Chemistry_BigH2 )
01164 {
01165 return h2.Average_collH_dissoc_s;
01166 }
01167 else
01168 {
01169 ASSERT( fp_equal( rate->a, 1. ) );
01170 return HMRATE(4.67e-7,-1.,5.5e4);
01171 }
01172 }
01173 class mole_reaction_rh2s_dis_h : public mole_reaction
01174 {
01175 typedef mole_reaction_rh2s_dis_h T;
01176 public:
01177 virtual T* Create() const {return new T;}
01178 virtual const char* name() {return "rh2s_dis_h";}
01179 double rk() const
01180 {
01181 return rh2s_dis_h(this);
01182 }
01183 };
01184
01185 double rh2g_dis_h2(const mole_reaction *rate)
01186 {
01187 if( h2.lgEnabled && h2.lgEvaluated && hmi.lgH2_Chemistry_BigH2 )
01188 {
01189 return h2.Average_collH2_dissoc_g;
01190 }
01191 else
01192 {
01193
01194 ASSERT( fp_equal( rate->a, 1. ) );
01195 return HMRATE(5.5e-29*0.5/(SAHA*3.634e-5)*sqrt(300.),0.5,5.195e4);
01196 }
01197 }
01198 class mole_reaction_rh2g_dis_h2 : public mole_reaction
01199 {
01200 typedef mole_reaction_rh2g_dis_h2 T;
01201 public:
01202 virtual T* Create() const {return new T;}
01203 virtual const char* name() {return "rh2g_dis_h2";}
01204 double rk() const
01205 {
01206 return rh2g_dis_h2(this);
01207 }
01208 };
01209
01210 double rh2s_dis_h2(const mole_reaction *rate)
01211 {
01212 if( h2.lgEnabled && h2.lgEvaluated && hmi.lgH2_Chemistry_BigH2 )
01213 {
01214 return h2.Average_collH2_dissoc_s;
01215 }
01216 else
01217 {
01218 ASSERT( fp_equal( rate->a, 1. ) );
01219 return HMRATE(1e-11,0.,0.);
01220 }
01221 }
01222 class mole_reaction_rh2s_dis_h2 : public mole_reaction
01223 {
01224 typedef mole_reaction_rh2s_dis_h2 T;
01225 public:
01226 virtual T* Create() const {return new T;}
01227 virtual const char* name() {return "rh2s_dis_h2";}
01228 double rk() const
01229 {
01230 return rh2s_dis_h2(this);
01231 }
01232 };
01233
01234 double rh2s_dis_h2_nodeexcit(const mole_reaction *rate)
01235 {
01236 if( h2.lgEnabled && h2.lgEvaluated && hmi.lgH2_Chemistry_BigH2 )
01237 {
01238 return h2.Average_collH2_dissoc_s;
01239 }
01240 else
01241 {
01242 ASSERT( fp_equal( rate->a, 1. ) );
01243 return HMRATE(1e-11,0.,2.18e4);
01244 }
01245 }
01246 class mole_reaction_rh2s_dis_h2_nodeexcit : public mole_reaction
01247 {
01248 typedef mole_reaction_rh2s_dis_h2_nodeexcit T;
01249 public:
01250 virtual T* Create() const {return new T;}
01251 virtual const char* name() {return "rh2s_dis_h2_nodeexcit";}
01252 double rk() const
01253 {
01254 return rh2s_dis_h2_nodeexcit(this);
01255 }
01256 };
01257
01258 double bh2g_dis_h(const mole_reaction *rate)
01259 {
01260 return rh2g_dis_h(rate)*hmi.rel_pop_LTE_H2g;
01261 }
01262 class mole_reaction_bh2g_dis_h : public mole_reaction
01263 {
01264 typedef mole_reaction_bh2g_dis_h T;
01265 public:
01266 virtual T* Create() const {return new T;}
01267 virtual const char* name() {return "bh2g_dis_h";}
01268 double rk() const
01269 {
01270 return bh2g_dis_h(this);
01271 }
01272 };
01273
01274 double bh2s_dis_h(const mole_reaction *rate)
01275 {
01276 return rh2s_dis_h(rate)*hmi.rel_pop_LTE_H2s;
01277 }
01278 class mole_reaction_bh2s_dis_h : public mole_reaction
01279 {
01280 typedef mole_reaction_bh2s_dis_h T;
01281 public:
01282 virtual T* Create() const {return new T;}
01283 virtual const char* name() {return "bh2s_dis_h";}
01284 double rk() const
01285 {
01286 return bh2s_dis_h(this);
01287 }
01288 };
01289
01290 double bh2g_dis_h2(const mole_reaction *rate)
01291 {
01292 return rh2g_dis_h2(rate)*hmi.rel_pop_LTE_H2g;
01293 }
01294 class mole_reaction_bh2g_dis_h2 : public mole_reaction
01295 {
01296 typedef mole_reaction_bh2g_dis_h2 T;
01297 public:
01298 virtual T* Create() const {return new T;}
01299 virtual const char* name() {return "bh2g_dis_h2";}
01300 double rk() const
01301 {
01302 return bh2g_dis_h2(this);
01303 }
01304 };
01305
01306 class mole_reaction_bh2s_dis_h2 : public mole_reaction
01307 {
01308 typedef mole_reaction_bh2s_dis_h2 T;
01309 public:
01310 virtual T* Create() const {return new T;}
01311 virtual const char* name() {return "bh2s_dis_h2";}
01312 double rk() const
01313 {
01314 return rh2s_dis_h2(this)*hmi.rel_pop_LTE_H2s;
01315 }
01316 };
01317
01318 class mole_reaction_h2photon : public mole_reaction
01319 {
01320 typedef mole_reaction_h2photon T;
01321 public:
01322 virtual T* Create() const {return new T;}
01323 virtual const char* name() {return "h2photon";}
01324 double rk() const
01325 {
01326 return h2.photoionize_rate;
01327 }
01328 };
01329
01330 class mole_reaction_h2crphot : public mole_reaction
01331 {
01332 typedef mole_reaction_h2crphot T;
01333 public:
01334 virtual T* Create() const {return new T;}
01335 virtual const char* name() {return "h2crphot";}
01336 double rk() const
01337 {
01338 return secondaries.csupra[ipHYDROGEN][0]*2.02;
01339 }
01340 };
01341
01342 class mole_reaction_h2crphh : public mole_reaction
01343 {
01344 typedef mole_reaction_h2crphh T;
01345 public:
01346 virtual T* Create() const {return new T;}
01347 virtual const char* name() {return "h2crphh";}
01348 double rk() const
01349 {
01350 if( h2.lgEnabled && h2.lgEvaluated && hmi.lgH2_Chemistry_BigH2 )
01351 {
01352 return secondaries.x12tot*10.;
01353 }
01354 else
01355 {
01356 return secondaries.x12tot*3.;
01357 }
01358 }
01359 };
01360
01361 class mole_reaction_h2scrphh : public mole_reaction
01362 {
01363 typedef mole_reaction_h2scrphh T;
01364 public:
01365 virtual T* Create() const {return new T;}
01366 virtual const char* name() {return "h2scrphh";}
01367 double rk() const
01368 {
01369 if( h2.lgEnabled && h2.lgEvaluated && hmi.lgH2_Chemistry_BigH2 )
01370 {
01371 return secondaries.x12tot*10.;
01372 }
01373 else
01374 {
01375 return secondaries.x12tot*3.;
01376 }
01377 }
01378 };
01379
01380 class mole_reaction_radath : public mole_reaction
01381 {
01382 typedef mole_reaction_radath T;
01383 public:
01384 virtual T* Create() const {return new T;}
01385 virtual const char* name() {return "radath";}
01386 double rk() const
01387 {
01388 return MAX2(0.,2.325*MIN2(5000.,phycon.te)-1375.)*1e-20;
01389 }
01390 };
01391
01392 class mole_reaction_gamtwo : public mole_reaction
01393 {
01394 typedef mole_reaction_gamtwo T;
01395 public:
01396 virtual T* Create() const {return new T;}
01397 virtual const char* name() {return "gamtwo";}
01398 double rk() const
01399 {
01400 t_phoHeat dummy;
01401 return GammaK(opac.ih2pnt[0],opac.ih2pnt[1],opac.ih2pof,1.,&dummy);
01402 }
01403 };
01404
01405 class mole_reaction_hlike_phot : public mole_reaction
01406 {
01407 typedef mole_reaction_hlike_phot T;
01408 public:
01409 virtual T* Create() const {return new T;}
01410 virtual const char* name() {return "hlike_phot";}
01411 double rk() const
01412 {
01413
01414
01415
01416
01417
01418 if( !conv.nTotalIoniz )
01419 iso_photo( ipH_LIKE, ipHYDROGEN );
01420 return iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].gamnc;
01421 }
01422 };
01423
01424 class mole_reaction_h2s_sp_decay : public mole_reaction
01425 {
01426 typedef mole_reaction_h2s_sp_decay T;
01427 public:
01428 virtual T* Create() const {return new T;}
01429 virtual const char* name() {return "h2s_sp_decay";}
01430 double rk() const
01431 {
01432
01433
01434 if( h2.lgEnabled && h2.lgEvaluated && hmi.lgH2_Chemistry_BigH2 )
01435 {
01436 return h2.Average_A;
01437 }
01438 else
01439 {
01440 return 2e-7;
01441 }
01442 }
01443 };
01444
01445 double h2_collh2_deexc(const mole_reaction *)
01446 {
01447 if( h2.lgEnabled && h2.lgEvaluated && hmi.lgH2_Chemistry_BigH2 )
01448 {
01449 return h2.Average_collH2_deexcit;
01450 }
01451 else
01452 {
01453 return ((1.4e-12*phycon.sqrte * sexp( 18100./(phycon.te + 1200.) ))/6.);
01454 }
01455 }
01456 class mole_reaction_h2_collh2_deexc : public mole_reaction
01457 {
01458 typedef mole_reaction_h2_collh2_deexc T;
01459 public:
01460 virtual T* Create() const {return new T;}
01461 virtual const char* name() {return "h2_collh2_deexc";}
01462 double rk() const
01463 {
01464 return h2_collh2_deexc(this);
01465 }
01466 };
01467
01468 double h2_collh_deexc(const mole_reaction *)
01469 {
01470 if ( h2.lgEnabled && h2.lgEvaluated && hmi.lgH2_Chemistry_BigH2 )
01471 {
01472 return h2.Average_collH_deexcit;
01473 }
01474 else
01475 {
01476 return ((1e-12*phycon.sqrte * sexp(1000./phycon.te ))/6.);
01477 }
01478 }
01479 class mole_reaction_h2_collh_deexc : public mole_reaction
01480 {
01481 typedef mole_reaction_h2_collh_deexc T;
01482 public:
01483 virtual T* Create() const {return new T;}
01484 virtual const char* name() {return "h2_collh_deexc";}
01485 double rk() const
01486 {
01487 return h2_collh_deexc(this);
01488 }
01489 };
01490
01491 double h2_collh2_excit(const mole_reaction *rate)
01492 {
01493
01494 if ( h2.lgEnabled && h2.lgEvaluated && hmi.lgH2_Chemistry_BigH2 )
01495 {
01496 return h2.Average_collH2_excit;
01497 }
01498 else
01499 {
01500 return h2_collh2_deexc(rate)*sexp( 30172./phycon.te);
01501 }
01502 }
01503 class mole_reaction_h2_collh2_excit : public mole_reaction
01504 {
01505 typedef mole_reaction_h2_collh2_excit T;
01506 public:
01507 virtual T* Create() const {return new T;}
01508 virtual const char* name() {return "h2_collh2_excit";}
01509 double rk() const
01510 {
01511 return h2_collh2_excit(this);
01512 }
01513 };
01514
01515 double h2_collh_excit(const mole_reaction *rate)
01516 {
01517
01518 if ( h2.lgEnabled && h2.lgEvaluated && hmi.lgH2_Chemistry_BigH2 )
01519 {
01520 return h2.Average_collH_excit;
01521 }
01522 else
01523 {
01524 return h2_collh_deexc(rate)*sexp( 30172./phycon.te);
01525 }
01526 }
01527 class mole_reaction_h2_collh_excit : public mole_reaction
01528 {
01529 typedef mole_reaction_h2_collh_excit T;
01530 public:
01531 virtual T* Create() const {return new T;}
01532 virtual const char* name() {return "h2_collh_excit";}
01533 double rk() const
01534 {
01535 return h2_collh_excit(this);
01536 }
01537 };
01538
01539 class mole_reaction_h2gexcit : public mole_reaction
01540 {
01541 typedef mole_reaction_h2gexcit T;
01542 public:
01543 virtual T* Create() const {return new T;}
01544 virtual const char* name() {return "h2gexcit";}
01545 double rk() const
01546 {
01547 return hmi.H2_H2g_to_H2s_rate_used;
01548 }
01549 };
01550
01551 class mole_reaction_h2sdissoc : public mole_reaction
01552 {
01553 typedef mole_reaction_h2sdissoc T;
01554 public:
01555 virtual T* Create() const {return new T;}
01556 virtual const char* name() {return "h2sdissoc";};
01557 double rk() const
01558 {
01559
01560
01561
01562
01563
01564
01565 return hmi.H2_photodissoc_used_H2s +hmi.H2_Solomon_dissoc_rate_used_H2s;
01566 }
01567 };
01568
01569 class mole_reaction_h2gdissoc : public mole_reaction
01570 {
01571 typedef mole_reaction_h2gdissoc T;
01572 public:
01573 virtual T* Create() const {return new T;}
01574 virtual const char* name() {return "h2gdissoc";}
01575 double rk() const
01576 {
01577
01578
01579
01580
01581
01582
01583 return hmi.H2_photodissoc_used_H2g+hmi.H2_Solomon_dissoc_rate_used_H2g;
01584 }
01585 };
01586
01587 class mole_reaction_hd_photodissoc : public mole_reaction
01588 {
01589 typedef mole_reaction_hd_photodissoc T;
01590 public:
01591 virtual T* Create() const {return new T;}
01592 virtual const char* name() {return "hd_photodissoc";}
01593 double rk() const
01594 {
01595 return hd.photodissoc_BigH2_H2g + hd.Solomon_dissoc_rate_g;
01596 }
01597 };
01598 }
01599
01600
01601 double hmirat(double te)
01602 {
01603 double hmirat_v;
01604
01605 DEBUG_ENTRY( "hmirat()" );
01606
01607
01608 if( te < 31.62 )
01609 {
01610 hmirat_v = 8.934e-18*phycon.sqrte*phycon.te003*phycon.te001*
01611 phycon.te001;
01612 }
01613 else if( te < 90. )
01614 {
01615 hmirat_v = 5.159e-18*phycon.sqrte*phycon.te10*phycon.te03*
01616 phycon.te03*phycon.te003*phycon.te001;
01617 }
01618 else if( te < 1200. )
01619 {
01620 hmirat_v = 2.042e-18*te/phycon.te10/phycon.te03;
01621 }
01622 else if( te < 3800. )
01623 {
01624 hmirat_v = 8.861e-18*phycon.te70/phycon.te03/phycon.te01*
01625 phycon.te003;
01626 }
01627 else if( te <= 1.4e4 )
01628 {
01629
01630 hmirat_v = 8.204e-17*phycon.sqrte/phycon.te10/phycon.te01*
01631 phycon.te003;
01632 }
01633 else
01634 {
01635
01636 hmirat_v = 5.44e-16*phycon.te20/phycon.te01;
01637 }
01638
01639 return( hmirat_v );
01640 }
01641 STATIC void mole_h2_grain_form(void);
01642
01644 STATIC void mole_h_reactions(void);
01645
01646
01647 namespace {
01648 class mole_reaction_gamheh : public mole_reaction
01649 {
01650 typedef mole_reaction_gamheh T;
01651 public:
01652 virtual T* Create() const {return new T;}
01653 virtual const char* name() {return "gamheh";}
01654 double rk() const
01655 {
01656 double retval = 0.;
01657 long int limit,
01658 i;
01659
01660
01661
01662
01663
01664 limit = MIN2(hmi.iheh2-1 , rfield.nflux );
01665 for( i=hmi.iheh1-1; i < limit; i++ )
01666 {
01667 retval += rfield.flux[0][i] + rfield.ConInterOut[i]+ rfield.outlin[0][i] + rfield.outlin_noplot[i];
01668 }
01669 retval *= 4e-18;
01670
01671
01672 retval += 3.*iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].gamnc;
01673
01674 return retval;
01675 }
01676 };
01677
01678 enum {exclude, base, umisthack, federman, lithium, deuterium, ti, misc};
01679 static int source;
01680
01681
01682 static bool lgReactInitialized = false;
01683 }
01684
01685 void mole_create_react( void )
01686 {
01687
01688
01689
01690
01691
01692 if( lgReactInitialized )
01693 return;
01694
01695 lgReactInitialized = true;
01696
01697 DEBUG_ENTRY("mole_create_react()");
01698
01699 if (UDFA)
01700 read_data("rate05_feb17.csv",parse_udfa);
01701
01702
01703
01704 newfunc<mole_reaction_null>();
01705 newfunc<mole_reaction_hmrate>();
01706 newfunc<mole_reaction_hmrate_exo>();
01707 newfunc<mole_reaction_constrate>();
01708 newfunc<mole_reaction_th85rate>();
01709 newfunc<mole_reaction_crnurate>();
01710 newfunc<mole_reaction_crnurate_noalbedo>();
01711 newfunc<mole_reaction_co_lnu_c_o_lnu>();
01712 newfunc<mole_reaction_vib_evap>();
01713 newfunc<mole_reaction_th85rate_co>();
01714 newfunc<mole_reaction_grn_abs>();
01715
01716 newfunc<mole_reaction_grn_react>();
01717
01718 newfunc<mole_reaction_grn_photo>();
01719 newfunc<mole_reaction_oh_c2h2_co_ch3>();
01720 newfunc<mole_reaction_h_hnc_hcn_h>();
01721
01722 newfunc<mole_reaction_gamheh>();
01723 newfunc<mole_reaction_hd_photodissoc>();
01724 newfunc<mole_reaction_h2gdissoc>();
01725 newfunc<mole_reaction_h2sdissoc>();
01726 newfunc<mole_reaction_h2gexcit>();
01727 newfunc<mole_reaction_h2_collh_excit>();
01728 newfunc<mole_reaction_h2_collh2_excit>();
01729 newfunc<mole_reaction_h2_collh_deexc>();
01730 newfunc<mole_reaction_h2_collh2_deexc>();
01731 newfunc<mole_reaction_h2s_sp_decay>();
01732 newfunc<mole_reaction_hlike_phot>();
01733 newfunc<mole_reaction_gamtwo>();
01734 newfunc<mole_reaction_h2ph3p>();
01735 newfunc<mole_reaction_radath>();
01736 newfunc<mole_reaction_cryden_ov_bg>();
01737 newfunc<mole_reaction_h2scrphh>();
01738 newfunc<mole_reaction_h2crphh>();
01739 newfunc<mole_reaction_h2photon>();
01740 newfunc<mole_reaction_h2crphot>();
01741 newfunc<mole_reaction_rh2g_dis_h>();
01742 newfunc<mole_reaction_rh2s_dis_h>();
01743 newfunc<mole_reaction_rh2g_dis_h2>();
01744 newfunc<mole_reaction_rh2s_dis_h2>();
01745 newfunc<mole_reaction_rh2s_dis_h2_nodeexcit>();
01746 newfunc<mole_reaction_bh2g_dis_h>();
01747 newfunc<mole_reaction_bh2s_dis_h>();
01748 newfunc<mole_reaction_bh2g_dis_h2>();
01749 newfunc<mole_reaction_bh2s_dis_h2>();
01750 newfunc<mole_reaction_radasc>();
01751 newfunc<mole_reaction_assoc_ion>();
01752 newfunc<mole_reaction_grnh2>();
01753 newfunc<mole_reaction_grnh2s>();
01754 newfunc<mole_reaction_h2_spon_diss>();
01755 newfunc<mole_reaction_bhneut>();
01756 newfunc<mole_reaction_hneut>();
01757 newfunc<mole_reaction_asdbs>();
01758 newfunc<mole_reaction_asdbg>();
01759 newfunc<mole_reaction_asdfs>();
01760 newfunc<mole_reaction_asdfg>();
01761 newfunc<mole_reaction_c3bod>();
01762 newfunc<mole_reaction_cionhm>();
01763 newfunc<mole_reaction_hmphoto>();
01764 newfunc<mole_reaction_hmihph2p>();
01765 newfunc<mole_reaction_hmattach>();
01766 newfunc<mole_reaction_hpoexch>();
01767 newfunc<mole_reaction_hopexch>();
01768
01769 source = base;
01770 read_data("mole_co_base.dat",parse_base);
01771 if (mole_global.lgFederman)
01772 {
01773 source = federman;
01774 read_data("mole_co_federman.dat",parse_base);
01775 }
01776 if (!mole_global.lgLeidenHack)
01777 {
01778 source = umisthack;
01779 read_data("mole_co_umisthack.dat",parse_base);
01780 }
01781
01782 source = lithium;
01783 read_data("mole_lithium.dat",parse_base);
01784
01785 source = deuterium;
01786 read_data("mole_deuterium.dat",parse_base);
01787
01788 #if 0
01789 source = ti;
01790 read_data("mole_ti.dat",parse_base);
01791 #endif
01792
01793 source = misc;
01794 read_data("mole_misc.dat",parse_base);
01795
01796
01797 if (!mole_global.lgProtElim)
01798 {
01799 source = exclude;
01800 newreact("C+,OH=>CO,H+","hmrate",0.,0.,0.);
01801 }
01802
01803 newreact("H,H,grn=>H2,grn","grnh2",1.,0.,0.);
01804 newreact("H,H,grn=>H2*,grn","grnh2s",1.,0.,0.);
01805 newreact("H-,PHOTON=>H,e-","hmphoto",1.,0.,0.);
01806
01807
01808
01809
01810 fixit();
01811
01812 for(ChemAtomList::iterator atom=unresolved_atom_list.begin(); atom != unresolved_atom_list.end(); ++atom)
01813 {
01814 long nelem = (*atom)->el->Z-1;
01815 if( nelem >= ipHELIUM && dense.lgElmtOn[nelem] )
01816 {
01817 char react[32];
01818 sprintf(react,"H-,%s+=>H,%s", (*atom)->label().c_str(), (*atom)->label().c_str() );
01819 newreact(react,"hmrate",4e-6/sqrt(300.),-0.5,0.);
01820 }
01821 }
01822
01823 newreact("H,e-=>H-,PHOTON","hmattach",1.,0.,0.);
01824 newreact("H-,H+=>H2+,e-","hmihph2p",1.,0.,0.);
01825 newreact("H-,e-=>H,e-,e-","cionhm",1.,0.,0.);
01826 newreact("H,e-,e-=>H-,e-","c3bod",1.,0.,0.);
01827 newreact("H,H-=>H2,e-","asdfg",1.,0.,0.);
01828 newreact("H,H-=>H2*,e-","asdfs",1.,0.,0.);
01829 newreact("H2,e-=>H,H-","asdbg",1.,0.,0.);
01830 newreact("H2*,e-=>H,H-","asdbs",1.,0.,0.);
01831 newreact("H-,H+=>H,H","hneut",1.,0.,0.);
01832 newreact("H,H=>H-,H+","bhneut",1.,0.,0.);
01833 newreact("H2*=>H,H","h2_spon_diss",1.,0.,0.);
01834 if (!mole_global.lgLeidenHack)
01835 {
01836 newreact("H,H=>H2,PHOTON","radasc",3e-14,0.,0.);
01837 newreact("H,H=>H2+,e-","assoc_ion",3.27e-10,-0.35,17829.);
01838 newreact("H2,H=>H,H,H","rh2g_dis_h",1.,0.,0.);
01839 newreact("H2,H2=>H,H,H2","rh2g_dis_h2",1.,0.,0.);
01840 newreact("H,H,H=>H2,H","bh2g_dis_h",1.,0.,0.);
01841 newreact("H,H,H2=>H2,H2","bh2g_dis_h2",1.,0.,0.);
01842
01843 newreact("H2,PHOTON=>H2+,e-","h2photon",1.,0.,0.);
01844 newreact("H2,CRPHOT=>H2+,e-","h2crphot",1.,0.,0.);
01845 newreact("H2*,PHOTON=>H2+,e-","h2photon",1.,0.,0.);
01846 newreact("H2*,CRPHOT=>H2+,e-","h2crphot",1.,0.,0.);
01847 newreact("H2,CRPHOT=>H,H","h2crphh",1.,0.,0.);
01848 newreact("H2,CRPHOT=>H+,H-","cryden_ov_bg",3.9e-21,0.,0.);
01849 newreact("H2*,CRPHOT=>H+,H-","cryden_ov_bg",3.9e-21,0.,0.);
01850 newreact("H2,CRPHOT=>H+,H,e-","cryden_ov_bg",2.2e-19,0.,0.);
01851 newreact("H2*,CRPHOT=>H+,H,e-","cryden_ov_bg",2.2e-19,0.,0.);
01852 newreact("H+,H=>H2+,PHOTON","radath",1.,0.,0.);
01853 newreact("H2+,PHOTON=>H,H+","gamtwo",1.,0.,0.);
01854 newreact("H2+,CRPHOT=>H,H+","hlike_phot",1.,0.,0.);
01855
01856
01857 newreact("H3+,CRPHOT=>H2+,H+,e-","hlike_phot",2.,0.,0.);
01858 newreact("H,H+,e-=>H2+,e-","hmrate",2e-7 * SAHA * (4./(2.*1.)) * 3.634e-5 * pow(300.,-1.50),-1.50,0.);
01859 }
01860 else
01861 {
01862 newreact("H2,CRPHOT=>H2+,e-","hmrate",4.4e-17,0.,0.);
01863 newreact("H2,CRPHOT=>H,H+,e-","hmrate",1e-19,0.,0.);
01864 newreact("H2*,CRPHOT=>H,H+,e-","hmrate",1e-19,0.,0.);
01865 newreact("H2*,CRPHOT=>H2+,e-","hmrate",4.4e-17,0.,0.);
01866 newreact("H2,CRPHOT=>H,H","hmrate",5e-18,0.,0.);
01867 newreact("e-,H3+=>H2,H","hmrate",2.5e-8,-0.3,0.);
01868 newreact("e-,H3+=>H,H,H","hmrate",7.5e-8,-0.3,0.);
01869 newreact("H+,H=>H2+,PHOTON","hmrate",5.3e-19,1.85,0);
01870
01871
01872
01873 newreact("H2+,e-=>H,H","hmrate",1.6e-8,-0.43,0.);
01874 newreact("H2+,PHOTON=>H,H+","th85rate",5.7e-10,0.,1.9);
01875 }
01876
01877 newreact("H2*,CRPHOT=>H,H","h2scrphh",1.,0.,0.);
01878 newreact("H2,H2+=>H,H3+","h2ph3p",1.,0.,0.);
01879 newreact("H2*=>H2,PHOTON","h2s_sp_decay",1.,0.,0.);
01880 newreact("H2*,H2=>H2,H2","h2_collh2_deexc",1.,0.,0.);
01881 newreact("H2*,H=>H2,H","h2_collh_deexc",1.,0.,0.);
01882 newreact("H2,H2=>H2*,H2","h2_collh2_excit",1.,0.,0.);
01883 newreact("H2,H=>H2*,H","h2_collh_excit",1.,0.,0.);
01884
01885
01886 newreact("H2*,H=>H,H,H","rh2s_dis_h",1.,0.,0.);
01887 newreact("H2*,H2=>H2,H,H","rh2s_dis_h2",1.,0.,0.);
01888 newreact("H2*,H2*=>H2,H,H","rh2s_dis_h2",1.,0.,0.);
01889 newreact("H2*,H2*=>H2*,H,H","rh2s_dis_h2_nodeexcit",1.,0.,0.);
01890
01891 #if 0
01892
01893 newreact("H2,He=>He,H,H","rh2g_dis_h",1.,0.,0.);
01894 newreact("H2,H+=>H+,H,H","rh2g_dis_h",1.,0.,0.);
01895 newreact("H2,H3+=>H3+,H,H","rh2g_dis_h",1.,0.,0.);
01896 newreact("H2,e-=>e-,H,H","rh2g_dis_h",1.,0.,0.);
01897 newreact("H2*,He=>He,H,H","rh2s_dis_h",1.,0.,0.);
01898 newreact("H2*,H+=>H+,H,H","rh2s_dis_h",1.,0.,0.);
01899 newreact("H2*,H3+=>H3+,H,H","rh2s_dis_h",1.,0.,0.);
01900 newreact("H2*,e-=>e-,H,H","rh2s_dis_h",1.,0.,0.);
01901
01902
01903 newreact("He,H,H=>H2,He","bh2g_dis_h",1.,0.,0.);
01904 newreact("H+,H,H=>H2,H+","bh2g_dis_h",1.,0.,0.);
01905 newreact("H3+,H,H=>H2,H3+","bh2g_dis_h",1.,0.,0.);
01906 newreact("e-,H,H=>H2,e-","bh2g_dis_h",1.,0.,0.);
01907 newreact("H,H,H=>H2*,H","bh2s_dis_h",1.,0.,0.);
01908 newreact("H2,H,H=>H2*,H2","bh2s_dis_h2",1.,0.,0.);
01909 newreact("H2,H,H=>H2*,H2*","bh2s_dis_h2",1.,0.,0.);
01910 newreact("H2*,H,H=>H2*,H2*","bh2s_dis_h2",1.,0.,0.);
01911 newreact("He,H,H=>H2*,He","bh2s_dis_h",1.,0.,0.);
01912 newreact("H+,H,H=>H2*,H+","bh2s_dis_h",1.,0.,0.);
01913 newreact("H3+,H,H=>H2*,H3+","bh2s_dis_h",1.,0.,0.);
01914 newreact("e-,H,H=>H2*,e-","bh2s_dis_h",1.,0.,0.);
01915 #endif
01916
01917 newreact("H2,PHOTON=>H2*","h2gexcit",1.,0.,0.);
01918 newreact("H2*,PHOTON=>H,H","h2sdissoc",1.,0.,0.);
01919 newreact("H2,PHOTON=>H,H","h2gdissoc",1.,0.,0.);
01920 newreact("HeH+,PHOTON=>H,He+","gamheh",1.,0.,0.);
01921
01922 if(0)
01923 {
01924
01925 newreact("OHgrn,Hgrn=>H2Ogrn","grn_react",1.,0.,0.);
01926 }
01927
01928 if (UDFA)
01929 {
01930 fprintf(stderr,"Finished testing against UDFA database\n");
01931 exit(-1);
01932 }
01933
01934 if (0)
01935 plot_sparsity();
01936
01937 mole_check_reverse_reactions();
01938
01939 if( deut.lgElmtOn )
01940 mole_generate_isotopologue_reactions( "H", "D" );
01941
01942 long index = 0;
01943 for( mole_reaction_i it = mole_priv::reactab.begin(); it != mole_priv::reactab.end(); ++it, ++index )
01944 it->second->index = index;
01945
01946 mole.reaction_rks.resize( index );
01947 mole.old_zone = -1;
01948 memset( &mole.reaction_rks[0], 0, (unsigned long)(index)*sizeof(double) );
01949
01950
01951 for( mole_reaction_i it = mole_priv::reactab.begin(); it != mole_priv::reactab.end(); ++it )
01952 register_reaction_vectors( it->second );
01953 }
01954
01955 STATIC void mole_generate_isotopologue_reactions( string atom_old, string atom_new )
01956 {
01957 DEBUG_ENTRY("mole_generate_isotopologue_reactions()");
01958
01959 bool lgDebug = false;
01960
01961
01962 vector<string> newReactionStrings;
01963 vector<mole_reaction*> oldReactions;
01964
01965
01966 for( mole_reaction_ci it = mole_priv::reactab.begin(); it != mole_priv::reactab.end(); ++it )
01967 {
01968 bool lgFound = false;
01969
01970 for( long i=0; i<it->second->nreactants; ++i )
01971 {
01972
01973 for( nAtoms_i atom=it->second->reactants[i]->nAtom.begin(); atom != it->second->reactants[i]->nAtom.end(); ++atom )
01974 {
01975 if( atom->first->label() == atom_old )
01976 {
01977 lgFound = true;
01978 continue;
01979 }
01980 }
01981 if( lgFound )
01982 continue;
01983 }
01984
01985 if( !lgFound )
01986 continue;
01987
01988 if( lgDebug )
01989 fprintf( ioQQQ, "DEBUGGG mole_generate_isotopologue_reactions %s ..........\n", it->first.c_str() );
01990
01991 for( long i=0; i<it->second->nreactants; ++i )
01992 {
01993
01994 if( it->second->reactants[i]->mole_mass < 0.99 * ATOMIC_MASS_UNIT )
01995 continue;
01996
01997 vector<string> react_iso_labels;
01998 ChemAtomList atomsLeftToRight;
01999 vector< int > numAtoms;
02000 string embellishments;
02001
02002 bool lgParseOK = parse_species_label( it->second->reactants[i]->label.c_str(), atomsLeftToRight, numAtoms, embellishments );
02003 ASSERT( lgParseOK == true );
02004
02005 create_isotopologues_one(
02006 atomsLeftToRight,
02007 numAtoms,
02008 atom_old,
02009 atom_new,
02010 embellishments,
02011 react_iso_labels );
02012 for( unsigned j=0; j<react_iso_labels.size(); ++j )
02013 {
02014 for( long k=0; k<it->second->nproducts; ++k )
02015 {
02016
02017 if( it->second->products[k]->mole_mass < 0.99 * ATOMIC_MASS_UNIT )
02018 continue;
02019
02020 vector<string> prod_iso_labels;
02021 atomsLeftToRight.clear();
02022 numAtoms.clear();
02023 embellishments.clear();
02024
02025 lgParseOK = parse_species_label( it->second->products[k]->label.c_str(), atomsLeftToRight, numAtoms, embellishments );
02026 ASSERT( lgParseOK == true );
02027
02028 create_isotopologues_one(
02029 atomsLeftToRight,
02030 numAtoms,
02031 atom_old,
02032 atom_new,
02033 embellishments,
02034 prod_iso_labels );
02035 for( unsigned l=0; l<prod_iso_labels.size(); ++l )
02036 {
02037 string react_string;
02038
02039
02040 for( long i1=0; i1<i; ++i1 )
02041 {
02042 react_string += it->second->reactants[i1]->label;
02043 if( i1 != it->second->nreactants-1 )
02044 react_string += ",";
02045 }
02046 react_string += react_iso_labels[j];
02047 if( i != it->second->nreactants-1 )
02048 react_string += ",";
02049 for( long i2=i+1; i2<it->second->nreactants; ++i2 )
02050 {
02051 react_string += it->second->reactants[i2]->label;
02052 if( i2 != it->second->nreactants-1 )
02053 react_string += ",";
02054 }
02055
02056 react_string += "=>";
02057
02058 for( long k1=0; k1<k; ++k1 )
02059 {
02060 react_string += it->second->products[k1]->label;
02061 if( k1 != it->second->nproducts-1 )
02062 react_string += ",";
02063 }
02064 react_string += prod_iso_labels[l];
02065 if( k != it->second->nproducts-1 )
02066 react_string += ",";
02067 for( long k2=k+1; k2<it->second->nproducts; ++k2 )
02068 {
02069 react_string += it->second->products[k2]->label;
02070 if( k2 != it->second->nproducts-1 )
02071 react_string += ",";
02072 }
02073
02074 if( lgDebug )
02075 fprintf( ioQQQ, "DEBUGGG mole_generate_isotopologue_reactions .................... %s\n",
02076 react_string.c_str() );
02077
02078 string canon_react_string = canonicalize_reaction_label( react_string.c_str() );
02079
02080 newReactionStrings.push_back( canon_react_string );
02081 oldReactions.push_back( it->second.get_ptr() );
02082 }
02083 }
02084 }
02085 }
02086 }
02087
02088 ASSERT( oldReactions.size() == newReactionStrings.size() );
02089
02090
02091 vector<mole_reaction*>::const_iterator it2 = oldReactions.begin();
02092 for( vector<string>::const_iterator it1 = newReactionStrings.begin(); it1 != newReactionStrings.end(); ++it1, ++it2 )
02093 {
02094 fixit();
02095
02096
02097 if( mole_priv::reactab.find( it1->c_str() ) == mole_priv::reactab.end() )
02098 newreact( it1->c_str(), (*it2)->name(), (*it2)->a, (*it2)->b, (*it2)->c );
02099 }
02100
02101 return;
02102 }
02103
02104 STATIC void mole_check_reverse_reactions(void)
02105 {
02106 DEBUG_ENTRY("mole_check_reverse_reactions()");
02107
02108 char chLabel[50], chLabelSave[50];
02109 int exists;
02110
02111 for(mole_reaction_i p=mole_priv::reactab.begin();
02112 p != mole_priv::reactab.end(); ++p)
02113 {
02114 mole_reaction_i q = p;
02115 strcpy( chLabel, p->second->label.c_str() );
02116 strcpy( chLabelSave, p->second->label.c_str() );
02117 char *str = chLabel;
02118 const char *delim = "=>";
02119 char *chNewProducts = strtok( str, delim );
02120 char *chNewReactants = strtok( NULL, delim );
02121 char chNewReaction[50];
02122
02123 strcpy( chNewReaction, chNewReactants );
02124 strcat( chNewReaction, "=>" );
02125 strcat( chNewReaction, chNewProducts );
02126
02127
02128 q = mole_priv::reactab.find(chNewReaction);
02129
02130 exists = (q != mole_priv::reactab.end());
02131 if ( !exists )
02132 {
02133 if( trace.lgTraceMole )
02134 {
02135 fprintf(ioQQQ,"Warning! No reverse reaction for %30s. ", p->second->label.c_str() );
02136 fprintf( ioQQQ,"\n" );
02137 }
02138
02139 fixit();
02140
02141 }
02142 }
02143
02144 return;
02145 }
02146
02147 STATIC double mole_get_equilibrium_constant( const char buf[] )
02148 {
02149 DEBUG_ENTRY("mole_get_equilibrium_constant()");
02150
02151 mole_reaction *rate = mole_findrate_s(buf);
02152
02153 if( !rate )
02154 return 0;
02155
02156
02157 double result = 1.;
02158 for( long i=0; i<rate->nreactants; ++i)
02159 result *= mole_partition_function(rate->reactants[i]);
02160 for( long i=0; i<rate->nproducts; ++i)
02161 result /= mole_partition_function(rate->products[i]);
02162
02163
02164 return result;
02165 }
02166
02167 STATIC double mole_partition_function( const molecule* const sp)
02168 {
02169 DEBUG_ENTRY("mole_partition_function()");
02170
02171 if( sp->label == "PHOTON" || sp->label == "CRPHOT" )
02172 {
02173 fixit();
02174 fixit();
02175 return 1.;
02176 }
02177 else if( sp->label == "CRP" || sp->label == "grn" )
02178 return 1.;
02179
02180 fixit();
02181 double q = 1.;
02182
02183 double deltaH = sp->form_enthalpy * (1e10/AVOGADRO/BOLTZMANN);
02184 ASSERT( sp->mole_mass > 0. );
02185 double part_fun = pow(phycon.te*sp->mole_mass/(HION_LTE_POP*ELECTRON_MASS), 1.5) * q * dsexp(deltaH/phycon.te);
02186 ASSERT( part_fun > 1e-100 );
02187 ASSERT( part_fun < BIGFLOAT );
02188
02189 return part_fun;
02190 }
02191
02192 void mole_cmp_num_in_out_reactions()
02193 {
02194 DEBUG_ENTRY("mole_cmp_num_in_out_reactions()");
02195
02196 vector<long> numForm, numDest;
02197 numForm.resize( mole_global.num_total );
02198 numDest.resize( mole_global.num_total );
02199
02200 for(mole_reaction_i p=mole_priv::reactab.begin(); p != mole_priv::reactab.end(); p++)
02201 {
02202 count_ptr<mole_reaction> rate = p->second;
02203 for( long i=0; i<rate->nreactants; ++i)
02204 {
02205 ++numDest[ rate->reactants[i]->index ];
02206 }
02207
02208 for( long i=0; i<rate->nproducts; ++i)
02209 {
02210 ++numForm[ rate->products[i]->index ];
02211 }
02212 }
02213
02214 for( unsigned i=0; i<numForm.size(); ++i )
02215 {
02216 if( numForm[i]==0 && numDest[i]==0 )
02217 continue;
02218 if( numForm[i]>1 && numDest[i]>1 )
02219 continue;
02220 if( mole_global.list[i]->isMonatomic() )
02221 continue;
02222 fprintf( ioQQQ, "DEBUGGG mole_cmp_num_in_out_reactions %*s: in %4li out %4li\n", CHARS_SPECIES, mole_global.list[i]->label.c_str(), numForm[i], numDest[i] );
02223 }
02224
02225 return;
02226 }
02227
02228 STATIC char *getcsvfield(char **s,char c);
02229 STATIC void parse_base(char *s)
02230 {
02231 char *label, *reactstr, *f;
02232 double a, b, c;
02233 label = getcsvfield(&s,':');
02234 reactstr = getcsvfield(&s,':');
02235 f = getcsvfield(&s,':');
02236 a = atof(f);
02237 f = getcsvfield(&s,':');
02238 b = atof(f);
02239 f = getcsvfield(&s,':');
02240 c = atof(f);
02241
02242 newreact(label,reactstr,a,b,c);
02243
02244 }
02245
02246 STATIC void newreact(const char label[], const char fun[], double a, double b, double c)
02247 {
02248 DEBUG_ENTRY("newreact()");
02249
02250 ratefunc rate = findfunc(fun);
02251 if (rate.get_ptr() == NULL)
02252 {
02253 fprintf(stderr,"Rate function %s not known for reaction %s. Aborting. Sorry.\n",fun,label);
02254 cdEXIT( EXIT_FAILURE );
02255 }
02256
02257 rate->label = label;
02258 rate->a = a;
02259 rate->b = b;
02260 rate->c = c;
02261
02262 rate->source = source;
02263
02264 rate->photon = 0;
02265
02266 if( parse_reaction( rate, label ) == 0 )
02267 return;
02268
02269 canonicalize_reaction( rate );
02270
02271 const char *rateLabelPtr = rate->label.c_str();
02272
02273 ASSERT(lgReactBalance(rate));
02274
02275 rate->udfastate = ABSENT;
02276 if (UDFA)
02277 {
02278 compare_udfa(rate);
02279 if (rate->udfastate == ABSENT)
02280 {
02281 fprintf(stderr,"Reaction %s not in UDFA\n",rateLabelPtr);
02282 }
02283 }
02284
02285
02286 if (rate->nreactants == 2 && rate->reactants[0]->mole_mass!=0.0 && rate->reactants[1]->mole_mass!=0.0 )
02287 {
02288 rate->reduced_mass = 1./(1./rate->reactants[0]->mole_mass+1./rate->reactants[1]->mole_mass);
02289 }
02290 else
02291 {
02292 rate->reduced_mass = 0.;
02293 }
02294
02295
02296 mole_reaction_i p = mole_priv::reactab.find(rateLabelPtr);
02297 int exists = (p != mole_priv::reactab.end());
02298
02299 if(exists && !t_version::Inst().lgRelease && t_version::Inst().nBetaVer == 0)
02300 {
02301
02302 fprintf(ioQQQ,"Warning: duplicate reaction %s -- using new version\n",rateLabelPtr);
02303 }
02304 mole_priv::reactab[rateLabelPtr] = rate;
02305 }
02306
02307 STATIC long parse_reaction( count_ptr<mole_reaction>& rate, const char label[] )
02308 {
02309 DEBUG_ENTRY("parse_reaction()");
02310
02311 for (int i=0;i<MAXREACTANTS;i++)
02312 {
02313 rate->reactants[i] = NULL;
02314 }
02315 rate->nreactants = 0;
02316
02317 for (int i=0;i<MAXPRODUCTS;i++)
02318 {
02319 rate->products[i] = NULL;
02320 }
02321 rate->nproducts = 0;
02322
02323 bool lgProd = false;
02324 string buf = "";
02325 for(int i=0;!i || label[i-1]!='\0';i++)
02326 {
02327 if(label[i] == ',' || label[i] == '=' || label[i] == '\0')
02328 {
02329 molecule *sp = findspecies(buf.c_str());
02330 if( sp == null_mole || !sp->isEnabled )
02331 {
02332 if( trace.lgTraceMole )
02333 fprintf(ioQQQ,"Mole_reactions: ignoring reaction %s (species %s not active)\n",label,buf.c_str());
02334 return 0;
02335 }
02336 buf = "";
02337 if(! lgProd)
02338 {
02339 if (rate->nreactants >= MAXREACTANTS)
02340 {
02341 fprintf(stderr,"Mole_reactions: Too many reactants in %s, only %d allowed\n",label,MAXREACTANTS);
02342 exit(-1);
02343 }
02344 rate->reactants[rate->nreactants] = sp;
02345 rate->nreactants++;
02346 }
02347 else
02348 {
02349 if (rate->nproducts >= MAXPRODUCTS)
02350 {
02351 fprintf(stderr,"Mole_reactions: Too many products in %s, only %d allowed\n",label,MAXPRODUCTS);
02352 exit(-1);
02353 }
02354 rate->products[rate->nproducts] = sp;
02355 rate->nproducts++;
02356 }
02357 if(label[i] == '=')
02358 {
02359 i++;
02360 if (label[i] != '>')
02361 {
02362 fprintf(ioQQQ,"Format error in reaction %s\n",label);
02363 cdEXIT(EXIT_FAILURE);
02364 }
02365 lgProd = true;
02366 }
02367 }
02368 else
02369 {
02370 buf += label[i];
02371 }
02372 }
02373
02374 ASSERT( rate->nreactants );
02375 ASSERT( rate->nproducts );
02376
02377 return 1;
02378 }
02379
02380 STATIC string canonicalize_reaction_label( const char label[] )
02381 {
02382 DEBUG_ENTRY("canonicalize_reaction_label()");
02383
02384
02385 ratefunc rate = findfunc("null");
02386 rate->label = label;
02387 parse_reaction( rate, label );
02388 canonicalize_reaction( rate );
02389
02390
02391
02392
02393 return rate->label;
02394 }
02395
02396 STATIC void canonicalize_reaction( count_ptr<mole_reaction>& rate )
02397 {
02398 DEBUG_ENTRY("canonicalize_reaction()");
02399
02400
02401
02402
02403 t_mole_global::sort(rate->reactants,rate->reactants+rate->nreactants);
02404 t_mole_global::sort(rate->products,rate->products+rate->nproducts);
02405
02406
02407 string newLabel;
02408 for( long i=0; i<rate->nreactants; ++i )
02409 {
02410 newLabel += rate->reactants[i]->label;
02411 if( i != rate->nreactants-1 )
02412 newLabel += ",";
02413 }
02414 newLabel += "=>";
02415 for( long i=0; i<rate->nproducts; ++i )
02416 {
02417 newLabel += rate->products[i]->label;
02418 if( i != rate->nproducts-1 )
02419 newLabel += ",";
02420 }
02421
02422
02423 rate->label = newLabel;
02424
02425 return;
02426 }
02427
02428 STATIC void register_reaction_vectors( count_ptr<mole_reaction> rate )
02429 {
02430 DEBUG_ENTRY("register_reaction_vectors()");
02431
02432 for (long k=0;k<rate->nreactants;k++)
02433 {
02434 rate->rvector[k] = NULL;
02435 rate->rvector_excit[k] = NULL;
02436 }
02437
02438 for (long k=0;k<rate->nproducts;k++)
02439 {
02440 rate->pvector[k] = NULL;
02441 rate->pvector_excit[k] = NULL;
02442 }
02443
02444
02445 for (long i=0;i<rate->nproducts;i++)
02446 {
02447 if (rate->pvector[i] == NULL)
02448 {
02449 for (long k=0;k<rate->nreactants;k++)
02450 {
02451 if (rate->rvector[k] == NULL)
02452 {
02453 if (rate->products[i] == rate->reactants[k])
02454 {
02455 rate->rvector[k] = rate->products[i];
02456 rate->pvector[i] = rate->reactants[k];
02457 break;
02458 }
02459 }
02460 }
02461 }
02462 }
02463
02464
02465 for (long i=0;i<rate->nproducts;i++)
02466 {
02467 if (rate->pvector[i] == NULL)
02468 {
02469 for (long k=0;k<rate->nreactants;k++)
02470 {
02471 if (rate->rvector[k] == NULL)
02472 {
02473 if (rate->products[i]->groupnum != -1 &&
02474 rate->products[i]->groupnum ==
02475 rate->reactants[k]->groupnum)
02476 {
02477 rate->rvector[k] = rate->products[i];
02478 rate->pvector[i] = rate->reactants[k];
02479 break;
02480 }
02481 }
02482 }
02483 }
02484 }
02485
02486
02487 for (long i=0;i<rate->nproducts;i++)
02488 {
02489 if (rate->pvector[i] == NULL && rate->pvector_excit[i] == NULL)
02490 {
02491 for (long k=0;k<rate->nreactants;k++)
02492 {
02493 if (rate->rvector[k] == NULL && rate->rvector_excit[k] == NULL )
02494 {
02495 if ( lgDifferByExcitation( *rate->products[i], *rate->reactants[k] ) )
02496 {
02497 rate->rvector_excit[k] = rate->products[i];
02498 rate->pvector_excit[i] = rate->reactants[k];
02499 break;
02500 }
02501 }
02502 }
02503 }
02504 }
02505
02506 return;
02507 }
02508
02509
02510 STATIC void plot_sparsity(void)
02511 {
02512 FILE *sparsefp;
02513 int i, j, nb, ch;
02514 long int ratej;
02515 double **c;
02516
02517 c = (double **)MALLOC((size_t)mole_global.num_total*sizeof(double *));
02518 c[0] = (double *)MALLOC((size_t)mole_global.num_total*
02519 mole_global.num_total*sizeof(double));
02520
02521 for(i=1;i<mole_global.num_total;i++)
02522 {
02523 c[i] = c[i-1]+mole_global.num_total;
02524 }
02525
02526 for ( j=0; j < mole_global.num_total; j++ )
02527 {
02528 for ( i=0; i < mole_global.num_total; i++ )
02529 {
02530 c[j][i] = 0.;
02531 }
02532 }
02533
02534 for(mole_reaction_i p=mole_priv::reactab.begin();
02535 p != mole_priv::reactab.end(); ++p)
02536 {
02537 mole_reaction &rate = *p->second;
02538
02539 for (j=0;j<rate.nreactants;j++)
02540 {
02541 ratej = rate.reactants[j]->index;
02542 for (i=0;i<rate.nreactants;i++)
02543 {
02544 if (rate.rvector[i] == NULL)
02545 c[ratej][rate.reactants[i]->index] = 1.0;
02546 }
02547 for (i=0;i<rate.nproducts;i++)
02548 {
02549 if (rate.pvector[i] == NULL)
02550 c[ratej][rate.products[i]->index] = 1.0;
02551 }
02552 }
02553 }
02554
02555 sparsefp = fopen("sparse.pbm","w");
02556 fprintf(sparsefp,"P4\n%d %d\n",
02557 mole_global.num_total,mole_global.num_total);
02558
02559 for ( j=0; j < mole_global.num_total; j++ )
02560 {
02561 nb = ch = 0;
02562 for ( i=0; i < mole_global.num_total; i++ )
02563 {
02564 ch = (ch << 1) | (c[i][j] != 0.0);
02565 nb++;
02566 if (nb == 8)
02567 {
02568 fputc(ch,sparsefp);
02569 nb = ch = 0;
02570 }
02571 }
02572 if (nb != 0)
02573 {
02574 ch <<= 8-nb;
02575 fputc(ch,sparsefp);
02576 }
02577 }
02578 fclose(sparsefp);
02579 free(c[0]);
02580 free(c);
02581 }
02582
02583 STATIC bool lgReactBalance(const count_ptr<mole_reaction> &rate)
02584 {
02585 molecule::nAtomsMap nel;
02586 int dcharge, n, sign;
02587 bool lgOK = true;
02588
02589 dcharge = 0;
02590 for (n=0;n<rate->nreactants;n++)
02591 {
02592 for( nAtoms_i it = rate->reactants[n]->nAtom.begin(); it != rate->reactants[n]->nAtom.end(); ++it )
02593 nel[it->first] += it->second;
02594 dcharge += rate->reactants[n]->charge;
02595 }
02596 for (n=0;n<rate->nproducts;n++)
02597 {
02598 for( nAtoms_i it = rate->products[n]->nAtom.begin(); it != rate->products[n]->nAtom.end(); ++it )
02599 nel[it->first] -= it->second;
02600 dcharge -= rate->products[n]->charge;
02601 }
02602 if (dcharge != 0)
02603 {
02604 fprintf(stderr,"Reaction %s charge out of balance by %d\n",
02605 rate->label.c_str(),dcharge);
02606 lgOK = false;
02607 }
02608
02609 for( nAtoms_i it = nel.begin(); it != nel.end(); ++it )
02610 {
02611 if(it->second != 0)
02612 {
02613 if(it->second > 0)
02614 sign = 1;
02615 else
02616 sign = -1;
02617 fprintf(stderr,"Error: reaction %s %s %d of element %s\n",
02618 rate->label.c_str(),sign==1?"destroys":"creates",
02619 sign*it->second,
02620 it->first->label().c_str() );
02621 lgOK = false;
02622 }
02623 }
02624 return lgOK;
02625 }
02626
02627 #ifndef ZLIB_H
02628 #define gzopen(file,mode) fopen(file,mode)
02629 #define gzgets(fp,buf,siz) fgets(buf,siz,fp)
02630 #define gzclose(fp) fclose(fp)
02631 #define gzFile FILE
02632 #endif
02633
02634 enum {BUFSIZE=256};
02635
02636 namespace
02637 {
02638 class formula_species {
02639 public:
02640 molecule *reactants[MAXREACTANTS], *products[MAXPRODUCTS];
02641 };
02642
02643 bool operator< (const formula_species &a, const formula_species &b)
02644 {
02645 int i;
02646 for (i=0;i<MAXREACTANTS;i++)
02647 {
02648 if (a.reactants[i]<b.reactants[i])
02649 return true;
02650 else if (a.reactants[i]>b.reactants[i])
02651 return false;
02652 }
02653 for (i=0;i<MAXPRODUCTS;i++)
02654 {
02655 if (a.products[i]<b.products[i])
02656 return true;
02657 else if (a.products[i]>b.products[i])
02658 return false;
02659 }
02660 return false;
02661 }
02662
02663 struct udfa_reaction_s {
02664 int index;
02665 formula_species l;
02666 char source;
02667 double a, b, c, trange[2];
02668 };
02669 }
02670
02671 static map <formula_species,struct udfa_reaction_s *> udfatab;
02672
02673 STATIC void read_data(const char file[], void (*parse)(char *s))
02674 {
02675 FILE *fp;
02676 char buf[BUFSIZE];
02677
02678 fp = open_data(file,"r");
02679 if (!fp)
02680 {
02681 fprintf(stderr,"Error, could not read %s\n",file);
02682 exit(-1);
02683 }
02684
02685 fixit();
02686 while(fgets(buf,BUFSIZE,fp))
02687 {
02688 if( buf[0] == '#' )
02689 continue;
02690 parse(buf);
02691 }
02692 fclose(fp);
02693 }
02694 #define FLTEQ(a,b) (fabs((a)-(b)) <= 1e-6*fabs((a)+(b)))
02695 STATIC void parse_udfa(char *s)
02696 {
02697 char *f;
02698 struct udfa_reaction_s *r;
02699 map <formula_species, struct udfa_reaction_s *>::iterator p;
02700 unsigned int havespecies = 1, i, n;
02701
02702
02703 int lgCRPHOT=0;
02704
02705 r = (struct udfa_reaction_s *)MALLOC(sizeof(struct udfa_reaction_s));
02706 f = getcsvfield(&s,',');
02707 r->index = atoi(f);
02708
02709
02710 for (n=0;n<MAXREACTANTS;n++)
02711 {
02712 r->l.reactants[n] = NULL;
02713 }
02714 i = 0;
02715 for (n=0;n<MAXREACTANTS;n++)
02716 {
02717 f = getcsvfield(&s,',');
02718 if (f[0] != '\0')
02719 {
02720 i++;
02721 r->l.reactants[n] = findspecies(f);
02722 if (r->l.reactants[n] == null_mole)
02723 havespecies = 0;
02724 if (!strncmp(f,"CRPHOT",6))
02725 lgCRPHOT = 1;
02726 }
02727 }
02728 t_mole_global::sort(r->l.reactants,r->l.reactants+i);
02729
02730
02731 for (n=0;n<MAXPRODUCTS;n++)
02732 {
02733 r->l.products[n] = NULL;
02734 }
02735 i = 0;
02736 for (n=0;n<MAXPRODUCTS;n++)
02737 {
02738 f = getcsvfield(&s,',');
02739 if (f[0] != '\0')
02740 {
02741 i++;
02742 r->l.products[n] = findspecies(f);
02743 if (r->l.products[n] == null_mole)
02744 havespecies = 0;
02745 }
02746 }
02747
02748 t_mole_global::sort(r->l.products,r->l.products+i);
02749
02750
02751 f = getcsvfield(&s,',');
02752 r->a = atof(f);
02753 f = getcsvfield(&s,',');
02754 r->b = atof(f);
02755 f = getcsvfield(&s,',');
02756 r->c = atof(f);
02757
02758 if (lgCRPHOT)
02759 {
02760
02761
02762
02763
02764 ASSERT (FLTEQ(r->a,1.3e-17));
02765 r->a = r->c;
02766 r->c = 0.;
02767 }
02768
02769
02770 f = getcsvfield(&s,',');
02771 r->source = f[0]?f[0]:'?';
02772 for (n=0;n<2;n++) {
02773 f = getcsvfield(&s,',');
02774 r->trange[n] = atof(f);
02775 }
02776
02777 if (!havespecies)
02778 {
02779 free(r);
02780 }
02781 else
02782 {
02783 p = udfatab.find(r->l);
02784 if (p != udfatab.end())
02785 {
02786 fprintf(stderr,"Duplicate reaction\n");
02787 }
02788 udfatab[r->l] = r;
02789 }
02790 }
02791 STATIC char *getcsvfield(char **s, char c)
02792 {
02793 char *sv, *f;
02794
02795 sv = strchr(*s,c);
02796 if (sv) {
02797 *sv++ = '\0';
02798 }
02799 f = *s;
02800 *s = sv;
02801 return f;
02802 }
02803 STATIC void compare_udfa(const count_ptr<mole_reaction> &rate)
02804 {
02805 formula_species s;
02806 unsigned int n;
02807 map<formula_species, struct udfa_reaction_s *>::iterator p;
02808 struct udfa_reaction_s *u;
02809
02810 for (n=0;n<MAXREACTANTS;n++)
02811 {
02812 s.reactants[n] = rate->reactants[n];
02813 }
02814 for (n=0;n<MAXPRODUCTS;n++)
02815 {
02816 s.products[n] = rate->products[n];
02817 }
02818 p = udfatab.find(s);
02819 if (p == udfatab.end() )
02820 {
02821
02822 return;
02823 }
02824 else
02825 {
02826 u = p->second;
02827 if (FLTEQ(rate->a,u->a) && FLTEQ(rate->b,u->b) && FLTEQ(rate->c,u->c))
02828 {
02829 rate->udfastate = CORRECT;
02830
02831 }
02832 else
02833 {
02834 rate->udfastate = CONFLICT;
02835
02836
02837 }
02838 }
02839 }
02840
02841 namespace
02842 {
02843 ratefunc findfunc (const char name[])
02844 {
02845 return count_ptr<mole_reaction>(mole_priv::functab[name]->Create());
02846 }
02847 }
02848
02849 void mole_update_rks(void)
02850 {
02851 enum { DEBUG_MOLE = false };
02852
02853 DEBUG_ENTRY("mole_update_rks()");
02854
02855 mole_h2_grain_form();
02856
02857 mole_h_reactions();
02858
02859 for (mole_reaction_i p
02860 =mole_priv::reactab.begin(); p != mole_priv::reactab.end(); ++p)
02861 {
02862 mole_reaction &rate = *p->second;
02863 long index = rate.index;
02864 realnum oldrk = (realnum)mole.reaction_rks[index];
02865 realnum newrk = rate.a*rate.rk();
02866 mole.reaction_rks[index] = newrk;
02867 if (DEBUG_MOLE)
02868 {
02869 if (fabs(newrk-oldrk) > 0.1*newrk)
02870 fprintf(ioQQQ,"%s: %15.8g => %15.8g\n",
02871 rate.label.c_str(),oldrk,newrk);
02872 }
02873 }
02874 }
02875 void mole_rk_bigchange(void)
02876 {
02877 enum { DEBUG_MOLE = false };
02878
02879 DEBUG_ENTRY("mole_rk_bigchange()");
02880
02881 if ( mole.old_reaction_rks.size() == 0 )
02882 {
02883 mole.old_zone = -1;
02884 mole.old_reaction_rks.resize(mole.reaction_rks.size());
02885 }
02886
02887 if (nzone > 1)
02888 {
02889 ASSERT(mole.old_zone == nzone - 1);
02890 double bigchange = 0.;
02891 unsigned long bigindex = ULONG_MAX;
02892 for (unsigned long index=0; index<mole.reaction_rks.size(); ++index)
02893 {
02894 double oldrk = mole.old_reaction_rks[index],
02895 newrk = mole.reaction_rks[index],
02896 sum = oldrk+newrk, diff = newrk-oldrk;
02897 if (sum > 0.)
02898 {
02899 double change = fabs(diff)/sum;
02900 if (change > bigchange)
02901 {
02902 bigchange = change;
02903 bigindex = index;
02904 }
02905 }
02906 }
02907
02908 for (mole_reaction_i p
02909 =mole_priv::reactab.begin(); p != mole_priv::reactab.end(); ++p)
02910 {
02911 mole_reaction &rate = *p->second;
02912 if (bigindex == (unsigned) rate.index)
02913 {
02914 double oldrk = mole.old_reaction_rks[bigindex],
02915 newrk = mole.reaction_rks[bigindex];
02916 double pct = 0.;
02917 if (oldrk > 0.)
02918 pct = 100.*(newrk-oldrk)/oldrk;
02919 fprintf(ioQQQ, "Zone %ld, big chemistry rate change %s:"
02920 " %15.8g => %15.8g (%6.2g%%)\n",
02921 nzone,rate.label.c_str(),oldrk,newrk,pct);
02922 break;
02923 }
02924 }
02925 }
02926
02927 mole.old_zone = nzone;
02928 for (unsigned long index=0; index<mole.reaction_rks.size(); ++index)
02929 {
02930 mole.old_reaction_rks[index] = mole.reaction_rks[index];
02931 }
02932 }
02933
02934 STATIC void mole_h2_grain_form(void)
02935 {
02936 double T_ortho_para_crit;
02937 realnum AveVelH = GetAveVelocity( dense.AtomicWeight[ipHYDROGEN] );
02938 realnum AveVelH2 = GetAveVelocity( 2.f * dense.AtomicWeight[ipHYDROGEN] );
02939
02940
02941
02942
02943 if( gv.lgDustOn() )
02944 {
02945
02946 # ifndef IGNORE_QUANTUM_HEATING
02947
02948 GrainDrive();
02949 # endif
02950
02951
02952
02953
02954 h2.rate_grain_op_conserve = 0.;
02955
02956 h2.rate_grain_J1_to_J0 = 0.;
02957
02958
02959 for( size_t nd=0; nd < gv.bin.size(); nd++ )
02960 {
02961 # ifndef IGNORE_QUANTUM_HEATING
02962 long k, qnbin;
02963 double *qtemp, *qprob;
02964 bool lgUseQHeat = gv.lgGrainPhysicsOn && gv.bin[nd]->lgQHeat;
02965 # endif
02966
02967
02968
02969
02970
02971
02972
02973 double sticking_probability_H = 1./(1. + 0.04*sqrt(gv.bin[nd]->tedust+phycon.te) +
02974 0.002*phycon.te + 8e-6*phycon.te*phycon.te);
02975
02976 # ifndef IGNORE_QUANTUM_HEATING
02977
02978 if( lgUseQHeat )
02979 {
02980 qtemp = (double*)MALLOC((size_t)(NQGRID*sizeof(double)));
02981 qprob = (double*)MALLOC((size_t)(NQGRID*sizeof(double)));
02982
02983 qheat(qtemp,qprob,&qnbin,nd);
02984
02985 if( gv.bin[nd]->lgUseQHeat )
02986 {
02987 ASSERT( qnbin > 0 );
02988 }
02989 else
02990 {
02991 qnbin = 1;
02992 qprob[0] = 1.;
02993 qtemp[0] = gv.bin[nd]->tedust;
02994 }
02995
02996 gv.bin[nd]->rate_h2_form_grains_HM79 = 0.;
02997
02998 for( k=0; k < qnbin; k++ )
02999 {
03000
03001
03002
03003
03004
03005
03006
03007
03008 double conversion_efficiency_HM79 = 1/(1. + 1e4*sexp(920./qtemp[k]));
03009 sticking_probability_H = 1./(1. + 0.04*sqrt(qtemp[k]+phycon.te) +
03010 0.002*phycon.te + 8e-6*phycon.te*phycon.te);
03011
03012 gv.bin[nd]->rate_h2_form_grains_HM79 += qprob[k] * sticking_probability_H *
03013 conversion_efficiency_HM79;
03014 }
03015
03016
03017
03018
03019 gv.bin[nd]->rate_h2_form_grains_HM79 *= 0.5 * AveVelH *
03020 gv.bin[nd]->IntArea/4. * gv.bin[nd]->cnv_H_pCM3;
03021
03022 ASSERT( gv.bin[nd]->rate_h2_form_grains_HM79 > 0. );
03023 }
03024 else
03025 # endif
03026 {
03027
03028
03029
03030
03031
03032
03033
03034
03035 double conversion_efficiency_HM79 = 1/(1. + 1e4*sexp(920./gv.bin[nd]->tedust));
03036
03037
03038
03039
03040 gv.bin[nd]->rate_h2_form_grains_HM79 = 0.5 * AveVelH * gv.bin[nd]->IntArea/4. *
03041
03042 gv.bin[nd]->cnv_H_pCM3 * sticking_probability_H * conversion_efficiency_HM79;
03043 ASSERT( gv.bin[nd]->rate_h2_form_grains_HM79 > 0. );
03044 }
03045
03046 # ifndef IGNORE_QUANTUM_HEATING
03047 if( lgUseQHeat )
03048 {
03049
03050
03051
03052 double f = 1e-10;
03053
03054
03055 double sqrt_term = 35.399494936611667;
03056
03057 gv.bin[nd]->rate_h2_form_grains_CT02 = 0.;
03058
03059 for( k=0; k < qnbin; k++ )
03060 {
03061 double beta_alpha = 0.25 * sqrt_term *sexp(200./qtemp[k] );
03062
03063 double xi = 1./ (1. + 1.3e13*sexp(1.5*1e4/qtemp[k])*sqrt_term/(2.*f) );
03064
03065 double beta = 3e12 * sexp( 320. / qtemp[k] );
03066
03067
03068 double recombination_efficiency_CT02 = xi / (1. + 0.005*f/2./SDIV(beta) + beta_alpha );
03069 sticking_probability_H = 1./(1. + 0.04*sqrt(qtemp[k]+phycon.te) +
03070 0.002*phycon.te + 8e-6*phycon.te*phycon.te);
03071
03072
03073
03074
03075 gv.bin[nd]->rate_h2_form_grains_CT02 += qprob[k] * sticking_probability_H *
03076 recombination_efficiency_CT02;
03077 }
03078
03079
03080
03081
03082
03083 gv.bin[nd]->rate_h2_form_grains_CT02 *= 0.5 * AveVelH *
03084 gv.bin[nd]->IntArea/4. * gv.bin[nd]->cnv_H_pCM3;
03085
03086 ASSERT( gv.bin[nd]->rate_h2_form_grains_CT02 > 0. );
03087
03088 free(qtemp);
03089 free(qprob);
03090 }
03091 else
03092 # endif
03093 {
03094
03095
03096
03097 double f = 1e-10;
03098
03099
03100 double sqrt_term = 35.399494936611667;
03101 double beta_alpha = 0.25 * sqrt_term *sexp(200./gv.bin[nd]->tedust );
03102
03103 double xi = 1./ (1. + 1.3e13*sexp(1.5*1e4/ gv.bin[nd]->tedust )*sqrt_term/(2.*f) );
03104
03105 double beta = 3e12 * sexp( 320. / gv.bin[nd]->tedust );
03106
03107
03108 double recombination_efficiency_CT02 = xi / (1. + 0.005*f/2./SDIV(beta) + beta_alpha );
03109
03110
03111 fixit();
03112 #if 0
03113
03114 if( gv.bin[nd]->matType == MAT_CAR || gv.bin[nd]->matType == MAT_CAR2 )
03115 {
03116 double Td = gv.bin[nd]->tedust;
03117 recombination_efficiency_CT02 = 1./(
03118 ( 1. + 4.609569629185726e24*sexp(45000./Td) ) *
03119 ( 1. + sexp(800./Td) / (0.5389970511202651 * sexp(540./Td) + 5.6333909478365e-14*sqrt(Td)) )
03120 );
03121
03122 }
03123
03124 else if( gv.bin[nd]->matType == MAT_SIL || gv.bin[nd]->matType == MAT_SIL2 )
03125 {
03126 double Td = gv.bin[nd]->tedust;
03127 recombination_efficiency_CT02 = 1./(
03128 ( 1. + 6.998161265697586e24*sexp(45000./Td) ) *
03129 ( 1. + sexp(450./Td) / (0.4266153643741715 * sexp(340./Td) + 2.5335919594255e-14*sqrt(Td)) )
03130 );
03131
03132 }
03133 #endif
03134
03135
03136
03137
03138
03139 gv.bin[nd]->rate_h2_form_grains_CT02 = 0.5 * AveVelH * gv.bin[nd]->IntArea/4. *
03140 gv.bin[nd]->cnv_H_pCM3 * sticking_probability_H * recombination_efficiency_CT02;
03141 ASSERT( gv.bin[nd]->rate_h2_form_grains_CT02 > 0. );
03142 }
03143
03144 # ifndef IGNORE_QUANTUM_HEATING
03145
03146 sticking_probability_H = 1./(1. + 0.04*sqrt(gv.bin[nd]->tedust+phycon.te) +
03147 0.002*phycon.te + 8e-6*phycon.te*phycon.te);
03148 # endif
03149
03150
03151
03152
03153
03154
03155
03156
03157
03158
03159
03160 h2.rate_grain_op_conserve += AveVelH2 * gv.bin[nd]->IntArea/4. *
03161 gv.bin[nd]->cnv_H_pCM3 * sticking_probability_H;
03162
03163
03164
03165
03166
03167
03168
03169
03170
03179
03180
03181
03182
03183
03184
03185
03186 T_ortho_para_crit = 2. * hmi.Tad / log( POW2(60. *1.1e11)*hmi.Tad);
03187 if( gv.bin[nd]->tedust < T_ortho_para_crit )
03188 {
03189 double efficiency_opr = sexp(60.*1.1e11*sqrt(hmi.Tad)*sexp(hmi.Tad/gv.bin[nd]->tedust));
03190
03191
03192 h2.rate_grain_J1_to_J0 += AveVelH2 * gv.bin[nd]->IntArea/4. *
03193 gv.bin[nd]->cnv_H_pCM3 * sticking_probability_H * efficiency_opr;
03194 }
03195 }
03196
03197
03198
03199
03200
03201
03202
03203
03204
03205 h2.rate_grain_op_conserve *= h2.lgH2_grain_deexcitation;
03206 h2.rate_grain_J1_to_J0 *= h2.lgH2_grain_deexcitation;
03207
03208 }
03209 else
03210 {
03211
03212 for( size_t nd=0; nd < gv.bin.size(); nd++ )
03213 {
03214 gv.bin[nd]->rate_h2_form_grains_CT02 = 0.;
03215 gv.bin[nd]->rate_h2_form_grains_HM79 = 0.;
03216 }
03217
03218 h2.rate_grain_op_conserve = 0.;
03219
03220 h2.rate_grain_J1_to_J0 = 0.;
03221 }
03222
03223
03224
03225
03226
03227 gv.rate_h2_form_grains_used_total = 0.;
03228 for( size_t nd=0; nd < gv.bin.size(); nd++ )
03229 {
03230 if( hmi.chJura == 'C' )
03231 {
03232
03233
03234
03235 gv.bin[nd]->rate_h2_form_grains_used =
03236 gv.bin[nd]->rate_h2_form_grains_CT02*hmi.ScaleJura;
03237 gv.rate_h2_form_grains_used_total += gv.bin[nd]->rate_h2_form_grains_used;
03238 }
03239 else if( hmi.chJura == 'T' )
03240 {
03241
03242 gv.bin[nd]->rate_h2_form_grains_used =
03243 gv.bin[nd]->rate_h2_form_grains_HM79*hmi.ScaleJura;
03244 gv.rate_h2_form_grains_used_total += gv.bin[nd]->rate_h2_form_grains_used;
03245 }
03246 else if( hmi.chJura == 'S' )
03247 {
03248
03249
03250 gv.bin[nd]->rate_h2_form_grains_used =
03251 3e-18 * phycon.sqrte / gv.bin.size() * dense.gas_phase[ipHYDROGEN]*hmi.ScaleJura;
03252
03253 gv.rate_h2_form_grains_used_total += gv.bin[nd]->rate_h2_form_grains_used;
03254 }
03255
03256
03257
03258 else if( hmi.chJura == 'F' )
03259 {
03260
03261
03262 gv.bin[nd]->rate_h2_form_grains_used = hmi.rate_h2_form_grains_set*dense.gas_phase[ipHYDROGEN] / gv.bin.size();
03263 gv.rate_h2_form_grains_used_total += gv.bin[nd]->rate_h2_form_grains_used;
03264 }
03265 }
03266 ASSERT( gv.rate_h2_form_grains_used_total >= 0. );
03267
03268 # ifndef IGNORE_QUANTUM_HEATING
03269 printf( " fnzone %.2f H2 rate %.4e\n", fnzone, gv.rate_h2_form_grains_used_total );
03270 # endif
03271
03272
03273
03274
03275 }
03276
03277 STATIC void mole_h_reactions( void )
03278 {
03279 static double teused=-1;
03280 double exph2,
03281 exph2s,
03282 exphp,
03283 ex3hp;
03284 long i;
03285 double h2esc,
03286 th2,
03287 cr_H2s ,
03288 cr_H2dis,
03289 cr_H2dis_ELWERT_H2g,
03290 cr_H2dis_ELWERT_H2s;
03291
03292 DEBUG_ENTRY( "hmole_reactions()" );
03293
03294
03295
03296
03297 bool need_update = ! fp_equal( phycon.te, teused );
03298
03299 if( need_update )
03300 {
03301 teused = phycon.te;
03302
03303
03304
03305
03306 hmi.exphmi = sexp(8.745e3/phycon.te);
03307 if( hmi.exphmi > 0. )
03308 {
03309
03310 hmi.rel_pop_LTE_Hmin = SAHA/(phycon.te32*hmi.exphmi)*(1./(2.*2.));
03311 }
03312 else
03313 {
03314 hmi.rel_pop_LTE_Hmin = 0.;
03315 }
03316
03317
03318
03319 exphp = sexp(3.072e4/phycon.te);
03320 if( exphp > 0. )
03321 {
03322
03323
03324 hmi.rel_pop_LTE_H2p = SAHA/(phycon.te32*exphp)*(4./(2.*1.))*3.634e-5;
03325 }
03326 else
03327 {
03328 hmi.rel_pop_LTE_H2p = 0.;
03329 }
03330
03331
03332
03333 ex3hp = sexp(1.882e4/phycon.te);
03334 if( ex3hp > 0. )
03335 {
03336
03337
03338 hmi.rel_pop_LTE_H3p = SAHA/(phycon.te32*ex3hp)*(4./(2.*1.))*3.634e-5;
03339 }
03340 else
03341 {
03342 hmi.rel_pop_LTE_H3p = 0.;
03343 }
03344 }
03345
03346
03347
03348
03349
03350
03351 if( h2.lgEnabled && h2.lgEvaluated && hmi.lgH2_Chemistry_BigH2 )
03352 {
03353
03354 hmi.rel_pop_LTE_H2g = h2.rel_pop_LTE_g;
03355 hmi.rel_pop_LTE_H2s = h2.rel_pop_LTE_s;
03356 }
03357 else
03358 {
03359 if (need_update)
03360 {
03361
03362 exph2 = sexp((5.195e4)/phycon.te);
03363
03364
03365 if( exph2 > 0. )
03366 {
03367
03368 hmi.rel_pop_LTE_H2g = SAHA/(phycon.te32*exph2)*(1./(2.*2.))*3.634e-5;
03369 }
03370 else
03371 {
03372 hmi.rel_pop_LTE_H2g = 0.;
03373 }
03374
03375
03376
03377
03378
03379 exph2s = sexp(2.178e4/phycon.te);
03380
03381 if( exph2s > 0. )
03382 {
03383
03384 hmi.rel_pop_LTE_H2s = SAHA/(phycon.te32*exph2s)*(1./(2.*2.))*3.634e-5;
03385 }
03386 else
03387 {
03388 hmi.rel_pop_LTE_H2s = 0.;
03389 }
03390 }
03391 }
03392 {
03393
03394
03395
03396
03397 enum {DEBUG_LOC=false};
03398
03399 if( DEBUG_LOC && nzone>187&& iteration > 1)
03400 {
03401
03402 fprintf(ioQQQ,"ph2lteee\t%.2e\t%.1e\t%.1e\n",
03403 hmi.rel_pop_LTE_H2g,
03404 sexp(2.178e4/phycon.te),
03405 phycon.te);
03406 }
03407 }
03408
03409
03410 hmi.hmicol = hmirat(phycon.te)*EN1RYD*phycon.te*1.15e-5;
03411
03412 fixit();
03413 for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
03414 {
03415 if( (*diatom)->lgEnabled && mole_global.lgStancil )
03416 (*diatom)->Mol_Photo_Diss_Rates();
03417 }
03418
03419
03420
03421
03422
03423 hmi.hmicol *= dense.eden*findspecieslocal("H")->den;
03424
03425
03426
03427
03428
03429
03432
03433
03434
03435 t_phoHeat photoHeat;
03436
03437
03438
03439 hmi.HMinus_photo_rate = GammaBn( hmi.iphmin-1 , iso_sp[ipHE_LIKE][ipHELIUM].fb[0].ipIsoLevNIonCon , opac.iphmop ,
03440 0.055502 , &hmi.HMinus_induc_rec_rate , &hmi.HMinus_induc_rec_cooling, &photoHeat );
03441
03442
03443 hmi.HMinus_photo_heat = photoHeat.HeatNet;
03444
03445 {
03446
03447
03448 enum {DEBUG_LOC=false};
03449
03450 if( DEBUG_LOC)
03451 {
03452 fprintf(ioQQQ,"hminphoto\t%li\t%li\t%.2e\n", nzone, conv.nPres2Ioniz , hmi.HMinus_photo_rate );
03453 }
03454 }
03455
03456
03457 hmi.HMinus_induc_rec_rate *= hmi.rel_pop_LTE_Hmin*dense.eden;
03458
03459
03461 hmi.HMinus_induc_rec_cooling *= hmi.rel_pop_LTE_Hmin*dense.eden*findspecieslocal("H")->den;
03462
03463 {
03464
03465
03466 enum {DEBUG_LOC=false};
03467
03468 if( DEBUG_LOC && nzone>400)
03469 {
03470 fprintf(ioQQQ,"hmoledebugg %.2f ",fnzone);
03471 GammaPrt(
03472 hmi.iphmin-1 , iso_sp[ipHE_LIKE][ipHELIUM].fb[0].ipIsoLevNIonCon , opac.iphmop ,
03473
03474 ioQQQ,
03475
03476 hmi.HMinus_photo_rate,
03477
03478 hmi.HMinus_photo_rate*0.05);
03479 }
03480 }
03481
03482
03483
03484
03485
03486
03487
03488
03489
03490
03491
03492 hmi.UV_Cont_rel2_Habing_TH85_depth = 0.;
03493
03494
03495
03496
03497
03498 for( i=rfield.ipG0_TH85_lo; i < rfield.ipG0_TH85_hi; ++i )
03499 {
03500 hmi.UV_Cont_rel2_Habing_TH85_depth += ((rfield.flux[0][i-1]) + (rfield.ConInterOut[i-1])+
03501 (rfield.outlin[0][i-1])+ (rfield.outlin_noplot[i-1]))*rfield.anu[i-1];
03502 }
03503
03504
03505
03506
03507
03508 hmi.UV_Cont_rel2_Habing_TH85_depth = (realnum)(hmi.UV_Cont_rel2_Habing_TH85_depth*EN1RYD/1.6e-3);
03509
03510 if( nzone<=1 )
03511 {
03512 hmi.UV_Cont_rel2_Habing_TH85_face = hmi.UV_Cont_rel2_Habing_TH85_depth;
03513 }
03514
03515
03516
03517 hmi.UV_Cont_rel2_Habing_spec_depth = 0.;
03518 for( i=rfield.ipG0_spec_lo; i < rfield.ipG0_spec_hi; ++i )
03519 {
03520 hmi.UV_Cont_rel2_Habing_spec_depth += ( rfield.flux[0][i-1] + rfield.ConInterOut[i-1]+
03521 rfield.outlin[0][i-1]+ rfield.outlin_noplot[i-1])*rfield.anu[i-1];
03522 }
03523 hmi.UV_Cont_rel2_Habing_spec_depth = (realnum)(hmi.UV_Cont_rel2_Habing_spec_depth*EN1RYD/1.6e-3);
03524
03525
03526
03527 if( hmi.UV_Cont_rel2_Draine_DB96_face ==0 )
03528 {
03529
03530 for( i=rfield.ipG0_DB96_lo; i < rfield.ipG0_DB96_hi; ++i )
03531 {
03532 hmi.UV_Cont_rel2_Draine_DB96_face += ( rfield.flux[0][i-1] + rfield.ConInterOut[i-1]+
03533 rfield.outlin[0][i-1]+ rfield.outlin_noplot[i-1]);
03534 }
03535
03536
03537
03538
03539
03540
03541
03542
03543 hmi.UV_Cont_rel2_Draine_DB96_face = hmi.UV_Cont_rel2_Draine_DB96_face/(1.232e7f * 1.71f);
03544 }
03545
03546
03547
03548 hmi.H2Opacity = (realnum)(1.2e-14*(1e5/GetDopplerWidth(2.f*dense.AtomicWeight[ipHYDROGEN])));
03549
03550 th2 = (colden.colden[ipCOL_H2g]+ colden.colden[ipCOL_H2s])*hmi.H2Opacity;
03551
03552
03553 h2esc = esc_PRD_1side(th2,1e-4);
03554
03555
03556
03557
03558
03559
03560
03561 hmi.H2_Solomon_dissoc_rate_TH85_H2g = 3.4e-11 * hmi.UV_Cont_rel2_Habing_TH85_depth * h2esc;
03562 hmi.H2_Solomon_dissoc_rate_TH85_H2s = 3.4e-11 * hmi.UV_Cont_rel2_Habing_TH85_depth * h2esc;
03563 hmi.H2_H2g_to_H2s_rate_TH85 = hmi.H2_Solomon_dissoc_rate_TH85_H2g*9.;
03564
03565
03566 hmi.H2_Solomon_dissoc_rate_BHT90_H2g = 3.4e-11 * hmi.UV_Cont_rel2_Habing_TH85_depth * h2esc;
03567 hmi.H2_Solomon_dissoc_rate_BHT90_H2s = 3.4e-11 * hmi.UV_Cont_rel2_Habing_TH85_depth * h2esc;
03568 hmi.H2_H2g_to_H2s_rate_BHT90 = hmi.H2_Solomon_dissoc_rate_BHT90_H2g*9.;
03569
03570 {
03571
03572
03573
03574
03575 double x = (colden.colden[ipCOL_H2g]+colden.colden[ipCOL_H2s]) / 5e14;
03576 double sqrtx = sqrt(1.+x);
03577
03578 double b5 = GetDopplerWidth(2.f*dense.AtomicWeight[ipHYDROGEN])/1e5;
03579
03580 double fshield = 0.965/POW2(1.+x/b5) + 0.035/sqrtx *
03581 sexp(8.5e-4*sqrtx);
03582
03583
03584
03585
03586
03587
03588 hmi.UV_Cont_rel2_Draine_DB96_depth = hmi.UV_Cont_rel2_Draine_DB96_face *
03589 (realnum)(sexp( opac.TauAbsFace[rfield.ip1000A-1] )/radius.r1r0sq);
03590
03591
03592
03593
03594 if( !mole_global.lgLeidenHack )
03595 {
03596
03597 hmi.H2_Solomon_dissoc_rate_BD96_H2g = 4.6e-11 * hmi.UV_Cont_rel2_Draine_DB96_depth * fshield;
03598 hmi.H2_Solomon_dissoc_rate_BD96_H2s = 4.6e-11 * hmi.UV_Cont_rel2_Draine_DB96_depth * fshield;
03599 }
03600 else
03601 {
03602
03603 hmi.H2_Solomon_dissoc_rate_BD96_H2g = 5.18e-11* (hmi.UV_Cont_rel2_Habing_TH85_face/1.66f)
03604 *sexp(3.02*rfield.extin_mag_V_point)* fshield;
03605 hmi.H2_Solomon_dissoc_rate_BD96_H2s = 5.18e-11* (hmi.UV_Cont_rel2_Habing_TH85_face/1.66f)
03606 *sexp(3.02*rfield.extin_mag_V_point)* fshield;
03607 }
03608
03609
03610
03611
03612
03613 hmi.H2_H2g_to_H2s_rate_BD96 = 5.67* hmi.H2_Solomon_dissoc_rate_BD96_H2g;
03614 }
03615
03616
03617 if( hmi.UV_Cont_rel2_Draine_DB96_face > SMALLFLOAT )
03618 {
03619
03620
03621
03622
03623
03624 double b5 = GetDopplerWidth(2.f*dense.AtomicWeight[ipHYDROGEN])/1e5;
03625
03626
03627
03628
03629
03630
03631
03632 double x_H2g, x_H2s,
03633 fshield_H2g, fshield_H2s,
03634 f_H2s;
03635 static double a_H2g, a_H2s,
03636 e1_H2g, e1_H2s,
03637 e2_H2g,
03638 b_H2g,
03639 sl_H2g, sl_H2s,
03640 k_f_H2s,
03641 k_H2g_to_H2s,
03642 log_G0_face = -1;
03643
03644
03645
03646
03647
03648
03649
03650
03651 {
03652 if(hmi.UV_Cont_rel2_Draine_DB96_face <= 1.)
03653 {
03654 log_G0_face = 0.;
03655 }
03656 else if(hmi.UV_Cont_rel2_Draine_DB96_face >= 1e7)
03657 {
03658 log_G0_face = 7.;
03659 }
03660 else
03661 {
03662 log_G0_face = log10(hmi.UV_Cont_rel2_Draine_DB96_face);
03663 }
03664
03665
03666
03667
03668 a_H2g = 0.06 * log_G0_face + 1.32;
03669 a_H2s = 0.375 * log_G0_face + 2.125;
03670
03671 e1_H2g = -0.05 * log_G0_face + 2.25;
03672 e1_H2s = -0.125 * log_G0_face + 2.625;
03673
03674 e2_H2g = -0.005 * log_G0_face + 0.625;
03675
03676 b_H2g = -4.0e-3 * log_G0_face + 3.2e-2;
03677
03678
03679 sl_H2g = 4.0e14;
03680 sl_H2s = 9.0e15;
03681
03682
03683 k_f_H2s = MAX2(0.1,2.375 * log_G0_face - 1.875 );
03684
03685
03686 k_H2g_to_H2s = MAX2(1.,-1.75 * log_G0_face + 11.25);
03687
03688
03689
03690
03691 }
03692
03693
03694 f_H2s = k_f_H2s * pow((double)hmi.UV_Cont_rel2_Draine_DB96_depth, 0.2 );
03695
03696
03697 x_H2g = (colden.colden[ipCOL_H2g]) / sl_H2g;
03698 x_H2s = (colden.colden[ipCOL_H2s]) / sl_H2s;
03699
03700
03701 fshield_H2g = 0.965/pow(1.+x_H2g/b5,e1_H2g) + b_H2g/pow(1.+x_H2g/b5,e2_H2g);
03702 fshield_H2s = 0.965/pow(1.+x_H2s/b5,e1_H2s);
03703
03704
03705 hmi.H2_Solomon_dissoc_rate_ELWERT_H2g = a_H2g * 4.6e-11 * fshield_H2g * hmi.UV_Cont_rel2_Draine_DB96_depth;
03706 hmi.H2_Solomon_dissoc_rate_ELWERT_H2s = a_H2s * 4.6e-11 * fshield_H2s * (hmi.UV_Cont_rel2_Draine_DB96_depth + f_H2s);
03707
03708
03709 hmi.H2_H2g_to_H2s_rate_ELWERT = k_H2g_to_H2s * hmi.H2_Solomon_dissoc_rate_ELWERT_H2g;
03710
03711
03712
03713
03714 hmi.H2_photodissoc_ELWERT_H2s = hmi.UV_Cont_rel2_Draine_DB96_depth*1e-11;
03715 hmi.H2_photodissoc_ELWERT_H2g = hmi.H2_photodissoc_ELWERT_H2s * 1.0e-10;
03716 }
03717 else
03718 {
03719 hmi.H2_Solomon_dissoc_rate_ELWERT_H2g = 0.;
03720 hmi.H2_Solomon_dissoc_rate_ELWERT_H2s = 0.;
03721 hmi.H2_photodissoc_ELWERT_H2s = 0.;
03722 hmi.H2_photodissoc_ELWERT_H2g = 0.;
03723 }
03724
03725
03726 hmi.H2_photodissoc_TH85 = hmi.UV_Cont_rel2_Habing_TH85_depth*1e-11;
03727
03728
03729 hmi.H2_photodissoc_BHT90 = hmi.UV_Cont_rel2_Habing_TH85_depth*1e-11;
03730
03731
03732
03733
03734
03735
03736 cr_H2s = secondaries.x12tot*0.9 / 2. * hmi.lgLeidenCRHack;
03737
03738
03739 cr_H2dis = secondaries.x12tot*0.1 / 2. * hmi.lgLeidenCRHack;
03740
03741
03742
03743
03744
03745 cr_H2dis_ELWERT_H2g = secondaries.x12tot*5e-8 * hmi.lgLeidenCRHack;
03746 cr_H2dis_ELWERT_H2s = secondaries.x12tot*4e-2 * hmi.lgLeidenCRHack;
03747
03748
03749
03750
03751
03752
03753
03754
03755
03756
03757
03758 if( h2.lgEnabled && h2.lgEvaluated && hmi.lgH2_Chemistry_BigH2 )
03759 {
03760
03761
03762
03763
03764
03765
03766 hmi.H2_Solomon_dissoc_rate_used_H2g = h2.Solomon_dissoc_rate_g;
03767 ASSERT( hmi.H2_Solomon_dissoc_rate_used_H2g >= 0. );
03768
03769 hmi.H2_Solomon_dissoc_rate_used_H2s = h2.Solomon_dissoc_rate_s;
03770 ASSERT( hmi.H2_Solomon_dissoc_rate_used_H2s >= 0. );
03771
03772
03773 hmi.H2_H2g_to_H2s_rate_used = h2.gs_rate();
03774 ASSERT( hmi.H2_H2g_to_H2s_rate_used >= 0. );
03775
03776
03777
03778
03779
03780 hmi.H2_photodissoc_used_H2s = h2.photodissoc_BigH2_H2s;
03781
03782 if( mole_global.lgStancil && h2.lgEnabled )
03783 hmi.H2_photodissoc_used_H2s = h2.Cont_Dissoc_Rate_H2s;
03784 ASSERT( hmi.H2_photodissoc_used_H2s >= 0. );
03785
03786
03787
03788 hmi.H2_photodissoc_used_H2g = h2.photodissoc_BigH2_H2g;
03789
03790 if( mole_global.lgStancil && h2.lgEnabled )
03791 hmi.H2_photodissoc_used_H2g = h2.Cont_Dissoc_Rate_H2g;
03792 ASSERT( hmi.H2_photodissoc_used_H2g >= 0. );
03793 }
03794 else if( hmi.chH2_small_model_type == 'T' )
03795 {
03796
03797
03798 hmi.H2_Solomon_dissoc_rate_used_H2g = hmi.H2_Solomon_dissoc_rate_TH85_H2g + cr_H2dis;
03799
03800 hmi.H2_Solomon_dissoc_rate_used_H2s = hmi.H2_Solomon_dissoc_rate_TH85_H2s + cr_H2dis;
03801 hmi.H2_H2g_to_H2s_rate_used = hmi.H2_H2g_to_H2s_rate_TH85 + cr_H2s;
03802
03803
03804
03805 hmi.H2_photodissoc_used_H2s = hmi.H2_photodissoc_TH85;
03806
03807
03808 hmi.H2_photodissoc_used_H2g = hmi.H2_photodissoc_TH85*1.0e-10f;
03809 }
03810
03811 else if( hmi.chH2_small_model_type == 'H' )
03812 {
03813
03814
03815
03816 hmi.H2_Solomon_dissoc_rate_used_H2g = hmi.H2_Solomon_dissoc_rate_BHT90_H2g + cr_H2dis;
03817
03818 hmi.H2_Solomon_dissoc_rate_used_H2s = hmi.H2_Solomon_dissoc_rate_BHT90_H2s + cr_H2dis;
03819 hmi.H2_H2g_to_H2s_rate_used = hmi.H2_H2g_to_H2s_rate_BHT90 + cr_H2s;
03820
03821
03822
03823 hmi.H2_photodissoc_used_H2s = hmi.H2_photodissoc_BHT90;
03824
03825
03826 hmi.H2_photodissoc_used_H2g = hmi.H2_photodissoc_BHT90*1.0e-10f;
03827 }
03828
03829 else if( hmi.chH2_small_model_type == 'B' )
03830 {
03831
03832
03833 hmi.H2_Solomon_dissoc_rate_used_H2g = hmi.H2_Solomon_dissoc_rate_BD96_H2g + cr_H2dis;
03834
03835 hmi.H2_Solomon_dissoc_rate_used_H2s = hmi.H2_Solomon_dissoc_rate_BD96_H2s + cr_H2dis;
03836
03837 hmi.H2_H2g_to_H2s_rate_used = hmi.H2_H2g_to_H2s_rate_BD96 + cr_H2s;
03838
03839
03840
03841
03842 hmi.H2_photodissoc_used_H2s = hmi.H2_photodissoc_TH85;
03843
03844
03845 hmi.H2_photodissoc_used_H2g = hmi.H2_photodissoc_TH85*1.0e-10f;
03846 }
03847 else if( hmi.chH2_small_model_type == 'E' )
03848 {
03849
03850
03851 hmi.H2_Solomon_dissoc_rate_used_H2g = hmi.H2_Solomon_dissoc_rate_ELWERT_H2g + cr_H2dis_ELWERT_H2g;
03852 hmi.H2_Solomon_dissoc_rate_used_H2s = hmi.H2_Solomon_dissoc_rate_ELWERT_H2s + cr_H2dis_ELWERT_H2s;
03853 hmi.H2_H2g_to_H2s_rate_used = hmi.H2_H2g_to_H2s_rate_ELWERT + cr_H2s;
03854
03855
03856
03857
03858 hmi.H2_photodissoc_used_H2s = hmi.H2_photodissoc_ELWERT_H2s;
03859 hmi.H2_photodissoc_used_H2g = hmi.H2_photodissoc_ELWERT_H2g;
03860 }
03861 else
03862 TotalInsanity();
03863
03864 {
03865
03866 enum {DEBUG_LOC=false};
03867
03868 if( DEBUG_LOC && h2.lgEnabled )
03869 {
03870 fprintf(ioQQQ," Solomon H2 dest rates: TH85 %.2e BD96 %.2e Big %.2e excit rates: TH85 %.2e Big %.2e\n",
03871 hmi.H2_Solomon_dissoc_rate_TH85_H2g,
03872 hmi.H2_Solomon_dissoc_rate_BD96_H2g,
03873 h2.Solomon_dissoc_rate_g,
03874 hmi.H2_H2g_to_H2s_rate_TH85 ,
03875 h2.gs_rate() );
03876 }
03877 }
03878 return;
03879 }
03880 mole_reaction *mole_findrate_s(const char buf[])
03881 {
03882 DEBUG_ENTRY("mole_findrate_s()");
03883
03884 string newbuf = canonicalize_reaction_label(buf);
03885
03886 mole_reaction_i p = mole_priv::reactab.find(newbuf);
03887
03888 if(p != mole_priv::reactab.end())
03889 return &(*p->second);
03890 else
03891 return NULL;
03892 }
03893
03894 double t_mole_local::findrk(const char buf[]) const
03895 {
03896 DEBUG_ENTRY("t_mole_local::findrk()");
03897
03898 mole_reaction *rate = mole_findrate_s(buf);
03899
03900 if(!rate)
03901 return 0.;
03902
03903
03904 ASSERT( !isnan( reaction_rks[ rate->index ] ) );
03905
03906 return reaction_rks[ rate->index ];
03907 }
03908 double t_mole_local::findrate(const char buf[]) const
03909 {
03910 double ret;
03911 int i;
03912
03913 DEBUG_ENTRY("t_mole_local::findrate()");
03914
03915 mole_reaction *rate = mole_findrate_s(buf);
03916 if(!rate)
03917 {
03918 return 0.;
03919 }
03920
03921 ret = reaction_rks[ rate->index ];
03922 for(i=0;i<rate->nreactants;i++)
03923 ret *= species[ rate->reactants[i]->index ].den;
03924
03925 return ret;
03926 }
03927
03928
03929
03930 double t_mole_local::sink_rate_tot(const char chSpecies[]) const
03931 {
03932 DEBUG_ENTRY("t_mole_local::sink_rate_tot()");
03933
03934 const molecule* const sp = findspecies(chSpecies);
03935 double ratev = sink_rate_tot(sp);
03936
03937 return ratev;
03938 }
03939 double t_mole_local::sink_rate_tot(const molecule* const sp) const
03940 {
03941 DEBUG_ENTRY("t_mole_local::sink_rate_tot()");
03942 double ratev = 0;
03943
03944 for(mole_reaction_i p=mole_priv::reactab.begin();
03945 p != mole_priv::reactab.end(); ++p)
03946 {
03947 mole_reaction &rate = *p->second;
03948 ratev += sink_rate( sp, rate );
03949 }
03950
03951 return ratev;
03952 }
03953
03954 double t_mole_local::sink_rate(const molecule* const sp, const char buf[]) const
03955 {
03956 const mole_reaction* const rate = mole_findrate_s(buf);
03957 return sink_rate( sp, *rate );
03958 }
03959
03960 double t_mole_local::sink_rate(const molecule* const sp, const mole_reaction& rate) const
03961 {
03962 DEBUG_ENTRY("t_mole_local::sink_rate()");
03963
03964 int ipthis = -1;
03965 for(int i=0;i<rate.nreactants && ipthis == -1;i++)
03966 {
03967 if(rate.reactants[i] == sp && rate.rvector[i]==NULL && rate.rvector_excit[i]==NULL )
03968 {
03969 ipthis = i;
03970 }
03971 }
03972 if(ipthis != -1)
03973 {
03974 double ratevi = rate.a * rate.rk();
03975 for(int i=0;i<rate.nreactants;i++)
03976 {
03977 if(i!=ipthis)
03978 {
03979 ratevi *= species[ rate.reactants[i]->index ].den;
03980 }
03981 }
03982 return ratevi;
03983 }
03984 else
03985 return 0.;
03986 }
03987
03988 #ifdef PRINTREACT
03989 STATIC void printreact(mole_reaction *rate)
03990 {
03991 int i;
03992
03993 DEBUG_ENTRY("printreact()");
03994
03995 for(i=0;i<rate->nreactants;i++) {
03996 fprintf(stderr,"%s,",rate->reactants[i]->label);
03997 }
03998 fprintf(stderr,"=>");
03999 for(i=0;i<rate->nproducts;i++) {
04000 fprintf(stderr,"%s,",rate->products[i]->label);
04001 }
04002 fprintf(stderr,"\n");
04003
04004 }
04005 #endif
04006
04007 double t_mole_local::dissoc_rate(const char chSpecies[]) const
04008 {
04009 DEBUG_ENTRY("t_mole_local::dissoc_rate()");
04010
04011 molecule *sp = findspecies(chSpecies);
04012 if (sp == null_mole)
04013 return 0.0;
04014 ASSERT(sp->isMonatomic());
04015 const chem_atom *tgt = sp->nAtom.begin()->first.get_ptr();
04016 molecule *ph = findspecies("PHOTON");
04017 double ratev = 0.0;
04018
04019 for (mole_reaction_i p
04020 =mole_priv::reactab.begin(); p != mole_priv::reactab.end(); ++p)
04021 {
04022 mole_reaction &rate = *p->second;
04023
04024
04025 int ipph = 0;
04026 for (int i=0;i<rate.nreactants;i++)
04027 {
04028 if (rate.reactants[i] == ph)
04029 ipph++;
04030 }
04031 if (!ipph)
04032 continue;
04033
04034
04035
04036 int ipspin = 0, ipfreein = 0;
04037 for (int i=0;i<rate.nreactants;i++)
04038 {
04039 if (rate.reactants[i] == sp)
04040 ++ipspin;
04041 if (rate.reactants[i]->isMonatomic() && tgt == sp->nAtom.begin()->first.get_ptr())
04042 ++ipfreein;
04043 }
04044 int ipspout = 0, ipfreeout = 0;
04045 for (int i=0;i<rate.nproducts;i++)
04046 {
04047 if (rate.products[i] == sp)
04048 ++ipspout;
04049 if (rate.products[i]->isMonatomic() && tgt == sp->nAtom.begin()->first.get_ptr())
04050 ++ipfreeout;
04051 }
04052
04053
04054 int newsp = ipspout-ipspin;
04055 if (newsp <= 0)
04056 continue;
04057
04058
04059 int nbondsbroken = ipfreeout-ipfreein;
04060 if (nbondsbroken <= 0)
04061 continue;
04062
04063 double fracbroken = nbondsbroken/((double)ipfreeout);
04064 ASSERT( fracbroken <= 1.0 );
04065
04066 double ratevi = reaction_rks[ rate.index ];
04067 for (int i=0;i<rate.nreactants;i++)
04068 {
04069 ratevi *= species[ rate.reactants[i]->index ].den;
04070 }
04071
04072
04073
04074
04075 double ratesp = ratevi*newsp;
04076
04077 ratesp *= fracbroken;
04078
04079
04080 ratev += ratesp;
04081 }
04082 return ratev;
04083 }
04084 double t_mole_local::source_rate_tot(const char chSpecies[]) const
04085 {
04086 DEBUG_ENTRY("t_mole_local::source_rate_tot()");
04087
04088 molecule *sp = findspecies(chSpecies);
04089 double ratev = source_rate_tot(sp);
04090
04091 return ratev;
04092 }
04093 double t_mole_local::source_rate_tot(const molecule* const sp) const
04094 {
04095 DEBUG_ENTRY("t_mole_local::source_rate_tot()");
04096 double ratev = 0;
04097
04098 for (mole_reaction_i p =mole_priv::reactab.begin(); p != mole_priv::reactab.end(); ++p)
04099 {
04100 mole_reaction &rate = *p->second;
04101 int ipthis = 0;
04102 for(int i=0;i<rate.nproducts;i++)
04103 {
04104 if( rate.products[i] == sp && rate.pvector[i]==NULL && rate.pvector_excit[i]==NULL )
04105 {
04106 ipthis++;
04107 }
04108 }
04109 if(ipthis)
04110 {
04111 double ratevi = rate.a * rate.rk();
04112 for(int i=0;i<rate.nreactants;i++)
04113 {
04114 ratevi *= species[ rate.reactants[i]->index ].den;
04115 }
04116 ratev += ipthis*ratevi;
04117 }
04118 }
04119
04120 return ratev;
04121 }
04122
04123 double t_mole_local::chem_heat(void) const
04124 {
04125
04126
04127
04128
04129
04130
04131
04132
04133
04134
04135
04136 DEBUG_ENTRY("t_mole_local::chem_heat()");
04137
04138 double heating = 0.;
04139 map<double,string> heatMap;
04140 molecule *ph = findspecies("PHOTON");
04141 molecule *crph = findspecies("CRPHOT");
04142 molecule *grn = findspecies("grn");
04143
04144
04145 for (mole_reaction_i p
04146 =mole_priv::reactab.begin(); p != mole_priv::reactab.end(); ++p)
04147 {
04148 mole_reaction &rate = *p->second;
04149
04150
04151 bool lgCanSkip = false;
04152 for (int i=0;i<rate.nproducts;i++)
04153 {
04154 if( rate.products[i] == ph || rate.products[i] == crph )
04155 lgCanSkip = true;
04156 }
04157 for (int i=0;i<rate.nreactants;i++)
04158 {
04159 if( rate.reactants[i] == ph || rate.reactants[i] == crph )
04160 lgCanSkip = true;
04161 }
04162
04163 for (int i=0;i<rate.nreactants;i++)
04164 {
04165 if( rate.reactants[i] == grn && rate.rvector[i] != NULL )
04166 lgCanSkip = true;
04167 }
04168
04169 if( lgCanSkip )
04170 continue;
04171
04172
04173 double rate_tot = reaction_rks[ rate.index ];
04174 for( long i=0; i < rate.nreactants; ++i )
04175 {
04176 rate_tot *= species[ rate.reactants[i]->index ].den;
04177 }
04178
04179 realnum reaction_enthalpy = 0.;
04180
04181
04182 for( long i=0; i < rate.nreactants; ++i )
04183 {
04184 reaction_enthalpy += rate.reactants[i]->form_enthalpy;
04185 }
04186
04187
04188 for( long i=0; i < rate.nproducts; ++i )
04189 {
04190 reaction_enthalpy -= rate.products[i]->form_enthalpy;
04191 }
04192
04193
04194
04195
04196
04197
04198 double heat = reaction_enthalpy*rate_tot*(1e10/AVOGADRO);
04199 heatMap[heat] = rate.label;
04200 heating += heat;
04201 }
04202
04203
04204 long index = 0;
04205
04206
04207
04208
04209 for( map<double,string>::reverse_iterator it = heatMap.rbegin(); it != heatMap.rend(); ++it, ++index )
04210 {
04211 fprintf( ioQQQ, "DEBUGGG heat %li\t%li\t%.6e\t%s\n", index, nzone, it->first, it->second.c_str() );
04212 if( index==2 )
04213 break;
04214 }
04215 index = 0;
04216 for( map<double,string>::iterator it = heatMap.begin(); it != heatMap.end(); ++it, ++index )
04217 {
04218 if( it->first >= 0. )
04219 break;
04220 fprintf( ioQQQ, "DEBUGGG cool %li\t%li\t%.6e\t%s\n", index, nzone, it->first, it->second.c_str() );
04221 if( index==2 )
04222 break;
04223 }
04224
04225 return heating;
04226 }
04227
04228
04229 void mole_punch(FILE *punit, const char speciesname[], const char args[], bool lgHeader, bool lgData, double depth)
04230 {
04231 int i, ipthis;
04232 double ratevi;
04233 molecule *sp;
04234
04235 DEBUG_ENTRY("mole_punch()");
04236
04237 sp = findspecies(speciesname);
04238
04239 if (lgHeader)
04240 {
04241 fprintf (punit,"#Depth");
04242 }
04243 if (lgData)
04244 {
04245 fprintf (punit,"%.5e",depth);
04246 }
04247
04248 for (mole_reaction_i p
04249 =mole_priv::reactab.begin(); p != mole_priv::reactab.end(); ++p)
04250 {
04251 mole_reaction &rate = *p->second;
04252 ipthis = 0;
04253
04254 for (i=0;i<rate.nreactants;i++)
04255 {
04256 if( rate.reactants[i] == sp )
04257 {
04258 if( ( strcmp( args, "DEST" )==0 && rate.rvector[i]==NULL && rate.rvector_excit[i]==NULL ) ||
04259 ( strcmp( args, "DFLT" )==0 && rate.rvector[i]==NULL && rate.rvector_excit[i]==NULL ) ||
04260 ( strcmp( args, "CATA" )==0 && rate.rvector[i]!=NULL ) ||
04261 strcmp( args, "ALL " )==0 )
04262 ipthis++;
04263 }
04264 }
04265
04266 for(i=0;i<rate.nproducts;i++)
04267 {
04268 if( rate.products[i] == sp )
04269 {
04270 if( ( strcmp( args, "CREA" )==0 && rate.pvector[i]==NULL && rate.pvector_excit[i]==NULL ) ||
04271 ( strcmp( args, "DFLT" )==0 && rate.pvector[i]==NULL && rate.pvector_excit[i]==NULL ) ||
04272 ( strcmp( args, "CATA" )==0 && rate.pvector[i]!=NULL ) ||
04273 strcmp( args, "ALL " )==0 )
04274 ipthis++;
04275 }
04276 }
04277
04278 if(ipthis)
04279 {
04280 if (lgHeader)
04281 {
04282 fprintf(punit,"\t%s",rate.label.c_str());
04283 }
04284 if (lgData)
04285 {
04286 ratevi = mole.reaction_rks[ rate.index ];
04287 for(i=0;i<rate.nreactants;i++)
04288 {
04289 ratevi *= mole.species[ rate.reactants[i]->index ].den;
04290 }
04291 fprintf(punit,"\t%.3e",ratevi);
04292 }
04293 }
04294 }
04295 fprintf(punit,"\n");
04296 }
04297