/////////////////////////////////////////////////////////////////////////////
//                                                                         //
// Ideal mixture of two EOSs                                               // 
//                                                                         //
// Burkhard Militzer                                  Berkeley, 09-01-19   //
//                                                                         //
/////////////////////////////////////////////////////////////////////////////

#ifndef _IDEALMIXTURE_
#define _IDEALMIXTURE_

#include "Standard.h"
#include "EOSTableSimple.h"
#include "FindRoot.h"

class IdealMixture {
 public:
  EOSTableDiffT eos1;
  double N1;
  EOSTableDiffT eos2;
  double N2;
  //  double f1,NAtomsFU1,massFraction1,N1;
  //  double f2,NAtomsFU2,massFraction2,N2;

  bool hugoniotFlag;
  double hugoniotE0;
  double hugoniotP0;
  double hugoniotnn0;

  IdealMixture(EOSTableDiffT & eos1_, EOSTableDiffT & eos2_):
    eos1(eos1_), eos2(eos2_), hugoniotFlag(false) {
  }
  IdealMixture(EOSTableDiffT & eos1_, const double N1_, EOSTableDiffT & eos2_, const double N2_):
    eos1(eos1_), eos2(eos2_), hugoniotFlag(false) {
    SetNumberFraction(N1_,N2_);
  }

  void SetNumberFraction( const double N1_) {
    SetNumberFraction(N1_,1.0-N1_);
  }
  void SetNumberFraction(const double N1_, const double N2_) {
    N1 = N1_;
    N2 = N2_;
    SetHugoniotParameters(); // should be called in  constructor and when N1 and N2 become available
    /*
    //    Write3(N1,N2,N1+N2);
    double mFU1 = eos1.GetMassFormulaUnit();
    double mFU2 = eos2.GetMassFormulaUnit();
    double m1 = N1*mFU1;
    double m2 = N2*mFU2;
    double m  = m1+m2;
    double massFraction1 = m1/m;
    double massFraction2 = m2/m;
    //    Write2(m1/(PC::u/PC::me),m2/(PC::u/PC::me));
    //    Write2(massFraction1,massFraction2);
    //    Write2(massFraction1,m1/m);
    //    Write2(massFraction2,m2/m);
    */
  }
  double GetMassFormulaUnit() const {
    return N1*eos1.GetMassFormulaUnit() + N2*eos2.GetMassFormulaUnit();
  }

  string GetName() const {
    return eos1.GetName() + "+" + eos2.GetName();
  }
  double GetNElectronsFormulaUnit() const {
    return N1*eos1.GetNElectronsFormulaUnit() + N2*eos2.GetNElectronsFormulaUnit();
  }
  double AverageIonicCharge() const {
    return GetNElectronsFormulaUnit() / ( N1*eos1.nFU + N2*eos2.nFU );
  }

  /*
  double EntropyGivenTemperatureAndPressure(const double T, const double P) {
    double S1 = eos1.EntropyGivenTemperatureAndPressure(T,P) * f1; // because of f1, we now have the entropy per H2
    double S2 = eos2.EntropyGivenTemperatureAndPressure(T,P) * f2; // because of f2, we now have the entropy per H2O
    //    Write(-T*N1*log(N1));
    //    Write(-T*N2*log(N2));
    return S1*N1 + S2*N2 - T*N1*log(N1) - T*N2*log(N2);
  }
  */

