Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::granularPressureModels::Chao::
granularPressureCoeffPrime
(
    const phaseModel& phase1,
    const phaseModel& phase2,
    const volScalarField& Theta1,
    const volScalarField& Theta2,
    const volScalarField& g0,
    const volScalarField& g0prime,
    const dimensionedScalar& e
) const
{
    volScalarField m0
    (
        constant::mathematical::pi/6.0
       *(phase1.rho()*pow3(phase1.d()) + phase2.rho()*pow3(phase2.d()))
    );
    tmp<volScalarField> pCoeff
    (
        2.0*phase2*constant::mathematical::pi/(3.0*m0)*phase1.rho()
       *phase2.rho()*(1.0 + e)*pow3((phase1.d() + phase2.d())/2.0)
       *(Theta1 + Theta2 + 0.2*magSqr(phase1.U() - phase2.U()))
    );

    if (&phase1 == &phase2)
    {
        return pCoeff*(2.0*g0 + phase1*g0prime);

    }
    else
    {
        return pCoeff*(g0 + phase1*g0prime);
    }
}

double dgrowth_factor_dt(double a, double omega_m, double omega_l)
{
  double o_m, o_l, do_m, do_l,  hubble_a, g0, g, extra_fac, dg, dDdt;
  o_m = omega_m * 1./pow3(a)/(omega_m * 1./pow3(a) + omega_l);
  o_l = omega_l / (omega_m * 1./pow3(a) + omega_l);
  do_m = -3 * omega_m * omega_l / (a * a * a * a) / (omega_m / (a * a * a) + omega_l) / (omega_m / (a * a * a) + omega_l);
  do_l = -do_m;

  hubble_a = sqrt(omega_m/(a*a*a) + omega_l);

  g0 = 2.5 * omega_m / (pow(omega_m, 4.0/7.0) - omega_l + (1.0 + 0.5*omega_m)*(1.0 + omega_l/70.0));
  g = 2.5 * o_m / (pow(o_m, 4.0/7.0) - o_l + (1.0 + 0.5*o_m)*(1.0 + o_l/70.0));

  //da_dtau = sqrt(1 + o_m * (1 / a - 1) + o_l * (a * a - 1));   //tau = H0*t
  extra_fac = - ( 4/7.* pow(o_m, -3./7) * do_m - do_l
		  +(do_m /2. *(1 + o_l / 70.) - (1 + o_m / 2.) * do_l / 70)
		    /(1 + o_l/70.)/(1 + o_l / 70))/(pow(o_m, 4.0/7.0) - o_l + (1.0 + 0.5*o_m)*(1.0 + o_l/70.0))/(pow(o_m, 4.0/7.0) - o_l + (1.0 + 0.5*o_m)*(1.0 + o_l/70.0));

  dg = 2.5 * do_m / (pow(o_m, 4.0/7.0) - o_l + (1.0 + 0.5*o_m)*(1.0 + o_l/70.0)) + 2.5 * o_m * extra_fac;

  dDdt = a * hubble_a * (dg * a +  g )/ g0;

  return dDdt;

}

Foam::NonSphereDragForce<CloudType>::NonSphereDragForce
(
    CloudType& owner,
    const fvMesh& mesh,
    const dictionary& dict
)
:
    ParticleForce<CloudType>(owner, mesh, dict, typeName, true),
    phi_(readScalar(this->coeffs().lookup("phi"))),
    a_(exp(2.3288 - 6.4581*phi_ + 2.4486*sqr(phi_))),
    b_(0.0964 + 0.5565*phi_),
    c_(exp(4.9050 - 13.8944*phi_ + 18.4222*sqr(phi_) - 10.2599*pow3(phi_))),
    d_(exp(1.4681 + 12.2584*phi_ - 20.7322*sqr(phi_) + 15.8855*pow3(phi_)))
{
    if (phi_ <= 0 || phi_ > 1)
    {
        FatalErrorIn
        (
            "NonSphereDrag<CloudType>::NonSphereDrag"
            "("
                "const dictionary&, "
                "CloudType&"
            ")"
        )   << "Ratio of surface of sphere having same volume as particle to "
            << "actual surface area of particle (phi) must be greater than 0 "
            << "and less than or equal to 1" << exit(FatalError);
    }
}

