/****************************************************************************/
/* This code was implemented for calculating nuclear EoS based on the       */
/* Walecka model with matter composed of n, p, e, muons, Lambda, Sigma      */
/* and Cascade                                                              */
/****************************************************************************/

#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_multiroots.h>

#define hc 197.327
#define nmax 250000
#define epsilon 1e-8
#define me 0.511/hc // electron mass
#define mmu 105.0/hc // muon mass
#define mp 938.3/hc // proton mass
#define mn 939.6/hc // neutron mass
#define ml 1116.0/hc // Lambda mass
#define msp 1189.0/hc // Sigma + mass
#define msm 1197.0/hc // Sigma - mass
#define ms0 1193.0/hc // Sigma 0 mass
#define mc0 1315.0/hc //Chi 0 mass
#define mcm 1321.0/hc //Chi - mass

struct rparams
{
  double gmw2; // (g_omega/m_omega)^2
  double gms2; // (g_sigma/m_sigma)^2
  double gmr2; // (g_rho/m_rho)^2
  double b;
  double c;
  double xs; // x_sigma
  double xw; // x_omega
  double xr; // x_rho
  double rho; //density
};

int Walecka_f (const gsl_vector * x, void *params, gsl_vector * f);     
double Eq1 (double x1, double x2, double x3, double x4, double x5, void *params);
double Eq2 (double x1, double x2, double x3, double x4, double x5, void *params);
double Eq3 (double x1, double x2, double x3, double x4, double x5, void *params);
double Eq4 (double x1, double x2, double x3, double x4, double x5, void *params);
double Eq5 (double x1, double x2, double x3, double x4, double x5, void *params);
double kf2(double mub, double mue, double m8, double Vw, double Vr,double qe, double I3, void *params);
double dens_s(double m8, double k);
double dens(double k);

int main (void)
{

  int status;
  double mue;
  double gs_s, gww, grr, b, c;
  double mub, Vw, Vr, Vs;
  double rhoaux;
  double gms2, gmw2, gmr2;
  double gss_i, mun_i, mue_i, gww_i, grr_i;

  //  double xs=1.0, xr=1.0, xw=1.0; //for G300
  double xs=0.9, xr=0.9, xw=0.9; //for GM4
  size_t iter = 0;
     
  const size_t n = 5;

  const gsl_multiroot_fsolver_type *T;
  gsl_multiroot_fsolver *s;

  rhoaux=0.01; //fm^-3
     
  /*  gmw2=4.733;  //for G300
  gms2=9.031;
  gmr2=4.825;
  b=0.003305;
  c=0.01529;
  */
  gmw2=7.149;  //for GM4
  gms2=11.79;
  gmr2=4.411;
  b=0.002947;
  c=-0.001070;

  gss_i=sqrt(300.0/(0.9999*939.0));
  gss_i=asin(gss_i);
  mun_i=mn;
  mue_i=0.12*938.0/hc*pow(rhoaux/0.17, 2.0/3.0);
  gww_i= sqrt(226.0/hc);
  grr_i=1000.0/hc*rhoaux;

  for (rhoaux=0.01; rhoaux<=1.7;){ 
 
    iter=0;

    struct rparams p = {gmw2, gms2, gmr2, b, c, xs, xw, xr, rhoaux};

    gsl_multiroot_function f = {&Walecka_f, n, &p};

    double x_init[5] = {gss_i, gww_i, grr_i, mun_i, mue_i};

    gsl_vector *x=gsl_vector_alloc(n);
     
    gsl_vector_set(x, 0, x_init[0]);
    gsl_vector_set(x, 1, x_init[1]);
    gsl_vector_set(x, 2, x_init[2]);
    gsl_vector_set(x, 3, x_init[3]);
    gsl_vector_set(x, 4, x_init[4]);
    
    T = gsl_multiroot_fsolver_hybrid;
    s = gsl_multiroot_fsolver_alloc(T, 5);
    gsl_multiroot_fsolver_set(s, &f, x);
    
    do
      {
	iter++;
	
	printf ("iter = %3zu x = %g %g %g %g %g " "f(x) = %g %g %g %g %g \n", iter,
		gsl_vector_get (s->x, 0),
		gsl_vector_get (s->x, 1),
		gsl_vector_get (s->x, 2),
		gsl_vector_get (s->x, 3),
		gsl_vector_get (s->x, 4),
		gsl_vector_get (s->f, 0),
		gsl_vector_get (s->f, 1),
		gsl_vector_get (s->f, 2),
		gsl_vector_get (s->f, 3),
		gsl_vector_get (s->f, 4));

	//      print_state (iter, s);
	status = gsl_multiroot_fsolver_iterate (s);
	
	if (status)   /* check if solver is stuck */
	  break;
	
	status = gsl_multiroot_test_residual (s->f, 1e-7);
      }
    while (status == GSL_CONTINUE && iter < 1000);
    
    printf ("status = %s\n", gsl_strerror (status));
    
    Vs=gsl_vector_get (s->x,0);
    Vw=gsl_vector_get (s->x,1);
    Vr=gsl_vector_get (s->x,2);
    mub=gsl_vector_get (s->x,3);
    mue=gsl_vector_get (s->x,4);
    
    //Mapeio:
    
    gs_s=0.9999*mn*sin(Vs)*sin(Vs);
    gww=Vw*Vw;
    grr=Vr;
        
    gsl_multiroot_fsolver_free (s);
    gsl_vector_free (x);

    gss_i=sqrt(gs_s/(0.9999*mn));
    gss_i=asin(gss_i);
    mun_i=mub;
    mue_i=mue;
    gww_i= sqrt(gww);
    grr_i=grr;
    
    if(rhoaux<0.1)
      rhoaux+=0.005;
    else
      rhoaux+=0.025;
  }

  return 0;
}