  double Pressure(const double nn, const double T);
  double PressureErrorBar(const double nn, const double T) { // only approximate
    double P = Pressure(nn,T);
    double dP1 = eos1.PressureErrorBarGivenTemperatureAndPressure(T,P); // * f1; // because of f1, we now have the entropy per H2
    double dP2 = eos2.PressureErrorBarGivenTemperatureAndPressure(T,P); // * f2; // because of f2, we now have the entropy per H2O
    double dP = dP1*N1 + dP2*N2;
    return dP;
  }
  double Energy(const double nn, const double T) {
    double P = Pressure(nn,T);
    return EnergyGivenTemperatureAndPressure(T,P);
  }
  double EnergyGivenTemperatureAndPressure(const double T, const double P) {
    double E1 = eos1.EnergyGivenTemperatureAndPressure(T,P); // * f1; // because of f1, we now have the entropy per H2
    double E2 = eos2.EnergyGivenTemperatureAndPressure(T,P); // * f2; // because of f2, we now have the entropy per H2O
    double E = E1*N1 + E2*N2;
    return E;
  }
  double EnergyErrorBar(const double nn, const double T) {
    double P = Pressure(nn,T);
    double dE1 = eos1.EnergyErrorBarGivenTemperatureAndPressure(T,P); // * f1; // because of f1, we now have the entropy per H2
    double dE2 = eos2.EnergyErrorBarGivenTemperatureAndPressure(T,P); // * f2; // because of f2, we now have the entropy per H2O
    double dE = dE1*N1 + dE2*N2;
    return dE;
  }

  double DensityGivenTemperatureAndPressure(const double T, const double P) {
    //    cout << " z(1): "; Write2(T*PC::AUToK,P*PC::AUToGPa);
    double nn1 = eos1.DensityGivenTemperatureAndPressure(T,P); // number of FU per unit volume
    //    cout << " z(2): "; Write(nn1);
    double nn2 = eos2.DensityGivenTemperatureAndPressure(T,P);
    //    cout << " z(3): "; Write(nn2);
    double V1  = 1.0/nn1; // volume per FU1
    double V2  = 1.0/nn2; // volume per FU2
    double VMix = N1*V1 + N2*V2; // (number of FU) times (volume per FU)
    double nn   = 1.0/VMix; // number of FU per unit volume
    return nn;
  }
  double DensitySIGivenTemperatureAndPressure(const double T, const double P) {
    return MassDensitySI(DensityGivenTemperatureAndPressure(T,P));
  }
  double DensityGCCGivenTemperatureAndPressure(const double T, const double P) {
    return MassDensityGCC(DensityGivenTemperatureAndPressure(T,P));
  }

  double MassDensity(const double nn) const {
    double VMix = 1.0/nn; // nn = number of formula unit per unit volume, Vmix = volume of one mixture formular unit
    double mFU1 = eos1.GetMassFormulaUnit();
    double mFU2 = eos2.GetMassFormulaUnit();
    double mMix = N1*mFU1 + N2*mFU2;
    return mMix/VMix; // in atomic units = electron masses per Bohr^3
  }
  double MassDensitySI(const double nn) const {
    return MassDensity(nn)*PC::me/PC::ab/PC::ab/PC::ab;
  }
  double MassDensityGCC(const double nn) const {
    return MassDensitySI(nn)/1000.0;
  }
  double NumberDensityFromMassDensity(const double rho) const { // rho is in atomic units = electron masses per Bohr^3
    return rho/MassDensity(1.0);
  }
  double NumberDensityFromMassDensitySI(const double rho) const {
    return NumberDensityFromMassDensity(rho/PC::me*PC::ab*PC::ab*PC::ab);
  }
  double NumberDensityFromMassDensityGCC(const double rho) const {
    return NumberDensityFromMassDensitySI(1000.0*rho);
  }
  double MassDensityGivenTemperatureAndPressure(const double T, const double P) {
    double nn1 = eos1.DensityGivenTemperatureAndPressure(T,P); // FU per unit volume
    double nn2 = eos2.DensityGivenTemperatureAndPressure(T,P);
    double V1  = 1.0/nn1;
    double V2  = 1.0/nn2;
    double VMix = N1*V1 + N2*V2;    
    double mFU1 = eos1.GetMassFormulaUnit();
    double mFU2 = eos2.GetMassFormulaUnit();
    double mMix = N1*mFU1 + N2*mFU2;
    return mMix/VMix;
  }
  double MassDensitySIGivenTemperatureAndPressure(const double T, const double P) {
    return MassDensityGivenTemperatureAndPressure(T,P)*PC::me/PC::ab/PC::ab/PC::ab;
  }
  double MassDensityGCCGivenTemperatureAndPressure(const double T, const double P) {
    return MassDensitySIGivenTemperatureAndPressure(T,P)/1000.0;
  }
  
