gwm17/catima
1
0
Fork 0
mirror of https://github.com/gwm17/catima.git synced 2024-11-26 12:08:52 -05:00
Code Issues Projects Releases Wiki Activity

Merge pull request #40 from hrosiak/reactions2

Browse Source 
reactions2
This commit is contained in:
Andrej Prochazka 2018-07-31 18:54:36 +02:00 committed by GitHub
parent 852ce689fe 3ddc3205fb
commit 4e1dd18a49
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
2 changed files with 14 additions and 5 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
reactions.cpp
View File

@ -22,16 +22,17 @@ double nonreaction_rate(Projectile &projectile, const Material &target, const Co
int ap = lround(projectile.A); int ap = lround(projectile.A);
int zp = lround(projectile.Z); int zp = lround(projectile.Z);
int zt = target.get_element(0).Z;
int at = abundance::get_isotope_a(zt,0); // most abundand natural isotope mass
auto data = _storage.Get(projectile,target,c); auto data = _storage.Get(projectile,target,c);
Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num); Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num);
if(energy_out(projectile.T, target.thickness(), range_spline) < emin_reaction)return -1.0; if(energy_out(projectile.T, target.thickness(), range_spline) < emin_reaction)return -1.0;
auto sigma_r = [&](double th){ auto sigma_r = [&](double th){
double stn_sum=0.0, sum=0.0; double stn_sum=0.0, sum=0.0;
double e = energy_out(projectile.T, th, range_spline); double e = energy_out(projectile.T, th, range_spline);
for(unsigned int i = 0;i<target.ncomponents();i++){ for(unsigned int i = 0;i<target.ncomponents();i++){
int zt = target.get_element(i).Z;
int at = abundance::get_isotope_a(zt,0); // most abundand natural isotope mass
stn_sum += target.molar_fraction(i); stn_sum += target.molar_fraction(i);
sum += target.molar_fraction(i)*SigmaR_Kox(ap,zp,e,at,zt); sum += target.molar_fraction(i)*SigmaR_Kox(ap,zp,e,at,zt);
} }
@ -50,9 +51,8 @@ double nonreaction_rate(Projectile &projectile, const Material &target, const Co
cs = target.number_density_cm2()*(cs0 + cs1)/2.0; cs = target.number_density_cm2()*(cs0 + cs1)/2.0;
} }
else{ else{
cs = catima::integrator.integrate(sigma_r,0,target.number_density_cm2()); cs = Avogadro*catima::integrator.integrate(sigma_r,0,target.thickness())/target.M();
} }
return exp(-cs*0.0001); return exp(-cs*0.0001);
} }

11
tests/test_reaction.cpp
View File

@ -11,11 +11,12 @@ const lest::test specification[] =
{ {
CASE("reaction"){ CASE("reaction"){
catima::Projectile proj{12,6,6,870}; catima::Projectile proj{12,6,6,870};
catima::Projectile proj2{238,92,92,500};
auto c = catima::get_material(6); auto c = catima::get_material(6);
auto h = catima::get_material(1); auto h = catima::get_material(1);
catima::Material water({{0,8,2},{0,1,1}},1.0); catima::Material water({{0,8,2},{0,1,1}},1.0);
c.thickness(2.0); c.thickness(2.0);
double r; double r,r2;
r= catima::nonreaction_rate(proj, c); r= catima::nonreaction_rate(proj, c);
EXPECT(r == approx(0.92,0.02)); EXPECT(r == approx(0.92,0.02));
@ -40,6 +41,14 @@ const lest::test specification[] =
r= catima::nonreaction_rate(proj, c); r= catima::nonreaction_rate(proj, c);
EXPECT(r == -1.0); EXPECT(r == -1.0);
proj.T=870;
water.thickness(1);
r= catima::nonreaction_rate(proj2, water);
r2= catima::nonreaction_rate(proj, water);
EXPECT( (r > 0 && r<1.0) );
EXPECT( (r2 > 0 && r2<1.0) );
EXPECT( r2>r );
} }
}; };