c-------------------------------------------------------------------------------
c     MAGMOD Simulator for NEOMOD
c 
c     Based on the parameters listed in the parameter file (neomod.par), 
c     the code generates absolute magnitudes for 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     Compile: ifort -o magmod_simulator magmod_simulator.f or gfortran -o magmod_simulator magmod_simulator.f  
c     Run: magmod_simulator < magmod.par > output_file.txt 
c     See magmod.par for input parameters
c------------------------------------------------------------------------------- 

      implicit NONE
      
      character*100 filename,modelname
      integer iseed,j
      integer nh,nseg,npoint,nneo,counter
      real*8 nneo_float,hmag1,hmag2,maxmod,ran3,hmag,wmod
      real*8 minh,maxh,href,nref,hpoint(100),minh2,maxh2,gamma(100)
      real*8 nlogp(100),y2(100),nlog0,yp1,ypn,nreal,nmag1,nmag2
      real*8 dh,hmagleft,hmagright

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(*,*)'Exiting...'
         stop
      end if

c     Read the header of the model file
      open(1,file=filename,status='old')
      read(1,'(a)')modelname
      read(1,*)nh,minh,maxh     ! absolute magnitude binning: # of H bins, Hmin,Hmax 
      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 >= than ',minh2
         write(*,*)'Exiting...'
         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 <= than ',maxh2
         write(*,*)'Exiting...'
         stop
      end if
      read(1,*)(gamma(j),j=1,nseg) ! segment slopes (-1.0 indicates N(<15)=50)
      close(1)

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

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.ge.0.d0) then
         nneo=nneo_float
      else
         nneo=nreal
      end if

c     determine maxmod (differential!)
      dh=0.01d0
      hmagleft=hmag2-dh
      hmagright=hmag2
      call splint(hpoint,nlogp,y2,npoint,hmagleft,nmag1)
      call splint(hpoint,nlogp,y2,npoint,hmagright,nmag2)
      maxmod=10.d0**nmag2-10.d0**nmag1
      
c     Generate NEOs
      counter=0
 300  continue
      hmag=hmag1+(hmag2-hmag1)*ran3(iseed)
      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)
      wmod=10.d0**nmag2-10.d0**nmag1
      if(ran3(iseed)*maxmod.lt.wmod) then   ! rejection method
         write(*,8848)hmag
 8848    format(f7.3)
         counter=counter+1
      end if
      if(counter.ge.nneo) stop  ! cycle until we generate the desired number of NEOs 
      goto 300

      end

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

      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