//*************
// tov_solver_polytrope.C
// This code generates a sequence of models for polytropes using the TOV
// equations for NS structure.
// Note that "polytrope" here implies P=kappa * rho^gamma
// where "rho" is the REST mass density, not the energy density
// Here, we assume kappa=1, since results may be scaled to any value
//
// The call is just ./executable gam   where "gam" is the value of gamma above
//
// The output is a file gam_n.nnn.dat, where "n.nnn" is the value of gamma.
//
// (c)2019 Josh Faber, please consider this either GPL
//  or CC-BY-SA as far as licensing goes.
//***************************************************8

#include <cstdlib>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <string>
#include <stdlib.h>
#include <math.h>

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, double k, double gam);
void derivs(double x, double* v, double* y, double k, double gam);

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


  double gam = atof(argv[1]);


  char filename[20];
  sprintf(filename,"gam_%.3f.dat",gam);
  ofstream outfile;
  outfile.open(filename);
  

  double step=1.0e-6;
  int np = (int) (100.0/step);

  double* rho = new double [np];
  double* mass = new double [np];
  double* massb = new double [np];

  int i,j;

  //initial guess
  double k=1;
  double compact;

  int nloop;
  int nloopmax=800;
  double rhoc=1.0e-3;
  
  for (nloop=0; nloop<nloopmax; nloop++) {
    
    for(i=0; i<np; i++) {
      rho[i] = 0.0;
      mass[i]=0.;
      massb[i]=0.;
    }
    
    rho[0]=rhoc;
    mass[0]=0.0;
    massb[0]=0.0;

    //    cout<<"Enthalpy check:"<<log(1.0+gam/(gam-1.0)*k*pow(rho[0],gam-1.0))<<endl;
    
    double y[3];
    
    int done=0;
    double xsurf;

     double rrk4,mrk4,mbrk4;
    
    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,k,gam);
      
      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]);
	cout<<"r:"<<xsurf<<" massg:"<<mass[i]<<" "<<" massb:"<<massb[i]<<" compact:"<<mass[i]/xsurf<<" rhoc:"<<rhoc<<endl;
	outfile<<xsurf<<" "<<mass[i]<<" "<<massb[i]<<" "<<mass[i]/xsurf<<" "<<rhoc<<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;

	//	cout<<"plot:"<<xsurf<<" "<<mass[i]/xsurf<<" "<<k*pow(rho[i],gam-1.0)<<endl;

      }
      
      compact = mass[i]/xsurf;

    }

    rhoc*=1.01;
  }

  outfile.close();
  
  return 0;

}

void rk4(double x, double* y, double step, double k, double gam) {

  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,k,gam);
  for (i=0; i<3; i++) {
    k1[i]=step*v[i];
    t1[i] = y[i]+0.5*k1[i];
  }
  if(t1[0]<0) {
    done=1;
    y[0]=0.;
    y[1]=0.;
  }
  xtry+=h;
  derivs(xtry,v,t1,k,gam);
  for (i=0; i<3; i++) {
    k2[i]=step*v[i];
    t2[i] = y[i]+0.5*k2[i];
  }
  if(t2[0]<0) {
    done=1;
    y[0]=0.;
    y[1]=0.;
  }
  derivs(xtry,v,t2,k,gam);	
  for (i=0; i<3; i++) {
    k3[i]=step*v[i];
    t3[i] = y[i]+k3[i];
  }
  if(t3[0]<0) {
    done=1;
    y[0]=0.;
    y[1]=0.;
  }
  xtry+=h;
  derivs(xtry,v,t3,k,gam);
  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, double k, double gam) {

  double press,dpdrho,eps;
  double rho = y[0];
  double mass = y[1];
  if(y[0]>0) {
    press = k*pow(rho,gam);
    eps = press/(gam-1.0);
    dpdrho = k*gam*pow(rho,gam-1.0);
  } else {
    press = 0;
    eps=0;
    dpdrho = 0;
  }

  if(x>0 && rho>0) {
      v[0] = -1.0/x/x*(rho+eps+press)*(mass+4.0*M_PI*x*x*x*press)/(1.0-2.0*mass/x)/dpdrho;
      v[1] = 4.0*M_PI*x*x*(rho+eps);
      v[2] = 4.0*M_PI*x*x*rho/sqrt(1.0-2.0*mass/x);
  } else {
      v[0]=0.;
      v[1]=0.;
      v[2]=0.;
  }
  
}