int Walecka_f (const gsl_vector * x, void *params, gsl_vector * f)
{  
  const double x0 = gsl_vector_get (x, 0);
  const double x1 = gsl_vector_get (x, 1);
  const double x2 = gsl_vector_get (x, 2);
  const double x3 = gsl_vector_get (x, 3);
  const double x4 = gsl_vector_get (x, 4);

  const double y0 = Eq1(x0, x1, x2, x3, x4, params);
  const double y1 = Eq2(x0, x1, x2, x3, x4, params);
  const double y2 = Eq3(x0, x1, x2, x3, x4, params);
  const double y3 = Eq4(x0, x1, x2, x3, x4, params);
  const double y4 = Eq5(x0, x1, x2, x3, x4, params);
     
  gsl_vector_set (f, 0, y0);
  gsl_vector_set (f, 1, y1);
  gsl_vector_set (f, 2, y2);
  gsl_vector_set (f, 3, y3);
  gsl_vector_set (f, 4, y4);
  
  return GSL_SUCCESS;
}

double Eq1 (double x1, double x2, double x3, double x4, double x5, void *params)  //equacao para o campo sigma
{
  double gms2 = ((struct rparams *) params)->gms2;
  double b = ((struct rparams *) params)->b;
  double c = ((struct rparams *) params)->c;
  double xs = ((struct rparams *) params)->xs;

  double Vs, Vr, Vw, mub, mue;
  double rho_s=0.0, m8, k2, k; 
  double f;

  //Mapeio:

  Vs=0.9999*mn*sin(x1)*sin(x1);
  Vw=x2*x2;
  Vr=x3;
  mub=x4;
  mue=x5;
  
  //proton:
  m8= mp-xs*Vs; //massa efetiva
  k2=kf2(mub,mue,m8,Vw,Vr,1.0,0.5,params); //qe, I3
  if(k2>0.0){
    k=sqrt(k2);
    rho_s+=dens_s(m8,k);
  }
      
  //Neutron
  m8= mn-xs*Vs;
  k2=kf2(mub,mue,m8,Vw,Vr,0.0,-0.5,params); //qe, I3
  if(k2>0.0){
    k=sqrt(k2);
    rho_s+=dens_s(m8,k);
  }

  //Lambda
  m8= ml-xs*Vs;
  k2=kf2(mub,mue,m8,Vw,Vr,0.0,0.0,params); //qe, I3
  if(k2>0.0){
    k=sqrt(k2);
    rho_s+=dens_s(m8,k);
  }

  //Sigma -
  m8= msm-xs*Vs;
  k2=kf2(mub,mue,m8,Vw,Vr,-1.0,-1.0,params); //qe, I3
  if(k2>0.0){
    k=sqrt(k2);
    rho_s+=dens_s(m8,k);
  }

  //Sigma +
  m8= msp-xs*Vs;
  k2=kf2(mub,mue,m8,Vw,Vr,1.0,1.0,params); //qe, I3
  if(k2>0.0){
    k=sqrt(k2);
    rho_s+=dens_s(m8,k);
  }

  //Sigma 0
  m8= ms0-xs*Vs;
  k2=kf2(mub,mue,m8,Vw,Vr,0.0,0.0,params); //qe, I3
  if(k2>0.0){
    k=sqrt(k2);
    rho_s+=dens_s(m8,k);
  }

  //Chi -
  m8= mcm-xs*Vs;
  k2=kf2(mub,mue,m8,Vw,Vr,-1.0,-0.5,params); //qe, I3
  if(k2>0.0){
    k=sqrt(k2);
    rho_s+=dens_s(m8,k);
  }

  //Chi 0
  m8= mc0-xs*Vs;
  k2=kf2(mub,mue,m8,Vw,Vr,0.0,0.5,params); //qe, I3
  if(k2>0.0){
    k=sqrt(k2);
    rho_s+=dens_s(m8,k);
  }

  f=b*mn*Vs*Vs;
  f=f+c*pow(Vs,3.0);
  f=f-xs*rho_s;
  f=Vs+gms2*f; 

  return f;
}     

