Skip to content

Commit

Permalink
1/1000 is not an error, start work on He,H (#25)
Browse files Browse the repository at this point in the history
  • Loading branch information
dschwen committed Feb 4, 2016
1 parent f6d1df1 commit 869452e
Show file tree
Hide file tree
Showing 5 changed files with 31 additions and 38 deletions.
3 changes: 3 additions & 0 deletions apps/mytrim_bobmsq.C
Original file line number Diff line number Diff line change
Expand Up @@ -149,6 +149,9 @@ int main(int argc, char *argv[])
material->prepare(); // all materials added
sample->material.push_back(material); // add material to sample
break;

default:
return 1;
}

const int nstep = 1000;
Expand Down
27 changes: 1 addition & 26 deletions apps/mytrim_solid.C
Original file line number Diff line number Diff line change
Expand Up @@ -92,13 +92,6 @@ int main(int argc, char *argv[])
// create a FIFO for recoils
std::queue<IonBase*> recoils;

Real norm;
MassInverter *m = new MassInverter;
EnergyInverter *e = new EnergyInverter;

Real A1, A2, Etot, E1, E2;
int Z1, Z2;

snprintf(fname, 199, "%s.Erec", argv[1]);
FILE *erec = fopen(fname, "wt");

Expand All @@ -118,29 +111,11 @@ int main(int argc, char *argv[])
ff1->_md = 0;
ff1->id = simconf->id++;

// generate fission fragment data
A1 = m->x(dr250());
//A1 = 131;

A2 = 235.0 - A1;
Etot = e->x(dr250());
E1 = Etot * A2 / (A1 + A2);
//E1 = 100;

E2 = Etot - E1;
Z1 = round((A1 * 92.0) / 235.0);
//Z1 = 54;

Z2 = 92 - Z1;

/* ff1->_Z = Z1;
ff1->_m = A1;
ff1->_E = E1 * 1.0e6; */

ff1->_Z = 53;
ff1->_m = 127;
ff1->_E = 70.0 * 1.0e6;

Real norm;
do
{
for (int i = 0; i < 3; ++i) ff1->_dir(i) = dr250() - 0.5;
Expand Down
35 changes: 25 additions & 10 deletions material.C
Original file line number Diff line number Diff line change
Expand Up @@ -119,25 +119,28 @@ Real
MaterialBase::rstop(const IonBase * ion, int z2)
{
Real e, vrmin, yrmin, v, vr, yr, vmin, m1;
Real a, b, q, q1, l, l0, l1;
Real a, b, q, /*q1,*/ l, l0, l1;
Real zeta;
int z1 = ion->_Z;
Real fz1 = Real(z1), fz2 = Real(z2);
const int z1 = ion->_Z;
const Real fz1 = Real(z1);
const Real fz2 = Real(z2);
Real eee, sp, power;
Real se;

// scoeff
Real lfctr = _simconf->scoef[z1-1].lfctr;
Real mm1 = _simconf->scoef[z1-1].mm1;
Real vfermi = _simconf->scoef[z2-1].vfermi;
const Real lfctr = _simconf->scoef[z1-1].lfctr;
const Real mm1 = _simconf->scoef[z1-1].mm1;
const Real vfermi = _simconf->scoef[z2-1].vfermi;
//Real atrho = _simconf->scoef[z2-1].atrho;

if (ion->_m == 0.0)
m1 = mm1;
else
m1 = ion->_m;

e = 0.001 * ion->_E / m1;
// we store ion energy in eV but ee is needed in keV
const Real ee = 0.001 * ion->_E;
e = ee / m1;

if (z1 == 1)
{
Expand All @@ -147,6 +150,7 @@ MaterialBase::rstop(const IonBase * ion, int z2)
std::cerr << "proton stopping not yet implemented!\n";
exit(1);
#endif
// rpstop(z2, e, pcoef, se(i))
}
else if (z1 == 2)
{
Expand All @@ -157,8 +161,19 @@ MaterialBase::rstop(const IonBase * ion, int z2)
exit(1);
#endif
// Helium electronic stopping powers [RST0820]
// Real heo = 10.0;
// Real he = std::max(heo, )
Real he0 = 10.0;
Real he = std::max(he0, e);
b = std::log(e);
a = 0.2865 + 0.1266 * b - 0.001429 * b*b + 0.02402 * b*b*b - 0.1135 * std::pow(b, 4.0) + 0.001475 * std::pow(b, 5.0);
Real heh = 1.0 - std::exp(-std::min(30.0, a));
// add z1^3 effect to He/H stopping power ratio heh
a = (1.0 + (0.007 + 0.00005 * z2) * std::exp(-std::pow((7.6 * std::max(0.0, std::log(he))), 2.0) ));
heh *= a * a;

// call rstop(z2, he, pcoeff, sp)
se = sp * heh * 4;
if (e <= he0)
se *= std::sqrt(e / he0);
}
else
{
Expand All @@ -185,7 +200,7 @@ MaterialBase::rstop(const IonBase * ion, int z2)
l1 = 0.0;
else if (q < std::max(0.0, 0.9 - 0.025 * fz1))
{//210
q1 = 0.2;
// q1 = 0.2; in the original code, but never used
l1 = b * (q - 0.2) / std::abs(std::max(0.0, 0.9 - 0.025 * fz1) - 0.2000001);
}
else if (q < std::max(0.0, 1.0 - 0.025 * std::min(16.0, fz1)))
Expand Down
2 changes: 1 addition & 1 deletion sample.C
Original file line number Diff line number Diff line change
Expand Up @@ -20,4 +20,4 @@ Real
SampleBase::rangeMaterial(Point & /* pos */, Point & /* dir */)
{
return 100000.0;
};
}
2 changes: 1 addition & 1 deletion trim.C
Original file line number Diff line number Diff line change
Expand Up @@ -361,7 +361,7 @@ TrimBase::followRecoil()
_recoil->_md = 0;
*/
return true;
};
}

void
TrimBase::vacancyCreation()
Expand Down

0 comments on commit 869452e

Please sign in to comment.