diff --git a/src/nulibtable.F90 b/src/nulibtable.F90 index 42122f9..ea7f42d 100644 --- a/src/nulibtable.F90 +++ b/src/nulibtable.F90 @@ -17,6 +17,8 @@ module nulibtable real*8, allocatable,save :: nulibtable_ewidths(:) real*8, allocatable,save :: nulibtable_ebottom(:) real*8, allocatable,save :: nulibtable_etop(:) + real*8, allocatable,save :: nulibtable_logenergies(:) + real*8, allocatable,save :: nulibtable_logetop(:) real*8, allocatable,save :: nulibtable_emissivities(:,:,:,:) real*8, allocatable,save :: nulibtable_absopacity(:,:,:,:) @@ -640,3 +642,331 @@ subroutine nulibtable_epannihil_single_species_range_energy(xtemp,xeta, & end subroutine nulibtable_epannihil_single_species_range_energy + +!!!! ROUTINES FOR COMPATIBILITY WITH GR1D !!!! + +!this takes temp,eta, and return phi0/1 over energy (both in and out) and species range +subroutine nulibtable_inelastic_range_species_range_energy2(xtemp,xeta, & + eas,eas_n1,eas_n2,eas_n3,eas_n4) + + use nulibtable + implicit none + + real*8, intent(in) :: xtemp, xeta !inputs + real*8 :: xltemp, xleta !log versions + integer, intent(in) :: eas_n1,eas_n2,eas_n3,eas_n4 + real*8, intent(out) :: eas(eas_n1,eas_n2,eas_n3,eas_n4) + integer :: ins,ing_in,ing_out + real*8 :: xeas(eas_n1*eas_n2*(eas_n2+1)/2) + integer :: index + real*8 :: energy_conversion = 1.60217733d-6*5.59424238d-55 + + + if(size(eas,1).ne.nulibtable_number_species) then + stop "nulibtable_inelastic_range_species_range_energy: supplied array dimensions (1) is not commensurate with table" + endif + if(size(eas,2).ne.nulibtable_number_groups) then + stop "nulibtable_inelastic_range_species_range_energy: supplied array dimensions (2) is not commensurate with table" + endif + if(size(eas,3).ne.nulibtable_number_groups) then + stop "nulibtable_inelastic_range_species_range_energy: supplied array dimensions (3) is not commensurate with table" + endif + if(size(eas,4).ne.2) then + stop "nulibtable_inelastic_range_species_range_energy: supplied array dimensions (4) is not commensurate with table" + endif + + xltemp = log10(xtemp) + xleta = log10(xeta) + + if (xltemp.lt.nulibtable_logItemp_min) stop "temp below nulib inelastic table minimum temp" + if (xltemp.gt.nulibtable_logItemp_max) stop "temp above nulib inelastic table maximum temp" + if (xleta.lt.nulibtable_logIeta_min) stop "eta below nulib inelastic table minimum eta" + if (xleta.gt.nulibtable_logIeta_max) stop "eta above nulib inelastic table maximum eta" + + xeas = 0.0d0 + call intp2d_many_mod(xltemp,xleta,xeas,nulibtable_Itable_Phi0,nulibtable_nItemp, & + nulibtable_nIeta,eas_n1*eas_n2*(eas_n2+1)/2,nulibtable_logItemp, & + nulibtable_logIeta) + + index = 0 + do ins=1,nulibtable_number_species + do ing_in=1,nulibtable_number_groups + do ing_out=1,ing_in + index = index + 1 + eas(ins,ing_in,ing_out,1) = 10.0d0**xeas(index) + enddo + enddo + do ing_in=1,nulibtable_number_groups + do ing_out=ing_in+1,nulibtable_number_groups + eas(ins,ing_in,ing_out,1) = & + exp(-(nulibtable_energies(ing_out)-nulibtable_energies(ing_in))/ & + (xtemp*energy_conversion))*eas(ins,ing_out,ing_in,1) !cernohorsky 94 + enddo + enddo + enddo + + xeas = 0.0d0 + call intp2d_many_mod(xltemp,xleta,xeas,nulibtable_Itable_Phi1,nulibtable_nItemp, & + nulibtable_nIeta,eas_n1*eas_n2*(eas_n2+1)/2,nulibtable_logItemp, & + nulibtable_logIeta) + + index = 0 + do ins=1,nulibtable_number_species + do ing_in=1,nulibtable_number_groups + do ing_out=1,ing_in + index = index + 1 + !this way because we interpolate \phi1/\phi0 + eas(ins,ing_in,ing_out,2) = xeas(index)*eas(ins,ing_in,ing_out,1) + enddo + enddo + do ing_in=1,nulibtable_number_groups + do ing_out=ing_in+1,nulibtable_number_groups + eas(ins,ing_in,ing_out,2) = & + exp(-(nulibtable_energies(ing_out)-nulibtable_energies(ing_in))/ & + (xtemp*energy_conversion))*eas(ins,ing_out,ing_in,2) !cernohorsky 94 + enddo + enddo + enddo + + + +end subroutine nulibtable_inelastic_range_species_range_energy2 + +!this takes temp,eta, and return phi0/1 over energy (both in and out) range for a single species +subroutine nulibtable_inelastic_single_species_range_energy2(xtemp,xeta, & + lns,eas,eas_n1,eas_n2,eas_n3) + + use nulibtable + implicit none + + real*8, intent(in) :: xtemp, xeta !inputs + real*8 :: xltemp, xleta !log versions + integer, intent(in) :: lns,eas_n1,eas_n2,eas_n3 + real*8, intent(out) :: eas(eas_n1,eas_n2,eas_n3) + integer :: ing_in,ing_out + real*8 :: xeas(eas_n1*(eas_n1+1)/2) + integer :: index + real*8 :: energy_conversion = 1.60217733d-6*5.59424238d-55 + integer :: startindex,endindex + + if(size(eas,1).ne.nulibtable_number_groups) then + stop "nulibtable_inelastic_single_species_range_energy: supplied array dimensions (1) is not commensurate with table" + endif + if(size(eas,2).ne.nulibtable_number_groups) then + stop "nulibtable_inelastic_single_species_range_energy: supplied array dimensions (2) is not commensurate with table" + endif + if(size(eas,3).ne.2) then + stop "nulibtable_inelastic_single_species_range_energy: supplied array dimensions (3) is not commensurate with table" + endif + + xltemp = log10(xtemp) + xleta = log10(xeta) + + if (xltemp.lt.nulibtable_logItemp_min) stop "temp below nulib inelastic table minimum temp" + if (xltemp.gt.nulibtable_logItemp_max) stop "temp above nulib inelastic table maximum temp" + if (xleta.lt.nulibtable_logIeta_min) stop "eta below nulib inelastic table minimum eta" + if (xleta.gt.nulibtable_logIeta_max) stop "eta above nulib inelastic table maximum eta" + + startindex = (lns-1)*nulibtable_number_groups*(nulibtable_number_groups+1)/2+1 + endindex = startindex + nulibtable_number_groups*(nulibtable_number_groups+1)/2 - 1 + + xeas = 0.0d0 + call intp2d_many_mod(xltemp,xleta,xeas,nulibtable_Itable_Phi0(:,:,startindex:endindex), & + nulibtable_nItemp,nulibtable_nIeta,eas_n1*(eas_n1+1)/2, & + nulibtable_logItemp,nulibtable_logIeta) + + index = 0 + do ing_in=1,nulibtable_number_groups + do ing_out=1,ing_in + index = index + 1 + eas(ing_in,ing_out,1) = 10.0d0**xeas(index) + enddo + enddo + do ing_in=1,nulibtable_number_groups + do ing_out=ing_in+1,nulibtable_number_groups + eas(ing_in,ing_out,1) = & + exp(-(nulibtable_energies(ing_out)-nulibtable_energies(ing_in))/ & + (xtemp*energy_conversion))*eas(ing_out,ing_in,1) !cernohorsky 94 + enddo + enddo + + xeas = 0.0d0 + call intp2d_many_mod(xltemp,xleta,xeas,nulibtable_Itable_Phi1(:,:,startindex:endindex), & + nulibtable_nItemp,nulibtable_nIeta,eas_n1*(eas_n1+1)/2, & + nulibtable_logItemp,nulibtable_logIeta) + + index = 0 + do ing_in=1,nulibtable_number_groups + do ing_out=1,ing_in + index = index + 1 + !this way because we interpolate \phi1/\phi0 + eas(ing_in,ing_out,2) = xeas(index)*eas(ing_in,ing_out,1) + enddo + enddo + do ing_in=1,nulibtable_number_groups + do ing_out=ing_in+1,nulibtable_number_groups + eas(ing_in,ing_out,2) = & + exp(-(nulibtable_energies(ing_out)-nulibtable_energies(ing_in))/ & + (xtemp*energy_conversion))*eas(ing_out,ing_in,2) !cernohorsky 94 + enddo + enddo + + + +end subroutine nulibtable_inelastic_single_species_range_energy2 + +!this takes temp,eta, and return phi0/1 for ep-annihilation over energy (both neutrinos) and species range +subroutine nulibtable_epannihil_range_species_range_energy2(xtemp,xeta, & + eas,eas_n1,eas_n2,eas_n3,eas_n4) + + use nulibtable + implicit none + + real*8, intent(in) :: xtemp, xeta !inputs + real*8 :: xltemp, xleta !log versions + integer, intent(in) :: eas_n1,eas_n2,eas_n3,eas_n4 + real*8, intent(out) :: eas(eas_n1,eas_n2,eas_n3,eas_n4) + integer :: ins,ing_this,ing_that + real*8 :: xeas(eas_n1*eas_n2*eas_n3*2) + integer :: index + real*8 :: energy_conversion = 1.60217733d-6*5.59424238d-55 + + + if(size(eas,1).ne.nulibtable_number_species) then + stop "nulibtable_epannihil_range_species_range_energy: supplied array dimensions (1) is not commensurate with table" + endif + if(size(eas,2).ne.nulibtable_number_groups) then + stop "nulibtable_epannihil_range_species_range_energy: supplied array dimensions (2) is not commensurate with table" + endif + if(size(eas,3).ne.nulibtable_number_groups) then + stop "nulibtable_enannihil_range_species_range_energy: supplied array dimensions (3) is not commensurate with table" + endif + if(size(eas,4).ne.4) then + stop "nulibtable_epannihil_range_species_range_energy: supplied array dimensions (4) is not commensurate with table" + endif + + xltemp = log10(xtemp) + xleta = log10(xeta) + + if (xltemp.lt.nulibtable_logItemp_min) stop "temp below nulib epannihil table minimum temp" + if (xltemp.gt.nulibtable_logItemp_max) stop "temp above nulib epannihil table maximum temp" + if (xleta.lt.nulibtable_logIeta_min) stop "eta below nulib epannihil table minimum eta" + if (xleta.gt.nulibtable_logIeta_max) stop "eta above nulib epannihil table maximum eta" + + xeas = 0.0d0 + call intp2d_many_mod(xltemp,xleta,xeas,nulibtable_epannihiltable_Phi0,nulibtable_nItemp, & + nulibtable_nIeta,eas_n1*eas_n2*eas_n2*2,nulibtable_logItemp, & + nulibtable_logIeta) + + index = 0 + do ins=1,nulibtable_number_species + do ing_this=1,nulibtable_number_groups + do ing_that=1,nulibtable_number_groups + index = index + 1 + eas(ins,ing_this,ing_that,1) = 10.0d0**xeas(index) + index = index + 1 + eas(ins,ing_this,ing_that,2) = 10.0d0**xeas(index) +! eas(ins,ing_this,ing_that,1) = eas(ins,ing_this,ing_that,2)* & +! exp(-(nulibtable_energies(ing_this)+nulibtable_energies(ing_that))/(xtemp*energy_conversion)) + + enddo + enddo + enddo + + xeas = 0.0d0 + call intp2d_many_mod(xltemp,xleta,xeas,nulibtable_epannihiltable_Phi1,nulibtable_nItemp, & + nulibtable_nIeta,eas_n1*eas_n2*eas_n2*2,nulibtable_logItemp, & + nulibtable_logIeta) + + index = 0 + do ins=1,nulibtable_number_species + do ing_this=1,nulibtable_number_groups + do ing_that=1,nulibtable_number_groups + !this way because we interpolate \phi1/\phi0 + index = index + 1 + eas(ins,ing_this,ing_that,3) = xeas(index)*eas(ins,ing_this,ing_that,1) + index = index + 1 + eas(ins,ing_this,ing_that,4) = xeas(index)*eas(ins,ing_this,ing_that,2) +! eas(ins,ing_this,ing_that,3) = eas(ins,ing_this,ing_that,4)* & +! exp(-(nulibtable_energies(ing_this)+nulibtable_energies(ing_that))/(xtemp*energy_conversion)) + enddo + enddo + enddo + + + +end subroutine nulibtable_epannihil_range_species_range_energy2 + +!this takes temp,eta, and return epannihilation phi0/1 over energy (both nus) range for a single species +subroutine nulibtable_epannihil_single_species_range_energy2(xtemp,xeta, & + lns,eas,eas_n1,eas_n2,eas_n3) + + use nulibtable + implicit none + + real*8, intent(in) :: xtemp, xeta !inputs + real*8 :: xltemp, xleta !log versions + integer, intent(in) :: lns,eas_n1,eas_n2,eas_n3 + real*8, intent(out) :: eas(eas_n1,eas_n2,eas_n3) + integer :: ing_this,ing_that + real*8 :: xeas(eas_n1*eas_n2*2) + integer :: index + real*8 :: energy_conversion = 1.60217733d-6*5.59424238d-55 + integer :: startindex,endindex + + if(size(eas,1).ne.nulibtable_number_groups) then + stop "nulibtable_epannihil_single_species_range_energy: supplied array dimensions (1) is not commensurate with table" + endif + if(size(eas,2).ne.nulibtable_number_groups) then + stop "nulibtable_epannihil_single_species_range_energy: supplied array dimensions (2) is not commensurate with table" + endif + if(size(eas,3).ne.4) then + stop "nulibtable_epannihil_single_species_range_energy: supplied array dimensions (3) is not commensurate with table" + endif + + xltemp = log10(xtemp) + xleta = log10(xeta) + + if (xltemp.lt.nulibtable_logItemp_min) stop "temp below nulib epannihil table minimum temp" + if (xltemp.gt.nulibtable_logItemp_max) stop "temp above nulib epannihil table maximum temp" + if (xleta.lt.nulibtable_logIeta_min) stop "eta below nulib epannihil table minimum eta" + if (xleta.gt.nulibtable_logIeta_max) stop "eta above nulib epannihil table maximum eta" + + startindex = (lns-1)*nulibtable_number_groups*nulibtable_number_groups*2+1 + endindex = startindex + nulibtable_number_groups*nulibtable_number_groups*2 - 1 + + xeas = 0.0d0 + call intp2d_many_mod(xltemp,xleta,xeas,nulibtable_epannihiltable_Phi0(:,:,startindex:endindex), & + nulibtable_nItemp,nulibtable_nIeta,eas_n1*eas_n2*2,nulibtable_logItemp,nulibtable_logIeta) + + index = 0 + do ing_this=1,nulibtable_number_groups + do ing_that=1,nulibtable_number_groups + index = index + 1 + eas(ing_this,ing_that,1) = 10.0d0**xeas(index) + index = index + 1 + eas(ing_this,ing_that,2) = 10.0d0**xeas(index) +! eas(ing_this,ing_that,1) = eas(ing_this,ing_that,2)* & +! (exp(-(nulibtable_energies(ing_this)+nulibtable_energies(ing_that))/(xtemp*energy_conversion))) + enddo + enddo + + xeas = 0.0d0 + call intp2d_many_mod(xltemp,xleta,xeas,nulibtable_epannihiltable_Phi1(:,:,startindex:endindex), & + nulibtable_nItemp,nulibtable_nIeta,eas_n1*eas_n2*2,nulibtable_logItemp,nulibtable_logIeta) + + index = 0 + do ing_this=1,nulibtable_number_groups + do ing_that=1,nulibtable_number_groups + index = index + 1 + !this way because we interpolate \phi1/\phi0 + eas(ing_this,ing_that,3) = xeas(index)*eas(ing_this,ing_that,1) + index = index + 1 + eas(ing_this,ing_that,4) = xeas(index)*eas(ing_this,ing_that,2) +! eas(ing_this,ing_that,3) = eas(ing_this,ing_that,4)*& +! (exp(-(nulibtable_energies(ing_this)+nulibtable_energies(ing_that))/(xtemp*energy_conversion))) + enddo + enddo + +end subroutine nulibtable_epannihil_single_species_range_energy2 + diff --git a/src/nulibtable_reader.F90 b/src/nulibtable_reader.F90 index 920a42a..b165c0e 100644 --- a/src/nulibtable_reader.F90 +++ b/src/nulibtable_reader.F90 @@ -85,7 +85,9 @@ subroutine nulibtable_reader(filename,include_Ielectron,include_epannihil_kernel allocate(nulibtable_ewidths(nulibtable_number_groups)) allocate(nulibtable_ebottom(nulibtable_number_groups)) allocate(nulibtable_etop(nulibtable_number_groups)) - + allocate(nulibtable_logenergies(nulibtable_number_groups)) + allocate(nulibtable_logetop(nulibtable_number_groups)) + call h5dopen_f(file_id, "nrho",dset_id, error) call h5dread_f(dset_id, H5T_NATIVE_INTEGER, nulibtable_nrho, dims1, error) call h5dclose_f(dset_id, error) diff --git a/src/scripts/interpolate_table.py b/src/scripts/interpolate_table.py index 9f5d21b..4202f11 100644 --- a/src/scripts/interpolate_table.py +++ b/src/scripts/interpolate_table.py @@ -25,29 +25,29 @@ def interp3d(x,y,z,xt,yt,zt,data): #------- determine location in (equidistant!!!) table - ix = 1 + int( (x - xt[0] - 1e-10) * dxi ) - iy = 1 + int( (y - yt[0] - 1e-10) * dyi ) - iz = 1 + int( (z - zt[0] - 1e-10) * dzi ) + ix = 1 + np.array( (x - xt[0] - 1e-10) * dxi ).astype(int) + iy = 1 + np.array( (y - yt[0] - 1e-10) * dyi ).astype(int) + iz = 1 + np.array( (z - zt[0] - 1e-10) * dzi ).astype(int) - ix = max( 1, min( ix, nx ) ) - iy = max( 1, min( iy, ny ) ) - iz = max( 1, min( iz, nz ) ) + ix = np.maximum( 1, np.minimum( ix, nx ) ) + iy = np.maximum( 1, np.minimum( iy, ny ) ) + iz = np.maximum( 1, np.minimum( iz, nz ) ) #------- set-up auxiliary arrays for Lagrange interpolation - delx = xt[ix] - x - dely = yt[iy] - y - delz = zt[iz] - z - - corners = np.zeros(8) - corners[0] = data[ix , iy , iz ] - corners[1] = data[ix-1, iy , iz ] - corners[2] = data[ix , iy-1, iz ] - corners[3] = data[ix , iy , iz-1] - corners[4] = data[ix-1, iy-1, iz ] - corners[5] = data[ix-1, iy , iz-1] - corners[6] = data[ix , iy-1, iz-1] - corners[7] = data[ix-1, iy-1, iz-1] + delx = xt[[ix]] - x + dely = yt[[iy]] - y + delz = zt[[iz]] - z + + corners = np.zeros((8,)+ix.shape) + corners[0] = data[[ix ], [iy ], [iz ]] + corners[1] = data[[ix-1], [iy ], [iz ]] + corners[2] = data[[ix ], [iy-1], [iz ]] + corners[3] = data[[ix ], [iy ], [iz-1]] + corners[4] = data[[ix-1], [iy-1], [iz ]] + corners[5] = data[[ix-1], [iy ], [iz-1]] + corners[6] = data[[ix ], [iy-1], [iz-1]] + corners[7] = data[[ix-1], [iy-1], [iz-1]] # coefficients @@ -87,22 +87,22 @@ def interp2d(x, y, xt, yt, data): #------- determine location in (equidistant!!!) table - ix = 1 + int( (x - xt[0] - 1e-10) * dxi ) - iy = 1 + int( (y - yt[0] - 1e-10) * dyi ) + ix = 1 + np.array( (x - xt[0] - 1e-10) * dxi ).astype(int) + iy = 1 + np.array( (y - yt[0] - 1e-10) * dyi ).astype(int) - ix = max( 1, min( ix, nx ) ) - iy = max( 1, min( iy, ny ) ) + ix = np.maximum( 1, np.minimum( ix, nx ) ) + iy = np.maximum( 1, np.minimum( iy, ny ) ) #------- set-up auxiliary arrays for Lagrange interpolation - delx = xt[ix] - x - dely = yt[iy] - y + delx = xt[[ix]] - x + dely = yt[[iy]] - y - corners = np.zeros(4) - corners[0] = data[ix , iy ] - corners[1] = data[ix-1, iy ] - corners[2] = data[ix , iy-1] - corners[3] = data[ix-1, iy-1] + corners = np.zeros((4,)+ix.shape) + corners[0] = data[[ix ], [iy ]] + corners[1] = data[[ix-1], [iy ]] + corners[2] = data[[ix ], [iy-1]] + corners[3] = data[[ix-1], [iy-1]] #------ set up coefficients of the interpolation polynomial and @@ -132,7 +132,7 @@ def interpolate_eas(ig, s, rho, T, Ye, table, datasetname): #! table h5py.File() #! datasetname {"absorption_opacity", "emissivities", "scattering_opacity", "scattering_delta"} #!--------------------------------------------------------------------- - data = table[datasetname] + data = np.array(table[datasetname]) zt = np.log10(table["rho_points"]) yt = np.log10(table["temp_points"]) @@ -168,7 +168,7 @@ def interpolate_kernel(ig_in, s, ig_out, eta, T, table, datasetname): #! table h5py.File() #! datasetname {"inelastic_phi0", "inelastic_phi1"} #!--------------------------------------------------------------------- - data = table[datasetname][ig_in,s,ig_out,:,:] + data = np.array(table[datasetname][ig_in,s,ig_out,:,:]) xt = np.log10(table["eta_Ipoints"]) yt = np.log10(table["temp_Ipoints"]) @@ -184,11 +184,11 @@ def interpolate_kernel(ig_in, s, ig_out, eta, T, table, datasetname): return interp2d(x, y, xt, yt, data) def interpolate_eos(rho, T, Ye, table, datasetname): - data = table[datasetname] + data = np.array(table[datasetname]) - zt = table["logrho"] - yt = table["logtemp"] - xt = table["ye"] + zt = np.array(table["logrho"]) + yt = np.array(table["logtemp"]) + xt = np.array(table["ye"]) z = np.log10(rho) y = np.log10(T)