c-------------------------------------------------------------------------------
c     NEOMOD Simulator
c
c     The code inputs a binned NEO model generated by MultiNest (input_model.dat). 
c     Based on the parameters listed in the parameter file (neomod.par), 
c     it generates a,e,i,H for any number of model NEOs.  
c     
c     The currently available magnitude range is 9<H<33 mag, but note that NEOMOD     
c     was calibrated for 15<H<28 mag, and H<15 and H>28 are extrapolations
c
c     The output distribution is smoothed by interpolation.
c     Columns 6-17 in the input file list relative contribution of sources to each bin 
c
c     Compile: ifort -o neomod_simulator neomod_simulator.f or gfortran -o neomod_simulator neomod_simulator.f  
c     Run: neomod_simulator < neomod.par > output_file.txt 
c     See neomod.par for input parameters
c-------------------------------------------------------------------------------

      implicit NONE

      integer mh,ma,me,mi
      parameter(mh=52,ma=42,me=25,mi=22)  ! needs to be >= than the number of bins in the model file

      integer ms
      parameter(ms=12)                    ! >= than the number of sources (nsour in the model file)
      
      character*100 filename,modelname
      integer ih,ia,ie,ii,nh,na,ne,ni
      real*8 model(mh,ma,me,mi),modelin
      real*8 minh,maxh,mina,maxa,mine,maxe,mini,maxi
      real*8 hmag,a,e,inc,da,de,di,dh
      real*8 alpha(mh,ma,me,mi,ms)
      
      integer j,npoint,nseg,nsour
      real*8 href,nref,gamma(100),hpoint(100),nlogp(100)
      real*8 yp1,ypn,y2(100),nmag1,nmag2,hmag1,hmag2,nlog0

      integer nneo,counter
      real*8 nreal,nneo_float,nmag

      integer iseed
      real*8 ran3,x,maxmod,ran3_new,ran2

      integer ia0,ia1,ie0,ie1,ii0,ii1,ih0,ih1
      real*8 bh,ba,be,bi,wmod
      real*8 xd,yd,zd,x0,x1,y0,y1,z0,z1
      real*8 c0,c1,c00,c10,c01,c11
      real*8 c000,c100,c010,c001,c111,c011,c101,c110
      real*8 v0,v1,vd,wmod0,wmod1

      real*8 minh2,maxh2,hmag_fake,cfac,hmagleft,hmagright,num1,num2
      
c     Read input parameters from neomod.par
      read(*,'(a)')filename ! NEO model file (only for bookkeeping)
      read(*,*)iseed        ! seed for the random number generator 
      read(*,*)nneo_float   ! number of NEOs to be generated, <0 for real number
      read(*,*)hmag1,hmag2  ! desired magnitude range
      if(hmag1.ge.hmag2) then
         write(*,*)'hmag2 in neomod.par must be greater than hmag1'
         write(*,*)'Exciting...'
         stop
      end if
            
c     Zero the orbit model
      do ih=1,mh
         do ia=1,ma
            do ie=1,me
               do ii=1,mi
                  model(ih,ia,ie,ii)=0.d0
               end do
            end do
         end do
      end do
         
