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 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(10),hpoint(10),nlogp(10) real*8 yp1,ypn,y2(10),nmag1,nmag2,hmag1,hmag2,nlog0 integer nneo,counter real*8 nreal,nneo_float,nmag integer iseed real*8 ran3,x,maxmod 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 c Read input parameters from param.in 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 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 if(hmag1.lt.minh) then write(*,*) &'The requested magnitude range is outside of model domain' write(*,*)'hmag1 in neomod.in must be greater than ',minh write(*,*)'Exciting...' stop end if if(hmag2.gt.maxh) then write(*,*) &'The requested magnitude range is outside of model domain' write(*,*)'hmag2 in neomod.in must be smaller than ',maxh write(*,*)'Exciting...' stop end if 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 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 PDF & (alpha(ih,ia,ie,ii,j),j=1,nsour) ! source contributions model(ih,ia,ie,ii)=modelin if(minh+ih*dh.gt.hmag1.and.minh+(ih-1)*dh.lt.hmag2)then if(modelin.gt.maxmod)maxmod=modelin 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(