  double GetMinimumTemperature() const {
    return max(eos1.GetMinimumTemperature(),eos2.GetMinimumTemperature());
  }
  double GetMaximumTemperature() const {
    return min(eos1.GetMaximumTemperature(),eos2.GetMaximumTemperature());
  }

  bool MinimumMaximumPressureForTemperatureSafely(const double T, double & PMin, double & PMax, const bool check=true) {
    double PMin1 = eos1.MinimumPressureForTemperature(T);
    double PMin2 = eos2.MinimumPressureForTemperature(T);
    PMin = max(PMin1,PMin2);
    double PMax1 = eos1.MaximumPressureForTemperature(T);
    double PMax2 = eos2.MaximumPressureForTemperature(T);
    PMax = min(PMax1,PMax2);
    double ok = (PMin*(1.0+1e-8)<=PMax);
    if (check && !ok) {
      cout << endl << " Error: MinimumMaximumPressureForTemperatureSafely for T[K]=" << T*PC::AUToK << endl;
      cout << " " << eos1.name << " " ; Write2(PMin1*PC::AUToGPa,PMax1*PC::AUToGPa);
      cout << " " << eos2.name << " " ; Write2(PMin2*PC::AUToGPa,PMax2*PC::AUToGPa);
      error("No overlap in pressure ranges at T=",T*PC::AUToK,PMin*PC::AUToGPa,PMax*PC::AUToGPa);
      return false;
    }
    return ok;
  }
  double MinimumPressureForTemperatureSafely(const double T) {
    double PMin,PMax;
    MinimumMaximumPressureForTemperatureSafely(T,PMin,PMax);
    return PMin;
  }
  double MaximumPressureForTemperatureSafely(const double T) {
    double PMin,PMax;
    MinimumMaximumPressureForTemperatureSafely(T,PMin,PMax);
    return PMax;
  }
  double MinimumPressureForTemperature(const double T) {
    return MinimumPressureForTemperatureSafely(T);
    double PMin1 = eos1.MinimumPressureForTemperature(T);
    cout << " y(1): "; Write(PMin1*PC::AUToGPa);
    double PMin2 = eos2.MinimumPressureForTemperature(T);
    cout << " y(2): "; Write(PMin2*PC::AUToGPa);
    double PMin = max(PMin1,PMin2);
    cout << " y(3): "; Write(PMin*PC::AUToGPa);
    return PMin;
  }
  double  MinimumDensityForTemperature(const double T) {
    double PMin = MinimumPressureForTemperature(T);
    //    cout << " yy(1): "; Write2(T*PC::AUToK,PMin*PC::AUToGPa);
    double nnMin = DensityGivenTemperatureAndPressure(T,PMin);
    //    cout << " yy(2): "; Write3(T*PC::AUToK,PMin*PC::AUToGPa,MassDensityGCC(nnMin));
    return nnMin;
  }
  double MaximumPressureForTemperature(const double T) {
    return MaximumPressureForTemperatureSafely(T);
    double PMax1 = eos1.MaximumPressureForTemperature(T);
    //    Write(PMax1*PC::AUToGPa);
    double PMax2 = eos2.MaximumPressureForTemperature(T);
    //    Write(PMax2*PC::AUToGPa);
    double PMax = min(PMax1,PMax2);
    return PMax;
  }
  double MaximumDensityForTemperature(const double T) {
    double PMax = MaximumPressureForTemperature(T);
    //    cout << " xx(1): "; Write(PMax*PC::AUToGPa);
    double nnMax = DensityGivenTemperatureAndPressure(T,PMax);
    //    cout << " xx(2): "; Write(nnMax);
    return nnMax;
  }
  bool CheckForOverlapAtTemperature(const double T) {
    double PMin,PMax;
    return MinimumMaximumPressureForTemperatureSafely(T,PMin,PMax,false);
  }