c     Read the file with binned PDF of orbits
      open(1,file=filename,status='old')
      read(1,'(a)')modelname
      read(1,*)nh,minh,maxh     ! absolute magnitude binning: # of H bins, Hmin,Hmax 
      dh=(maxh-minh)/float(nh)
      read(1,*)href,nref        ! reference magnitude, reference pop 
      read(1,*)nseg        ! number of spline segments for absolute magnitude distribution  
      npoint=nseg+1        ! number of boundaries between segments
      read(1,*)(hpoint(j),j=1,npoint)  ! segment boundaries
      minh2=hpoint(1)
      maxh2=hpoint(npoint)
      if(hmag1.lt.minh2) then
         write(*,*)
     &'The requested magnitude range is outside of model domain'
         write(*,*)
     &'hmag1 in neomod.par must be equal or greater than ',minh2
         write(*,*)'Exciting...'
         stop
      end if
      if(hmag1.ge.maxh) then
         write(*,*)
     &'The requested magnitude range is outside of model domain'
         write(*,*)'hmag1 in neomod.par must be smaller than ',maxh
         write(*,*)'Exciting...'
         stop
      end if
      if(hmag2.gt.maxh2) then
         write(*,*)
     &'The requested magnitude range is outside of model domain'
         write(*,*)
     &'hmag2 in neomod.par must be equal or smaller than ',maxh2
         write(*,*)'Exciting...'
         stop
      end if
      if(hmag2.lt.minh) then
         write(*,*)
     &'The requested magnitude range is outside of model domain'
         write(*,*)
     &'hmag2 in neomod.par must be greater than ',minh
         write(*,*)'Exciting...'
         stop
      end if 

      read(1,*)(gamma(j),j=1,nseg) ! segment slopes (-1.0 indicates N(<15)=50)
      read(1,*)na,mina,maxa ! semimajor axis binning: # of a bins, amin (au), amax(au)
      read(1,*)ne,mine,maxe ! eccentricity binning: # of e bins, emin, emax
      read(1,*)ni,mini,maxi ! orbital inclination binning: # of i bins, imin (deg), imax (deg)
      read(1,*)nsour            ! number of sources 
      read(1,*) ! ignore header 
      read(1,*) ! ih,ia,ie,ii,modelin, plus "nsour" columns with relative source contributions   
      read(1,*) !
      maxmod=0.d0
 100  continue
      read(1,*,end=200,err=200)ih,ia,ie,ii,modelin, ! binned differential distribution
     &     (alpha(ih,ia,ie,ii,j),j=1,nsour) ! source contributions
      model(ih,ia,ie,ii)=modelin
      if(minh+ih*dh.ge.hmag1.and.minh+(ih-1)*dh.le.hmag2)then  ! set maxmod for hmag1<=bin<=hmag2
!     this is not executed for bins with partial overlap with the desired range 
!     solution: set hmag1,hmag2 as multiples of 0.25 mag 
         if(modelin.gt.maxmod)maxmod=modelin  ! maxmod in the overlap of the model domain with (hmag1,hmag2) 
      end if
      goto 100
 200  continue
      close(1)

c     Bin sizes
      da=(maxa-mina)/float(na)
      de=(maxe-mine)/float(ne)
      di=(maxi-mini)/float(ni)

c     Define cumulative magnitude distribution logN(<H) from splines
      nlog0=log10(nref) 
      nlogp(5)=nlog0+(hpoint(5)-href)*gamma(4) ! Href falls into 4nd segment
      nlogp(4)=nlog0+(hpoint(4)-href)*gamma(4)
      nlogp(3)=nlogp(4)+(hpoint(3)-hpoint(4))*gamma(3)
      if(gamma(2).ge.0.d0) then
         nlogp(2)=nlogp(3)+(hpoint(2)-hpoint(3))*gamma(2)
      else 
         nlogp(2)=log10(50.d0)
         gamma(2)=(nlogp(2)-nlogp(3))/(hpoint(2)-hpoint(3))
      end if
      nlogp(1)=nlogp(2)+(hpoint(1)-hpoint(2))*gamma(1)
      do j=6,npoint
         nlogp(j)=nlogp(j-1)+(hpoint(j)-hpoint(j-1))*gamma(j-1)
      end do
      yp1=1.d30
      ypn=1.d30
      call spline(hpoint,nlogp,npoint,yp1,ypn,y2)

c     Correct maxmod if hmag2>maxh - magnitude range in the extrapolated domain
      if(hmag2.gt.maxh) then
         hmagleft=hmag2-0.5d0*dh
         hmagright=hmag2+0.5d0*dh
         call splint(hpoint,nlogp,y2,npoint,hmagleft,nmag1)
         call splint(hpoint,nlogp,y2,npoint,hmagright,nmag2)
         num1=10.d0**nmag2-10.d0**nmag1
         hmagleft=maxh-dh
         hmagright=maxh         
         call splint(hpoint,nlogp,y2,npoint,hmagleft,nmag1)
         call splint(hpoint,nlogp,y2,npoint,hmagright,nmag2)
         num2=10.d0**nmag2-10.d0**nmag1
         maxmod=maxmod*num1/num2
      end if

c      hmag=17.75d0
c      call splint(hpoint,nlogp,y2,npoint,hmag,nmag)
c      write(*,*)'Number of H<',hmag,10.d0**nmag
c      stop

c     Rescale distribution to the desired number of model NEOs (nneo)
      call splint(hpoint,nlogp,y2,npoint,hmag1,nmag1)
      call splint(hpoint,nlogp,y2,npoint,hmag2,nmag2)
      nreal=10.d0**nmag2-10.d0**nmag1

      if(nneo_float.gt.0.d0) then
         nneo=nneo_float
      else
         nneo=nreal
      end if

c     Generate model NEOs
      counter=0
