program nstar
      include 'implno.dek'
      include 'const.dek'

c..
c..this program generates non-rotating neutron star models
c..

c..declare
      character*10     eosfile
      character*80     string
      integer          jj,nstep
      double precision denc,mass,den,dlo,dhi,dstp,presc,
     1                 con1,third
      parameter        (con1  = 4.0d0/3.0d0 * pi,
     1                  third = 1.0d0/3.0d0)


c..declare the integration technology
      external          derv,rkqc
      integer           xdim,ydim,kount,nok,nbad,iprint,i
      parameter         (xdim=600, ydim=4)
      double precision  xrk(xdim),yrk(ydim,xdim),bc(ydim),dydx,stptry,
     1                  stpmin,tol,odescal,start,sstop,dxsav



c..stop integrating when the density or pressure drop below these numbers
      double precision stopden,stop_pres
      common  /sbyiol/ stopden,stop_pres



c..communicate if general relativistic corrections are included
      logical          genrel
      common /grcm1/   genrel



c..for communicating the nuclear eos
      integer          np,npmax
      parameter        (npmax=500)
      double precision eray(npmax),pray(npmax),xnray(npmax)
      common /nccf1/   eray,pray,xnray,np




c..formats
01    format(1x,i4,1p8e12.4)
02    format(1x,1p5e12.4,i5)
03    format(a) 
04    format(1x,a,t10,a,t22,a,t34,a,t46,a,t58,a,t70,a,t82,a)
05    format(1x,t5,a,t18,a,t30,a,t42,a,t54,a,t66,a,t78,a)




c..set a flag if general relativistic corrections are to be used
      genrel = .true.
c      genrel = .false.


c..open and rad in one of the nuclear eos tables
c      eosfile = 'eosC'
      eosfile = '3.d'
      call nucread(eosfile,xnray,pray,eray,np,npmax)



c..get the central density and find the central pressure
c      write(6,*) 'give the central density =>'
c      read(5,*) denc
c      call nuceos(eray,pray,np,denc,presc)


     
c..get the output file name
      write(6,*) 'give output file name =>'
      read(5,03) string

      open(unit=3, file=string, status='unknown')



c..for comparing the central density to total mass 
c..nstep evenly spaced points in log space

      dlo   = 14.0d0
      dhi   = 16.0d0
      nstep = 41
      dstp = 0.0d0
      if (nstep .ne. 1) dstp = (dhi - dlo)/(nstep - 1)


c.sweep over the central density

      write(3,05) 'density','massg','massp','mass0','radius'

      do jj = 1,nstep
       denc = dlo + float(jj-1) * dstp
       denc = 10.0d0**denc
       call nuceos(eray,pray,np,denc,presc)



c..integrate in radius (in cm) from start to sstop
      start   = 1.0e2
      stptry  = start
      stpmin  = 1.0d0
      sstop   = rsol


c..set the initial conditions: 
c..bc(1) is the gravitational mass
c..bc(2) is the proper mass 
c..bc(3) is the amu or rest mass
c..bc(4) is the central pressure
c..these are expansions about zero radius

      bc(1)   = con1 * start**3 * denc
      bc(2)   = bc(1)
      bc(3)   = bc(1)
      bc(4)   = presc - 0.5d0 * con1 * g * start**2 * denc**2



c..get the stopping density and pressure
      stopden   = max(1.0d0, 1.0e-9 * denc)
      call nuceos(eray,pray,np,stopden,stop_pres)

      tol     = 1.0e-6
      dxsav   = 0.0d0
      odescal = 1.0
      iprint  = 0


c..and integrate the odes in routine derv

      call fermint(start,stptry,stpmin,sstop,bc,
     1            tol,dxsav,xdim,
     2            xrk,yrk,xdim,ydim,xdim,ydim,
     3            nok,nbad,kount,odescal,iprint,
     4            derv,rkqc)




c..write out the profile
c      write(3,04) '  i','radius','mass_grav','mass_prop','mass_rest',
c     1            'pressure','density'
c      do i=1,kount
c       call invert_nuceos(eray,pray,np,den,yrk(4,i))
c       write(3,01) i,xrk(i)/1.0e5,yrk(1,i)/msol,yrk(2,i)/msol,
c     1             yrk(3,i)/msol,yrk(4,i),den
c      enddo



      write(3,02) denc,yrk(1,kount)/msol,yrk(2,kount)/msol,
     1            yrk(3,kount)/msol,xrk(kount)/1.0e5,kount

      write(6,02) denc,yrk(1,kount)/msol,yrk(2,kount)/msol,
     1            yrk(3,kount)/msol,xrk(kount)/1.0e5,kount