double Eq2 (double x1, double x2, double x3, double x4, double x5, void *params)  //equacao para omega
{
  double gmw2 = ((struct rparams *) params)->gmw2;
  double xw = ((struct rparams *) params)->xw;
  double xs = ((struct rparams *) params)->xs;

  double Vs, Vr, Vw, mub, mue;
  double n=0.0;
  double f,m8, k2,k; 

  //Mapeio:

  Vs=0.9999*mn*sin(x1)*sin(x1);
  Vw=x2*x2;
  Vr=x3;
  mub=x4;
  mue=x5;

  //proton:
  m8= mp-xs*Vs; //massa efetiva
  k2=kf2(mub,mue,m8,Vw,Vr,1.0,0.5,params); //qe, I3
  if(k2>0.0){
    k=sqrt(k2);
    n+=dens(k);
  }

  //Neutron
  m8= mn-xs*Vs;
  k2=kf2(mub,mue,m8,Vw,Vr,0.0,-0.5,params); //qe, I3
  if(k2>0.0){
    k=sqrt(k2);
    n+=dens(k);
  }

  //Lambda
  m8= ml-xs*Vs;
  k2=kf2(mub,mue,m8,Vw,Vr,0.0,0.0,params); //qe, I3
  if(k2>0.0){
    k=sqrt(k2);
    n+=dens(k);
  }

  //Sigma -
  m8= msm-xs*Vs;
  k2=kf2(mub,mue,m8,Vw,Vr,-1.0,-1.0,params); //qe, I3
  if(k2>0.0){
    k=sqrt(k2);
    n+=dens(k);
  }

  //Sigma +
  m8= msp-xs*Vs;
  k2=kf2(mub,mue,m8,Vw,Vr,1.0,1.0,params); //qe, I3
  if(k2>0.0){
    k=sqrt(k2);
    n+=dens(k);
  }

  //Sigma 0
  m8= ms0-xs*Vs;
  k2=kf2(mub,mue,m8,Vw,Vr,0.0,0.0,params); //qe, I3
  if(k2>0.0){
    k=sqrt(k2);
    n+=dens(k);
  }

  //Chi -
  m8= mcm-xs*Vs;
  k2=kf2(mub,mue,m8,Vw,Vr,-1.0,-0.5,params); //qe, I3
  if(k2>0.0){
    k=sqrt(k2);
    n+=dens(k);
  }

  //Chi 0
  m8= mc0-xs*Vs;
  k2=kf2(mub,mue,m8,Vw,Vr,0.0,0.5,params); //qe, I3
  if(k2>0.0){
    k=sqrt(k2);
    n+=dens(k);
  }

  f=(xw*gmw2)*n;
  f=Vw-f;

  return f;
}

