program tkap implicit none save include 'vector_eos.dek' include 'vector_kap.dek' c.. c..tests the opacity routines c.. c..ionmax = number of isotopes in the network c..xmass = mass fractions c..ymass = molar fractions c..aion = number of nucleons c..zion = number of protons integer ionmax parameter (ionmax=3) double precision xmass(ionmax),ymass(ionmax), 1 aion(ionmax),zion(ionmax), 2 temp,den,abar,zbar, 3 pin,xin,ein,xh,xhe,xz,w4,w5,w6,zalf, 4 orad,ocond,opac c..a format statement 01 format(1x,a6,t10,a6,1pe10.3, 1 t30,a6,1pe10.3,t50,a6,1pe10.3) c..set the mass fractions, z's and a's of the composition c..hydrogen xmass(1) = 0.75d0 aion(1) = 1.0d0 zion(1) = 1.0d0 c..helium xmass(2) = 0.23d0 aion(2) = 4.0d0 zion(2) = 2.0d0 c..carbon 12 xmass(3) = 0.02d0 aion(3) = 12.0d0 zion(3) = 6.0d0 c..get abar, zbar and a few other composition variables c..this routine is included with any of the eos routines from fxt's webpage call azbar(xmass,aion,zion,ionmax, 1 ymass,abar,zbar) c..set the input vector temp_row(1) = 1.0d8 den_row(1) = 1.0d6 abar_row(1) = abar zbar_row(1) = zbar c..here the eos pipeline is only 1 element long jlo_eos = 1 jhi_eos = 1 c..call an equation of state c..you can download this eos from fxt's eos code webpage c..see the end of this file for how its brought in call eosfxt c..one *must* always call the eos before calling the opacity routine c..since some of the required input must come from an eos: c..load the electron-positron pressure and its derivatives pin_row(1) = pele_row(1) + ppos_row(1) dpind_row(1) = dpepd_row(1) dpint_row(1) = dpept_row(1) c..load the electron-positron number densities and its derivatives xin_row(1) = xne_row(1) + xnp_row(1) dxind_row(1) = dxned_row(1) dxint_row(1) = dxnet_row(1) c..load the electron chemical potential and its derivatives ein_row(1) = etaele_row(1) deind_row(1) = detad_row(1) deint_row(1) = detat_row(1) c..get the composition variables call kapbar(xmass,aion,zion,ionmax, 1 xh,xhe,xz,w4,w5,w6,zalf) c..load the composition variables xh_row(1) = xh xhe_row(1) = xhe xz_row(1) = xz w4_row(1) = w4 w5_row(1) = w5 w6_row(1) = w6 zalf_row(1) = zalf c..here the opacity pipeline is only 1 element long jlo_kap = 1 jhi_kap = 1 c..get the opacities and conductivities call kapfxt c..write out what we got write(6,01) 'kapfxt','orad =',orad_row(1), 1 'dord =',dord_row(1), 2 'dort =',dort_row(1) write(6,01) ' ','ocond=',ocond_row(1), 1 'docd =',docd_row(1), 2 'doct =',doct_row(1) write(6,01) ' ','opac =',opac_row(1), 1 'dopd =',dopd_row(1), 2 'dopt =',dopt_row(1) write(6,*) write(6,01) ' ','srad =',s2rad_row(1), 1 'dsrd =',dsrd_row(1), 2 'dsrt =',dsrt_row(1) write(6,01) ' ','scond=',scond_row(1), 1 'dscd =',dscd_row(1), 2 'dsct =',dsct_row(1) write(6,01) ' ','sigma=',sigma_row(1), 1 'dsigd=',dsigd_row(1), 2 'dsigt=',dsigt_row(1) write(6,*) ' ' stop 'normal termination' end subroutine kapbar(xmass,aion,zion,ionmax, 1 xh,xhe,xz,w4,w5,w6,zalf) implicit none save c..this routine calculates composition variables for an opacity routine c..input: c..mass fractions = xmass(1:ionmax) c..number of nucleons = aion(1:ionmax) c..charge of nucleus = zion(1:ionmax) c..number of isotopes = ionmax c.. c..output: c..mass fraction of hydrogen = xh c..mass fraction of helium = xhe c..mass fraction of metals = xz c..charge squared weighted mass fraction of hydrogen = w4 c..charge squared weighted mass fraction of helium = w5 c..charge squared weighted mass fraction of metals = w6 c..average charge squared aion**(2/3) = zalf c..declare integer ionmax,i,iz double precision xmass(ionmax),aion(ionmax),zion(ionmax), 1 xh,xhe,xz,w4,w5,w6,zalf, 2 w(6),ztemp,twoth parameter (twoth = 2.0d0/3.0d0) zalf = 0.0d0 do i=1,6 w(i) = 0.0d0 enddo do i = 1,ionmax iz = min(3,max(1,int(zion(i)))) ztemp = zion(i) * zion(i) * xmass(i)/aion(i) w(iz) = w(iz) + xmass(i) w(iz+3) = w(iz+3) + ztemp zalf = zalf + ztemp * aion(i)**twoth enddo xh = w(1) xhe = w(2) xz = w(3) w4 = w(4) w5 = w(5) w6 = w(6) return end subroutine kapfxt include 'implno.dek' include 'const.dek' include 'vector_eos.dek' include 'vector_kap.dek' c..this routine calculates an opacity/thermal conductivity c..input through common block: c..jlokap = lower limit of pipeline c..jhikap = upper limit of pipeline c..temp_row = temperature temp (in K) c..den_row = density den (in g/cm**3) c..input through common block and returned by, for example, routine azbar: c..abar_row = average number of nucleons in nucleus c..zbar_row = average charge on each nucleus c..input through common block and returned by, for example, routine kapbar: c..xh_row = hydrogen mass fraction c..xhe_row = helium mass fraction c..xz_row = metal mass fraction c..w4_row = charge weighted hydrogen mass fraction c..w5_row = charge weighted helium mass fraction c..w6_row = charge weighted metal mass fraction c..zalf_row = charge squared aion^2/3 weighted molar fractions c..input through common block and returned by an equation of state: c..pin_row = electron+positron pressure (in erg/cm**3) c..dpind_row = derivative of pin_row with density c..dpint_row = derivative of pin_row with temperature c..xin_row = electron+positron number density (in #/cm**3) c..dxind_row = derivative of xin_row with density, from an eos c..dxint_row = derivative of xin_row with temperature, from an eos c..ein_row = electron degeneracy parameter (dimensionless, eta/kt) c..deind_row = derivative of ein_row with density c..deind_row = derivative of ein_row with temperature c..main output through common block: c..orad_row = radiative opacity (in cm**2/g) c..dord_row = derivative of the radiative opacity with density c..dort_row = derivative of the radiative opacity with temperature c..ocond_row = thermal conductivity turned into an opacity (in cm**2/g) c..docd_row = derivative of the thermal opacity with density c..doct_row = derivative of the thermal opacity with temperature c..opac_row = total opacity (in cm**2/g) c..dopd_row = derivative of total opacity with density c..dopt_row = derivative of the total opacity with temperature c..s2rad_row = radiative opacity turned into a conductivity (in erg/cm/K/sec) c..dsrd_row = derivative of the radiative conductivity with density c..dsrt_row = derivative of the radiative conductivity with temperature c..scond_row = thermal conductivity (in erg/cm/K/sec) c..dscd_row = derivative of the thermal conductivity with density c..dsct_row = derivative of the thermal conductivity with temperature c..sigma_row = total thermal conductivity (in erg/cm/K/sec) c..dsigd_row = derivative of total conductivity with density c..dsigt_row = derivative of the total conductivity with temperature c..lots of the various pieces are also avaliable through common block c..declare c..for the thermodynamic state double precision temp,den,abar,zbar,pep,dpepdd,dpepdt,xne, 1 dxnedd,dxnedt,eta,detadd,detadt c..for the composition variables double precision xh,xhe,xz,w4,w5,w6,zalf, 1 xxmue,xxmuebd,xxmuebt,zcz,zbeta,xmu,xmu12,zalpha c..for the iben radiative opacity for small hydrogen compositions double precision oiben1,oiben1bd,oiben1bt, 1 xkc,xkap,xkb,xka,dbar,dbardt, 2 xkcbd,xkcbt,xkapbd,xkapbt,xkrho,xcon1,zbound parameter (xkrho = 0.67d0, 1 xcon1 = 1.0d0/24.55d0, 2 zbound = 0.1d0) c..for the iben radiative opacity for large hydrogen high temp composition double precision oiben2,oiben2bd,oiben2bt, 1 d0log,xka1,xkw,xkaz,dbar1log,dbar2log,dbar1, 2 d0logbt,xka1bt,xkwbt,xkazbd,xkazbt, 3 dbar1logbt,dbar2logbt,dbar1bt,dzfac,xln10 parameter (xln10 = 2.3025850929940459d0) c..for the christy radiative opacity for large hydro low temps composition double precision ochrsty,ochrstybd,ochrstybt, 1 t4,t4r,t43,t44,t45,t46, 2 ck1,ck3,ck2,ck4,ck5,ck6, 3 xkcx,xkcy,xkcz, 4 ck1bt,ck3bd,ck3bt,ck2bd,ck2bt,ck4bt,ck5bt,ck6bt, 5 xkcxbd,xkcxbt,xkcybt,xkczbt c..for easy toggeling of the iben/christy regimes double precision t6_switch1,t6_switch2 parameter (t6_switch1 = 0.5d0, t6_switch2 = 0.9d0) c..orginal limits c parameter (t6_switch1 = 1.0d0, t6_switch2 = 1.5d0) c..for the weaver compton scattering opacity double precision ocompt,ocomptbt,ocomptbd, 1 th,thbt,fact,factbt,facetax,faceta, 2 facetaxbt,facetaxbd,facetabt,facetabd c..for the plasma frequency opacity cutoff double precision tcutoff,cutfac,cutfacbt,cutfacbd c..for fudging the molecular opacity double precision xkf,xkfbt,xkfbd,t7,t73,t74,t7peek parameter (t7peek = 1.0d20) c..for the total radiative opacity and its derivatives double precision orad,doraddd,doraddt,oradmlt c..for the total radiative opacity, converted to a thermal conductivity double precision s2rad,ds2raddd,ds2raddt c..for the iben low density thermal conductivity oh double precision oh,ohbt,ohbd, 1 dlog10,zdel,zdell10,eta0,eta02,thpl,thpla,thplb, 2 cfac1,cfac2,zdelbd,zdelbt,zdell10bd,zdell10bt, 3 eta0bd,eta0bt,eta02bd,eta02bt,thplbd,thplbt, 4 thplabd,thplabt,thplbbd,thplbbt,cfac1bd,cfac1bt, 5 cfac2bd,cfac2bt,pefac,pefacbd,pefacbt,pefaca, 6 pefacabd,pefacabt,pefacal,pefacalbd,pefacalbt, 7 pefacb,pefacbbd,pefacbbt,pefacbl,pefacblbd, 8 pefacblbt,dnefac,dnefacbd,dnefacbt,wpar2, 9 wpar2bd,wpar2bt,walf,walfbd,walfbt,walf10,thx, & thxbd,thxbt,thy,thybd,thybt,thc,thcbd,thcbt c..for the iben high density thermal conductivity ov double precision ov,ovbd,ovbt, 1 d6,d613,d623,d623i,xne6,xne613,xne613i,xne623, 2 ef,efbd,efbt,em,embd,embt,zalphal,eg1,eg1bd, 3 eg1bt,eg2,eg2bd,eg2bt,eg,egbd,egbt,ef15,ef25, 4 ef35,gam,gambd,gambt c..for the transition between the low and high density thermal conductivities double precision farg,fargbd,fargbt,ffac,ffacbd,ffacbt, 1 drel,drel10,drelim,drelbt, 2 x,x1,x2,alfa,beta c..for the yakovlev high density thermal conductivity ov c..meff = hbar/(me*c) * (3.0 * pi * pi)**(1.0d0/3.0d0) c..weid = pi**2 * kerg**2 /(3.0 * me) c..iec = 4.0 * e**4 * me/(3.0 * pi * hbar**3) c..xec = hbar/kerg * e * sqrt(4.0 * pi/me) c..rt3 = sqrt(3) double precision xmas,ymas,wfac,cint,vie,cie,tpe,jy,vee,cee, 1 xmasbd,xmasbt,ymasbd,ymasbt,wfacbd,wfacbt, 2 viebd,viebt,ciebd,ciebt,tpebd,tpebt,xg,xgbd, 3 xgbt,jy1,jy1bd,jy1bt,jy2,jy2bd,jy2bt,jybd,jybt, 4 veebd,veebt,ceebd,ceebt,ov1,ov1bd,ov1bt, 5 xrel,xrelbd,xrelbt,beta2,beta2bd,beta2bt, 6 xrel2,xrel3,xrel4,jy3,jy3bd,jy3bt,jy4,jy4bd, 7 jy4bt,meff,weid,iec,xec,rt3 parameter (meff = 1.194648642401440e-10, 1 weid = 6.884326138694269e-5, 2 iec = 1.754582332329132e16, 3 xec = 4.309054377592449e-7, 4 rt3 = 1.7320508075688772d0) c..for easy toggling between the iben high density fits (whichk=1) c..and the yakovlev high density fits (whichk=2) integer whichk c parameter (whichk = 1) parameter (whichk = 2) c..for the total thermal conductivity, converted to an opacity double precision ocond,doconddd,doconddt,ocondmlt c..for the total thermal conductivity double precision scond,dsconddd,dsconddt c..for the total opacity and its derivatives double precision opac,dopacdd,dopacdt c..for the total thermal conductivity and its derivatives double precision sigma,dsigdd,dsigdt c..miscellaneous variables and constants c..con2 = con1 * sqrt(4.0*pi*qe*qe/me) c..k2c = conversion factor between opacities anc conductivities integer iz,i,j double precision t6,t6r,t6log,t6log10,t6inv,xdum,ydum,zdum, 1 third,twoth,twofive,sixfive,sixpi5,con1,con2,k2c parameter (third = 1.0d0/3.0d0, 1 twoth = 2.0d0/3.0d0, 2 twofive = 2.0d0/5.0d0, 3 sixfive = 6.0d0/5.0d0, 4 sixpi5 = pi*pi*pi*pi*pi/6.0d0, 5 con1 = hbar/(4.0d0 * kerg), 6 con2 = 1.07726359439811217d-7, 7 k2c = 4.0d0/3.0d0 * asol * clight) c..format statement for debugging 111 format(1x,1p3e14.6) c..set the multipliers oradmlt = 1.0d0 ocondmlt = 1.0d0 c..start of pipeline loop do j=jlo_kap,jhi_kap temp = temp_row(j) den = den_row(j) abar = abar_row(j) zbar = zbar_row(j) pep = pin_row(j) dpepdd = dpind_row(j) dpepdt = dpint_row(j) xne = xin_row(j) dxnedd = dxind_row(j) dxnedt = dxint_row(j) eta = ein_row(j) detadd = deind_row(j) detadt = deint_row(j) xh = xh_row(j) xhe = xhe_row(j) xz = xz_row(j) w4 = w4_row(j) w5 = w5_row(j) w6 = w6_row(j) zalpha = zalf_row(j) c..initialize opac = 0.0d0 dopacdd = 0.0d0 dopacdt = 0.0d0 orad = 0.0d0 doraddd = 0.0d0 doraddt = 0.0d0 oiben1 = 0.0d0 oiben1bd = 0.0d0 oiben1bt = 0.0d0 oiben2 = 0.0d0 oiben2bd = 0.0d0 oiben2bt = 0.0d0 ochrsty = 0.0d0 ochrstybd = 0.0d0 ochrstybt = 0.0d0 ocompt = 0.0d0 ocomptbd = 0.0d0 ocomptbt = 0.0d0 oh = 0.0d0 ohbd = 0.0d0 ohbt = 0.0d0 ov = 0.0d0 ovbd = 0.0d0 ovbt = 0.0d0 ocond = 0.0d0 doconddd = 0.0d0 doconddt = 0.0d0 c..some common factors zcz = w6 * third zbeta = w4 + w5 + w6 xmu = max(1.0d-99,zbeta-1.0d0) xmu12 = sqrt(xmu) xxmue = den*avo/xne xxmuebd = avo/xne - xxmue/xne * dxnedd xxmuebt = -xxmue/xne * dxnedt t6 = temp * 1.0d-6 t6r = sqrt(t6) t6log = log(t6) t6log10 = log10(t6) t6inv = 1.0d0/t6 c..iben's radiative opacity fit for small hydrogen (apj 196 525 1975) if (xh .lt. 1.0e-5) then xkc = (2.019e-4*den/t6**1.7d0)**(2.425d0) xkap = 1.0d0 + xkc * (1.0d0 + xkc * xcon1) xkb = 3.86d0 + 0.252d0*xmu12 + 0.018d0*xmu xka = (1.25d0 + 0.488d0*xmu12 + 0.092d0*xmu)*3.437d0 dbar = exp(-xka + xkb*t6log) oiben1 = xkap * (den/dbar)**xkrho xkcbd = 2.425d0 * xkc/den xkcbt = -4.1225d0 * xkc * t6inv xdum = 1.0d0 + 2.0d0 * xkc * xcon1 xkapbd = xkcbd * xdum xkapbt = xkcbt * xdum dbardt = dbar * xkb * t6inv oiben1bd = oiben1 * (xkapbd/xkap + xkrho/den) oiben1bt = oiben1 * (xkapbt/xkap - xkrho*dbardt/dbar) * 1.0d-6 end if c..iben's radiative opacity fit for large hydrogen, high temp if ( .not.((xh.ge.1.0e-5) .and. (t6.lt.t6_switch1)) .and. 1 .not.((xh.lt.1.0e-5) .and. (xz.gt.zbound)) ) then if (t6 .gt. 1.0) then d0log = -(3.868d0 + 0.806d0*xh) + 1.8d0*t6log else d0log = -(3.868d0 + 0.806d0*xh) + (3.42d0 - 0.52d0*xh)*t6log endif xka1 = 2.809d0 * exp(-(1.74d0 - 0.755d0*xh) 1 * (t6log10 - 0.22d0 + 0.1375d0*xh)**2) xkw = 4.05d0 * exp(-(0.306d0 - 0.04125d0*xh) 1 * (t6log10 - 0.18d0 + 0.1625d0*xh)**2) xkaz = xka1 * exp(-0.5206d0 * ((log(den)-d0log)/xkw)**2) 1 * xz/0.02d0 dbar2log = -(4.283d0 + 0.7196d0*xh) + 3.86d0*t6log dbar1log = -5.296d0 + 4.833d0*t6log if (dbar2log .lt. dbar1log) dbar1log = dbar2log dbar1 = exp(dbar1log) oiben2 = (den/dbar1)**xkrho * exp(xkaz) if (t6 .gt. t6_switch1) then d0logbt = 1.8d0 * t6inv else d0logbt = (3.42d0 - 0.52d0*xh) * t6inv endif xka1bt = -xka1 * 2.0d0 * (1.74d0 - 0.755d0*xh) 1 * (t6log10 - 0.22d0 + 0.1375d0*xh)/(t6*xln10) xkwbt = -xkw * 2.0d0 * (0.306d0 - 0.04125d0*xh) 1 * (t6log10 - 0.18d0 + 0.1625d0*xh)/(t6*xln10) xdum = (log(den)-d0log)/xkw xkazbd = -xkaz * 1.0412d0 * xdum/(den*xkw) xkazbt = xkaz * (1.0412d0 * xdum/xkw * (xdum*xkwbt + d0logbt) 1 + xka1bt/xka1) dbar2logbt = 3.86d0 * t6inv dbar1logbt = 4.833d0 * t6inv if (dbar2log .eq. dbar1log) dbar1logbt = dbar2logbt dbar1bt = dbar1 * dbar1logbt oiben2bd = oiben2 * (xkrho/den + xkazbd) oiben2bt = oiben2 * (xkazbt - xkrho/dbar1 * dbar1bt) * 1.0d-6 end if c..christy's radiative opacity for large hydro low temp (apj 144 108 1966) if ((t6.lt. t6_switch2) .and. (xh .ge. 1.0e-5)) then t4 = temp * 1.0d-4 t4r = sqrt(t4) t43 = t4 * t4 * t4 t44 = t43 * t4 t45 = t44 * t4 t46 = t45 * t4 ck1 = 2.0d6/t44 + 2.1d0*t46 ck3 = 4.0d-3/t44 + 2.0d-4/den**(0.25d0) ck2 = 4.5d0*t46 + 1.0d0/(t4*ck3) ck4 = 1.4d3*t4 + t46 ck5 = 1.0d6 + 0.1d0*t46 ck6 = 20.0d0*t4 + 5.0d0*t44 + t45 xkcx = xh*(t4r/ck1 + 1.0d0/ck2) xkcy = xhe*(1.0d0/ck4 + 1.5d0/ck5) xkcz = xz*(t4r/ck6) xdum = xkcx + xkcy + xkcz ochrsty = pep * xdum ck1bt = -8.0d6/t45 + 12.6d0*t45 ck3bd = -5.0d-5/den**(1.25d0) ck3bt = -16.0d-3/t45 ck2bd = -ck3bd/(t4*ck3**2) ck2bt = 27.0d0*t45 - (ck3 + t4*ck3bt)/(t4*ck3)**2 ck4bt = 1.4d3 + 6.0d0 * t45 ck5bt = 0.6d0 * t45 ck6bt = 20.0d0 + 20.0d0*t43 + 5.0d0*t44 xkcxbd = -xh/ck2**2 * ck2bd xkcxbt = xh*(0.5d0/(ck1*t4r) - t4r*ck1bt/ck1**2 - ck2bt/ck2**2) xkcybt = -xhe*(ck4bt/ck4**2 + 1.5d0*ck5bt/ck5**2) xkczbt = xz*(0.5d0/(t4r*ck6) - t4r*ck6bt/ck6**2) ochrstybd = dpepdd * xdum + pep * xkcxbd ochrstybt = dpepdt*xdum + pep*(xkcxbt+xkcybt+xkczbt) * 1.0d-4 end if c..calculate radiative opacity in presence of hydrogen if (xh .ge. 1.0e-5) then if (t6 .lt. t6_switch1) then orad = ochrsty doraddd = ochrstybd doraddt = ochrstybt else if (t6 .le. t6_switch2) then xdum = 1.5d0 - t6 ydum = t6 - 1.0d0 orad = 2.0d0*(ochrsty*xdum + oiben2*ydum) doraddd = 2.0d0*(ochrstybd*xdum + oiben2bd*ydum) doraddt = 2.0d0*(ochrstybt*xdum - ochrsty*1.0d-6 1 + oiben2bt*ydum + oiben2*1.0d-6) else orad = oiben2 doraddd = oiben2bd doraddt = oiben2bt end if c..calculate radiative opacity in the absence of hydrogen else if (xz .gt. zbound) then orad = oiben1 doraddd = oiben1bd doraddt = oiben1bt else dzfac = 1.0d0/zbound xdum = xz*dzfac ydum = (zbound-xz)*dzfac orad = oiben1*xdum + oiben2*ydum doraddd = oiben1bd*xdum + oiben2bd*ydum doraddt = oiben1bt*xdum + oiben2bt*ydum end if end if c..add in compton scattering opacity including effects of relativity and c..degeneracy. see buechler & yueh (1976), sampson (ap.j.129,734,1959) c..and chin(ap.j.142,1481,1965). this fit formula from wzw, valid from c..nondegenerate to partial degeneracy th = min(511.0d0, temp * 8.617d-8) fact = 1.0d0 + 2.75d-2*th - 4.88d-5*th*th facetax = 1.0d100 if (eta .le. 500.0) facetax = exp(0.522d0*eta - 1.563d0) faceta = 1.0d0 + facetax xdum = 1.0d0/(fact * faceta) ocompt = 6.65205d-25 * xdum * xne/den thbt = 8.617d-8 if (th .eq. 511.0) thbt = 0.0d0 factbt = thbt * (2.75d-2 - 9.76d-5 * th) facetaxbt = 0.0d0 facetaxbd = 0.0d0 if (eta .le. 500.0) then facetaxbt = facetax * 0.522d0 * detadt facetaxbd = facetax * 0.522d0 * detadd end if facetabt = facetaxbt facetabd = facetaxbd ydum = -xdum * (factbt * faceta + fact * facetabt) zdum = -xdum * fact * facetabd ocomptbt = ocompt * (ydum + dxnedt/xne) ocomptbd = ocompt * (zdum + dxnedd/xne - 1.0d0/den) orad = orad + ocompt doraddt = doraddt + ocomptbt doraddd = doraddd + ocomptbd c..cutoff radiative opacity when 4kt/hbar is less than the plasma frequency tcutoff = con2 * sqrt(xne) if (temp .lt. tcutoff) then if (tcutoff .gt. 200.0*temp) then orad = orad * 2.658d86 doraddd = doraddd * 2.658d86 doraddt = doraddt * 2.658d86 else cutfac = exp(tcutoff/temp - 1.0d0) cutfacbt = cutfac * tcutoff 1 * (0.5d0 * dxnedt/(xne*temp) - 1.0d0/temp**2) cutfacbd = cutfac * tcutoff * 0.5d0 * dxnedd/(xne*temp) xdum = orad orad = orad * cutfac doraddd = doraddd * cutfac + xdum * cutfacbd doraddt = doraddt * cutfac + xdum * cutfacbt end if end if c..fudge radiative opacity for low t; molecules and stuff zdum = orad t7 = temp * 1.0d-7 t73 = t7 * t7 * t7 t74 = t73 * t7 xkf = t7peek * den * t74 xdum = 1.0d0/(xkf + orad) orad = xkf * orad * xdum xkfbt = 4.0d0 * t7peek * den * t73 xkfbd = t7peek * t74 ydum = -xdum*xdum * (xkfbt + doraddt) doraddt = xkfbt*zdum*xdum + xkf*doraddt*xdum +xkf*zdum*ydum ydum = -xdum*xdum * (xkfbd + doraddd) doraddd = xkfbd*zdum*xdum + xkf*doraddd*xdum + xkf*zdum*ydum c..this section computes the thermal conductivity c..drel is the dividing line between nondegenerate and degenerate regions, c..taken from clayton eq. 2-34. if the density is larger than drel, then c..use the simpler degenerate expressions. if the density is smaller than c..drelim, use the ugly non degenerate formulas. in between drel and drelim, c..apply a smooth blending of the two. dlog10 = log10(den) drel = 2.4d-7 * zbar/abar * temp * sqrt(temp) if (temp .le. 1.0d5) drel = drel * 15.0d0 drelbt = 1.5d0 * drel/temp drel10 = log10(drel) c drelim = drel10 + 0.3 drelim = drel10 + 1.0d0 c..original iben limits c drel = 1.0d0 c drelbt = 0.0d0 c drel10 = 6.0d0 c drelim = 6.3d0 c drelim = 7.0d0 c..for the iben low density thermal conductivity c..iben's (apj 196 525 1975)fits to hubbard & lampe (apj s18,279,1969) c..note: iben's paper uses log10, a lot of the constants here c..were flipped to log, a legacy effect. log(x) = log(10) * log10(x) if (dlog10 .lt. drelim) then zdel = xne/(avo*t6*t6r) zdell10 = log10(zdel) eta0 = exp(-1.20322d0 + twoth * log(zdel)) eta02 = eta0*eta0 zdelbd = zdel/xne * dxnedd zdelbt = zdel * (dxnedt/xne - 1.5d-6/t6) zdell10bd = zdelbd/(zdel * xln10) zdell10bt = zdelbt/(zdel * xln10) eta0bd = eta0 * twoth/zdel * zdelbd eta0bt = eta0 * twoth/zdel * zdelbt eta02bd = 2.0d0 * eta0 * eta0bd eta02bt = 2.0d0 * eta0 * eta0bt c..thpl factor if (zdell10 .lt. 0.645) then xdum = zdel * (1.0d0 + 0.024417d0*zdel) thpl = -7.5668d0 + log(xdum) ydum = (1.0d0 + 0.048834d0*zdel)/xdum thplbd = zdelbd * ydum thplbt = zdelbt * ydum else if (zdell10 .lt. 2.5) then xdum = zdel*(1.0d0 + 0.02804d0*zdel) thpl = -7.58110d0 + log(xdum) ydum = (1.0d0 + 0.05608d0*zdel)/xdum thplbd = zdelbd * ydum thplbt = zdelbt * ydum if (zdell10 .ge. 2.0) then thpla = thpl thplabd = thplbd thplabt = thplbt ydum = 1.0d0 + 9.376d0/eta02 xdum = zdel * zdel * ydum thplb = -11.0742d0 + log(xdum) zdum = -zdel*9.376d0/eta02**2 thplbbd = zdel/xdum * (2.0d0*zdelbd*ydum + zdum * eta02bd) thplbbt = zdel/xdum * (2.0d0*zdelbt*ydum + zdum * eta02bt) cfac1 = 2.5d0 - zdell10 cfac1bd = -zdell10bd cfac1bt = -zdell10bt cfac2 = zdell10 - 2.0d0 cfac2bd = zdell10bd cfac2bt = zdell10bt thpl = 2.0d0 * (cfac1*thpla + cfac2*thplb) thplbd = 2.0d0 * (cfac1bd*thpla + cfac1*thplabd + 1 cfac2bd*thplb + cfac2*thplbbd) thplbt = 2.0d0 * (cfac1bt*thpla + cfac1*thplabt + 1 cfac2bt*thplb + cfac2*thplbbt) end if else ydum = 1.0d0 + 9.376d0/eta02 xdum = zdel * zdel * ydum thpl = -11.0742d0 + log(xdum) zdum = -zdel*9.376d0/eta02**2 thplbd = zdel/xdum * (2.0d0*zdelbd*ydum + zdum * eta02bd) thplbt = zdel/xdum * (2.0d0*zdelbt*ydum + zdum * eta02bt) end if end if c..pefac factor if (zdell10 .lt. 2.0) then pefac = 1.0d0 + 0.021876d0 * zdel pefacbd = 0.021876d0 * zdelbd pefacbt = 0.021876d0 * zdelbt if (zdell10 .gt. 1.5) then pefaca = pefac pefacabd = pefacbd pefacabt = pefacbt pefacal = log(pefaca) pefacalbd = pefacabd/pefaca pefacalbt = pefacabt/pefaca pefacb = 0.4d0 * eta0 + 1.64496d0/eta0 pefacbbd = 0.4d0 * eta0bd - 1.64496d0/eta02 * eta0bd pefacbbt = 0.4d0 * eta0bt - 1.64496d0/eta02 * eta0bt pefacbl = log(pefacb) pefacblbd = pefacbbd/pefacb pefacblbt = pefacbbt/pefacb cfac1 = 2.0d0 - zdell10 cfac1bd = -zdell10bd cfac1bt = -zdell10bt cfac2 = zdell10 - 1.5d0 cfac2bd = zdell10bd cfac2bt = zdell10bt pefac = exp(2.0d0 * (cfac1*pefacal + cfac2*pefacbl)) pefacbd = pefac * 2.0d0 * 1 (cfac1bd*pefacal + cfac1*pefacalbd + 2 cfac2bd*pefacbl + cfac2*pefacblbd) pefacbt = pefac * 2.0d0 * 1 (cfac1bt*pefacal + cfac1*pefacalbt + 2 cfac2bt*pefacbl + cfac2*pefacblbt) end if else pefac = 0.4d0 * eta0 + 1.64496d0/eta0 pefacbd = 0.4d0 * eta0bd - 1.64496d0/eta02 * eta0bd pefacbt = 0.4d0 * eta0bt - 1.64496d0/eta02 * eta0bt end if c..walf factor if (zdel.lt.40.0) then dnefac = 1.0d0 + zdel * (3.4838d-4 * zdel - 2.8966d-2) dnefacbd = zdelbd * (6.9676d-4 * zdel - 2.8966d-2) dnefacbt = zdelbt * (6.9676d-4 * zdel - 2.8966d-2) else xdum = 1.0d0 - 0.8225d0/eta02 ydum = 0.8225d0/eta02**2 dnefac = 1.5d0/eta0 * xdum dnefacbd = 1.5d0*(ydum/eta0 * eta02bd - xdum/eta02 * eta0bd) dnefacbt = 1.5d0*(ydum/eta0 * eta02bt - xdum/eta02 * eta0bt) endif xdum = xxmue*zbeta + dnefac zdum = t6r*pefac wpar2 = 9.24735e-3 * zdel * xdum/zdum ydum = xxmuebd*zbeta + dnefacbd wpar2bd = wpar2 * (zdelbd/zdel + ydum/xdum - t6r/zdum*pefacbd) ydum = xxmuebt*zbeta + dnefacbt wpar2bt = wpar2 * (zdelbt/zdel + ydum/xdum 1 - (t6r * pefacbt + pefac * 5.0d-7/t6r)/zdum) walf = 0.5d0 * log(wpar2) walfbd = 0.5d0/wpar2 * wpar2bd walfbt = 0.5d0/wpar2 * wpar2bt walf10 = 0.5d0 * log10(wpar2) c..thx factor if (walf10 .le. -3.0) then thx = exp(2.413d0 - 0.124d0*walf) thxbd = -0.124d0*thx*walfbd thxbt = -0.124d0*thx*walfbt else if (walf10 .le. -1.0) then thx = exp(0.299d0 - walf*(0.745d0 + 0.0456d0*walf)) thxbd = -thx*walfbd*(0.745d0 + 0.0912d0*walf) thxbt = -thx*walfbt*(0.745d0 + 0.0912d0*walf) else thx = exp(0.426d0 - 0.558d0*walf) thxbd = -0.558d0*thx*walfbd thxbt = -0.558d0*thx*walfbt end if c..thy factor if (walf10 .le. -3.0) then thy = exp(2.158d0 - 0.111d0*walf) thybd = -0.111d0*thy*walfbd thybt = -0.111d0*thy*walfbt else if (walf10 .le. 0.0) then thy = exp(0.553d0 - walf*(0.55d0 + 0.0299d0*walf)) thybd = -thy*walfbd*(0.55d0 + 0.0598d0*walf) thybt = -thy*walfbt*(0.55d0 + 0.0598d0*walf) else thy = exp(0.553d0 - 0.6d0*walf) thybd = -0.6d0*thy*walfbd thybt = -0.6d0*thy*walfbt end if c..thc factor if (walf10 .le. -2.5) then thc = exp(2.924d0 - 0.1d0*walf) thcbd = -0.1d0*thc*walfbd thcbt = -0.1d0*thc*walfbt else if (walf10 .le. 0.5) then thc = exp(1.6740d0 - walf*(0.511d0 + 0.0338d0*walf)) thcbd = -thc*walfbd*(0.511d0 + 0.0676d0*walf) thcbt = -thc*walfbt*(0.511d0 + 0.0676d0*walf) else thc = exp(1.941d0 - 0.785d0*walf) thcbd = -0.785d0*thc*walfbd thcbt = -0.785d0*thc*walfbt end if c..the non-degenerate, non-relativistic electron thermal opacity xdum = 1.0d0 / (t6*exp(thpl)) ydum = xh*thx + xhe*thy + zcz*thc oh = ydum * xdum ohbd = (xh*thxbd + xhe*thybd + zcz*thcbd)*xdum 1 - ydum*xdum*thplbd zdum = -xdum * xdum * exp(thpl) * (t6 * thplbt + 1.0d-6) ohbt = (xh*thxbt + xhe*thybt + zcz*thcbt)*xdum + ydum*zdum end if c..iben's fits to canuto for high density thermal conductivity ov if (dlog10 .gt. drel10) then if (whichk .eq. 1) then d6 = den * 1.0d-6 d613 = d6**third d623 = d6**twoth d623i = 1.0d0/d623 xdum = 1.0d-6/avo xne6 = xne * xdum xne613 = xne6**third xne613i = 1.0d0/xne613 xne623 = xne6**twoth ydum = sqrt(1.0d0 + xne623) ef = ydum - 1.0d0 zdum = 0.5d0/ydum * twoth * xne613i * xdum efbd = zdum * dxnedd efbt = zdum * dxnedt gam = 22.76d0 * d613 * zalpha * t6inv gambd = gam/d6 * third * 1.0d-6 gambt = -gam * t6inv * 1.0d-6 em = min(1.0d0, 0.5d0 + log10(ef)) embd = 0.0d0 embt = 0.0d0 if (em .lt. 1.0) then embd = efbd/(ef*xln10) embt = efbt/(ef*xln10) end if zalphal = log(zalpha) xdum = 0.767d0-0.168d0*zalphal eg1 = 2.010d0-0.298d0*zalphal + xdum * em eg1bd = xdum * embd eg1bt = xdum * embt eg2 = 1.0d0 - (1.0d0 + gam)**(-0.85d0) xdum = 0.85d0 * (1.0d0 + gam)**(-1.85d0) eg2bd = xdum * gambd eg2bt = xdum * gambt eg = exp(eg1*eg2) egbd = eg * (eg1bd*eg2 + eg1*eg2bd) egbt = eg * (eg1bt*eg2 + eg1*eg2bt) xdum = 6.753e-8 * t6 * t6 * zbeta ef15 = ef**(1.5d0) ef25 = ef15 * ef ef35 = ef25 * ef ydum = eg * (ef25 + ef35) ov = xdum/ydum ovbd = -ov/ydum * ( egbd * (ef25 + ef35) + 1 efbd * (2.5d0 * ef15 + 3.5d0 * ef25)) ovbt = -ov/ydum * ( egbt * (ef25 + ef35) + 1 efbt * (2.5d0 * ef15 + 3.5d0 * ef25)) 2 + 1.35060d-13 * t6 * zbeta/ydum c..here is the yakovlev and urpin (soviet astro 24 303 1980) c..fit to for a high density thermal conductivity ov else if (whichk .eq. 2) then xmas = meff * xne**third xmasbd = third * xmas/xne * dxnedd xmasbt = third * xmas/xne * dxnedt ymas = sqrt(1.0d0 + xmas*xmas) ymasbd = xmas/ymas * xmasbd ymasbt = xmas/ymas * xmasbt wfac = weid * temp/ymas * xne wfacbd = wfac * (dxnedd/xne - ymasbd/ymas) wfacbt = wfac * (dxnedt/xne - ymasbt/ymas + 1.0d0/temp) c..the yakovlev and urpin coulomb integral c cint = log(sqrt(1.5 + 3.0/plasg) * (2.0/3.0*pi*zbar)**third) - c 1 0.5e-4 * yakfac/(1.0d0 + 1.0e-4*yakfac) c c..the nandkumar and pethick (mnras 209 511 1984) coulomb integral c cint = sqrt(pi/3.0d0) * log(zbar**third) + c 1 twoth * log(1.32 + 2.33/sqrt(plasg)) - c 2 0.484e-4 * yakfac/(1.0d0 + yakc2*yakfac) c..nobodies coulomb integral cint = 1.0d0 c..the ion-electron collision frequency and the thermal conductivity xdum = iec * zbar * cint vie = xdum * ymas viebd = xdum * ymasbd viebt = xdum * ymasbt cie = wfac/vie ciebd = wfacbd/vie - cie/vie * viebd ciebt = wfacbt/vie - cie/vie * viebt c..start of the electron-electron collision frequency ydum = xne/ymas xdum = sqrt(ydum) tpe = xec * xdum tpebd = 0.5d0 * tpe/(ydum*ymas) * (dxnedd - xne/ymas * ymasbd) tpebt = 0.5d0 * tpe/(ydum*ymas) * (dxnedt - xne/ymas * ymasbt) xg = rt3 * tpe/temp xgbd = rt3 * tpebd/temp xgbt = (rt3*tpebt - xg)/temp c..from the function j(y) from potekhin, chabier & yakovlev a&a 1997 xdum = zbar/abar * den * 1.0d-6 xrel = 1.009d0 * xdum**third xrelbd = third * xrel/xdum * zbar/abar * 1.0d-6 xrelbt = 0.0d0 xdum = 1.0d0 + xrel**2 beta2 = xrel**2/xdum beta2bd = 2.0d0 * beta2 * xrelbd * (1.0d0/xrel - xrel/xdum) beta2bt = 0.0d0 xrel2 = xrel * xrel xrel3 = xrel2 * xrel xrel4 = xrel3 * xrel jy1 = 1.0d0 + sixfive/xrel2 + twofive/xrel4 jy1bd = -xrelbd/xrel3 * (2.0d0*sixfive + 4.0d0*twofive/xrel2) jy1bt = 0.0d0 xdum = 1.0d0 + 0.07414d0 * xg ydum = 3.0d0 * xdum**3 jy2 = xg**3 / ydum zdum = jy2 * (3.0d0/xg - 0.22242d0/xdum) jy2bd = xgbd * zdum jy2bt = xgbt * zdum xdum = 2.81d0 - 0.810d0*beta2 + xg ydum = xgbd - 0.810d0 * beta2bd zdum = xgbt - 0.810d0 * beta2bt jy3 = log(xdum/xg) jy3bd = ydum/xdum - xgbd/xg jy3bt = zdum/xdum - xgbt/xg xdum = 13.91d0 + xg ydum = xg/xdum jy4 = sixpi5 * ydum**4 zdum = 4.0d0 * jy4/ydum * (1.0d0 - xg/xdum)/xdum jy4bd = zdum * xgbd jy4bt = zdum * xgbt jy = jy1 * (jy2 * jy3 + jy4) jybd = jy1bd*jy2*jy3 + jy1*jy2bd*jy3 + jy1*jy2*jy3bd + 1 jy1bd*jy4 + jy1*jy4bd jybt = jy1bt*jy2*jy3 + jy1*jy2bt*jy3 + jy1*jy2*jy3bt + 1 jy1bt*jy4 + jy1*jy4bt c..finish of the electron-electron collicion frequency and the c..elecrton-electron thermal conductivity xdum = xmas/ymas ydum = 0.511d0 * temp**2 zdum = sqrt(xdum) vee = ydum * xdum/ymas * zdum * jy veebd = vee * (xmasbd/xmas - 2.0d0 * ymasbd/ymas 1 + 0.5d0/xdum * (xmasbd/ymas - xdum/ymas * ymasbd) 2 + jybd/jy) veebt = vee * (xmasbt/xmas - 2.0d0 * ymasbt/ymas 1 + 0.5d0/xdum * (xmasbt/ymas - xdum/ymas * ymasbt) 2 + jybt/jy + temp/ydum * 1.022d0) cee = wfac/vee ceebd = wfacbd/vee - cee/vee * veebd ceebt = wfacbt/vee - cee/vee * veebt c..sum of electron-ion and electron-electron thermal conductivity c..and conversion to an opacity xdum = cie * cee ydum = cee + cie ov1 = xdum/ydum ov1bd = (ciebd*cee + cie*ceebd - ov1 * (ceebd + ciebd))/ydum ov1bt = (ciebt*cee + cie*ceebt - ov1 * (ceebt + ciebt))/ydum ov = k2c/(ov1*den)*temp**3 ovbd = -ov * (1.0d0/den + ov1bd/ov1) ovbt = ov * (3.0d0/temp - ov1bt/ov1) end if end if c..interpolate between the low density and high density regimes if need be c..form from iben's 1975 paper if (dlog10 .le. drel10) then ocond = oh doconddd = ohbd doconddt = ohbt else if (dlog10 .gt. drel10 .and. dlog10 .lt. drelim) then x = den x1 = 10.0d0**drel10 x2 = 10.0d0**drelim alfa = (x-x2)/(x1-x2) beta = (x-x1)/(x2-x1) ocond = alfa*oh + beta*ov doconddd = alfa*ohbd + beta*ovbd doconddt = alfa*ohbt + beta*ovbt else if (dlog10 .ge. drelim) then ocond = ov doconddd = ovbd doconddt = ovbt end if c..the total opacity xdum = orad * ocond ydum = ocond + orad opac = xdum/ydum dopacdd = (doraddd*ocond + orad*doconddd)/ydum 1 - xdum/ydum**2 * (doconddd + doraddd) dopacdt = (doraddt*ocond + orad*doconddt)/ydum 1 - xdum/ydum**2 * (doconddt + doraddt) c..the total conductivities s2rad = k2c/(orad*den)*temp**3 ds2raddd = -s2rad * (1.0d0/den + doraddd/orad) ds2raddt = s2rad * (3.0d0/temp - doraddt/orad) scond = k2c/(ocond*den)*temp**3 dsconddd = -scond * (1.0d0/den + doconddd/ocond) dsconddt = scond * (3.0d0/temp - doconddt/ocond) sigma = k2c/(opac*den)*temp**3 dsigdd = -sigma * (1.0d0/den + dopacdd/opac) dsigdt = sigma * (3.0d0/temp - dopacdt/opac) c..end of pipeline loop, fill with the results orad_row(j) = orad dord_row(j) = doraddd dort_row(j) = doraddt ocond_row(j) = ocond docd_row(j) = doconddd doct_row(j) = doconddt opac_row(j) = opac dopd_row(j) = dopacdd dopt_row(j) = dopacdt s2rad_row(j) = s2rad dsrd_row(j) = ds2raddd dsrt_row(j) = ds2raddt scond_row(j) = scond dscd_row(j) = dsconddd dsct_row(j) = dsconddt sigma_row(j) = sigma dsigd_row(j) = dsigdd dsigt_row(j) = dsigdt enddo return end include 'public_timmes.f'