Foam::tmp<Foam::volScalarField> Foam::dragModels::Tenneti::CdRe() const
{
    volScalarField alpha1
    (
        max(pair_.dispersed(), pair_.continuous().residualAlpha())
    );

    volScalarField alpha2
    (
        max(scalar(1) - pair_.dispersed(), pair_.continuous().residualAlpha())
    );

    volScalarField F0
    (
        5.81*alpha1/pow3(alpha2) + 0.48*pow(alpha1, 1.0/3.0)/pow4(alpha2)
    );

    volScalarField F1
    (
        pow(alpha1, 3)*max(pair_.Re(), residualRe_)
       *(0.95 + 0.61*pow3(alpha1)/sqr(alpha2))
    );

    // Tenneti et al. correlation includes the mean pressure drag.
    // This was removed here by multiplying F by alpha2 for consistency with
    // the formulation used in OpenFOAM
    return
        SchillerNaumann_->CdRe()/(alpha2*max(pair_.Re(), residualRe_)) +
        24.0*sqr(alpha2)*(F0 + F1);
}

/*===========================================================================
* create a timeline for current cosmological model
*===========================================================================*/
void create_timeline2(double a_init, double a_fin, tlptr timeline)
{
   int    iloop;
   double a,t,omega,lambda,rhoc,hubble;
   
   for(iloop=0; iloop<MAXTIME2; iloop++)
     {
      a      = ((double)iloop+1.0)/(double)MAXTIME2 * (a_fin-a_init) + a_init;
      t      = calc_t(a);
      omega  = calc_omega(a);
      lambda = calc_lambda(a);
      hubble = H0    * sqrt(
                            simu.lambda0 * (1.-1./pow2(a))         +
                            simu.omega0  * (1./pow3(a)-1./pow2(a)) + 1./pow2(a));
      rhoc   = rhoc0 *     (
                            simu.lambda0*(1.-1./pow2(a))           +
                            simu.omega0*(1./pow3(a)-1./pow2(a))    + 1./pow2(a));
      
      timeline->a[iloop]      = a;
      timeline->t[iloop]      = t;
      timeline->omega[iloop]  = omega;
      timeline->lambda[iloop] = lambda;
      timeline->hubble[iloop] = hubble;
      timeline->age[iloop]    = t*Mpc/H0/Gyr;
      timeline->virial[iloop] = calc_virial(a);
     }
}

Foam::tmp<Foam::volScalarField> Foam::dragModels::Syamlal::K() const
{
    const phaseModel& phase1 = pair_.phase1();
    const phaseModel& phase2 = pair_.phase2();

    const volScalarField& Pfric = kineticTheorySystem_.frictionalPressure();
    scalar e = kineticTheorySystem_.es()[pair_];
    scalar Cf = kineticTheorySystem_.Cfs()[pair_];
    scalar pi = Foam::constant::mathematical::pi;

    tmp<volScalarField> g0 = kineticTheorySystem_.gs0(phase1, phase2);
    return
        (
            3.0*(1.0 + e)*(pi/2.0 + Cf*sqr(pi)/8.0)
           *phase1*phase1.rho()*phase2*phase2.rho()
           *sqr(phase1.d() + phase2.d())*g0*pair_.magUr()
        )
       /(
            2.0*pi
           *(
                phase1.rho()*pow3(phase1.d())
              + phase2.rho()*pow3(phase2.d())
            )
        )
      + C1_*Pfric;
}

standardRadiation::standardRadiation
(
    const surfaceFilmModel& owner,
    const dictionary& dict
)
:
    filmRadiationModel(typeName, owner, dict),
    QrPrimary_
    (
        IOobject
        (
            //pqr "Qr", // (net radiation) same name as Qr on primary region to enable mapping
            "Qin", // (incident radiation) same name as Qr on primary region to enable mapping
            owner.time().timeName(),
            owner.regionMesh(),
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        owner.regionMesh(),
        dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0),
        owner.mappedPushedFieldPatchTypes<scalar>()
    ),
    QrNet_
    (
        IOobject
        (
            "QrNet",
            owner.time().timeName(),
            owner.regionMesh(),
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        owner.regionMesh(),
        dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0),
        zeroGradientFvPatchScalarField::typeName
    ),
    /*QrFilm_ is used for diagnostics*/
    QrFilm_
    (
        IOobject
        (
            "QrFilm",
            owner.time().timeName(),
            owner.regionMesh(),
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        owner.regionMesh(),
        dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0),
        zeroGradientFvPatchScalarField::typeName
    ),
    delta_(owner.delta()),
    deltaMin_(readScalar(coeffs_.lookup("deltaMin"))),
    beta_(readScalar(coeffs_.lookup("beta"))),
    kappaBar_(readScalar(coeffs_.lookup("kappaBar")))
{}