double Eq3 (double x1, double x2, double x3, double x4, double x5, void *params)  //equacao para o rho
{
  double gmr2 = ((struct rparams *) params)->gmr2;
  double xs = ((struct rparams *) params)->xs;
  double xr = ((struct rparams *) params)->xr;

  double Vs, Vr, Vw, mub, mue;
  double n=0.0;
  double f,m8,k2,k; 

  //Mapeio:

  Vs=0.9999*mn*sin(x1)*sin(x1);
  Vw=x2*x2;
  Vr=x3;
  mub=x4;
  mue=x5;

  //proton:
  m8= mp-xs*Vs;
  k2=kf2(mub,mue,m8,Vw,Vr,1.0,0.5,params); //qe, I3
  if(k2>0.0){
    k=sqrt(k2);
    n+=0.5*dens(k);
  }

  //Neutron
  m8= mn-xs*Vs;
  k2=kf2(mub,mue,m8,Vw,Vr,0.0,-0.5,params); //qe, I3
  if(k2>0.0){
    k=sqrt(k2);
    n+=-0.5*dens(k);
  }

  //Sigma -
  m8= msm-xs*Vs;
  k2=kf2(mub,mue,m8,Vw,Vr,-1.0,-1.0,params); //qe, I3
  if(k2>0.0){
    k=sqrt(k2);
    n+=-1.0*dens(k);
  }

  //Sigma +
  m8= msp-xs*Vs;
  k2=kf2(mub,mue,m8,Vw,Vr,1.0,1.0,params); //qe, I3
  if(k2>0.0){
    k=sqrt(k2);
    n+=dens(k);
  }

  //Chi -
  m8= mcm-xs*Vs;
  k2=kf2(mub,mue,m8,Vw,Vr,-1.0,-0.5,params); //qe, I3
  if(k2>0.0){
    k=sqrt(k2);
    n+=-0.5*dens(k);
  }

  //Chi 0
  m8= mc0-xs*Vs;
  k2=kf2(mub,mue,m8,Vw,Vr,0.0,0.5,params); //qe, I3
  if(k2>0.0){
    k=sqrt(k2);
    n+=0.5*dens(k);
  }

  f=n*xr*gmr2;
  f=f-Vr;

  return f;
}
 
double Eq4 (double x1, double x2, double x3, double x4, double x5, void *params)  //neutralidade da carga eletrica (sum(ni*Qi)=0)
{
  double xs = ((struct rparams *) params)->xs;

  double Vs, Vr, Vw, mub, mue;
  double n=0.0, ke=0.0, kmu=0.0;
  double m8, k2, k; 

  //Mapeio:

  Vs=0.9999*mn*sin(x1)*sin(x1);
  Vw=x2*x2;
  Vr=x3;
  mub=x4;
  mue=x5;

  //barions

  //proton:
  m8= mp-xs*Vs;
  k2=kf2(mub,mue,m8,Vw,Vr,1.0,0.5,params); //qe, I3
  if(k2>0.0){
    k=sqrt(k2);
    n+=dens(k);
  }

  //Sigma -
  m8= msm-xs*Vs;
  k2=kf2(mub,mue,m8,Vw,Vr,-1.0,-1.0,params); //qe, I3
  if(k2>0.0){
    k=sqrt(k2);
    n+=-1.0*dens(k);
  }

  //Sigma +
  m8= msp-xs*Vs;
  k2=kf2(mub,mue,m8,Vw,Vr,1.0,1.0,params); //qe, I3
  if(k2>0.0){
    k=sqrt(k2);
    n+=dens(k);
  }

  //Chi -
  m8= mcm-xs*Vs;
  k2=kf2(mub,mue,m8,Vw,Vr,-1.0,-0.5,params); //qe, I3
  if(k2>0.0){
    k=sqrt(k2);
    n+=-1.0*dens(k);
  }

  //leptons

  ke=mue*mue-me*me;
  if(ke>0.0){
    ke=sqrt(ke);
    n+=-1.0*dens(ke);
  }

  kmu=mue*mue-mmu*mmu;
  if(kmu>0.0){
    kmu=sqrt(kmu);
    n+=-1.0*dens(kmu);
  }

  return n;
}     

