Commit 832f0c68 authored by Alain O' Miniussi's avatar Alain O' Miniussi
Browse files

Merge commit '2112f94e' into 597-compile-with-cuda

parents 23425d46 2112f94e
......@@ -759,6 +759,8 @@ namespace fargOCA {
real
oplin(real temp, real rho)
{
using std::pow;
using std::cbrt;
// data t234,t456,t678/1.6e3, 5.7e3, 2.28e6/
real t234 = 1.6e3;
real t456 = 5.7e3;
......@@ -780,20 +782,20 @@ namespace fargOCA {
/* test T against (T_23 * T_34 * T_34)**0.333333333 */
real opacity = std::numeric_limits<real>::quiet_NaN();
if( temp < t234*std::pow(rho,4.44444444e-2)){
if( temp < t234*pow(rho,4.44444444e-2)){
// disjoint opacity laws
real o1 = ak1*std::pow(temp, 2);
real o2 = ak2*temp/std::pow(temp, 8);
real o1 = ak1*pow(temp, 2);
real o2 = ak2*temp/pow(temp, 8);
real o3 = ak3*temp;
// parameters used for smoothing
real o1an = std::pow(o1, 2);
real o2an = std::pow(o2, 2);
real o1an = pow(o1, 2);
real o2an = pow(o2, 2);
// smoothed and continuous opacity law for regions 1, 2, and 3.
opacity=std::pow(std::pow(o1an*o2an/(o1an+o2an),2)+std::pow(o3/(1+1.e22/std::pow(temp,10)),4),0.25);
} else if ( temp < t456*std::pow(rho,2.381e-2)){
opacity=pow(pow(o1an*o2an/(o1an+o2an),2)+pow(o3/(1+1.e22/pow(temp,10)),4),0.25);
} else if ( temp < t456*pow(rho,2.381e-2)){
// to avoid overflow
real ts4 = 1.e-4*temp;
real rho13 = std::cbrt(rho);
real rho13 = cbrt(rho);
real rho23 = rho13*rho13;
real ts42 = ts4*ts4;
real ts44 = ts42*ts42;
......@@ -804,28 +806,28 @@ namespace fargOCA {
real o4=bk4*rho23/(ts48*ts4);
real o5=bk5*rho23*ts42*ts4;
// parameters used for smoothing
real o4an = std::pow(o4,4);
real o3an = std::pow(o3,4);
real o4an = pow(o4,4);
real o3an = pow(o3,4);
// smoothed and continuous opacity law for regions 3, 4, and 5.
// FortranForm opacty=((o4an*o3an/(o4an+o3an))+(o5/(1+6.561e-5/ts48))**4)**0.25
opacity = std::pow((o4an*o3an/(o4an+o3an))+std::pow(o5/(1.0+6.561e-5/ts48),4.0),0.25);
} else if (temp < t678*std::pow(rho,2.267e-1) || rho <= 1.0e-10){
opacity = pow((o4an*o3an/(o4an+o3an))+pow(o5/(1.0+6.561e-5/ts48),4.0),0.25);
} else if (temp < t678*pow(rho,2.267e-1) || rho <= 1.0e-10){
// to avoid overflow
real ts4 = 1.e-4*temp;
real rho13 = std::cbrt(rho);
real rho13 = cbrt(rho);
real rho23 = rho13*rho13;
// disjoint opacity laws for 5, 6, and 7.
real o5 = bk5*rho23*std::pow(ts4,3);
real o6 = bk6*rho13*std::pow(ts4,10);
real o7 = bk7*rho/(std::pow(ts4, 2)*sqrt(ts4));
real o5 = bk5*rho23*pow(ts4,3);
real o6 = bk6*rho13*pow(ts4,10);
real o7 = bk7*rho/(pow(ts4, 2)*sqrt(ts4));
// parameters used for smoothing
real o6an = std::pow(o6,2);
real o7an = std::pow(o7,2);
real o6an = pow(o6,2);
real o7an = pow(o7,2);
// smoothed and continuous opacity law for regions 5, 6, and 7.
/* FortranForm
opacty=((o6an*o7an/(o6an+o7an))**2+(o5/(1+(ts4/(1.1*rho**0.04762))**10)**4)**0.25 */
opacity = std::pow(std::pow(o6an*o7an/(o6an+o7an), 2)+std::pow(o5/(1+std::pow(ts4/(1.1*std::pow(rho,0.04762)),10)),4.0),0.25);
opacity = pow(pow(o6an*o7an/(o6an+o7an), 2)+pow(o5/(1+pow(ts4/(1.1*pow(rho,0.04762)),10)),4.0),0.25);
} else {
// no scattering! (So undefined behavior is ok ?)
opacity = std::numeric_limits<real>::infinity();
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment