
       program myOHmodel
       character*64 fname,gname
c load traj file
       write(*,'("enter name of trajectory file")')
       read(*,101)fname
       open(30,file=fname,status='old')
c load output file
       write(*,'("enter name of output file")')
       read(*,101)gname
       open(31,file=gname,status='new')
101    format(a)
c initialize functions used
       fnz=1.e-10
1      continue
c read traj point
       read (30,*,end=99)x,y,z,day
       write(*,*)x,y,z,day
       rho=sqrt(x**2+y**2)
       r=rho
       fnr=1.e-10
c calculate density from analytic approx
       if(r.lt.1.75)go to 12	      
       chk=(r-4.)**2/2. 
       if(chk.gt.10.)chk=10.
       if(r.lt.4.0)h=.37*exp(-chk)
       if(r.ge.4.0)h=.37
       chk=(z/h)**2
       if(chk.gt.10.)chk=10.
c       fnz=2.30*exp(-(z/h)**2)
        fnz=2.30*exp(-chk)
       if(r.lt.2.0)h2=.5
       if(r.ge.2.0)h2=.5*(r/2.)**.5
       chk=(z/h2)**2
       if(chk.gt.10.)chk=10.
c       fnz=fnz+.70*exp(-(z/h2)**2)
       fnz=fnz+.70*exp(-chk)
       rpk=3.8
       if(r.lt.rpk)then
       chk=(r-rpk)**2/(1.**2)
       if(chk.gt.10.)chk=10.
c      fnr=3.00*exp(-(r-rpk)**2/(1.**2))
       fnr=3.00*exp(-chk)
       else
       chk=(r-rpk)**2/(.7**2)
       if(chk.gt.10.)chk=10.
       fnr=2.30*exp(-chk)
c      fnr=2.30*exp(-(r-rpk)**2/(.7**2))
       chk=(r-rpk)**2/(1.)
       if(chk.gt.10.)chk=10.
       fnr=fnr+.70*exp(-chk)
c      fnr=fnr+.70*exp(-(r-rpk)**2/(1.))
       endif
12     continue
       dens=fnr*fnz
c write to output file
       write(31,*)x,y,z,day,dens
c      write(*,*)x,y,z,day,dens
       go to 1
99     continue
       stop
       end
