/////////////////////////////////////////////////////////////////////////////
//                                                                         //
//  CMS Method                                                             //
//                                                                         //
//  B. Militzer                                     Berkeley 04-18-17      //
//                                                                         //
//  based on William Hubbard's methodology and                             //
//  based on Sean Wahl's tidal CMS code                                    //
//                                                                         //
/////////////////////////////////////////////////////////////////////////////

#ifndef _CMS_
#define _CMS_

#include "Array.h"
#include "Form.h"
#include "Spline.h"
#include "LSpline.h"
#include "MCParameters.h"
#include "MatrixAlgebra.h"
#include "Point.h"

class Units;

class CMS {
 public:
  //  static const int nq = 48;          // decent value for Jupiter unless TWO are needed
  //  static const int nq = 48*4;        // needed for improved Jupiter's TWE calculations, see **5891**
  //  static const int nq = 48*8;            // tried on 1/13/22
  static const int nq = 24;              // our standard value but is probably too many for U+N
  //  static const int nq = 12;          // probably good enough for U+N

  static const int np = 16;

  //  int nl, nlCore, nlEnv, nInt, nIntCore;
  int nl, nInt;

  //  normalized gravitational harmonics
  Array2 <double> jTilde,jTildeP;
  Array1 <double> jPP;

  // unnormalized gravitational harmonics
  Array2 <double> jReg,jRegP;

  // total (observed) gravitational moments
  Array1 <double> jTotal,jTotalPrev;

  // gaussian quadrature grid (mu,phi) and weights
  Array1 <double> xq, wq; // xqtheta 
  //  Array1 <double> xq2, wq2; // xq2phi
  
  // half of the nu = cos(phi) grid: phi = 0,pi
  //       real(dp), allocatable, dimension(:) :: xq2half, wq2half,&
  //         xq2phihalf 

  Array1 <double> lambda;
  Array1 <double> lambdaOrg;
  double scaleCoreRadius;
  LSpline sRhoLambda2; // spline representing lambda^2(rho) named 's'  earlier
  LSpline sLambda2Rho; // spline representing rho(lambda^2) named 'ss' earlier

  double PTarget; // 1.0 bar for Jupiter. 0.1 or 1.0 bar for Saturn
  Array1 <double> delta, pressure, potential; // 'press' is in planetary units (not CGS)
  Array1 <double> rho, rhoPrev, rhoZ1Derivative, rhoZ2Derivative, PAverage;
  Array2 <double> xi;
  Array2 <double> zeta, zetaPrev, zetaPrevPrev;
  Array2 <double> equiPot,equiPotPrev;
  Array1 <double> potOnEquator;
  Array1 <double> volume;      // volume enclosed by the surface of layer 'j'. See Moments.C
  Array1 <double> j2Hat;       // M*J2Layer, See Moments.C

  // precalculated values 
  // Legendre and associated Legendre polynomials
  Array2 <double> Pn_mu; // index order changed from (i,n) to (n,i). n = polynomial order, i = quadrature point index
  Array1 <double> Pn_0;

  double GM_SI,  a_SI,  M_SI;  // CGS->SI, multiply values 10^-3, 10^-6, 10^-2, 10^-3
  /*
  double dASigma;
  double  J2Juno,  J4Juno,  J6Juno,  J8Juno,  J10Juno;
  double dJ2Juno, dJ4Juno, dJ6Juno, dJ8Juno, dJ10Juno; 
  double  J3Juno,  J5Juno,  J7Juno,  J9Juno; // J11's error bar similar to J11 --> do not yet use
  double dJ3Juno, dJ5Juno, dJ7Juno, dJ9Juno;

  double  J2Target,  J4Target,  J6Target,  J8Target,  J10Target;
  double dJ2Target, dJ4Target, dJ6Target, dJ8Target, dJ10Target; 
  double  J3Target,  J5Target,  J7Target,  J9Target; // J11's error bar similar to J11 --> do not yet use
  double dJ3Target, dJ5Target, dJ7Target, dJ9Target;
  */
  double  J2Target,  J4Target;
  double dJ2Target, dJ4Target;

  double J2Tol, qRot, rotationPeriod, omega_SI;
  double MEarth_SI;

  double newtonStepSize;

  //  bool hasCore;
  //  double rCore; // mCore, rhoCore; // we no longer keep mCore up to date during the iterations
  //  double mCoreTarget;           
  //  double Z1MinFactor;           // increases minimum value of Z1 to Z1MinFactor*p.GetZ0ProtoSolar() - store inside CMS so that it is same everywhere

  // true if we use an EOS table, false if you a polytrope
  //  bool flagWorkWithEOSTable; 
  // flag to check whether EOS file has been read in
  //  bool flagEOSTableAlreadyReadInOrGenerated;

  //  LSpline PRhoSpline;
  //  double PMinAU; // in AU, marks surface pressure. Do not try to calculated rho(P) for a lower pressure
  //  double PMaxAU; // in AU, marks pressure above which we need to extrpolate
  bool densitiesAlreadyInitialized;
  //  Array1 <double> Log10PTable;
  //  Array1 <double> Log10rhoTableOriginal; // deprecated: used in CMS::SetUpEOSTable(const string & filename, const MCParameters & p) {
  //  Array1 <double> Log10rhoTable;         // new:        used in CMS::SetUpEOSTable(const MCParameters & p) {
  //  Array1 <double> Log10rhoTableModified;

  bool uranusFlag;
  
  int iWarning,nWarning;

  LSpline PCumMassSpline;          // log(P_AU),cum_mass
  LSpline rCumMassSpline;          // rE,       cum_mass
  LSpline PVolumetricRadiusSpline; // log(P_AU),rS
  LSpline PEquatorialRadiusSpline; // log(P_AU),rE
  LSpline EquatorialRadiusPSpline; // rE,log(P_AU)
 //  bool changeZ1Flag;