c..back for another central density
      enddo

      close(unit=3)
      stop 'normal termination'
      end









      subroutine derv(x,y,dydx)
      include 'implno.dek'
      include 'const.dek'

c..this routine sets up the continuity and hydrostatic equilibrium ode's. 
c..x is the radial coordinate, y(1) is the gravitational mass,
c..y(2) is the proper mass, y(3) is the rest mass, y(4) is the pressure

c..declare the pass
      double precision x,y(*),dydx(*)


c..communicate if general relativistic corrections are included
      logical          genrel
      common /grcm1/   genrel



c..for communicating the nuclear eos
      integer          np,npmax
      parameter        (npmax=500)
      double precision eray(npmax),pray(npmax),xnray(npmax)
      common /nccf1/   eray,pray,xnray,np


c..local variables
      integer          i,niter
      double precision massg,massp,mass0,den,pres,cor,xn,con1,c2
      parameter        (con1   = 4.0d0 * pi,
     1                  c2     = clight*clight)


c..map the input vector
      massg = y(1)
      massp = y(2)
      mass0 = y(3) 
      pres  = y(4)


c..a cold nuclear eos
      call invert_nuceos(eray,pray,np,den,pres)
      call invert_nuceos(xnray,pray,np,xn,pres)


c..the gravitational mass
      dydx(1) = con1 * x**2 * den 



c..the proper mass
      if (genrel) then
       cor = 1.0d0/sqrt(1.0d0 - (2.0d0*g*massg)/(x*c2)) 
      else
       cor = 1.0d0
      end if

      dydx(2) = dydx(1) * cor


c..the rest mass
      dydx(3) = con1 * x**2 * xn/avo * cor


c..here is d(press)/dr
      if (genrel) then
       cor = (1.0d0 + pres/(den*c2)) * 
     1       (1.0d0 + (con1*pres*x**3)/(massg*c2)) /
     2       (1.0d0 - (2.0d0*g*massg)/(x*c2))
      else
       cor = 1.0d0
      end if

      dydx(4) = -g * massg/x**2 * den * cor 

      return
      end







      subroutine nucread(file,xnray,pray,eray,np,npmax)
      include 'implno.dek'

c..opens one of the nuclear eos files
c..and returns the number density and pressure arrays

c..declare the pass
      character*(*)    file
      integer          np,npmax
      double precision xnray(npmax),pray(npmax),eray(npmax) 


c..local variables
      integer          i
      double precision xx

c..open the file and read it
      open(unit=22,file=file,status='old')
      read(22,*) np
      if (np .gt. npmax) stop 'np > npmax in routine teos'
      do i=1,np
       read(22,*) eray(i),pray(i),xx,xnray(i)
      enddo
      close(unit=22)


c..work in log space
      do i=1,np
       eray(i)  = log10(eray(i))
       pray(i)  = log10(pray(i))
       xnray(i) = log10(xnray(i))
      enddo 


      return
      end




      subroutine nuceos(xnray,pray,np,xn,pres)
      include 'implno.dek'

c..given the number density, this routine returns the pressure
c..by interpolating one of the nuclear eos files

c..declare the pass
      integer          np
      double precision xnray(np),pray(np),xn,pres


c..local variables
c..order = 2 is a linear interpolation
c..order = 3 is a quadratic interpolant
c..order = 4 is a cubic interpolant
c..stick with a linear interpolant in case there are sharpish phase transitions

      integer          jat,order
      parameter        (order = 2)
      double precision xnl,pl,dy


      xnl = log10(xn)
      call locate(xnray,np,xnl,jat)

c..comment these out if you don't mind linear extrapolation off the table
c      if (jat .eq. 0) stop 'xnl off lower table bound in nuceos'
c      if (jat .eq. np) stop 'xnl off upper table bound in nuceos'

      jat = max(1, min(jat - order/2 + 1, np - order + 1))
      call polint(xnray(jat),pray(jat),order,xnl,pl,dy)
      pres = 10.0d0**pl

      return
      end
 




      subroutine invert_nuceos(xnray,pray,np,xn,pres)
      include 'implno.dek'

c..given the pressure, this routine returns the number density
c..by interpolating one of the nuclear eos files

c..declare the pass
      integer          np
      double precision xnray(np),pray(np),xn,pres


c..local variables
c..order = 2 is a linear interpolation
c..order = 3 is a quadratic interpolant
c..order = 4 is a cubic interpolant
c..stick with a linear interpolant in case there are sharpish phase transitions

      integer          jat,order
      parameter        (order = 2)
      double precision xnl,pl,dy

      pl = log10(pres)
      call locate(pray,np,pl,jat)

