**********************************************************************
*** neutrinoyields. This program calculates yields in neutrino
*** telescopes for a generic WIMP.
*** Author: Joakim Edsjo, edsjo@fysik.su.se
*** Date: June 24, 2021
**********************************************************************

      program neutrinoyields
      implicit none
      include 'dsio.h'


      real*8 mwimp,svann,sigsi,sigsd
      logical selfconj
      real*8 eth,thmax,rho0in,rateea,ratesu
      integer kind,rtype,ptype,pdgann,istat

c...Initialize DarkSUSY and set default parameters
      call dsinit
c      call dssenu_set('tabmed')  ! default capture calculation
c      call dssem_sunset('default') ! default solar model

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

c...Set up astrophysics (vel. distribution is currently the default)      
      rho0in=0.3d0                 ! local halo density in GeV/cm^3
      
c...Set up particle physics model
      mwimp=500
      selfconj=.true.           ! self-conjugated DM (it own anti-particle)
      svann=3d-26               ! ann cross sect, cm^3/s
      pdgann=24                 ! annihilation to W+ W- (pdg code 24)
      sigsi=1.0d-7              ! SI scattering cross section, pb
      sigsd=1.0d-5              ! SD scattering cross section, pb

      call dsgivemodel_generic_wimp(mwimp,selfconj,svann,pdgann,sigsi)

c...Note, sigsd is not part of a typcial generic WIMP model, give it separately
      call dsgivemodel_generic_wimp_opt('sigsdp',sigsd)
      call dsgivemodel_generic_wimp_opt('sigsdn',sigsd)

c...Set up neutrino-telescope parameters
      eth=1.0d0      ! energy threshold (of neutrino/muon), GeV
      thmax=30.0d0   ! the maximum half-aperture angle, degrees
      kind=1         ! 1=integrated
                     ! 2=differential
                     ! 3=mixed (diff in E, integrated in theta)
      rtype=3        ! 1=neutrino flux
                     ! 2=neutrino to muon conversion rate
                     ! 3=muon flux
      ptype=3        ! 1=particles only
                     ! 2=anti-particles only
                     ! 3=sum of particle and anti-particle rates

c... NB: Assume here that the active DM particle makes up 100% of the local
c...     DM density. Otherwise, rescale rho0in by a corresponding fraction
      call dssenu_rates(eth,thmax,1,rtype,ptype,rho0in,rateea,ratesu,
     &  istat)

      write(*,*) '  Flux from the Earth = ',rateea, ' km^-2 yr^-1'
      write(*,*) '  Flux from the Sun =   ',ratesu, ' km^-2 yr^-1'


      end