  void PrintOnePoint(ostream & of, const double T, const double P) {
    double rho = DensityGCCGivenTemperatureAndPressure(T,P);
    double V = 1.0/DensityGivenTemperatureAndPressure(T,P);
    of << " rho[g/cc]= " << fix(11,8,rho);
    of << " V[A^3]= " << fix(16,8,V*PC::AUToA3);
    
    of << " T[K]= " << fix(9,0,T*PC::AUToK);
    of << " P[GPa]= " << fix(16,6,P*PC::AUToGPa);
    
    double E = EnergyGivenTemperatureAndPressure(T,P);
    of << " E[Ha]= " << fix(16,8,E);
    
    of << endl;
  }

  void SetHugoniotParameters() {
    hugoniotFlag = ( eos1.GetHugoniotFlag() & eos2.GetHugoniotFlag() );
    if (hugoniotFlag) {
      hugoniotE0 =  eos1.GetHugoniotE0()*N1 + eos2.GetHugoniotE0()*N2;
      hugoniotP0 = (eos1.GetHugoniotP0()*N1 + eos2.GetHugoniotP0()*N2)/(N1+N2); // this is not exact but we may not have E0 and V0 for a common P0 for both materials
      double nn1 = eos1.GetHugoniotnn0();
      double nn2 = eos2.GetHugoniotnn0();
      double V1  = 1.0/nn1; // volume per FU1
      double V2  = 1.0/nn2; // volume per FU2
      double VMix = N1*V1 + N2*V2; // (number of FU) times (volume per FU)
      hugoniotnn0 = 1.0/VMix; // number of FU per unit volume
    }
  }

  bool GetHugoniotFlag() const {
    //    return eos1.GetHugoniotFlag() & eos2.GetHugoniotFlag();
    return hugoniotFlag;
  }
  double GetHugoniotE0() const {
    //    return eos1.GetHugoniotE0()*N1 + eos2.GetHugoniotE0()*N2;
    return hugoniotE0;
  }
  double GetHugoniotP0() const  {// this is not exact but we may not have E0 and V0 for a common P0 for both materials
    //    return (eos1.GetHugoniotP0()*N1 + eos2.GetHugoniotP0()*N2)/(N1+N2); 
    return hugoniotP0;
  }
  double GetHugoniotnn0() const { // again if the two P0 values are too different this will cause problems even for heterogenous mixtures
    //    double nn1 = eos1.GetHugoniotnn0();
    //    double nn2 = eos2.GetHugoniotnn0();
    //    double V1  = 1.0/nn1; // volume per FU1
    //    double V2  = 1.0/nn2; // volume per FU2
    //    double VMix = N1*V1 + N2*V2; // (number of FU) times (volume per FU)
    //    double nn   = 1.0/VMix; // number of FU per unit volume
    //    return nn;
    return hugoniotnn0;
  }
  void OverwriteHugoniotnn0(const double nn0) {
    hugoniotnn0 = nn0;
  }
  void OverwriteHugoniotRho0GCC(const double rho0) {
    cout << "Setting a new Hugoniot rho0[g/cc]= " << rho0 << endl;
    OverwriteHugoniotnn0(NumberDensityFromMassDensityGCC(rho0));
  }
  void OverwriteHugoniotE0(const double E0) {
    cout << "Setting a new Hugoniot E0[Ha]= " << E0 << endl;
    hugoniotE0 = E0;
  }

};

class FindPressureGivenTemperatureAndDensity {
  IdealMixture & im;
  const double nn,T;

public:
  FindPressureGivenTemperatureAndDensity(IdealMixture & im_, const double nn_, const double T_):
    im(im_),nn(nn_),T(T_) {}

  double operator()(const double P) {
    double nnn = im.DensityGivenTemperatureAndPressure(T,P);
    return nnn - nn;
  }

};

double IdealMixture::Pressure(const double nn, const double T) {
  FindPressureGivenTemperatureAndDensity f(*this,nn,T);
  double PMin = MinimumPressureForTemperature(T);
  double PMax = MaximumPressureForTemperature(T);
  double PTol = 1e-8*(PMax-PMin);
  double P = FindRootRegulaFalsiSafely(PMin,PMax,PTol,f);
  return P;
}

#endif // _IDEALMIXTURE_