c..comment these out if you don't mind linear extrapolation off the table
c      if (jat .eq. 0) stop 'xnl off lower table bound in invert_nuceos'
c      if (jat .eq. np) stop 'xnl off upper table bound in invert_nuceos'

      jat = max(1, min(jat - order/2 + 1, np - order + 1))
      call polint(pray(jat),xnray(jat),order,pl,xnl,dy)
      xn = 10.0d0**xnl

      return
      end








      subroutine locate(xx,n,x,j) 
      include 'implno.dek' 

c..given an array xx of length n, and a value of x, this routine returns 
c..a value j such that x is between xx(j) and xx(j+1). the array xx must be  
c..monotonic. j=0 or j=n indicates that x is out of range. bisection is used  
c..to find the entry 

c..declare 
      integer           n,j,jl,ju,jm 
      double precision  xx(n),x 


c..initialize 
      jl = 0 
      ju = n+1 


c..compute a midpoint, and replace either the upper or lower limit 
 10    if (ju-jl .gt. 1) then 
       jm = (ju+jl)/2 
       if ( (xx(n) .ge. xx(1)) .eqv. (x .ge. xx(jm)) ) then 
        jl = jm 
       else 
        ju = jm 
       end if 
       goto 10 
      end if 
      if (x .eq. xx(1))then 
        j = 1 
      else if(x .eq. xx(n))then 
        j = n - 1 
      else 
        j = jl 
      end if 
      return 
      end 








      subroutine polint(xa,ya,n,x,y,dy)
      include 'implno.dek'
c..
c..given arrays xa and ya of length n and a value x, this routine returns a 
c..value y and an error estimate dy. if p(x) is the polynomial of degree n-1
c..such that ya = p(xa) ya then the returned value is y = p(x) 
c..
c..declare
      integer          n,nmax,ns,i,m
      parameter        (nmax=10)
      double precision xa(n),ya(n),x,y,dy,c(nmax),d(nmax),dif,dift,
     1                 ho,hp,w,den

c..find the index ns of the closest table entry; initialize the c and d tables
      ns  = 1
      dif = abs(x - xa(1))
      do i=1,n
       dift = abs(x - xa(i))
       if (dift .lt. dif) then
        ns  = i
        dif = dift
       end if
       c(i)  = ya(i)
       d(i)  = ya(i)
      enddo

c..first guess for y
      y = ya(ns)

c..for each column of the table, loop over the c's and d's and update them
      ns = ns - 1
      do m=1,n-1
       do i=1,n-m
        ho   = xa(i) - x
        hp   = xa(i+m) - x
        w    = c(i+1) - d(i)
        den  = ho - hp
        if (den .eq. 0.0) stop ' 2 xa entries are the same in polint'
        den  = w/den
        d(i) = hp * den
        c(i) = ho * den
       enddo

c..after each column is completed, decide which correction c or d, to add
c..to the accumulating value of y, that is, which path to take in the table
c..by forking up or down. ns is updated as we go to keep track of where we
c..are. the last dy added is the error indicator.
       if (2*ns .lt. n-m) then
        dy = c(ns+1)
       else
        dy = d(ns)
        ns = ns - 1
       end if
       y = y + dy
      enddo
      return
      end










      subroutine fermint(start,stptry,stpmin,stopp,bc,
     1                   eps,dxsav,kmax, 
     2                   xrk,yrk,xphys,yphys,xlogi,ylogi,
     3                   nok,nbad,kount,odescal,iprint,
     4                   derivs,steper)  
      include 'implno.dek'

c..generic ode driver, tuned a bit for polytropes


c..declare the pass
      external         derivs,steper
      integer          xphys,yphys,xlogi,ylogi,nok,nbad,
     1                 kmax,kount,iprint
      double precision start,stptry,stpmin,stopp,bc(yphys),eps,
     2                 dxsav,xrk(xphys),yrk(yphys,xphys),odescal 


c..local variables
      integer          nmax,stpmax,i,j,nstp
      parameter        (nmax = 20, stpmax=10000)  
      double precision yscal(nmax),y(nmax),dydx(nmax),  
     1                 x,xsav,h,hdid,hnext,zero,one,tiny
      parameter        (zero=0.0d0, one=1.0d0, tiny=1.0d-15)



c..stop integrating when the density or pressure drop below these numbers
      double precision stopden,stop_pres
      common  /sbyiol/ stopden,stop_pres



c..here are the format statements for printouts as we integrate
100   format(1x,i4,1p12e10.2)



