program Jfactors
      implicit none
      include 'dsmpconst.h'                         ! to have access to pi
      integer i, npoints
      real*8 cospsi0, psi0, theta, annwidth         ! direction of l.o.s. integration
      real*8 rhos,rs,distance,rinn,rmax,alpha,rho0  ! halo parameters
      real*8 jpsi0, jpsi, dsjfactor
      character*80 halolabel, outfile

      call dsinit
      outfile='data/Jfactors_Einasto_017_20kpc_core_1kpc.dat'
      open(unit=10,file=outfile)
      write(10,*) 'Theta [deg]          | J [GeV**2 cm**-5]  '
     
      distance=8.5d0 ! kpc, distance of halo to earth
      rinn=1.d0      ! kpc, inner cutoff (min 1.0d-3 or so)
      rmax=50.d0     ! kpc,  radial size of DM halo
      rho0=0.4d0     ! GeV/cm^3, local DM density
      psi0=0.0d0     ! Direction to GC [in degree]
      alpha=0.17     ! Einasto parameter alpha
      rs=20.0d0      ! kpc, scale radius
      rhos = rho0*exp(2./alpha*((distance/rs)**alpha-1.))
      halolabel='dummyein'
      call sethalo_ein(halolabel,rhos,rs,alpha,distance,rinn,rmax)
     
      npoints=200    ! # points + 1
      do i=0,npoints
        theta =15.d0*i/(1.0d0*npoints) ! (Half) opening angle of observation cone [deg]
        annwidth = 0.1d0 ! annwidth/2 in both directions from given theta
        jpsi0=0.0d0
        if (theta.gt.annwidth/2.) jpsi0=dsjfactor(halolabel,psi0,(theta-annwidth/2)*pi/180.)
        jpsi=dsjfactor(halolabel,psi0,(theta+annwidth/2)*pi/180.) - jpsi0
        write(10,*) theta, jpsi*kpc*1.0d21
      enddo

      close(10)
      end ! end of main program



********************************************************************
***   profile setting routine from DMhalo_predef.                ***
***   See there for detailed comments!                           ***   
********************************************************************
      subroutine sethalo_ein(tag,rhos,rs,alpha,distance,rinn,rmax)
      implicit none
      character*12 tag
      real*8 rhos,rs,alpha,distance,rinn,rmax
      logical match
      integer nin
      real*8 rein(6)
      character*10 chin(6)
      call dslabcheck(12,tag,3,'ein',match)
      if(.not.match) then
        write(*,*) 'DS: call to sethalo_ein with tag = ',tag
        write(*,*) 'DS: not containg the string - ein'
        write(*,*) 'DS: program stopped'
        stop
      endif  
      chin(1)='objdist'
      rein(1)=distance
      chin(2)='radintr'
      rein(2)=rinn
      chin(3)='radouttr'
      rein(3)=rmax
      chin(4)='rhosein'
      rein(4)=rhos
      chin(5)='rsein'
      rein(5)=rs
      chin(6)='alphaein'
      rein(6)=alpha
      nin=6
      call dsdmdset_halomodel(tag,nin,chin,rein)
      return
      end