  Array1 <double> zetaPolar;
  Array1 <double> xiPolar;
  Array1 <Spline> xiMuSpline; // for every layer 'l', one spline xi as function of mu from equator to either pole

  Spline rhoSpline;
  Spline potentialSpline;
  Spline pressureSpline;

  //  bool filterSet, northFilter;
  //  double lFilter,sigmaFilter;

  int iterMin;

  int mLayers;
  double eta;
  Array1 <double> lambdaBoundaries;
  Array1 <double> lambdaBoundariesOrg;
  Array1 <int>    indexBoundaries;
  Array1 <double> dlambdaBoundaries_dEta;
  Array1 <double> dLa_dEta;

  double oxygenMass;
  double carbonNitrogenMass;
  double fO;  // oxygenMass         /(oxygenMass+carbonNitrogenMass);
  double fCN; // carbonNitrogenMass /(oxygenMass+carbonNitrogenMass);
  double fOExpected;  //  7.0*barotrope.mO                      / (7.0*barotrope.mO + 4.0*barotrope.mC + 1.0*barotrope.mN);
  double fCNExpected; //  (4.0*barotrope.mC + 1.0*barotrope.mN) / (7.0*barotrope.mO + 4.0*barotrope.mC + 1.0*barotrope.mN);
  bool   mCNFlag;
  bool   H2H3Flag;

  bool   hydrogenMassFlag;
  double hydrogenMassAbsorbed;
  double hydrogenMassReleased;
  bool   diffRhoFlag;
  double diffRhoCNHRhoOH; // diff between top of CNH layer (rho2) and bottom of OH later (rho1) 

  CMS(bool uranusFlag_): 
    uranusFlag(uranusFlag_),
    iWarning(0),
    nWarning(1000),
      //    changeZ1Flag(false),
      //    filterSet(false),
    iterMin(2), // replaces iter>1 in Iterate()
    eta(1.0),
    mCNFlag(false),
    H2H3Flag(false),
    hydrogenMassFlag(false),
    diffRhoFlag(false)
  {
  } 

  void SetUpHHeImmiscibilityPenalty(const string & HHeFileName);
  double HHeImmiscibilityPenaltyOnePoint(const double PP, const double TT) const;
  double HHeImmiscibilityPenalty(const double PSwitch1,const double PSwitch2) const;