c..initialize   
      if (ylogi .gt. yphys) stop 'ylogi > yphys in routine fermint'
      if (yphys .gt. nmax)  stop 'yphys > nmax in routine fermint'
      x     = start   
      h     = sign(stptry,stopp-start) 
      nok   = 0 
      nbad  = 0
      kount = 0   


c..store the first step 
      do i=1,ylogi
       y(i) = bc(i)  
      enddo
      xsav = x - 2.0d0 * dxsav



c..take at most stpmax steps
      do nstp=1,stpmax
       call derivs(x,y,dydx)


c..scaling vector used to monitor accuracy  
       do i=1,ylogi

c..constant fractional accuracy        
        yscal(i) = abs(y(i)) + tiny

c..const. frac. cept near zero  
c        yscal(i)=abs(y(i)) + abs(h*dydx(i)) + tiny  

c..step size dependent accuracy   
c         yscal(i) = abs(odescal * h * dydx(i)) + tiny  

c..for stiffs
c        yscal(i) = max(odescal,y(i))

c..strait scaling (decrease to get more time steps)
c        yscal(i) = odescal
       enddo



c..store intermediate results   
       if (kmax .gt. 0) then
        if ( (abs(dxsav) - abs(x-xsav)) .le. tiny) then 
         if ( kount .lt. (kmax-1) ) then  
          kount         = kount+1  
          xrk(kount)    = x   
          do i=1,ylogi 
           yrk(i,kount) = y(i)
          enddo
          if (iprint .eq. 1) then
           write(6,100) kount,xrk(kount),(yrk(j,kount), j=1,ylogi)
          end if
          xsav=x 
         end if
        end if  
       end if



c..if the step can overshoot the stop point or the dxsav increment then cut it
       if ((x+h-stopp)*(x+h-start) .gt. zero) h = stopp - x  
       if (dxsav.ne.zero .and. h.gt.(xsav-x+dxsav)) h = xsav + dxsav-x


c..do an integration step
       call steper(y,dydx,ylogi,x,h,eps,yscal,hdid,hnext,derivs)   
       if (hdid.eq.h) then
        nok = nok+1   
       else 
        nbad = nbad+1 
       end if



c..this is the normal exit point, save the final step   
       if (nstp .eq. stpmax .or. (x-stopp)*(stopp-start) .ge. zero  .or.
     1     y(4) .le. stop_pres) then
        do i=1,ylogi  
         bc(i) = y(i) 
        enddo
        if (kmax.ne.0) then   
         kount         = kount+1  
         xrk(kount)    = x   
         do i=1,ylogi 
          yrk(i,kount) = y(i) 
         enddo
         if (iprint .eq. 1) then
           write(6,100) kount,xrk(kount),(yrk(j,kount), j=1,ylogi)
         end if
        end if
        return  
       end if


c..set the step size for the next iteration; stay above stpmin
       h=hnext
       if (abs(hnext).lt.stpmin) return
c       if (abs(hnext).lt.stpmin) stop 'hnext < stpmin in fermint'


c..back for another iteration or death
      enddo
      write(6,*) '> than stpmax steps required in fermint' 
      return
      end







c..explicit: 
c..routine rk4 is a fourth order runge-kutta stepper 
c..routine rkqc is a step doubling rk4 integrator; plug into odeint 




       subroutine rk4(y,dydx,n,x,h,yout,derivs)  
       include 'implno.dek' 
c..  
c..given values for the variables y(1:n) and their derivatives dydx(1:n) known 
c..at x, use the fourth order runge-kutta method to advance the solution over 
c..an interval h and return the incremented variables in yout(1:n) (which need 
c..not be a distinct array from y). one supplies the routine derivs which  
c..evaluates the right hand side of the ode's.    
c..  
c..declare   
       external          derivs  
       integer           n,nmax,i    
       parameter         (nmax = 2000)   
       double precision  x,h,y(n),dydx(n),yout(n), 
     1                   yt(nmax),dyt(nmax),dym(nmax), 
     2                   hh,h6,xh    

c..initialize the step sizes and weightings  
       hh = h*0.5d0 
       h6 = h/6.0d0  
       xh = x + hh   

c..the first step    
       do i=1,n   
        yt(i) = y(i) + hh*dydx(i)   
       enddo

c..the second step   
       call derivs(xh,yt,dyt)    
       do i=1,n   
        yt(i) = y(i) + hh*dyt(i)    
       enddo

c..the third step    
       call derivs(xh,yt,dym)    
       do i=1,n   
        yt(i)  = y(i) + h*dym(i) 
        dym(i) = dyt(i) + dym(i)    
       enddo

