program mssm_flux
      implicit none
      real*8 egam, psi0, theta0, gaflux, dsgafluxsph, mx, dsmwimp
      integer ierr, istat
      character*50 refhalolab, filename, slhafile

      call dsinit

      slhafile = 'stau_coannihilation.slha2'
      call dsgive_model_SLHA(slhafile,0)
      call dsmodelsetup(ierr,istat)
      mx=dsmwimp()
      write(*,*) 'mass = ', mx
c      call dsibset('off')
c      call dsib2set('standard')

      refhalolab = 'mwnfwdef'  ! 'default' NFW profile
      psi0       = 0.0d0       ! direction: Galactic Center
      theta0     = 0.1d0       ! opening angle

      filename='data/flux_MSSM.dat'
      open(unit=47,file=filename,status='unknown',form='formatted')
      write(47,'(A)') '   E [GeV] |  flux [ph cm^-2 s^-1 GeV^-1]'

      egam=mx/1d2               ! start energy
 10   gaflux=dsgafluxsph(egam,1,1.d0,refhalolab,psi0,theta0,istat)
         write(47,*) egam, gaflux
         egam = egam*1.05       ! increase energy by 5% per bin
      if (egam.lt.mx) goto 10

      close(47)
      end