  void Init(const int nl_, const int np_, const int nInt_) {

    nl = nl_;
    if (nl<256) warning("Number of layers is low for debugging",nl);

    //    np = np_;
    if (np!=np_) error("np != np_",np_);

    nInt = nInt_;
    if ( (nl-1)%nInt != 0 ) error("nl-1 must a multiple of nInt",nl,nInt);

    cout << " CMS: Init() called with "; Write4(nl,nInt,np,nq);

    jTilde.Resize(nl,np+1);
    jTildeP.Resize(nl,np+1);
    jPP.Resize(nl);
    jReg.Resize(nl,np+1);
    jRegP.Resize(nl,np+1);
    jTotal.Resize(np+1);
    jTotalPrev.Resize(np+1);

    lambda.Resize(nl);    
    lambda  = 0.0;
    delta.Resize(nl);     
    rho.Resize(nl);       
    rhoPrev.Resize(nl);   
    rhoZ1Derivative.Resize(nl); // added on 4/10/20 to enable faster M and J2 matching
    rhoZ2Derivative.Resize(nl); // added on 3/8/20 to allow Z2 scaling of core less models
    pressure.Resize(nl);     
    PAverage.Resize(nl);     
    potential.Resize(nl);     

    xi.Resize(nl,nq);       
    zeta.Resize(nl,nq);      
    zetaPrev.Resize(nl,nq); 
    zetaPrevPrev.Resize(nl,nq); 
    equiPot.Resize(nl,nq); 
    equiPotPrev.Resize(nl,nq); 
    potOnEquator.Resize(nl);

    volume.Resize(nl+1);
    j2Hat.Resize(nl+1);

    // precalculated values 
    // legendre and associated legendre polynomials
    Pn_mu.Resize(np+1,nq); // index order change from (i,n) to (n,i) 
    Pn_0.Resize(np+1);     
    Pn_mu=0.0;
    Pn_0 =0.0;
    Precompute();
    Reset();
  }
  void Reset(const bool print=true) {
    //    cout << "Reseting all variables" << endl;
    jTilde  = 0.0;
    jTildeP = 0.0;
    jReg  = 0.0;
    jRegP = 0.0;
    jTotal     = 0.0;
    jTotalPrev = 0.0;
   
    delta   = 0.0;
    rho     = 0.0;
    rhoPrev = 0.0;
    pressure  = 0.0;
    potential = 0.0;

    zeta     = 0.0;
    zetaPrev = 0.0;
    xi       = 0.0;

    MEarth_SI  = 5.9722e+24;      // kg

    if (uranusFlag) {
      // https://en.wikipedia.org/wiki/Uranus
      //      SetGM_SI(8.6810e25 * PC::G); // Uranus
      //      SetReferenceRadiusInKMAndPressureInBar(25559.0,1.0); 
      SetM_SI(14.536 * MEarth_SI);
      SetReferenceRadiusInKMAndPressureInBar(25559.0,1.0); 
      J2Target  = 0.351099e-2;  // J2*10^6 = 3510.99 +- 0.72
      J4Target  = -0.3361e-4;   // J4*10^6 = -33.61  +- 1
      dJ2Target = 0.000072e-2;
      dJ4Target = 0.0100e-4;
      SetRotationPeriod(17,14,40);
      //      SetRotationPeriod(16,34,24);
    } else {
      // https://en.wikipedia.org/wiki/Neptune
      //      SetGM_SI(1.02413e26 * PC::G); // Neptune
      //      SetReferenceRadiusInKMAndPressureInBar(24764,1.0); 
      SetM_SI(17.148 * MEarth_SI);
      SetReferenceRadiusInKMAndPressureInBar(24766.0,1.0); 
      J2Target  = 0.35294e-2;
      J4Target  =-0.358e-4;
      dJ2Target = 0.0045e-2;
      dJ4Target = 0.029e-4;
      SetRotationPeriod(16,6,40);
      //      SetRotationPeriod(17,27,29);
    }
    
    /*
    double ratio = aStandardKm / aEmployedKm;    
    bool printFlag = (abs(ratio-1.0)>1e-6);
    if (printFlag) cout << "Scaling J_n by ratio (aOrg/aNew)^n with aOrg[km]= " << aStandardKm << " aNew[km]= " << aEmployedKm << " (aOrg/aNew)=" << ratio << endl;
    if (printFlag) PrintJnTarget();
    J2Juno   *= pow(ratio,2.0);
    J4Juno   *= pow(ratio,4.0);
    J6Juno   *= pow(ratio,6.0);
    J8Juno   *= pow(ratio,8.0);
    J10Juno  *= pow(ratio,10.0);
    dJ2Juno  *= pow(ratio,2.0);
    dJ4Juno  *= pow(ratio,4.0);
    dJ6Juno  *= pow(ratio,6.0);
    dJ8Juno  *= pow(ratio,8.0);
    dJ10Juno *= pow(ratio,10.0);
    J3Juno   *= pow(ratio,3.0);
    J5Juno   *= pow(ratio,5.0);
    J7Juno   *= pow(ratio,7.0);
    J9Juno   *= pow(ratio,9.0);
    dJ3Juno  *= pow(ratio,3.0);
    dJ5Juno  *= pow(ratio,5.0);
    dJ7Juno  *= pow(ratio,7.0);
    dJ9Juno  *= pow(ratio,9.0);
    */
    //    if (printFlag) PrintJnTarget();
    PrintJnTarget();

    /*
    // by default we try to match the Juno measurements with our CMS models
    // but later we may want to our CMS models to match JnTarget=JnJuno-JnWind
    J2Target  = J2Juno  ;
    J4Target  = J4Juno  ; 
    J6Target  = J6Juno  ;
    J8Target  = J8Juno  ;
    J10Target = J10Juno ;
    dJ2Target = dJ2Juno ; 
    dJ4Target = dJ4Juno ; 
    dJ6Target = dJ6Juno ;
    dJ8Target = dJ8Juno ;
    dJ10Target= dJ10Juno;

    J3Target  = J3Juno ;
    J5Target  = J5Juno ; 
    J7Target  = J7Juno ;
    J9Target  = J9Juno ;
    dJ3Target = dJ3Juno; 
    dJ5Target = dJ5Juno; 
    dJ7Target = dJ7Juno;
    dJ9Target = dJ9Juno;
 
    //    warning("Overwriting uncertainties of J8 and J10 with that of J6 to have equal optimization weights.");
    //    dJ8Target  = dJ6Target;
    //    dJ10Target = dJ6Target;
    if (print) cout << " Using the 1-sigma uncertainties from Durante (2020) as they are." << endl;
    */  

    SetNewtonStepSizeToStandardValue(print);
    Precompute();

    // define radii of every layer on the grid
    double eps = 1.e-6;
    zeta = 1.0-eps;
    CalculateAndSetAllSpheroidVolumes(); // added on 5/24/21

    densitiesAlreadyInitialized = false;

    //    hasCore=false;
    //    rCore  = 0.0;
  }

  void SetRotationPeriod(const double h, const double m, const double s) {
    double P = 3600.0*h + 60.0*m + s;
    SetRotationPeriod(P);
  }
  void SetRotationPeriod(const double P) { // in seconds
    rotationPeriod = P;
    double hours   = floor(rotationPeriod / 3600.0);
    double minutes = floor((rotationPeriod - 3600.0*hours) / 60.0);
    double seconds = rotationPeriod - 3600.0*hours - 60.0*minutes;
    cout << " Rotation period= " << hours << ":" << minutes << ":" << seconds << " = " << rotationPeriod << " seconds." << endl;
    //    Reset(); // required to update qRot for example

    omega_SI       = 2.0*pi/rotationPeriod;
    qRot = sqr(omega_SI)*cube(a_SI) / GM_SI; // same when we switch from cgs to SI
    //    if (print) cout << "Rotational parameter: "; Write(qRot);
    cout << " Rotational parameter: "; Write(qRot);
  }

  void PrintRotationPeriod(const double P) {
    double hours   = floor(P / 3600.0);
    double minutes = floor((P - 3600.0*hours) / 60.0);
    double seconds = P - 3600.0*hours - 60.0*minutes;
    cout << " Rotation period= " << hours << ":" << minutes << ":" << seconds << " = " << P << " seconds." << endl;
  }
  void PrintRotationPeriod() {
    PrintRotationPeriod(rotationPeriod);
  }

  void SetReferenceRadiusInKMAndPressureInBar(const double a_km, const double PReferenceInBar) {
    a_SI = a_km * 1.0e3;                  // 1 km = 10^3 m
    PTarget = PressureBarToPU(PReferenceInBar);
    pressure[0] = PTarget;
    //    Write(PTarget);
    cout << " R_eq= " << a_km << " km." << endl;
    cout << " P_target= " << PTarget << " PU = " << PressurePUToBar(PTarget) << " bar." << endl;
  }
  double SurfacePressurePU() const {
    return pressure[0];
  }
  double SurfacePressureAU() const {
    return PressurePUToAU( SurfacePressurePU() );
  }

