/////////////////////////////////////////////////////////////////////////////
//                                                                         //
// Ideal Fermi gas in 3D                                                   // 
//                                                                         //
// Burkhard Militzer                                Washington, 03-23-07   //
//                                                                         //
/////////////////////////////////////////////////////////////////////////////

#include "Standard.h"
#include "GKIntegration.h"
#include "FindRoot.h"
#include "FermiGas.h"
#include "Spline.h"

inline double FermiDirac(const double z, const double a) {
  return 1.0/(exp(z-a)+1.0);
}

class FermiIntegral {
public:
  FermiIntegral(const double aa, const double ee):a(aa),e(ee) {}

  const double a,e;
  double operator() (const double z) {
    return pow(z,e)*FermiDirac(z,a);
  }
};

class FermiIntegralInfinity {
public:
  FermiIntegralInfinity(const double aa, const double ee):a(aa),e(ee) {}

  const double a,e;
  double operator() (const double z1) {
    double z=1.0/z1;
    return z*z * pow(z,e)*FermiDirac(z,a);
  }
};

double CalculateFermiIntegral(const double a, const double ee) {
  FermiIntegral fi1(a,ee);
  GKIntegration<FermiIntegral,GK15> fii1(fi1,"Fermi integral I");
  double f1 = fii1.Integrate(0.0,1.0,1e-8);

  FermiIntegralInfinity fi2(a,ee);
  GKIntegration<FermiIntegralInfinity,GK15> fii2(fi2,"Fermi integral II");
  double f2 = fii2.Integrate(0.0,1.0,1e-8);
  
  return f1+f2;
}

double CC(const double lambda, const double beta) {
  return 4.0*pi*pi*pow(lambda*beta,3.0/2.0);
}

class FindMuFermiGas {
public:
  FindMuFermiGas(const double la, const double be, const double n, const double gg):
    lambda(la),beta(be),ne(n),g(gg),c(CC(lambda,beta)) {}

  const double lambda,beta,ne,g,c;

  double operator()(const double a) const {
    double diff = ne/g * c - CalculateFermiIntegral(a,1.0/2.0);
    //    Write2(a,diff);
    return diff;
  }
};

double MuFermiGas(const double lambda, const double beta, const double ne, const double g) {
  FindMuFermiGas fmfg(lambda,beta,ne,g);

  double betaMu = FindRootRegulaFalsi(5.0,10.0,1e-8,fmfg);
  double mu = betaMu/beta;
  return mu;
}

// U per particle
double EnergyFermiGas(const double lambda, const double beta, const double ne, const double g) {
  double mu = MuFermiGas(lambda,beta,ne,g);
  //  Write(mu);

  double fi12 = CalculateFermiIntegral(beta*mu,1.0/2.0); 
  //  Write2(mu,fi12);

  double fi32 = CalculateFermiIntegral(beta*mu,3.0/2.0);
  //  Write3(mu,fi12,fi32);
  
  double c = CC(lambda,beta);

  double ne2 = g/c * fi12;
  
  if (fabs(ne-ne2)>1e-8*ne)
    error("Problem with Fermi integrals");

  double uFermi = g*fi32 / (c*beta*ne);

  return uFermi;
}

// 3 PV = 2 U
double PressureFermiGas(const double lambda, const double beta, const double ne, const double g) {
  double P = 2.0/3.0 * EnergyFermiGas(lambda,beta,ne,g) * ne;
  return P;
}

// regular E_F = highest occupied level at T=0
double FermiEnergy(const double lambda, const double ne, const double g) {
  return lambda * pow(6.0*pi*pi*ne/g,2.0/3.0);
}

double FermiDegeneracy(const double lambda, const double ne, const double g, const double T) {
  return T/FermiEnergy(lambda,ne,g);
}

// F   = -PV        + mu*N
// 3PV =  2 U = 2 E_kin        (no pot. E)
// F   = -2/3 * U   + mu*N
// F/N = -2/3 * U/N + mu
// return F/N
double FreeEnergyFermiGas(const double lambda, const double beta, const double ne, const double g) {
  double mu     = MuFermiGas(lambda,beta,ne,g);
  double uFermi = EnergyFermiGas(lambda,beta,ne,g);
  double fFermi = -2.0/3.0*uFermi + mu;
  return fFermi;
}

