c======================================================================c
c
c      program readeos
c
	subroutine readeos(tin,yin,rin,fer,pre,ent,ene,cnu,cpu,ceu,cnv,anu
     &  ,pnu,xnu,xpu,xau,xiu,emm)
c
c======================================================================c
c
c******************** usage explanation   *****************************c

c	Default setup for it is to use as subroutine. 
c     Call load_table once in the main routine you have.
c
c	if used as main routine:
c	   uncomment call load_table in this subroutine,
c	   and uncomment the test input parameters,
c	   and comment the return command.
c	
c	last update:
c     G. Shen. 3/5/10/
c**********************************************************************c
c
c---- input parameters for interpolation
c
c	tin: T     [MeV]
c	yin: Yp    
c     rin: n     [fm^-3]
c
c---- output
c
c	tin: T     [MeV]
c	yin: Yp    
c     rin: n     [fm^-3]
c
c	fer: free energy per baryon [MeV/baryon] relative to nucleon mass: 939 MeV.
c	pre: pressure		        [MeV/fm**3]
c	ent: entropy                [kb/baryon]
c	ene: internal energy        [MeV/baryon] relative to nucleon mass
c	cnu: n chemical potential   [MeV]
c	cpu: p chemical potential   [MeV]
c	ceu: e chemical potential   [MeV]
c	cnv: equilibrium electron neutrino chemical potential [MeV]
c	anu: mass number of heavies (it deos not include alpha particles.)
c	pnu: proton number of heavies (same as above.)
c	xnu: neutron mass fraction
c	xpu: proton mass fraction
c	xau: alpha mass fraction
c	xiu: ion mass fraction
c	emm: effective mass of nucleon.         [MeV]
c
c----------------------------------------------------------------------c
      implicit real*8(a-h,o-z)

	common /readin/ t1(700,240,53),y1(700,240,53),d1(700,240,53),
     &f1(700,240,53),p1(700,240,53),s1(700,240,53),cn1(700,240,53),
     &cp1(700,240,53),an1(700,240,53),pn1(700,240,53),xn1(700,240,53),
     &xp1(700,240,53),xa1(700,240,53),xi1(700,240,53),em1(700,240,53),
     &ce1(700,240,53),e1(700,240,53)

	common /work/   t2(700,240,53),y2(700,240,53),d2(700,240,53)
		data int/218/, inb/657/
c
c
c---- comment this line if used as subroutine and call load_table in the
c	main subroutine.
c
c	call load_table
c
c---- input parameters: comment these lines if used as subroutine.
c	tin = 0.500
c	yin = 0.065
c	rin = 2.d-3
c
c---- bracket the location of input point in 3-d phase space
	do i = 2, int
         if(tin.ge.t1(1,i-1,1).and.tin.lt.t1(1,i,1)) then
	     iff = i-1
	     goto 10
         endif
	enddo
	   write(*,*) 'tin =  ', tin
	stop  'need T: [0,80] MeV'
c
   10 do j = 2, 53
         if(yin.ge.y1(1,1,j-1).and.yin.lt.y1(1,1,j)) then
	     jff = j-1
	     goto 20
	   endif
	enddo
   	stop  'need Yp: [0,0.56]'
c
   20 do k = 2, inb
         if(rin.ge.d1(k-1,1,1).and.rin.lt.d1(k,1,1)) then
	     kff = k-1
	     goto 30
	   endif
      enddo	
	  write(*,*) tin,yin,rin
	stop 'need nb: [10^-8,1.6]'
c
c   30 write(*,*) 'bracket T:  ', iff, iff+1
c	write(*,*) 'bracket Yp: ', jff, jff+1
c	write(*,*) 'bracket nb: ', kff, kff+1
c
c calculating the weight factors
  30	grho = (log10(rin) - d2(kff,iff,jff))/(d2(kff+1,iff,jff)
     &	-d2(kff,iff,jff))
        
	If(iff.eq.1) then
	  gtem = (tin - t1(kff,iff,jff))/(t1(kff,iff+1,jff)-
     &	  t1(kff,iff,jff))
        else
          gtem = (log10(tin) - t2(kff,iff,jff))/(t2(kff,iff+1,jff)-
     &	  t2(kff,iff,jff))
	endif

	gye  = (yin-y1(kff,iff,jff))/(y1(kff,iff,jff+1)-y1(kff,iff,jff))
	