c      write(*,'(a)')'----------------------------------------'
c      write(*,'(a)')'     H(mag)         e   inc(deg)   a(au)'
c      write(*,'(a)')'----------------------------------------'
 300  continue
      hmag=hmag1+ran3(iseed)*(hmag2-hmag1)
      hmag_fake=hmag ! this is used for the orbital distribution, no info outside model domain 
      if(hmag.le.minh) hmag_fake=minh+1.d-10 ! with safety margin
      if(hmag.ge.maxh) hmag_fake=maxh-1.d-10

      a=mina+ran3(iseed)*(maxa-mina)
      e=mine+ran3(iseed)*(maxe-mine)
      inc=mini+ran3(iseed)*(maxi-mini)
      if(a*(1.d0-e).gt.1.3d0) goto 300
      
      ih=int((hmag_fake-minh)/dh)+1  ! here we use the last available bin for extrapolation 
      ia=int((a-mina)/da)+1
      ie=int((e-mine)/de)+1
      ii=int((inc-mini)/di)+1      

c---------------------------------------------------------------

c     Smooth orbital elements and absolute magnitudes 

c     Power law (linear in log) interpolation for absolute magnitudes
      bh=minh+(ih-0.5d0)*dh
      if(hmag_fake.gt.bh) then
         ih0=ih
         ih1=ih+1
      else
         ih0=ih-1
         ih1=ih
      end if
      v0=minh+(ih0-0.5d0)*dh
      v1=minh+(ih1-0.5d0)*dh
      if(ih0.lt.1)ih0=1    ! border effects
      if(ih1.gt.nh)ih1=nh
            
c     Trilinear interpolation for orbital elements
      ba=mina+(ia-0.5d0)*da
      be=mine+(ie-0.5d0)*de
      bi=mini+(ii-0.5d0)*di
      if(a.gt.ba) then
         ia0=ia
         ia1=ia+1
      else
         ia0=ia-1
         ia1=ia
      end if
      if(e.gt.be) then
         ie0=ie
         ie1=ie+1
      else
         ie0=ie-1
         ie1=ie
      end if
      if(inc.gt.bi) then
         ii0=ii
         ii1=ii+1
      else
         ii0=ii-1
         ii1=ii
      end if
      x0=mina+(ia0-0.5d0)*da
      x1=mina+(ia1-0.5d0)*da
      y0=mine+(ie0-0.5d0)*de
      y1=mine+(ie1-0.5d0)*de
      z0=mini+(ii0-0.5d0)*di
      z1=mini+(ii1-0.5d0)*di
      xd=(a-x0)/(x1-x0)
      yd=(e-y0)/(y1-y0)
      zd=(inc-z0)/(z1-z0)
      if(ia0.lt.1)ia0=1
      if(ia1.gt.na)ia1=na
      if(ie0.lt.1)ie0=1
      if(ie1.gt.ne)ie1=ne
      if(ii0.lt.1)ii0=1
      if(ii1.gt.ni)ii1=ni

c     Step 1: trilinear for ih0
      c000=model(ih0,ia0,ie0,ii0)
      c100=model(ih0,ia1,ie0,ii0)
      c010=model(ih0,ia0,ie1,ii0)
      c001=model(ih0,ia0,ie0,ii1)
      c111=model(ih0,ia1,ie1,ii1)
      c011=model(ih0,ia0,ie1,ii1)
      c101=model(ih0,ia1,ie0,ii1)
      c110=model(ih0,ia1,ie1,ii0)
      c00=c000*(1.d0-xd)+c100*xd
      c01=c001*(1.d0-xd)+c101*xd
      c10=c010*(1.d0-xd)+c110*xd
      c11=c011*(1.d0-xd)+c111*xd
      c0=c00*(1.d0-yd)+c10*yd
      c1=c01*(1.d0-yd)+c11*yd
      wmod0=c0*(1.d0-zd)+c1*zd

c     Step 2: trilinear for ih1
      c000=model(ih1,ia0,ie0,ii0)
      c100=model(ih1,ia1,ie0,ii0)
      c010=model(ih1,ia0,ie1,ii0)
      c001=model(ih1,ia0,ie0,ii1)
      c111=model(ih1,ia1,ie1,ii1)
      c011=model(ih1,ia0,ie1,ii1)
      c101=model(ih1,ia1,ie0,ii1)
      c110=model(ih1,ia1,ie1,ii0)
      c00=c000*(1.d0-xd)+c100*xd
      c01=c001*(1.d0-xd)+c101*xd
      c10=c010*(1.d0-xd)+c110*xd
      c11=c011*(1.d0-xd)+c111*xd
      c0=c00*(1.d0-yd)+c10*yd
      c1=c01*(1.d0-yd)+c11*yd
      wmod1=c0*(1.d0-zd)+c1*zd

