C   mirrorHepPh0308254 - make a simple calculation for testing formula (5)
C   of http://de.arXiv.org/abs/hep-ph/0308254
C
C   Copyright (C) 2004 Boud Roukema
C
C   This program is free software; you can redistribute it and/or modify
C   it under the terms of the GNU General Public License as published by
C   the Free Software Foundation; either version 2, or (at your option)
C   any later version.
C
C   This program is distributed in the hope that it will be useful,
C   but WITHOUT ANY WARRANTY; without even the implied warranty of
C   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C   GNU General Public License for more details.
C
C   You should have received a copy of the GNU General Public License
C   along with this program; if not, write to the Free Software Foundation,
C   Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.  

C     To compile and run and sort the results by velocity:
C     g77 mirrorvelocity.f ; ./a.out | sort -n +3

      program  mirror
      
      parameter(mxmirror=4)
      parameter(mxnorm= 2)
      real*4   xmmirror(mxmirror)
      real*4   xmnorm(mxnorm)
      character*2  chmirr(mxmirror),chnorm(mxnorm)

      parameter(mxthresh=3)
      real*4   thresh(mxthresh)

c     H, He, O, Fe
      data  chmirr /'H', 'He', 'O', 'Fe'/
      data  xmmirror / 1.0, 4.0, 16.0, 55.845/ ! atomic units
c     Na,Ge,I,
      data chnorm /'Na','Ge'/ !,'I'/
      data  xmnorm   / 23.0, 72.64/ !, 127.0 / ! atomic units

      data thresh / 10.0, 
     _     4.0,                ! Edeleweiss  astro-ph/0206271
     _     1.0/  ! DAMA/NaI


c     units = sqrt ( 1keV  / (0.938 GeV/c^2) ) 
      units = sqrt ( 1.0 / (0.938e6 / (3e5*3e5) ) ) ! km/s
      print*,'units=',units
      units = sqrt ( 1.783e-33 /1.660539e-27 ) * 3e5

      print*,'units=',units
      
      do imirr=1,mxmirror
         do inorm=1,mxnorm
            do ithresh=1,mxthresh
               xx= xmnorm(inorm) + xmmirror(imirr) 
               vmin = sqrt( ( xx*xx * thresh(ithresh) )
     _              / ( 2.0*xmnorm(inorm)
     _              * xmmirror(imirr) *xmmirror(imirr) ) )
     _              * units ! convert to km/s
               write(6,'(3f8.1,a,f8.1)')xmmirror(imirr),xmnorm(inorm),
     _              thresh(ithresh),'   ',vmin
            enddo
         enddo
      enddo

      end
