subroutine optics(xsize,ysize,xdive,ydive,flength,llength)
      implicit none
      real xsize,ysize,xdive,ydive,flength,llength
      include 'optics.inc'
      data plotrays/.false./
      data idnt/10/
      data lunnt/85/
      data is_open/.false./
      data nament,filent/'optics','optics.hbook'/
      data ntform/'event:i,r0(2):r,dcos0(2):r,r1(2):r,dcos1(2):r'/
      data dist_laser_lens1/0.20/   !m
      data lens1_thickness/0.0142/  !m
      data lens1_radius/0.0250/     !m
      data lens1_rcurve/0.0345/     !m
      data dist_lens1_lens2/0.225/  !m
c      data lens2_thickness/0.0078/  !m
c      data lens2_radius/0.0250/     !m
c      data lens2_rcurve/0.0690/     !m  (15 cm f.l.)
      data lens2_thickness/0.0066/  !m
      data lens2_radius/0.0250/     !m
      data lens2_rcurve/0.0920/     !m  (20 cm f.l.)
      data dist_lens2_lens3/0.175/  !m
      data lens3_thickness/0.0142/  !m
      data lens3_radius/0.0254/     !m
      data lens3_rcurve/0.0345/     !m
      data dist_laser_sample/0.675/ !m
      data lens_ref_index/1.566/
      vector ioptics(1)
      if (ioptics(1).eq.7654321) then
        is_open=.true.
      else
        ioptics(1)=7654321
      endif
      ssizefwhm(1)=xsize
      ssizefwhm(2)=ysize
      ssizerms(1)=3e-3
      ssizerms(2)=3e-3
      sdiverms(1)=xdive
      sdiverms(2)=ydive
      dist_laser_sample=flength
      dist_lens1_lens2=llength
      print *, 'Welcome to the UV laser ablation optics simulation'
      print *, 'Source size (x * y):',ssizefwhm(1),' *',ssizefwhm(2)
      print *, 'Source divergence (x * y):',sdiverms(1),' *',sdiverms(2)
      print *, 'Distance from source to first lens:',dist_laser_lens1
      print *, 'Distance from first lens to second lens:',
     &          dist_lens1_lens2
      print *, 'Distance from second lens to third lens:',
     &          dist_lens2_lens3
      print *, 'Radius of curvature of first lens:',lens1_rcurve
      print *, 'Radius of curvature of second lens:',lens2_rcurve
      print *, 'Radius of curvature of third lens:',lens3_rcurve
      print *, 'Diameter of first lens:',2*lens1_radius
      print *, 'Diameter of second lens:',2*lens2_radius
      print *, 'Diameter of third lens:',2*lens3_radius
      print *, 'Refractive index of lenses:',lens_ref_index
      print *, 'Distance from laser to sample:',dist_laser_sample
      print *, '(all distances above are in meters)'
      end

      subroutine simulate(nsim)
      implicit none
      integer nsim
      include 'optics.inc'
      integer icycle
      call makehists
      do event=1,nsim
        call genevent(event)
        call trackevent(event)
        call hfnt(idnt)
      enddo
      call hrout(0,icycle,' ')
      end

      subroutine makehists
      implicit none
      include 'optics.inc'
      character*80 title
      integer istat
      logical hexist
      external hexist
      if (.not.is_open) then 
        call hropen(lunnt,nament,filent,'n',1024,istat)
        if (istat.ne.0) then
          print *, 'Error opening ntuple output file:',istat
          return
        else
          is_open=.true.
        endif
      endif
      call hrin(0,99999,0)
      do while (hexist(idnt))
        idnt=idnt+1
      enddo
      write(title,*) 'Ablation laser optics simulation with two lenses'
      call hbnt(idnt,title,' ')
      call hbname(idnt,'ray',event,ntform)
      end

      subroutine genevent(ievent)
      implicit none
      integer ievent
      include 'optics.inc'
      real runif(2),rgauss(4)
      external ranlux
      call ranlux(runif,2)
      call rnormx(rgauss,4,ranlux)
      r0(1)=(runif(1)-0.5)*ssizefwhm(1)+rgauss(1)*ssizerms(1)
      r0(2)=(runif(2)-0.5)*ssizefwhm(2)+rgauss(2)*ssizerms(2)
      dcos0(1)=rgauss(3)*sdiverms(1)
      dcos0(2)=rgauss(4)*sdiverms(2)
      r1(1)=9999.
      r1(2)=9999.
      dcos1(1)=9999.
      dcos1(2)=9999.
      end

      subroutine trackevent(ievent)
      implicit none
      integer ievent
      include 'optics.inc'
      real nhat(3),khat(3),r(3)
      real t,rabs,rperp,rdotk,dist
      real track1(10000),track2(10000),track3(1000)
      integer nstep,n
      real discrim