c     Step 3: power law for magnitudes
      vd=(hmag_fake-v0)/(v1-v0)
      wmod=log10(wmod0+1.d-30)*(1.d0-vd)+log10(wmod1+1.d-30)*vd
      wmod=10.d0**wmod

c     Step 4: correct hmod if outside the (hmin,hmax) interval
      cfac=1.d0
      if(hmag.lt.minh) then
         hmagleft=hmag-0.5d0*dh
         hmagright=hmag+0.5d0*dh         
         call splint(hpoint,nlogp,y2,npoint,hmagleft,nmag1)
         call splint(hpoint,nlogp,y2,npoint,hmagright,nmag2)
         num1=10.d0**nmag2-10.d0**nmag1
         hmagleft=minh
         hmagright=minh+dh         
         call splint(hpoint,nlogp,y2,npoint,hmagleft,nmag1)
         call splint(hpoint,nlogp,y2,npoint,hmagright,nmag2)
         num2=10.d0**nmag2-10.d0**nmag1
         cfac=num1/num2
      end if
      if(hmag.gt.maxh) then
         hmagleft=hmag-0.5d0*dh
         hmagright=hmag+0.5d0*dh         
         call splint(hpoint,nlogp,y2,npoint,hmagleft,nmag1)
         call splint(hpoint,nlogp,y2,npoint,hmagright,nmag2)
         num1=10.d0**nmag2-10.d0**nmag1
         hmagleft=maxh-dh
         hmagright=maxh         
         call splint(hpoint,nlogp,y2,npoint,hmagleft,nmag1)
         call splint(hpoint,nlogp,y2,npoint,hmagright,nmag2)
         num2=10.d0**nmag2-10.d0**nmag1
         cfac=num1/num2
      end if
      wmod=wmod*cfac

      if(wmod.gt.maxmod) then
         write(*,*)'Error: wmod>maxmod, exiting...'
         stop
      end if

c----------------------------------------------------------------
      
c     Output model NEOs
      if(ran3(iseed)*maxmod.lt.wmod) then
         write(*,8848)hmag,a,e,inc
 8848    format(f7.3,f6.3,f7.4,f6.2)