  void SetNewtonStepSizeInitiallyToReducedValue(const bool print=true) {
    //    newtonStepSize=min(1.0,512.0/nl);
    newtonStepSize=0.25;
    if (print) Write(newtonStepSize);
  }
  void SetNewtonStepSizeToStandardValue(const bool print=true) {
    newtonStepSize=1.0;
    if (print) Write(newtonStepSize);
  }

  double ZetaDifference(const Array2 <double> & z1, const Array2 <double> & z2) const;
  double ZetaMin() const {
    return 0.8;
  }
  double ZetaMax() const {
    return 1.2;
  }
  void CheckZetaLimits() const;
  void PerformRegulaFalsiStep(const int iter, const Array2 <double> & zetaPrev, Array2 <double> & zeta, 
			      const Array2 <double> & equiPotPrev, const Array2 <double> & equiPot) const;

  int MaxSpheroidIntersections(const Array1 <double> & r) const;
  bool CheckForSpheroidIntersections(const Array1 <double> & r) const;
  bool CheckAndRepairSpheroidIntersections(Array2 <double> & zeta, const Array2 <double> & zetaPrev) const;
  //  void RepairSpheroidIntersections(Array1 <double> & r) const;
  //  void RepairSpheroidIntersections(const Array1 <double> & r, Array1 <double> & rr, const int i0, 
  //				   bool & flag, double & d2) const;

  void Precompute();
  void CopyMoments() {
    jTotalPrev = jTotal;
  }
  void UpdateZeta(const Array2 <double> & zetaOld, Array2 <double> & zetaNew, Array2 <double> & equiPotNew);
  void UpdateZetaAccelerated(const Array2 <double> & zetaOld, Array2 <double> & zetaNew, Array2 <double> & equiPotNew);

  //  double GetZChange(const double M, const double dMdZ) const;
  //  bool ChangeZ1(const double M, const double dMdZ1, MCParameters & p);
  //  bool ChangeZ2(const double M, const double dMdZ2, MCParameters & p);
  //  bool ChangeZ12(const double M, const double dMdZ1, const double dMdZ2, MCParameters & p); 
  void CheckTotalPlanetMass(const bool print);
  void NormalizeRhoToMatchTotalPlanetMass(MCParameters & p);

  double CalculateMass() const;
  double CalculateSpheroidVolume(const int j)             const; // Calculate the volume enclosed by the surface of layer 'j'. See Moments.C
  void CalculateAndSetAllSpheroidVolumes() { // added only on 3/3/2020
    for(int j=0; j<nl; j++) 
      volume[j] = CalculateSpheroidVolume(j);
    volume[nl] = 0.0; // so hat we can determine layerVolume[l=nl-1]=volume[l]-volume[l+1] as usual
    //    layerVolume.Resize(nl<);
    //    for(int j=0; j<nl-1; j++) 
    //      layerVolume[j] = volume[j] - volume[j+1];
    //    layerVolume[nl-1] = volume[nl-1];
  }
  double GetSpheroidVolume(const int j) const { // may be called for j=nl (a point inside core). The core is spheroid l=nl-1
    return volume[j];
  }
  //  double CalculateCoreVolume() const { // removed on 5/29/21 because we have more core layers now
  //    return CalculateSpheroidVolume(nl-1);
  //  }
  double CalculateSphericalizedSpheroidRadius(const int j) const {
    double V = GetSpheroidVolume(j); // 4/3 Pi r^3 = V   
    double rS = pow( 3.0*V / (4.0*pi), 1.0/3.0 );
    return rS;
  }
  //  double GetCoreVolume() const { // removed on 5/29/21 because we have more core layers now
  //    return volume[nl-1]; // requires volumes[] to be set
  //  }
  double CalculateLayerVolume(const int j) const {
    if (j==nl-1) return CalculateSpheroidVolume(nl-1);
    return CalculateSpheroidVolume(j) - CalculateSpheroidVolume(j+1); 
  }
  double GetLayerVolume(const int j) const { // add only on 3/3/2020, requires volumes[] to be set
    return GetSpheroidVolume(j) - GetSpheroidVolume(j+1);  // no also works for core with j=nl-1
  }
  double GetJ2HatDifference(const int j) const {
    return j2Hat[j] - j2Hat[j+1];
  }

  double CalculateLayerMass(const int j) const {
    return CalculateLayerVolume(j)*rho[j];
  }
  double GetLayerMass(const int j) const {  // add only on 3/3/2020, requires volumes[] to be set
    return GetLayerVolume(j)*rho[j];
  }
  double CalculateSpheroidMass(const int j) const { // Calculate the mass enclosed by the surface of layer 'j'. Could be made more efficient
    double m=0.0;
    for(int l=nl-1; l>=j; l--) {
      m += CalculateLayerMass(l);
    }
    return m;
  }
  double GetSpheroidMass(const int j) const { // requires volumes[] to be set. Yields zero for j==nl
    double m=0.0;
    for(int l=nl-1; l>=j; l--) {
      m += GetLayerMass(l);
    }
    return m;
  }
  //  void UpdateCoreMass() {
  //    mCore = CalculateCoreMass();
  //  }
  void CalculateMoments();                        // in Moments.C
  void CalculateUnnormalizedMoments();            // in Moments.C
  void CalculateNormalizedMoments();              // in Moments.C
  double CalculateVolumeIntegral(const int mode);         // in CMS.C
  double CalculateMomentOfInertia() const;                // in CMS.C
  double CalculateMomentOfInertiaAgain() const;           // in CMC.C
  double CalculateMomentOfInertiaAgainAndAgain() const;   // in CMC.C