// F/N = -2/3 * U/N + mu
// S/N = (U/N-F/N) / T
// S/N = (U/N+2/3*U/N- mu) /T
// S/N = (5/3 * U/N  - mu) /T
double EntropyFermiGas(const double lambda, const double beta, const double ne, const double g) {
  double mu     = MuFermiGas(lambda,beta,ne,g);
  double uFermi = EnergyFermiGas(lambda,beta,ne,g);
  double sFermi = (5.0/3.0*uFermi - mu)*beta;
  return sFermi;
}

// move inside SahaFermi.C orherwise it will get initialized all the time!!!
// FermiFunctionTable fermiFunctionTable;


/////////////////////////////////////////////////////////////////////////////////////

// PV = N kb T     = N / beta
// P  = N/V / beta = ne/beta
double PressureBoltzmannGas(const double beta, const double ne) {
  return ne/beta;
}

// per particle
double EnergyBoltzmannGas(const double beta) {
  return 3.0/2.0/beta;
}

// Wrong name
// double LambdaDebyeLambdaBeta(const double lambda, const double beta) {

// corrected name
// input: lambda = hb^2/(2*m) in A.U.
//        beta   = 1/(kb T)   in A.U.
double LambdaDeBroglieLambdaBeta(const double lambda, const double beta) {
  double laDB = sqrt(4.0*pi*lambda*beta);
  return laDB;
}


// input: mass in m_e and T in [Ha/kb]
double LambdaDeBroglieMassTemperature(const double m, const double T) {
  //  return PC::h/sqrt(2.0*pi*PC::kb*T*m) / PC::ab;
  //  return 2.0*pi/sqrt(2.0*pi*T*m); // in atmomic units
  return LambdaDeBroglieLambdaBeta(0.5/m,1.0/T);
}

// input: nn = N/V with V=ab^3 
// Now    returns  f = F/N 
// Before returned f = F/V
//
// REMEMBER DEFAULT ARGUMENTS: g=1 and e=0.0
double FreeEnergyBoltzmannGas(const double lambda, const double beta, const double nn, const double g, const double e) {
  if (nn<0.0) error("nn",nn);
  if (nn==0.0) return 0.0;
  if (beta<0.0) error("beta",beta);
  double laDB = LambdaDeBroglieLambdaBeta(lambda,beta);
  double f = ( log ( laDB*laDB*laDB /g*nn ) - 1.0 ) / beta + e; // F/N
  //  f *= nn;  // F/N -> F/V
  return f;
}

// return also per particle: S/N
double EntropyBoltzmannGas(const double lambda, const double beta, const double nn, const double g, const double e) {
  double F = FreeEnergyBoltzmannGas(lambda, beta, nn, g, e);
  double E = EnergyBoltzmannGas(beta);
  double S = (E-F)*beta;    // F = E -TS
  return S;
}

/////////////////////////////////////////////////////////////////////////////////////

