!23456789a123456789b123456789c123456789d123456789e123456789f123456789g12 ! ! Copyright D.J. Jeffery, 2006jan01. ! ! rde.f for preparing radioactive decay data. ! ! References: ! \bibitem[Browne \& Firestone(1986)]{browne1986} ! Browne, E., \& Firestone, R. B. 1986, ! Table of Radioactive Isotopes (New York: ! John Wiley \& Sons, Inc.) ! % Radioacitve isotopes and decay chains and ! % radiations including all the complex stuff like ! % X-ray, atomic electron, internal bremstrahlung X-ray. ! ! \bibitem[Huo(1992)]{huo1992} Huo, J. 1992, ! Nuclear Data Sheets, 67, 523. ! % all about isobars of atomic mass 56, ! % especially Ni56 and Co56. ! % An update of Huo~\etal 1987. ! ! \bibitem[Jeffery(1998)]{jeffery1998} Jeffery, D. J. 1998, ! astro-ph/9811356 ! % the longer LS grey $\gamma$-ray transfer procedure ! % for supernovae. ! % Effectively rejected A\&SS, but just for being ! % prolix and (in the referee's opinion) worthless. ! % But the referee did not say anything was wrong. ! % So unpublished rather than preprint is best. ! % The later implies eventual publication which ! % may not be forthcoming. ! % There is a good review of radioactive decay energy ! % quantities in Section 5. ! ! LBNL Isotopes ! Project - LUNDS Universitet ! WWW Table of Radioactive Isotopes ! http://ie.lbl.gov/toi/nuclide.asp?iZA=270056 ! Co-56 ! ! \bibitem[Metcalf et al.(2004)]{metcalf2004} Metcalf, M., Reid, J., ! \&~Cohen, M. 2004, ! Fortran 95/2003 Explained ! (Oxford: Oxford University Press) (MRC) ! % Pretty good reference book. ! !----------------------------------------------------------------------- ! program rde use numerical use physical use astronomical use rde_mod use rde_mod2 implicit none ! integer, parameter :: nelevel=200 integer, parameter :: nxray=200 ! real (kind=nprecision) :: beta_speedf real (kind=nprecision) :: con_mev_erg real (kind=nprecision) :: con_t_half_t_efold integer :: i,j,k,l,m,n integer :: idata integer :: ielevel real (kind=nprecision) :: small real (kind=nprecision) :: x real (kind=nprecision) :: xmass_np_diff ! real (kind=nprecision) :: elevel(10,nelevel) ! Energy level ! ! and beta/EC data. integer :: ixray ! Number of X-ray emissions. real (kind=nprecision) :: xray(2,nxray) ! X-ray data. ! namelist/param2/elevel,xray ! beta_speedf(x)=sqrt(1.d0-1.d0/(1.d0+x/emass_mev)**2) ! !----------------------------------------------------------------------- write(*,*) 'Before preliminary results calculated.' !----------------------------------------------------------------------- ! con_mev_erg=1.d+6*echarge*1.d+7 ! 1.d+6 for MeV to eV, 1.d+7 to ergs. con_t_half_t_efold=1.d0/log(2.d0) ! Convert half-life to e-folding time. small=1.d-10 xmass_np_diff=x_neutron_mass_mev-proton_mass_mev ! write(*,*) write(*,*) 'con_mev_erg,con_t_half_t_efold' write(*,*) con_mev_erg,con_t_half_t_efold ! !----------------------------------------------------------------------- write(*,*) 'Before data-read and evaluation of the quantities.' !----------------------------------------------------------------------- ! open(unit=1,file='rde.dat',status='old',action='read') ! call readcopy(1,0,0,'',idata) open(unit=2,file='rde_out.dat',status='unknown',action='write') open(unit=3,file='rde_out2.dat',status='unknown',action='write') i=0 do410 : do elevel(:,:)=0.d0 read(1,param,end=210) read(1,param2,end=210) i=i+1 write(*,*) write(*,*) 'Nuclide ',i ! t_efold=t_half*con_t_half_t_efold ! !----------------------------------------------------------------------- ! write(*,*) 'Before counting the energy levels and xray ', ! write(*,*) 'emissions (IB plus lines).' !----------------------------------------------------------------------- ! do420 : do j=2,nelevel ! There is always 1 level and the ground state ! ! level can fool the level counter. if(elevel(1,j) .lt. small) exit do420 end do do420 ielevel=min(j-1,nelevel)! The min function is redundant: j-1 would do. do430 : do j=1,nxray if(xray(1,j) .lt. small) exit do430 end do do430 ixray=min(j-1,nxray) ! The min function is redundant: j-1 would do. ! elevel(1,:ielevel)=elevel(1,:ielevel)*1.d-3 ! Conversion to MeV. elevel(3,:ielevel)=elevel(3,:ielevel)*1.d-2 ! Conversion to fraction. elevel(5,:ielevel)=elevel(5,:ielevel)*1.d-2 ! Conversion to fraction. qval_gamma_direct=sum( elevel(1,:ielevel) & & *(elevel(3,:ielevel)+elevel(5,:ielevel)) ) betafrac=sum(elevel(3,:ielevel)) qval_xray=sum( .001d0*xray(1,:ixray) & & *xray(2,:ixray)*.01d0) ! With conversion to MeV and ! ! fraction. ! !----------------------------------------------------------------------- ! ! The tricky part is the neutrino energy. We are assuming it escapes, ! and so we don't want to count it. ! But it is always included in the standard total disintegration energy. ! ! So we get the direct gamma fraction as above by counting up the energy ! That must come from decays from the excited daughter. ! ! We add to this, the mean beta kinetic energy per decay (not per ! beta particle) that we get from a table. ! ! We add to this, the pair annihilation energy if a beta+ is created. ! We assume that the positron is always annihilated. ! ! We neglect all atomic binding energies. See Enge-304--305. ! They are not always negligible, but they are at another order of ! accuracy than here. ! ! But we could account for them using the X-ray and Auger electron ! emission. But the Auger electron emission can be clumped with the ! internal conversion electron emission which competes with gamma-ray ! emission. Thus, one is in danger of double counting ! some energy both as gamma-rays and as internval conversion ! electrons. It's a bloody mess trying to sort it all out. ! There is also internal bremstrahlung X-ray emission. ! Better leave all this to anther day. ! But can get the data to do it in browne1986. It may be oline too. ! One should probably have xray and atelectron and photon and lepton ! categories. ! ! We also neglect nucleus recoil energy. It's probably pretty darn ! negligible since no one ever talks about it. ! ! Note in beta+ and EC, one goes from nucleus 1 and Z electrons in the ! universe to nucleus 2 and Z-1 electrons in the universe after annihilation. ! ! The total disintegration energy is atom 1 to atom 2 which also is ! Z to Z-1 electrons in the universe. So aside from binding eneryg ! everything is accurate. ! ! Now in beta- decay, one goes from nucleus 1 and Z electrons in the universe ! to nucleus 2 and Z+1 electrons in the universe. ! ! The total disintegration energy is atom 1 to atom 2 which also is ! Z to Z+1 electrons in the universe. So aside from binding eneryg ! everything is accurate. ! ! Recall the reactions: ! EC e**(-) + p**(+) --> n + neutrino ! beta+ decay p**(+) --> n + e**(+)+neutrino (Enge-302) ! beta- decay n --> p**(+) + e**(-)+antineutrino (Enge-303) !----------------------------------------------------------------------- ! qval_nu_omit=qval_gamma_direct & & +qval_beta & & +2.d0*emass_mev*betafrac*real(max(ibeta,0),nprecision) & & +qval_atomic & & +qval_xray fatomic=qval_atomic/qval_nu_omit ! Fraction of KE energy in atomic electrons. fbeta=qval_beta/qval_nu_omit ! Fraction of KE energy in positrons. flepton=fatomic+fbeta ! Fraction of KE energy in leptons fxray=qval_xray/qval_nu_omit ! Fraction of KE energy in x_rays. fgamma=1.d0-flepton-fxray ! Fraction of energy in gamma-rays fphoton=1.d0-flepton ! Fraction of energy in photons. qval_gamma=fgamma*qval_nu_omit if(betafrac .ge. small) then beta_ke=qval_beta/betafrac ! This is the mean beta KE per beta. beta_speed=beta_speedf(beta_ke) else beta_ke=0.d0 beta_speed=0.d0 end if ! c_coef=(qval_nu_omit*con_mev_erg) & & /(amu*real(ia,nprecision)*t_efold*daysec) c_coef_lepton=c_coef*(fatomic+fbeta) c_coef_photon=c_coef*(fgamma+fxray) c_coef_sun=c_coef*xmass_sun if(ia_parent .ne. 0) then t_efold_parent=t_half_parent*con_t_half_t_efold b_coef=(qval_nu_omit*con_mev_erg) & & /(amu*real(ia,nprecision) & & *(t_efold-t_efold_parent)*daysec) else t_efold_parent=0.d0 b_coef=0.d0 end if b_coef_lepton=b_coef*flepton b_coef_photon=b_coef*fphoton b_coef_sun=b_coef*xmass_sun ! write(2,*) write(2,*) 'Nuclide ',i write(2,param) ! radio(i)%inuclide=i radio(i)%aanuclide=aanuclide radio(i)%flepton=flepton radio(i)%fphoton=fphoton radio(i)%gfactor=(qval_nu_omit/(real(ia,nprecision) & & *t_efold) )*rw7 ! ! Note for hfactor it should be ia_parent, but that is ia to the ! level of accuracy I am working. Note hfactor will be negative ! if the parent e-folding time is longer than the daughter e-folding ! time. ! radio(i)%hfactor=(qval_nu_omit/(real(ia,nprecision) & & *(t_efold-t_efold_parent)) )*rw7_parent radio(i)%ia=ia radio(i)%ia_parent=ia_parent radio(i)%iz=iz radio(i)%iz_parent=iz_parent radio(i)%qval_nu_omit=qval_nu_omit radio(i)%t_efold=t_efold radio(i)%t_efold_parent=t_efold_parent write(3,*) radio(i) ! end do do410 210 continue close(unit=3) close(unit=2) close(unit=1) ! end program rde ! !23456789a123456789b123456789c123456789d123456789e123456789f123456789g12 ! function decay(t) use numerical use physical use astronomical implicit none ! real (kind=nprecision) :: decay real (kind=nprecision), parameter :: ccco=6.70d9*xmass_sun ! Convert to per solar masses. real (kind=nprecision), parameter :: ccni=3.94d10*xmass_sun ! Convert to per solar masses. real (kind=nprecision), parameter :: fpe=3.20d-2 real (kind=nprecision), parameter :: fph=.968d0 ! real (kind=nprecision), parameter :: ftot=fpe+fph ! In principle 1 exactly. real (kind=nprecision), parameter :: ftot=1.d0 real (kind=nprecision), parameter :: gg=0.184641 real (kind=nprecision), parameter :: teco=111.48 real (kind=nprecision), parameter :: teni=8.767 ! real (kind=nprecision) :: decay real (kind=nprecision) :: t ! print*,'ccco,ccni' print*,ccco,ccni ! ! decay=6.31d43*exp(-t/8.8d0) + 1.43d43*exp(-t/111.d0) ! This is wrong in H2006. ! decay=ccni*( exp(-t/teni) + gg*(exp(-t/teco)-exp(-t/teni))*ftot ) ! end function decay