51 STATIC void newreact(
const char label[],
const char fun[],
double a,
double b,
double c);
83 template<
class T>
void newfunc()
89 ratefunc findfunc(
const char name[]);
92 class mole_reaction_hmrate;
94 class mole_reaction_th85rate;
96 class mole_reaction_crnurate_noalbedo;
100 class mole_reaction_grn_react;
102 class mole_reaction_h2_collh_excit;
104 class mole_reaction_h2_collh2_excit;
106 class mole_reaction_h2_collh_deexc;
108 class mole_reaction_h2_collh2_deexc;
110 class mole_reaction_rh2g_dis_h;
112 class mole_reaction_rh2s_dis_h;
114 class mole_reaction_rh2g_dis_h2;
116 class mole_reaction_rh2s_dis_h2;
118 class mole_reaction_rh2s_dis_h2_nodeexcit;
120 class mole_reaction_bh2g_dis_h;
122 class mole_reaction_bh2s_dis_h;
124 class mole_reaction_bh2g_dis_h2;
126 class mole_reaction_hneut;
128 class mole_reaction_cionhm;
145 typedef mole_reaction_null T;
147 virtual T*
Create()
const {
return new T;}
148 virtual const char*
name() {
return "null";}
186 for(n=0;n<nreact;n++)
213 return pow(te/300.,rate->
b)*exp(-rate->
c/te);
218 typedef mole_reaction_hmrate_exo T;
220 virtual T*
Create()
const {
return new T;}
222 virtual const char*
name() {
return "hmrate_exo";}
233 te =
MIN2( te, -10. * this->c );
235 return pow(te/300.,this->b)*exp(-te/this->c);
242 typedef mole_reaction_hmrate T;
244 virtual T*
Create()
const {
return new T;}
245 virtual const char*
name() {
return "hmrate";}
254 typedef mole_reaction_constrate T;
256 virtual T*
Create()
const {
return new T;}
257 virtual const char*
name() {
return "constrate";}
298 typedef mole_reaction_th85rate T;
300 virtual T*
Create()
const {
return new T;}
301 virtual const char*
name() {
return "th85rate";}
304 return th85rate(
this);
330 typedef mole_reaction_crnurate_noalbedo T;
332 virtual T*
Create()
const {
return new T;}
334 virtual const char*
name() {
return "crnurate_noalbedo";}
337 return crnurate_noalbedo(
this);
343 typedef mole_reaction_crnurate T;
345 virtual T*
Create()
const {
return new T;}
346 virtual const char*
name() {
return "crnurate";}
349 return crnurate_noalbedo(
this)/(1.0-
albedo);
358 typedef mole_reaction_cryden_ov_bg T;
360 virtual T*
Create()
const {
return new T;}
361 virtual const char*
name() {
return "cryden_ov_bg";}
370 typedef mole_reaction_co_lnu_c_o_lnu T;
372 virtual T*
Create()
const {
return new T;}
373 virtual const char*
name() {
return "co_lnu_c_o_lnu";}
385 for( ns=0; ns<2; ++ns )
431 typedef mole_reaction_vib_evap T;
433 virtual T*
Create()
const {
return new T;}
434 virtual const char*
name() {
return "vib_evap";}
437 double binding_energy, exponent, vib_freq;
443 binding_energy = this->
b;
444 double bin_total=0.0;
445 for(
size_t nd=0; nd <
gv.
bin.size() ; nd++ )
447 double bin_area =
gv.
bin[nd]->IntArea*
gv.
bin[nd]->cnv_H_pCM3;
448 exponent += exp(-binding_energy/
gv.
bin[nd]->tedust)*bin_area;
449 bin_total += bin_area;
451 exponent /= bin_total;
452 const double surface_density_of_sites = 1.5e15;
455 vib_freq = sqrt(2*surface_density_of_sites*(0.3*
BOLTZMANN)*binding_energy/(
PI*
PI*this->reactants[0]->mole_mass));
475 typedef mole_reaction_grn_abs T;
477 virtual T*
Create()
const {
return new T;}
478 virtual const char*
name() {
return "grn_abs";}
486 if (this->reactants[0]->n_nuclei() != 0)
487 mass = this->reactants[0]->mole_mass;
489 mass = this->reactants[1]->mole_mass;
497 typedef mole_reaction_grn_react T;
499 virtual T*
Create()
const {
return new T;}
500 virtual const char*
name() {
return "grn_react";}
503 return grn_react(
this);
527 double quant_barrier = 1e-8;
528 double surface_density_of_sites = 1.5e15;
533 double activ_barrier = rate->
c;
542 double dust_density = 0.0;
547 for(
size_t nd=0; nd <
gv.
bin.size(); nd++ )
549 double bin_density =
gv.
bin[nd]->IntArea*
gv.
bin[nd]->cnv_H_pCM3;
550 Exp_i += exp(-E_i/
gv.
bin[nd]->tedust)*bin_density;
551 Exp_j += exp(-E_j/
gv.
bin[nd]->tedust)*bin_density;
552 dust_density += bin_density/(4*1e-10);
560 double diff_i = vib_freq_i*Exp_i/total_sites;
561 double diff_j = vib_freq_j*Exp_j/total_sites;
568 return scale*(diff_i + diff_j)/
SDIV(dust_density);
573 typedef mole_reaction_grn_photo T;
575 virtual T*
Create()
const {
return new T;}
576 virtual const char*
name() {
return "grn_photo";}
604 typedef mole_reaction_th85rate_co T;
606 virtual T*
Create()
const {
return new T;}
608 virtual const char*
name() {
return "th85rate_co";}
612 double esc_co, column;
632 if (this->reactants[0]->n_nuclei() != 0)
633 column =
mole.
species[ this->reactants[0]->index ].column;
635 column =
mole.
species[ this->reactants[1]->index ].column;
637 esc_co = 4.4e-15 * column /
642 return esca0k2(esc_co)*th85rate(
this);
648 typedef mole_reaction_oh_c2h2_co_ch3 T;
650 virtual T*
Create()
const {
return new T;}
651 virtual const char*
name() {
return "oh_c2h2_co_ch3";}
670 typedef mole_reaction_h_hnc_hcn_h T;
672 virtual T*
Create()
const {
return new T;}
674 virtual const char*
name() {
return "h_hnc_hcn_h";}
710 return 1. - 4.938e-6;
715 double frac_H2star_grains(
void)
741 typedef mole_reaction_h2ph3p T;
743 virtual T*
Create()
const {
return new T;}
744 virtual const char*
name() {
return "h2ph3p";}
766 typedef mole_reaction_hopexch T;
768 virtual T*
Create()
const {
return new T;}
769 virtual const char*
name() {
return "hopexch";}
778 typedef mole_reaction_hpoexch T;
780 virtual T*
Create()
const {
return new T;}
781 virtual const char*
name() {
return "hpoexch";}
790 typedef mole_reaction_hmattach T;
792 virtual T*
Create()
const {
return new T;}
793 virtual const char*
name() {
return "hmattach";}
802 typedef mole_reaction_hmihph2p T;
804 virtual T*
Create()
const {
return new T;}
805 virtual const char*
name() {
return "hmihph2p";}
829 typedef mole_reaction_hmphoto T;
831 virtual T*
Create()
const {
return new T;}
832 virtual const char*
name() {
return "hmphoto";}
862 typedef mole_reaction_cionhm T;
864 virtual T*
Create()
const {
return new T;}
865 virtual const char*
name() {
return "cionhm";}
872 double assoc_detach(
void)
887 y=545969508.1323510+x*71239.23653059864;
893 typedef mole_reaction_c3bod T;
895 virtual T*
Create()
const {
return new T;}
896 virtual const char*
name() {
return "c3bod";}
905 typedef mole_reaction_asdfg T;
907 virtual T*
Create()
const {
return new T;}
908 virtual const char*
name() {
return "asdfg";}
917 typedef mole_reaction_asdfs T;
919 virtual T*
Create()
const {
return new T;}
920 virtual const char*
name() {
return "asdfs";}
929 typedef mole_reaction_asdbg T;
931 virtual T*
Create()
const {
return new T;}
932 virtual const char*
name() {
return "asdbg";}
942 typedef mole_reaction_asdbs T;
944 virtual T*
Create()
const {
return new T;}
945 virtual const char*
name() {
return "asdbs";}
955 typedef mole_reaction_bhneut T;
957 virtual T*
Create()
const {
return new T;}
958 virtual const char*
name() {
return "bhneut";}
969 return (hneut(
this)*ratio*
992 return 3.4738192887404660e-008;
997 typedef mole_reaction_hneut T;
999 virtual T*
Create()
const {
return new T;}
1000 virtual const char*
name() {
return "hneut";}
1010 typedef mole_reaction_h2_spon_diss T;
1012 virtual T*
Create()
const {
return new T;}
1013 virtual const char*
name() {
return "h2_spon_diss";}
1034 virtual const char*
name() {
return "grnh2tot";}
1037 return grnh2tot(
this);
1043 typedef mole_reaction_grnh2 T;
1045 virtual T*
Create()
const {
return new T;}
1046 virtual const char*
name() {
return "grnh2";}
1049 return grnh2tot(
this)*(1.-frac_H2star_grains());
1055 typedef mole_reaction_grnh2s T;
1057 virtual T*
Create()
const {
return new T;}
1058 virtual const char*
name() {
return "grnh2s";}
1061 return grnh2tot(
this)*frac_H2star_grains();
1067 typedef mole_reaction_radasc T;
1069 virtual T*
Create()
const {
return new T;}
1070 virtual const char*
name() {
return "radasc";}
1078 return hmrate(
this) *
1090 typedef mole_reaction_assoc_ion T;
1092 virtual T*
Create()
const {
return new T;}
1093 virtual const char*
name() {
return "assoc_ion";}
1100 return hmrate(
this) *
1131 typedef mole_reaction_rh2g_dis_h T;
1133 virtual T*
Create()
const {
return new T;}
1134 virtual const char*
name() {
return "rh2g_dis_h";}
1137 return rh2g_dis_h(
this);
1150 return HMRATE(4.67e-7,-1.,5.5e4);
1155 typedef mole_reaction_rh2s_dis_h T;
1157 virtual T*
Create()
const {
return new T;}
1158 virtual const char*
name() {
return "rh2s_dis_h";}
1161 return rh2s_dis_h(
this);
1175 return HMRATE(5.5e-29*0.5/(
SAHA*3.634e-5)*sqrt(300.),0.5,5.195e4);
1180 typedef mole_reaction_rh2g_dis_h2 T;
1182 virtual T*
Create()
const {
return new T;}
1183 virtual const char*
name() {
return "rh2g_dis_h2";}
1186 return rh2g_dis_h2(
this);
1199 return HMRATE(1e-11,0.,0.);
1204 typedef mole_reaction_rh2s_dis_h2 T;
1206 virtual T*
Create()
const {
return new T;}
1207 virtual const char*
name() {
return "rh2s_dis_h2";}
1210 return rh2s_dis_h2(
this);
1223 return HMRATE(1e-11,0.,2.18e4);
1226 class mole_reaction_rh2s_dis_h2_nodeexcit :
public mole_reaction
1228 typedef mole_reaction_rh2s_dis_h2_nodeexcit T;
1230 virtual T*
Create()
const {
return new T;}
1231 virtual const char*
name() {
return "rh2s_dis_h2_nodeexcit";}
1234 return rh2s_dis_h2_nodeexcit(
this);
1244 typedef mole_reaction_bh2g_dis_h T;
1246 virtual T*
Create()
const {
return new T;}
1247 virtual const char*
name() {
return "bh2g_dis_h";}
1250 return bh2g_dis_h(
this);
1260 typedef mole_reaction_bh2s_dis_h T;
1262 virtual T*
Create()
const {
return new T;}
1263 virtual const char*
name() {
return "bh2s_dis_h";}
1266 return bh2s_dis_h(
this);
1276 typedef mole_reaction_bh2g_dis_h2 T;
1278 virtual T*
Create()
const {
return new T;}
1279 virtual const char*
name() {
return "bh2g_dis_h2";}
1282 return bh2g_dis_h2(
this);
1288 typedef mole_reaction_bh2s_dis_h2 T;
1290 virtual T*
Create()
const {
return new T;}
1291 virtual const char*
name() {
return "bh2s_dis_h2";}
1300 typedef mole_reaction_h2photon T;
1302 virtual T*
Create()
const {
return new T;}
1303 virtual const char*
name() {
return "h2photon";}
1312 typedef mole_reaction_h2crphot T;
1314 virtual T*
Create()
const {
return new T;}
1315 virtual const char*
name() {
return "h2crphot";}
1324 typedef mole_reaction_h2crphh T;
1326 virtual T*
Create()
const {
return new T;}
1327 virtual const char*
name() {
return "h2crphh";}
1343 typedef mole_reaction_h2scrphh T;
1345 virtual T*
Create()
const {
return new T;}
1346 virtual const char*
name() {
return "h2scrphh";}
1362 typedef mole_reaction_radath T;
1364 virtual T*
Create()
const {
return new T;}
1365 virtual const char*
name() {
return "radath";}
1374 typedef mole_reaction_gamtwo T;
1376 virtual T*
Create()
const {
return new T;}
1377 virtual const char*
name() {
return "gamtwo";}
1387 typedef mole_reaction_hlike_phot T;
1389 virtual T*
Create()
const {
return new T;}
1390 virtual const char*
name() {
return "hlike_phot";}
1406 typedef mole_reaction_h2s_sp_decay T;
1408 virtual T*
Create()
const {
return new T;}
1409 virtual const char*
name() {
return "h2s_sp_decay";}
1438 typedef mole_reaction_h2_collh2_deexc T;
1440 virtual T*
Create()
const {
return new T;}
1441 virtual const char*
name() {
return "h2_collh2_deexc";}
1444 return h2_collh2_deexc(
this);
1461 typedef mole_reaction_h2_collh_deexc T;
1463 virtual T*
Create()
const {
return new T;}
1464 virtual const char*
name() {
return "h2_collh_deexc";}
1467 return h2_collh_deexc(
this);
1485 typedef mole_reaction_h2_collh2_excit T;
1487 virtual T*
Create()
const {
return new T;}
1488 virtual const char*
name() {
return "h2_collh2_excit";}
1491 return h2_collh2_excit(
this);
1509 typedef mole_reaction_h2_collh_excit T;
1511 virtual T*
Create()
const {
return new T;}
1512 virtual const char*
name() {
return "h2_collh_excit";}
1515 return h2_collh_excit(
this);
1521 typedef mole_reaction_h2gexcit T;
1523 virtual T*
Create()
const {
return new T;}
1524 virtual const char*
name() {
return "h2gexcit";}
1533 typedef mole_reaction_h2sdissoc T;
1535 virtual T*
Create()
const {
return new T;}
1536 virtual const char*
name() {
return "h2sdissoc";};
1551 typedef mole_reaction_h2gdissoc T;
1553 virtual T*
Create()
const {
return new T;}
1554 virtual const char*
name() {
return "h2gdissoc";}
1569 typedef mole_reaction_hd_photodissoc T;
1571 virtual T*
Create()
const {
return new T;}
1572 virtual const char*
name() {
return "hd_photodissoc";}
1598 else if( te < 1200. )
1602 else if( te < 3800. )
1607 else if( te <= 1.4e4 )
1630 typedef mole_reaction_gamheh T;
1632 virtual T*
Create()
const {
return new T;}
1633 virtual const char*
name() {
return "gamheh";}
1645 for( i=
hmi.
iheh1-1; i < limit; i++ )
1658 enum {exclude, base, umisthack, federman, lithium, deuterium, ti, misc};
1662 static bool lgReactInitialized =
false;
1672 if( lgReactInitialized )
1675 lgReactInitialized =
true;
1684 newfunc<mole_reaction_null>();
1685 newfunc<mole_reaction_hmrate>();
1686 newfunc<mole_reaction_hmrate_exo>();
1687 newfunc<mole_reaction_constrate>();
1688 newfunc<mole_reaction_th85rate>();
1689 newfunc<mole_reaction_crnurate>();
1690 newfunc<mole_reaction_crnurate_noalbedo>();
1691 newfunc<mole_reaction_co_lnu_c_o_lnu>();
1692 newfunc<mole_reaction_vib_evap>();
1693 newfunc<mole_reaction_th85rate_co>();
1694 newfunc<mole_reaction_grn_abs>();
1696 newfunc<mole_reaction_grn_react>();
1698 newfunc<mole_reaction_grn_photo>();
1699 newfunc<mole_reaction_oh_c2h2_co_ch3>();
1700 newfunc<mole_reaction_h_hnc_hcn_h>();
1702 newfunc<mole_reaction_gamheh>();
1703 newfunc<mole_reaction_hd_photodissoc>();
1704 newfunc<mole_reaction_h2gdissoc>();
1705 newfunc<mole_reaction_h2sdissoc>();
1706 newfunc<mole_reaction_h2gexcit>();
1707 newfunc<mole_reaction_h2_collh_excit>();
1708 newfunc<mole_reaction_h2_collh2_excit>();
1709 newfunc<mole_reaction_h2_collh_deexc>();
1710 newfunc<mole_reaction_h2_collh2_deexc>();
1711 newfunc<mole_reaction_h2s_sp_decay>();
1712 newfunc<mole_reaction_hlike_phot>();
1713 newfunc<mole_reaction_gamtwo>();
1714 newfunc<mole_reaction_h2ph3p>();
1715 newfunc<mole_reaction_radath>();
1716 newfunc<mole_reaction_cryden_ov_bg>();
1717 newfunc<mole_reaction_h2scrphh>();
1718 newfunc<mole_reaction_h2crphh>();
1719 newfunc<mole_reaction_h2photon>();
1720 newfunc<mole_reaction_h2crphot>();
1721 newfunc<mole_reaction_rh2g_dis_h>();
1722 newfunc<mole_reaction_rh2s_dis_h>();
1723 newfunc<mole_reaction_rh2g_dis_h2>();
1724 newfunc<mole_reaction_rh2s_dis_h2>();
1725 newfunc<mole_reaction_rh2s_dis_h2_nodeexcit>();
1726 newfunc<mole_reaction_bh2g_dis_h>();
1727 newfunc<mole_reaction_bh2s_dis_h>();
1728 newfunc<mole_reaction_bh2g_dis_h2>();
1729 newfunc<mole_reaction_bh2s_dis_h2>();
1730 newfunc<mole_reaction_radasc>();
1731 newfunc<mole_reaction_assoc_ion>();
1732 newfunc<mole_reaction_grnh2>();
1733 newfunc<mole_reaction_grnh2s>();
1734 newfunc<mole_reaction_h2_spon_diss>();
1735 newfunc<mole_reaction_bhneut>();
1736 newfunc<mole_reaction_hneut>();
1737 newfunc<mole_reaction_asdbs>();
1738 newfunc<mole_reaction_asdbg>();
1739 newfunc<mole_reaction_asdfs>();
1740 newfunc<mole_reaction_asdfg>();
1741 newfunc<mole_reaction_c3bod>();
1742 newfunc<mole_reaction_cionhm>();
1743 newfunc<mole_reaction_hmphoto>();
1744 newfunc<mole_reaction_hmihph2p>();
1745 newfunc<mole_reaction_hmattach>();
1780 newreact(
"C+,OH=>CO,H+",
"hmrate",0.,0.,0.);
1783 newreact(
"H,H,grn=>H2,grn",
"grnh2",1.,0.,0.);
1784 newreact(
"H,H,grn=>H2*,grn",
"grnh2s",1.,0.,0.);
1785 newreact(
"H-,PHOTON=>H,e-",
"hmphoto",1.,0.,0.);
1794 long nelem = (*atom)->el->Z-1;
1798 sprintf(react,
"H-,%s+=>H,%s", (*atom)->label().c_str(), (*atom)->label().c_str() );
1799 newreact(react,
"hmrate",4e-6/sqrt(300.),-0.5,0.);
1803 newreact(
"H,e-=>H-,PHOTON",
"hmattach",1.,0.,0.);
1804 newreact(
"H-,H+=>H2+,e-",
"hmihph2p",1.,0.,0.);
1805 newreact(
"H-,e-=>H,e-,e-",
"cionhm",1.,0.,0.);
1806 newreact(
"H,e-,e-=>H-,e-",
"c3bod",1.,0.,0.);
1807 newreact(
"H,H-=>H2,e-",
"asdfg",1.,0.,0.);
1808 newreact(
"H,H-=>H2*,e-",
"asdfs",1.,0.,0.);
1809 newreact(
"H2,e-=>H,H-",
"asdbg",1.,0.,0.);
1810 newreact(
"H2*,e-=>H,H-",
"asdbs",1.,0.,0.);
1811 newreact(
"H-,H+=>H,H",
"hneut",1.,0.,0.);
1812 newreact(
"H,H=>H-,H+",
"bhneut",1.,0.,0.);
1813 newreact(
"H2*=>H,H",
"h2_spon_diss",1.,0.,0.);
1816 newreact(
"H,H=>H2,PHOTON",
"radasc",3e-14,0.,0.);
1817 newreact(
"H,H=>H2+,e-",
"assoc_ion",3.27e-10,-0.35,17829.);
1818 newreact(
"H2,H=>H,H,H",
"rh2g_dis_h",1.,0.,0.);
1819 newreact(
"H2,H2=>H,H,H2",
"rh2g_dis_h2",1.,0.,0.);
1820 newreact(
"H,H,H=>H2,H",
"bh2g_dis_h",1.,0.,0.);
1821 newreact(
"H,H,H2=>H2,H2",
"bh2g_dis_h2",1.,0.,0.);
1823 newreact(
"H2,PHOTON=>H2+,e-",
"h2photon",1.,0.,0.);
1824 newreact(
"H2,CRPHOT=>H2+,e-",
"h2crphot",1.,0.,0.);
1825 newreact(
"H2*,PHOTON=>H2+,e-",
"h2photon",1.,0.,0.);
1826 newreact(
"H2*,CRPHOT=>H2+,e-",
"h2crphot",1.,0.,0.);
1827 newreact(
"H2,CRPHOT=>H,H",
"h2crphh",1.,0.,0.);
1828 newreact(
"H2,CRPHOT=>H+,H-",
"cryden_ov_bg",3.9e-21,0.,0.);
1829 newreact(
"H2*,CRPHOT=>H+,H-",
"cryden_ov_bg",3.9e-21,0.,0.);
1830 newreact(
"H2,CRPHOT=>H+,H,e-",
"cryden_ov_bg",2.2e-19,0.,0.);
1831 newreact(
"H2*,CRPHOT=>H+,H,e-",
"cryden_ov_bg",2.2e-19,0.,0.);
1832 newreact(
"H+,H=>H2+,PHOTON",
"radath",1.,0.,0.);
1833 newreact(
"H2+,PHOTON=>H,H+",
"gamtwo",1.,0.,0.);
1834 newreact(
"H2+,CRPHOT=>H,H+",
"hlike_phot",1.,0.,0.);
1837 newreact(
"H3+,CRPHOT=>H2+,H+,e-",
"hlike_phot",2.,0.,0.);
1838 newreact(
"H,H+,e-=>H2+,e-",
"hmrate",2e-7 *
SAHA * (4./(2.*1.)) * 3.634e-5 * pow(300.,-1.50),-1.50,0.);
1842 newreact(
"H2,CRPHOT=>H2+,e-",
"hmrate",4.4e-17,0.,0.);
1843 newreact(
"H2,CRPHOT=>H,H+,e-",
"hmrate",1e-19,0.,0.);
1844 newreact(
"H2*,CRPHOT=>H,H+,e-",
"hmrate",1e-19,0.,0.);
1845 newreact(
"H2*,CRPHOT=>H2+,e-",
"hmrate",4.4e-17,0.,0.);
1846 newreact(
"H2,CRPHOT=>H,H",
"hmrate",5e-18,0.,0.);
1847 newreact(
"e-,H3+=>H2,H",
"hmrate",2.5e-8,-0.3,0.);
1848 newreact(
"e-,H3+=>H,H,H",
"hmrate",7.5e-8,-0.3,0.);
1849 newreact(
"H+,H=>H2+,PHOTON",
"hmrate",5.3e-19,1.85,0);
1853 newreact(
"H2+,e-=>H,H",
"hmrate",1.6e-8,-0.43,0.);
1854 newreact(
"H2+,PHOTON=>H,H+",
"th85rate",5.7e-10,0.,1.9);
1857 newreact(
"H2*,CRPHOT=>H,H",
"h2scrphh",1.,0.,0.);
1858 newreact(
"H2,H2+=>H,H3+",
"h2ph3p",1.,0.,0.);
1859 newreact(
"H2*=>H2,PHOTON",
"h2s_sp_decay",1.,0.,0.);
1860 newreact(
"H2*,H2=>H2,H2",
"h2_collh2_deexc",1.,0.,0.);
1861 newreact(
"H2*,H=>H2,H",
"h2_collh_deexc",1.,0.,0.);
1862 newreact(
"H2,H2=>H2*,H2",
"h2_collh2_excit",1.,0.,0.);
1863 newreact(
"H2,H=>H2*,H",
"h2_collh_excit",1.,0.,0.);
1866 newreact(
"H2*,H=>H,H,H",
"rh2s_dis_h",1.,0.,0.);
1867 newreact(
"H2*,H2=>H2,H,H",
"rh2s_dis_h2",1.,0.,0.);
1868 newreact(
"H2*,H2*=>H2,H,H",
"rh2s_dis_h2",1.,0.,0.);
1869 newreact(
"H2*,H2*=>H2*,H,H",
"rh2s_dis_h2_nodeexcit",1.,0.,0.);
1873 newreact(
"H2,He=>He,H,H",
"rh2g_dis_h",1.,0.,0.);
1874 newreact(
"H2,H+=>H+,H,H",
"rh2g_dis_h",1.,0.,0.);
1875 newreact(
"H2,H3+=>H3+,H,H",
"rh2g_dis_h",1.,0.,0.);
1876 newreact(
"H2,e-=>e-,H,H",
"rh2g_dis_h",1.,0.,0.);
1877 newreact(
"H2*,He=>He,H,H",
"rh2s_dis_h",1.,0.,0.);
1878 newreact(
"H2*,H+=>H+,H,H",
"rh2s_dis_h",1.,0.,0.);
1879 newreact(
"H2*,H3+=>H3+,H,H",
"rh2s_dis_h",1.,0.,0.);
1880 newreact(
"H2*,e-=>e-,H,H",
"rh2s_dis_h",1.,0.,0.);
1883 newreact(
"He,H,H=>H2,He",
"bh2g_dis_h",1.,0.,0.);
1884 newreact(
"H+,H,H=>H2,H+",
"bh2g_dis_h",1.,0.,0.);
1885 newreact(
"H3+,H,H=>H2,H3+",
"bh2g_dis_h",1.,0.,0.);
1886 newreact(
"e-,H,H=>H2,e-",
"bh2g_dis_h",1.,0.,0.);
1887 newreact(
"H,H,H=>H2*,H",
"bh2s_dis_h",1.,0.,0.);
1888 newreact(
"H2,H,H=>H2*,H2",
"bh2s_dis_h2",1.,0.,0.);
1889 newreact(
"H2,H,H=>H2*,H2*",
"bh2s_dis_h2",1.,0.,0.);
1890 newreact(
"H2*,H,H=>H2*,H2*",
"bh2s_dis_h2",1.,0.,0.);
1891 newreact(
"He,H,H=>H2*,He",
"bh2s_dis_h",1.,0.,0.);
1892 newreact(
"H+,H,H=>H2*,H+",
"bh2s_dis_h",1.,0.,0.);
1893 newreact(
"H3+,H,H=>H2*,H3+",
"bh2s_dis_h",1.,0.,0.);
1894 newreact(
"e-,H,H=>H2*,e-",
"bh2s_dis_h",1.,0.,0.);
1897 newreact(
"H2,PHOTON=>H2*",
"h2gexcit",1.,0.,0.);
1898 newreact(
"H2*,PHOTON=>H,H",
"h2sdissoc",1.,0.,0.);
1899 newreact(
"H2,PHOTON=>H,H",
"h2gdissoc",1.,0.,0.);
1900 newreact(
"HeH+,PHOTON=>H,He+",
"gamheh",1.,0.,0.);
1905 newreact(
"OHgrn,Hgrn=>H2Ogrn",
"grn_react",1.,0.,0.);
1910 fprintf(stderr,
"Finished testing against UDFA database\n");
1924 it->second->index = index;
1937 DEBUG_ENTRY(
"mole_generate_isotopologue_reactions()");
1939 bool lgDebug =
false;
1942 vector<string> newReactionStrings;
1943 vector<mole_reaction*> oldReactions;
1948 bool lgFound =
false;
1950 for(
long i=0; i<it->second->nreactants; ++i )
1953 for(
nAtoms_i atom=it->second->reactants[i]->nAtom.begin(); atom != it->second->reactants[i]->nAtom.end(); ++atom )
1955 if( atom->first->label() == atom_old )
1969 fprintf(
ioQQQ,
"DEBUGGG mole_generate_isotopologue_reactions %s ..........\n", it->first.c_str() );
1971 for(
long i=0; i<it->second->nreactants; ++i )
1977 vector<string> react_iso_labels;
1979 vector< int > numAtoms;
1980 string embellishments;
1982 bool lgParseOK =
parse_species_label( it->second->reactants[i]->label.c_str(), atomsLeftToRight, numAtoms, embellishments );
1983 ASSERT( lgParseOK ==
true );
1992 for(
unsigned j=0; j<react_iso_labels.size(); ++j )
1994 for(
long k=0; k<it->second->nproducts; ++k )
2000 vector<string> prod_iso_labels;
2001 atomsLeftToRight.clear();
2003 embellishments.clear();
2005 lgParseOK =
parse_species_label( it->second->products[k]->label.c_str(), atomsLeftToRight, numAtoms, embellishments );
2006 ASSERT( lgParseOK ==
true );
2015 for(
unsigned l=0; l<prod_iso_labels.size(); ++l )
2017 string react_string;
2020 for(
long i1=0; i1<i; ++i1 )
2022 react_string += it->second->reactants[i1]->label;
2023 if( i1 != it->second->nreactants-1 )
2024 react_string +=
",";
2026 react_string += react_iso_labels[j];
2027 if( i != it->second->nreactants-1 )
2028 react_string +=
",";
2029 for(
long i2=i+1; i2<it->second->nreactants; ++i2 )
2031 react_string += it->second->reactants[i2]->label;
2032 if( i2 != it->second->nreactants-1 )
2033 react_string +=
",";
2036 react_string +=
"=>";
2038 for(
long k1=0; k1<k; ++k1 )
2040 react_string += it->second->products[k1]->label;
2041 if( k1 != it->second->nproducts-1 )
2042 react_string +=
",";
2044 react_string += prod_iso_labels[l];
2045 if( k != it->second->nproducts-1 )
2046 react_string +=
",";
2047 for(
long k2=k+1; k2<it->second->nproducts; ++k2 )
2049 react_string += it->second->products[k2]->label;
2050 if( k2 != it->second->nproducts-1 )
2051 react_string +=
",";
2055 fprintf(
ioQQQ,
"DEBUGGG mole_generate_isotopologue_reactions .................... %s\n",
2056 react_string.c_str() );
2060 newReactionStrings.push_back( canon_react_string );
2061 oldReactions.push_back( it->second.get_ptr() );
2068 ASSERT( oldReactions.size() == newReactionStrings.size() );
2071 vector<mole_reaction*>::const_iterator it2 = oldReactions.begin();
2072 for( vector<string>::const_iterator it1 = newReactionStrings.begin(); it1 != newReactionStrings.end(); ++it1, ++it2 )
2078 newreact( it1->c_str(), (*it2)->name(), (*it2)->a, (*it2)->b, (*it2)->c );
2088 char chLabel[50], chLabelSave[50];
2095 strcpy( chLabel, p->second->label.c_str() );
2096 strcpy( chLabelSave, p->second->label.c_str() );
2097 char *str = chLabel;
2098 const char *delim =
"=>";
2099 char *chNewProducts = strtok( str, delim );
2100 char *chNewReactants = strtok( NULL, delim );
2101 char chNewReaction[50];
2103 strcpy( chNewReaction, chNewReactants );
2104 strcat( chNewReaction,
"=>" );
2105 strcat( chNewReaction, chNewProducts );
2115 fprintf(
ioQQQ,
"Warning! No reverse reaction for %30s. ", p->second->label.c_str() );
2116 fprintf(
ioQQQ,
"\n" );
2137 double ln_result = 0.;
2143 ln_result += log(fac);
2151 ln_result -= log(fac);
2153 double result = exp( ln_result );
2163 if( sp->
label ==
"PHOTON" || sp->
label ==
"CRPHOT" )
2169 else if( sp->
label ==
"CRP" || sp->
label ==
"grn" )
2187 vector<long> numForm, numDest;
2205 for(
unsigned i=0; i<numForm.size(); ++i )
2207 if( numForm[i]==0 && numDest[i]==0 )
2209 if( numForm[i]>1 && numDest[i]>1 )
2222 char *label, *reactstr, *f;
2237 STATIC void newreact(
const char label[],
const char fun[],
double a,
double b,
double c)
2241 ratefunc rate = findfunc(fun);
2242 if (rate.get_ptr() == NULL)
2244 fprintf(stderr,
"Rate function %s not known for reaction %s. Aborting. Sorry.\n",fun,label);
2248 rate->label = label;
2262 const char *rateLabelPtr = rate->label.c_str();
2266 rate->udfastate =
ABSENT;
2270 if (rate->udfastate ==
ABSENT)
2272 fprintf(stderr,
"Reaction %s not in UDFA\n",rateLabelPtr);
2277 if (rate->nreactants == 2 && rate->reactants[0]->mole_mass!=0.0 && rate->reactants[1]->mole_mass!=0.0 )
2279 rate->reduced_mass = 1./(1./rate->reactants[0]->mole_mass+1./rate->reactants[1]->mole_mass);
2283 rate->reduced_mass = 0.;
2293 fprintf(
ioQQQ,
"Warning: duplicate reaction %s -- using new version\n",rateLabelPtr);
2314 bool lgProd =
false;
2316 for(
int i=0;!i || label[i-1]!=
'\0';i++)
2318 if(label[i] ==
',' || label[i] ==
'=' || label[i] ==
'\0')
2324 fprintf(
ioQQQ,
"Mole_reactions: ignoring reaction %s (species %s not active)\n",label,buf.c_str());
2332 fprintf(stderr,
"Mole_reactions: Too many reactants in %s, only %d allowed\n",label,
MAXREACTANTS);
2342 fprintf(stderr,
"Mole_reactions: Too many products in %s, only %d allowed\n",label,
MAXPRODUCTS);
2351 if (label[i] !=
'>')
2353 fprintf(
ioQQQ,
"Format error in reaction %s\n",label);
2376 ratefunc rate = findfunc(
"null");
2377 rate->label = label;
2414 rate->
label = newLabel;
2546 sparsefp = fopen(
"sparse.pbm",
"w");
2547 fprintf(sparsefp,
"P4\n%d %d\n",
2555 ch = (ch << 1) | (c[i][j] != 0.0);
2577 int dcharge, n,
sign;
2584 nel[it->first] += it->second;
2590 nel[it->first] -= it->second;
2595 fprintf(stderr,
"Reaction %s charge out of balance by %d\n",
2596 rate->
label.c_str(),dcharge);
2600 for(
nAtoms_i it = nel.begin(); it != nel.end(); ++it )
2608 fprintf(stderr,
"Error: reaction %s %s %d of element %s\n",
2609 rate->
label.c_str(),
sign==1?
"destroys":
"creates",
2611 it->first->label().c_str() );
2619 #define gzopen(file,mode) fopen(file,mode)
2620 #define gzgets(fp,buf,siz) fgets(buf,siz,fp)
2621 #define gzclose(fp) fclose(fp)
2629 class formula_species {
2634 bool operator< (
const formula_species &a,
const formula_species &b)
2639 if (a.reactants[i]<b.reactants[i])
2641 else if (a.reactants[i]>b.reactants[i])
2646 if (a.products[i]<b.products[i])
2648 else if (a.products[i]>b.products[i])
2654 struct udfa_reaction_s {
2658 double a, b, c, trange[2];
2662 static map <formula_species,struct udfa_reaction_s *>
udfatab;
2672 fprintf(stderr,
"Error, could not read %s\n",file);
2685 #define FLTEQ(a,b) (fabs((a)-(b)) <= 1e-6*fabs((a)+(b)))
2689 struct udfa_reaction_s *r;
2690 map <formula_species, struct udfa_reaction_s *>::iterator p;
2691 unsigned int havespecies = 1, i, n;
2696 r = (
struct udfa_reaction_s *)
MALLOC(
sizeof(
struct udfa_reaction_s));
2703 r->l.reactants[n] = NULL;
2715 if (!strncmp(f,
"CRPHOT",6))
2724 r->l.products[n] = NULL;
2762 r->source = f[0]?f[0]:
'?';
2765 r->trange[n] = atof(f);
2777 fprintf(stderr,
"Duplicate reaction\n");
2798 map<formula_species, struct udfa_reaction_s *>::iterator p;
2799 struct udfa_reaction_s *u;
2834 ratefunc findfunc (
const char name[])
2842 enum { DEBUG_MOLE =
false };
2854 long index = rate.
index;
2860 if (fabs(newrk-oldrk) > 0.1*newrk)
2861 fprintf(
ioQQQ,
"%s: %15.8g => %15.8g\n",
2862 rate.
label.c_str(),oldrk,newrk);
2868 enum { DEBUG_MOLE =
false };
2881 double bigchange = 0.;
2882 unsigned long bigindex = ULONG_MAX;
2887 sum = oldrk+newrk, diff = newrk-oldrk;
2890 double change = fabs(diff)/sum;
2891 if (change > bigchange)
2903 if (bigindex == (
unsigned) rate.
index)
2909 pct = 100.*(newrk-oldrk)/oldrk;
2910 fprintf(
ioQQQ,
"Zone %ld, big chemistry rate change %s:"
2911 " %15.8g => %15.8g (%6.2g%%)\n",
2927 double T_ortho_para_crit;
2937 # ifndef IGNORE_QUANTUM_HEATING
2950 for(
size_t nd=0; nd <
gv.
bin.size(); nd++ )
2952 # ifndef IGNORE_QUANTUM_HEATING
2954 double *qtemp, *qprob;
2964 double sticking_probability_H = 1./(1. + 0.04*sqrt(
gv.
bin[nd]->tedust+
phycon.
te) +
2967 # ifndef IGNORE_QUANTUM_HEATING
2971 qtemp = (
double*)
MALLOC((
size_t)(
NQGRID*
sizeof(double)));
2972 qprob = (
double*)
MALLOC((
size_t)(
NQGRID*
sizeof(double)));
2974 qheat(qtemp,qprob,&qnbin,nd);
2976 if(
gv.
bin[nd]->lgUseQHeat )
2984 qtemp[0] =
gv.
bin[nd]->tedust;
2987 gv.
bin[nd]->rate_h2_form_grains_HM79 = 0.;
2989 for( k=0; k < qnbin; k++ )
2999 double conversion_efficiency_HM79 = 1/(1. + 1e4*
sexp(920./qtemp[k]));
3000 sticking_probability_H = 1./(1. + 0.04*sqrt(qtemp[k]+
phycon.
te) +
3003 gv.
bin[nd]->rate_h2_form_grains_HM79 += qprob[k] * sticking_probability_H *
3004 conversion_efficiency_HM79;
3010 gv.
bin[nd]->rate_h2_form_grains_HM79 *= 0.5 * AveVelH *
3011 gv.
bin[nd]->IntArea/4. *
gv.
bin[nd]->cnv_H_pCM3;
3013 ASSERT(
gv.
bin[nd]->rate_h2_form_grains_HM79 > 0. );
3026 double conversion_efficiency_HM79 = 1/(1. + 1e4*
sexp(920./
gv.
bin[nd]->tedust));
3031 gv.
bin[nd]->rate_h2_form_grains_HM79 = 0.5 * AveVelH *
gv.
bin[nd]->IntArea/4. *
3033 gv.
bin[nd]->cnv_H_pCM3 * sticking_probability_H * conversion_efficiency_HM79;
3034 ASSERT(
gv.
bin[nd]->rate_h2_form_grains_HM79 > 0. );
3037 # ifndef IGNORE_QUANTUM_HEATING
3046 double sqrt_term = 35.399494936611667;
3048 gv.
bin[nd]->rate_h2_form_grains_CT02 = 0.;
3050 for( k=0; k < qnbin; k++ )
3052 double beta_alpha = 0.25 * sqrt_term *
sexp(200./qtemp[k] );
3054 double xi = 1./ (1. + 1.3e13*
sexp(1.5*1e4/qtemp[k])*sqrt_term/(2.*f) );
3056 double beta = 3e12 *
sexp( 320. / qtemp[k] );
3059 double recombination_efficiency_CT02 = xi / (1. + 0.005*f/2./
SDIV(beta) + beta_alpha );
3060 sticking_probability_H = 1./(1. + 0.04*sqrt(qtemp[k]+
phycon.
te) +
3066 gv.
bin[nd]->rate_h2_form_grains_CT02 += qprob[k] * sticking_probability_H *
3067 recombination_efficiency_CT02;
3074 gv.
bin[nd]->rate_h2_form_grains_CT02 *= 0.5 * AveVelH *
3075 gv.
bin[nd]->IntArea/4. *
gv.
bin[nd]->cnv_H_pCM3;
3077 ASSERT(
gv.
bin[nd]->rate_h2_form_grains_CT02 > 0. );
3091 double sqrt_term = 35.399494936611667;
3092 double beta_alpha = 0.25 * sqrt_term *
sexp(200./
gv.
bin[nd]->tedust );
3094 double xi = 1./ (1. + 1.3e13*
sexp(1.5*1e4/
gv.
bin[nd]->tedust )*sqrt_term/(2.*f) );
3096 double beta = 3e12 *
sexp( 320. /
gv.
bin[nd]->tedust );
3099 double recombination_efficiency_CT02 = xi / (1. + 0.005*f/2./
SDIV(beta) + beta_alpha );
3107 double Td =
gv.
bin[nd]->tedust;
3108 recombination_efficiency_CT02 = 1./(
3109 ( 1. + 4.609569629185726e24*
sexp(45000./Td) ) *
3110 ( 1. +
sexp(800./Td) / (0.5389970511202651 *
sexp(540./Td) + 5.6333909478365e-14*sqrt(Td)) )
3117 double Td =
gv.
bin[nd]->tedust;
3118 recombination_efficiency_CT02 = 1./(
3119 ( 1. + 6.998161265697586e24*
sexp(45000./Td) ) *
3120 ( 1. +
sexp(450./Td) / (0.4266153643741715 *
sexp(340./Td) + 2.5335919594255e-14*sqrt(Td)) )
3130 gv.
bin[nd]->rate_h2_form_grains_CT02 = 0.5 * AveVelH *
gv.
bin[nd]->IntArea/4. *
3131 gv.
bin[nd]->cnv_H_pCM3 * sticking_probability_H * recombination_efficiency_CT02;
3132 ASSERT(
gv.
bin[nd]->rate_h2_form_grains_CT02 > 0. );
3135 # ifndef IGNORE_QUANTUM_HEATING
3137 sticking_probability_H = 1./(1. + 0.04*sqrt(
gv.
bin[nd]->tedust+
phycon.
te) +
3152 gv.
bin[nd]->cnv_H_pCM3 * sticking_probability_H;
3178 if(
gv.
bin[nd]->tedust < T_ortho_para_crit )
3184 gv.
bin[nd]->cnv_H_pCM3 * sticking_probability_H * efficiency_opr;
3203 for(
size_t nd=0; nd <
gv.
bin.size(); nd++ )
3205 gv.
bin[nd]->rate_h2_form_grains_CT02 = 0.;
3206 gv.
bin[nd]->rate_h2_form_grains_HM79 = 0.;
3219 for(
size_t nd=0; nd <
gv.
bin.size(); nd++ )
3226 gv.
bin[nd]->rate_h2_form_grains_used =
3233 gv.
bin[nd]->rate_h2_form_grains_used =
3241 gv.
bin[nd]->rate_h2_form_grains_used =
3259 # ifndef IGNORE_QUANTUM_HEATING
3270 static double teused=-1;
3280 cr_H2dis_ELWERT_H2g,
3281 cr_H2dis_ELWERT_H2s;
3388 enum {DEBUG_LOC=
false};
3393 fprintf(
ioQQQ,
"ph2lteee\t%.2e\t%.1e\t%.1e\n",
3407 (*diatom)->Mol_Photo_Diss_Rates();
3439 enum {DEBUG_LOC=
false};
3457 enum {DEBUG_LOC=
false};
3459 if( DEBUG_LOC &&
nzone>400)
3568 double sqrtx = sqrt(1.+x);
3572 double fshield = 0.965/
POW2(1.+x/b5) + 0.035/sqrtx *
3624 double x_H2g, x_H2s,
3625 fshield_H2g, fshield_H2s,
3627 static double a_H2g, a_H2s,
3660 a_H2g = 0.06 * log_G0_face + 1.32;
3661 a_H2s = 0.375 * log_G0_face + 2.125;
3663 e1_H2g = -0.05 * log_G0_face + 2.25;
3664 e1_H2s = -0.125 * log_G0_face + 2.625;
3666 e2_H2g = -0.005 * log_G0_face + 0.625;
3668 b_H2g = -4.0e-3 * log_G0_face + 3.2e-2;
3675 k_f_H2s =
MAX2(0.1,2.375 * log_G0_face - 1.875 );
3678 k_H2g_to_H2s =
MAX2(1.,-1.75 * log_G0_face + 11.25);
3693 fshield_H2g = 0.965/pow(1.+x_H2g/b5,e1_H2g) + b_H2g/pow(1.+x_H2g/b5,e2_H2g);
3694 fshield_H2s = 0.965/pow(1.+x_H2s/b5,e1_H2s);
3858 enum {DEBUG_LOC=
false};
3862 fprintf(
ioQQQ,
" Solomon H2 dest rates: TH85 %.2e BD96 %.2e Big %.2e excit rates: TH85 %.2e Big %.2e\n",
3881 return &(*p->second);
3957 for(
int i=0;i<rate.
nreactants && ipthis == -1;i++)
3966 double ratevi = rate.
a * rate.
rk();
3990 fprintf(stderr,
"=>");
3994 fprintf(stderr,
"\n");
4028 int ipspin = 0, ipfreein = 0;
4036 int ipspout = 0, ipfreeout = 0;
4046 int newsp = ipspout-ipspin;
4051 int nbondsbroken = ipfreeout-ipfreein;
4052 if (nbondsbroken <= 0)
4055 double fracbroken = nbondsbroken/((double)ipfreeout);
4056 ASSERT( fracbroken <= 1.0 );
4067 double ratesp = ratevi*newsp;
4069 ratesp *= fracbroken;
4103 double ratevi = rate.
a * rate.
rk();
4108 ratev += ipthis*ratevi;
4130 double heating = 0.;
4131 map<double,string> heatMap;
4143 bool lgCanSkip =
false;
4171 realnum reaction_enthalpy = 0.;
4180 for(
long i=0; i < rate.
nproducts; ++i )
4190 double heat = reaction_enthalpy*rate_tot*(1e10/
AVOGADRO);
4191 heatMap[heat] = rate.
label;
4201 for( map<double,string>::reverse_iterator it = heatMap.rbegin(); it != heatMap.rend(); ++it, ++index )
4203 fprintf(
ioQQQ,
"DEBUGGG heat %li\t%li\t%.6e\t%s\n", index,
nzone, it->first, it->second.c_str() );
4208 for( map<double,string>::iterator it = heatMap.begin(); it != heatMap.end(); ++it, ++index )
4210 if( it->first >= 0. )
4212 fprintf(
ioQQQ,
"DEBUGGG cool %li\t%li\t%.6e\t%s\n", index,
nzone, it->first, it->second.c_str() );
4221 void mole_punch(FILE *punit,
const char speciesname[],
const char args[],
bool lgHeader,
bool lgData,
double depth)
4233 fprintf (punit,
"#Depth");
4237 fprintf (punit,
"%.5e",depth);
4252 ( strcmp( args,
"CATA" )==0 && rate.
rvector[i]!=NULL ) ||
4253 strcmp( args,
"ALL " )==0 )
4264 ( strcmp( args,
"CATA" )==0 && rate.
pvector[i]!=NULL ) ||
4265 strcmp( args,
"ALL " )==0 )
4274 fprintf(punit,
"\t%s",rate.
label.c_str());
4283 fprintf(punit,
"\t%.3e",ratevi);
4287 fprintf(punit,
"\n");
