!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.
!
!
!  <a href="http://ie.lbl.gov/toi/"> LBNL Isotopes 
!      Project - LUNDS Universitet
!      WWW Table of Radioactive Isotopes </a>
!         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
      implicit none
!
      integer, parameter :: nelevel=200
      integer, parameter :: nxray=200
!
      real (kind=nprecision) :: con_mev_erg
      real (kind=nprecision) :: con_thalf_tefold
      integer :: i,j,k,l,m,n 
      integer :: idata
      integer :: ielevel
      real (kind=nprecision) :: small 
      real (kind=nprecision) :: xmass_np_diff 
!
      character :: aanuclide*10           ! Nuclide name. 
      real (kind=nprecision) :: b_coef    ! Energy generation B coefficient 
!                                         !   (erg/s/g).
      real (kind=nprecision) :: b_coef_lepton  ! For betas and atomic el.
      real (kind=nprecision) :: b_coef_photon ! For gammas and X-rays. 
      real (kind=nprecision) :: b_coef_sun  ! Energy generation B coefficient 
!                                           !   (erg/s/mass_sun).
      real (kind=nprecision) :: betamax   ! Maximum beta particle energy (MeV).
      real (kind=nprecision) :: betafrac  ! Fraction of time there is a 
!                                         !   beta particle. 
      real (kind=nprecision) :: beta_ke   ! Average beta particle energy 
!                                         !  per beta (MeV).
      real (kind=nprecision) :: c_coef    ! Energy generation C 
!                                         !  coefficient (erg/s/g).
      real (kind=nprecision) :: c_coef_lepton  ! For betas and atomic el. 
      real (kind=nprecision) :: c_coef_photon ! For gammas and X-rays. 
      real (kind=nprecision) :: c_coef_sun  ! Energy generation C 
!                                           !   coefficient (erg/s/mass_sun).
      real (kind=nprecision) :: elevel(10,nelevel) ! Energy level 
!                                                  !  and beta/EC data. 
      real (kind=nprecision) :: fatomic   ! Fraction of energy in atomic 
!                                         !   electron KE.
      real (kind=nprecision) :: fbeta     ! Fraction of energy in positron KE.
      real (kind=nprecision) :: fgamma    ! Fraction of energy in gamma-rays.
      real (kind=nprecision) :: flepton   ! Fraction of energy in leptons. 
      real (kind=nprecision) :: fphoton   ! Fraction of energy in photons. 
      real (kind=nprecision) :: fxray     ! Fraction of energy in X-rays. 
      integer :: ibeta                    ! The sign of the beta particle. 
      integer :: iz                       ! Atomic number Z.
      integer :: ia                       ! Atomic mass number A.
      integer :: ia_parent                ! Atomic mass number A of parent.
      integer :: ixray                    ! Number of X-ray emissions. 
      real (kind=nprecision) :: qval      ! Total disintegration energy (MeV).
      real (kind=nprecision) :: qval_nu_omit  ! The energy that doesn't escape 
!                                             !   as neutrinos (MeV).
      real (kind=nprecision) :: qvalerr   ! Error in qval.
      real (kind=nprecision) :: qval_atomic  ! The energy in atomic electrons. 
      real (kind=nprecision) :: qval_beta    ! The energy in betas. 
      real (kind=nprecision) :: qval_gamma   ! The energy in gammas. 
      real (kind=nprecision) :: qval_gamma_direct ! The energy in gamma lines. 
      real (kind=nprecision) :: qval_xray    ! The energy in gammas. 
      real (kind=nprecision) :: tefold    ! e-folding time (days). 
      real (kind=nprecision) :: tefold_parent    ! e-folding time (days) 
!                                                !   of parent. 
      real (kind=nprecision) :: thalf     ! Half-life (days).
      real (kind=nprecision) :: thalferr  ! Half-life error (days).
      real (kind=nprecision) :: thalf_parent     ! Half-life (days) of parent.
      real (kind=nprecision) :: xray(2,nxray)    ! X-ray data. 
      real (kind=nprecision) :: zzfiction ! Just a place-holder variable. 
!
!-----------------------------------------------------------------------
!
      namelist/param/aanuclide                                              &
     & ,b_coef                                                              &
     & ,b_coef_lepton                                                       &
     & ,b_coef_photon                                                       &
     & ,b_coef_sun                                                          &
     & ,betamax                                                             &
     & ,betafrac                                                            &
     & ,beta_ke                                                             &
     & ,c_coef                                                              &
     & ,c_coef_lepton                                                       &
     & ,c_coef_photon                                                       &
     & ,c_coef_sun                                                          &
     & ,fatomic                                                             &
     & ,fbeta                                                               &
     & ,fgamma                                                              &
     & ,flepton                                                             &
     & ,fphoton                                                             &
     & ,fxray                                                               &
     & ,ibeta                                                               &
     & ,ia                                                                  &
     & ,ia_parent                                                           &
     & ,iz                                                                  &
     & ,iz_parent                                                           &
     & ,qval                                                                &
     & ,qval_nu_omit                                                        &
     & ,qvalerr                                                             &
     & ,qval_atomic                                                         &
     & ,qval_beta                                                           &
     & ,qval_gamma                                                          &
     & ,qval_xray                                                           &
     & ,tefold                                                              &
     & ,tefold_parent                                                       &
     & ,thalf                                                               &
     & ,thalferr                                                            &
     & ,thalf_parent                                                        &
     & ,zzfiction
!
      namelist/param2/elevel,xray
!
!-----------------------------------------------------------------------
      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_thalf_tefold=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_thalf_tefold'
      write(*,*) con_mev_erg,con_thalf_tefold
!
!-----------------------------------------------------------------------
      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)
        i=0
        do410 : do
          elevel(:,:)=0.d0
          read(1,param,end=210)
          read(1,param2,end=210)
          i=i+1
          write(*,*)
          write(*,*) 'Nuclide ',i
!
          tefold=thalf*con_thalf_tefold
!
!-----------------------------------------------------------------------
!          write(*,*) 'Before counting the energy levels and xray ',
!          write(*,*) 'emissions (IB plus lines).'
!-----------------------------------------------------------------------
!
          j=1
          do420 : do j=1,nelevel
            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 fashion as above by counting up the energy
!   That must come from decays from the excited child. 
!
!   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.
            else
              beta_ke=0.d0
          end if
!
          c_coef=(qval_nu_omit*con_mev_erg)                                &
     &     /(amu*real(ia,nprecision)*tefold*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
              tefold_parent=thalf_parent*con_thalf_tefold
              b_coef=(qval_nu_omit*con_mev_erg)                            &
     &         /(amu*real(ia,nprecision)*(tefold-tefold_parent)*daysec)
            else
              tefold_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(*,param)
!
        end do do410
  210   continue
      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
