//*************
// tov_lorenepars_polytrope.C
// This interactive code calculates the maximum mass for the sequence of
// polytropic models calculated by tov_solver_polytrope.
// Next, it prompts for a mass and radius, determines if such a configuration
// is on the stable branch for that EOS (below the maximum baryon mass),
// and outputs parameters for use with lorene using simple linear interpolation.
//
// The call is just ./executable, but you will need a file in the form
// gam_n.nnn.dat, where "n.nnn" is the value of gamma.
//
// The output is a file par_eos.d.gam_n.nnn that may be renamed and used
// with the Lorene code
//
// (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;

int nlines = 800;

using namespace std;

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

  cout<<"Enter gamma value"<<endl;

  double gamval;
  cin >> gamval;


  char filename[20];
  sprintf(filename,"gam_%.3f.dat",gamval);

  cout<<"filename? "<<filename<<endl;
  ifstream tovdatfile;
  tovdatfile.open(filename);

  char outfilename[30];
  sprintf(outfilename,"par_eos.d.poly_%.3f.dat",gamval);
  ofstream loreosfile;
  loreosfile.open(outfilename);

  double* r = new double[nlines];
  double* massg = new double[nlines];
  double* massb = new double[nlines];
  double* rhoc = new double[nlines];
  double compact,rval,massgval,massbval,rhoval;
  double massold=0.0;
  int i,imax;
  int done=0;
  for(i=0; i<nlines && done==0; i++) {
    tovdatfile >> rval >> massgval >> massbval >> compact >> rhoval;
    if(massbval > massold) {
      r[i]=rval;
      massg[i]=massgval;
      massb[i]=massbval;
      rhoc[i]=rhoval;
      massold=massbval;
    } else {
      done=1;
      imax=i-1;
    }
  }
  tovdatfile.close();

  cout<<"Maximum mass configuration for k=1:"<<endl;
  cout<<"r:"<<r[imax]<<" mg and mb:"<<massg[imax]<<" "<<massb[imax]<<" rhoc:"<<rhoc[imax]<<" --> Compactness:"<<massg[imax]/r[imax]<<endl;

  double Cmax = massg[imax]/r[imax];

  cout<<"What (gravitational) mass do you want (in Msun)?"<<endl;
  double mns,rns;
  cin>>mns;
  cout<<"What radius do you want (in km)?"<<endl;
  cin>>rns;

  double Cns = ggrav*mns*msun/rns/kilometer/speedlight/speedlight;
  cout<<"That yields a compactness of "<<Cns<<endl;

  if(Cns>Cmax) {
    cout<<" That is larger than the maximum allowed compactness, Cmax="<<Cmax<<"!!!  Better luck next time!"<<endl;
    return -2;
  } else {
    cout<<"We can do that!"<<endl;

    done=0;
    double Cnew,Cold = massg[0]/r[0];
    for (i=1; i<imax && done==0; i++) {
      Cnew=massg[i]/r[i];
      if(Cnew<Cns) {
	Cold=Cnew;
      } else {
	//Interpolate!
	double interp = (Cns-Cold)/(Cnew-Cold);
	double radns = r[i-1]+interp*(r[i]-r[i-1]);
	double mgns = massg[i-1]+interp*(massg[i]-massg[i-1]);
	double mbns = massb[i-1]+interp*(massb[i]-massb[i-1]);
	double rhons = rhoc[i-1]+interp*(rhoc[i]-rhoc[i-1]);

	cout<<"Dimensionless values:"<<endl;
	cout<<"R="<<radns<<" Massg:"<<mgns<<" Massb:"<<mbns<<" rhoc:"<<rhons<<endl;

	cout<<"Compactness check! "<<Cns<<" "<<mgns/radns<<endl;

	//Scale to solar masses
	double scale = mns/mgns;
	cout<<"Scaling factor:"<<scale<<endl;
	double rhocscale = rhons/scale/scale;
	double kappa = pow(scale,2.0*gamval-2.0);

	cout<<"Rescaled density and kappa:"<<rhocscale<<" "<<kappa<<" porcheck:"<<pow(rhons,gamval-1.0)<<" "<<kappa*pow(rhocscale,gamval-1.0)<<endl;
	
	//Rescale for Lorene
	double radlor = rns/10;
	double mlor = 11.9789*mns;
	double kaplor = kappa*pow(3720.8,1.0-gamval);
	double denslor = rhocscale*3720.8;
	cout<<"In Lorene units, we want a radius and grav. mass of "<<radlor<<" "<<mlor<<endl;

	cout<<"Lorene compactness check:"<<0.0123265*mlor/radlor<<endl;
	cout<<"Lorene por check:"<<pow(rhons,gamval-1.0)<<" "<<kappa*pow(rhocscale,gamval-1.0)<<" "<<kaplor*pow(denslor,gamval-1.0)<<endl;

	cout<<"Lorene kappa value:"<<kaplor<<endl;
	cout<<"Lorene log enth:"<<log(1.0+gamval/(gamval-1.0)*kaplor*pow(denslor,gamval-1.0))<<endl;
	
	
	done=1;

	loreosfile << "1    Type of the EOS (cf. documentation of Eos::eos_from_file)" <<endl;
	loreosfile << "Star 1  EOS " <<endl;
	loreosfile << gamval<<"        adiabatic index gamma" <<endl;
	loreosfile << kaplor << "      pressure coefficient kappa [rho_nuc c^2 / n_nuc^gamma]"    <<endl;
	loreosfile << "1.              mean particle mass [m_b]"<<endl; 
      }
    }

  }
  tovdatfile.close();
  loreosfile.close();
  return 0;
}