  void CalculateRhoFromDeltaRho() {
    double rho_cur = 0.0;
    for(int j=0; j<nl; j++) {
      rho_cur += delta[j];
      rho[j] = rho_cur;
    }
  }
  void CalculateDeltaRhoFromRho() {
    delta[0] = rho[0]; 
    for(int j=1; j<nl; j++) {
      delta[j] = rho[j] - rho[j-1];
    }
  }
  void UpdatePressure() {
    Check("q1");
    CalculateSurfacePotentials(); // sets potential[i] by calculating the potential in the equatorial plane
    Check("q2");
    CalculateSurfacePressures();  // pressure[i] = pressure[i-1] + rho[i-1] * ( potential[i] - potential[i-1] ) -- takes no time
    Check("q3");
  }
  void CalculateSurfacePressures() {
    for(int i=1; i<nl; i++) {
      pressure[i] = pressure[i-1] + rho[i-1] * ( potential[i] - potential[i-1] );
    }
  }
  void DetermineSurfacePressures(const Array1 <double> & pressureInFile) {
    double dPMax = 0.0;
    for(int j=0; j<nl; j++) { 
      double scale =  (pressureInFile[j] + pressure[j]) / 2.0;
      double dP =  abs(pressureInFile[j] - pressure[j]) / scale;
      if (j==0 || dP>dPMax) dPMax = dP;
      //    Write2(j,dP);
    }
    cout << " Pressure accuracy= " << dPMax << endl;
    //  Write3(iter,PressurePUToBar(pressure[0]),PressurePUToBar(pressure[jReference]));
  }
  void CalculateSurfacePotentials() { // no acceleration employed
    potential[0] = ExternalPotential(1.0,0.0); // call for equator zeta0=1.0 and mu=0
    for(int i=1; i<nl; i++) {
      potential[i] = InternalPotential(i,1.0,0.0); // call for equatorial plane and mu=0
      //    Write3(i,InternalPotentialSlow(i,1.0,0.0),InternalPotential(i,1.0,0.0));
    }
  }
  void CheckSurfacePotentials();
  double CentralPot() const {
    double U = 0.0;
    for(int i=0; i<nl; i++) {
      U -= jTildeP(i,0) / lambda[i];
    }
    return U;
  }
  double CentralPressure() const {
    double PC = pressure[nl-1] + rho[nl-1] * ( CentralPot() - potential[nl-1] );
    return PC;
  }
  double AveragePressure(const int i) const {
    double nextPressure = (i<nl-1) ? pressure[i+1] : CentralPressure();
    double averagePressure = 0.5 * ( pressure[i] + nextPressure );
    return averagePressure;
  }
  double AveragePressureGPa(const int i) const {
    return PressurePUToGPa(AveragePressure(i));
  }
  double AverageLambda(const int i) const {
    double nextLambda    = (i<nl-1) ? lambda[i+1] : 0.0;
    double averageLambda = 0.5 * ( lambda[i] + nextLambda );
    return averageLambda;
  }

  void SetLambdas();
  void SetLambdasWeighted();
  void SetLambdasSmoothlyWeighted();
  void SetLambdasBasedOnReferenceCalculation();
  void SetLambdasBasedOnSplinesFromReferenceCalculation();
  void ScaleLambdaGrid();
  void PrepareAdjustmentOfLambdaGrid(const Array1 <double> & rr);
  int LayerBeginningIndex(const int layer) const {
    return indexBoundaries[layer];
  }
  int LayerEndingIndex(const int layer) const {
    return indexBoundaries[layer+1];
  }
  int LayerSpheroidNumvber(const int layer) const {
    return LayerEndingIndex(layer) - LayerBeginningIndex(layer);
  }
  int GetLayerIndexFromSpheroidIndex(const int j) const {
    int m=0;
    for( ; m<mLayers; m++) {
      if (j<indexBoundaries[m+1]) break;
    }
    return m;
  }
  void PrintLambdaGrid() const {
    for(int j=0; j<nl; j++) {
      Write2(j,lambda[j]);
    }
  }
  void PrintDeltaLambda() const {
    for(int j=0; j<nl-1; j++) {
      Write3(j,lambda[j],lambda[j]-lambda[j+1]);
    }
  }
  void CheckLambdaGrid() const {
    if (CheckForNan(lambda)) error("NaNs found in lambda grid",lambda);
    for(int j=0; j<nl-1; j++) {
      if (lambda[j]<=lambda[j+1]) error("problem in lambda grid",j,lambda[j],lambda[j+1]);
    }
  }
  //  void SetUpEOSTable(const string & filename);
  //  void SetUpEOSTable(const string & filename, const MCParameters & p);
  //  void SetUpEOSTable(MCParameters & p);

  void SetUpOmegaSplinesAndRotationalPotential(const MCParameters & p, const bool print=false);
  double Evaluate(MCParameters & p, const bool matchPlanetMass, const bool print);
  double Chi2(const bool print) const;

  void WriteLayerCMSDataToFile(const string & fileName) const;
  //  void WriteLayerEOSDataToFile(const string & fileName, const bool printFlag=false);
  void WriteHarmonicsToMachineReadableFile(const string & fileName);
  void WriteBarotropeToFile(const string & fileName) const;
  void WriteMuGridToFile(const string & fileName) const;
  void WriteXiDataToFile(const string & fileName) const;
  void WriteGraphicsFileForOneQuarter(const string & fileName);
  void WriteGraphicsFileForOneHemisphere(const string & fileName);
  //  void PrintReferenceConditions() const;
  void ReadFilesAndPerformCMSCalculations(const string & prefix="");