/*==========================================================================
 * cubic:
 c
 c     Function for evaluating a cubic equation.  In the case where
 c     there is one real root, this routine returns it.  If there are three
 c     roots, it returns the smallest positive one.
 c
 c     The cubic equation is a1*x**3 + a2*x**2 + a3*x + a4 = 0.
 c
 c     Ref: Tallarida, R. J., "Pocket Book of Integrals and Mathematical
 c     Formulas," 2nd ed., Boca Raton: CRC Press 1992, pp. 8-9.
 c
 *==========================================================================*/
double cubic(double a1, double a2, double a3, double a4)
{
  double a,b,d,r,theta,pi,ar,ai,y,y1,y2,y3,p,q;
  
  if (fabs(a1) < ZERO)
    {
      fprintf(io.logfile,"ERROR: Quadratic/linear equation passed to CUBIC.");
      return(0.0);
    }
  
  
  p = a3/a1 - pow2(a2/a1)/3.;
  q = a4/a1 - a2*a3/pow2(a1)/3. + 2.*pow3(a2/a1)/27.;
  
  d = pow3(p)/27. + pow2(q)/4.;
  
  if (d > 0.)
    {
      a = -0.5*q + sqrt(d);
      b = -0.5*q - sqrt(d);
      
      if (a > 0.)
        a =   pow( a, 1./3.);
      else
        a = - pow(-a, 1./3.);
      
      if (b > 0.)
        b =   pow( b, 1./3.);
      else
        b = - pow(-b, 1./3.);
      
      y = a + b;
    }
  else
    {
      ar    = -0.5*q;
      ai    = sqrt(-d);
      r     = pow( pow2(ar)+pow2(ai), 1./6.);
      theta = atan2(ai,ar);
      y1    = 2. * r * cos(theta/3.) - a2/a1/3.;
      y     = y1;
      pi    = 4.*atan(1.);
      y2    = 2. * r * cos(theta/3.+2.*pi/3.) - a2/a1/3.;
      
      if (y < 0. || (y2 > 0. && y2 < y))
        y = y2;
      
      y3 = 2. * r * cos(theta/3.-2.*pi/3.) - a2/a1/3.;
      
      if (y < 0. || (y3 > 0. && y3 < y))
        y = y3;
    }
  
  return(y);
}

 static Complexd spharmGround(int l, int m, Complexd th, Complexd ph) {
     switch(l) {
     case 0:
         return 0.5*sqrt(1./M_PI);
     case 1:
         switch(m) {
         case  0:
             return  0.5*sqrt(3/ M_PI) * ph.r;
         case  1:
             return -0.5*sqrt(3/M_2PI) * ph.i * th;
         case -1:
             return  0.5*sqrt(3/M_2PI) * ph.i * th.conj();
         default:
             ;
         }
     case 2:
         switch(m) {
         case  0:
             return  0.25*sqrt( 5/ M_PI) * (3*ph.r*ph.r-1);
         case  1:
             return -0.50*sqrt(15/M_2PI) * ph.r*ph.i * th;
         case -1:
             return  0.50*sqrt(15/M_2PI) * ph.r*ph.i * th.conj();
         case  2:
             return  0.25*sqrt(15/M_2PI) * pow2(th*ph.i);
         case -2:
             return  0.25*sqrt(15/M_2PI) * pow2(th.conj()*ph.i);
         default:
             ;
         }
     case 3:
         switch(m) {
         case  0:
             return  0.250*sqrt(  7/ M_PI) * (5*pow3(ph.r) - 3*ph.r);
         case  1:
             return -0.125*sqrt( 21/ M_PI) * ph.i*(5*pow2(ph.r)-1) * th;
         case -1:
             return  0.125*sqrt( 21/ M_PI) * ph.i*(5*pow2(ph.r)-1) * th.conj();
         case  2:
             return  0.250*sqrt(105/M_2PI) * pow2(ph.i)*ph.r * pow2(th);
         case -2:
             return  0.250*sqrt(105/M_2PI) * pow2(ph.i)*ph.r * pow2(th.conj());
         case  3:
             return -0.125*sqrt( 35/ M_PI) * pow3(th*ph.i);
         case -3:
             return  0.125*sqrt( 35/ M_PI) * pow3(th.conj()*ph.i);
         default:
             ;
         }
     default:
         ;
     }
     return 0;
 }