c Interpolating:

c free energy:
	 fer           = (1-grho)*(1-gtem)*(1-gye) *
     &                    f1(kff,iff,jff)              +
     &                   (1-grho)*(gtem)*(1-gye)   *
     &                     f1(kff,iff+1,jff)           +
     &                   (grho)*(1-gtem)*(1-gye)   *
     &                     f1(kff+1,iff,jff)           +
     &                   (grho)*(gtem)*(1-gye)     *
     &                     f1(kff+1,iff+1,jff)         +
     &                   (1-grho)*(1-gtem)*gye     *
     &                     f1(kff,iff,jff+1)           +
     &                   (1-grho)*(gtem)*(gye)     *
     &                     f1(kff,iff+1,jff+1)         +
     &                   (grho)*(1-gtem)*(gye)     *
     &                     f1(kff+1,iff,jff+1)         +
     &                   (grho)*(gtem)*(gye)       * 
     &                     f1(kff+1,iff+1,jff+1)
c
c pressure:
	 pre           = (1-grho)*(1-gtem)*(1-gye) *
     &                     p1(kff,iff,jff)             +
     &                   (1-grho)*(gtem)*(1-gye)   *
     &                     p1(kff,iff+1,jff)           +
     &                   (grho)*(1-gtem)*(1-gye)   *
     &                     p1(kff+1,iff,jff)           +
     &                   (grho)*(gtem)*(1-gye)     *
     &                     p1(kff+1,iff+1,jff)         +
     &                   (1-grho)*(1-gtem)*gye     *
     &                     p1(kff,iff,jff+1)           +
     &                   (1-grho)*(gtem)*(gye)     *
     &                     p1(kff,iff+1,jff+1)         +
     &                   (grho)*(1-gtem)*(gye)     *
     &                     p1(kff+1,iff,jff+1)         +
     &                   (grho)*(gtem)*(gye)       * 
     &                     p1(kff+1,iff+1,jff+1)
c
c entropy:
	 ent           = (1-grho)*(1-gtem)*(1-gye) *
     &                    s1(kff,iff,jff)              +
     &                   (1-grho)*(gtem)*(1-gye)   *
     &                     s1(kff,iff+1,jff)           +
     &                   (grho)*(1-gtem)*(1-gye)   *
     &                     s1(kff+1,iff,jff)           +
     &                   (grho)*(gtem)*(1-gye)     *
     &                     s1(kff+1,iff+1,jff)         +
     &                   (1-grho)*(1-gtem)*gye     *
     &                     s1(kff,iff,jff+1)           +
     &                   (1-grho)*(gtem)*(gye)     *
     &                     s1(kff,iff+1,jff+1)         +
     &                   (grho)*(1-gtem)*(gye)     *
     &                     s1(kff+1,iff,jff+1)         +
     &                   (grho)*(gtem)*(gye)       * 
     &                     s1(kff+1,iff+1,jff+1)
c
c internal energy relative to nucleon mass
	ene = fer + tin*ent
c
c 
c	 ene           = (1-grho)*(1-gtem)*(1-gye) *
c     &                    e1(kff,iff,jff)              +
c     &                   (1-grho)*(gtem)*(1-gye)   *
c     &                     e1(kff,iff+1,jff)           +
c     &                   (grho)*(1-gtem)*(1-gye)   *
c     &                     e1(kff+1,iff,jff)           +
c     &                   (grho)*(gtem)*(1-gye)     *
c     &                     e1(kff+1,iff+1,jff)         +
c     &                   (1-grho)*(1-gtem)*gye     *
c     &                     e1(kff,iff,jff+1)           +
cc     &                   (1-grho)*(gtem)*(gye)     *
c     &                     e1(kff,iff+1,jff+1)         +
c     &                   (grho)*(1-gtem)*(gye)     *
c     &                     e1(kff+1,iff,jff+1)         +
c     &                   (grho)*(gtem)*(gye)       * 
c     &                     e1(kff+1,iff+1,jff+1)

