program test_chempot implicit none save c..tests the chemical potential interpolations c..declare integer i,j,imax,jmax parameter (imax = 200, jmax=200) double precision temp,den,ye,etaele,etapos, 1 tlo,thi,tstp,dlo,dhi,dstp,eta,etamax,etamin, 2 switch tlo = 4.0d0 thi = 12.0d0 tstp = (thi - tlo)/float(jmax-1) dlo = -11.0d0 dhi = 14.0d0 dstp = (dhi - dlo)/float(imax-1) ye = 0.5d0 etamax = -1.0e20 etamin = 1.0e20 switch = 1.0d-30 open(unit=2,file='eta.dat',status='unknown') do j=1,jmax temp = tlo + (j-1)*tstp temp = 10.0d0**temp do i=1,imax den = dlo + (i-1)*dstp den = 10.0d0**den call chempot(temp,den,ye,etaele,etapos) if (etaele .lt. -switch) then eta = -log10(abs(etaele)) else if (etaele .gt. switch) then eta = log10(abs(etaele)) else eta = abs(etaele) end if write(2,10) log10(temp),log10(den),etaele 10 format(1x,1p3e14.6) etamax = max(etamax,etaele) etamin = min(etamin,etaele) enddo write(2,*) enddo close(unit=2) write(6,10) etamin,etamax end subroutine chempot(temp,den,ye,etaele,etapos) implicit none save c..given a temperature temp [K], density den [g/cm**3], and a composition c..characterized by y_e, this routine returns the electron and positron c..chemical potentials. derivatives with respect to temperature and c,,density are computed, but not returned. c.. c..references: timmes & swesty apj 1999 c..declare the pass double precision temp,den,ye,etaele,etapos c..declare local variables c..for the tables integer i,j,imax,jmax parameter (imax = 271, jmax = 101) c..for the density and temperature table double precision d(imax),t(jmax) c..for chemical potential table double precision ef(imax,jmax),efd(imax,jmax), 1 eft(imax,jmax),efdt(imax,jmax) c..for storing the differences double precision dt_sav(jmax), 1 dti_sav(jmax), 2 dd_sav(imax), 3 ddi_sav(imax) c..for the interpolations integer iat,jat double precision tlo,thi,tstp,tstpi,dlo,dhi,dstp,dstpi, 1 tsav,dsav,deni,tempi,kt,ktinv,beta,dbetadt double precision dth,dti,dd,ddi,xt,xd,mxt,mxd, 1 si0t,si1t,si0mt,si1mt, 2 si0d,si1d,si0md,si1md, 3 dsi0t,dsi1t,dsi0mt,dsi1mt, 4 dsi0d,dsi1d,dsi0md,dsi1md, 5 x,z,din,fi(16), 6 xpsi0,xdpsi0,xpsi1,xdpsi1,h3, 7 w0t,w1t,w0mt,w1mt,w0d,w1d,w0md,w1md, 8 detadd,detadt,detapdd,detapdt c..physical constants and parameters double precision clight,kerg,me,mecc,positron_start parameter (clight = 2.99792458d10, 1 kerg = 1.380658d-16, 2 me = 9.1093897d-28, 3 mecc = me * clight * clight, 4 positron_start = 0.02d0) c..for initialization integer ifirst data ifirst/0/ c..cubic hermite polynomial statement functions c..psi0 & derivatives xpsi0(z) = z * z * (2.0d0*z - 3.0d0) + 1.0 xdpsi0(z) = z * (6.0d0*z - 6.0d0) c..psi1 & derivatives xpsi1(z) = z * ( z * (z - 2.0d0) + 1.0d0) xdpsi1(z) = z * (3.0d0*z - 4.0d0) + 1.0d0 c..bicubic hermite polynomial statement function h3(i,j,w0t,w1t,w0mt,w1mt,w0d,w1d,w0md,w1md) = 1 fi(1) *w0d*w0t + fi(2) *w0md*w0t 2 + fi(3) *w0d*w0mt + fi(4) *w0md*w0mt 3 + fi(5) *w0d*w1t + fi(6) *w0md*w1t 4 + fi(7) *w0d*w1mt + fi(8) *w0md*w1mt 5 + fi(9) *w1d*w0t + fi(10) *w1md*w0t 6 + fi(11) *w1d*w0mt + fi(12) *w1md*w0mt 7 + fi(13) *w1d*w1t + fi(14) *w1md*w1t 8 + fi(15) *w1d*w1mt + fi(16) *w1md*w1mt c..read the table and construct the deltas only once if (ifirst .eq. 0) then ifirst = 1 open(unit=19,file='chem_table.dat',status='old') tlo = 3.0d0 thi = 13.0d0 tstp = (thi - tlo)/float(jmax-1) tstpi = 1.0d0/tstp dlo = -12.0d0 dhi = 15.0d0 dstp = (dhi - dlo)/float(imax-1) dstpi = 1.0d0/dstp do j=1,jmax tsav = tlo + (j-1)*tstp t(j) = 10.0d0**(tsav) do i=1,imax dsav = dlo + (i-1)*dstp d(i) = 10.0d0**(dsav) read(19,*) ef(i,j),efd(i,j),eft(i,j),efdt(i,j) enddo enddo c..construct the temperature and density deltas and their inverses do j=1,jmax-1 dth = t(j+1) - t(j) dti = 1.0d0/dth dt_sav(j) = dth dti_sav(j) = dti end do do i=1,imax-1 dd = d(i+1) - d(i) ddi = 1.0d0/dd dd_sav(i) = dd ddi_sav(i) = ddi enddo c..close the file and write a diagnostic close(unit=19) write(6,*) write(6,*) 'finished reading chemical potential table' write(6,04) 'imax=',imax,' jmax=',jmax 04 format(1x,4(a,i4)) write(6,03) 'temp(1) =',t(1),' temp(jmax) =',t(jmax) write(6,03) 'ye*den(1) =',d(1),' ye*den(imax) =',d(imax) 03 format(1x,4(a,1pe11.3)) write(6,*) end if c..normal execution starts here c..enter the table with ye*den din = ye*den c..bomb proof the input if (temp .gt. t(jmax)) then write(6,01) 'temp=',temp,' t(jmax)=',t(jmax) 01 format(1x,5(a,1pe11.3)) write(6,*) 'temp too hot, off grid, returning' return end if if (temp .lt. t(1)) then write(6,01) 'temp=',temp,' t(1)=',t(1) write(6,*) 'temp too cold, off grid, returning' return end if if (din .gt. d(imax)) then write(6,01) 'den*ye=',din,' d(imax)=',d(imax) write(6,*) 'ye*den too big, off grid, returning' return end if if (din .lt. d(1)) then write(6,01) 'ye*den=',din,' d(1)=',d(1) write(6,*) 'ye*den too small, off grid, returning' return end if c..initialize deni = 1.0d0/den tempi = 1.0d0/temp kt = kerg * temp ktinv = 1.0d0/kt beta = kt/mecc dbetadt = kerg/mecc c..hash locate this temperature and density jat = int((log10(temp) - tlo)*tstpi) + 1 jat = max(1,min(jat,jmax-1)) iat = int((log10(din) - dlo)*dstpi) + 1 iat = max(1,min(iat,imax-1)) c..look in the electron chemical potential table only once fi(1) = ef(iat,jat) fi(2) = ef(iat+1,jat) fi(3) = ef(iat,jat+1) fi(4) = ef(iat+1,jat+1) fi(5) = eft(iat,jat) fi(6) = eft(iat+1,jat) fi(7) = eft(iat,jat+1) fi(8) = eft(iat+1,jat+1) fi(9) = efd(iat,jat) fi(10) = efd(iat+1,jat) fi(11) = efd(iat,jat+1) fi(12) = efd(iat+1,jat+1) fi(13) = efdt(iat,jat) fi(14) = efdt(iat+1,jat) fi(15) = efdt(iat,jat+1) fi(16) = efdt(iat+1,jat+1) c..get the interpolation weight functions xt = max( (temp - t(jat))*dti_sav(jat), 0.0d0) xd = max( (din - d(iat))*ddi_sav(iat), 0.0d0) mxt = 1.0d0 - xt mxd = 1.0d0 - xd si0t = xpsi0(xt) si1t = xpsi1(xt)*dt_sav(jat) si0mt = xpsi0(mxt) si1mt = -xpsi1(mxt)*dt_sav(jat) si0d = xpsi0(xd) si1d = xpsi1(xd)*dd_sav(iat) si0md = xpsi0(mxd) si1md = -xpsi1(mxd)*dd_sav(iat) c..derivatives of weight functions dsi0t = xdpsi0(xt)*dti_sav(jat) dsi1t = xdpsi1(xt) dsi0mt = -xdpsi0(mxt)*dti_sav(jat) dsi1mt = xdpsi1(mxt) dsi0d = xdpsi0(xd)*ddi_sav(iat) dsi1d = xdpsi1(xd) dsi0md = -xdpsi0(mxd)*ddi_sav(iat) dsi1md = xdpsi1(mxd) c..electron chemical potential etaele = h3(iat,jat, 1 si0t, si1t, si0mt, si1mt, 2 si0d, si1d, si0md, si1md) c..derivative with respect to density x = h3(iat,jat, 1 si0t, si1t, si0mt, si1mt, 2 dsi0d, dsi1d, dsi0md, dsi1md) detadd = ye * x c..derivative with respect to temperature detadt = h3(iat,jat, 1 dsi0t, dsi1t, dsi0mt, dsi1mt, 2 si0d, si1d, si0md, si1md) c..positron chemical potential etapos = 0.0d0 detapdd = 0.0d0 detapdt = 0.0d0 if (beta .gt. positron_start) then etapos = -etaele - 2.0d0/beta detapdd = -detadd detapdt = -detadt + 2.0d0/beta**2 * dbetadt end if return end