  double DeltaLambda(const int j) const {
    double dLambda = (j==nl-1) ? lambda[j] : lambda[j] - lambda[j+1]; 
    return dLambda;
  }
  double DeltaJnHat(const int j, const int n) const {
    double dJnHat = (j==nl-1) ? jReg(j,n)/delta[j] : jReg(j,n)/delta[j] - jReg(j+1,n)/delta[j+1];     
    return dJnHat;
  }
  double WnFunction(const int j, const int n) const {
    return DeltaJnHat(j,n) / DeltaLambda(j);
  }

  void SuperAdiabaticityTest();

  void PrintPolarRadiusShort();
  void PrintPolarRadiusLong();

  void EOSPolytrope() {
    SetInitialDensitiesToPolytropeEOS();
  }
  void SetInitialDensitiesToPolytropeEOS();
  void ReadInitialDensitiesFromFile();

  void UpdateDensity() {
    UpdateDensityBarotrope();
    CalculateDeltaRhoFromRho();
  }
  void UpdateDensityBarotrope(); // see EOS.C
  void UpdateDensityPolytrope();
  void UpdateDensityEOSTable();
  void UpdateDensityEOSTableAndMatchCoreMass();

  double KPolytrope() const;

  double ExternalPotential(const double zeta0, const double mu);
  double ExternalV(const double zeta0, const double mu);
  double ExternalQ(const double zeta0, const double mu);
  double InternalPotentialSlow(const int jLayer, const double zetaJ, const double mu);
  double InternalPotential(const int jLayer, const double zetaJ, const double mu);

  void DetermineAllEquiPotentials();

  void SetGM_SI(const double GM_SI_) {
    GM_SI = GM_SI_;
    M_SI  = GM_SI / PC::G;
  }
   void SetM_SI(const double M_SI_) {
    M_SI = M_SI_;
    GM_SI = M_SI * PC::G;
  }
 
  double DensitySIToPU(const double rhoSI) const {
    double rhoPU = rhoSI * cube(a_SI) / M_SI;
    return rhoPU;
  }
  double DensityPUToSI(const double rhoPU) const {
    double rhoSI = rhoPU / cube(a_SI) * M_SI;
    return rhoSI; // kg/m^3
  }
  double DensityPUToGCC(const double rhoPU) const { 
    return 0.001*DensityPUToSI(rhoPU); // g/cc = 10^-3 kg/m^3
  }
  double DensityGCCToPU(const double rhoGCC) const { // g/cc = 10^-3 kg/m^3
    return DensitySIToPU(1000.0*rhoGCC);
  }

  double PressurePUToSI(const double pressurePU) const {
    double pressureSI = pressurePU * PC::G * sqr(M_SI) / sqr(sqr(a_SI));
    return pressureSI;
  }
  double PressureSIToPU(const double pressureSI) const {
    double pressurePU = pressureSI / PC::G / sqr(M_SI) * sqr(sqr(a_SI));
    return pressurePU;
  }
  double PressureBarToSI(const double pressureBar) const {
    double P_SI = pressureBar * 1.0e5; // 1 bar = 10^5 Pa
    //    double P_CGS  = P_Pa   * 10.0;     // 1 P_CGS = 0.1 Pa = 10^-10 GPa
    return P_SI;
  }
  double PressureSIToBar(const double pressureSI) const {
    return pressureSI / PressureBarToSI(1.0);
  }

  double PressureBarToPU(const double pressureBar) const {
    return PressureSIToPU(PressureBarToSI(pressureBar));
  }
  double PressurePUToBar(const double pressurePU) const {
    return pressurePU / PressureBarToPU(1.0);
  }
  double PressureAUToPU(const double pressureAU) const {
    return PressureBarToPU(pressureAU*PC::AUToBar);
  }
  double PressurePUToAU(const double pressurePU) const {
    return pressurePU / PressureAUToPU(1.0);
  }
  double PressurePaToPU(const double pressurePa) const {
    return PressureSIToPU(pressurePa);
  }
  double PressurePUToPa(const double pressurePU) const {
    return pressurePU / PressurePaToPU(1.0);
  }
  double PressureGPaToPU(const double pressureGPa) const {
    return PressureSIToPU(1e9*pressureGPa);
  }
  double PressurePUToGPa(const double pressurePU) const {
    return 1e-9*PressurePUToPa(pressurePU);
  }

  double MassPUToEarthMasses(const double m) const {
    double mME = m * M_SI / MEarth_SI;
    return mME;
  }
  double MassEarthMassesToPU(const double m) const {
    return m / MassPUToEarthMasses(1.0);
  }

  double GetJn(const int n) const {
    return jTotal[n];
  }
  double GetJnTarget(const int n) const {
    if (n== 2) return J2Target;
    if (n== 4) return J4Target;
    error("In GetJnTarget(n), n is out of range",n);
    return -1.0;
  }
  double GetJnDiscrepancyCMSMinusTarget(const int n) const { // positive result means the model value is larger than the observations
    return GetJn(n) - GetJnTarget(n); 
  }

