**********************************************************************
*** wimpyields. This program calculates various stable particle yields
*** (positrons, neutrinos, gamma-rays,...) resulting from final states
*** of DM annihilation or decay. It is thus an example of how to use 
*** DarkSUSY without relying on a specific SUSY model (instead linking 
*** to the generic WIMP model).
*** Author: Joakim Edsjo, edsjo@fysik.su.se
*** Date: October 25, 2017
*** mod: May 27, 2020 -- added example for how to compare different
***                      versions of tabulated yields [TB]
**********************************************************************

      program wimpyields
      implicit none
      include 'dsio.h'

      real*8 dsanyield_sim
      integer annpdg,yieldpdg,diff,ndec,ntot
      real*8 yield,mwimp,decmin,e
      character*1 hel
      integer i,istat
      character*80 filename


      call dsinit
      prtlevel = 2 ! This increases the level of output,
                   ! e.g. information about which tables are loaded

c...Setup, see header of src/an_yield/dsanyield_sim.f for details
      annpdg=3                  ! PDG code of annihilation products
      yieldpdg=22               ! PDG code of yield of interest
      mwimp=15d0                ! WIMP mass (GeV)
      hel='0'                   ! helicity state, '0'=unpolarized
      diff=1                    ! differential yields
      decmin=-4.d0              ! lowest energy, 10^-decmin * mwimp
      ndec=30                   ! number of bins per decade
                                ! total number of bins = ndec*(-decmin)
      filename='data/yield_DS_15GeV.dat'
c      filename='yield_Amoroso.dat'
c      filename='yield_Plehn.dat'

c...Change default yield tables
      call dsanyield_set('yieldtables','default')
c      call dsanyield_set('yieldtables','Amoroso') ! tables from S. Amoroso et al.
                                                   ! JCAP05(2019)007
c      call dsanyield_set('yieldtables','Plehn')   ! tables from T. Plehn et al. for
                                                   ! sub-GeV DM, arXiv: 1911.11147 [hep-ph]

c... this is an example of how to change default settings for contributed yield tables
c      call dsanyield_set('Plehn_binning','fine')  ! default: 'coarse'
c      call dsanyield_set('Plehn_model','SM') ! default: 'B-L'

c... the default is to return zero yields outside the range of contributed tables
c... this behaviour can be changed by a call to dsanyield_set('fallback',VALUE).
c... For VALUE='mass', e.g., the standard DS yields are returned outside the
c... tabulated mass range
c      call dsanyield_set('fallback','mass')


c...Preliminary calculations
      ntot=ndec*int(-decmin)

c...Open file and perform calculation      
      open(unit=47,file=filename,status='unknown',form='formatted')
      if (diff.eq.1) then
         write(47,'(A)') '#  i        E [GeV]   yield [GeV^-1 ann^-1]'
      else
         write(47,'(A)') '#  i        E [GeV]   yield [ann^-1]'
      endif
      
      do i=1,ntot
         e=mwimp*10**(decmin-(dble(i)-0.5d0)/abs(ntot)*decmin)
         yield=dsanyield_sim(mwimp,e,annpdg,hel,yieldpdg,diff,istat)
         write(47,100) i,e,yield
 100     format(I4,1x,E14.6,1x,E14.6)
      enddo

      end