c
c neutron chemical potential:
	 cnu           = (1-grho)*(1-gtem)*(1-gye) *
     &                    cn1(kff,iff,jff)              +
     &                   (1-grho)*(gtem)*(1-gye)   *
     &                     cn1(kff,iff+1,jff)           +
     &                   (grho)*(1-gtem)*(1-gye)   *
     &                     cn1(kff+1,iff,jff)           +
     &                   (grho)*(gtem)*(1-gye)     *
     &                     cn1(kff+1,iff+1,jff)         +
     &                   (1-grho)*(1-gtem)*gye     *
     &                     cn1(kff,iff,jff+1)           +
     &                   (1-grho)*(gtem)*(gye)     *
     &                     cn1(kff,iff+1,jff+1)         +
     &                   (grho)*(1-gtem)*(gye)     *
     &                     cn1(kff+1,iff,jff+1)         +
     &                   (grho)*(gtem)*(gye)       * 
     &                     cn1(kff+1,iff+1,jff+1)
c
c proton chemical potential:
	 cpu           = (1-grho)*(1-gtem)*(1-gye) *
     &                    cp1(kff,iff,jff)              +
     &                   (1-grho)*(gtem)*(1-gye)   *
     &                     cp1(kff,iff+1,jff)           +
     &                   (grho)*(1-gtem)*(1-gye)   *
     &                     cp1(kff+1,iff,jff)           +
     &                   (grho)*(gtem)*(1-gye)     *
     &                     cp1(kff+1,iff+1,jff)         +
     &                   (1-grho)*(1-gtem)*gye     *
     &                     cp1(kff,iff,jff+1)           +
     &                   (1-grho)*(gtem)*(gye)     *
     &                     cp1(kff,iff+1,jff+1)         +
     &                   (grho)*(1-gtem)*(gye)     *
     &                     cp1(kff+1,iff,jff+1)         +
     &                   (grho)*(gtem)*(gye)       * 
     &                     cp1(kff+1,iff+1,jff+1)
c
c electron chemical potential:
	 ceu           = (1-grho)*(1-gtem)*(1-gye) *
     &                    ce1(kff,iff,jff)              +
     &                   (1-grho)*(gtem)*(1-gye)   *
     &                     ce1(kff,iff+1,jff)           +
     &                   (grho)*(1-gtem)*(1-gye)   *
     &                     ce1(kff+1,iff,jff)           +
     &                   (grho)*(gtem)*(1-gye)     *
     &                     ce1(kff+1,iff+1,jff)         +
     &                   (1-grho)*(1-gtem)*gye     *
     &                     ce1(kff,iff,jff+1)           +
     &                   (1-grho)*(gtem)*(gye)     *
     &                     ce1(kff,iff+1,jff+1)         +
     &                   (grho)*(1-gtem)*(gye)     *
     &                     ce1(kff+1,iff,jff+1)         +
     &                   (grho)*(gtem)*(gye)       * 
     &                     ce1(kff+1,iff+1,jff+1)
c
c equilibrium electron neutrino chemical potential
	 cnv = cpu + ceu - cnu
c
c mass number:
	 anu           = (1-grho)*(1-gtem)*(1-gye) *
     &                    an1(kff,iff,jff)              +
     &                   (1-grho)*(gtem)*(1-gye)   *
     &                     an1(kff,iff+1,jff)           +
     &                   (grho)*(1-gtem)*(1-gye)   *
     &                     an1(kff+1,iff,jff)           +
     &                   (grho)*(gtem)*(1-gye)     *
     &                     an1(kff+1,iff+1,jff)         +
     &                   (1-grho)*(1-gtem)*gye     *
     &                     an1(kff,iff,jff+1)           +
     &                   (1-grho)*(gtem)*(gye)     *
     &                     an1(kff,iff+1,jff+1)         +
     &                   (grho)*(1-gtem)*(gye)     *
     &                     an1(kff+1,iff,jff+1)         +
     &                   (grho)*(gtem)*(gye)       * 
     &                     an1(kff+1,iff+1,jff+1)
c
c proton number:
	 pnu           = (1-grho)*(1-gtem)*(1-gye) *
     &                    pn1(kff,iff,jff)              +
     &                   (1-grho)*(gtem)*(1-gye)   *
     &                     pn1(kff,iff+1,jff)           +
     &                   (grho)*(1-gtem)*(1-gye)   *
     &                     pn1(kff+1,iff,jff)           +
     &                   (grho)*(gtem)*(1-gye)     *
     &                     pn1(kff+1,iff+1,jff)         +
     &                   (1-grho)*(1-gtem)*gye     *
     &                     pn1(kff,iff,jff+1)           +
     &                   (1-grho)*(gtem)*(gye)     *
     &                     pn1(kff,iff+1,jff+1)         +
     &                   (grho)*(1-gtem)*(gye)     *
     &                     pn1(kff+1,iff,jff+1)         +
     &                   (grho)*(gtem)*(gye)       * 
     &                     pn1(kff+1,iff+1,jff+1)