double Eq5 (double x1, double x2, double x3, double x4, double x5,void *params)  //densidade (de número) barionica  (MeV^3)
{
  double xs = ((struct rparams *) params)->xs;
  double rho = ((struct rparams *) params)->rho;

  double Vs, Vr, Vw, mub, mue;
  double n=0.0;
  double f,m8,k2,k; 

  //Mapeio:

  Vs=0.9999*mn*sin(x1)*sin(x1);
  Vw=x2*x2;
  Vr=x3;
  mub=x4;
  mue=x5;

  //proton:
  m8= mp-xs*Vs;
  k2=kf2(mub,mue,m8,Vw,Vr,1.0,0.5,params); //qe, I3
  if(k2>0.0){
    k=sqrt(k2);
    n+=dens(k);
  }

  //Neutron
  m8= mn-xs*Vs;
  k2=kf2(mub,mue,m8,Vw,Vr,0.0,-0.5,params); //qe, I3
  if(k2>0.0){
    k=sqrt(k2);
    n+=dens(k);
  }

  //Lambda
  m8= ml-xs*Vs;
  k2=kf2(mub,mue,m8,Vw,Vr,0.0,0.0,params); //qe, I3
  if(k2>0.0){
    k=sqrt(k2);
    n+=dens(k);
  }

  //Sigma -
  m8= msm-xs*Vs;
  k2=kf2(mub,mue,m8,Vw,Vr,-1.0,-1.0,params); //qe, I3
  if(k2>0.0){
    k=sqrt(k2);
    n+=dens(k);
  }

  //Sigma +
  m8= msp-xs*Vs;
  k2=kf2(mub,mue,m8,Vw,Vr,1.0,1.0,params); //qe, I3
  if(k2>0.0){
    k=sqrt(k2);
    n+=dens(k);
  }

  //Sigma 0
  m8= ms0-xs*Vs;
  k2=kf2(mub,mue,m8,Vw,Vr,0.0,0.0,params); //qe, I3
  if(k2>0.0){
    k=sqrt(k2);
    n+=dens(k);
  }

  //Chi -
  m8= mcm-xs*Vs;
  k2=kf2(mub,mue,m8,Vw,Vr,-1.0,-0.5,params); //qe, I3
  if(k2>0.0){
    k=sqrt(k2);
    n+=dens(k);
  }

  //Chi 0
  m8= mc0-xs*Vs;
  k2=kf2(mub,mue,m8,Vw,Vr,0.0,0.5,params); //qe, I3
  if(k2>0.0){
    k=sqrt(k2);
    n+=dens(k);
  }

  f=n-rho;
  return f;
}     

double kf2(double mub, double mue, double m8, double Vw, double Vr, double qe, double I3, void *params)
{
  double xr = ((struct rparams *) params)->xr;
  double xw = ((struct rparams *) params)->xw;

  double k2, mu;

  mu=mub-qe*mue;
  k2=mu-xw*Vw-I3*xr*Vr; //momento de Fermi ao quadrado
  k2=k2*k2;
  k2=k2-m8*m8;
  //  printf("k2=%9.9lf\n",k2);

  return k2;
}

double dens_s(double m8, double k)
{
  double Pi=M_PI;
  double ns;

  ns=k+sqrt(k*k+m8*m8);
  ns=m8*m8*log(ns/m8);
  ns=k*sqrt(k*k+m8*m8)-ns;
  ns=m8/(2.0*Pi*Pi)*ns;

  return ns;
}

double dens(double k)
{
  double Pi=M_PI;
  double n;

  n=pow(k,3.0);
  n=n/(3.0*Pi*Pi);

  return n;
}