/* GROUND STATE COMPOSITION AND STRUCTURE OF THE OUTER CRUST OF MAGNETARS      */

/* THE CODE COMPUTES THE FOLLOWING QUANTITIES:                                 */
/* - EQUILIBRIUM NUCLEI CONSTITUTING EACH CRUSTAL LAYER,                       */
/* - TRANSITION PRESSURES,                                                     */
/* - DENSITIES OF ADJACENT LAYERS,                                             */
/* - RELATIVE NUCLEAR ABUNDANCES,                                              */
/* - DEPTHS RELATIVE TO NEUTRON-DRIP TRANSITION.                               */

/* REMARKS:                                                                    */
/* -ONLY LANDAU-RABI QUANTIZATION OF ELECTRON MOTION INCLUDED,                 */
/* -ELECTRONS ARE ASSUMED TO BE ULTRARELATIVISTIC,                             */
/* -ELECTRONS ARE SUPPOSED TO BE IN THE LOWEST LANDAU-RABI LEVEL.              */

/* execution: smagcrust <atomic_mass_data_file>                                */
/* options:                                                                    */
/* -w Wigner-Seitz approximation                                               */
/* The format of <atomic_mass_data_file> should be the following:              */
/*     atomic number / mass number / mass excess in MeV                        */
/* help: smagcrust -h                                                          */

/* code written by Nicolas Chamel                                              */
/* version 22.05.2020                                                          */


/* C headers */
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>

/* --- fundamental constants from NIST CODATA 2018 --- */
/* https://physics.nist.gov/cuu/Constants/             */

/* atomic mass unit in MeV/c^2 */
const double ATOMIC_MASS=931.49410242 ;

/* electron mass in MeV/c^2 */
const double ELECTRON_MASS=0.51099895 ;

/* neutron mass in MeV/c^2 */
const double NEUTRON_MASS=939.56542052 ;

/* neutron magnetic moment to nuclear magneton ratio */
const double MUN=-1.91304273 ;

/* dielectric permittivity of vacuum in units of 10^12 F/m */
const double EPS0=8.8541878128 ;

/* elementary electric charge in units of 10^-19  C */
const double E1=1.602176634 ;

/* hbar c in units MeV.fm */
const double HBARC=197.3269804 ;