/*
int main() {
  cout.precision(8);
  const double AUToK = 315774.660873;     // old and wrong 315774.9076860828;
  const double AUToA = 0.529177208607388; // old and wrong 0.52917679499658208;
  const double AUToM = AUToA*1e-10;
  const double UToAU = 1.66053873e-27 / 9.10938188e-31; // nuclear mass unit 'u' / mass pf electron

  // electrons
  double m      = 1.0;
  double g      = 2.0;

  // protons
  //  m      = 1836.152786976163; // p/me
  //  g      = 2.0; // spin state should not matter, only add constant to free energy
  //  g = 1.0; // compare with Hugh

  // helium atoms
  //  m      = 4.00260*UToAU;
  //  g      = 2.0;

  // neon atoms
  //  m      = 20.1797*UToAU;
  //  g      = 1.0;
  //  g      = 2.0;

  // argon atoms
  //  m      = 39.948*UToAU;
  //  g      = 2.0;

  double lambda = 0.5/m;
  //  double TinK   =  125000.0;
  //  double TinK   = 1250000.000;
  //  double TinK   = 10000.000;
  //  double TinK   = 5000.000;
  //  double T      = TinK / AUToK;
  double T = 1.0/0.0010375984; // Kevin value at 304332275K
  double TinK = T*AUToK;

  double beta   = 1.0 / T;

  Write3(m,TinK,LambdaDeBroglieMassTemperature(m,T));

  //  int N = 110; // this number should not matter below if rs is kept constant
  //  int N = 64; // this number should not matter below if rs is kept constant
  //  int N = 1; // this number should not matter below if rs is kept constant
  int N = 8*6;

  double rs = 0.5441299520956544;
  //  double rs = 0.9;
  //  double rs = 1.1;
  //  double rs = 1.86;

  double L = rs*pow(4.0/3.0*pi*N , 1.0/3.0);

  // Comment section out to start with density in cm^-3
  //  double nPercc = 9.01631e+23; // particles per cm^3
  //  //  double nPercc = 1.82389e+24; // particles per cm^3
  //  double nPerM3 = nPercc*1e6;
  //  double nAU    = nPerM3*AUToM*AUToM*AUToM;
  //  rs = pow(3.0/(4.0*pi*nAU),1.0/3.0);

  //  L = 9.74871962;
  //  L = 11.74275818937414587033; // Hugh H runs
  //  L = 11.05111842476021967327; // Hugh He runs
  //  L = 11.63426901208; // Hugh's H128 run at 220 GPa
  //  L = 11.7630916425 ; // Hugh's H128+Ne1 run at 220 GPa
  //  L = 11.79820275   ;   // Hugh's H128+Ar1 run at 220 GPa
  //  L = 10.9358073361 ; // Hugh's He64 run at 220 GPa
  //  L = 11.0509105549 ; // Hugh's He64+Ne1 run at 220 GPa
  //  L = 11.1311861210 ; // Hugh's He64+Ar1 run at 220 GPa
  //  L = 8.9363207;
  L = 6.365627; // Kevin's Carbon run at 1Mbar
  rs = L / pow(4.0/3.0*pi*N , 1.0/3.0);

  double V = L*L*L;
  double nn= double(N)/V;
  double V1= V/double(N);
  Write6(rs,TinK,V,nn,V1,m);

  double EF = FermiEnergy(lambda,nn,g);
  cout << endl;
  cout << "Check Fermi degeneracy (electron gas only, no ions, all in AU per electron]:" << endl;
  Write3(T,EF,T/EF);
  double E0 = 3.0/5.0*EF; // T=0 limit of energy of Fermi gas
  cout << " E(T->0) = " << E0 << endl << endl;

  double mu = MuFermiGas(lambda,beta,nn,g); // mu < EF, mu(T->0)=EF
  if (TinK==125000.0) Write2("0.224743",mu);
  cout << " Check mu<EF: ";
  Write2(mu,mu/EF);
  cout << endl;

  double P      = PressureFermiGas(lambda,beta,nn,g);
  double PBoltz = PressureBoltzmannGas(beta,nn);
  cout << " Check P>PBoltz: "; 
  Write4(T/EF,P,PBoltz,P/PBoltz);
  cout << endl;

  double E = EnergyFermiGas(lambda,beta,nn,g);
  double EBoltz = EnergyBoltzmannGas(beta);
  if (TinK==125000.0) Write2("0.714316",E);
  cout << " Compare energies: "; 
  Write4(E,EBoltz,E/EBoltz,E0);

  cout << endl;

  double F = FreeEnergyFermiGas(lambda,beta,nn,g);
  double TS = E - F;
  double S  = TS/T;
  //  Write2(F,-2.0/3.0*E+mu); 
  Write4(F,E,TS,S); // F = E-TS < E   ,   F(T->0) = E
  cout << " Fermi entropy in AU = 1/kb : "; Write(S);
  cout << " Fermi entropy in Ha/kb     : "; Write(S/AUToK);
  cout << endl;

  double FBoltz = FreeEnergyBoltzmannGas(lambda,beta,nn,g);
  double TSBoltz = EBoltz - FBoltz;
  double SBoltz  = TSBoltz/T;
  double PVBoltz = PBoltz*V1;
  double GBoltz  = FBoltz+PVBoltz;

  Write4(FBoltz,EBoltz,TSBoltz,SBoltz);
  Write3(GBoltz,FBoltz,PVBoltz);
  Write4(FBoltz*N,EBoltz*N,TSBoltz*N,SBoltz*N);
  Quit("F0");
  Write3(GBoltz*N,FBoltz*N,PVBoltz*N);
  cout << " Note that SBoltz can go negative for small T." << endl;
  cout << " Boltzmann entropy in AU = 1/kb : "; Write(SBoltz);
  cout << " Boltzmann entropy in Ha/kb     : "; Write(SBoltz/AUToK);
  cout << endl;

  cout << " Check E = d(F*beta)/d beta" << endl;
  double dBeta = 1e-6*beta;
  double beta1 = beta+dBeta;
  double F1 = FreeEnergyFermiGas(lambda,beta1,nn,g);
  double EE = (F1*beta1-F*beta)/dBeta;
  Write4(T/EF,EF,E,EE);

  double FB1 = FreeEnergyBoltzmannGas(lambda,beta1,nn,g);
  double EEB = (FB1*beta1-FBoltz*beta)/dBeta;
  Write2(EBoltz,EEB);
}
*/