  void PrintEvenJnCompact() const {
    cout << "Jn*10^6:";
    //    int p = cout.precision();
    //    cout.precision(6);
    for(int n=2; n<=np; n+=2) {
      if (n<=12) {
	cout << " J" << n << "= " << fix(10,4,1e6*GetJn(n));
	//	cout << " J" << n << "= " << fix(14,8,1e6*GetJn(n)); // extra many digits
      } else {
	cout << " J" << n << "= " << fix(10,6,1e6*GetJn(n));
	//	cout << " J" << n << "= " << fix(14,8,1e6*GetJn(n)); // extra many digits
      }
      //      int p = (n<=12) ? 4 : 6;
      //      cout << " J" << n << "= " << 1e6*GetJn(n);
      //      cout << " J" << n << "= " << fix(p+6,p,1e6*GetJn(n));
    }
    //    cout.precision(p);
    cout << endl;
  }
  void PrintEvenJnTarget() const {
    int n1 = 10;
    int n2 = 4;
    //    int n1 = 14;
    //    int n2 = 8;
    cout << " Target:        ";
    cout << " J" << 2 << "= " << fix(n1,n2,1e6*J2Target);
    cout << " J" << 4 << "= " << fix(n1,n2,1e6*J4Target);
    cout << endl;
  }

  void PrintJnTarget() const {
    PrintEvenJnTarget();
  }
  void PrintCMSResults() const;
  void PrintCMSModelInfo() {
  }

  bool CheckForNan(const double x) const { // if x==NaN return true
    // if x!=NaN then x!=x will be false
    // if x==NaN then x!=x will be true
    return (x != x); 
  }
  bool CheckForNan(const Array1 <double> & a) const { // if x==NaN return true
    for(int i=0; i<a.Size(); i++) 
      if (CheckForNan(a[i])) return true;
    return false;
  }
  bool CheckForNan(const Array2 <double> & a) const { // if x==NaN return true
    for(int i=0; i<a.dim[0]; i++) 
      for(int j=0; j<a.dim[1]; j++) 
	if (CheckForNan(a(i,j))) return true;
    return false;
  }
  /*
  bool CheckJnForNan() const {
    for(int n=0; n<=np; n+=2) {
      if (CheckForNan(GetJn(n))) return true;
    }
    return false;
  }
  bool CheckDeltaForNan() const {
    return CheckForNan(delta);
  }
  */
  void CheckJnForNan() const {
    //    cout << "Encountered NaNs in Jn." << endl;
    //    cout.flush();
    //    throw("Encountered NaNs in Jn.");
    for(int n=0; n<=np; n+=2) {
      if (CheckForNan(GetJn(n))) {
	//	error("Encountered NaN in Jn",n,GetJn(n));
	warning("Encountered NaN in Jn",n,GetJn(n));
	throw("Encountered NaN in Jn");
      }
    }
  }
  void CheckDeltaForNan() const {
    //    if (CheckForNan(delta)) error("Encountered NaNs in delta array",delta);
    if (CheckForNan(delta)) throw("Encountered NaNs in delta.");
  }

  void Check(const string & s) const {
    return;
    //    cout << s << ":" << endl;
    //    cout << s << ":"; Write7(zeta.Sum(),potential.Sum(),pressure.Sum(),rho.Sum(),delta.Sum(),volume.Sum(),lambda.Sum());
    //    cout << s << ":"; Write2(MassPUToEarthMasses(GetCoreMass()),MassPUToEarthMasses(CalculateCoreMass()));
    //    cout << s << ":"; Write2(MassPUToEarthMasses(GetCoreMass()),CalculateMass());
    string ss;
    //    if (CheckForNan(Pn_mu))     ss += " Pn_mu";
    //    if (CheckForNan(zeta))      ss += " zeta";
    if (CheckForNan(jTotal))    ss += " jTotal";
    //    if (CheckForNan(jTilde))    ss += " jTilde";
    //    if (CheckForNan(jReg))      ss += " jReg";
    //    if (CheckForNan(jRegP))     ss += " jRegP";
    //    if (CheckForNan(zeta))      ss += " zeta";
    if (CheckForNan(potential)) ss += " potential";
    if (CheckForNan(pressure))  ss += " pressure";
    if (CheckForNan(delta))     ss += " delta";
    if (CheckForNan(rho))       ss += " rho";
    //    if (CheckForNan(volume))    ss += " volume";
    //    if (CheckForNan(lambda))    ss += " lambda";
    if (ss.length()>0) {
      //      Write(delta);
      //      Write(rho);
      error("Encountered NaNs in",ss);
      warning("Encountered NaNs in",ss);
      throw("Encountered NaNs in"+ss);
    }
  }

  friend ostream& operator<<(ostream & os, const CMS & cms) {
    os << endl;
    os << " Model output: " << endl;
    os << endl;
    
    Form f(sci,16,10);
    f.SetShowPos();
    for(int i=0; i<=cms.np; i+=2) {
      os << " J " << fix(2,0,i) << " = " << f(cms.jTotal[i]) << endl;
    }
    return os;
  }

  void SaveVector(const string & fileName, const Array1 <double> & a, const bool txt) const {
    int p = 10;
    if (txt) ::SaveVectorText(fileName+".txt",a,false,p);
    else ::SaveVector(fileName+".dat",a);
  }

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

  // pass in l/a0 with a0 equatorial radius of planet
  // return omega / omega0 with omega is stored in qRot = a0^3 omega0^2 / GM
  double Omega(const double l) const { 
#if defined DIFF_ROT && !defined TWE
    return omegaSpline.Interpolate(l);
#else
    return 1.0;
#endif
  }

  void Omega(const double l, double & w, double & dw) const { 
#if defined DIFF_ROT && !defined TWE
    w  = omegaSpline.Interpolate(l);
    dw = omegaSpline.Derivative(l);
#else
    w  = 1.0;
    dw = 0.0;
#endif
  }

  double RotationalPotential(const double l) const { 
#if defined DIFF_ROT && !defined TWE
    return qRot * rotPotSpline.Interpolate(l);
#else
    return qRot * 0.5*l*l;
#endif
 }

