//-------------------------------------------------
// tov_solver.C:  Copyright 2019, Joshua Faber

// Licensing terms: please use similar terms to either the GPL or Creative Commons-BY-SA.  Thanks!

// To call this routine, look at the definitions in Read et al. (2008)-- PRD 79, 124032 (2009); arXiv:0812.2163
// Note: you likely have at least a moral obligation to cite their work if using their EOS model
//-------------------------------------------------
// The call is <executable> name_string P_1 gam_1 gam_2 gam_3

// where name_string is a string that is appended to the end of the filename for ID purposes
// P_1 is the logarithm of the cgs pressure at log(rho_cgs)=14.7
// and the 3 gamma's are the 3-part polytrope adiabatic indicies

// There are two files created:

// par_eos.d.name_string is a Lorene-compatible par_eos.d file

// densenthmassradius.name_string is a tabulation of central density, log(enthalpy), gravitational mass.
//                                                  radius in km, compactness, and baryon mass

// Finally the code will output to screen the configuration details for the model with the largest baryon mass, then stop.
//-------------------------------------------------
#include <cstdlib>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <string>
#include <stdlib.h>
#include <math.h>

double lowdens=1.0;

// Physical constants
double msun = 1.98847e33;
double kilometer=1.0e5;
double ggrav=6.67408e-8;
double speedlight=2.9979e10;

using namespace std;

void rk4(double x, double* y, double step, int npoly, double* kpoly, double* gam, double* rhob, double* aread);
void derivs(double x, double* v, double* y, int npoly, double* kpoly, double* gam, double* rhob, double* aread);
double enth(double rho, int npoly, double* kpoly, double* gam, double* rhob, double* aread);