c
c neutron mass fraction:
	 xnu           = (1-grho)*(1-gtem)*(1-gye) *
     &                    xn1(kff,iff,jff)              +
     &                   (1-grho)*(gtem)*(1-gye)   *
     &                     xn1(kff,iff+1,jff)           +
     &                   (grho)*(1-gtem)*(1-gye)   *
     &                     xn1(kff+1,iff,jff)           +
     &                   (grho)*(gtem)*(1-gye)     *
     &                     xn1(kff+1,iff+1,jff)         +
     &                   (1-grho)*(1-gtem)*gye     *
     &                     xn1(kff,iff,jff+1)           +
     &                   (1-grho)*(gtem)*(gye)     *
     &                     xn1(kff,iff+1,jff+1)         +
     &                   (grho)*(1-gtem)*(gye)     *
     &                     xn1(kff+1,iff,jff+1)         +
     &                   (grho)*(gtem)*(gye)       * 
     &                     xn1(kff+1,iff+1,jff+1)
c
c proton mass fraction:
	 xpu           = (1-grho)*(1-gtem)*(1-gye) *
     &                    xp1(kff,iff,jff)              +
     &                   (1-grho)*(gtem)*(1-gye)   *
     &                     xp1(kff,iff+1,jff)           +
     &                   (grho)*(1-gtem)*(1-gye)   *
     &                     xp1(kff+1,iff,jff)           +
     &                   (grho)*(gtem)*(1-gye)     *
     &                     xp1(kff+1,iff+1,jff)         +
     &                   (1-grho)*(1-gtem)*gye     *
     &                     xp1(kff,iff,jff+1)           +
     &                   (1-grho)*(gtem)*(gye)     *
     &                     xp1(kff,iff+1,jff+1)         +
     &                   (grho)*(1-gtem)*(gye)     *
     &                     xp1(kff+1,iff,jff+1)         +
     &                   (grho)*(gtem)*(gye)       * 
     &                     xp1(kff+1,iff+1,jff+1)
c
c alpha mass fraction:
	 xau           = (1-grho)*(1-gtem)*(1-gye) *
     &                    xa1(kff,iff,jff)              +
     &                   (1-grho)*(gtem)*(1-gye)   *
     &                     xa1(kff,iff+1,jff)           +
     &                   (grho)*(1-gtem)*(1-gye)   *
     &                     xa1(kff+1,iff,jff)           +
     &                   (grho)*(gtem)*(1-gye)     *
     &                     xa1(kff+1,iff+1,jff)         +
     &                   (1-grho)*(1-gtem)*gye     *
     &                     xa1(kff,iff,jff+1)           +
     &                   (1-grho)*(gtem)*(gye)     *
     &                     xa1(kff,iff+1,jff+1)         +
     &                   (grho)*(1-gtem)*(gye)     *
     &                     xa1(kff+1,iff,jff+1)         +
     &                   (grho)*(gtem)*(gye)       * 
     &                     xa1(kff+1,iff+1,jff+1)
c
c ion mass fraction:
	 xiu           = (1-grho)*(1-gtem)*(1-gye) *
     &                    xi1(kff,iff,jff)              +
     &                   (1-grho)*(gtem)*(1-gye)   *
     &                     xi1(kff,iff+1,jff)           +
     &                   (grho)*(1-gtem)*(1-gye)   *
     &                     xi1(kff+1,iff,jff)           +
     &                   (grho)*(gtem)*(1-gye)     *
     &                     xi1(kff+1,iff+1,jff)         +
     &                   (1-grho)*(1-gtem)*gye     *
     &                     xi1(kff,iff,jff+1)           +
     &                   (1-grho)*(gtem)*(gye)     *
     &                     xi1(kff,iff+1,jff+1)         +
     &                   (grho)*(1-gtem)*(gye)     *
     &                     xi1(kff+1,iff,jff+1)         +
     &                   (grho)*(gtem)*(gye)       * 
     &                     xi1(kff+1,iff+1,jff+1)