c         write(*,8848)hmag,a,e,inc,(alpha(ih,ia,ie,ii,j),j=1,nsour)   
c 8848    format(4(1x,f9.5),12(1x,f6.3))
        counter=counter+1
      end if
      
      if(counter.ge.nneo) stop  ! cycle until we generate the desired number of NEOs 
      
      goto 300

      end
      
      SUBROUTINE spline(x,y,n,yp1,ypn,y2)
      INTEGER n,NMAX
      REAL*8 yp1,ypn,x(n),y(n),y2(n)
      PARAMETER (NMAX=500)
      INTEGER i,k
      REAL*8 p,qn,sig,un,u(NMAX)
      if (yp1.gt.0.99d30) then
         y2(1)=0.d0
         u(1)=0.d0
      else
         y2(1)=-0.5d0
         u(1)=(3.d0/(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
      endif
      do i=2,n-1
         sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
         p=sig*y2(i-1)+2.d0
         y2(i)=(sig-1.d0)/p
         u(i)=(6.d0*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1))
     &        /(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*u(i-1))/p
      enddo
      if (ypn.gt.0.99d30) then
         qn=0.d0
         un=0.d0
      else
         qn=0.5d0
         un=(3.d0/(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
      endif
      y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.d0)
      do k=n-1,1,-1
         y2(k)=y2(k)*y2(k+1)+u(k)
      enddo 
      return
      END

      SUBROUTINE splint(xa,ya,y2a,n,x,y)
      INTEGER n
      REAL*8 x,y,xa(n),y2a(n),ya(n)
      INTEGER k,khi,klo
      REAL*8 a,b,h
      klo=1
      khi=n
 1    if(khi-klo.gt.1) then
         k=(khi+klo)/2
         if(xa(k).gt.x)then
            khi=k
         else
            klo=k
         endif
         goto 1
      endif
      h=xa(khi)-xa(klo)
c      if (h.eq.0.) pause ’bad xa input in splint’
      a=(xa(khi)-x)/h
      b=(x-xa(klo))/h
      y=a*ya(klo)+b*ya(khi)+
     &     ((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h**2)/6.d0
      return
      END

      function ran3(idum)
      implicit real*8(a-h,o-z)
      parameter (mbig=1000000000,mseed=161803398,mz=0,fac=1.d-9)
      dimension ma(55)
      data iff /0/
      save
c --------------------------------------------------------------------
      if((idum.ge.0).and.(iff.ne.0)) goto 20
        iff=1
        mj=mseed-iabs(idum)
        mj=mod(mj,mbig)
        ma(55)=mj
        mk=1
        do 11 i=1,54
          ii=mod(21*i,55)
          ma(ii)=mk
          mk=mj-mk
          if(mk.lt.mz)mk=mk+mbig
          mj=ma(ii)
11      continue
        do 13 k=1,4
          do 12 i=1,55
            ma(i)=ma(i)-ma(1+mod(i+30,55))
            if(ma(i).lt.mz)ma(i)=ma(i)+mbig
12        continue
13      continue
        inext=0
        inextp=31
        idum=1
20    continue
      inext=inext+1
      if(inext.eq.56)inext=1
      inextp=inextp+1
      if(inextp.eq.56)inextp=1
      mj=ma(inext)-ma(inextp)
      if(mj.lt.mz)mj=mj+mbig
      ma(inext)=mj
      ran3=mj*fac
      return
      end

      FUNCTION ran3_new(idum)
      INTEGER idum
      INTEGER MBIG,MSEED,MZ
C     REAL MBIG,MSEED,MZ
      REAL*8 ran3,FAC
      PARAMETER (MBIG=1000000000,MSEED=161803398,MZ=0,FAC=1.d0/MBIG)
C     PARAMETER (MBIG=4000000.,MSEED=1618033.,MZ=0.,FAC=1./MBIG)
      INTEGER i,iff,ii,inext,inextp,k
      INTEGER mj,mk,ma(55)
C     REAL mj,mk,ma(55)
      SAVE iff,inext,inextp,ma
      DATA iff /0/
      if(idum.lt.0.or.iff.eq.0)then
        iff=1
        mj=MSEED-iabs(idum)
        mj=mod(mj,MBIG)
        ma(55)=mj
        mk=1
        do 11 i=1,54
          ii=mod(21*i,55)
          ma(ii)=mk
          mk=mj-mk
          if(mk.lt.MZ)mk=mk+MBIG
          mj=ma(ii)
11      continue
        do 13 k=1,4
          do 12 i=1,55
            ma(i)=ma(i)-ma(1+mod(i+30,55))
            if(ma(i).lt.MZ)ma(i)=ma(i)+MBIG
12        continue
13      continue
        inext=0
        inextp=31
        idum=1
      endif
      inext=inext+1
      if(inext.eq.56)inext=1
      inextp=inextp+1
      if(inextp.eq.56)inextp=1
      mj=ma(inext)-ma(inextp)
      if(mj.lt.MZ)mj=mj+MBIG
      ma(inext)=mj
      ran3_new=mj*FAC
      return
      END

      FUNCTION ran2(idum)
      INTEGER idum,IM1,IM2,IMM1,IA1,IA2,IQ1,IQ2,IR1,IR2,NTAB,NDIV
      REAL*8 ran2,AM,EPS,RNMX
      PARAMETER (IM1=2147483563,IM2=2147483399,AM=1.d0/IM1,IMM1=IM1-1,
     *IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211,IR2=3791,
     *NTAB=32,NDIV=1+IMM1/NTAB,EPS=1.2d-7,RNMX=1.d0-EPS)
      INTEGER idum2,j,k,iv(NTAB),iy
      SAVE iv,iy,idum2
      DATA idum2/123456789/, iv/NTAB*0/, iy/0/
      if (idum.le.0) then
        idum=max(-idum,1)
        idum2=idum
        do 11 j=NTAB+8,1,-1
          k=idum/IQ1
          idum=IA1*(idum-k*IQ1)-k*IR1
          if (idum.lt.0) idum=idum+IM1
          if (j.le.NTAB) iv(j)=idum
11      continue
        iy=iv(1)
      endif
      k=idum/IQ1
      idum=IA1*(idum-k*IQ1)-k*IR1
      if (idum.lt.0) idum=idum+IM1
      k=idum2/IQ2
      idum2=IA2*(idum2-k*IQ2)-k*IR2
      if (idum2.lt.0) idum2=idum2+IM2
      j=1+iy/NDIV
      iy=iv(j)-idum2
      iv(j)=idum
      if(iy.lt.1)iy=iy+IMM1
      ran2=min(AM*iy,RNMX)
      return
      END