int main(int argc, char** argv) {

  int i,j,npoly;
  double gam[7],kpoly[7],rhobreak[6],rhob[6],pressbreak[6],pressb[6],pcgs[6];

  string eosname=argv[1];
  string filename="par_eos.d."+eosname;
  string filename2="densenthmassradius."+eosname;

  ofstream outfile,outfile2;
  outfile.open(filename.c_str());
  outfile2.open(filename2.c_str());

  double presspoly=atof(argv[2]);

  //See read et al. for the underlying piecewise polytropic EOS model and low-density crust parameters
  
  npoly=7;
  
  gam[0]=1.58425;
  gam[1]=1.28733;
  gam[2]=0.62223;
  gam[3]=1.35692;

  gam[4]=atof(argv[3]);
  gam[5]=atof(argv[4]);
  gam[6]=atof(argv[5]);
  
  kpoly[0]=6.8011e-9;

  rhobreak[0]=7.3875;
  rhobreak[1]=11.5779;
  rhobreak[2]=12.4196;
  rhobreak[4]=14.7;
  rhobreak[5]=15.0;
  
  pressbreak[0]=log10(kpoly[0])+gam[0]*rhobreak[0];
  pressbreak[1]=pressbreak[0]+gam[1]*(rhobreak[1]-rhobreak[0]);
  pressbreak[2]=pressbreak[1]+gam[2]*(rhobreak[2]-rhobreak[1]);
  for(i=0; i<=2; i++)pcgs[i]=pressbreak[i]+2.0*log10(speedlight);
  
  double presscrust=pcgs[2];
  
  rhobreak[3]=(presspoly-presscrust+rhobreak[2]*gam[3]-rhobreak[4]*gam[4])/(gam[3]-gam[4]);
  pressbreak[3]=pressbreak[2]+gam[3]*(rhobreak[3]-rhobreak[2]);
  pressbreak[4]=pressbreak[3]+gam[4]*(rhobreak[4]-rhobreak[3]);
  pressbreak[5]=pressbreak[4]+gam[5]*(rhobreak[5]-rhobreak[4]);

  for(i=0; i<=5; i++)rhob[i]=pow(10.0,rhobreak[i]);
  for(i=0; i<=5; i++)pressb[i]=pow(10.0,pressbreak[i]);
  for(i=3; i<=5; i++)pcgs[i]=pressbreak[i]+2.0*log10(speedlight);

  for(i=1; i<=6; i++)kpoly[i]=kpoly[i-1]*pow(rhob[i-1],gam[i-1]-gam[i]);
  //  for(i=0;i<=5;i++)cout<<"gam,k,rho,p,pcgs,rhob:"<<gam[i]<<" "<<kpoly[i]<<" "<<rhobreak[i]<<" "<<pressbreak[i]<<" "<<pcgs[i]<<" "<<rhob[i]<<endl;

  double epsoverrho[6],aread[7];
  aread[0]=0.0;
  epsoverrho[0]=1.0+kpoly[0]*pow(rhob[0],gam[0]-1.0)/(gam[0]-1.0);
  for(i=1; i<=5; i++)epsoverrho[i]=epsoverrho[i-1]+kpoly[i]/(gam[i]-1.0)*(pow(rhob[i],gam[i]-1)-pow(rhob[i-1],gam[i]-1));
  for(i=1; i<=6; i++)aread[i]=epsoverrho[i-1]-1.0-kpoly[i]/(gam[i]-1.0)*pow(rhob[i-1],gam[i]-1.0);

  // The following create a par_eosN.d file for Lorene, where "N" should be 1 or 2
  outfile << "110    Type of the EOS (cf. documentation of Eos::eos_from_file)" <<endl;
  outfile << "Star 1  EOS " <<endl;
  outfile << "7                       number of polytropes" <<endl;
  outfile << gam[0] << "                 array of adiabatic index" <<endl;
  outfile << gam[1] <<endl;
  outfile << gam[2] <<endl;
  outfile << gam[3] <<endl;
  outfile << gam[4] <<endl;
  outfile << gam[5] <<endl;
  outfile << gam[6] <<endl;
  outfile << kpoly[0] <<endl; // k0
  outfile << pressbreak[0]<<endl;// log(P0)
  outfile << rhobreak[0] << "                 array of the exponent of fiducial densities logRho" <<endl;
  outfile << rhobreak[1] << endl;
  outfile << rhobreak[2] <<endl;
  outfile << rhobreak[3] <<endl;
  outfile << rhobreak[4] <<endl;
  outfile << rhobreak[5] <<endl;
  outfile << "0.                      array of percentage"<<endl;
  outfile << "0." <<endl;
  outfile << "0." <<endl;
  outfile << "0." <<endl;
  outfile << "0." <<endl;
  outfile << "0." <<endl;

  outfile.close();

  //Simple RK4 integrator
  double step=1.0e-5;
  int np = (int) (100.0/step);
  
  double* rho = new double [np];
  double* mass = new double [np];
  double* massb = new double [np];

  //Here is where we loop over central densities -- this is a good starting point as radii are still unphysically large
  // and masses too small, but they get better fast
  int iloop,nloop=300;
  double rhoc=3e14;


  //Initialize the check for increasing (baryon) mass with central density, necessary for stability
  double oldmass=-1.0;
  double newmass=0.0;
    
  for(iloop=1; iloop<=nloop && newmass>oldmass; iloop++) {
    for(i=0; i<np; i++) {
      rho[i] = 0.0;
      mass[i]=0.0;
      massb[i]=0.0;
    }
    rho[0]=rhoc;
    //    cout<<"check-- rhoc, h:"<<rhoc<<" "<<enth(rhoc,npoly,kpoly,gam,rhob,aread)<<endl;
    //     cout<<"Enthalpy check:"<<log(1.0+gam/(gam-1.0)*k*pow(rho[0],gam-1.0))<<endl;
    //  1 + gam P / (gam-1) rho    
    
    double y[3];
    
    int done=0;
    double xsurf,mrk4,mbrk4,rrk4,prevgmass,prevradius,prevrhoc;
    double kmfact=ggrav*msun/pow(speedlight,2)/kilometer;

    for (i=1; i<np && done==0; i++) {
      double xpt=(i-1)*step;
      y[0]=rho[i-1];
      y[1]=mass[i-1];
      y[2]=massb[i-1];
      
      rk4(xpt,y,step,npoly,kpoly,gam,rhob,aread);

      if(y[0]<=0) {
	done=1;
	rho[i] = 0.0;
	mass[i] = mass[i-1]+rho[i-1]*(mass[i-1]-mass[i-2])/(rho[i-2]-rho[i-1]);
	massb[i] = massb[i-1]+rho[i-1]*(massb[i-1]-massb[i-2])/(rho[i-2]-rho[i-1]);
	xsurf = xsurf+rho[i-1]*step/(rho[i-2]-rho[i-1]);

	// mass in Msun, radius in km, compactness, baryon mass in Msun, central density, log(enthalpy)
	double logenth=log(enth(rhoc,npoly,kpoly,gam,rhob,aread));
	if(massb[i]>newmass)outfile2<<rhoc<<" "<<logenth<<" "<<mass[i]<<" "<<xsurf*kmfact<<" "<<mass[i]/xsurf<<" "<<massb[i]<<endl;

	if(massb[i]<newmass)cout<<eosname<<" EOS: Maximum baryon mass:"<<newmass<<" grav mass:"<<prevgmass<<
			      " radius in km:"<<prevradius<<" rhoc:"<<prevrhoc<<endl;
	  
	mrk4=mass[i];
	mbrk4=massb[i];
	rrk4=xsurf;
	
      } else {
	rho[i]=y[0];
	mass[i]=y[1];
	massb[i]=y[2];
	xsurf = xpt+step;
	
      }
    }
    double compact = mass[i]/xsurf;


    //Updates
    oldmass=newmass;
    newmass=mbrk4;
    prevgmass=mrk4;
    prevradius = xsurf*kmfact;
    prevrhoc=rhoc;
    
    //Increase the central density; geometric spacing, factor of 1% increase step-to-step
    rhoc *= 1.01;

  }

  outfile2.close();
  
}