c
c effective mass:
	 emm           = (1-grho)*(1-gtem)*(1-gye) *
     &                    em1(kff,iff,jff)              +
     &                   (1-grho)*(gtem)*(1-gye)   *
     &                     em1(kff,iff+1,jff)           +
     &                   (grho)*(1-gtem)*(1-gye)   *
     &                     em1(kff+1,iff,jff)           +
     &                   (grho)*(gtem)*(1-gye)     *
     &                     em1(kff+1,iff+1,jff)         +
     &                   (1-grho)*(1-gtem)*gye     *
     &                     em1(kff,iff,jff+1)           +
     &                   (1-grho)*(gtem)*(gye)     *
     &                     em1(kff,iff+1,jff+1)         +
     &                   (grho)*(1-gtem)*(gye)     *
     &                     em1(kff+1,iff,jff+1)         +
     &                   (grho)*(gtem)*(gye)       * 
     &                     em1(kff+1,iff+1,jff+1)
c
c	write(*,*) 'outputs '
c	write(*,'(21(1pe15.5))') tin,yin,rin,fer,pre,ent,ene,cnu,cpu,ceu,
c     &  cnv,anu,pnu,xnu,xpu,xau,xiu,emm
c
c      write(*,*) 'interpolation completed'
c
c---- uncomment this line if used as subroutine.
	return
c
      end
c
c========================================================================c
c
	subroutine load_table
c
c     read table by 1D-1D-1D.
c	For each temperature T, proton fraction Yp runs from 0 to 0.56;
c	for each Yp, baryon density nb runs from 10**-8 to 1.6 fm**-3.
c
c	11/12/09
c	modified to readin eos_b.dat only.
c
c========================================================================c
c
	implicit real*8(a-h,o-z)
c		
	common /readin/ t1(700,240,53),y1(700,240,53),d1(700,240,53),
     &f1(700,240,53),p1(700,240,53),s1(700,240,53),cn1(700,240,53),
     &cp1(700,240,53),an1(700,240,53),pn1(700,240,53),xn1(700,240,53),
     &xp1(700,240,53),xa1(700,240,53),xi1(700,240,53),em1(700,240,53),
     &ce1(700,240,53),e1(700,240,53)
	common /work/   t2(700,240,53),y2(700,240,53),d2(700,240,53)
c
		data int/218/, inb/657/
c
c
	m = 1
c---- big table
c
	write(*,*) 'input choice: 1 for Baryon; 2 for all'
c
	read(*,*) is
	if(is.eq.1) then
      	open(90,file='eosb_64.dat',status='unknown')      ! read in baryon only
	    write(*,*) 'loading eosb_64.dat'
	elseif(is.eq.2) then
		open(90,file='eos_64.dat',status='unknown')        ! read in all
	    write(*,*) 'loading eos_64.dat'
	endif
c
	do i = 1, int

	   do j = 1, 53

	       do k = 1, inb
c
	read(90,*) t1(k,i,j),               ! temperature                    [MeV]
     &	       y1(k,i,j),				! proton fraction                  
     &		   d1(k,i,j),				! baryon number density          [fm**-3]
     &		   f1(k,i,j),				! free energy per baryon         [MeV]
     &		   p1(k,i,j),				! pressure                       [MeV fm**-3]
     &		   s1(k,i,j)				! energy per baryon              [MeV]
     &		   ,cn1(k,i,j),				! chemical potential for neutron [MeV]
     &		   cp1(k,i,j),				! chemical potential for proton  [MeV]
     &		   ce1(k,i,j),				! chemical potential for electron[MeV] 
     &		   an1(k,i,j),				! average mass number       
     &		   pn1(k,i,j),				! average proton number
     &		   xn1(k,i,j),				! fraction of free neutron
     &		   xp1(k,i,j),				! fraction of free proton
     &		   xa1(k,i,j),				! fraction of alphas
     &		   xi1(k,i,j),				! fraction of heavy nuclei
     &		   em1(k,i,j)				! effective mass                 [MeV]

c
	d2(k,i,j) = log10(d1(k,i,j))
        if(i.ne.1) t2(k,i,j) = log10(t1(k,i,j))
	       enddo
   	   enddo
	     write(*,*) i, t1(k-1,i,j-1)
	enddo
	close(90)
c
	write(*,*) 'loading completed: number of points =   ', m
c
	return
c	
	end
c======================================================================c
