diff --git a/src/Interpolated_data.h5 b/src/Interpolated_data.h5 new file mode 100644 index 0000000..a26ce55 Binary files /dev/null and b/src/Interpolated_data.h5 differ diff --git a/src/Makefile b/src/Makefile index 33413e9..f2173ba 100644 --- a/src/Makefile +++ b/src/Makefile @@ -13,7 +13,8 @@ SOURCES = fermi.F90 \ scattering.F90 \ weak_magnetism_correction.F90 \ inelastic_electron_scattering.F90 \ - eos_interface.F90 + eos_interface.F90\ + bremsstrahlung.F90 NT_SOURCES=nulibtable.F90 \ nulibtable_reader.F90 \ @@ -92,9 +93,9 @@ EXTRAINCS += $(HDF5INCS) EXTRAOBJECTS += $(HDF5LIBS) ifeq ($(WEAK_RATES),1) -all: nulib.a ../nulibtable_driver ../point_example ../make_table_weakrates +all: nulib.a ../nulibtable_driver ../point_example ../bremsstrahlung_point_example ../make_table_weakrates else -all: nulib.a ../nulibtable_driver ../point_example ../make_table_example +all: nulib.a ../nulibtable_driver ../point_example ../bremsstrahlung_point_example ../make_table_example endif ../nulibtable_driver: $(EXTRADEPS) $(NT_OBJECTS) $(F_OBJECTS) $(f90_OBJECTS) @@ -103,6 +104,9 @@ endif ../point_example: $(EXTRADEPS) $(F_OBJECTS) $(MOD_OBJECTS) $(OBJECTS) $(f90_OBJECTS)point_example.F90 $(F90) $(F90FLAGS) $(DEFS) $(MODINC) $(EXTRAINCS) -o $@ point_example.F90 $(MOD_OBJECTS) $(OBJECTS) $(F_OBJECTS) $(EXTRAOBJECTS) $(f90_OBJECTS) +../bremsstrahlung_point_example: $(EXTRADEPS) $(F_OBJECTS) $(MOD_OBJECTS) $(OBJECTS) $(f90_OBJECTS)bremsstrahlung_point_example.F90 + $(F90) $(F90FLAGS) $(DEFS) $(MODINC) $(EXTRAINCS) -o $@ bremsstrahlung_point_example.F90 $(MOD_OBJECTS) $(OBJECTS) $(F_OBJECTS) $(EXTRAOBJECTS) $(f90_OBJECTS) + ../make_table_example: $(EXTRADEPS) $(F_OBJECTS) $(MOD_OBJECTS) $(OBJECTS) $(f90_OBJECTS)make_table_example.F90 $(F90) $(F90FLAGS) $(DEFS) $(MODINC) $(EXTRAINCS) -o $@ make_table_example.F90 $(MOD_OBJECTS) $(OBJECTS) $(F_OBJECTS) $(EXTRAOBJECTS) $(f90_OBJECTS) @@ -138,6 +142,7 @@ clean: rm -rf ../make_table_example rm -rf ../make_table_weakrates rm -rf ../point_example + rm -rf ../bremsstrahlung_point_example rm -rf ../nulibtable_driver rm -rf test rm -rf *.o diff --git a/src/bremsstrahlung.F90 b/src/bremsstrahlung.F90 new file mode 100644 index 0000000..f6a9c1f --- /dev/null +++ b/src/bremsstrahlung.F90 @@ -0,0 +1,411 @@ +!-*-f90-*- + + + +function Bremsstrahlung_Phi0_Hannestad(nu_energy_x,nubar_energy_x& + ,matter_temperature,n_N,neutrino_species,pro_ann) result(Phi0_nn) + +! Rates from Hannestad & Raffelt 1998 +use nulib +implicit none +real*8,intent(in) :: nu_energy_x,nubar_energy_x ! dimensionless +real*8, intent(in) :: matter_temperature ! MeV +real*8, intent(in) :: n_N ! n/cm3 +integer, intent(in) :: neutrino_species ! 1 to 6 +integer,intent(in) :: pro_ann ! 0 or 1 + +real*8 :: eta_star +real*8 :: Gamma_sigma +real*8 :: S_sig +real*8 :: p_f +real*8 :: s_func +real*8 :: g +real*8 :: x +real*8 :: y +real*8 :: G_f +real*8 :: G_rampp +real*8 :: GC2_raffelt +real*8,parameter :: C_a = gA/2.0d0 +!~ real*8,parameter :: C_a = 0.5d0 + + +!functions +real*8 :: find_g +real*8 :: find_s +real*8 :: find_s2 + +!output +real*8 :: Phi0_nn +real*8 :: Phi0_nn_old +real*8 :: Phi0_nn_Raff + + + +!Rampp +p_f = hbarc_mevcm/clight * (3.0d0 * pi **2 * n_N)**(1.0d0/3.0d0) ! Mev s cm-1 + +eta_star = (p_f)**2 /(2.0d0 *m_amu/clight**2 * matter_temperature) !adimensional + + + +Gamma_sigma = 8.0d0 * SQRT(2.0d0*pi) * 15.0d0**2 / (3.0d0*pi**2) * & + eta_star**(1.5d0) * & + matter_temperature**2/m_amu ! MeV !alpha_pi = 15 + + +G_f = Gfermi*(hbarc_mevcm) ! Mev-1 cm +G_rampp = (G_f*hbarc_mevcm**2)**2/(pi*(hbarc_mevcm)**4) !cm2 MeV^-2 +GC2_raffelt = 2.1d0 * 1d-44 !cm2 Mev-2 + +!Raffelt + +!~ y = 1.94d0* 10.0d0/matter_temperature ! adim +y = m_pi**2/(matter_temperature*m_amu) !adim +x = nu_energy_x + nubar_energy_x ! adimensionnal +s_func = find_s(eta_star,y,x) ! adimenssional +g = find_g(eta_star,y) ! adimenssional + + + +!Phi 0 + +S_sig = Gamma_sigma / & + ( ( x *matter_temperature)**2 + & + 0.25d0 * ( Gamma_sigma * g) **2 ) & + * s_func ! MeV^-1 + + +Phi0_nn = 3.0d0*G_rampp * C_a**2 * n_N * S_sig & + *hbarc_mevcm**3 *clight/(2.0d0*pi)**2 ! cm3 s-1 + +Phi0_nn_old = 6.0d0*G_f ** 2 * C_a**2 * n_N * S_sig & + *(hbarc_mevcm)**3 *clight ! cm3 s-1 + + +Phi0_nn_Raff = GC2_raffelt * n_N * S_sig & + *hbarc_mevcm**3 *clight/(2.0d0*pi) ! cm3 s-1 + + +if (pro_ann .EQ. 1) then + Phi0_nn = Phi0_nn_old ! absorption +else if (pro_ann .EQ. 0) then + Phi0_nn = Phi0_nn_old * exp(-x) ! production +else + write(*,*) "brem : annihilation or production wrongly inserted" +end if + + +if ( Phi0_nn .NE. Phi0_nn) then + write(*,*) "NaNing", s_func,x,y,matter_temperature,eta_star + stop +endif +if ( Phi0_nn .LT. 0.0d0) then + write(*,*) "negative Han brem", s_func,x,y,matter_temperature,eta_star + stop +endif + +!~ write(*,*) S_sig +end function Bremsstrahlung_Phi0_Hannestad + + + +function Bremsstrahlung_Phi0_Kuroda(nu_energy_x,nubar_energy_x& + ,matter_temperature,n_N,neutrino_species,pro_ann) result(Phi0_nn) + +! Kuroda rates from Kuroda et al. 2016 ( Raffelt 1996) modified with factor from Raffelt & Seckel 1995 +use nulib +implicit none +real*8,intent(in) :: nu_energy_x,nubar_energy_x ! dimensionless +real*8, intent(in) :: matter_temperature ! MeV +real*8, intent(in) :: n_N ! n/cm3 +integer, intent(in) :: neutrino_species +integer,intent(in) :: pro_ann ! 0 or 1 + +real*8 :: p_f +real*8 :: x +real*8 :: Phi_ND +real*8 :: Phi_ND_2 +real*8 :: Phi_D +real*8 :: Phi_D_2 +real*8 :: constante +real*8,parameter :: G = 1.55d-33 ! G**2 in Kuroda Mev-2 cm3 s-1 +real*8,parameter :: C_A = gA/2.0d0 + + +!output +real*8 :: Phi0_nn + +!function +real*8 :: find_s2 + +x = nubar_energy_x + nu_energy_x ! dimensionless +p_f = hbarc_mevcm/clight * ( 3.0d0 * pi**2 * n_N)**(1.0d0/3.0d0) ! cm2 s-1 + + +Phi_D_2 = G / 4.0d0 * 15.0d0**2*(4.0d0*pi)**2 * pi**4 & + / (12.0d0 * pi**9) * matter_temperature/x * (4.0d0 * pi**2 + x**2) & + / ( 1.0d0 - exp(-x)) * 3.0d0* C_A**2.0d0 * 2.0d0 * hbarc_mevcm & + * ( 3.0d0 *pi**2 * n_N) ** (1.0d0/3.0d0) + + +Phi_ND_2 = hbarc_mevcm**6 * G/4.0d0 * 15.0d0**2 *(4.0d0*pi)**2/m_amu**4& + *pi**4 * 2.0d0 * m_amu**(1.5d0) & + / pi**(5.5d0) /( matter_temperature**(1.5d0) *x**2) & + *3.0d0 * C_A**2.0d0 * n_N ** 2 * find_s2(0.0d0,0.0d0,x)! using the s_kl part as well + + +if ( pro_ann .EQ. 0 ) then ! pro + Phi0_nn = exp(-x)* MIN(Phi_D_2,Phi_ND_2) !cm3 s-1 +else if (pro_ann .EQ. 1 ) then !abs + Phi0_nn = MIN(Phi_D_2,Phi_ND_2) ! cm3 s-1 +else + stop " pro_ann not allowed" +endif + + + +end function Bremsstrahlung_Phi0_Kuroda + + + +function find_g(eta_star,y) result (g) + implicit none + real*8, intent(in) :: eta_star + real*8, intent(in) :: y + + !output + real*8 :: g + + real*8 :: alpha1 + real*8 :: alpha2 + real*8 :: alpha3 + real*8 :: p1 + real*8 :: p2 + + alpha1 = (0.5d0 + 1.0d0/eta_star)/ & + ((1.0d0 + 1.0d0/eta_star) * (25.0d0 * y**2 + 1.0d0)) + & + (0.5d0 + eta_star/15.6d0)* (25.0d0 * y**2 )/(25.0d0*y**2 + 1.0d0) + + + alpha2 = (0.63d0 + 0.04d0 * eta_star**(1.45d0) )/& + (1.0d0 + 0.02d0 * eta_star**(2.5d0)) + + alpha3 = 1.2d0 * exp(0.6d0*eta_star - 0.4d0 * eta_star**(1.5d0)) + + p1= (1.8d0 + 0.45d0*eta_star) / ( 1.0d0 + 0.15d0 *eta_star**(1.5d0)) + + p2 = 2.3d0 - 0.05d0*eta_star/(1.0d0 + 0.025*eta_star) + + g = (alpha1 + alpha2 *y **(p1))/ & + (1.0d0 + alpha3 * y **(p2) + alpha2 * y**(p1+2.0d0)/13.75d0) + + if (g .LT. 0.0d0 ) stop "brem : g < 0" + +end function find_g + + + +function find_s(eta, y,x) result(s) +use nulib +implicit none + +real*8, intent(in) :: x +real*8,intent(in) :: y +real*8,intent(in) :: eta + +real*8 :: s_ND +real*8 :: s_D +real*8 :: u +real*8 :: f +real*8 :: h +real*8 :: p +real*8 :: C +real*8 :: G +real*8 :: big_F + +!result +real*8 :: s + + +s_ND = 2.0d0 * sqrt(pi) * ( x + 2.0d0 - exp(-y/12.0d0))**1.5d0 * & + (x**2 + 2.0d0*x*y + 5.0d0/3.0d0 * y**2 +1.0d0) / ( sqrt(pi) & + + ( pi **(1.0d0/8.0d0) + x + y)**4) + + +u = sqrt(y/(2.0d0*eta)) + +f = 1.0d0 - 5.0d0*u/6.0d0*atan(2.0d0/u) + u**2 /(3.0d0*(u**2 +4.0d0)) & + + u**2 /(6.0d0*sqrt(2.*u**2 +4.0d0)) * atan( 2.0d0*sqrt(2.0d0*u**2+4.0d0)/u**2) + +s_D =3.0d0*(pi/2.0d0)**2.5d0 * eta**(-2.5d0)* (x**2 + 4.0d0*pi**2)*x & + /(4.0d0*pi**2 *(1.0d0 - exp(-x))) * f + + + +h = 0.1d0 * eta /(2.39d0 + 0.1d0 * eta**1.1) + +p = 0.67d0 + 0.18d0*y**0.4d0 + +C = 1.1d0* x**1.1d0 * h /(2.3 +h * x**0.93d0+0.0001d0*x**1.2)*30.0d0 & + /(30.0d0+0.005d0*x**2.8d0) + +G = 1.0d0 - 0.0044*x**1.1 * y/(0.8d0+0.06d0*y*1.05d0)* sqrt(eta)/(eta+0.2d0) + +big_F = 1.0d0 + 1.0d0/(( 3.0d0 +( x-1.2d0)**2 + x**(-4.0d0)) * & + (1.0d0+ eta**2.0d0)*(1+y**4.0d0)) + +s = ( s_ND**(-p) + ABS(s_D)**(-p)) ** (-1.0d0/p) * big_F *(1.00d0+ C*G) ! adim + +if ( s.NE.s) then + write(*,*) "find_s NaNing",s_ND,s_D, big_F,C,G,eta,y,x + stop +endif + +end function find_s + +function find_s2(eta, y,x) result(s) +use nulib +implicit none + +real*8, intent(in) :: x +real*8,intent(in) :: y +real*8,intent(in) :: eta + +real*8 :: s_ND +real*8 :: s_D +real*8 :: u +real*8 :: f +real*8 :: h +real*8 :: p +real*8 :: C +real*8 :: G +real*8 :: big_F + + +!result +real*8 :: s + + +s_ND = sqrt(1.0d0+x*pi/4.0d0) - (1.0d0+(pi*x/4.0d0)**(5.0d0/4.0d0))**0.4d0& + + (1+ (64.0d0/(169.0d0*pi)*x)**(5.0d0/4.0d0))**(-0.4d0) ! (s_0 - s_k.l) from B.12 & B.16 in Raffelt & Seckel 1995 + + + + +s_D = (x**2 + 4.0d0 * pi**2 ) * x/(1.0d0-exp(-abs(x)))/(4.0d0*pi**2) & + *3.0d0*(pi/2.0d0)**2.5d0/(eta**2.5d0) + +s= s_ND ! used in Kuroda rates + + +if ( s.NE.s) then + write(*,*) "find_s2 NaNing",s_ND,s_D, big_F,C,G,eta,y,x + stop +endif + +end function find_s2 + +function Bremsstrahlung_Phi0_gang(nu_energy_x,nubar_energy_x& + ,matter_temperature,n_N,Ye,neutrino_species,pro_ann) result(Phi0_nn) +use nulib + +implicit none +real*8,intent(in) :: nu_energy_x,nubar_energy_x ! dimensionless +real*8, intent(in) :: matter_temperature ! MeV +real*8, intent(in) :: n_N ! n/cm3 +real*8, intent(in) :: Ye ! n/cm3 +integer, intent(in) :: neutrino_species ! 1 to 6 +integer,intent(in) :: pro_ann ! 0 or 1 + + +real*8 :: Phi0_nn + +! internal variables +integer :: in_N,i +real*8 :: Itable_n_N(25) +real*8 :: temp_array(25) = (/(i, i=2,50, 2)/) +real*8 :: Ye_array(26) = (/(i, i=0,50,2)/)/100.0d0 +real*8 :: om_array(40) = (/(i,i = 1,40,1)/)/10.0d0 -1.4d0 +real*8 :: omega +real*8 :: S,n_N_fm +real*8 :: G_f +real*8,parameter :: C_a = gA/2.0d0 +real*8 :: eas_brem(25) +real*8 :: dx,delt +integer :: indx +real*8 :: t_try + + do in_N=1,25 + Itable_n_N(in_N) = & + (-4.0d0+dble(in_N-1)/dble(24)*(4.0d0)) + enddo + +n_N_fm= log10(n_N*(1.0d-13)**3) + +omega= log10((nu_energy_x+nubar_energy_x)*matter_temperature) +S=0.0d0 +indx=0 + +if ( matter_temperature .GT. maxval(temp_array) .or. matter_temperature .LT. minval(temp_array) & + .or. Ye .GT. maxval(Ye_array) .or. Ye .LT. minval(Ye_array) & + .or. abs(n_N_fm) .GE. maxval(abs(Itable_n_N)) .or. abs(n_N_fm) .LE. minval(abs(Itable_n_N)) & + .or. omega .GT. maxval(om_array) .or. omega .LT. minval(om_array) ) then + +else + + call intp3d_many (omega, Ye,n_N_fm, eas_brem, 1, gang_table,40,26,25,25& + ,om_array, Ye_array,Itable_n_N ) + + + + dx = 0.5d0 + indx= 1+ INT((matter_temperature-temp_array(1))*dx) + + + delt = (matter_temperature-temp_array(indx)) / 2.0d0 + + S = eas_brem(indx) + (eas_brem(indx+1)-eas_brem(indx))*delt + + + +endif + + + + +!Rampp +G_f = Gfermi*(hbarc_mevcm) ! Mev-1 cm + +Phi0_nn = 6.0d0*G_f ** 2 * C_a**2 * n_N * S & + *(hbarc_mevcm)**3 *clight ! cm3 s-1 + +if (pro_ann .EQ. 1) then + Phi0_nn = Phi0_nn ! absorption +else if (pro_ann .EQ. 0) then + Phi0_nn = Phi0_nn * exp(-(nu_energy_x+nubar_energy_x)) ! production +else + write(*,*) "brem : annihilation or production wrongly inserted" +end if + + + +if (Phi0_nn .LT. 0.0d0) then + write(*,*) " brem negative : problem" + write(*,*) "S",S,"Phi0_nn :", Phi0_nn,"Ye :", Ye, "omega :",omega,"n :"& + ,n_N_fm,"n_cm : ",n_N,"T :", matter_temperature,"pro_ann :",pro_ann + write(*,*) eas_brem + write(*,*) "arg", shape(eas_brem), shape(gang_table) + write(*,*) + write(*,*) "nu,nubar,temp",nu_energy_x,nubar_energy_x& + ,matter_temperature + write(*,*) neutrino_species,pro_ann + write(*,*) + write(*,*) om_array + write(*,*) temp_array + write(*,*) Ye_array + write(*,*) Itable_n_N + stop +endif + + +end function Bremsstrahlung_Phi0_gang diff --git a/src/bremsstrahlung_point_example.F90 b/src/bremsstrahlung_point_example.F90 new file mode 100644 index 0000000..6d10587 --- /dev/null +++ b/src/bremsstrahlung_point_example.F90 @@ -0,0 +1,658 @@ +!-*-f90-*- +program bremsstrahlung_point_example + + use nulib + use inputparser + use eosmodule, ONLY : energy_shift +#if NUCLEI_HEMPEL + use nuclei_hempel, only : set_up_Hempel +#endif + implicit none + + !many people use different number of species, this is to denote how they are devided up. + ! mytable_neutrino_scheme = 1 (three output species) + ! species #1: electron neutrino #2 electron antineutrino + ! #3: muon+tau neutrino+antineutrino + ! neutrino_scheme = 2 (four output species) + ! species #1: electron neutrino #2 electron antineutrino + ! #3: muon+tau neutrino #4 mu and tau antineutrino + ! neutrino_scheme = 3 (six output species) + ! species #1: electron neutrino #2 electron antineutrino + ! #3: muon neutrino #4 mu antineutrino + ! #5: tau neutrino #6 tau antineutrino + integer :: mypoint_neutrino_scheme = 1 + + !number of species for nulib to calculate interactions for, must + !be six currently, average is done via above parameter + integer :: mypoint_number_species = 6 + + !number of species for nulib to output interactions for, must + !be commensurate with neutrino_scheme above + integer :: mypoint_number_output_species = 3 + + !number of energy groups + integer :: mypoint_number_groups = 50 !18 + + character*200 :: parameters = "./parameters" + + !local variables + real*8, allocatable,dimension(:,:) :: local_emissivity + real*8, allocatable,dimension(:,:) :: local_absopacity + real*8, allocatable,dimension(:,:) :: local_scatopacity + real*8, allocatable,dimension(:,:) :: local_delta + real*8, allocatable,dimension(:,:) :: local_Phi0, local_Phi1 + real*8, allocatable,dimension(:,:,:) :: local_Phi0_epannihil, local_Phi1_epannihil + real*8, allocatable,dimension(:,:,:) :: local_Phi0_bremsstrahlung + real*8, allocatable,dimension(:,:,:) :: local_Phi0_bremsstrahlung2 + real*8, allocatable,dimension(:,:) :: blackbody_spectra + real*8, allocatable,dimension(:) :: eos_variables + real*8,allocatable, dimension(:) :: Phi0_brem,Phi0_brem2 + real*8,allocatable, dimension(:) :: Phi0_brem_ann,Phi0_brem2_ann + real*8 :: matter_prs,matter_ent,matter_cs2,matter_dedt,matter_dpderho,matter_dpdrhoe + real*8 :: single_neutrino_emissivity_from_NNBrem_given_energypoint + real*8 :: single_neutrino_emissivity_from_epannihil_given_energypoint + integer :: keytemp,keyerr + real*8 :: precision = 1.0d-10 + real*8 :: xrho, xtemp, xye + real*8 :: n_N,eta_star + integer :: i,j,inde,in_N,zone + real*8 :: dxfac,mindx + real*8 :: Q,Q_ann,Q2,Q2_ann,M,dE_dEdt,dE_dEdt2,dE_dEdt3,dE_dEdt1_ann,dE_dEdt2_ann,dE_dEdt3_ann,k,N,M2 + real*8 :: Kuro_inv_lamb,Hann_inv_lamb,epan_inv_lamb + real*8 :: dn_dEdt,dn_dEdt2,dn_dEdt3,ratio,ratio_2,dn_dEdt_burrows + real*8 :: dE_dEdt_2,dE_dEdt2_2,dE_dEdt3_2,dE_dEdt_burrows,dn_dEdt_bruenn + real*8 :: Bremsstrahlung_Phi0_Hannestad + real*8 :: Bremsstrahlung_Phi0_Kuroda + real*8 :: Bremsstrahlung_Phi0_gang + real*8 :: epannihil_Phi_Bruenn + real*8 :: S_sigma_static,Itable_n_N(25) + real*8 :: om_array(40) = (/(i,i = 1,40,1)/)/10.0d0 -1.4d0 + real*8 :: fermi +!~ real*8 :: find_s + real*8,dimension(100) :: integral1,integral2,integral3,temp_array,dens_array,ener_array + real*8:: integral + real*8 :: J_1,J_1_bar,phi_a,phi_p,energy,energy_2 + real*8,dimension(100) :: M1_mom,M1_mom_inv + real*8 :: ye_array(13) = (/(i,i = 43,559,43)/)/1000.0d0 + real*8 :: rho_array(143) = (/(i,i = 6,148,1)/)/10.0d0 + real*8,dimension(4) :: rho_figure +!~ real*8,dimension(100) :: rho_array + real*8,dimension(13) :: eps_array + real*8,dimension(796) :: rho_test,R_test,T_test,Ye_test + character(len=13),dimension(4) :: name_list + + + !fermi function + real*8 :: fermidirac,fermidirac_dimensionless + real*8 :: x + real*8 :: find_s,find_s2 + real*8 :: find_g + real*8 :: nue_absorption_on_n + real*8 :: anue_absorption_on_p + + + real*8,parameter :: nulib_emissivity_gf = 5.59424238d-55/(6.77140812d-06**3*2.03001708d+05) + real*8,parameter :: nulib_opacity_gf = 1.0d0/6.77140812d-06 + real*8,parameter :: nulib_energy_gf = 1.60217733d-6*5.59424238d-55 + real*8,parameter :: nulib_kernel_gf = 6.77140812d-06**3/2.03001708d+05 + + !allocate the arrays for the point values + allocate(local_emissivity(mypoint_number_output_species,mypoint_number_groups)) + allocate(local_absopacity(mypoint_number_output_species,mypoint_number_groups)) + allocate(local_scatopacity(mypoint_number_output_species,mypoint_number_groups)) + allocate(local_delta(mypoint_number_output_species,mypoint_number_groups)) + allocate(blackbody_spectra(mypoint_number_output_species,mypoint_number_groups)) + + call input_parser(parameters) + + !this sets up many cooefficients and creates the energy grid (one + !zone + log spacing) see nulib.F90:initialize_nulib + call initialize_nulib(mypoint_neutrino_scheme,mypoint_number_species,mypoint_number_groups) +#if NUCLEI_HEMPEL + call set_up_Hempel ! set's up EOS for nuclear abundances +#endif + !read in EOS table & set reference mass + call read_eos_table(eos_filename) + m_ref = m_amu !for SFHo_EOS (Hempel) + ! m_ref = m_n !for LS220 +!~ name_list(1)='brem_4e11.txt' +!~ name_list(2)=('brem_4e12.txt') +!~ name_list(3)=('brem_4e13.txt') +!~ name_list(4)=('brem_4e14.txt') + + +!~ rho_figure(1) = 4.0d11 +!~ rho_figure(2) = 4.0d12 +!~ rho_figure(3) = 4.0d13 +!~ rho_figure(4) = 4.0d14 + + + !example point + xrho = 8.36008541071749d12!(1*10**(-1.00d0)) * m_n*mev_to_gram* 1.0d39 ! 1.665d13!! 4.435d13 !g/cm^3 + xtemp = 9.494951990570902!3.0d0*(xrho/1.0d11)**(1.0d0/3.0d0)!14.511804010728726d0! !MeV + xye = 0.08210467543349459 + !set up energies bins + do_integrated_BB_and_emissivity = .true. + mindx = 2.0d0 + bin_bottom(1) = 0.0d0 !MeV + bin_bottom(2) = 2.0d0 !MeV + bin_bottom(3) = bin_bottom(2)+mindx + bin_bottom(number_groups) = 400.0d0 ! MeV + + call nulib_series2(number_groups-1,bin_bottom(2),bin_bottom(number_groups),mindx,dxfac) + do i=4,number_groups + bin_bottom(i) = bin_bottom(i-1)+(bin_bottom(i-1)-bin_bottom(i-2))*dxfac + enddo + + !calculate bin widths & energies from the bottom of the bin & energy at top on bin + do i=1,number_groups-1 + energies(i) = (bin_bottom(i)+bin_bottom(i+1))/2.0d0 + bin_widths(i) = bin_bottom(i+1)-bin_bottom(i) + bin_top(i) = bin_bottom(i+1) + enddo + + energies(number_groups) = bin_bottom(number_groups)+bin_widths(number_groups-1)*dxfac/2.0d0 + bin_widths(number_groups) = 2.0*(energies(number_groups)-bin_bottom(number_groups)) + bin_top(number_groups) = bin_bottom(number_groups)+bin_widths(number_groups) + + + + + + allocate(eos_variables(total_eos_variables)) + eos_variables(:) = 0.0d0 + eos_variables(rhoindex) = xrho + eos_variables(tempindex) = xtemp + eos_variables(yeindex) =xye + + + !! EOS stuff + call set_eos_variables(eos_variables) + write(*,*) "Rho: ",eos_variables(rhoindex)," g/ccm" + write(*,*) "T: ",eos_variables(tempindex)," MeV" + write(*,*) "Ye: ",eos_variables(yeindex) + write(*,*) "X_n: ",eos_variables(xnindex) + write(*,*) "X_p: ",eos_variables(xpindex) + write(*,*) "X_alpha: ",eos_variables(xaindex) + +!~ write(*,*) energy_shift + +!~ write(*,*) eos_variables(energyindex) +!~ write(*,*) eos_variables(xpindex) +!~ eos_variables(:) = 0.0d0 +!~ eos_variables(tempindex) = 0.012d0 + +!~ open(unit=100, file="EOS_table_T=0.dat") +!~ do i = 1,size(rho_array) +!~ eos_variables(rhoindex) = 10**rho_array(i) +!~ do j = 1,13 +!~ eos_variables(yeindex) = ye_array(j) +!~ call set_eos_variables(eos_variables) +!~ eps_array(j) = eos_variables(energyindex) +!~ enddo +!~ write(100,*) eps_array +!~ enddo + +!~ close(100) +!~ stop + + +!~ open(unit=100,file='test_temp_rho.dat') + +!~ do i = 1,size(R_test) +!~ read(100,*,IOSTAT=j) R_test(i),rho_test(i),T_test(i),Ye_test(i) + +!~ enddo +!~ close(100) + if (add_nue_kernel_bremsstrahlung.or.add_anue_kernel_bremsstrahlung.or. & + add_numu_kernel_bremsstrahlung.or.add_anumu_kernel_bremsstrahlung.or. & + add_nutau_kernel_bremsstrahlung.or.add_anutau_kernel_bremsstrahlung .or. & + add_numu_kernel_gangguo .or. add_anumu_kernel_gangguo .or. add_nutau_kernel_gangguo & + .or. add_anutau_kernel_gangguo) then + + write(*,*) + write(*,*) + write(*,*) "brem call" + + + allocate(local_Phi0_bremsstrahlung(mypoint_number_output_species,mypoint_number_groups,2)) + allocate(local_Phi0_bremsstrahlung2(mypoint_number_output_species,mypoint_number_groups,2)) + allocate(local_Phi0_epannihil(mypoint_number_output_species,mypoint_number_groups,2)) + allocate(local_Phi1_epannihil(mypoint_number_output_species,mypoint_number_groups,2)) + allocate(Phi0_brem(mypoint_number_groups)) + allocate(Phi0_brem2(mypoint_number_groups)) + allocate(Phi0_brem_ann(mypoint_number_groups)) + allocate(Phi0_brem2_ann(mypoint_number_groups)) + + + write(*,*) "brem emi", single_neutrino_emissivity_from_NNBrem_given_energypoint( & + 3,energies(mypoint_number_groups),eos_variables) + + write(*,*) "echem",eos_variables(mueindex),energies(11) + write(*,*) "echem fermi",fermidirac_dimensionless( & + energies(11),-eos_variables(mueindex))& + ,exp(-energies(11)-eos_variables(mueindex)) + + + + + +!~ do zone= 1,4!size(R_test) +!~ write(*,*) name_list(zone),rho_figure(zone) +!~ write(*,*) 'test' +!~ open(unit = 404, file = name_list(zone) ) + open(unit = 404, file = "brem.txt" ) + +!~ xrho = rho_figure(zone)!(1*10**(-1.00d0)) * m_n*mev_to_gram* 1.0d39 ! 1.665d13!! 4.435d13 !g/cm^3 +!~ xtemp = 3.0d0*(xrho/1.0d11)**(1.0d0/3.0d0)!14.511804010728726d0! !MeV +!~ xye = 0.2 +!~ write(*,*) xrho,xtemp,xye +!~ stop + +!~ eos_variables(rhoindex) = xrho +!~ eos_variables(tempindex) = xtemp +!~ eos_variables(yeindex) = xye + + + call set_eos_variables(eos_variables) + + ratio_2 = 0.0d0 + S_sigma_static = 0.0d0 + do i=1,mypoint_number_groups !neutrino energ + + +!~ dE_dEdt = 0.0d0 +!~ dE_dEdt2 = 0.0d0 +!~ dE_dEdt3 = 0.0d0 +!~ dn_dEdt = 0.0d0 +!~ dn_dEdt2 = 0.0d0 +!~ dn_dEdt3 = 0.0d0 + Phi0_brem = 0.0d0 + Phi0_brem2 = 0.0d0 + Phi0_brem_ann = 0.0d0 + Phi0_brem2_ann = 0.0d0 + call single_epannihil_kernel_point_return_all(i, & + eos_variables(mueindex)/eos_variables(tempindex),& +!~ 0.0d0/eos_variables(tempindex),& + eos_variables(tempindex),local_Phi0_epannihil & + ,local_Phi1_epannihil,mypoint_neutrino_scheme) + + + + do inde =1,3 + local_Phi0_bremsstrahlung = 0.0d0 + local_Phi0_bremsstrahlung2 = 0.0d0 +!~ n_N =0.0d0 + + + ! neutron-neutron n_n + if (inde .EQ. 1) then + + n_N = eos_variables(rhoindex)& + *eos_variables(xnindex)& + /(m_n*mev_to_gram) + + + ! proton-proton n_n + else if (inde .EQ. 2 ) then + + n_N = eos_variables(rhoindex)& + *eos_variables(xpindex)& + /(m_p*mev_to_gram) + + ! neutron-proton n_n + else if (inde .EQ. 3 ) then + n_N = eos_variables(rhoindex) & + *sqrt(eos_variables(xnindex)*eos_variables(xpindex))& + /(sqrt(m_p*m_n)*mev_to_gram) + else + write(*,*) "Wrong number of species for brem test" + endif + + + + ! Hannestad version + call single_bremsstrahlung_kernel_point_return_all_Hannestad(i, & + n_N, & + eos_variables(tempindex) & + ,local_Phi0_bremsstrahlung,mypoint_neutrino_scheme) + +!~ ! Kuroda version + call single_bremsstrahlung_kernel_point_return_all_gang(i, & + eos_variables(rhoindex)/(m_n*mev_to_gram)& +!~ n_N & + ,eos_variables(yeindex), & + eos_variables(tempindex)& + ,local_Phi0_bremsstrahlung2,mypoint_neutrino_scheme) + +!~ write(*,*) eos_variables(rhoindex)/(m_amu*mev_to_gram),n_N + ! Hannestad production kernels + if (inde .LT. 3) then + Phi0_brem = Phi0_brem + local_Phi0_bremsstrahlung(3,:,1) + + + else if ( inde .EQ. 3) then ! factor for neutron-proton process + Phi0_brem = Phi0_brem + & + 28.0d0/3.0d0*local_Phi0_bremsstrahlung(3,:,1) + endif + + ! Hannestad annihilation kernels + if (inde .LT. 3) then + Phi0_brem_ann = Phi0_brem_ann + local_Phi0_bremsstrahlung(3,:,2) + + + else if ( inde .EQ. 3) then ! factor for neutron-proton process + Phi0_brem_ann = Phi0_brem_ann + & + 28.0d0/3.0d0*local_Phi0_bremsstrahlung(3,:,2) + endif + + + + ! Kuroda production kernels +!~ if (inde .LT. 3) then +!~ Phi0_brem2 = Phi0_brem2 + local_Phi0_bremsstrahlung2(3,:,1) + +!~ else if ( inde .EQ. 3) then ! factor for neutron-proton process +!~ Phi0_brem2 = Phi0_brem2 + & +!~ 28.0d0/3.0d0*local_Phi0_bremsstrahlung2(3,:,1) +!~ endif + if ( inde .EQ. 1) then + Phi0_brem2 = Phi0_brem2 + local_Phi0_bremsstrahlung2(3,:,1) + endif + + + + ! Kuroda annihilation kernels +!~ if (inde .LT. 3) then +!~ Phi0_brem2_ann = Phi0_brem2_ann + local_Phi0_bremsstrahlung2(3,:,2) + + if ( inde .EQ. 1) then ! factor for neutron-proton process + Phi0_brem2_ann = Phi0_brem2_ann + & + local_Phi0_bremsstrahlung2(3,:,2) + endif + enddo + + + + Q = 0.0d0 + Q_ann = 0.0d0 + Q2 = 0.0d0 + Q2_ann = 0.0d0 + M = 0.0d0 + M2 = 0.0d0 + N = 0.0d0 + ratio = 0.0d0 + ratio_2 = 0.0d0 + + + ! Loop over antineutrino energies + do j=1,mypoint_number_groups +!~ write(*,*) j + ! Kotake pair production rates + +!~ ! Hannestad + Q = Q + energies(j)**2 /(2.0d0* pi * hbarc_mevcm)**3& + *2.0d0*pi *bin_widths(j) & + * Phi0_brem(j) +!~ * local_Phi0_bremsstrahlung(3,j,1) + +!~ ! Hannestad +!~ Q = Q + (dble(j)/10.0d0)**2 /(2.0d0* pi * hbarc_mevcm)**3& +!~ *2.0d0*pi *((dble(j+1)-dble(j))/10.0d0)&!*bin_widths(j) & +!~ * Bremsstrahlung_Phi0_Hannestad((dble(i))/xtemp,(dble(j)/10.0d0)/xtemp,& +!~ xtemp,xrho/(m_n*mev_to_gram),3,0) + + Q_ann = Q_ann + energies(j)**2 /(2.0d0* pi * hbarc_mevcm)**3& + * bin_widths(j) *2.0d0*pi& + * Phi0_brem_ann(j) + + + ! Kuroda/gang + Q2 = Q2 + energies(j)**2 /(2.0d0* pi * hbarc_mevcm)**3 & + * bin_widths(j) *2.0d0*pi & +!~ * Phi0_brem2(j) + * bremsstrahlung_Phi0_gang(energies(i)/eos_variables(tempindex) & + ,energies(j)/eos_variables(tempindex) & + ,eos_variables(tempindex), eos_variables(rhoindex)/(m_n*mev_to_gram),xye,3,0) + +!~ ! Kuroda/gang +!~ Q2 = Q2 + (dble(j)/10.0d0)**2 /(2.0d0* pi * hbarc_mevcm)**3 & +!~ *2.0d0*pi *((dble(j+1)-dble(j))/10.0d0)&!* bin_widths(j) & +!~ * Phi0_brem2(j) +!~ * bremsstrahlung_Phi0_gang((dble(i))/xtemp,(dble(j)/10.0d0)/xtemp & +!~ ,xtemp, xrho/(m_n*mev_to_gram),xye,3,0) +!~ stop + Q2_ann = Q2_ann + energies(j)**2 /(2.0d0* pi * hbarc_mevcm)**3 & + * bin_widths(j)*2.0d0*pi * & + bremsstrahlung_Phi0_gang(energies(i)/eos_variables(tempindex), & + energies(j)/eos_variables(tempindex) & + ,eos_variables(tempindex), eos_variables(rhoindex)/(m_n*mev_to_gram),xye,3,1) + + ! electron-positron + M = M + energies(j)**2 /(2.0d0* pi * hbarc_mevcm)**3& + *2.0d0*pi* local_Phi0_epannihil(3,j,1) * bin_widths(j) + M2 = M2 + energies(j)**2 /(2.0d0* pi * hbarc_mevcm)**3& + *2.0d0*pi* local_Phi0_epannihil(3,j,2) * bin_widths(j) + + + enddo + + + !!!number production spectra + + ! Hannestad + dn_dEdt =1.0d0/(2.0d0 * pi *hbarc_mevcm)**3 *energies(i)**2 & + * Q + ! Kuroda + dn_dEdt2 =1.0d0/(2.0d0 * pi * hbarc_mevcm)**3 * energies(i)**2 & + * Q2 + ! positron-electron pair + dn_dEdt3 =1.0d0/(2.0d0 * pi * hbarc_mevcm)**3 * energies(i)**2 & + *M + ! Burrows + dn_dEdt_burrows =single_neutrino_emissivity_from_NNBrem_given_energypoint(3,& + energies(i),eos_variables)/mev_to_erg*(4.0d0*pi)/energies(i)& + /(4.0d0 *pi) + + ! Bruenn epan + dn_dEdt_bruenn =single_neutrino_emissivity_from_epannihil_given_energypoint(3,& + energies(i),eos_variables)/mev_to_erg*(4.0d0*pi)/energies(i)& + /(4.0d0 *pi) + + + !!!energy production spectra + + ! Hannestad + dE_dEdt =1.0d0/(2.0d0 * pi * hbarc_mevcm)**3 * energies(i)**3 & + * (4.0d0 * pi) * (Q ) + + ! Kuroda + dE_dEdt2 =1.0d0/(2.0d0 * pi * hbarc_mevcm)**3 * energies(i)**3 & + * 4.0d0 * pi * Q2 +!~ write(*,*) "1",dE_dEdt +!~ write(*,*) "2",dE_dEdt2 + ! positron-electron pair + dE_dEdt3 =1.0d0/(2.0d0 * pi *hbarc_mevcm)**3 * energies(i)**3 & + * (4.0d0 * pi) * M + + dE_dEdt1_ann =1.0d0/(2.0d0 * pi *hbarc_mevcm)**3 * energies(i)**3 & + * (4.0d0 * pi) * Q_ann + dE_dEdt2_ann =1.0d0/(2.0d0 * pi *hbarc_mevcm)**3 * energies(i)**3 & + * (4.0d0 * pi) * Q2_ann + dE_dEdt3_ann =1.0d0/(2.0d0 * pi *hbarc_mevcm)**3 * energies(i)**3 & + * (4.0d0 * pi) * M2 + + ! Energy production spectra Burrows + dE_dEdt_burrows = single_neutrino_emissivity_from_NNBrem_given_energypoint(1,& + energies(i),eos_variables)/mev_to_erg*(4.0d0*pi) + + !!! inverse mean free path for comparaison with Bruenn 2018 + + + ratio = (dE_dEdt-dE_dEdt_burrows) /dE_dEdt *100.0d0 + + + write(404,*) energies(i) & + ,dn_dEdt,dn_dEdt2,dn_dEdt3 & + ,dn_dEdt_burrows & + ,dn_dEdt_bruenn & + ,dE_dEdt &!/n_N & ! + ,dE_dEdt2 &!/n_N & ! + ,dE_dEdt3 &!/n_N & ! + ,dE_dEdt_burrows &!/n_N & ! + ,dE_dEdt-dE_dEdt2,dE_dEdt-dE_dEdt_burrows & + ,ratio,dE_dEdt2-dE_dEdt_burrows + enddo +!~ enddo +eta_star = (hbarc_mevcm/clight * (3.0d0 * pi **2 * xrho/(m_n*mev_to_gram))& + **(1.0d0/3.0d0))**2 & + /(2.0d0 *m_amu/clight**2 * xtemp) !adimensional + +!~ write(*,*) eta_star +close(404) + + + +stop + + +write(*,*) "S_sigma_static" + + +open(unit=100, file="S_sigma_static.dat") + do in_N=1,25 + Itable_n_N(in_N) = & + (-4.0d0+dble(in_N-1)/dble(24)*(4.0d0)) + enddo +!~ write(*,*) Itable_n_N +do inde = 100,400 + xrho = 10**(-dble(inde)/100.0d0) * m_n*mev_to_gram* 1.0d39 !!g/cm^3 + xtemp = 3.0d0*(xrho/1.0d11)**(1.0d0/3.0d0)!14.511804010728726d0! !MeV + N=0.0d0 + M=0.0d0 + Q2 = 0.0d0 + do i =2,40000 + Q = dble(i)/100.0d0!10**(om_array(i)) + M = M+ & + Bremsstrahlung_Phi0_Hannestad(Q/xtemp/2.0d0,Q/xtemp/2.0d0& + ,xtemp,xrho/(m_n*mev_to_gram),3,1)& + /(6.0d0 * (gA/2.0d0)**2 * (Gfermi*hbarc_mevcm)**2 & + * xrho/(m_n*mev_to_gram)*(hbarc_mevcm)**3 *clight) & + /(2.0d0*pi)*(dble(i)-dble(i-1))/100.0d0& + * ( 1.0d0 + exp(-((Q)/xtemp)))! /10**(0.1) + N = N+ & + bremsstrahlung_Phi0_gang(Q/xtemp/2.0d0,Q/xtemp/2.0d0 & + ,xtemp, xrho/(m_n*mev_to_gram),xye,3,1) & + /(6.0d0*(Gfermi*hbarc_mevcm) ** 2 * (gA/2.0d0)**2 & + * xrho/(m_n*mev_to_gram) & + *(hbarc_mevcm)**3 *clight) & + /(2.0d0*pi)* (dble(i)-dble(i-1))/100.0d0& + * ( 1.0d0 +exp(-((Q)/xtemp)) ) !*(10**(om_array(i))-10**(om_array(i-1))) + + Q2 = Q2 + bremsstrahlung_Phi0_gang(Q/xtemp/2.0d0,Q/xtemp/2.0d0 & + ,xtemp, xrho/(m_n*mev_to_gram),xye,3,1) & + /(6.0d0*(Gfermi*hbarc_mevcm) ** 2 * (gA/2.0d0)**2 & + * xrho/(m_n*mev_to_gram) & + *(hbarc_mevcm)**3 *clight) & + *(dble(i)-dble(i-1))/100.0d0 & + *(Gfermi*hbarc_mevcm*1.0d13)**2*gA**2/(160.0d0*xtemp**6) * Q**5 *exp(-((Q)/xtemp)) & + * 1.0d13 * (hbarc_mevcm*1.0d13)**3 + enddo + +!~ enddo + eta_star = (hbarc_mevcm/clight * (3.0d0 * pi **2 * n_N)**(1.0d0/3.0d0))**2 & + /(2.0d0 *m_amu/clight**2 * (real(j)/10.0d0)) !adimensional + write(100,*) xrho/(m_n*mev_to_gram) * 1d-39, N,M,& + 1.5d0/((hbarc_mevcm/clight * (3.0d0 * pi **2 * & + xrho/(m_n*mev_to_gram))**(1.0d0/3.0d0))**2 & + /(2.0d0 *m_n/clight**2 * xtemp)),Q2 !adimensional + + +enddo + + +close(100) + +xrho = 1e-2 * m_n*mev_to_gram* 1.0d39 ! 1.665d13!! 4.435d13 !g/cm^3 + +xtemp = 3.0d0*(xrho/1.0d11)**(1.0d0/3.0d0)!14.511804010728726d0! !MeV +xye = 0.2d0!0.28092030837767307d0 !dimensionless + + + open(unit=60,file="test_phi.txt") + +do i=1,400 + Q = 0.0d0 + Q2 = 0.0d0 + do j =1,40000 + Q = Q + (dble(j)/100.0d0)**2 /(2.0d0* pi * hbarc_mevcm)**3& + *2.0d0*pi *((dble(j+1)-dble(j))/100.0d0)&!*bin_widths(j) & + * Bremsstrahlung_Phi0_Hannestad((dble(i)/1.0d0)/xtemp,(dble(j)/100.0d0)/xtemp,& + xtemp,xrho/(m_n*mev_to_gram),3,0) + + Q2 = Q2 + (dble(j)/100.0d0)**2 /(2.0d0* pi * hbarc_mevcm)**3 & + *2.0d0*pi *((dble(j+1)-dble(j))/100.0d0)&!* bin_widths(j) & +!~ * Phi0_brem2(j) + * bremsstrahlung_Phi0_gang((dble(i)/1.0d0)/xtemp,(dble(j)/100.0d0)/xtemp & + ,xtemp, xrho/(m_n*mev_to_gram),xye,3,0) + enddo + ! Hannestad + dn_dEdt =1.0d0/(2.0d0 * pi *hbarc_mevcm)**3 * (dble(i)/1.0d0)**2&!energies(i)**2 & + * Q + ! Kuroda + dn_dEdt2 =1.0d0/(2.0d0 * pi * hbarc_mevcm)**3 * (dble(i)/1.d0)**2&!energies(i)**2 & + * Q2 + write(60,*) dble(i)/1.0d0,dn_dEdt,dn_dEdt2 +enddo +close(60) +stop + + + !! reproduce Fig. 4.2 in Raffelt 1985 + + open(unit=5,file="find_s.txt") + do j = 10,200,10 + + eta_star = (hbarc_mevcm/clight * (3.0d0 * pi **2 * n_N)**(1.0d0/3.0d0))**2 & + /(2.0d0 *m_amu/clight**2 * (real(j)/10.0d0)) !adimensional + + do inde = 1,mypoint_number_groups + Q = 0.0d0 + Q2 = 0.0d0 + do i = 1,mypoint_number_groups + Q = Q + find_s(eta_star,0.0d0,& + (energies(i)+energies(inde))/(real(j)/10.0d0))& + * bin_widths(i) + Q2 = Q2 + find_s2(eta_star,0.0d0,& + (energies(i)+energies(inde))/(real(j)/10.0d0))& + * bin_widths(i) + enddo + write(5,*) (energies(inde)+energies(i))/(real(j)/10.0d0),& + eta_star,find_s(eta_star,0.0d0,energies(inde)/(real(j)/10.0d0)),& + ABS(find_s(eta_star,0.0d0,energies(inde)/(real(j)/10.0d0))& + -find_s2(eta_star,0.0d0,energies(inde)/(real(j)/10.0d0)))& + /find_s(eta_star,0.0d0,energies(inde)/(real(j)/10.0d0)),& + real(j)/10.0d0,& + (energies(inde)+energies(i))/(real(j)/10.0d0)*5.0d0/8.0d0& + *exp(-(energies(inde)+energies(i))/(real(j)/10.0d0)) & + *find_s(eta_star,0.0d0,(energies(inde)+energies(i))/(real(j)/10.0d0)),& + (energies(inde)+energies(i))/(real(j)/10.0d0)*5.0d0/8.0d0& + *exp(-(energies(inde)+energies(i))/(real(j)/10.0d0)) & + *find_s2(eta_star,0.0d0,(energies(inde)+energies(i))/(real(j)/10.0d0)) + enddo + enddo + close(5) + + + + + + deallocate(Phi0_brem) + deallocate(local_Phi0_bremsstrahlung) + deallocate(local_Phi0_epannihil) + deallocate(local_Phi1_epannihil) + endif + + +end program bremsstrahlung_point_example + diff --git a/src/constants.inc b/src/constants.inc index 523069e..484d58c 100644 --- a/src/constants.inc +++ b/src/constants.inc @@ -8,6 +8,7 @@ real*8, parameter :: m_n = 939.565346d0 !MeV real*8, parameter :: m_p = 938.272013d0 !MeV real*8, parameter :: m_amu = 931.494061d0 !MeV + real*8, parameter :: m_pi = 139.57018d0 !MeV real*8, parameter :: m_alpha = 3727.379240 !MeV real*8, parameter :: sin2thetaW = 0.23d0 !dimensionless real*8, parameter :: clight = 29979245800.0d0 !cm/s @@ -19,6 +20,7 @@ real*8, parameter :: planck = 6.626176d-27 !h (i.e. hbar*2*pi), in cgs real*8, parameter :: avo = 6.0221367d23 !N_A, avogadro's number + !conversions real*8, parameter :: mev_to_gram = 1.782661758d-27 !one MeV is # grams real*8, parameter :: mev_to_erg = 1.60217733d-6 !one MeV is # ergs diff --git a/src/electron_positron_annihilation.F90 b/src/electron_positron_annihilation.F90 index d29ee5b..61f9494 100644 --- a/src/electron_positron_annihilation.F90 +++ b/src/electron_positron_annihilation.F90 @@ -426,3 +426,4 @@ function d1(E1,E2,E3) result(answer) end function d1 end function epannihil_H1_Bruenn + diff --git a/src/emissivities.F90 b/src/emissivities.F90 index 8a02b41..9de69f5 100644 --- a/src/emissivities.F90 +++ b/src/emissivities.F90 @@ -73,6 +73,7 @@ function single_neutrino_emissivity_from_epannihil_given_energypoint( & Gfermi**2*hbarc_mevcm**2*clight/pi !units (MeV*cm)^-6 * MeV^9 * erg*MeV^-1 * MeV^-4 (MeV*cm)^2 * cm*s^-1 = erg/cm^3/s eta = eos_variables(mueindex)/eos_variables(tempindex) +!~ eta = 0.0d0 nu_energy_x = nu_energy/eos_variables(tempindex) emissivity = epannihil_dQdenu_BRT06(nu_energy_x,eta,neutrino_species)*nu_energy_x**3 @@ -204,7 +205,7 @@ subroutine total_emissivities(neutrino_species,energy_point,energy_bottom,energy endif !add in the electron neutrino emission from Nucleon-Nucleon bremsstrahlung - if (add_nue_emission_NNBrems) then + if (add_nue_emission_bremsstrahlung) then if (do_integrated_BB_and_emissivity) then total_emissivity = total_emissivity + & !total emmissivity, dimensions ergs/cm^3/s/MeV/srad single_neutrino_emissivity_from_NNBrem_given_energyrange( & @@ -232,7 +233,7 @@ subroutine total_emissivities(neutrino_species,energy_point,energy_bottom,energy endif !add in the electron antineutrino emission from Nucleon-Nucleon bremsstrahlung - if (add_anue_emission_NNBrems) then + if (add_anue_emission_bremsstrahlung) then if (do_integrated_BB_and_emissivity) then total_emissivity = total_emissivity + & !total emmissivity, dimensions ergs/cm^3/s/MeV/srad single_neutrino_emissivity_from_NNBrem_given_energyrange( & @@ -260,7 +261,7 @@ subroutine total_emissivities(neutrino_species,energy_point,energy_bottom,energy endif !add in the mu neutrino emission from Nucleon-Nucleon bremsstrahlung - if (add_numu_emission_NNBrems) then + if (add_numu_emission_bremsstrahlung) then if (do_integrated_BB_and_emissivity) then total_emissivity = total_emissivity + & !total emmissivity, dimensions ergs/cm^3/s/MeV/srad single_neutrino_emissivity_from_NNBrem_given_energyrange( & @@ -288,7 +289,7 @@ subroutine total_emissivities(neutrino_species,energy_point,energy_bottom,energy endif !add in the mu antineutrino emission from Nucleon-Nucleon bremsstrahlung - if (add_anumu_emission_NNBrems) then + if (add_anumu_emission_bremsstrahlung) then if (do_integrated_BB_and_emissivity) then total_emissivity = total_emissivity + & !total emmissivity, dimensions ergs/cm^3/s/MeV/srad single_neutrino_emissivity_from_NNBrem_given_energyrange( & @@ -298,6 +299,8 @@ subroutine total_emissivities(neutrino_species,energy_point,energy_bottom,energy single_neutrino_emissivity_from_NNBrem_given_energypoint( & neutrino_species,energy_point,eos_variables) endif +!~ write(*,*) single_neutrino_emissivity_from_NNBrem_given_energyrange( & +!~ neutrino_species,energy_bottom,energy_top,eos_variables) endif endif @@ -316,7 +319,7 @@ subroutine total_emissivities(neutrino_species,energy_point,energy_bottom,energy endif !add in the tau neutrino emission from Nucleon-Nucleon bremsstrahlung - if (add_nutau_emission_NNBrems) then + if (add_nutau_emission_bremsstrahlung) then if (do_integrated_BB_and_emissivity) then total_emissivity = total_emissivity + & !total emmissivity, dimensions ergs/cm^3/s/MeV/srad single_neutrino_emissivity_from_NNBrem_given_energyrange( & @@ -332,6 +335,7 @@ subroutine total_emissivities(neutrino_species,energy_point,energy_bottom,energy if (neutrino_species.eq.6) then !add in the tau antineutrino emission from ep annihilation if (add_anutau_emission_epannihil) then + if (do_integrated_BB_and_emissivity) then total_emissivity = total_emissivity + & !total emmissivity, dimensions ergs/cm^3/s/MeV/srad single_neutrino_emissivity_from_epannihil_given_energyrange( & @@ -344,7 +348,7 @@ subroutine total_emissivities(neutrino_species,energy_point,energy_bottom,energy endif !add in the tau antineutrino emission from Nucleon-Nucleon bremsstrahlung - if (add_anutau_emission_NNBrems) then + if (add_anutau_emission_bremsstrahlung) then if (do_integrated_BB_and_emissivity) then total_emissivity = total_emissivity + & !total emmissivity, dimensions ergs/cm^3/s/MeV/srad single_neutrino_emissivity_from_NNBrem_given_energyrange( & @@ -406,6 +410,7 @@ subroutine return_emissivity_spectra_given_neutrino_scheme(emissivity_spectra,eo emissivity_spectra(ns,:) = emissivity_spectra(ns,:) + ec_emissivity(:) !erg/cm^3/s/MeV/srad end if if (add_anue_emission_weakinteraction_poscap.and.ns.eq.2) then + stop "emissivities :: anue weak rates are not yet implemented" call microphysical_electron_capture(ns,eos_variables,ec_emissivity) emissivity_spectra(ns,:) = emissivity_spectra(ns,:) + ec_emissivity(:) !erg/cm^3/s/MeV/srad end if diff --git a/src/make_table_example.F90 b/src/make_table_example.F90 index 6adcbe1..cae8ccf 100644 --- a/src/make_table_example.F90 +++ b/src/make_table_example.F90 @@ -31,7 +31,7 @@ program make_table_example integer :: number_output_species = 3 !number of energy groups - integer :: mytable_number_groups = 18 + integer :: mytable_number_groups = 12 !NuLib parameters file (weak rates and EOS) character*200 :: parameters_filename = "./parameters" @@ -53,17 +53,25 @@ program make_table_example real*8, allocatable,dimension(:,:,:,:,:) :: table_delta + !final Itable parameters integer :: final_Itable_size_temp, final_Itable_size_eta, final_Itable_size_inE + integer :: final_Itable_size_n_N,final_Itable_size_Ye real*8 :: Imin_logtemp,Imax_logtemp real*8 :: Imin_logeta,Imax_logeta + real*8 :: Imin_logn_N,Imax_logn_N + real*8 :: Imin_logYe,Imax_logYe real*8, allocatable,dimension(:) :: Itable_temp real*8, allocatable,dimension(:) :: Itable_eta real*8, allocatable,dimension(:) :: Itable_inE + real*8, allocatable,dimension(:) :: Itable_n_N + real*8, allocatable,dimension(:) :: Itable_Ye real*8, allocatable,dimension(:,:,:,:,:) :: Itable_Phi0 real*8, allocatable,dimension(:,:,:,:,:) :: Itable_Phi1 real*8, allocatable,dimension(:,:,:,:,:,:) :: epannihiltable_Phi0 !for ep-annihilation kernels, need both production and destruction kernels real*8, allocatable,dimension(:,:,:,:,:,:) :: epannihiltable_Phi1 !for ep-annihilation kernels, need both production and destruction kernels + real*8, allocatable,dimension(:,:,:,:,:,:) :: bremsstrahlungtable_Phi0 !for ep-annihilation kernels, need both production and destruction kernels + real*8, allocatable,dimension(:,:,:,:,:,:,:) :: bremsstrahlungtable_gang_Phi0 !for ep-annihilation kernels, need both production and destruction kernels !versioning @@ -74,7 +82,7 @@ program make_table_example !local variables to help in making tables integer :: irho,itemp,iye,ns,ng - integer :: ieta,iinE + integer :: ieta,iinE,in_N real*8, allocatable,dimension(:,:) :: local_emissivity real*8, allocatable,dimension(:,:) :: local_absopacity real*8, allocatable,dimension(:,:) :: local_scatopacity @@ -83,13 +91,16 @@ program make_table_example real*8, allocatable,dimension(:,:) :: local_Phi1 real*8, allocatable,dimension(:,:,:) :: local_Phi0_epannihil real*8, allocatable,dimension(:,:,:) :: local_Phi1_epannihil + real*8, allocatable,dimension(:,:,:) :: local_Phi0_bremsstrahlung + real*8, allocatable,dimension(:,:,:) :: local_Phi0_bremsstrahlung_gang real*8, allocatable,dimension(:) :: eos_variables real*8 :: matter_prs,matter_ent,matter_cs2,matter_dedt,matter_dpderho,matter_dpdrhoe integer :: keytemp,keyerr real*8 :: precision = 1.0d-10 integer :: i real*8 dxfac,mindx - logical :: doing_inelastic, doing_epannihil + logical :: doing_inelastic, doing_epannihil,doing_bremsstrahlung,doing_bremsstrahlung_gang + logical :: skip_emission #ifdef __MPI__ !MPI variables @@ -107,6 +118,7 @@ program make_table_example real*8, allocatable,dimension(:,:,:,:,:) :: Itable_Phi0_node real*8, allocatable,dimension(:,:,:,:,:) :: Itable_Phi1_node integer :: mpi_final_Itable_size_temp + !MPI initialization @@ -129,28 +141,36 @@ program make_table_example base="NuLib" vnum="1.0" + skip_emission = .false. + adhoc_nux_factor = 0.0d0 !increase for adhoc nux heating (also set !add_nux_absorption_on_n_and_p to true) !set up table - final_table_size_ye = 51 + final_table_size_ye = 56 final_table_size_rho = 82 - final_table_size_temp = 65 + final_table_size_temp = 100 - final_Itable_size_temp = 65 - final_Itable_size_eta = 61 + final_Itable_size_temp = 100 + final_Itable_size_eta = 120 + final_Itable_size_n_N = 80 + final_Itable_size_Ye =50 final_Itable_size_inE = mytable_number_groups - min_ye = 0.035d0 - max_ye = 0.55d0 - min_logrho = 6.0d0 - max_logrho = 15.5d0 - min_logtemp = log10(0.05d0) - max_logtemp = log10(150.0d0) - Imin_logtemp = log10(0.05d0) - Imax_logtemp = log10(150.0d0) + min_ye = 0.01d0 + max_ye = 0.6d0 + min_logrho = -2.22024 + max_logrho = 15.5001 + min_logtemp = -2.99d0!log10(0.05d0) + max_logtemp = 2.19d0!log10(150.0d0) + Imin_logtemp = -2.99d0!log10(0.005d0) + Imax_logtemp = 2.19d0! log10(150.0d0) Imin_logeta = log10(0.1d0) Imax_logeta = log10(100.0d0) - + Imin_logn_N = 20.0d0! 35.0d0! + Imax_logn_N = log10(1.0d39)!log10(1.0d39)!45.0d0! + Imax_logYe = log10(0.6d0) + Imin_logYe = log10(0.0100d0) + !set up energies bins do_integrated_BB_and_emissivity = .false. mindx = 2.0d0 @@ -202,6 +222,9 @@ program make_table_example bin_widths(number_groups) = 2.0*(energies(number_groups)-bin_bottom(number_groups)) bin_top(number_groups) = bin_bottom(number_groups)+bin_widths(number_groups) +if ( skip_emission) then + go to 395 +endif allocate(table_ye(final_table_size_ye)) allocate(table_rho(final_table_size_rho)) @@ -273,6 +296,7 @@ program make_table_example write(*,*) "Rho:", 100.0*dble(irho-1)/dble(final_table_size_rho),"%" #endif do itemp=1,final_table_size_temp +!~ write(*,*) "Temp:", 100.0*dble(itemp-1)/dble(final_table_size_temp),"%" do iye=1,final_table_size_ye eos_variables = 0.0d0 @@ -310,12 +334,10 @@ program make_table_example eos_variables(rhoindex),eos_variables(tempindex),eos_variables(yeindex),ns,ng stop endif - if (.not. do_transport_opacities) then - if (local_delta(ns,ng).ne.local_delta(ns,ng)) then - write(*,"(a,1P3E18.9,i6,i6)") "We have a NaN in scat delta", & - eos_variables(rhoindex),eos_variables(tempindex),eos_variables(yeindex),ns,ng - stop - endif + if (local_delta(ns,ng).ne.local_delta(ns,ng)) then + write(*,"(a,1P3E18.9,i6,i6)") "We have a NaN in scat delta", & + eos_variables(rhoindex),eos_variables(tempindex),eos_variables(yeindex),ns,ng + stop endif if (log10(local_emissivity(ns,ng)).ge.300.0d0) then @@ -336,23 +358,21 @@ program make_table_example eos_variables(tempindex),eos_variables(yeindex),ns,ng stop endif - if (.not. do_transport_opacities) then - if (local_delta(ns,ng).gt.1.0d0) then - write(*,"(a,1P4E18.9,i6,i6)") "delta > 1", & - local_delta(ns,ng),eos_variables(rhoindex), & - eos_variables(tempindex),eos_variables(yeindex),ns,ng - stop - endif - if (local_delta(ns,ng).lt.-1.0d0) then - write(*,"(a,1P4E18.9,i6,i6)") "delta < -1", & - local_delta(ns,ng),eos_variables(rhoindex), & - eos_variables(tempindex),eos_variables(yeindex),ns,ng - stop - endif + if (local_delta(ns,ng).gt.1.0d0) then + write(*,"(a,1P4E18.9,i6,i6)") "delta > 1", & + local_delta(ns,ng),eos_variables(rhoindex), & + eos_variables(tempindex),eos_variables(yeindex),ns,ng + stop + endif + if (local_delta(ns,ng).lt.-1.0d0) then + write(*,"(a,1P4E18.9,i6,i6)") "delta < -1", & + local_delta(ns,ng),eos_variables(rhoindex), & + eos_variables(tempindex),eos_variables(yeindex),ns,ng + stop endif enddo !do ng=1,mytable_number_groups enddo !do ns=1,number_output_species - + !set global table do ns=1,number_output_species do ng=1,mytable_number_groups @@ -366,6 +386,7 @@ program make_table_example table_absopacity(irho,itemp,iye,ns,ng) = local_absopacity(ns,ng) !cm^-1 table_scatopacity(irho,itemp,iye,ns,ng) = local_scatopacity(ns,ng) !cm^-1 table_delta(irho,itemp,iye,ns,ng) = local_delta(ns,ng) !dimensionless + #endif enddo !do ns=1,number_output_species enddo !do ng=1,mytable_number_groups @@ -380,7 +401,7 @@ program make_table_example deallocate(eos_variables) enddo!do irho=1,final_table_size_rho !$OMP END PARALLEL DO! end do - +395 CONTINUE #ifdef __MPI__ call mpi_barrier(mpi_comm_world, ierror) if(mpirank.eq.0)write(*,*) "Finished Opacity Table" @@ -431,8 +452,27 @@ program make_table_example else doing_epannihil = .false. endif + + if (add_nue_kernel_bremsstrahlung.or.add_anue_kernel_bremsstrahlung.or. & + add_numu_kernel_bremsstrahlung.or.add_anumu_kernel_bremsstrahlung.or. & + add_nutau_kernel_bremsstrahlung.or.add_anutau_kernel_bremsstrahlung) then + + doing_bremsstrahlung = .true. + else + doing_bremsstrahlung = .false. + endif + + if ( add_numu_kernel_gangguo.or.add_anumu_kernel_gangguo.or. & + add_nutau_kernel_gangguo.or.add_anutau_kernel_gangguo) then - if (doing_inelastic.or.doing_epannihil) then + doing_bremsstrahlung_gang = .true. + else + doing_bremsstrahlung_gang = .false. + endif + + + if (doing_inelastic.or.doing_epannihil.or.doing_bremsstrahlung & + .or.doing_bremsstrahlung_gang) then write(*,*) "Making Inelastic Table Opacity Table" @@ -475,6 +515,8 @@ program make_table_example allocate(Itable_temp(final_Itable_size_temp)) allocate(Itable_eta(final_Itable_size_eta)) allocate(Itable_inE(final_Itable_size_inE)) + allocate(Itable_n_N(final_Itable_size_n_N)) + allocate(Itable_Ye(final_Itable_size_Ye)) allocate(Itable_Phi0(final_Itable_size_temp,final_Itable_size_eta, & final_Itable_size_inE,number_output_species,mytable_number_groups)) @@ -484,6 +526,11 @@ program make_table_example final_Itable_size_inE,number_output_species,mytable_number_groups,2)) allocate(epannihiltable_Phi1(final_Itable_size_temp,final_Itable_size_eta, & final_Itable_size_inE,number_output_species,mytable_number_groups,2)) + allocate(bremsstrahlungtable_Phi0(final_Itable_size_temp,final_Itable_size_n_N, & + final_Itable_size_inE,number_output_species,mytable_number_groups,2)) + allocate(bremsstrahlungtable_gang_Phi0(final_Itable_size_temp,final_Itable_size_n_N, & + final_Itable_size_Ye,final_Itable_size_inE,number_output_species, & + mytable_number_groups,2)) #ifdef __MPI__ !mpi node tables for inelastic @@ -502,14 +549,29 @@ program make_table_example Itable_eta(ieta) = & 10.0d0**(Imin_logeta+dble(ieta-1)/dble(final_Itable_size_eta-1)*(Imax_logeta-Imin_logeta)) enddo - + + do in_N=1,final_Itable_size_n_N + Itable_n_N(in_N) = & + 10.0d0**(Imin_logn_N+dble(in_N-1)/dble(final_Itable_size_n_N-1)*(Imax_logn_N-Imin_logn_N)) + enddo + + do iye=1,final_Itable_size_Ye + Itable_Ye(iye) = & + 10.0d0**(Imin_logYe+dble(iye-1)/dble(final_Itable_size_Ye-1)*(Imax_logYe-Imin_logYe)) + enddo + +!~ write(*,*) "T", Itable_temp, size(Itable_temp) +!~ write(*,*) +!~ write(*,*) "n", Itable_n_N,size(Itable_n_N) + #ifdef __MPI__ !mpi_scatterv sends portions of Itable_temp to different nodes call mpi_scatterv(Itable_temp,sendcounts,displs,mpi_double,Itable_temp_subset,& recvcount,mpi_double,0,mpi_comm_world,ierror) mpi_final_Itable_size_temp = recvcount #endif - !$OMP PARALLEL DO PRIVATE(local_Phi0,local_Phi1,local_Phi0_epannihil,local_Phi1_epannihil,ieta,iinE,ns,ng) + !$OMP PARALLEL DO PRIVATE(local_Phi0,local_Phi1,local_Phi0_epannihil,local_Phi1_epannihil & + !$OMP ,local_Phi0_bremsstrahlung,local_Phi0_bremsstrahlung_gang,in_N,ieta,iinE,ns,ng) !loop over temp,eta,inE of table, do each point #ifdef __MPI__ do itemp=1,mpi_final_Itable_size_temp @@ -521,15 +583,18 @@ program make_table_example allocate(local_Phi1(number_output_species,mytable_number_groups)) allocate(local_Phi0_epannihil(number_output_species,mytable_number_groups,2)) allocate(local_Phi1_epannihil(number_output_species,mytable_number_groups,2)) + allocate(local_Phi0_bremsstrahlung(number_output_species,mytable_number_groups,2)) + allocate(local_Phi0_bremsstrahlung_gang(number_output_species,mytable_number_groups,2)) #ifdef __MPI__ write(*,*) "Temp:", 100.0*dble(displs(mpirank)+itemp-1)/dble(final_Itable_size_temp),"%" #else - write(*,*) "Temp:", 100.0*dble(itemp-1)/dble(final_Itable_size_temp),"%" + write(*,*) "Temp:", 100.0*dble(itemp-1)/dble(final_Itable_size_temp),"%",Itable_temp(itemp) #endif - do ieta=1,final_Itable_size_eta - do iinE=final_Itable_size_inE,1,-1 + do iinE=final_Itable_size_inE,1,-1 +!~ write(*,*) "Eta:", 100.0*dble(ieta-1)/dble(final_Itable_size_eta),"%" + do ieta=1,final_Itable_size_eta #ifdef __MPI__ call single_Ipoint_return_all(iinE,Itable_eta(ieta), & @@ -691,18 +756,120 @@ program make_table_example epannihiltable_Phi1(itemp,ieta,iinE,ns,ng,1) = local_Phi1_epannihil(ns,ng,1) !cm^3/s epannihiltable_Phi1(itemp,ieta,iinE,ns,ng,2) = local_Phi1_epannihil(ns,ng,2) !cm^3/s enddo !do ns=1,number_output_species - enddo !do ng=1,mytable_number_groups + enddo !do ng=1,mytable_number_groups + enddo!do ieta=1,final_Itable_size_eta + + if (doing_bremsstrahlung .or. doing_bremsstrahlung_gang) then + do in_N=1,final_Itable_size_n_N + call single_bremsstrahlung_kernel_point_return_all_Hannestad(iinE,Itable_n_N(in_N), & + Itable_temp(itemp),local_Phi0_bremsstrahlung,mytable_neutrino_scheme) - enddo!do iinE=1,final_Itable_size_inE - enddo!do ieta=1,final_Itable_size_eta + !calculate and check that the number is not NaN or Inf + !(.gt.1.0d300) + do ns=1,number_output_species + do ng=1,mytable_number_groups + + if (local_Phi0_bremsstrahlung(ns,ng,1).ne.local_Phi0_bremsstrahlung(ns,ng,1)) then + write(*,"(a,1P2E18.9,i6,i6,i6)") "We have a NaN in Phi0_bremsstrahlung production", & + Itable_temp(itemp),Itable_n_N(in_N),iinE,ns,ng + stop + endif + if (local_Phi0_bremsstrahlung(ns,ng,2).ne.local_Phi0_bremsstrahlung(ns,ng,2)) then + write(*,"(a,1P2E18.9,i6,i6,i6)") "We have a NaN in Phi0_bremsstrahlung annihilation", & + Itable_temp(itemp),Itable_n_N(in_N),iinE,ns,ng + stop + endif + + if (log10(local_Phi0_bremsstrahlung(ns,ng,1)).ge.300.0d0) then + write(*,"(a,1P3E18.9,i6,i6,i6)") "We have a Inf in Phi0_bremsstrahlung production", & + local_Phi0_bremsstrahlung(ns,ng,1),Itable_temp(itemp),Itable_n_N(in_N),iinE,ns,ng + stop + endif + if (log10(local_Phi0_bremsstrahlung(ns,ng,2)).ge.300.0d0) then + write(*,"(a,1P3E18.9,i6,i6,i6)") "We have a Inf in Phi0_bremsstrahlung annihilation", & + local_Phi0_bremsstrahlung(ns,ng,2),Itable_temp(itemp),Itable_n_N(in_N),iinE,ns,ng + stop + endif + + enddo !do ng=1,mytable_number_groups + enddo !do ns=1,number_output_species + + !set global table + do ns=1,number_output_species + do ng=1,mytable_number_groups + bremsstrahlungtable_Phi0(itemp,in_N,iinE,ns,ng,1) = local_Phi0_bremsstrahlung(ns,ng,1) !cm^3/s + bremsstrahlungtable_Phi0(itemp,in_N,iinE,ns,ng,2) = local_Phi0_bremsstrahlung(ns,ng,2) !cm^3/s + enddo !do ns=1,number_output_species + enddo !do ng=1,mytable_number_groups + + if (doing_bremsstrahlung_gang) then + do iye = 1,final_Itable_size_Ye + call single_bremsstrahlung_kernel_point_return_all_gang(iinE,Itable_n_N(in_N),Itable_Ye(iye),& + Itable_temp(itemp),local_Phi0_bremsstrahlung_gang,mytable_neutrino_scheme) +!~ if (maxval(local_Phi0_bremsstrahlung_gang) .GT. 0.0d0) then +!~ write(*,*) Itable_n_N(in_N),Itable_Ye(iye),Itable_temp(itemp),local_Phi0_bremsstrahlung_gang +!~ endif + !calculate and check that the number is not NaN or Inf + !(.gt.1.0d300) + do ns=1,number_output_species + do ng=1,mytable_number_groups + + if (local_Phi0_bremsstrahlung_gang(ns,ng,1).ne.local_Phi0_bremsstrahlung_gang(ns,ng,1)) then + write(*,"(a,1P2E18.9,i6,i6,i6)") "We have a NaN in Phi0_bremsstrahlung gang production", & + Itable_temp(itemp),Itable_n_N(in_N),iinE,ns,ng + stop + endif + if (local_Phi0_bremsstrahlung_gang(ns,ng,2).ne.local_Phi0_bremsstrahlung_gang(ns,ng,2)) then + write(*,"(a,1P2E18.9,i6,i6,i6)") "We have a NaN in Phi0_bremsstrahlung gang annihilation", & + Itable_temp(itemp),Itable_n_N(in_N),iinE,ns,ng + stop + endif + + if (log10(local_Phi0_bremsstrahlung_gang(ns,ng,1)).ge.300.0d0) then + write(*,"(a,1P3E18.9,i6,i6,i6)") "We have a Inf in Phi0_bremsstrahlung gang production", & + local_Phi0_bremsstrahlung_gang(ns,ng,1),Itable_temp(itemp),Itable_n_N(in_N),iinE,ns,ng + stop + endif + if (local_Phi0_bremsstrahlung_gang(ns,ng,1).LT. 0.0d0) then + write(*,"(a,1P3E18.9,i6,i6,i6)") "We have a n in Phi0_bremsstrahlung gang production", & + local_Phi0_bremsstrahlung_gang(ns,ng,1),Itable_temp(itemp),Itable_n_N(in_N),iinE,ns,ng + stop + endif + if (log10(local_Phi0_bremsstrahlung_gang(ns,ng,2)).ge.300.0d0) then + write(*,"(a,1P3E18.9,i6,i6,i6)") "We have a Inf in Phi0_bremsstrahlung gang annihilation", & + local_Phi0_bremsstrahlung_gang(ns,ng,2),Itable_temp(itemp),Itable_n_N(in_N),iinE,ns,ng + stop + endif + + enddo !do ng=1,mytable_number_groups + enddo !do ns=1,number_output_species + + !set global table + do ns=1,number_output_species + do ng=1,mytable_number_groups + bremsstrahlungtable_gang_Phi0(itemp,in_N,iye,iinE,ns,ng,1) = local_Phi0_bremsstrahlung_gang(ns,ng,1) !cm^3/s + bremsstrahlungtable_gang_Phi0(itemp,in_N,iye,iinE,ns,ng,2) = local_Phi0_bremsstrahlung_gang(ns,ng,2) !cm^3/s + enddo !do ns=1,number_output_species + enddo !do ng=1,mytable_number_groups + + enddo ! do iye = 1,final_Itable_size_Ye + endif + enddo!do in_N=1,final_Itable_size_n_N + endif + + enddo!do iinE=1,final_Itable_size_inE deallocate(local_Phi0) deallocate(local_Phi1) deallocate(local_Phi0_epannihil) deallocate(local_Phi1_epannihil) + deallocate(local_Phi0_bremsstrahlung) + deallocate(local_Phi0_bremsstrahlung_gang) enddo!do itemp=1,final_Itable_size_temp !$OMP END PARALLEL DO! end do - +!~ write(*,*) bremsstrahlungtable_Phi0(51,41,11,3,11,1) +!~ write(*,*) epannihiltable_Phi0(51,41,11,3,11,1) +!~ write(*,*) "T",Itable_temp(51),"n_N",Itable_n_N(41),"eta",Itable_eta(41),"En",energies(11) #ifdef __MPI__ call mpi_barrier(mpi_comm_world, ierror) if(mpirank.eq.0)write(*,*) "Finished Inelastic Table" @@ -735,7 +902,7 @@ program make_table_example timestamp = dble(values(1))*10000.0d0+dble(values(2))*100.0+dble(values(3)) + & (dble(values(5))+dble(values(6))/60.0d0 + dble(values(7))/3600.0d0 )/24.0 - if (doing_inelastic.or.doing_epannihil) then + if (doing_inelastic.or.doing_epannihil.or.doing_bremsstrahlung) then finaltable_filename = trim(adjustl(outdir))//trim(adjustl(base))//"_rho"//trim(adjustl(srho))// & "_temp"//trim(adjustl(stemp))//"_ye"//trim(adjustl(sye))// & "_ng"//trim(adjustl(sng))//"_ns"//trim(adjustl(sns))// & @@ -769,7 +936,7 @@ subroutine write_h5(filename,timestamp) !H5 stuff integer :: error,rank,cerror integer(HID_T) :: file_id,dset_id,dspace_id - integer(HSIZE_T) :: dims1(1), dims2(2), dims3(3), dims4(4), dims5(5), dims6(6)!, etc.... + integer(HSIZE_T) :: dims1(1), dims2(2), dims3(3), dims4(4), dims5(5), dims6(6),dims7(7)!, etc.... real*8 :: timestamp character(8) :: date @@ -932,6 +1099,8 @@ subroutine write_h5(filename,timestamp) call h5sclose_f(dspace_id, error) cerror = cerror + error + + call h5screate_simple_f(rank, dims5, dspace_id, error) call h5dcreate_f(file_id, "scattering_opacity", H5T_NATIVE_DOUBLE, & dspace_id, dset_id, error) @@ -949,8 +1118,9 @@ subroutine write_h5(filename,timestamp) call h5sclose_f(dspace_id, error) cerror = cerror + error endif - - if (doing_inelastic.or.doing_epannihil) then + + if (doing_inelastic.or.doing_epannihil.or.doing_bremsstrahlung& + .or.doing_bremsstrahlung_gang) then rank = 1 dims1(1) = 1 @@ -971,6 +1141,26 @@ subroutine write_h5(filename,timestamp) call h5dclose_f(dset_id, error) call h5sclose_f(dspace_id, error) cerror = cerror + error + + rank = 1 + dims1(1) = 1 + call h5screate_simple_f(rank, dims1, dspace_id, error) + call h5dcreate_f(file_id, "In_N", H5T_NATIVE_INTEGER, & + dspace_id,dset_id, error) + call h5dwrite_f(dset_id, H5T_NATIVE_INTEGER, final_Itable_size_n_N, dims1, error) + call h5dclose_f(dset_id, error) + call h5sclose_f(dspace_id, error) + cerror = cerror + error + + rank = 1 + dims1(1) = 1 + call h5screate_simple_f(rank, dims1, dspace_id, error) + call h5dcreate_f(file_id, "IYe", H5T_NATIVE_INTEGER, & + dspace_id,dset_id, error) + call h5dwrite_f(dset_id, H5T_NATIVE_INTEGER, final_Itable_size_Ye, dims1, error) + call h5dclose_f(dset_id, error) + call h5sclose_f(dspace_id, error) + cerror = cerror + error rank = 1 dims1(1) = final_Itable_size_temp @@ -990,6 +1180,26 @@ subroutine write_h5(filename,timestamp) call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE,Itable_eta, dims1, error) call h5dclose_f(dset_id, error) call h5sclose_f(dspace_id, error) + cerror = cerror + error + + rank = 1 + dims1(1) = final_Itable_size_n_N + call h5screate_simple_f(rank, dims1, dspace_id, error) + call h5dcreate_f(file_id, "n_N_Ipoints", H5T_NATIVE_DOUBLE, & + dspace_id, dset_id, error) + call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE,Itable_n_N, dims1, error) + call h5dclose_f(dset_id, error) + call h5sclose_f(dspace_id, error) + cerror = cerror + error + + rank = 1 + dims1(1) = final_Itable_size_Ye + call h5screate_simple_f(rank, dims1, dspace_id, error) + call h5dcreate_f(file_id, "Ye_Ipoints", H5T_NATIVE_DOUBLE, & + dspace_id, dset_id, error) + call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE,Itable_Ye, dims1, error) + call h5dclose_f(dset_id, error) + call h5sclose_f(dspace_id, error) cerror = cerror + error endif @@ -1046,7 +1256,59 @@ subroutine write_h5(filename,timestamp) cerror = cerror + error endif + + + if (doing_bremsstrahlung) then + rank = 6 + dims6(1) = final_Itable_size_temp + dims6(2) = final_Itable_size_n_N + dims6(3) = final_Itable_size_inE + dims6(4) = number_output_species + dims6(5) = number_groups + dims6(6) = 2 + + call h5screate_simple_f(rank, dims6, dspace_id, error) + call h5dcreate_f(file_id, "bremsstrahlung_phi0", H5T_NATIVE_DOUBLE, & + dspace_id, dset_id, error) + call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE,bremsstrahlungtable_Phi0, dims6, error) + call h5dclose_f(dset_id, error) + call h5sclose_f(dspace_id, error) + cerror = cerror + error + write(*,*) MAXVAL(bremsstrahlungtable_Phi0) + write(*,*) MINVAL(bremsstrahlungtable_Phi0) + endif + + if (doing_bremsstrahlung_gang) then +!~ if (.false.) then + rank = 7 + dims7(1) = final_Itable_size_temp + dims7(2) = final_Itable_size_n_N + dims7(3) = final_Itable_size_Ye + dims7(4) = final_Itable_size_inE + dims7(5) = number_output_species + dims7(6) = number_groups + dims7(7) = 2 + + call h5screate_simple_f(rank, dims7, dspace_id, error) + call h5dcreate_f(file_id, "bremsstrahlung_phi0_gang", H5T_NATIVE_DOUBLE, & + dspace_id, dset_id, error) + call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE,bremsstrahlungtable_gang_Phi0, dims7, error) + call h5dclose_f(dset_id, error) + call h5sclose_f(dspace_id, error) + cerror = cerror + error + write(*,*) "gang max : ",MAXVAL(bremsstrahlungtable_gang_Phi0),MAXLOC(bremsstrahlungtable_gang_Phi0) + + write(*,*) "gang min : ",MINVAL(bremsstrahlungtable_gang_Phi0),MINLOC(bremsstrahlungtable_gang_Phi0) +!~ write(*,*) bremsstrahlungtable_gang_Phi0(55,80,:,8,3,1,2) +!~ write(*,*) +!~ write(*,*) bremsstrahlungtable_gang_Phi0(55,79,:,8,3,1,2) +!~ write(*,*) +!~ write(*,*) bremsstrahlungtable_gang_Phi0(55,1,:,8,3,1,2) + + endif + + !must close h5 files, check for error if (cerror.ne.0) then diff --git a/src/nulib.F90 b/src/nulib.F90 index f4efece..5554d24 100644 --- a/src/nulib.F90 +++ b/src/nulib.F90 @@ -1,6 +1,6 @@ !-*-f90-*- module nulib - + use hdf5 implicit none !variables initialization @@ -24,6 +24,7 @@ module nulib real*8, allocatable,dimension(:) :: bin_widths ! MeV, energy width of bin real*8, allocatable,dimension(:) :: bin_bottom ! MeV, energy at bottom of bin real*8, allocatable,dimension(:) :: bin_top ! MeV, energy at top of bin + real*8, allocatable,dimension(:,:,:,:) :: gang_table ! MeV, energy at top of bin !EOS table character*200 :: eos_filename @@ -76,7 +77,14 @@ module nulib !tabulated weak rate table bounds real*8, dimension(4) :: table_bounds - + + !Gang guo table + integer(HSIZE_T) :: dims4_gang(4) + integer :: error_gang,rank_gang,cerror_gang + integer(HID_T) :: file_id_gang,dset_id_gang,dspace_id_gang + INTEGER(HID_T) :: dataspace ! + INTEGER(HSIZE_T), DIMENSION(4) :: dimsr, maxdimsr + !special terms real*8 :: adhoc_nux_factor = 0.0d0 @@ -160,6 +168,43 @@ subroutine initialize_nulib(neutrino_scheme_in,number_species_in,number_groups_i if (add_anutau_kernel_epannihil.and.add_anutau_emission_epannihil) then stop "initialize_nulib: you can't both thermally emit nu with/without a kernel: anutau" endif + if (add_nue_kernel_bremsstrahlung.and.add_nue_emission_bremsstrahlung) then + stop "initialize_nulib: you can't both thermally emit nu with/without a kernel: nue" + endif + if (add_anue_kernel_bremsstrahlung.and.add_anue_emission_bremsstrahlung) then + stop "initialize_nulib: you can't both thermally emit nu with/without a kernel: anue" + endif + if (add_numu_kernel_bremsstrahlung.and.add_numu_emission_bremsstrahlung) then + stop "initialize_nulib: you can't both thermally emit nu with/without a kernel: numu" + endif + if (add_anumu_kernel_bremsstrahlung.and.add_anumu_emission_bremsstrahlung) then + stop "initialize_nulib: you can't both thermally emit nu with/without a kernel: anumu" + endif + if (add_nutau_kernel_bremsstrahlung.and.add_nutau_emission_bremsstrahlung) then + stop "initialize_nulib: you can't both thermally emit nu with/without a kernel: nutau" + endif + if (add_anutau_kernel_bremsstrahlung.and.add_anutau_emission_bremsstrahlung) then + stop "initialize_nulib: you can't both thermally emit nu with/without a kernel: anutau" + endif + if (add_nue_kernel_gangguo.and.add_nue_emission_bremsstrahlung) then + stop "initialize_nulib: you can't both thermally emit nu with/without a kernel: nue" + endif + if (add_anue_kernel_gangguo.and.add_anue_emission_bremsstrahlung) then + stop "initialize_nulib: you can't both thermally emit nu with/without a kernel: anue" + endif + if (add_numu_kernel_gangguo.and.add_numu_emission_bremsstrahlung) then + stop "initialize_nulib: you can't both thermally emit nu with/without a kernel: numu" + endif + if (add_anumu_kernel_gangguo.and.add_anumu_emission_bremsstrahlung) then + stop "initialize_nulib: you can't both thermally emit nu with/without a kernel: anumu" + endif + if (add_nutau_kernel_gangguo.and.add_nutau_emission_bremsstrahlung) then + stop "initialize_nulib: you can't both thermally emit nu with/without a kernel: nutau" + endif + if (add_anutau_kernel_gangguo.and.add_anutau_emission_bremsstrahlung) then + stop "initialize_nulib: you can't both thermally emit nu with/without a kernel: anutau" + endif + if (neutrino_scheme.eq.1) then if ((add_numu_Iscattering_electrons.neqv.add_anumu_Iscattering_electrons).or. & @@ -197,9 +242,9 @@ subroutine initialize_nulib(neutrino_scheme_in,number_species_in,number_groups_i (add_numu_emission_epannihil.neqv.add_anutau_emission_epannihil)) then stop "With neutrino scheme 1, all 4 heavy lepton epannihils must be the same" endif - if ((add_numu_emission_NNBrems.neqv.add_anumu_emission_NNBrems).or. & - (add_numu_emission_NNBrems.neqv.add_nutau_emission_NNBrems).or. & - (add_numu_emission_NNBrems.neqv.add_anutau_emission_NNBrems)) then + if ((add_numu_emission_bremsstrahlung.neqv.add_anumu_emission_bremsstrahlung).or. & + (add_numu_emission_bremsstrahlung.neqv.add_nutau_emission_bremsstrahlung).or. & + (add_numu_emission_bremsstrahlung.neqv.add_anutau_emission_bremsstrahlung)) then stop "With neutrino scheme 1, all 4 heavy lepton NBrems must be the same" endif if ((add_numu_kernel_epannihil.neqv.add_anumu_kernel_epannihil).or. & @@ -236,8 +281,8 @@ subroutine initialize_nulib(neutrino_scheme_in,number_species_in,number_groups_i (add_anumu_emission_epannihil.neqv.add_anutau_emission_epannihil)) then stop "With neutrino scheme 2, each heavy lepton nu/anu epannihils must be the same" endif - if ((add_numu_emission_NNBrems.neqv.add_nutau_emission_NNBrems).or. & - (add_anumu_emission_NNBrems.neqv.add_anutau_emission_NNBrems)) then + if ((add_numu_emission_bremsstrahlung.neqv.add_nutau_emission_bremsstrahlung).or. & + (add_anumu_emission_bremsstrahlung.neqv.add_anutau_emission_bremsstrahlung)) then stop "With neutrino scheme 2, each heavy lepton nu/anu NBrems must be the same" endif if ((add_numu_kernel_epannihil.neqv.add_nutau_kernel_epannihil).or. & @@ -313,6 +358,64 @@ subroutine initialize_nulib(neutrino_scheme_in,number_species_in,number_groups_i call GaussLegendreQuadrature_weights_and_roots(64,GPQ_n64_roots,GPQ_n64_weights) call GaussLegendreQuadrature_weights_and_roots(128,GPQ_n128_roots,GPQ_n128_weights) + if (add_anutau_kernel_gangguo .or. add_nutau_kernel_gangguo .or. & + add_anumu_kernel_gangguo .or. add_numu_kernel_gangguo) then + + + !open HDF5 file, given filename + call h5open_f(error_gang) + + if (error_gang.ne.0) then + stop "Error reading in gang file" + endif + + + call h5fopen_f("Interpolated_data.h5", & + H5F_ACC_RDONLY_F,file_id_gang,error_gang) + + if (error_gang.ne.0) then + write(*,*) trim(adjustl("Interpolated_data.h5")) + stop "Error reading in gang table" + endif + rank_gang =4 + dims4_gang(1) = 25 + dims4_gang(2) = 25 + dims4_gang(3) = 26 + dims4_gang(4) = 40 + + call h5dopen_f(file_id_gang, "kernel", dset_id_gang, error_gang) + write(*,*) file_id_gang,dset_id_gang,dims4_gang + + CALL h5dget_space_f(dset_id_gang, dataspace, error_gang) + CALL h5sget_simple_extent_dims_f(dataspace, dimsr, maxdimsr, error_gang) + write(*,*) "arg", dimsr + + allocate(gang_table(dims4_gang(4),dims4_gang(3),dims4_gang(2),dims4_gang(1))) +!~ allocate(gang_table(dimsr(1),dimsr(2),dimsr(3),dimsr(4))) +!~ + + + call h5dread_f(dset_id_gang, H5T_NATIVE_DOUBLE,gang_table, dimsr, error_gang) + + call h5dclose_f(dset_id_gang, error_gang) + + cerror_gang = cerror_gang + error_gang + + if (cerror_gang.ne.0) then + write (*,*) "We have errors on reading Gang file", cerror_gang + stop + endif + + + call h5fclose_f(file_id_gang,error_gang) + call h5close_f(error_gang) + + + + endif + + + end subroutine initialize_nulib subroutine single_point_return_all(eos_variables, & @@ -527,9 +630,13 @@ subroutine single_point_return_all(eos_variables, & endif delta = delta / scattering_opacity !pull out the weights to get the opacity-weighted average + + !!! now calculations for non thermal processes temporary_spectra = 0.0d0 call return_emissivity_spectra_given_neutrino_scheme(temporary_spectra,eos_variables) + + if(debug) then write(*,*) "debug #3: emissivities and black bodies for species 1:", & temporary_spectra(1,1:number_groups),"bb",blackbody_spectra(1,1:number_groups) @@ -635,7 +742,8 @@ subroutine single_point_return_all(eos_variables, & end subroutine single_point_return_all - + + !calcualates the expansion of the scattering kernal, This is !\Phi_{out}_{0/1}, it does not contain the !exp(-\beta(\omega-\omega^\prime)), see make_table_example for symmetries @@ -998,6 +1106,439 @@ subroutine single_epannihil_kernel_point_return_all(iin,eta,temperature, & end subroutine single_epannihil_kernel_point_return_all + !calcualates the expansion of the bremsstrahlung kernal, These are + !|phi^{prod/annihl}_{0/1} from Hannestad & Raffelt 1998 + subroutine single_bremsstrahlung_kernel_point_return_all_Hannestad(iin,n_N,temperature, & + Phi0s,neutrino_local_scheme) + + implicit none + + !input + real*8, dimension(:,:,:), intent(out) :: Phi0s !species,other neutrino's energy, prod/annihilation + integer, intent(in) :: iin !this neutrinos energy + real*8, intent(in) :: n_N ! cm-3 + real*8, intent(in) :: temperature ! MeV + integer, intent(in) :: neutrino_local_scheme + + !neutrino species + !1 = electron neutrino !local scheme 1,2,3 + !2 = electron anti-neutrino !scheme 1,2,3 + !3 = muon neutrino !scheme 1,2,3 + !4 = muon anti-neutrino !scheme 2,3 + !5 = tau neutrino !scheme 3 + !6 = tau anti-neutrino !scheme 3 + + !x neutrino = (3+4+5+6) !scheme 1 + !y neutrino = (3+5) !scheme 2 + !z anti-neutrino = (4+6) !scheme 2 + + !local + integer :: ns,ng + integer :: number_local_species + real*8 :: nu_energy_x ! adim + real*8 :: nuother_energy_x ! adim + + !functions + real*8 :: bremsstrahlung_Phi0_Hannestad + + if (neutrino_local_scheme.ne.neutrino_scheme) then + write(*,*) neutrino_local_scheme, neutrino_scheme + stop "you are requesting different schemes" + endif + + if (neutrino_scheme.eq.1) then + number_local_species = 3 + else if (neutrino_scheme.eq.2) then + number_local_species = 4 + else if (neutrino_scheme.eq.3) then + number_local_species = 6 + else + stop "single_bremsstrahlung_kernel_return_all:incorrect neutrino scheme" + endif + + if (size(Phi0s,1).ne.number_local_species) then + stop "single_bremsstrahlung_kernel_return_all:provided array has wrong number of species " + endif + if (size(Phi0s,2).ne.number_groups) then + stop "single_bremsstrahlung_kernel_return_all:provided array has wrong number of groups" + endif + if (size(Phi0s,3).ne.2) then + stop "single_bremsstrahlung_kernel_return_all:provided array has wrong number of kernels" + endif + + nu_energy_x = energies(iin)/temperature ! adim + + Phi0s = 0.0d0 + + do ng=1,number_groups + nuother_energy_x = energies(ng)/temperature !adim + !electron neutrinos + if (add_nue_kernel_bremsstrahlung) then + Phi0s(1,ng,1) = bremsstrahlung_Phi0_Hannestad(nu_energy_x,nuother_energy_x,temperature,n_N,1,0)!production of Phi_0 + Phi0s(1,ng,2) = bremsstrahlung_Phi0_Hannestad(nu_energy_x,nuother_energy_x,temperature,n_N,1,1)!annihilation of Phi_0 + + endif + + !electron antineutrinos + if (add_anue_kernel_bremsstrahlung) then + Phi0s(2,ng,1) = bremsstrahlung_Phi0_Hannestad(nu_energy_x,nuother_energy_x,temperature,n_N,2,0)!production of Phi_0 + Phi0s(2,ng,2) = bremsstrahlung_Phi0_Hannestad(nu_energy_x,nuother_energy_x,temperature,n_N,2,1)!annihilation of Phi_0 + + endif + + if (number_local_species.eq.3) then + !already ensure that all 4 are the all included or not + if (add_numu_kernel_bremsstrahlung) then + Phi0s(3,ng,1) = (bremsstrahlung_Phi0_Hannestad(nu_energy_x,nuother_energy_x,temperature,n_N,3,0) + & + bremsstrahlung_Phi0_Hannestad(nu_energy_x,nuother_energy_x,temperature,n_N,4,0) + & + bremsstrahlung_Phi0_Hannestad(nu_energy_x,nuother_energy_x,temperature,n_N,5,0) + & + bremsstrahlung_Phi0_Hannestad(nu_energy_x,nuother_energy_x,temperature,n_N,6,0))/4.0d0 !production of Phi_0 + Phi0s(3,ng,2) = (bremsstrahlung_Phi0_Hannestad(nu_energy_x,nuother_energy_x,temperature,n_N,3,1) + & + bremsstrahlung_Phi0_Hannestad(nu_energy_x,nuother_energy_x,temperature,n_N,4,1) + & + bremsstrahlung_Phi0_Hannestad(nu_energy_x,nuother_energy_x,temperature,n_N,5,1) + & + bremsstrahlung_Phi0_Hannestad(nu_energy_x,nuother_energy_x,temperature,n_N,6,1))/4.0d0 !annihilation of Phi_0 + endif + + else if (number_local_species.eq.4) then + + !already ensure that mu and tau the same + if (add_numu_kernel_bremsstrahlung) then + Phi0s(3,ng,1) = (bremsstrahlung_Phi0_Hannestad(nu_energy_x,nuother_energy_x,temperature,n_N,3,0) + & + bremsstrahlung_Phi0_Hannestad(nu_energy_x,nuother_energy_x,temperature,n_N,5,0))/2.0d0 !production of Phi_0 + Phi0s(3,ng,2) = (bremsstrahlung_Phi0_Hannestad(nu_energy_x,nuother_energy_x,temperature,n_N,3,1) + & + bremsstrahlung_Phi0_Hannestad(nu_energy_x,nuother_energy_x,temperature,n_N,5,1))/2.0d0 !annihilation of Phi_0 + + endif + + !already ensure that amu and atau the same + if (add_anumu_kernel_bremsstrahlung) then + Phi0s(3,ng,1) = (bremsstrahlung_Phi0_Hannestad(nu_energy_x,nuother_energy_x,temperature,n_N,4,0) + & + bremsstrahlung_Phi0_Hannestad(nu_energy_x,nuother_energy_x,temperature,n_N,6,0))/2.0d0 !production of Phi_0 + Phi0s(3,ng,2) = (bremsstrahlung_Phi0_Hannestad(nu_energy_x,nuother_energy_x,temperature,n_N,4,1) + & + bremsstrahlung_Phi0_Hannestad(nu_energy_x,nuother_energy_x,temperature,n_N,6,1))/2.0d0 !annihilation of Phi_0 + endif + + else if (number_local_species.eq.6) then + + if (add_numu_Iscattering_electrons) then + Phi0s(3,ng,1) = bremsstrahlung_Phi0_Hannestad(nu_energy_x,nuother_energy_x,temperature,n_N,3,0) !production of Phi_0 + Phi0s(3,ng,2) = bremsstrahlung_Phi0_Hannestad(nu_energy_x,nuother_energy_x,temperature,n_N,3,1) !annihilation of Phi_0 + endif + + if (add_anumu_Iscattering_electrons) then + Phi0s(4,ng,1) = bremsstrahlung_Phi0_Hannestad(nu_energy_x,nuother_energy_x,temperature,n_N,4,0) !production of Phi_0 + Phi0s(4,ng,2) = bremsstrahlung_Phi0_Hannestad(nu_energy_x,nuother_energy_x,temperature,n_N,4,1) !annihilation of Phi_0 + endif + + if (add_nutau_Iscattering_electrons) then + Phi0s(5,ng,1) = bremsstrahlung_Phi0_Hannestad(nu_energy_x,nuother_energy_x,temperature,n_N,5,0) !production of Phi_0 + Phi0s(5,ng,2) = bremsstrahlung_Phi0_Hannestad(nu_energy_x,nuother_energy_x,temperature,n_N,5,1) !annihilation of Phi_0 + endif + + if (add_anutau_Iscattering_electrons) then + Phi0s(6,ng,1) = bremsstrahlung_Phi0_Hannestad(nu_energy_x,nuother_energy_x,temperature,n_N,6,0) !production of Phi_0 + Phi0s(6,ng,2) = bremsstrahlung_Phi0_Hannestad(nu_energy_x,nuother_energy_x,temperature,n_N,6,1) !annihilation of Phi_0 + endif + + else + stop "shouldn't be here" + endif + + enddo + + + end subroutine single_bremsstrahlung_kernel_point_return_all_Hannestad + + !calcualates the expansion of the bremsstrahlung kernal, These are + !|phi^{prod/annihl}_{0/1} from Kuroda et al. 2016 modified via s_kl from Raffelt & Seckel 1995 + subroutine single_bremsstrahlung_kernel_point_return_all_Kuroda(iin,n_N,temperature, & + Phi0s,neutrino_local_scheme) + + implicit none + + !input + real*8, dimension(:,:,:), intent(out) :: Phi0s !species,other neutrino's energy, prod/annihilation + integer, intent(in) :: iin !this neutrinos energy + real*8, intent(in) :: n_N ! cm-3 + real*8, intent(in) :: temperature ! MeV + integer, intent(in) :: neutrino_local_scheme + + !neutrino species + !1 = electron neutrino !local scheme 1,2,3 + !2 = electron anti-neutrino !scheme 1,2,3 + !3 = muon neutrino !scheme 1,2,3 + !4 = muon anti-neutrino !scheme 2,3 + !5 = tau neutrino !scheme 3 + !6 = tau anti-neutrino !scheme 3 + + !x neutrino = (3+4+5+6) !scheme 1 + !y neutrino = (3+5) !scheme 2 + !z anti-neutrino = (4+6) !scheme 2 + + !local + integer :: ns,ng + integer :: number_local_species + real*8 :: nu_energy_x ! adim + real*8 :: nuother_energy_x ! adim + + !functions + real*8 :: bremsstrahlung_Phi0_Kuroda + + if (neutrino_local_scheme.ne.neutrino_scheme) then + write(*,*) neutrino_local_scheme, neutrino_scheme + stop "you are requesting different schemes" + endif + + if (neutrino_scheme.eq.1) then + number_local_species = 3 + else if (neutrino_scheme.eq.2) then + number_local_species = 4 + else if (neutrino_scheme.eq.3) then + number_local_species = 6 + else + stop "single_bremsstrahlung_kernel_return_all_Kuro:incorrect neutrino scheme" + endif + + if (size(Phi0s,1).ne.number_local_species) then + stop "single_bremsstrahlung_kernel_return_all_Kuro:provided array has wrong number of species" + endif + if (size(Phi0s,2).ne.number_groups) then + stop "single_bremsstrahlung_kernel_return_all_Kuro:provided array has wrong number of groups" + endif + if (size(Phi0s,3).ne.2) then + stop "single_bremsstrahlung_kernel_return_all_Kuro:provided array has wrong number of kernels" + endif + + nu_energy_x = energies(iin)/temperature ! adim + + Phi0s = 0.0d0 + + do ng=1,number_groups + nuother_energy_x = energies(ng)/temperature !adim + !electron neutrinos + if (add_nue_kernel_bremsstrahlung) then + Phi0s(1,ng,1) = bremsstrahlung_Phi0_Kuroda(nu_energy_x,nuother_energy_x,temperature,n_N,1,0)!production of Phi_0 + Phi0s(1,ng,2) = bremsstrahlung_Phi0_Kuroda(nu_energy_x,nuother_energy_x,temperature,n_N,1,1)!annihilation of Phi_0 + + endif + + !electron antineutrinos + if (add_anue_kernel_bremsstrahlung) then + Phi0s(2,ng,1) = bremsstrahlung_Phi0_Kuroda(nu_energy_x,nuother_energy_x,temperature,n_N,2,0)!production of Phi_0 + Phi0s(2,ng,2) = bremsstrahlung_Phi0_Kuroda(nu_energy_x,nuother_energy_x,temperature,n_N,2,1)!annihilation of Phi_0 + + endif + + if (number_local_species.eq.3) then + !already ensure that all 4 are the all included or not + if (add_numu_kernel_bremsstrahlung) then + Phi0s(3,ng,1) = (bremsstrahlung_Phi0_Kuroda(nu_energy_x,nuother_energy_x,temperature,n_N,3,0) + & + bremsstrahlung_Phi0_Kuroda(nu_energy_x,nuother_energy_x,temperature,n_N,4,0) + & + bremsstrahlung_Phi0_Kuroda(nu_energy_x,nuother_energy_x,temperature,n_N,5,0) + & + bremsstrahlung_Phi0_Kuroda(nu_energy_x,nuother_energy_x,temperature,n_N,6,0))/4.0d0 !production of Phi_0 + Phi0s(3,ng,2) = (bremsstrahlung_Phi0_Kuroda(nu_energy_x,nuother_energy_x,temperature,n_N,3,1) + & + bremsstrahlung_Phi0_Kuroda(nu_energy_x,nuother_energy_x,temperature,n_N,4,1) + & + bremsstrahlung_Phi0_Kuroda(nu_energy_x,nuother_energy_x,temperature,n_N,5,1) + & + bremsstrahlung_Phi0_Kuroda(nu_energy_x,nuother_energy_x,temperature,n_N,6,1))/4.0d0 !annihilation of Phi_0 + endif + + else if (number_local_species.eq.4) then + + !already ensure that mu and tau the same + if (add_numu_kernel_bremsstrahlung) then + Phi0s(3,ng,1) = (bremsstrahlung_Phi0_Kuroda(nu_energy_x,nuother_energy_x,temperature,n_N,3,0) + & + bremsstrahlung_Phi0_Kuroda(nu_energy_x,nuother_energy_x,temperature,n_N,5,0))/2.0d0 !production of Phi_0 + Phi0s(3,ng,2) = (bremsstrahlung_Phi0_Kuroda(nu_energy_x,nuother_energy_x,temperature,n_N,3,1) + & + bremsstrahlung_Phi0_Kuroda(nu_energy_x,nuother_energy_x,temperature,n_N,5,1))/2.0d0 !annihilation of Phi_0 + + endif + + !already ensure that amu and atau the same + if (add_anumu_kernel_bremsstrahlung) then + Phi0s(3,ng,1) = (bremsstrahlung_Phi0_Kuroda(nu_energy_x,nuother_energy_x,temperature,n_N,4,0) + & + bremsstrahlung_Phi0_Kuroda(nu_energy_x,nuother_energy_x,temperature,n_N,6,0))/2.0d0 !production of Phi_0 + Phi0s(3,ng,2) = (bremsstrahlung_Phi0_Kuroda(nu_energy_x,nuother_energy_x,temperature,n_N,4,1) + & + bremsstrahlung_Phi0_Kuroda(nu_energy_x,nuother_energy_x,temperature,n_N,6,1))/2.0d0 !annihilation of Phi_0 + endif + + else if (number_local_species.eq.6) then + + if (add_numu_kernel_bremsstrahlung) then + Phi0s(3,ng,1) = bremsstrahlung_Phi0_Kuroda(nu_energy_x,nuother_energy_x,temperature,n_N,3,0) !production of Phi_0 + Phi0s(3,ng,2) = bremsstrahlung_Phi0_Kuroda(nu_energy_x,nuother_energy_x,temperature,n_N,3,1) !annihilation of Phi_0 + endif + + if (add_anumu_kernel_bremsstrahlung) then + Phi0s(4,ng,1) = bremsstrahlung_Phi0_Kuroda(nu_energy_x,nuother_energy_x,temperature,n_N,4,0) !production of Phi_0 + Phi0s(4,ng,2) = bremsstrahlung_Phi0_Kuroda(nu_energy_x,nuother_energy_x,temperature,n_N,4,1) !annihilation of Phi_0 + endif + + if (add_nutau_kernel_bremsstrahlung) then + Phi0s(5,ng,1) = bremsstrahlung_Phi0_Kuroda(nu_energy_x,nuother_energy_x,temperature,n_N,5,0) !production of Phi_0 + Phi0s(5,ng,2) = bremsstrahlung_Phi0_Kuroda(nu_energy_x,nuother_energy_x,temperature,n_N,5,1) !annihilation of Phi_0 + endif + + if (add_anutau_kernel_bremsstrahlung) then + Phi0s(6,ng,1) = bremsstrahlung_Phi0_Kuroda(nu_energy_x,nuother_energy_x,temperature,n_N,6,0) !production of Phi_0 + Phi0s(6,ng,2) = bremsstrahlung_Phi0_Kuroda(nu_energy_x,nuother_energy_x,temperature,n_N,6,1) !annihilation of Phi_0 + endif + + else + stop "shouldn't be here" + endif + + enddo + + + + end subroutine single_bremsstrahlung_kernel_point_return_all_Kuroda + + + !calcualates the expansion of the bremsstrahlung kernal, These are + !|phi^{prod/annihl}_{0/1} from gang et al. 2019 + subroutine single_bremsstrahlung_kernel_point_return_all_gang(iin,n_N,Ye,temperature, & + Phi0s,neutrino_local_scheme) + + implicit none + + !input + real*8, dimension(:,:,:), intent(out) :: Phi0s !species,other neutrino's energy, prod/annihilation + integer, intent(in) :: iin !this neutrinos energy + real*8, intent(in) :: n_N ! cm-3 + real*8, intent(in) :: Ye ! cm-3 + real*8, intent(in) :: temperature ! MeV + integer, intent(in) :: neutrino_local_scheme + + !neutrino species + !1 = electron neutrino !local scheme 1,2,3 + !2 = electron anti-neutrino !scheme 1,2,3 + !3 = muon neutrino !scheme 1,2,3 + !4 = muon anti-neutrino !scheme 2,3 + !5 = tau neutrino !scheme 3 + !6 = tau anti-neutrino !scheme 3 + + !x neutrino = (3+4+5+6) !scheme 1 + !y neutrino = (3+5) !scheme 2 + !z anti-neutrino = (4+6) !scheme 2 + + !local + integer :: ns,ng + integer :: number_local_species + real*8 :: nu_energy_x ! adim + real*8 :: nuother_energy_x ! adim + + !functions + real*8 :: bremsstrahlung_Phi0_gang + + if (neutrino_local_scheme.ne.neutrino_scheme) then + write(*,*) neutrino_local_scheme, neutrino_scheme + stop "you are requesting different schemes" + endif + + if (neutrino_scheme.eq.1) then + number_local_species = 3 + else if (neutrino_scheme.eq.2) then + number_local_species = 4 + else if (neutrino_scheme.eq.3) then + number_local_species = 6 + else + stop "single_bremsstrahlung_kernel_return_all_Kuro:incorrect neutrino scheme" + endif + + if (size(Phi0s,1).ne.number_local_species) then + stop "single_bremsstrahlung_kernel_return_all_Kuro:provided array has wrong number of species" + endif + if (size(Phi0s,2).ne.number_groups) then + stop "single_bremsstrahlung_kernel_return_all_Kuro:provided array has wrong number of groups" + endif + if (size(Phi0s,3).ne.2) then + stop "single_bremsstrahlung_kernel_return_all_Kuro:provided array has wrong number of kernels" + endif + + nu_energy_x = energies(iin)/temperature ! adim +!~ write(*,*) number_groups + Phi0s = 0.0d0 + + do ng=1,number_groups + nuother_energy_x = energies(ng)/temperature !adim + !electron neutrinos + if (add_nue_kernel_gangguo) then + Phi0s(1,ng,1) = bremsstrahlung_Phi0_gang(nu_energy_x,nuother_energy_x,temperature,n_N,Ye,1,0)!production of Phi_0 + Phi0s(1,ng,2) = bremsstrahlung_Phi0_gang(nu_energy_x,nuother_energy_x,temperature,n_N,Ye,1,1)!annihilation of Phi_0 + + endif + + !electron antineutrinos + if (add_anue_kernel_gangguo) then + Phi0s(2,ng,1) = bremsstrahlung_Phi0_gang(nu_energy_x,nuother_energy_x,temperature,n_N,Ye,2,0)!production of Phi_0 + Phi0s(2,ng,2) = bremsstrahlung_Phi0_gang(nu_energy_x,nuother_energy_x,temperature,n_N,Ye,2,1)!annihilation of Phi_0 + + endif + + if (number_local_species.eq.3) then + !already ensure that all 4 are the all included or not + if (add_numu_kernel_gangguo) then + Phi0s(3,ng,1) = (bremsstrahlung_Phi0_gang(nu_energy_x,nuother_energy_x,temperature,n_N,Ye,3,0) + & + bremsstrahlung_Phi0_gang(nu_energy_x,nuother_energy_x,temperature,n_N,Ye,4,0) + & + bremsstrahlung_Phi0_gang(nu_energy_x,nuother_energy_x,temperature,n_N,Ye,5,0) + & + bremsstrahlung_Phi0_gang(nu_energy_x,nuother_energy_x,temperature,n_N,Ye,6,0))/4.0d0 !production of Phi_0 + Phi0s(3,ng,2) = (bremsstrahlung_Phi0_gang(nu_energy_x,nuother_energy_x,temperature,n_N,Ye,3,1) + & + bremsstrahlung_Phi0_gang(nu_energy_x,nuother_energy_x,temperature,n_N,Ye,4,1) + & + bremsstrahlung_Phi0_gang(nu_energy_x,nuother_energy_x,temperature,n_N,Ye,5,1) + & + bremsstrahlung_Phi0_gang(nu_energy_x,nuother_energy_x,temperature,n_N,Ye,6,1))/4.0d0 !annihilation of Phi_0 + endif + + else if (number_local_species.eq.4) then + + !already ensure that mu and tau the same + if (add_numu_kernel_gangguo) then + Phi0s(3,ng,1) = (bremsstrahlung_Phi0_gang(nu_energy_x,nuother_energy_x,temperature,n_N,Ye,3,0) + & + bremsstrahlung_Phi0_gang(nu_energy_x,nuother_energy_x,temperature,n_N,Ye,5,0))/2.0d0 !production of Phi_0 + Phi0s(3,ng,2) = (bremsstrahlung_Phi0_gang(nu_energy_x,nuother_energy_x,temperature,n_N,Ye,3,1) + & + bremsstrahlung_Phi0_gang(nu_energy_x,nuother_energy_x,temperature,n_N,Ye,5,1))/2.0d0 !annihilation of Phi_0 + + endif + + !already ensure that amu and atau the same + if (add_anumu_kernel_gangguo) then + Phi0s(3,ng,1) = (bremsstrahlung_Phi0_gang(nu_energy_x,nuother_energy_x,temperature,n_N,Ye,4,0) + & + bremsstrahlung_Phi0_gang(nu_energy_x,nuother_energy_x,temperature,n_N,Ye,6,0))/2.0d0 !production of Phi_0 + Phi0s(3,ng,2) = (bremsstrahlung_Phi0_gang(nu_energy_x,nuother_energy_x,temperature,n_N,Ye,4,1) + & + bremsstrahlung_Phi0_gang(nu_energy_x,nuother_energy_x,temperature,n_N,Ye,6,1))/2.0d0 !annihilation of Phi_0 + endif + + else if (number_local_species.eq.6) then + + if (add_numu_kernel_gangguo) then + Phi0s(3,ng,1) = bremsstrahlung_Phi0_gang(nu_energy_x,nuother_energy_x,temperature,n_N,Ye,3,0) !production of Phi_0 + Phi0s(3,ng,2) = bremsstrahlung_Phi0_gang(nu_energy_x,nuother_energy_x,temperature,n_N,Ye,3,1) !annihilation of Phi_0 + endif + + if (add_anumu_kernel_gangguo) then + Phi0s(4,ng,1) = bremsstrahlung_Phi0_gang(nu_energy_x,nuother_energy_x,temperature,n_N,Ye,4,0) !production of Phi_0 + Phi0s(4,ng,2) = bremsstrahlung_Phi0_gang(nu_energy_x,nuother_energy_x,temperature,n_N,Ye,4,1) !annihilation of Phi_0 + endif + + if (add_nutau_kernel_gangguo) then + Phi0s(5,ng,1) = bremsstrahlung_Phi0_gang(nu_energy_x,nuother_energy_x,temperature,n_N,Ye,5,0) !production of Phi_0 + Phi0s(5,ng,2) = bremsstrahlung_Phi0_gang(nu_energy_x,nuother_energy_x,temperature,n_N,Ye,5,1) !annihilation of Phi_0 + endif + + if (add_anutau_kernel_gangguo) then + Phi0s(6,ng,1) = bremsstrahlung_Phi0_gang(nu_energy_x,nuother_energy_x,temperature,n_N,Ye,6,0) !production of Phi_0 + Phi0s(6,ng,2) = bremsstrahlung_Phi0_gang(nu_energy_x,nuother_energy_x,temperature,n_N,Ye,6,1) !annihilation of Phi_0 + endif + + else + stop "shouldn't be here" + endif + + enddo + + + + end subroutine single_bremsstrahlung_kernel_point_return_all_gang + end module nulib diff --git a/src/point_example.F90 b/src/point_example.F90 index 934964e..1b186a6 100644 --- a/src/point_example.F90 +++ b/src/point_example.F90 @@ -3,6 +3,7 @@ program point_example use nulib use inputparser + #if NUCLEI_HEMPEL use nuclei_hempel, only : set_up_Hempel #endif @@ -30,9 +31,9 @@ program point_example integer :: mypoint_number_output_species = 3 !number of energy groups - integer :: mypoint_number_groups = 24 + integer :: mypoint_number_groups = 18 !18 - character*200 :: parameters = "/projects/ceclub/gr1dnulib/GitHub/NuLib/parameters" + character*200 :: parameters = "./parameters" !local variables real*8, allocatable,dimension(:,:) :: local_emissivity @@ -40,15 +41,50 @@ program point_example real*8, allocatable,dimension(:,:) :: local_scatopacity real*8, allocatable,dimension(:,:) :: local_delta real*8, allocatable,dimension(:,:) :: local_Phi0, local_Phi1 + real*8, allocatable,dimension(:,:,:) :: local_Phi0_epannihil, local_Phi1_epannihil + real*8, allocatable,dimension(:,:,:) :: local_Phi0_bremsstrahlung + real*8, allocatable,dimension(:,:,:) :: local_Phi0_bremsstrahlung2 real*8, allocatable,dimension(:,:) :: blackbody_spectra real*8, allocatable,dimension(:) :: eos_variables + real*8,allocatable, dimension(:) :: Phi0_brem,Phi0_brem2 + real*8,allocatable, dimension(:) :: Phi0_brem_ann,Phi0_brem2_ann real*8 :: matter_prs,matter_ent,matter_cs2,matter_dedt,matter_dpderho,matter_dpdrhoe + real*8 :: single_neutrino_emissivity_from_NNBrem_given_energypoint + real*8 :: single_neutrino_emissivity_from_epannihil_given_energypoint integer :: keytemp,keyerr real*8 :: precision = 1.0d-10 real*8 :: xrho, xtemp, xye - integer i,j - real*8 dxfac,mindx + real*8 :: n_N,eta_star + integer :: i,j,inde + real*8 :: dxfac,mindx + real*8 :: Q,Q_ann,Q2,Q2_ann,M,dE_dEdt,dE_dEdt2,dE_dEdt3,k + real*8 :: Kuro_inv_lamb,Hann_inv_lamb,epan_inv_lamb + real*8 :: dn_dEdt,dn_dEdt2,dn_dEdt3,ratio,ratio_2,dn_dEdt_burrows + real*8 :: dE_dEdt_2,dE_dEdt2_2,dE_dEdt3_2,dE_dEdt_burrows + real*8 :: Bremsstrahlung_Phi0_Hannestad + real*8 :: Bremsstrahlung_Phi0_Kuroda + real*8 :: epannihil_Phi_Bruenn +!~ real*8 :: find_s + real*8,dimension(100) :: integral1,integral2,integral3,temp_array,dens_array,ener_array + real*8:: integral + real*8 :: J_1,J_1_bar,phi_a,phi_p,energy,energy_2 + real*8,dimension(100) :: M1_mom,M1_mom_inv + + + !fermi function + real*8 :: fermidirac,fermidirac_dimensionless + real*8 :: x + real*8 :: find_s,find_s2 + real*8 :: find_g + real*8 :: nue_absorption_on_n + real*8 :: anue_absorption_on_p + + real*8,parameter :: nulib_emissivity_gf = 5.59424238d-55/(6.77140812d-06**3*2.03001708d+05) + real*8,parameter :: nulib_opacity_gf = 1.0d0/6.77140812d-06 + real*8,parameter :: nulib_energy_gf = 1.60217733d-6*5.59424238d-55 + real*8,parameter :: nulib_kernel_gf = 6.77140812d-06**3/2.03001708d+05 + !allocate the arrays for the point values allocate(local_emissivity(mypoint_number_output_species,mypoint_number_groups)) allocate(local_absopacity(mypoint_number_output_species,mypoint_number_groups)) @@ -70,17 +106,18 @@ program point_example ! m_ref = m_n !for LS220 !example point - xrho = 1.0d12 !g/cm^3 - xtemp = 1.5d0 !MeV - xye = 0.35d0 !dimensionless - + + xrho = 4.0d14 !g/cm^3 + xtemp = 12.0d0 !MeV + xye = 0.1d0! !dimensionless + !set up energies bins - do_integrated_BB_and_emissivity = .false. - mindx = 1.0d0 + do_integrated_BB_and_emissivity = .true. + mindx = 2.0d0 bin_bottom(1) = 0.0d0 !MeV - bin_bottom(2) = 1.0d0 !MeV + bin_bottom(2) = 2.0d0 !MeV bin_bottom(3) = bin_bottom(2)+mindx - bin_bottom(number_groups) = 250.0d0 + bin_bottom(number_groups) = 200.0d0 ! MeV call nulib_series2(number_groups-1,bin_bottom(2),bin_bottom(number_groups),mindx,dxfac) do i=4,number_groups @@ -93,16 +130,17 @@ program point_example bin_widths(i) = bin_bottom(i+1)-bin_bottom(i) bin_top(i) = bin_bottom(i+1) enddo + energies(number_groups) = bin_bottom(number_groups)+bin_widths(number_groups-1)*dxfac/2.0d0 bin_widths(number_groups) = 2.0*(energies(number_groups)-bin_bottom(number_groups)) bin_top(number_groups) = bin_bottom(number_groups)+bin_widths(number_groups) - + write(*,*) energies + stop allocate(eos_variables(total_eos_variables)) eos_variables(:) = 0.0d0 eos_variables(rhoindex) = xrho eos_variables(tempindex) = xtemp eos_variables(yeindex) = xye - !! EOS stuff call set_eos_variables(eos_variables) write(*,*) "Rho: ",eos_variables(rhoindex)," g/ccm" @@ -112,6 +150,16 @@ program point_example write(*,*) "X_p: ",eos_variables(xpindex) write(*,*) "X_alpha: ",eos_variables(xaindex) + write(*,*) "eos", "xp",eos_variables(xpindex),"xn",eos_variables(xnindex), & + "mue",eos_variables(mueindex),"xa",eos_variables(xaindex), "xh",& + eos_variables(xhindex) + + write(*,*) eos_variables(xpindex)+eos_variables(xnindex)+eos_variables(xaindex)& + +eos_variables(xhindex) + + write(*,*) "emi pair 10 bin nue",single_neutrino_emissivity_from_epannihil_given_energypoint( & + 1,energies(10),eos_variables) + !calculate the full emissivities and opacities, this averages the !values based on the mypoint_neutrino_scheme this assumes detailed !balance and uses the opacities to fullying calculate the @@ -119,15 +167,20 @@ program point_example call single_point_return_all(eos_variables, & local_emissivity,local_absopacity,local_scatopacity,local_delta, & mypoint_neutrino_scheme) - - write(*,*) "Example of a single point call with returning all emissivity, absorptive opacity, and scattering opacity" - write(*,*) - do i=1,mypoint_number_output_species - do j=1,mypoint_number_groups - write(*,"(i4,i4,1P10E18.9)") i,j,energies(j),local_emissivity(i,j),local_absopacity(i,j), & - local_scatopacity(i,j),local_delta(i,j) - enddo - enddo + +write(*,*) local_emissivity(1,6), energies(6) +write(*,*) local_absopacity(1,6), energies(6) +write(*,*) nue_absorption_on_n(energies(6),eos_variables), energies(6) +write(*,*) anue_absorption_on_p(energies(6),eos_variables), energies(6) +stop +!~ write(*,*) "Example of a single point call with returning all emissivity, absorptive opacity, and scattering opacity" +!~ write(*,*) +!~ do i=1,mypoint_number_output_species +!~ do j=1,mypoint_number_groups +!~ write(*,"(i4,i4,1P10E18.9)") i,j,energies(j),eos_variables(tempindex),local_emissivity(i,j),local_absopacity(i,j), & +!~ local_scatopacity(i,j),local_delta(i,j) +!~ enddo +!~ enddo write(*,*) write(*,*) @@ -193,13 +246,58 @@ program point_example call single_Ipoint_return_all(mypoint_number_groups,eos_variables(mueindex)/eos_variables(tempindex), & eos_variables(tempindex),local_Phi0,local_Phi1,mypoint_neutrino_scheme) - write(*,*) "sample scattering kernel for energy",energies(mypoint_number_groups) +!~ write(*,*) "sample scattering kernel for energy",energies(mypoint_number_groups) do i=1,mypoint_number_output_species do j=1,mypoint_number_groups - write(*,*) energies(mypoint_number_groups),energies(j),local_Phi0(i,j),local_Phi1(i,j) +!~ write(*,*) energies(mypoint_number_groups),energies(j),local_Phi0(i,j),local_Phi1(i,j) enddo enddo endif + + + if (add_nue_kernel_epannihil.or.add_anue_kernel_epannihil.or. & + add_numu_kernel_epannihil.or.add_anumu_kernel_epannihil.or. & + add_nutau_kernel_epannihil.or.add_anutau_kernel_epannihil) then + + write(*,*) + write(*,*) + write(*,*) "epan call" + + allocate(local_Phi0_epannihil(mypoint_number_output_species,mypoint_number_groups,2)) + allocate(local_Phi1_epannihil(mypoint_number_output_species,mypoint_number_groups,2)) + + call single_epannihil_kernel_point_return_all(mypoint_number_groups, & + eos_variables(mueindex)/eos_variables(tempindex), & + eos_variables(tempindex),local_Phi0_epannihil,local_Phi1_epannihil,mypoint_neutrino_scheme) + + write(*,*) "epan absorption for nu_x",local_Phi0_epannihil(3,:,1) + + + deallocate(local_Phi0_epannihil) + deallocate(local_Phi1_epannihil) + + endif + + if (add_nue_kernel_bremsstrahlung.or.add_anue_kernel_bremsstrahlung.or. & + add_numu_kernel_bremsstrahlung.or.add_anumu_kernel_bremsstrahlung.or. & + add_nutau_kernel_bremsstrahlung.or.add_anutau_kernel_bremsstrahlung) then + + write(*,*) + write(*,*) + write(*,*) "epan call" + + allocate(local_Phi0_bremsstrahlung(mypoint_number_output_species,mypoint_number_groups,2)) + + call single_bremsstrahlung_kernel_point_return_all_Hannestad(i, & + n_N,eos_variables(tempindex),local_Phi0_bremsstrahlung,mypoint_neutrino_scheme) + + + write(*,*) "bremsstrahlung kernel absorption for nu_x", local_Phi0_bremsstrahlung + + deallocate(local_Phi0_bremsstrahlung) + endif + + end program point_example diff --git a/src/requested_interactions.inc b/src/requested_interactions.inc index 1a771be..20a0171 100644 --- a/src/requested_interactions.inc +++ b/src/requested_interactions.inc @@ -13,7 +13,7 @@ !absorptions logical :: add_nue_absorption_on_n = .true. logical :: add_anue_absorption_on_p = .true. - logical :: add_nue_absorption_on_A = .false. + logical :: add_nue_absorption_on_A = .true. logical :: add_nux_absorption_on_n_and_p = .false. !scatterings @@ -64,22 +64,42 @@ logical :: add_nutau_emission_epannihil = .true. logical :: add_anutau_emission_epannihil = .true. - logical :: add_nue_emission_NNBrems = .false. - logical :: add_anue_emission_NNBrems = .false. - logical :: add_numu_emission_NNBrems = .true. - logical :: add_anumu_emission_NNBrems = .true. - logical :: add_nutau_emission_NNBrems = .true. - logical :: add_anutau_emission_NNBrems = .true. + logical :: add_nue_emission_bremsstrahlung = .false. + logical :: add_anue_emission_bremsstrahlung = .false. + logical :: add_numu_emission_bremsstrahlung = .true. + logical :: add_anumu_emission_bremsstrahlung = .true. + logical :: add_nutau_emission_bremsstrahlung = .true. + logical :: add_anutau_emission_bremsstrahlung = .true. - logical :: add_nue_emission_weakinteraction_ecap = .true. - logical :: add_anue_emission_weakinteraction_poscap = .true. + logical :: add_nue_emission_weakinteraction_ecap = .false. + logical :: add_anue_emission_weakinteraction_poscap = .false. logical :: apply_kirchoff_to_pair_creation = .true. + logical :: separated_pair_processes = .false. + + + !kernels for full calculation of thermal processes logical :: add_nue_kernel_epannihil = .false. logical :: add_anue_kernel_epannihil = .false. logical :: add_numu_kernel_epannihil = .false. logical :: add_anumu_kernel_epannihil = .false. logical :: add_nutau_kernel_epannihil = .false. - logical :: add_anutau_kernel_epannihil = .false. \ No newline at end of file + logical :: add_anutau_kernel_epannihil = .false. + + !kernels for full calculation of thermal processes + logical :: add_nue_kernel_bremsstrahlung = .false. + logical :: add_anue_kernel_bremsstrahlung = .false. + logical :: add_numu_kernel_bremsstrahlung = .false. + logical :: add_anumu_kernel_bremsstrahlung = .false. + logical :: add_nutau_kernel_bremsstrahlung = .false. + logical :: add_anutau_kernel_bremsstrahlung = .false. + + !kernels from gang + logical :: add_nue_kernel_gangguo = .false. + logical :: add_anue_kernel_gangguo = .false. + logical :: add_numu_kernel_gangguo = .false. + logical :: add_anumu_kernel_gangguo = .false. + logical :: add_nutau_kernel_gangguo = .false. + logical :: add_anutau_kernel_gangguo = .false.