c Propagate the ray from the laser aperture to the first lens entrance surface.
c The lens has a planar entrance surface and a spherical exit surface.
      khat(1)=dcos0(1)
      khat(2)=dcos0(2)
      khat(3)=sqrt(1-khat(1)**2-khat(2)**2)
      r(1)=r0(1)
      r(2)=r0(2)
      r(3)=0
      if (plotrays) then
        do nstep=1,99999
          track1(nstep)=r(1)+khat(1)*nstep/1000.
          track2(nstep)=r(2)+khat(2)*nstep/1000.
          track3(nstep)=r(3)+khat(3)*nstep/1000.
          if (track3(nstep).gt.dist_laser_lens1) goto 1
        enddo
    1   continue
        print *, 'hit first lens at step',nstep
      endif
      r(1)=r(1)+dist_laser_lens1*khat(1)/khat(3)
      r(2)=r(2)+dist_laser_lens1*khat(2)/khat(3)
      r(3)=lens1_rcurve-lens1_thickness
c Check if the ray hits the first lens, otherwise let it go straight
      rperp=sqrt(r(1)**2+r(2)**2)
      if (rperp.lt.lens1_radius) then
c Refract the ray through the entrance surface into the first lens
        nhat(1)=0
        nhat(2)=0
        nhat(3)=1.0
        call refractray(khat,nhat,1.,lens_ref_index)
c Propagate the ray inside the lens
        rdotk=r(1)*khat(1)+r(2)*khat(2)+r(3)*khat(3)
        t=sqrt(lens1_rcurve**2-rperp**2-r(3)**2+rdotk**2)-rdotk
        if (plotrays) then
          do n=1,99999
            track1(nstep+n)=r(1)+khat(1)*n/1000.
            track2(nstep+n)=r(2)+khat(2)*n/1000.
            track3(nstep+n)=track3(nstep)+khat(3)*n/1000.
            if (track3(nstep+n).gt.track3(nstep)+t*khat(3)) goto 2
          enddo
    2     nstep=nstep+n
          print *, 'exit first lens at step',nstep
        endif
        r(1)=r(1)+t*khat(1)
        r(2)=r(2)+t*khat(2)
        r(3)=r(3)+t*khat(3)
c Refract the ray through the exit surface of the lens
        rabs=sqrt(r(1)**2+r(2)**2+r(3)**2)
        nhat(1)=r(1)/rabs
        nhat(2)=r(2)/rabs
        nhat(3)=r(3)/rabs
        call refractray(khat,nhat,lens_ref_index,1.)
c Check for special case of total internal reflection
        if (khat(3).lt.0) then
          return
        endif
      endif
      r(3)=r(3)-lens1_rcurve+lens1_thickness

c Propagate the ray from the first to the second lens
      r(3)=r(3)-dist_lens1_lens2
      r(3)=r(3)-lens2_rcurve+lens2_thickness
      rdotk=r(1)*khat(1)+r(2)*khat(2)+r(3)*khat(3)
      discrim=lens2_rcurve**2-r(1)**2-r(2)**2-r(3)**2+rdotk**2
      if (discrim.gt.0) then
        t=-sqrt(discrim)-rdotk
      else
        return
      endif 
      if (plotrays) then
          do n=1,99999
            track1(nstep+n)=r(1)+khat(1)*n/1000.
            track2(nstep+n)=r(2)+khat(2)*n/1000.
            track3(nstep+n)=track3(nstep)+khat(3)*n/1000.
            if (track3(nstep+n).gt.track3(nstep)+t*khat(3)) goto 3
          enddo
    3     nstep=nstep+n
          print *, 'hit second lens at step',nstep
      endif
      r(1)=r(1)+t*khat(1)
      r(2)=r(2)+t*khat(2)
      r(3)=r(3)+t*khat(3)
      rperp=sqrt(r(1)**2+r(2)**2)
      if (rperp.lt.lens2_radius) then
c Propagate the ray into the second lens.
c The lens has a spherical entrance surface and
c a planar exit surface, opposite to the first lens.
        rabs=sqrt(r(1)**2+r(2)**2+r(3)**2)
        nhat(1)=-r(1)/rabs
        nhat(2)=-r(2)/rabs
        nhat(3)=-r(3)/rabs
        call refractray(khat,nhat,1.,lens_ref_index)
c Propagate the ray inside the second lens
        dist=lens2_thickness-lens2_rcurve-r(3)
        if (plotrays) then
          do n=1,99999
            track1(nstep+n)=r(1)+khat(1)*n/1000.
            track2(nstep+n)=r(2)+khat(2)*n/1000.
            track3(nstep+n)=track3(nstep)+khat(3)*n/1000.
            if (track3(nstep+n).gt.track3(nstep)+dist) goto 4
          enddo
    4     nstep=nstep+n
          print *, 'exit second lens at step',nstep
        endif
        r(1)=r(1)+dist*khat(1)/khat(3)
        r(2)=r(2)+dist*khat(2)/khat(3)
c Refract the ray through the exit surface of the second lens
        nhat(1)=0
        nhat(2)=0
        nhat(3)=1
        call refractray(khat,nhat,lens_ref_index,1.)