//--------------------------------------------------------------------------------------------------
// getKernelValue
inline
double BasicSPH::getKernelValue(double dist) const
{
    double q = dist/_halfH;
    if (q<1.0)
    {
        return _kernelValueCoeff * ( pow3(2.0-q)-4*pow3(1.0-q) );
    }
    else
    {
        return _kernelValueCoeff * pow3(2.0-q);
    }
}

double c_correction(float mass, int snapnum)
{
  double c_original, c_new, Omega_new, Omega_original, ratio;

  c_original = 5 * pow(mass/10000., -0.1);
  Omega_new = Omega * 1./pow3(AA[snapnum]) / (Omega * 1./pow3(AA[snapnum]) + OmegaLambda);
  Omega_original = Omega_OriginalCosm * 1./pow3(AA_OriginalCosm[snapnum]) /
                  (Omega_OriginalCosm * 1./pow3(AA_OriginalCosm[snapnum]) + OmegaLambda_OriginalCosm);
  ratio = Omega_original/ Omega_new;
  c_new = find_c(c_original, ratio);

  return func_c(c_new) / func_c(c_original);
}

int main(void)
{
	for (int count_num = 100; count_num < 1000; count_num++)
	{
		int total = pow3(count_num % 10) + pow3(count_num / 100) + pow3(count_num / 10 - count_num / 100 * 10);
		if (total == count_num)
		{
			printf("%d\n", count_num);
		}
	}

	return 0;
}

Foam::tmp<Foam::volScalarField>
Foam::populationBalanceSubModels::aggregationKernels::sumAggregation::Ka
(
    const volScalarField& abscissa1,
    const volScalarField& abscissa2
) const
{   

    tmp<volScalarField> aggregationKernel(Ca_*(pow3(abscissa1) + pow3(abscissa2)));
        
    //Info << "aggregationKernel" << aggregationKernel.ref().dimensions() << endl;
    
    return aggregationKernel;
}

/*==========================================================================
* calc_virial: calculates the virial overdensity
*==========================================================================*/
double calc_virial(double a)
{
  double virial, age, omega_ta, eta, reduce, t1, t2;
  
  /* Check for a manual overdensity and return it if specified. Nothing else to do here then. */
  if(simu.UserDvir > 0)
    virial = (double)(simu.UserDvir);
  else
   {
#ifdef DARK_ENERGY
	/* Calculation of STH collapse is NOT implemented for Dark Energy models. Print a warning message, here or somewhere else */
	fprintf(stderr, "Warning: the calculation of the virial overdensity in dynamical dark energy cosmologies is NOT implemented in AHF.\n");
#else
    /* Figure out how old the universe is */
    age      = calc_t(a);
    
    /*  Figure out how much overdensity we need to reach maximum expansion by
     half the age of the universe.  The overdense expansion factor is defined
     to be 1 at maximum. */
    omega_ta = collapse(age/(double)2., (double)1.e-8);
    
    
    /* Figure out how far an object collapses to Virial equilibrium. */
    eta      = (double)2.*simu.lambda0/omega_ta; ///pow3(a); this obscure a^3 factor will prevent virial to approach the SCDM value of 178 at high redshift; not sure why it was there in the first place!
    
    if (eta < ZERO)
     {
      reduce = 0.5;
     }
    else
     {
      t1      = 2.*eta;
      t2      = -(2. + eta);
      reduce  = cubic(t1, (double)0., t2, (double)1.);
     }
    
    /* as found in M.A.K. Gross' original version */
    virial = pow3(a) * omega_ta/simu.omega0/pow3(reduce);

    /* check whether we want virial in terms of RhoCrit */
    if(simu.UseRhoBack == FALSE)
      virial *= calc_omega(a);
    
#endif
   }

  return (virial);
}