  void RotationalPotential(const double l, double & pot, double & dPot) const { 
#if defined DIFF_ROT && !defined TWE
    pot  = qRot * rotPotSpline.Interpolate(l);
    dPot = qRot * rotPotSpline.Derivative(l);
#else
    pot  = qRot * 0.5*l*l;
    dPot = qRot * l;
#endif
  }

  void SetUpPressureCumulativeMassSpline(); // needed to determine helium mass fraction Y2 as function of Y1.
  void SetUpRadiusCumulativeMassSpline();   // needed to determine helium mass fraction Y2 as function of Y1.
  void SetUpPressureRadiusSplines();
  double AverageZ2MetallicLayer(const MCParameters & p);

  /////////////////////////// terms for thermal wind equation /////////////////////////

  /*
  double KmToPlanetaryUnits(const double l_km) const {
    double l_cgs = l_km * 1000.0 * 100.0;
    return l_cgs / a_cgs;
  }
  double PlanetaryUnitsToKm(const double l_PU) const {
    return l_PU / KmToPlanetaryUnits(1.0); 
  }
  */
  double LengthPUToSI(const double l_PU) const {
    return l_PU*a_SI;
  }
  double LengthPUToCGS(const double l_PU) const {
    return 100.0*LengthPUToSI(l_PU);
  }
  double LengthSIToPU(const double l_SI) const {
    return l_SI/a_SI;
  }
  double LengthKmToPU(const double l_km) const {
    return 1000.0*l_km/a_SI;
  }
  double LengthPUToKm(const double l_PU) const {
    return 0.001*LengthPUToSI(l_PU);
  }

  double PotentialSIToPU(const double uSI) const {
    //    double GM_SI = GM_cgs * 1e-6;
    //    double a_SI  = 0.01*a_cgs;
    double uPU   = a_SI / GM_SI * uSI;
    return uPU;
  }
  double PotentialPUToSI(const double uPU) const {
    return uPU / PotentialSIToPU(1.0);
  }

  inline Point XiVector(const int l, const int imu) const;
  inline double DensityDerivativeInPolarDirection(const int l, const bool north);
  inline double DensityDerivativeOnGridPoint(const int l, const int imu);
  inline double DensityDerivativeOnEquator(const int l);
  Spline BuildDRhoDsSpline(const int l, const bool north);
  void SetXiSurfaceSplines();
  void DeriveDensityCorrectionOnAllGridPoints();
  void BuildRhoAndXiSplineForArbitryMu(const double mu, Spline & xiSpline, Spline & rhoSplineS, Spline & rhoSplineN) const;
  void BuildRhoAndXiSplineForArbitryMu(const double mu, Spline & xiSpline, Spline & rhoSplineS, Spline & rhoSplineN, Spline & dRhoSplineS, Spline & dRhoSplineN) const;

  double DeriveMuGivenL(const double L) const; // xiTWE must have been set by CMS::DeriveDensityCorrectionOnAllGridPoints() 
  double DeriveThetaGivenL(const double L) const {
    double mu = DeriveMuGivenL(L);
    double theta = acos(mu);
    return theta;
  }

  void SetIterMin(const int n) {
    iterMin = n;
    Write(iterMin);
  }

  void BuildXiMuSplines();

  void WriteNodesFileHeader(ofstream & of, const Units & units);
  double WriteNodeInfoOnePoint(const double x, const double y, const double z, const int material, ofstream & of, Units & u, const bool print=false);
  double WriteNodeInfoOnePoint(const Point & r, const int material, ofstream & of, Units & u, const bool print=false) {
    return WriteNodeInfoOnePoint(r[0],r[1],r[2],material,of,u);
  }

  void TestGravityGradientsOnePoint(const double x, const double y, const double z, const int material, ofstream & of, Units & u);
  void WriteNodesFiles(const string & nodesOutputFile, const bool nodesXDirections, const bool nodesZDirections, const bool SIUnits);

};

class Units {
 public:
  bool SIFlag;
  double fL; // factor to convert from PU to desired units (PU or SI)
  double fRho;
  double fP;
  double fPot;
  double fPotGrad;
  double fOmega;

  double uL; // current unit in SI units
  double uRho;
  double uP;
  double uPot;
  double uPotGrad;
  double uOmega;

  Units(CMS & cms, const bool SIUnits_) {
    SIFlag   = SIUnits_;
    fL       = (SIFlag) ? cms.LengthPUToSI(1.0)                          : 1.0; 
    //  fRho = (SIFlag) ? cms.DensityPUToCGS(1.0)                        : 1.0; // use g/cc not kg/m^3
    fRho     = (SIFlag) ? cms.DensityPUToSI(1.0)                         : 1.0; // use kg/m^3
    fP       = (SIFlag) ? cms.PressurePUToPa(1.0)                        : 1.0; // use Pa
    fPot     = (SIFlag) ? cms.PotentialPUToSI(1.0)                       : 1.0;
    fPotGrad = (SIFlag) ? cms.PotentialPUToSI(1.0)/cms.LengthPUToSI(1.0) : 1.0;
    fOmega   = (SIFlag) ? cms.omega_SI                                   : 1.0;

    uL       = (SIFlag) ? 1.0 : cms.LengthPUToSI(1.0);
    uRho     = (SIFlag) ? 1.0 : cms.DensityPUToSI(1.0);
    uP       = (SIFlag) ? 1.0 : cms.PressurePUToPa(1.0);
    uPot     = (SIFlag) ? 1.0 : cms.PotentialPUToSI(1.0);
    uPotGrad = (SIFlag) ? 1.0 : cms.PotentialPUToSI(1.0)/cms.LengthPUToSI(1.0);
    uOmega   = (SIFlag) ? 1.0 : cms.omega_SI; 
  }
};

#endif // _CMS_