/* subroutine to empty the buffer */
void empty_buffer(void)
{
  int c;
  while((c=getchar()) != EOF && c != '\n');
}


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

    /* useful combinations of constants */
    const double LAMBDAE=HBARC/ELECTRON_MASS ; /* electron Compton wave length */
    const double E2=E1*100.0/(4.0*M_PI*EPS0) ; 
    const double ALPHA=E2/HBARC ; /* fine structure constant */
    const double BREL=sqrt(ELECTRON_MASS/ALPHA/(LAMBDAE*LAMBDAE*LAMBDAE)*E1*1e33) ; /* characteristic magnetic field in G */
    
    double PRESS_CONST ;
    double DENS_CONST ;
    
    /* body-centered cubic lattice constant (1) or Wigner-Seitz approximation (0) */
    char BCC=1 ;
    
	/* variables */
    short NDATA ;
    short *Z, *A ;
    double *EXCESS ;
    double MASS1,MASS2,NUCLEON_MASS0 ;
    double xe,sqtxe,mue12,gam12,F12,F12bar ;
    double P120,n1max0,n2min0,P122,n1max2,n2min2,P12,n1max,n2min ;
    double u, gameq, gam0, gam2, gammae ;
    short Z1, A1, Z2, A2 ;
    double Pmin, g12, Peq, gs, gdrip ;
	short i, imin ;
    double y1,y2 ;
    FILE *data_file ;
    int read ;
    char nuc=1 ;
    char sol=0 ;
    double xi=0.0 ;
    double sum=0.0 ;
    double BSTAR ;
    char reliability=0 ;
    int input ;
    double Cion ;
    double depth ;

    
	/* check the number of arguments */
    if((argc<2)||(argc>3))
	{
		fprintf(stderr,"The number of arguments is invalid. Use option -h for help.\a\n" );
		return EXIT_FAILURE ;
	}
    
	/* print help */
	if(strcmp(argv[1],"-h") == 0 || strcmp(argv[1],"-help") == 0) 
	{
		printf("Usage:\n");
		printf("%s <filename> options\n",argv[0]);
		printf("where <filename> is the name of a file containing the atomic mass table.\n");
        printf("Options:\n");
        printf("-w Wigner-Seitz approximation\n") ;
        printf("The format of the mass table should be the following:\n") ;
        printf("atomic number / mass number / mass excess in MeV.\n") ;
		return EXIT_SUCCESS;
	}
    
    /* read arguments */
    if (argc==3)
    {
        if(strcmp(argv[2],"-w") == 0)
            BCC=0 ;
    }
    
    /* open file */
    data_file=fopen(argv[1], "r") ;
    if(data_file==NULL)
    {
        fprintf(stderr,"The file does not exist!\a\n") ;
        return EXIT_FAILURE ;
    }
    
    /* determine number of nuclei and check for errors */
    i=0 ;
    read=fscanf(data_file, "%hd %hd %lf",&Z1,&A1,&MASS1) ;
    while(read!=EOF)
    {
        if((read!=3)||(Z1<1)||(A1<1))
        {
            fprintf(stderr,"Error when reading file at line %d!\a\n",i) ;
            exit(EXIT_FAILURE) ;
        }
        read=fscanf(data_file, "%hd %hd %lf",&Z1,&A1,&MASS1) ;
        i++ ;
    }
    NDATA=i ;
    rewind(data_file) ;
    
    /* allocate memory for the arrays */
    Z=(short *)malloc(NDATA*sizeof(short)) ;
    A=(short *)malloc(NDATA*sizeof(short)) ;
    EXCESS=(double *)malloc(NDATA*sizeof(double)) ;

    /* read mass excesses */
    for (i=0; i<NDATA; i++)
    {
        if (fscanf(data_file, "%hd %hd %lf",&Z[i],&A[i],&EXCESS[i])!=3)
        {
            fprintf(stderr,"Error when reading mass excess at line %d!\a\n",i) ;
            exit(EXIT_FAILURE) ;
        }
    }
    fclose(data_file) ;
    
    printf("Magnetic field strength in units of critical field %4.3le G?\n",BREL) ;
    input=scanf("%lf",&BSTAR) ;
    /* check the parameter */
    while((input!=1)||(BSTAR<0))
    {
        fprintf(stderr,"Invalid argument!\a\n") ;
        empty_buffer() ;
        input=scanf("%lf",&BSTAR) ;
    }
    
    /* useful constants */
    DENS_CONST=BSTAR/(2.0*M_PI*M_PI*LAMBDAE*LAMBDAE*LAMBDAE) ;
    PRESS_CONST=0.5*DENS_CONST*ELECTRON_MASS ;

    if (BCC==0)
        /* Wigner-Seitz approximation                      */
        Cion=-0.9*cbrt(4.0*M_PI/3.0) ;
    else
        /* body-centered cubic lattice constant            */
        /* from Baiko et al., Phys. Rev. E64,057402 (2001) */
        Cion=-0.895929255682*cbrt(4.0*M_PI/3.0) ;
    
    printf("\f") ;
	printf("----------------------------------------------------------------------------------\n") ;	
	printf("GROUND-STATE COMPOSITION OF THE OUTER CRUST OF A COLD NONACCRETING MAGNETAR \n") ;
    printf("WITH A MAGNETIC FIELD  B=%4.3le G\n", BSTAR*BREL) ;
    printf("STRONGLY QUANTIZING REGIME FOR ELECTRONS IS ASSUMED.\n") ;
	printf("----------------------------------------------------------------------------------\n") ;	

    printf("Atomic mass table: %s\n",argv[1]) ;
    printf("%hd nuclei included.\n",NDATA) ;
    
    if (BCC==0)
        printf("Wigner-Seitz approximation adopted.\n") ;
    else
        printf("Body-centered cubic lattice assumed.\n") ;
    
    
    printf("\n") ;
    
    /* composition of the surface layers */
    Z1=26 ;
    A1=56 ;
    y1=(double)Z1/(double)A1 ;
    MASS1=-60.607+(double)A1*ATOMIC_MASS+0.0000144381*pow(Z1,2.39)+1.55468e-12*pow(Z1,5.35) ; /* subtract electron binding energy */
    NUCLEON_MASS0=MASS1/A1 ; /* mean nucleon mass in iron-56 */
    Pmin=1e-2 ;
    xe=pow(Z1*Z1*pow(-Cion*ALPHA,3)*BSTAR/(2.0*M_PI*M_PI),1.0/5.0) ;
    gameq=sqrt(1.0+xe*xe) ;
    gs=MASS1/(double)A1+y1*((gameq-1.0)*ELECTRON_MASS+4.0/3.0*Cion*E2*cbrt(Z1*Z1*DENS_CONST*xe)) ;
    g12=gs ;
    Peq=0.0 ;
    gdrip=NEUTRON_MASS ;
    
    /* composition of the layers below neutron drip */
    printf("Z1 A1 Z2 A2 \t xe \t n1max[fm-3] \t n2min[fm-3] \t P12[MeV fm-3] \t mue12[MeV] \t g12[MeV]\t abundance \t depth\n") ;
    while((g12<gdrip)&&(gameq<sqrt(1.0+2.0*BSTAR)))
    {
        for (i=0; i<NDATA; i++)
        {
            sol=0 ;
            
            Z2=Z[i] ;
            A2=A[i] ;
            y2=(double)Z2/(double)A2 ;
            MASS2=EXCESS[i]+(double)A2*ATOMIC_MASS+0.0000144381*pow(Z2,2.39)+1.55468e-12*pow(Z2,5.35) ;  /* subtract electron binding energy */
            
            /* rule out transitions such that Z1/A1=Z2/A2  */
            /* see Appendix A of Phys.Rev.C94,065802(2016) */
            if (y1!=y2)
            {
                /* analytical solution for the transition between (A1,Z1) and (A2,Z2) */
                mue12=(MASS2/(double)A2-MASS1/(double)A1)/(y1-y2)+ELECTRON_MASS ;
                gam12=mue12/ELECTRON_MASS ;
                F12=((4.0*y1/3.0-y2/3.0)*cbrt(Z1*Z1)-y2*cbrt(Z2*Z2))/(y1-y2) ;
                F12bar=F12*Cion/3.0*ALPHA*cbrt(BSTAR/(2.0*M_PI*M_PI)) ;
            
                if ((F12bar>0.0)&&(gam12>0.0))
                {
                    sol=1 ;
                    u=0.5*gam12/pow(sqrt(F12bar),3.0) ;
                    gammae=8.0*pow(sqrt(F12bar),3.0)*pow(sinh(asinh(u)/3.0),3.0) ;
                }
                if (F12bar<0.0)
                {
                
                    u=0.5*gam12/pow(sqrt(fabs(F12bar)),3.0) ;
                    if (u>=1.0)
                    {
                        gammae=8.0*pow(sqrt(fabs(F12bar)),3.0)*pow(cosh(acosh(u)/3.0),3.0) ;
                        sol=1 ;
                    }
                    if ((u<1.0)&&(u>=0.0))
                    {
                        gammae=8.0*pow(sqrt(fabs(F12bar)),3.0)*pow(cos(acos(u)/3.0),3.0) ;
                        sol=1 ;
                    }
                    if ((u>=-1.0)&&(u<0.0))
                    {
                        gam0=8.0*pow(sqrt(fabs(F12bar)),3.0)*pow(cos(acos(u)/3.0),3.0) ;
                        gam2=8.0*pow(sqrt(fabs(F12bar)),3.0)*pow(cos(acos(u)/3.0+4.0*M_PI/3.0),3.0) ;
                        sol=2 ;
                    }
                }
                /* handle special cases for which two real positive mathematical solutions exist */
                if (sol==2)
                {
                    sol=0 ;
                    xe=sqrt(gam0*gam0-1.0) ;
                    n1max0=xe/y1 ;
                    n2min0=xe/y2*(1.0+Cion*ALPHA/3.0*(cbrt(Z1*Z1)-cbrt(Z2*Z2))*cbrt(BSTAR/(2.0*M_PI*M_PI))*gam0/pow(xe,5.0/3.0)) ;
                    xe=sqrt(gam2*gam2-1.0) ;
                    n1max2=xe/y1 ;
                    n2min2=xe/y2*(1.0+Cion*ALPHA/3.0*(cbrt(Z1*Z1)-cbrt(Z2*Z2))*cbrt(BSTAR/(2.0*M_PI*M_PI))*gam2/pow(xe,5.0/3.0)) ;
                
                    /* check mechanical stability of each solution separately */
                    if ((n2min0>n1max0)&&(n1max0>0.0)&&(gam0>=1.0))
                    {
                        sol=1 ;
                        /*printf("k=0 solution\n") ;*/
                        gammae=gam0 ;
                    }
                    if ((n2min2>n1max2)&&(n1max2>0.0)&&(gam2>=1.0))
                    {
                        /*printf("k=2 solution\n") ;*/
                        gammae=gam2 ;
                        sol++ ;
                    }
                    /* if the two solutions are still admissible, select the one with the lowest transition pressure */
                    if (sol==2)
                    {
                        /*printf("Two real solutions found!\n") ;*/
                        xe=sqrt(gam0*gam0-1.0) ;
                        P120=PRESS_CONST*(xe*gam0-log(xe+gam0))+Cion*E2/3.0*cbrt(pow(DENS_CONST*xe,4)*Z1*Z1) ;
                        xe=sqrt(gam2*gam2-1.0) ;
                        P122=PRESS_CONST*(xe*gam2-log(xe+gam2))+Cion*E2/3.0*cbrt(pow(DENS_CONST*xe,4)*Z1*Z1) ;
                        if (P120<P122)
                            gammae=gam0 ;
                        else
                            gammae=gam2 ;
                        sol=1 ;
                    }
                }
                
                /* if a solution to the transition between (A1,Z1) and (A2,Z2) exists */
                if (sol==1)
                {
                    xe=sqrt(gammae*gammae-1.0) ;
                    n1max=DENS_CONST*xe/y1 ;
                    n2min=DENS_CONST*xe/y2*(1.0+Cion*ALPHA/3.0*(cbrt(Z1*Z1)-cbrt(Z2*Z2))*cbrt(BSTAR/(2.0*M_PI*M_PI))*gammae/pow(xe,5.0/3.0)) ;
                
                    /* check mechanical stability */
                    if ((n2min>n1max)&&(n1max>0.0)&&(gammae>=1.0))
                    {
                        P12=PRESS_CONST*(xe*gammae-log(xe+gammae))+Cion*E2/3.0*cbrt(pow(DENS_CONST*xe,4)*Z1*Z1) ;
                    
                        /* find lowest transition pressure              */
                        /* ensure the pressure increases with density   */
                        if ((P12<Pmin)&&(P12>=Peq))
                        {
                            Pmin=P12 ;
                            gameq=gammae ;
                            imin=i ;
                        }
                    }
                }
            }
        }
        Z2=Z[imin] ;
        A2=A[imin] ;
        y2=(double)Z2/(double)A2 ;
        MASS2=EXCESS[imin]+(double)A2*ATOMIC_MASS+0.0000144381*pow(Z2,2.39)+1.55468e-12*pow(Z2,5.35) ; 
        mue12=(MASS2/(double)A2-MASS1/(double)A1)/(y1-y2)+ELECTRON_MASS ;
        gam12=mue12/ELECTRON_MASS ;
        xe=sqrt(gameq*gameq-1.0) ;
        n1max=DENS_CONST*xe/y1 ;
        n2min=DENS_CONST*xe/y2*(1.0+Cion*ALPHA/3.0*(cbrt(Z1*Z1)-cbrt(Z2*Z2))*cbrt(BSTAR/(2.0*M_PI*M_PI))*gameq/pow(xe,5.0/3.0)) ;
        P12=PRESS_CONST*(xe*gameq-log(xe+gameq))+Cion*E2/3.0*cbrt(pow(DENS_CONST*xe,4)*Z1*Z1) ;
        g12=MASS1/(double)A1+y1*((gameq-1.0)*ELECTRON_MASS+4.0/3.0*Cion*E2*cbrt(Z1*Z1*DENS_CONST*xe)) ;
        if ((g12<gdrip)&&(gameq<sqrt(1.0+2.0*BSTAR)))
        {
            nuc++ ;
            sum+=xi ;
            xi=P12-sum ;
            Peq=P12 ;
            depth=(pow(g12/gs,2)-1.0)/(pow(NEUTRON_MASS/gs,2)-1.0) ;
            /* check reliability of the solution */
            if (gam12<=0.0)
            {
                reliability=1 ;
                printf("%hd %hd %hd %hd \t %4.2f \t %4.3le \t %4.3le \t %4.3le \t %4.2f* \t %6.3f \t %3.2le \t %5.3f\n",Z1,A1,Z2,A2,xe,n1max,n2min,P12,mue12,g12,xi,depth) ;
            }
            /* normal case */
            if (gam12>0.0)
                printf("%hd %hd %hd %hd \t %4.2f \t %4.3le \t %4.3le \t %4.3le \t %4.2f \t\t %6.3f \t %3.2le \t %5.3f\n",Z1,A1,Z2,A2,xe,n1max,n2min,P12,mue12,g12,xi,depth) ;

            Z1=Z2 ;
            A1=A2 ;
            y1=y2 ;
            MASS1=MASS2 ;
            Pmin=1e-2 ;
        }
    }
    /* check strongly quantizing regime */
    if (gameq>=sqrt(1.0+2.0*BSTAR))
    {
        printf("Calculation stopped: lowest Landau-Rabi level completely filled.\a\n") ;
        printf("Abundances cannot be normalized!\n") ;
        
    }
    if (gameq<sqrt(1.0+2.0*BSTAR))
    {
        /* properties of the neutron-drip transition */
        mue12=(-MASS1+(double)A1*gdrip)/(double)Z1+ELECTRON_MASS ; /* always positive ?? */
        gam12=mue12/ELECTRON_MASS ;
        F12bar=4.0/3.0*Cion/3.0*ALPHA*cbrt(Z1*Z1*BSTAR/(2.0*M_PI*M_PI)) ;  /* always negative */
        u=0.5*gam12/pow(sqrt(fabs(F12bar)),3.0) ;
        if (u>=1.0)
            gameq=8.0*pow(sqrt(fabs(F12bar)),3.0)*pow(cosh(acosh(u)/3.0),3.0) ;
        if ((u<1.0)&&(u>=0.0))
            gameq=8.0*pow(sqrt(fabs(F12bar)),3.0)*pow(cos(acos(u)/3.0),3.0) ;
        xe=sqrt(gameq*gameq-1.0) ;
        n1max=DENS_CONST*xe/y1 ;
        P12=PRESS_CONST*(xe*gameq-log(xe+gameq))+Cion*E2/3.0*cbrt(pow(DENS_CONST*xe,4)*Z1*Z1) ;
        sum+=xi ;
        xi=P12-sum ;
        depth=1.0 ;
        printf("%hd %hd -  - \t %4.2f \t %4.3le \t - \t\t %4.3le \t %4.2f \t\t %6.3f \t %3.2le \t %5.3f\n",Z1,A1,xe,n1max,P12,mue12,gdrip,xi,depth) ;
        printf("The outer crust is stratified into %d layers.\n",nuc) ;
        printf("Warning: nuclear abundances are not normalized. Their sum is not 1 but %4.3le.\n", sum+xi) ;
        printf("Depths are relative to the neutron-drip point.\n") ;
    }

    if(reliability==1)
        printf("* Warning: this solution is not expected to be accurate.\a\n") ;
    
    free(Z) ;
    Z=NULL ;
    free(A) ;
    A=NULL ;
    free(EXCESS) ;
    EXCESS=NULL ; 
    
	return EXIT_SUCCESS ;
}