c..the fourth step and accumulate the increments with the proper weights 
       call derivs(x+h,yt,dyt)   
       do i=1,n   
        yout(i) = y(i) + h6*(dydx(i) +dyt(i) + 2.0d0*dym(i))  
       enddo
       return    
       end   





      subroutine rkdumb(start,step,stopp,bc,xrk,yrk,xdim,ydim,derivs) 
      include 'implno.dek' 
c.. 
c..starting from initial values bc(1:ydim) known as start, use fourth order 
c..runge-kutta to advance the solution to stop using step as the increment. 
c..deriv supplies the right hand side of the ode's. 
c.. 
c..declare 
      external          derivs   
      integer           xdim,ydim,i,k,numstp,nmax   
      parameter         (nmax = 2000)   
      double precision  start,step,stopp,bc(ydim),xrk(ydim), 
     1                  yrk(ydim,xdim),x,y(nmax),dy(nmax)   

c..get the number of steps to take   
      numstp = int((stopp-start)/step)
      if (numstp .gt. xdim) stop ' numstp > xdim in rkdumb'    
      if (ydim   .gt. nmax) stop ' ydim   > nmax in rkdumb'    

c..set the initial conditions   
      x      = start 
      xrk(1) = x 
      do i=1,ydim 
       y(i)     = bc(i)  
       yrk(i,1) = y(i) 
      enddo

c..step on through 
      do i=1,numstp   
       call derivs(x,y,dy)   
       call rk4(y,dy,ydim,x,step,y,derivs)   
       if (x+step .eq. x) stop 'step size too small in rkdumb' 
       x = x + step 
       xrk(i+1) = x    
       do k=1,ydim    
        yrk(k,i+1) = y(k)  
       enddo
      enddo
      return 
      end    






      subroutine rkqc(y,dydx,n,x,htry,eps,yscal,hdid,hnext,derivs)   
      include 'implno.dek' 
c..  
c..fifth order, step doubling, runge-kutta ode integrator with monitering of  
c..local truncation errors. input are the vector y of length n, which has a  
c..known the derivative dydx at the point x, the step size to be attempted  
c..htry, the required accuracy eps, and the vector yscal against which the  
c..error is to be scaled. on output, y and x are replaced by their new values,
c..hdid is the step size that was actually accomplished, and hnext is the  
c..estimated next step size. derivs is a user supplied routine that computes  
c..the right hand side of the first order system of ode's. plug into odeint. 
c..  
c..declare   
      external         derivs    
      integer          n,nmax,i  
      parameter        (nmax = 2000)   
      double precision x,htry,eps,hdid,hnext,y(n),dydx(n),yscal(n),  
     1                 ytemp(nmax),ysav(nmax),dysav(nmax),fcor,safety, 
     2                 errcon,pgrow,pshrnk,xsav,h,hh,errmax  
      parameter        (fcor=1.0d0/15.0d0, pgrow = -0.2d0,  
     1                  pshrnk = -0.25d0,  safety=0.9d0,  errcon=6.0e-4) 

c..note errcon = (4/safety)**(1/pgrow)   
c..nmax is the maximum number of differential equations 


c..save the initial values  
      h      = htry 
      xsav   =  x 
      do i=1,n    
       ysav(i)  = y(i) 
       dysav(i) = dydx(i) 
      enddo

c..take two half steps   
1     hh = 0.5d0*h   
      call rk4(ysav,dysav,n,xsav,hh,ytemp,derivs)    
      x  = xsav + hh  
      call derivs(x,ytemp,dydx)  
      call rk4(ytemp,dydx,n,x,hh,y,derivs)   
      x  = xsav + h   
      if (x .eq. xsav) stop 'stepsize not significant in rkqc'  

c..now take the large step   
      call rk4(ysav,dysav,n,xsav,h,ytemp,derivs) 

c..ytemp is the error estimate   
      errmax = 0.0 
      do i=1,n    
       ytemp(i) = y(i) - ytemp(i)   
       errmax   = max(errmax,abs(ytemp(i)/yscal(i)))    
      enddo
      errmax     = errmax/eps  

c..truncation error too big, reduce the step size and try again  
      if (errmax .gt. 1.0) then 
       h = safety * h * (errmax**pshrnk)  
       go to  1   

c..truncation within limits, compute the size of the next step   
      else   
       hdid = h   
       if (errmax.gt.errcon) then 
        hnext = safety * h * (errmax**pgrow) 
       else 
        hnext = 4.0d0 * h 
       end if    
      end if  

c..mop up the fifth order truncation error   
      do i=1,n    
       y(i) = y(i) + ytemp(i)*fcor  
      enddo
      return 
      end