Commit cc138d7a authored by Alain O' Miniussi's avatar Alain O' Miniussi
Browse files

Bell Lin cleanup

parent 703d93cc
......@@ -758,11 +758,6 @@ namespace fargOCA {
real
oplin(real temp, real rho)
{
// data power1,power2,power3/4.44444444e-2, 2.381e-2, 2.267e-1/
real power1 = 4.44444444e-2;
real power2 = 2.381e-2;
real power3 = 2.267e-1;
// data t234,t456,t678/1.6e3, 5.7e3, 2.28e6/
real t234 = 1.6e3;
real t456 = 5.7e3;
......@@ -784,66 +779,55 @@ 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,power1)){
// different powers of temperature
real t2=temp*temp;
real t4=t2*t2;
real t8=t4*t4;
real t10=t8*t2;
if( temp < t234*std::pow(rho,4.44444444e-2)){
// disjoint opacity laws
real o1=ak1*t2;
real o2=ak2*temp/t8;
real o3=ak3*temp;
real o1 = ak1*std::pow(temp, 2);
real o2 = ak2*temp/std::pow(temp, 8);
real o3 = ak3*temp;
// parameters used for smoothing
real o1an=o1*o1;
real o2an=o2*o2;
real o1an = std::pow(o1, 2);
real o2an = std::pow(o2, 2);
// smoothed and continuous opacity law for regions 1, 2, and 3.
opacity=std::pow(std::pow(o1an*o2an/(o1an+o2an),2.0)+std::pow(o3/(1.0+1.e22/t10),4.0),0.25);
} else if ( temp < t456*std::pow(rho,power2)){
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)){
// to avoid overflow
real ts4=1.e-4*temp;
real rho13=std::pow(rho,0.333333333);
real rho23=rho13*rho13;
real ts42=ts4*ts4;
real ts44=ts42*ts42;
real ts48=ts44*ts44;
real ts4 = 1.e-4*temp;
real rho13 = std::pow(rho,0.333333333);
real rho23 = rho13*rho13;
real ts42 = ts4*ts4;
real ts44 = ts42*ts42;
real ts48 = ts44*ts44;
// disjoint opacity laws for 3, 4, and 5.
real o3=bk3*ts4;
real o4=bk4*rho23/(ts48*ts4);
real o5=bk5*rho23*ts42*ts4;
// parameters used for smoothing
real o4an=std::pow(o4,4.0);
real o3an=std::pow(o3,4.0);
real o4an = std::pow(o4,4);
real o3an = std::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,power3) || rho <= 1.0e-10){
} else if (temp < t678*std::pow(rho,2.267e-1) || rho <= 1.0e-10){
// to avoid overflow
real ts4=1.e-4*temp;
real rho13=std::pow(rho,0.333333333);
real rho23=rho13*rho13;
real ts42=ts4*ts4;
real ts44=ts42*ts42;
real ts4 = 1.e-4*temp;
real rho13 = std::pow(rho,0.333333333);
real rho23 = rho13*rho13;
// disjoint opacity laws for 5, 6, and 7.
real o5=bk5*rho23*ts42*ts4;
real o6=bk6*rho13*ts44*ts44*ts42;
real o7=bk7*rho/(ts42*sqrt(ts4));
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));
// parameters used for smoothing
real o6an=o6*o6;
real o7an=o7*o7;
real o6an = std::pow(o6,2);
real o7an = std::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.0)+std::pow(o5/(1.0+std::pow(ts4/(1.1*std::pow(rho,0.04762)),10.0)),4.0),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);
} else {
// no scattering! // -> actually, that's called an undefined result at best.
// and dividing with 0 is not such a good idea...
real ts42 = 0;
real ts4 = 0;
opacity=bk7*rho/(ts42*sqrt(ts4));
// no scattering! (So undefined behavior is ok ?)
opacity = std::numeric_limits<real>::infinity();
}
return opacity;
}
......
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