singleStepCombustion<CombThermoType, ThermoType>
::singleStepCombustion
(
    const word& modelType, const fvMesh& mesh
)
:
    CombThermoType(modelType, mesh),
    singleMixture_
    (
        dynamic_cast<singleStepReactingMixture<ThermoType>&>(this->thermo())
    ),
    wFuel_
    (
        IOobject
        (
            "wFuel",
            this->mesh().time().timeName(),
            this->mesh(),
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        this->mesh(),
        dimensionedScalar("zero", dimMass/pow3(dimLength)/dimTime, 0.0)
    )
{}

tmp<volScalarField> SpalartAllmarasIDDES<BasicTurbulenceModel>::fd
(
    const volScalarField& magGradU
) const
{
    return 1 - tanh(pow3(8*rd(this->nuEff(), magGradU)));
}

Foam::tmp<Foam::volScalarField>
Foam::radiation::greyMeanAbsorptionEmissionSoot::ECont(const label bandI) const
{
    tmp<volScalarField> E
    (
        new volScalarField
        (
            IOobject
            (
                "E",
                mesh_.time().timeName(),
                mesh_,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            mesh_,
            dimensionedScalar("E", dimMass/dimLength/pow3(dimTime), 0.0)
        )
    );

    if (mesh_.foundObject<volScalarField>("qrflamelet"))	    // PC added this line
    {
        const volScalarField& qrsoot =                              // PC added this line
            mesh_.lookupObject<volScalarField>("qrflamelet");       
        E().internalField() = qrsoot;       // PC modified this line
    }


    return E;
}

Foam::forceSuSp Foam::PlessisMasliyahDragForce<CloudType>::calcCoupled
(
    const typename CloudType::parcelType& p,
    const scalar dt,
    const scalar mass,
    const scalar Re,
    const scalar muc
) const
{
    scalar alphac(alphac_[p.cell()]);

    scalar cbrtAlphap(pow(1.0 - alphac, 1.0/3.0));

    scalar A =
        26.8*pow3(alphac)
       /(
            sqr(cbrtAlphap)
           *(1.0 - cbrtAlphap)
           *sqr(1.0 - sqr(cbrtAlphap))
          + SMALL
        );

    scalar B =
        sqr(alphac)
       /sqr(1.0 - sqr(cbrtAlphap));

    return forceSuSp
    (
        vector::zero,
        (mass/p.rho())
       *(A*(1.0 - alphac) + B*Re)*muc/(alphac*sqr(p.d()))
    );
}

tmp<volScalarField> standardRadiation::Shs()
{
    tmp<volScalarField> tShs
    (
        new volScalarField
        (
            IOobject
            (
                typeName + "_Shs",
                owner().time().timeName(),
                owner().regionMesh(),
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            owner().regionMesh(),
            dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0),
            zeroGradientFvPatchScalarField::typeName
        )
    );

    scalarField& Shs = tShs();
    const scalarField& QinP = QinPrimary_.internalField();
    const scalarField& delta = owner_.delta().internalField();
    const scalarField& alpha = owner_.alpha().internalField();

    // kvm, for now assume if surface is wet, then radiation is absorbed by film
    Shs = QinP*alpha;

    return tShs;
}

tmp<volScalarField> standardRadiation::Shs()
{
    tmp<volScalarField> tShs
    (
        new volScalarField
        (
            IOobject
            (
                typeName + "::Shs",
                owner().time().timeName(),
                owner().regionMesh(),
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            owner().regionMesh(),
            dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0),
            zeroGradientFvPatchScalarField::typeName
        )
    );

    scalarField& Shs = tShs();
    const scalarField& QrP = QrPrimary_.internalField();
    const scalarField& delta = delta_.internalField();

    Shs = beta_*(QrP*pos(delta - deltaMin_))*(1.0 - exp(-kappaBar_*delta));

    // Update net Qr on local region
    QrNet_.internalField() = QrP - Shs;
    QrNet_.correctBoundaryConditions();

    return tShs;
}