48 STATIC void newreact(
const char label[],
const char fun[],
double a,
double b,
double c);
88 template<
class T>
void newfunc()
94 ratefunc findfunc(
const char name[]);
97 class mole_reaction_hmrate;
99 class mole_reaction_th85rate;
101 class mole_reaction_crnurate_noalbedo;
105 class mole_reaction_grn_react;
107 class mole_reaction_h2_collh_excit;
109 class mole_reaction_h2_collh2_excit;
111 class mole_reaction_h2_collh_deexc;
113 class mole_reaction_h2_collh2_deexc;
115 class mole_reaction_rh2g_dis_h;
117 class mole_reaction_rh2s_dis_h;
119 class mole_reaction_rh2g_dis_h2;
121 class mole_reaction_rh2s_dis_h2;
123 class mole_reaction_rh2s_dis_h2_nodeexcit;
125 class mole_reaction_bh2g_dis_h;
127 class mole_reaction_bh2s_dis_h;
129 class mole_reaction_bh2g_dis_h2;
131 class mole_reaction_hneut;
133 class mole_reaction_cionhm;
150 typedef mole_reaction_null T;
152 virtual T*
Create()
const {
return new T;}
153 virtual const char*
name() {
return "null";}
191 for(n=0;n<nreact;n++)
220 r *=
pow(te/300.,rate->
b);
222 r *= exp(-rate->
c/te);
228 typedef mole_reaction_hmrate_exo T;
230 virtual T*
Create()
const {
return new T;}
232 virtual const char*
name() {
return "hmrate_exo";}
243 te =
MIN2( te, -10. * this->c );
245 return pow(te/300.,this->b)*exp(-te/this->c);
251 typedef mole_reaction_hmrate T;
253 virtual T*
Create()
const {
return new T;}
254 virtual const char*
name() {
return "hmrate";}
263 typedef mole_reaction_constrate T;
265 virtual T*
Create()
const {
return new T;}
266 virtual const char*
name() {
return "constrate";}
307 typedef mole_reaction_th85rate T;
309 virtual T*
Create()
const {
return new T;}
310 virtual const char*
name() {
return "th85rate";}
313 return th85rate(
this);
339 typedef mole_reaction_crnurate_noalbedo T;
341 virtual T*
Create()
const {
return new T;}
343 virtual const char*
name() {
return "crnurate_noalbedo";}
346 return crnurate_noalbedo(
this);
352 typedef mole_reaction_crnurate T;
354 virtual T*
Create()
const {
return new T;}
355 virtual const char*
name() {
return "crnurate";}
358 return crnurate_noalbedo(
this)/(1.0-
albedo);
367 typedef mole_reaction_cryden_ov_bg T;
369 virtual T*
Create()
const {
return new T;}
370 virtual const char*
name() {
return "cryden_ov_bg";}
379 typedef mole_reaction_co_lnu_c_o_lnu T;
381 virtual T*
Create()
const {
return new T;}
382 virtual const char*
name() {
return "co_lnu_c_o_lnu";}
394 for( ns=0; ns<2; ++ns )
440 typedef mole_reaction_vib_evap T;
442 virtual T*
Create()
const {
return new T;}
443 virtual const char*
name() {
return "vib_evap";}
446 double binding_energy, exponent, vib_freq;
452 binding_energy = this->
b;
453 double bin_total=0.0;
454 for(
size_t nd=0; nd <
gv.
bin.size() ; nd++ )
456 double bin_area =
gv.
bin[nd]->IntArea*
gv.
bin[nd]->cnv_H_pCM3;
457 exponent += exp(-binding_energy/
gv.
bin[nd]->tedust)*bin_area;
458 bin_total += bin_area;
460 exponent /= bin_total;
461 const double surface_density_of_sites = 1.5e15;
464 vib_freq = sqrt(2*surface_density_of_sites*(0.3*BOLTZMANN)*binding_energy/(PI*PI*this->reactants[0]->mole_mass));
484 typedef mole_reaction_grn_abs T;
486 virtual T*
Create()
const {
return new T;}
487 virtual const char*
name() {
return "grn_abs";}
495 if (this->reactants[0]->n_nuclei() != 0)
496 mass = this->reactants[0]->mole_mass;
498 mass = this->reactants[1]->mole_mass;
500 return sqrt(8.*BOLTZMANN*
phycon.
te/(PI*mass));
506 typedef mole_reaction_grn_react T;
508 virtual T*
Create()
const {
return new T;}
509 virtual const char*
name() {
return "grn_react";}
512 return grn_react(
this);
536 double quant_barrier = 1e-8;
537 double surface_density_of_sites = 1.5e15;
538 fixit(
"rate->a is _always_ overall scaling factor");
542 double activ_barrier = rate->
c;
545 fixit(
"Ordering of reactants may have switched");
546 double vib_freq_i = sqrt(2*surface_density_of_sites*(0.3*BOLTZMANN)*E_i/(PI*PI*rate->
reactants[0]->
mole_mass));
547 double vib_freq_j = sqrt(2*surface_density_of_sites*(0.3*BOLTZMANN)*E_j/(PI*PI*rate->
reactants[1]->
mole_mass));
551 double dust_density = 0.0;
554 fixit(
"should there be an area weighting to exponential terms?");
556 for(
size_t nd=0; nd <
gv.
bin.size(); nd++ )
558 double bin_density =
gv.
bin[nd]->IntArea*
gv.
bin[nd]->cnv_H_pCM3;
559 Exp_i += exp(-E_i/
gv.
bin[nd]->tedust)*bin_density;
560 Exp_j += exp(-E_j/
gv.
bin[nd]->tedust)*bin_density;
561 dust_density += bin_density/(4*1e-10);
569 double diff_i = vib_freq_i*Exp_i/total_sites;
570 double diff_j = vib_freq_j*Exp_j/total_sites;
574 double scale = exp(-2*(quant_barrier/(HBAReV*EN1EV))*sqrt(2.0*rate->
reduced_mass*0.3*BOLTZMANN*activ_barrier));
577 return scale*(diff_i + diff_j)/
SDIV(dust_density);
582 typedef mole_reaction_grn_photo T;
584 virtual T*
Create()
const {
return new T;}
585 virtual const char*
name() {
return "grn_photo";}
600 fixit(
"grainn photoionization");
613 typedef mole_reaction_th85rate_co T;
615 virtual T*
Create()
const {
return new T;}
617 virtual const char*
name() {
return "th85rate_co";}
641 if (this->reactants[0]->n_nuclei() != 0)
642 column =
mole.
species[ this->reactants[0]->index ].column;
644 column =
mole.
species[ this->reactants[1]->index ].column;
646 esc_co = 4.4e-15 * column /
651 return esca0k2(esc_co)*th85rate(
this);
657 typedef mole_reaction_oh_c2h2_co_ch3 T;
659 virtual T*
Create()
const {
return new T;}
660 virtual const char*
name() {
return "oh_c2h2_co_ch3";}
679 typedef mole_reaction_h_hnc_hcn_h T;
681 virtual T*
Create()
const {
return new T;}
683 virtual const char*
name() {
return "h_hnc_hcn_h";}
719 return 1. - 4.938e-6;
724 double frac_H2star_grains(
void)
750 typedef mole_reaction_h2ph3p T;
752 virtual T*
Create()
const {
return new T;}
753 virtual const char*
name() {
return "h2ph3p";}
772 typedef mole_reaction_hopexch T;
774 virtual T*
Create()
const {
return new T;}
775 virtual const char*
name() {
return "hopexch";}
784 typedef mole_reaction_hpoexch T;
786 virtual T*
Create()
const {
return new T;}
787 virtual const char*
name() {
return "hpoexch";}
796 typedef mole_reaction_hmattach T;
798 virtual T*
Create()
const {
return new T;}
799 virtual const char*
name() {
return "hmattach";}
808 typedef mole_reaction_hmihph2p T;
810 virtual T*
Create()
const {
return new T;}
811 virtual const char*
name() {
return "hmihph2p";}
835 typedef mole_reaction_hmphoto T;
837 virtual T*
Create()
const {
return new T;}
838 virtual const char*
name() {
return "hmphoto";}
870 typedef mole_reaction_cionhm T;
872 virtual T*
Create()
const {
return new T;}
873 virtual const char*
name() {
return "cionhm";}
880 double assoc_detach(
void)
895 y=545969508.1323510+x*71239.23653059864;
901 typedef mole_reaction_c3bod T;
903 virtual T*
Create()
const {
return new T;}
904 virtual const char*
name() {
return "c3bod";}
913 typedef mole_reaction_asdfg T;
915 virtual T*
Create()
const {
return new T;}
916 virtual const char*
name() {
return "asdfg";}
925 typedef mole_reaction_asdfs T;
927 virtual T*
Create()
const {
return new T;}
928 virtual const char*
name() {
return "asdfs";}
937 typedef mole_reaction_asdbg T;
939 virtual T*
Create()
const {
return new T;}
940 virtual const char*
name() {
return "asdbg";}
950 typedef mole_reaction_asdbs T;
952 virtual T*
Create()
const {
return new T;}
953 virtual const char*
name() {
return "asdbs";}
963 typedef mole_reaction_bhneut T;
965 virtual T*
Create()
const {
return new T;}
966 virtual const char*
name() {
return "bhneut";}
977 return (hneut(
this)*ratio*
1000 return 3.4738192887404660e-008;
1005 typedef mole_reaction_hneut T;
1007 virtual T*
Create()
const {
return new T;}
1008 virtual const char*
name() {
return "hneut";}
1018 typedef mole_reaction_h2_spon_diss T;
1020 virtual T*
Create()
const {
return new T;}
1021 virtual const char*
name() {
return "h2_spon_diss";}
1030 typedef mole_reaction_h2_ind_rad_diss T;
1032 virtual T*
Create()
const {
return new T;}
1033 virtual const char*
name() {
return "h2_ind_rad_diss";}
1043 fixit(
"Remove factor of dense.gas_phase[ipHYDROGEN] factor from"
1044 "derivation of rate_h2_form_grains_used to avoid"
1054 virtual const char*
name() {
return "grnh2tot";}
1057 return grnh2tot(
this);
1063 typedef mole_reaction_grnh2 T;
1065 virtual T*
Create()
const {
return new T;}
1066 virtual const char*
name() {
return "grnh2";}
1069 return grnh2tot(
this)*(1.-frac_H2star_grains());
1075 typedef mole_reaction_grnh2s T;
1077 virtual T*
Create()
const {
return new T;}
1078 virtual const char*
name() {
return "grnh2s";}
1081 return grnh2tot(
this)*frac_H2star_grains();
1087 typedef mole_reaction_radasc T;
1089 virtual T*
Create()
const {
return new T;}
1090 virtual const char*
name() {
return "radasc";}
1098 return hmrate(
this) *
1110 typedef mole_reaction_assoc_ion T;
1112 virtual T*
Create()
const {
return new T;}
1113 virtual const char*
name() {
return "assoc_ion";}
1120 return hmrate(
this) *
1151 typedef mole_reaction_rh2g_dis_h T;
1153 virtual T*
Create()
const {
return new T;}
1154 virtual const char*
name() {
return "rh2g_dis_h";}
1157 return rh2g_dis_h(
this);
1172 fprintf(
ioQQQ,
"invalid parameter for rh2s_dis_h\n" );
1180 typedef mole_reaction_rh2s_dis_h T;
1182 virtual T*
Create()
const {
return new T;}
1183 virtual const char*
name() {
return "rh2s_dis_h";}
1186 return rh2s_dis_h(
this);
1202 fprintf(
ioQQQ,
"invalid parameter for rh2g_dis_h2\n" );
1205 return hmrate4(5.5e-29*0.5/(SAHA*3.634e-5)*sqrt(300.),0.5,5.195e4,
phycon.
te);
1210 typedef mole_reaction_rh2g_dis_h2 T;
1212 virtual T*
Create()
const {
return new T;}
1213 virtual const char*
name() {
return "rh2g_dis_h2";}
1216 return rh2g_dis_h2(
this);
1231 fprintf(
ioQQQ,
"invalid parameter for rh2s_dis_h2\n" );
1239 typedef mole_reaction_rh2s_dis_h2 T;
1241 virtual T*
Create()
const {
return new T;}
1242 virtual const char*
name() {
return "rh2s_dis_h2";}
1245 return rh2s_dis_h2(
this);
1260 fprintf(
ioQQQ,
"invalid parameter for rh2s_dis_h2_nodeexcit\n" );
1266 class mole_reaction_rh2s_dis_h2_nodeexcit :
public mole_reaction
1268 typedef mole_reaction_rh2s_dis_h2_nodeexcit T;
1270 virtual T*
Create()
const {
return new T;}
1271 virtual const char*
name() {
return "rh2s_dis_h2_nodeexcit";}
1274 return rh2s_dis_h2_nodeexcit(
this);
1284 typedef mole_reaction_bh2g_dis_h T;
1286 virtual T*
Create()
const {
return new T;}
1287 virtual const char*
name() {
return "bh2g_dis_h";}
1290 return bh2g_dis_h(
this);
1300 typedef mole_reaction_bh2s_dis_h T;
1302 virtual T*
Create()
const {
return new T;}
1303 virtual const char*
name() {
return "bh2s_dis_h";}
1306 return bh2s_dis_h(
this);
1316 typedef mole_reaction_bh2g_dis_h2 T;
1318 virtual T*
Create()
const {
return new T;}
1319 virtual const char*
name() {
return "bh2g_dis_h2";}
1322 return bh2g_dis_h2(
this);
1328 typedef mole_reaction_bh2s_dis_h2 T;
1330 virtual T*
Create()
const {
return new T;}
1331 virtual const char*
name() {
return "bh2s_dis_h2";}
1340 typedef mole_reaction_h2photon T;
1342 virtual T*
Create()
const {
return new T;}
1343 virtual const char*
name() {
return "h2photon";}
1352 typedef mole_reaction_h2crphot T;
1354 virtual T*
Create()
const {
return new T;}
1355 virtual const char*
name() {
return "h2crphot";}
1364 typedef mole_reaction_h2crphh T;
1366 virtual T*
Create()
const {
return new T;}
1367 virtual const char*
name() {
return "h2crphh";}
1388 typedef mole_reaction_h2scrphh T;
1390 virtual T*
Create()
const {
return new T;}
1391 virtual const char*
name() {
return "h2scrphh";}
1407 typedef mole_reaction_radath T;
1409 virtual T*
Create()
const {
return new T;}
1410 virtual const char*
name() {
return "radath";}
1419 typedef mole_reaction_gamtwo T;
1421 virtual T*
Create()
const {
return new T;}
1422 virtual const char*
name() {
return "gamtwo";}
1438 typedef mole_reaction_hlike_phot T;
1440 virtual T*
Create()
const {
return new T;}
1441 virtual const char*
name() {
return "hlike_phot";}
1457 typedef mole_reaction_h2s_sp_decay T;
1459 virtual T*
Create()
const {
return new T;}
1460 virtual const char*
name() {
return "h2s_sp_decay";}
1489 typedef mole_reaction_h2_collh2_deexc T;
1491 virtual T*
Create()
const {
return new T;}
1492 virtual const char*
name() {
return "h2_collh2_deexc";}
1495 return h2_collh2_deexc(
this);
1512 typedef mole_reaction_h2_collh_deexc T;
1514 virtual T*
Create()
const {
return new T;}
1515 virtual const char*
name() {
return "h2_collh_deexc";}
1518 return h2_collh_deexc(
this);
1536 typedef mole_reaction_h2_collh2_excit T;
1538 virtual T*
Create()
const {
return new T;}
1539 virtual const char*
name() {
return "h2_collh2_excit";}
1542 return h2_collh2_excit(
this);
1560 typedef mole_reaction_h2_collh_excit T;
1562 virtual T*
Create()
const {
return new T;}
1563 virtual const char*
name() {
return "h2_collh_excit";}
1566 return h2_collh_excit(
this);
1572 typedef mole_reaction_h2gexcit T;
1574 virtual T*
Create()
const {
return new T;}
1575 virtual const char*
name() {
return "h2gexcit";}
1584 typedef mole_reaction_h2sdissoc T;
1586 virtual T*
Create()
const {
return new T;}
1587 virtual const char*
name() {
return "h2sdissoc";};
1602 typedef mole_reaction_h2gdissoc T;
1604 virtual T*
Create()
const {
return new T;}
1605 virtual const char*
name() {
return "h2gdissoc";}
1620 typedef mole_reaction_hd_photodissoc T;
1622 virtual T*
Create()
const {
return new T;}
1623 virtual const char*
name() {
return "hd_photodissoc";}
1654 else if( te < 1200. )
1658 else if( te < 3800. )
1663 else if( te <= 1.4e4 )
1687 typedef mole_reaction_gamheh T;
1689 virtual T*
Create()
const {
return new T;}
1690 virtual const char*
name() {
return "gamheh";}
1702 for( i=
hmi.
iheh1-1; i < limit; i++ )
1715 enum {exclude, base, umisthack, federman, lithium, deuterium, ti, misc, in_code, generated};
1719 static bool lgReactInitialized =
false;
1729 if( lgReactInitialized )
1732 lgReactInitialized =
true;
1741 newfunc<mole_reaction_null>();
1743 map<string,bool> canonical;
1752 newfunc<mole_reaction_hmrate>();
1753 newfunc<mole_reaction_hmrate_exo>();
1754 newfunc<mole_reaction_constrate>();
1755 newfunc<mole_reaction_th85rate>();
1756 newfunc<mole_reaction_crnurate>();
1757 newfunc<mole_reaction_crnurate_noalbedo>();
1758 newfunc<mole_reaction_co_lnu_c_o_lnu>();
1759 newfunc<mole_reaction_vib_evap>();
1760 newfunc<mole_reaction_th85rate_co>();
1761 newfunc<mole_reaction_grn_abs>();
1763 newfunc<mole_reaction_grn_react>();
1765 newfunc<mole_reaction_grn_photo>();
1766 newfunc<mole_reaction_oh_c2h2_co_ch3>();
1767 newfunc<mole_reaction_h_hnc_hcn_h>();
1769 newfunc<mole_reaction_gamheh>();
1770 newfunc<mole_reaction_hd_photodissoc>();
1771 newfunc<mole_reaction_h2gdissoc>();
1772 newfunc<mole_reaction_h2sdissoc>();
1773 newfunc<mole_reaction_h2gexcit>();
1774 newfunc<mole_reaction_h2_collh_excit>();
1775 newfunc<mole_reaction_h2_collh2_excit>();
1776 newfunc<mole_reaction_h2_collh_deexc>();
1777 newfunc<mole_reaction_h2_collh2_deexc>();
1778 newfunc<mole_reaction_h2s_sp_decay>();
1779 newfunc<mole_reaction_hlike_phot>();
1780 newfunc<mole_reaction_gamtwo>();
1781 newfunc<mole_reaction_h2ph3p>();
1782 newfunc<mole_reaction_radath>();
1783 newfunc<mole_reaction_cryden_ov_bg>();
1784 newfunc<mole_reaction_h2scrphh>();
1785 newfunc<mole_reaction_h2crphh>();
1786 newfunc<mole_reaction_h2photon>();
1787 newfunc<mole_reaction_h2crphot>();
1788 newfunc<mole_reaction_rh2g_dis_h>();
1789 newfunc<mole_reaction_rh2s_dis_h>();
1790 newfunc<mole_reaction_rh2g_dis_h2>();
1791 newfunc<mole_reaction_rh2s_dis_h2>();
1792 newfunc<mole_reaction_rh2s_dis_h2_nodeexcit>();
1793 newfunc<mole_reaction_bh2g_dis_h>();
1794 newfunc<mole_reaction_bh2s_dis_h>();
1795 newfunc<mole_reaction_bh2g_dis_h2>();
1796 newfunc<mole_reaction_bh2s_dis_h2>();
1797 newfunc<mole_reaction_radasc>();
1798 newfunc<mole_reaction_assoc_ion>();
1799 newfunc<mole_reaction_grnh2>();
1800 newfunc<mole_reaction_grnh2s>();
1801 newfunc<mole_reaction_h2_spon_diss>();
1802 newfunc<mole_reaction_h2_ind_rad_diss>();
1803 newfunc<mole_reaction_bhneut>();
1804 newfunc<mole_reaction_hneut>();
1805 newfunc<mole_reaction_asdbs>();
1806 newfunc<mole_reaction_asdbg>();
1807 newfunc<mole_reaction_asdfs>();
1808 newfunc<mole_reaction_asdfg>();
1809 newfunc<mole_reaction_c3bod>();
1810 newfunc<mole_reaction_cionhm>();
1811 newfunc<mole_reaction_hmphoto>();
1812 newfunc<mole_reaction_hmihph2p>();
1813 newfunc<mole_reaction_hmattach>();
1848 newreact(
"C+,OH=>CO,H+",
"hmrate",0.,0.,0.);
1853 newreact(
"H,H,grn=>H2,grn",
"grnh2",1.,0.,0.);
1854 newreact(
"H,H,grn=>H2*,grn",
"grnh2s",1.,0.,0.);
1855 newreact(
"H-,PHOTON=>H,e-",
"hmphoto",1.,0.,0.);
1860 fixit(
"this should be atom_list instead of unresolved_element_list, but we have not defined ionized species of all isotopes yet!!!");
1864 if( !(*atom)->lgHasLinkedIon() )
1866 long nelem = (*atom)->el->Z-1;
1870 sprintf(react,
"H-,%s+=>H,%s", (*atom)->label().c_str(), (*atom)->label().c_str() );
1871 newreact(react,
"hmrate",4e-6/sqrt(300.),-0.5,0.);
1875 newreact(
"H,e-=>H-,PHOTON",
"hmattach",1.,0.,0.);
1876 newreact(
"H-,H+=>H2+,e-",
"hmihph2p",1.,0.,0.);
1877 newreact(
"H-,e-=>H,e-,e-",
"cionhm",1.,0.,0.);
1878 newreact(
"H,e-,e-=>H-,e-",
"c3bod",1.,0.,0.);
1879 newreact(
"H,H-=>H2,e-",
"asdfg",1.,0.,0.);
1880 newreact(
"H,H-=>H2*,e-",
"asdfs",1.,0.,0.);
1881 newreact(
"H2,e-=>H,H-",
"asdbg",1.,0.,0.);
1882 newreact(
"H2*,e-=>H,H-",
"asdbs",1.,0.,0.);
1883 newreact(
"H-,H+=>H,H",
"hneut",1.,0.,0.);
1884 newreact(
"H,H=>H-,H+",
"bhneut",1.,0.,0.);
1889 newreact(
"H2*,PHOTON=>H,H,PHOTON",
"h2_ind_rad_diss",1.,0.,0.);
1890 newreact(
"H,H2+=>H+,H2",
"hmrate",0.,0.,0.);
1895 newreact(
"H2,PHOTON=>H,H,PHOTON",
"h2_ind_rad_diss",1.,0.,0.);
1896 newreact(
"H,H2+=>H+,H2",
"hmrate",6.4e-10f,0.,0.);
1901 newreact(
"H,H=>H2,PHOTON",
"radasc",3e-14,0.,0.);
1902 newreact(
"H,H=>H2+,e-",
"assoc_ion",3.27e-10,-0.35,17829.);
1903 newreact(
"H2,H=>H,H,H",
"rh2g_dis_h",1.,0.,0.);
1904 newreact(
"H2,H2=>H,H,H2",
"rh2g_dis_h2",1.,0.,0.);
1905 newreact(
"H,H,H=>H2,H",
"bh2g_dis_h",1.,0.,0.);
1906 newreact(
"H,H,H2=>H2,H2",
"bh2g_dis_h2",1.,0.,0.);
1908 newreact(
"H2,PHOTON=>H2+,e-",
"h2photon",1.,0.,0.);
1909 newreact(
"H2,CRPHOT=>H2+,e-",
"h2crphot",1.,0.,0.);
1910 newreact(
"H2*,PHOTON=>H2+,e-",
"h2photon",1.,0.,0.);
1911 newreact(
"H2*,CRPHOT=>H2+,e-",
"h2crphot",1.,0.,0.);
1912 newreact(
"H2,CRPHOT=>H,H",
"h2crphh",1.,0.,0.);
1913 newreact(
"H2,CRPHOT=>H+,H-",
"cryden_ov_bg",3.9e-21,0.,0.);
1914 newreact(
"H2*,CRPHOT=>H+,H-",
"cryden_ov_bg",3.9e-21,0.,0.);
1915 newreact(
"H2,CRPHOT=>H+,H,e-",
"cryden_ov_bg",2.2e-19,0.,0.);
1916 newreact(
"H2*,CRPHOT=>H+,H,e-",
"cryden_ov_bg",2.2e-19,0.,0.);
1917 newreact(
"H+,H=>H2+,PHOTON",
"radath",1.,0.,0.);
1918 newreact(
"H2+,PHOTON=>H,H+",
"gamtwo",1.,0.,0.);
1919 newreact(
"H2+,CRPHOT=>H,H+",
"hlike_phot",1.,0.,0.);
1922 newreact(
"H3+,CRPHOT=>H2+,H+,e-",
"hlike_phot",2.,0.,0.);
1923 newreact(
"H,H+,e-=>H2+,e-",
"hmrate",2e-7 * SAHA * (4./(2.*1.)) * 3.634e-5 *
pow(300.,-1.50),-1.50,0.);
1927 newreact(
"H2,CRPHOT=>H2+,e-",
"hmrate",4.4e-17,0.,0.);
1928 newreact(
"H2,CRPHOT=>H,H+,e-",
"hmrate",1e-19,0.,0.);
1929 newreact(
"H2*,CRPHOT=>H,H+,e-",
"hmrate",1e-19,0.,0.);
1930 newreact(
"H2*,CRPHOT=>H2+,e-",
"hmrate",4.4e-17,0.,0.);
1931 newreact(
"H2,CRPHOT=>H,H",
"hmrate",5e-18,0.,0.);
1932 newreact(
"e-,H3+=>H2,H",
"hmrate",2.5e-8,-0.3,0.);
1933 newreact(
"e-,H3+=>H,H,H",
"hmrate",7.5e-8,-0.3,0.);
1934 newreact(
"H+,H=>H2+,PHOTON",
"hmrate",5.3e-19,1.85,0);
1938 newreact(
"H2+,e-=>H,H",
"hmrate",1.6e-8,-0.43,0.);
1939 newreact(
"H2+,PHOTON=>H,H+",
"th85rate",5.7e-10,0.,1.9);
1942 newreact(
"H2*,CRPHOT=>H,H",
"h2scrphh",1.,0.,0.);
1943 newreact(
"H2,H2+=>H,H3+",
"h2ph3p",1.,0.,0.);
1944 newreact(
"H2*=>H2,PHOTON",
"h2s_sp_decay",1.,0.,0.);
1945 newreact(
"H2*,H2=>H2,H2",
"h2_collh2_deexc",1.,0.,0.);
1946 newreact(
"H2*,H=>H2,H",
"h2_collh_deexc",1.,0.,0.);
1947 newreact(
"H2,H2=>H2*,H2",
"h2_collh2_excit",1.,0.,0.);
1948 newreact(
"H2,H=>H2*,H",
"h2_collh_excit",1.,0.,0.);
1951 newreact(
"H2*,H=>H,H,H",
"rh2s_dis_h",1.,0.,0.);
1952 newreact(
"H2*,H2=>H2,H,H",
"rh2s_dis_h2",1.,0.,0.);
1953 newreact(
"H2*,H2*=>H2,H,H",
"rh2s_dis_h2",1.,0.,0.);
1954 newreact(
"H2*,H2*=>H2*,H,H",
"rh2s_dis_h2_nodeexcit",1.,0.,0.);
1958 newreact(
"H2,He=>He,H,H",
"rh2g_dis_h",1.,0.,0.);
1959 newreact(
"H2,H+=>H+,H,H",
"rh2g_dis_h",1.,0.,0.);
1960 newreact(
"H2,H3+=>H3+,H,H",
"rh2g_dis_h",1.,0.,0.);
1961 newreact(
"H2,e-=>e-,H,H",
"rh2g_dis_h",1.,0.,0.);
1962 newreact(
"H2*,He=>He,H,H",
"rh2s_dis_h",1.,0.,0.);
1963 newreact(
"H2*,H+=>H+,H,H",
"rh2s_dis_h",1.,0.,0.);
1964 newreact(
"H2*,H3+=>H3+,H,H",
"rh2s_dis_h",1.,0.,0.);
1965 newreact(
"H2*,e-=>e-,H,H",
"rh2s_dis_h",1.,0.,0.);
1968 newreact(
"He,H,H=>H2,He",
"bh2g_dis_h",1.,0.,0.);
1969 newreact(
"H+,H,H=>H2,H+",
"bh2g_dis_h",1.,0.,0.);
1970 newreact(
"H3+,H,H=>H2,H3+",
"bh2g_dis_h",1.,0.,0.);
1971 newreact(
"e-,H,H=>H2,e-",
"bh2g_dis_h",1.,0.,0.);
1972 newreact(
"H,H,H=>H2*,H",
"bh2s_dis_h",1.,0.,0.);
1973 newreact(
"H2,H,H=>H2*,H2",
"bh2s_dis_h2",1.,0.,0.);
1974 newreact(
"H2,H,H=>H2*,H2*",
"bh2s_dis_h2",1.,0.,0.);
1975 newreact(
"H2*,H,H=>H2*,H2*",
"bh2s_dis_h2",1.,0.,0.);
1976 newreact(
"He,H,H=>H2*,He",
"bh2s_dis_h",1.,0.,0.);
1977 newreact(
"H+,H,H=>H2*,H+",
"bh2s_dis_h",1.,0.,0.);
1978 newreact(
"H3+,H,H=>H2*,H3+",
"bh2s_dis_h",1.,0.,0.);
1979 newreact(
"e-,H,H=>H2*,e-",
"bh2s_dis_h",1.,0.,0.);
1982 newreact(
"H2,PHOTON=>H2*",
"h2gexcit",1.,0.,0.);
1983 newreact(
"H2*,PHOTON=>H,H",
"h2sdissoc",1.,0.,0.);
1984 newreact(
"H2,PHOTON=>H,H",
"h2gdissoc",1.,0.,0.);
1985 newreact(
"HeH+,PHOTON=>H,He+",
"gamheh",1.,0.,0.);
1990 newreact(
"OHgrn,Hgrn=>H2Ogrn",
"grn_react",1.,0.,0.);
1995 fprintf(stderr,
"Finished testing against UDFA database\n");
2011 it->second->index = index;
2025 DEBUG_ENTRY(
"mole_generate_isotopologue_reactions()" );
2027 bool lgDebug =
false;
2031 vector<string> newReactionStrings;
2032 vector<mole_reaction*> oldReactions;
2033 vector<realnum> branchingRatios;
2038 bool lgFound =
false;
2042 for(
long j=0; j<it->second->nproducts; ++j )
2045 if( it->second->products[j]->nNuclide.find( atomOld ) != it->second->products[j]->nNuclide.end() )
2046 numSites += it->second->products[j]->nNuclide[ atomOld ];
2049 for(
long i=0; i<it->second->nreactants; ++i )
2052 for(
nNucs_i atom=it->second->reactants[i]->nNuclide.begin(); atom != it->second->reactants[i]->nNuclide.end(); ++atom )
2054 if( atom->first->label() == atom_old )
2068 fprintf(
ioQQQ,
"DEBUGGG mole_generate_isotopologue_reactions %s ..........\n", it->first.c_str() );
2070 for(
long i=0; i<it->second->nreactants; ++i )
2073 if( it->second->reactants[i]->mole_mass < 0.99 * ATOMIC_MASS_UNIT )
2076 vector<string> react_iso_labels;
2078 vector< int > numAtoms;
2079 string embellishments;
2081 bool lgParseOK =
parse_species_label( it->second->reactants[i]->label.c_str(), atomsLeftToRight, numAtoms, embellishments );
2092 for(
unsigned j=0; j<react_iso_labels.size(); ++j )
2094 int numAtomsTot = 0;
2095 for(
long k=0; k<it->second->nproducts; ++k )
2098 if( it->second->products[k]->mole_mass < 0.99 * ATOMIC_MASS_UNIT )
2101 atomsLeftToRight.clear();
2103 embellishments.clear();
2105 lgParseOK =
parse_species_label( it->second->products[k]->label.c_str(), atomsLeftToRight, numAtoms, embellishments );
2106 ASSERT( lgParseOK ==
true );
2108 for(
unsigned position = 0; position < atomsLeftToRight.size(); ++position )
2110 string prod_iso_label;
2121 if( prod_iso_label.empty() )
2125 string react_string;
2127 for(
long i1=0; i1<i; ++i1 )
2129 react_string += it->second->reactants[i1]->label;
2130 if( i1 != it->second->nreactants-1 )
2131 react_string +=
",";
2133 react_string += react_iso_labels[j];
2134 if( i != it->second->nreactants-1 )
2135 react_string +=
",";
2136 for(
long i2=i+1; i2<it->second->nreactants; ++i2 )
2138 react_string += it->second->reactants[i2]->label;
2139 if( i2 != it->second->nreactants-1 )
2140 react_string +=
",";
2143 react_string +=
"=>";
2145 for(
long k1=0; k1<k; ++k1 )
2147 react_string += it->second->products[k1]->label;
2148 if( k1 != it->second->nproducts-1 )
2149 react_string +=
",";
2151 react_string += prod_iso_label;
2152 if( k != it->second->nproducts-1 )
2153 react_string +=
",";
2154 for(
long k2=k+1; k2<it->second->nproducts; ++k2 )
2156 react_string += it->second->products[k2]->label;
2157 if( k2 != it->second->nproducts-1 )
2158 react_string +=
",";
2163 newReactionStrings.push_back( canon_react_string );
2164 oldReactions.push_back( it->second.get_ptr() );
2166 branchingRatios.push_back( numAtoms[position]/numSites );
2169 fprintf(
ioQQQ,
"DEBUGGG mole_generate_isotopologue_reactions .................... %s\t\t(%2i/%2i).\n",
2170 canon_react_string.c_str(), numAtoms[position], numSites );
2172 numAtomsTot += numAtoms[position];
2175 ASSERT( numAtomsTot == numSites );
2180 ASSERT( oldReactions.size() == newReactionStrings.size() );
2183 vector<mole_reaction*>::const_iterator it2 = oldReactions.begin();
2184 vector<realnum>::const_iterator it3 = branchingRatios.begin();
2185 for( vector<string>::const_iterator it1 = newReactionStrings.begin(); it1 != newReactionStrings.end(); ++it1, ++it2, ++it3 )
2187 fixit(
"make adjustments to a for mass?");
2192 ASSERT( *it3 <= 1. + FLT_EPSILON );
2193 fixit(
"multiply by *it3 here. ASSERT at mole_reactions.cpp:1190 will blow currently.");
2194 newreact( it1->c_str(), (*it2)->name(), (*it2)->a , (*it2)->b, (*it2)->c );
2205 char chLabel[50], chLabelSave[50];
2212 strcpy( chLabel, p->second->label.c_str() );
2213 strcpy( chLabelSave, p->second->label.c_str() );
2214 char *str = chLabel;
2215 const char *delim =
"=>";
2216 char *chNewProducts = strtok( str, delim );
2217 char *chNewReactants = strtok( NULL, delim );
2218 char chNewReaction[50];
2220 strcpy( chNewReaction, chNewReactants );
2221 strcat( chNewReaction,
"=>" );
2222 strcat( chNewReaction, chNewProducts );
2232 fprintf(
ioQQQ,
"Warning! No reverse reaction for %30s. ", p->second->label.c_str() );
2236 fixit(
"NB reverse reactions should be generated here");
2245 DEBUG_ENTRY(
"mole_get_equilibrium_condition()" );
2254 DEBUG_ENTRY(
"mole_get_equilibrium_condition()" );
2260 double ln_result = 0.;
2266 ln_result += log(fac);
2273 ln_result -= log(fac);
2286 if( sp->
label ==
"PHOTON" || sp->
label ==
"CRPHOT" )
2288 fixit(
"How can we adapt existing structures to have a photon energy or range available here?");
2289 fixit(
"include 2hnu^3/c^2. Understand HNU3C2 macro!");
2292 else if( sp->
label ==
"CRP" || sp->
label ==
"grn" )
2295 fixit(
"need to figure out stat weight for any given particle");
2298 double deltaH = sp->
form_enthalpy * (1e10/AVOGADRO/BOLTZMANN);
2302 ASSERT( part_fun >= 0. );
2311 vector<long> numForm, numDest;
2329 for(
unsigned i=0; i<numForm.size(); ++i )
2331 if( numForm[i]==0 && numDest[i]==0 )
2333 if( numForm[i]>1 && numDest[i]>1 )
2346 char *label, *reactstr, *f;
2361 STATIC void newreact(
const char label[],
const char fun[],
double a,
double b,
double c)
2365 ratefunc rate = findfunc(fun);
2366 if (rate.get_ptr() == NULL)
2368 fprintf(stderr,
"Rate function %s not known for reaction %s. Aborting. Sorry.\n",fun,label);
2374 fprintf(
ioQQQ,
"\n PROBLEM Reaction %s has negative pre-coefficient. Aborting. Sorry.\n", label);
2378 rate->label = label;
2382 rate->source = source;
2397 fprintf(
ioQQQ,
" W-reaction %s disabled\n",rate->label.c_str());
2402 const char *rateLabelPtr = rate->label.c_str();
2406 rate->udfastate =
ABSENT;
2410 if (rate->udfastate ==
ABSENT)
2412 fprintf(stderr,
"Reaction %s not in UDFA\n",rateLabelPtr);
2417 if (rate->nreactants == 2 && rate->reactants[0]->mole_mass!=0.0 && rate->reactants[1]->mole_mass!=0.0 )
2419 rate->reduced_mass = 1./(1./rate->reactants[0]->mole_mass+1./rate->reactants[1]->mole_mass);
2423 rate->reduced_mass = 0.;
2434 fprintf(
ioQQQ,
"Warning: duplicate reaction %s -- using new version\n",rateLabelPtr);
2455 bool lgProd =
false;
2457 for(
int i=0;!i || label[i-1]!=
'\0';i++)
2459 if(label[i] ==
',' || label[i] ==
'=' || label[i] ==
'\0')
2465 fprintf(
ioQQQ,
"Mole_reactions: ignoring reaction %s (species %s not active)\n",label,buf.c_str());
2473 fprintf(stderr,
"Mole_reactions: Too many reactants in %s, only %d allowed\n",label,MAXREACTANTS);
2483 fprintf(stderr,
"Mole_reactions: Too many products in %s, only %d allowed\n",label,MAXPRODUCTS);
2492 if (label[i] !=
'>')
2524 ratefunc rate = findfunc(
"null");
2525 rate->label = label;
2562 rate->
label = newLabel;
2572 bool lgTrivial =
false;
2717 fprintf(sparsefp,
"P4\n%d %d\n",
2725 ch = (ch << 1) | (c[i][j] != 0.0);
2748 int dcharge, n,
sign;
2755 nel[it->first] += it->second;
2761 nel[it->first] -= it->second;
2766 fprintf(stderr,
"Reaction %s charge out of balance by %d\n",
2767 rate->
label.c_str(),dcharge);
2771 for(
nNucs_i it = nel.begin(); it != nel.end(); ++it )
2779 fprintf(stderr,
"Error: reaction %s %s %d of element %s\n",
2780 rate->
label.c_str(),sign==1?
"destroys":
"creates",
2782 it->first->label().c_str() );
2794 class formula_species {
2799 bool operator< (
const formula_species &a,
const formula_species &b)
2804 if (a.reactants[i]<b.reactants[i])
2806 else if (a.reactants[i]>b.reactants[i])
2811 if (a.products[i]<b.products[i])
2813 else if (a.products[i]>b.products[i])
2819 class udfa_reaction {
2824 double a, b, c, trange[2];
2828 static map <formula_species,count_ptr<udfa_reaction> >
udfatab;
2838 fprintf(stderr,
"Error, could not read %s\n",file);
2842 fixit(
"this seg-faults if file ends in blank line!");
2851 #define FLTEQ(a,b) (fabs((a)-(b)) <= 1e-6*fabs((a)+(b)))
2855 unsigned int havespecies = 1, i, n;
2867 r->l.reactants[n] = NULL;
2879 if (!strncmp(f,
"CRPHOT",6))
2888 r->l.products[n] = NULL;
2926 r->source = f[0]?f[0]:
'?';
2929 r->trange[n] = atof(f);
2934 if (udfatab.find(r->l) != udfatab.end())
2936 fprintf(stderr,
"Duplicate reaction\n");
2965 map<formula_species, count_ptr<udfa_reaction> >::iterator p
2967 if (p == udfatab.end() )
2991 ratefunc findfunc (
const char name[])
2999 enum { DEBUG_MOLE =
false };
3011 long index = rate.
index;
3016 if (fabs(newrk-oldrk) > 0.1*newrk)
3018 rate.
label.c_str(),oldrk,newrk);
3026 enum { DEBUG_MOLE =
false };
3037 double bigchange = 0.;
3038 unsigned long bigindex = ULONG_MAX;
3043 sum = oldrk+newrk, diff = newrk-oldrk;
3046 double change = fabs(diff)/sum;
3047 if (change > bigchange)
3059 if (bigindex == (
unsigned) rate.
index)
3065 pct = 100.*(newrk-oldrk)/oldrk;
3066 fprintf(
ioQQQ,
"Zone %ld, big chemistry rate change %s:"
3067 " %15.8g => %15.8g (%6.2g%%)\n",
3085 double T_ortho_para_crit,xi_ELRD, beta_alpha_ELRD, recombination_efficiency_ELRD, Td;
3097 fixit(
"Is this still necessary?");
3110 for(
size_t nd=0; nd <
gv.
bin.size(); nd++ )
3120 vector<double> qtemp, qprob;
3127 qheat(qtemp,qprob,&qnbin,nd);
3129 if(
gv.
bin[nd]->lgUseQHeat )
3137 qtemp[0] =
gv.
bin[nd]->tedust;
3140 gv.
bin[nd]->rate_h2_form_grains_HM79 = 0.;
3142 for( k=0; k < qnbin; k++ )
3152 double conversion_efficiency_HM79 = 1/(1. + 1e4*
sexp(920./qtemp[k]));
3155 gv.
bin[nd]->rate_h2_form_grains_HM79 += qprob[k] * sticking_probability_H *
3156 conversion_efficiency_HM79;
3162 gv.
bin[nd]->rate_h2_form_grains_HM79 *= 0.5 * AveVelH *
3163 gv.
bin[nd]->IntArea/4. *
gv.
bin[nd]->cnv_H_pCM3;
3165 ASSERT(
gv.
bin[nd]->rate_h2_form_grains_HM79 > 0. );
3177 double conversion_efficiency_HM79 = 1/(1. + 1e4*
sexp(920./
gv.
bin[nd]->tedust));
3182 gv.
bin[nd]->rate_h2_form_grains_HM79 = 0.5 * AveVelH *
gv.
bin[nd]->IntArea/4. *
3184 gv.
bin[nd]->cnv_H_pCM3 * sticking_probability_H * conversion_efficiency_HM79;
3185 ASSERT(
gv.
bin[nd]->rate_h2_form_grains_HM79 > 0. );
3196 double sqrt_term = 35.399494936611667;
3198 gv.
bin[nd]->rate_h2_form_grains_CT02 = 0.;
3200 for( k=0; k < qnbin; k++ )
3202 double beta_alpha = 0.25 * sqrt_term *
sexp(200./qtemp[k] );
3204 double xi = 1./ (1. + 1.3e13*
sexp(1.5*1e4/qtemp[k])*sqrt_term/(2.*f) );
3206 double beta = 3e12 *
sexp( 320. / qtemp[k] );
3209 double recombination_efficiency_CT02 = xi / (1. + 0.005*f/2./
SDIV(beta) + beta_alpha );
3215 gv.
bin[nd]->rate_h2_form_grains_CT02 += qprob[k] * sticking_probability_H *
3216 recombination_efficiency_CT02;
3223 gv.
bin[nd]->rate_h2_form_grains_CT02 *= 0.5 * AveVelH *
3224 gv.
bin[nd]->IntArea/4. *
gv.
bin[nd]->cnv_H_pCM3;
3226 ASSERT(
gv.
bin[nd]->rate_h2_form_grains_CT02 > 0. );
3236 double sqrt_term = 35.399494936611667;
3237 double beta_alpha = 0.25 * sqrt_term * exp(-200./
gv.
bin[nd]->tedust);
3239 double xi = 1./ (1. + 1.3e13*exp(-1.5e4/
gv.
bin[nd]->tedust)*sqrt_term/(2.*f) );
3241 double beta = 3e12 * exp(-320./
gv.
bin[nd]->tedust);
3244 double recombination_efficiency_CT02 = beta*xi / (beta + 0.0025*f + beta*beta_alpha);
3250 gv.
bin[nd]->rate_h2_form_grains_CT02 = 0.5 * AveVelH *
gv.
bin[nd]->IntArea/4. *
3251 gv.
bin[nd]->cnv_H_pCM3 * sticking_probability_H * recombination_efficiency_CT02;
3252 ASSERT(
gv.
bin[nd]->rate_h2_form_grains_CT02 >= 0. );
3265 gv.
bin[nd]->rate_h2_form_grains_ELRD= 0.;
3271 for( k=0; k < qnbin; k++ )
3275 beta_alpha_ELRD = exp(-800./Td) / (0.5389970511202651 * exp(-540./Td) + 5.6333909478365e-14*sqrt(Td) ) ;
3276 xi_ELRD= 1./(1. + 4.61e24*
sexp(45000./Td));
3277 recombination_efficiency_ELRD= (1. / (1. + beta_alpha_ELRD))*xi_ELRD;
3278 sticking_probability_H = 1./(1. + 0.04*sqrt(Td+
phycon.
te) +
3281 gv.
bin[nd]->rate_h2_form_grains_ELRD+= qprob[k] * sticking_probability_H *
3282 recombination_efficiency_ELRD;
3290 gv.
bin[nd]->rate_h2_form_grains_ELRD*= 0.5 * AveVelH *
3291 gv.
bin[nd]->IntArea/4. *
gv.
bin[nd]->cnv_H_pCM3;
3293 ASSERT(
gv.
bin[nd]->rate_h2_form_grains_ELRD> 0. );
3298 for( k=0; k < qnbin; k++ )
3302 beta_alpha_ELRD = exp(-450./Td) / (0.4266153643741715*exp(-340./Td) + 2.5335919594255e-14*sqrt(Td) ) ;
3303 xi_ELRD= 1./(1. + 7.00e24*
sexp(45000./Td));
3304 recombination_efficiency_ELRD= (1. / (1. + beta_alpha_ELRD))*xi_ELRD;
3305 sticking_probability_H = 1./(1. + 0.04*sqrt(Td+
phycon.
te) +
3308 gv.
bin[nd]->rate_h2_form_grains_ELRD+= qprob[k] * sticking_probability_H *
3309 recombination_efficiency_ELRD;
3316 gv.
bin[nd]->rate_h2_form_grains_ELRD*= 0.5 * AveVelH *
3317 gv.
bin[nd]->IntArea/4. *
gv.
bin[nd]->cnv_H_pCM3;
3319 ASSERT(
gv.
bin[nd]->rate_h2_form_grains_ELRD> 0. );
3329 gv.
bin[nd]->rate_h2_form_grains_ELRD= 0.;
3335 Td =
gv.
bin[nd]->tedust;
3337 beta_alpha_ELRD = exp(-800./Td) / (0.5389970511202651 * exp(-540./Td) + 5.6333909478365e-14*sqrt(Td) ) ;
3338 xi_ELRD= 1./(1. + 4.61e24*
sexp(45000./Td));
3339 recombination_efficiency_ELRD = (1. / (1. + beta_alpha_ELRD))*xi_ELRD;
3345 gv.
bin[nd]->rate_h2_form_grains_ELRD = 0.5 * AveVelH *
gv.
bin[nd]->IntArea/4. *
3346 gv.
bin[nd]->cnv_H_pCM3 * sticking_probability_H * recombination_efficiency_ELRD;
3348 ASSERT(
gv.
bin[nd]->rate_h2_form_grains_ELRD > 0. );
3353 Td =
gv.
bin[nd]->tedust;
3355 beta_alpha_ELRD = exp(-450./Td) / (0.4266153643741715*exp(-340./Td) + 2.5335919594255e-14*sqrt(Td) ) ;
3356 xi_ELRD= 1./(1. + 7.00e24*
sexp(45000./Td));
3357 recombination_efficiency_ELRD = (1. / (1. + beta_alpha_ELRD))*xi_ELRD;
3363 gv.
bin[nd]->rate_h2_form_grains_ELRD = 0.5 * AveVelH *
gv.
bin[nd]->IntArea/4. *
3364 gv.
bin[nd]->cnv_H_pCM3 * sticking_probability_H * recombination_efficiency_ELRD;
3366 ASSERT(
gv.
bin[nd]->rate_h2_form_grains_ELRD > 0. );
3387 gv.
bin[nd]->cnv_H_pCM3 * sticking_probability_H;
3413 if(
gv.
bin[nd]->tedust < T_ortho_para_crit )
3419 gv.
bin[nd]->cnv_H_pCM3 * sticking_probability_H * efficiency_opr;
3438 for(
size_t nd=0; nd <
gv.
bin.size(); nd++ )
3440 gv.
bin[nd]->rate_h2_form_grains_CT02 = 0.;
3441 gv.
bin[nd]->rate_h2_form_grains_HM79 = 0.;
3454 for(
size_t nd=0; nd <
gv.
bin.size(); nd++ )
3461 gv.
bin[nd]->rate_h2_form_grains_used =
3471 gv.
bin[nd]->rate_h2_form_grains_used =
3479 gv.
bin[nd]->rate_h2_form_grains_used =
3487 gv.
bin[nd]->rate_h2_form_grains_used =
3517 static double teused=-1;
3527 cr_H2dis_ELWERT_H2g,
3528 cr_H2dis_ELWERT_H2s;
3635 enum {DEBUG_LOC=
false};
3650 fixit(
"Wasted cycles if we don't use Stancil's rates below? Why not put this down there/if used?");
3654 (*diatom)->Mol_Photo_Diss_Rates();
3686 enum {DEBUG_LOC=
false};
3704 enum {DEBUG_LOC=
false};
3706 if( DEBUG_LOC &&
nzone>400)
3815 double sqrtx = sqrt(1.+x);
3819 double fshield = 0.965/
POW2(1.+x/b5) + 0.035/sqrtx *
3871 double x_H2g, x_H2s,
3872 fshield_H2g, fshield_H2s,
3874 static double a_H2g, a_H2s,
3907 a_H2g = 0.06 * log_G0_face + 1.32;
3908 a_H2s = 0.375 * log_G0_face + 2.125;
3910 e1_H2g = -0.05 * log_G0_face + 2.25;
3911 e1_H2s = -0.125 * log_G0_face + 2.625;
3913 e2_H2g = -0.005 * log_G0_face + 0.625;
3915 b_H2g = -4.0e-3 * log_G0_face + 3.2e-2;
3922 k_f_H2s =
MAX2(0.1,2.375 * log_G0_face - 1.875 );
3925 k_H2g_to_H2s =
MAX2(1.,-1.75 * log_G0_face + 11.25);
3940 fshield_H2g = 0.965/
pow(1.+x_H2g/b5,e1_H2g) + b_H2g/
pow(1.+x_H2g/b5,e2_H2g);
3941 fshield_H2s = 0.965/
pow(1.+x_H2s/b5,e1_H2s);
4105 enum {DEBUG_LOC=
false};
4109 fprintf(
ioQQQ,
" Solomon H2 dest rates: TH85 %.2e BD96 %.2e Big %.2e excit rates: TH85 %.2e Big %.2e\n",
4132 return &(*p->second);
4208 for(
int i=0;i<rate.
nreactants && ipthis == -1;i++)
4217 double ratevi = rate.
a * rate.
rk();
4261 int ipspin = 0, ipfreein = 0;
4269 int ipspout = 0, ipfreeout = 0;
4279 int newsp = ipspout-ipspin;
4284 int nbondsbroken = ipfreeout-ipfreein;
4285 if (nbondsbroken <= 0)
4288 double fracbroken = nbondsbroken/((double)ipfreeout);
4289 ASSERT( fracbroken <= 1.0 );
4300 double ratesp = ratevi*newsp;
4302 ratesp *= fracbroken;
4336 double ratevi = rate.
a * rate.
rk();
4341 ratev += ipthis*ratevi;
4363 double heating = 0.;
4364 map<double,string> heatMap;
4376 bool lgCanSkip =
false;
4404 realnum reaction_enthalpy = 0.;
4413 for(
long i=0; i < rate.
nproducts; ++i )
4423 double heat = reaction_enthalpy*rate_tot*(1e10/AVOGADRO);
4424 heatMap[heat] = rate.
label;
4434 for( map<double,string>::reverse_iterator it = heatMap.rbegin(); it != heatMap.rend(); ++it, ++index )
4436 fprintf(
ioQQQ,
"DEBUGGG heat %li\t%li\t%.6e\t%s\n", index,
nzone, it->first, it->second.c_str() );
4441 for( map<double,string>::iterator it = heatMap.begin(); it != heatMap.end(); ++it, ++index )
4443 if( it->first >= 0. )
4445 fprintf(
ioQQQ,
"DEBUGGG cool %li\t%li\t%.6e\t%s\n", index,
nzone, it->first, it->second.c_str() );
4468 double T2 = T_gas / 100.;
4469 double Tg2 = T_grain / 100.;
4470 double S = 1./(1. + 0.4*sqrt(Tg2 + T2) + 0.2*T2 + 0.08*T2*T2);
molecule * reactants[MAXREACTANTS]
double sink_rate_tot(const char chSpecies[]) const
t_mole_global mole_global
STATIC void register_reaction_vectors(count_ptr< mole_reaction > rate)
double H2_Solomon_dissoc_rate_used_H2g
static void sort(MoleculeList::iterator start, MoleculeList::iterator end)
vector< double > reaction_rks
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
STATIC double sticking_probability_H_func(double T_gas, double T_grain)
molecule::nNucsMap::iterator nNucs_i
void mole_create_react(void)
double H2_photodissoc_ELWERT_H2s
double GammaK(long int ipLoEnr, long int ipHiEnr, long int ipOpac, double yield1, t_phoHeat *photoHeat)
NORETURN void TotalInsanity(void)
bool lgH2_Chemistry_BigH2
double Average_collH2_excit
STATIC void mole_check_reverse_reactions(void)
double H2_H2g_to_H2s_rate_used
realnum UV_Cont_rel2_Draine_DB96_face
virtual const char * name()=0
bool lgLeiden_Keep_ipMH2s
double H2star_forms_grains
bool operator<(const count_ptr< T > &a, const count_ptr< T > &b)
double H2_H2g_to_H2s_rate_BHT90
double Cont_Dissoc_Rate_H2g
double H2_photodissoc_used_H2s
double H2_Solomon_dissoc_rate_TH85_H2s
double H2_photodissoc_BHT90
STATIC void mole_h2_grain_form(void)
double findrk(const char buf[]) const
map< string, count_ptr< mole_reaction > > reactab
double rate_grain_op_conserve
double H2_photodissoc_TH85
double H2_Solomon_dissoc_rate_BD96_H2g
sys_float sexp(sys_float x)
double H2_Solomon_dissoc_rate_ELWERT_H2s
double H2_H2g_to_H2s_rate_ELWERT
molezone * findspecieslocal(const char buf[])
double source_rate_tot(const char chSpecies[]) const
STATIC char * getcsvfield(char **s, char c)
double rate_h2_form_grains_set
void iso_photo(long ipISO, long nelem)
ChemNuclideList nuclide_list
double H2_Solomon_dissoc_rate_BD96_H2s
bool parse_species_label(const char label[], ChemNuclideList &atomsLeftToRight, vector< int > &numAtoms, string &embellishments)
double anu(size_t i) const
double Average_collH2_dissoc_g
double Average_collH2_deexcit
map< string, bool > offReactions
double H2_Solomon_dissoc_rate_TH85_H2g
STATIC void newreact(const char label[], const char fun[], double a, double b, double c)
static t_version & Inst()
double Solomon_dissoc_rate_g
t_iso_sp iso_sp[NISO][LIMELM]
double H2_Solomon_dissoc_rate_BHT90_H2g
double xIonDense[LIMELM][LIMELM+1]
double Average_collH2_dissoc_s
bool fp_equal(sys_float x, sys_float y, int n=3)
double H2_Solomon_dissoc_rate_used_H2s
STATIC long parse_reaction(count_ptr< mole_reaction > &rate, const char label[])
molecule * products[MAXPRODUCTS]
void mole_cmp_num_in_out_reactions(void)
STATIC void mole_h_reactions(void)
double chem_heat(void) const
double H2_photodissoc_ELWERT_H2g
STATIC bool lgReactionTrivial(count_ptr< mole_reaction > &rate)
vector< double > old_reaction_rks
molecule * rvector[MAXREACTANTS]
const bool ENABLE_QUANTUM_HEATING
void qheat(vector< double > &, vector< double > &, long *, size_t)
double photodissoc_BigH2_H2s
double H2_H2g_to_H2s_rate_BD96
map< string, count_ptr< mole_reaction > >::iterator mole_reaction_i
molecule * findspecies(const char buf[])
double Average_collH_deexcit
void create_isotopologues_one_position(unsigned position, ChemNuclideList &atoms, vector< int > &numAtoms, string atom_old, string atom_new, string embellishments, string &newLabel)
valarray< class molezone > species
virtual double rk() const =0
double sink_rate(const molecule *const sp, const mole_reaction &rate) const
double photodissoc_BigH2_H2g
void mole_rk_bigchange(void)
count_ptr< chem_nuclide > findnuclide(const char buf[])
vector< diatomics * > diatoms
double powi(double, long int)
STATIC bool lgReactBalance(const count_ptr< mole_reaction > &rate)
map< string, count_ptr< mole_reaction > > functab
double H2_Solomon_dissoc_rate_BHT90_H2s
diatomics h2("h2", 4100.,&hmi.H2_total, Yan_H2_CS)
realnum GetDopplerWidth(realnum massAMU)
bool exists(const molecule *m)
double column(const genericState &gs)
map< const count_ptr< chem_nuclide >, int, element_pointer_value_less > nNucsMap
realnum GetAveVelocity(realnum massAMU)
realnum UV_Cont_rel2_Habing_TH85_face
double H2_photodissoc_used_H2g
double esc_PRD_1side(double tau, double a)
double **** PhotoRate_Shell
double rate_grain_J1_to_J0
realnum gas_phase[LIMELM]
realnum UV_Cont_rel2_Draine_DB96_depth
vector< count_ptr< chem_nuclide > > ChemNuclideList
realnum AtomicWeight[LIMELM]
STATIC void mole_generate_isotopologue_reactions(string atom_old, string atom_new)
double frac_H2star_hminus()
double findrate(const char buf[]) const
molecule * pvector[MAXPRODUCTS]
map< string, count_ptr< mole_reaction > >::const_iterator mole_reaction_ci
STATIC void compare_udfa(const count_ptr< mole_reaction > &rate)
bool lgH2_grain_deexcitation
char chH2_small_model_type
bool isMonatomic(void) const
double Average_collH_excit
#define DEBUG_ENTRY(funcname)
realnum UV_Cont_rel2_Habing_spec_depth
double powpq(double x, int p, int q)
double Solomon_dissoc_rate_s
STATIC double sticking_probability_H_HM79(double T_gas, double T_grain)
bool lgDifferByExcitation(const molecule &mol1, const molecule &mol2)
STATIC void plot_sparsity(void)
double Cont_Dissoc_Rate_H2s
double GammaBn(long int ipLoEnr, long int ipHiEnr, long int ipOpac, double thresh, double *ainduc, double *rcool, t_phoHeat *photoHeat)
void mole_update_rks(void)
diatomics hd("hd", 4100.,&hmi.HD_total, Yan_H2_CS)
STATIC double mole_get_equilibrium_condition(const char buf[])
STATIC void read_data(const char file[], void(*parse)(char *s))
STATIC double mole_partition_function(const molecule *const sp)
realnum UV_Cont_rel2_Habing_TH85_depth
int fprintf(const Output &stream, const char *format,...)
double dissoc_rate(const char chSpecies[]) const
mole_reaction * mole_findrate_s(const char buf[])
sys_float SDIV(sys_float x)
double pow(double x, int i)
molecule * rvector_excit[MAXREACTANTS]
double H2_H2g_to_H2s_rate_TH85
STATIC void parse_udfa(char *s)
virtual mole_reaction * Create() const =0
double Average_collH_dissoc_g
double hmrate4(double a, double b, double c, double te)
double esca0k2(double taume)
void GammaPrt(long int ipLoEnr, long int ipHiEnr, long int ipOpac, FILE *ioFILE, double total, double threshold)
STATIC void canonicalize_reaction(count_ptr< mole_reaction > &rate)
STATIC string canonicalize_reaction_label(const char label[])
double Average_collH_dissoc_s
t_secondaries secondaries
double CharExcIonOf[NCX][LIMELM][LIMELM+1]
double H2_Solomon_dissoc_rate_ELWERT_H2g
double HMinus_induc_rec_rate
STATIC void parse_base(char *s)
double CharExcRecTo[NCX][LIMELM][LIMELM+1]
void create_isotopologues(ChemNuclideList &atoms, vector< int > &numAtoms, string atom_old, string atom_new, string embellishments, vector< string > &newLabels)
vector< diatomics * >::iterator diatom_iter
double rate_h2_form_grains_used_total
static map< formula_species, count_ptr< udfa_reaction > > udfatab
molecule * pvector_excit[MAXPRODUCTS]
double HMinus_induc_rec_cooling
double H2star_forms_hminus
1.8.5