c--------------------------------------------------------------------
c     Given the orbital elements and absolute magnitude of an object 
c     on input the program finds the respective bin in the model file
c     and outputs probabilities that the object was produced from 
c     NEOMOD sources, individually for each source
c
c     Compile: ifort -o give_me_prob give_me_prob.f or gfortran -o give_me_prob give_me_prob.f  
c     Run: give_me_prob 
c--------------------------------------------------------------------

      implicit NONE

      integer ia,ie,ii,iH,na,ne,ni,nH,nsource
      integer ia2,ie2,ii2,iH2,j
      real*8 prob(100),pdf 
      real*8 a,e,inc,hmag
      real*8 amin,amax,emin,emax,imin,imax,Hmin,Hmax
      real*8 da,de,di,dH
      character*100 modelfilename
      character*10 sourcename(100)

c     names of sources, must correspond to the model file header
      sourcename(1)='nu6'
      sourcename(2)='3:1'
      sourcename(3)='5:2'
      sourcename(4)='7:3'
      sourcename(5)='8:3'
      sourcename(6)='9:4'
      sourcename(7)='11:5'
      sourcename(8)='2:1'
      sourcename(9)='inner'
      sourcename(10)='Hungarias'
      sourcename(11)='Phocaeas'
      sourcename(12)='JFCs'

      write(*,*)'a (au), e, i (deg), H (mag)'
      read(*,*)a,e,inc,hmag

      open(1,file='neomod.par',status='old')
      read(1,'(a)') modelfilename
      close(1)

      open(1,file=modelfilename,status='old')
      read(1,*)
      read(1,*)nH,Hmin,Hmax
      read(1,*)
      read(1,*)
      read(1,*)
      read(1,*)
      read(1,*)na,amin,amax
      read(1,*)ne,emin,emax
      read(1,*)ni,imin,imax
      read(1,*)nsource
      da=(amax-amin)/dfloat(na)
      de=(emax-emin)/dfloat(ne)
      di=(imax-imin)/dfloat(ni)
      dH=(Hmax-Hmin)/dfloat(nH)
      read(1,*)
      read(1,*)
      read(1,*)
      if(a.lt.amin.or.a.gt.amax.or.
     &     e.lt.emin.or.e.gt.emax.or.
     &     inc.lt.imin.or.inc.gt.imax.or.
     &     hmag.lt.Hmin.or.hmag.gt.Hmax) then
         goto 200
      end if
      ia=int((a-amin)/da)+1
      ie=int((e-emin)/de)+1
      ii=int((inc-imin)/di)+1
      iH=int((hmag-hmin)/dH)+1
      if(ia.gt.na)ia=na
      if(ie.gt.ne)ie=ne
      if(ii.gt.ni)ii=ni
      if(iH.gt.nH)iH=nH
c     write(*,*)ia,ie,ii,iH
 100  continue
      read(1,*,end=200)iH2,ia2,ie2,ii2,pdf,(prob(j),j=1,nsource)
      if(ia2.eq.ia.and.ie2.eq.ie.and.ii2.eq.ii.and.iH2.eq.iH)then
         goto 300
      end if         
      goto 100
 200  continue
      write(*,*)'Orbit or absolute magnitude out of range'
      write(*,*)'Exiting...'
      stop
 300  continue
      write(*,'(a)')'-----------------'
      write(*,'(a)')'  Probabilities  '
      write(*,'(a)')'-----------------'
      do j=1,nsource
         write(*,8848)sourcename(j),prob(j)
 8848    format(a10,f7.3)
      end do
      write(*,'(a)')'-----------------'
c      write(*,*)'For source names, see the header of',modelfilename
c      write(*,*)pdf

      close(1)
      end