void rk4(double x, double* y, double step, int npoly, double* kpoly, double* gam, double* rhob, double* aread) {

  double h=step/2.0;    
  double t1[3], t2[3], t3[3], k1[3], k2[3], k3[3], k4[3], v[3], dv[3];
  
  int i;

  int done = 0;

  double xtry=x;

  derivs(xtry,v,y,npoly,kpoly,gam,rhob,aread);
  for (i=0; i<3; i++) {
    k1[i]=step*v[i];
    t1[i] = y[i]+0.5*k1[i];
  }
  if(t1[0]<lowdens) {
    done=1;
    y[0]=0.;
    y[1]=0.;
  }
  xtry+=h;
  derivs(xtry,v,t1,npoly,kpoly,gam,rhob,aread);
  for (i=0; i<3; i++) {
    k2[i]=step*v[i];
    t2[i] = y[i]+0.5*k2[i];
  }
  if(t2[0]<lowdens) {
    done=1;
    y[0]=0.;
    y[1]=0.;
  }
  derivs(xtry,v,t2,npoly,kpoly,gam,rhob,aread);
  for (i=0; i<3; i++) {
    k3[i]=step*v[i];
    t3[i] = y[i]+k3[i];
  }
  if(t3[0]<lowdens) {
    done=1;
    y[0]=0.;
    y[1]=0.;
  }
  xtry+=h;
  derivs(xtry,v,t3,npoly,kpoly,gam,rhob,aread);
  for (i=0; i<3; i++) {
    k4[i] = step*v[i];
  }

  for (i=0; i<3; i++) y[i] += (k1[i]+2*k2[i]+2*k3[i]+k4[i])/6.0;
  if(done==1)y[0]=0.;
}

void derivs(double x, double* v, double* y, int npoly, double* kpoly, double* gam, double* rhob, double* aread) {

  double press,dpdrho,eps;
  double rho = y[0];
  double mass = y[1];

  int i,nsegment=0;
  if(rho>0) {
    for (i=0; i<npoly-1; i++) {
      if(rho>rhob[i])nsegment+=1;
    }
    press=kpoly[nsegment]*pow(rho,gam[nsegment]);
    dpdrho=gam[nsegment]*press/rho;
    if(nsegment==0) {
      eps=rho+press/(gam[0]-1.0);
    } else {
      eps=rho*(1.0+aread[nsegment])+press/(gam[nsegment]-1);
    }
  } else {
    press = 0;
    eps=0;
    dpdrho = 0;
  }


  //We'll want to rescale away from density-based cgs units here.
  //If P and rho are DIVIDED by a factor f, then r and m must each be multiplied by a factor f^{1/2} to balance
  //dm/dr is dimensionless
  //dP/dr divides by f^{-3/2}
  double refdens=msun/pow(ggrav*msun/speedlight/speedlight,3.0);
  double rhoref=rho/refdens;
  double epsref=eps/refdens;
  double pressref=press/refdens;
  
  if(x>0 && rho>0) {
      v[0] = -1.0/x/x*(epsref+pressref)*(mass+4.0*M_PI*x*x*x*pressref)/(1.0-2.0*mass/x)/dpdrho;
      v[1] = 4.0*M_PI*x*x*(epsref);
      v[2] = 4.0*M_PI*x*x*(rhoref)/sqrt(1.0-2.0*mass/x);
  } else {
      v[0]=0.;
      v[1]=0.;
      v[2]=0.;
  }
  v[0]*=pow(refdens,1.0);
  
}

double enth(double rho, int npoly, double* kpoly, double* gam, double* rhob, double* aread) {

  int i,nsegment=0;
  double enthval;
  if(rho>0) {
    for (i=0; i<npoly-1; i++) {
      if(rho>rhob[i])nsegment+=1;
    }
    enthval=1.0+aread[nsegment]+kpoly[nsegment]*gam[nsegment]*pow(rho,gam[nsegment]-1.0)/(gam[nsegment]-1.0);
  } else {
    enthval=0.0;
  }
  return enthval;

}