c Check for special case of total internal reflection in second lens
        if (khat(3).lt.0) then
          return
        endif
      endif

c Propagate the ray from the second lens to the third
      if (plotrays) then
        do n=1,99999
          track1(nstep+n)=r(1)+khat(1)*n/1000.
          track2(nstep+n)=r(2)+khat(2)*n/1000.
          track3(nstep+n)=track3(nstep)+khat(3)*n/1000.
          if (track3(nstep+n).gt.track3(nstep)+dist_lens2_lens3) goto 5
        enddo
    5   nstep=nstep+n
        print *, 'hit third lens at step',nstep
      endif
      r(1)=r(1)+dist_lens2_lens3*khat(1)/khat(3)
      r(2)=r(2)+dist_lens2_lens3*khat(2)/khat(3)
      r(3)=lens3_rcurve-lens3_thickness
c Check if the ray hits the third lens, otherwise let it go straight
      rperp=sqrt(r(1)**2+r(2)**2)
      if (rperp.lt.lens3_radius) then
c Refract the ray through the entrance surface into the third lens
        nhat(1)=0
        nhat(2)=0
        nhat(3)=1.0
        call refractray(khat,nhat,1.,lens_ref_index)
c Propagate the ray inside the lens
        rdotk=r(1)*khat(1)+r(2)*khat(2)+r(3)*khat(3)
        t=sqrt(lens3_rcurve**2-rperp**2-r(3)**2+rdotk**2)-rdotk
        if (plotrays) then
          do n=1,99999
            track1(nstep+n)=r(1)+khat(1)*n/1000.
            track2(nstep+n)=r(2)+khat(2)*n/1000.
            track3(nstep+n)=track3(nstep)+khat(3)*n/1000.
            if (track3(nstep+n).gt.track3(nstep)+t*khat(3)) goto 6
          enddo
    6     nstep=nstep+n
          print *, 'exit third lens at step',nstep
        endif
        r(1)=r(1)+t*khat(1)
        r(2)=r(2)+t*khat(2)
        r(3)=r(3)+t*khat(3)
c Refract the ray through the exit surface of the lens
        rabs=sqrt(r(1)**2+r(2)**2+r(3)**2)
        nhat(1)=r(1)/rabs
        nhat(2)=r(2)/rabs
        nhat(3)=r(3)/rabs
        call refractray(khat,nhat,lens_ref_index,1.)
c Check for special case of total internal reflection
        if (khat(3).lt.0) then
          return
        endif
      endif

c Propagate the ray from third lens to the sample plane
      r(3)=r(3)-lens3_rcurve+lens3_thickness
      dist=dist_laser_sample-dist_laser_lens1
     +    -dist_lens1_lens2-dist_lens2_lens3-r(3)
      if (plotrays) then
        do n=1,99999
          track1(nstep+n)=r(1)+khat(1)*n/1000.
          track2(nstep+n)=r(2)+khat(2)*n/1000.
          track3(nstep+n)=track3(nstep)+khat(3)*n/1000.
          if (track3(nstep+n).gt.track3(nstep)+dist) goto 7
        enddo
    7   nstep=nstep+n
        print *, 'hit the sample at step',nstep
      endif
      r1(1)=r(1)+dist*khat(1)/khat(3)
      r1(2)=r(2)+dist*khat(2)/khat(3)
      dcos1(1)=khat(1)
      dcos1(2)=khat(2)
      if (plotrays) then
        call iselnt(10)
        call ipl(nstep,track3,track1)
      endif
      end

      subroutine refractray(khat,nhat,rindex1,rindex2)
      implicit none
      real khat(3),nhat(3)
      real rindex1,rindex2
c     include 'optics.inc'
      real phat(3),qhat(3)
      real sintheta1,sintheta2,costheta2
      phat(1)=khat(2)*nhat(3)-khat(3)*nhat(2)
      phat(2)=khat(3)*nhat(1)-khat(1)*nhat(3)
      phat(3)=khat(1)*nhat(2)-khat(2)*nhat(1)
      sintheta1=sqrt(phat(1)**2+phat(2)**2+phat(3)**2)
      phat(1)=phat(1)/sintheta1
      phat(2)=phat(2)/sintheta1
      phat(3)=phat(3)/sintheta1
      qhat(1)=nhat(2)*phat(3)-nhat(3)*phat(2)
      qhat(2)=nhat(3)*phat(1)-nhat(1)*phat(3)
      qhat(3)=nhat(1)*phat(2)-nhat(2)*phat(1)
      sintheta2=rindex1*sintheta1/rindex2
      if (sintheta2.lt.1) then
         costheta2=sqrt(1-sintheta2**2)
         khat(1)=nhat(1)*costheta2+qhat(1)*sintheta2
         khat(2)=nhat(2)*costheta2+qhat(2)*sintheta2
         khat(3)=nhat(3)*costheta2+qhat(3)*sintheta2
      else
         khat(1)=0
         khat(2)=0
         khat(3)=-1
      endif
      end