subroutine huygens(nbins,radius)
      real nbins,radius
      include 'huygens.inc'
      lambda = 532e-9      ! wavelength (m)
      apradius = 5e-6      ! radius of circular aperture (m)
      sdistance = 1e-2     ! distance from aperture to screen (m)
      if (radius.gt.0) then
        sradius = radius
      else
        sradius = 1.0e-2 
      endif
      idN=10
      call hbook2(idN,'integral stats',nbins,-sradius*1e3,sradius*1e3,
     +                            nbins,-sradius*1e3,sradius*1e3,0)
      call hbook2(idN+1,'real part sum',nbins,-sradius*1e3,sradius*1e3,
     +                            nbins,-sradius*1e3,sradius*1e3,0)
      call hbook2(idN+2,'real part var',nbins,-sradius*1e3,sradius*1e3,
     +                            nbins,-sradius*1e3,sradius*1e3,0)
      call hbook2(idN+3,'imag part sum',nbins,-sradius*1e3,sradius*1e3,
     +                            nbins,-sradius*1e3,sradius*1e3,0)
      call hbook2(idN+4,'imag part var',nbins,-sradius*1e3,sradius*1e3,
     +                            nbins,-sradius*1e3,sradius*1e3,0)
      end

      subroutine huygen(ncounts)
      include 'huygens.inc'
      double precision x(2),xx(2)
      double precision dxmag,dxmag2
      real k,omega
      real strength
      real Rpart,Ipart
      real rndm(4)
      k = 2*pi/lambda
      omega = k*clight
      strength = sdistance/(2*pi)*(pi*apradius**2)
      do n=1,ncounts
    1   call ranlux(rndm,4)
        x(1) = rndm(1)*2-1
        x(2) = rndm(2)*2-1
        if (x(1)*x(1)+x(2)*x(2).gt.1) goto 1
        x(1) = x(1)*apradius
        x(2) = x(2)*apradius
        xx(1) = (rndm(3)*2-1)*sradius
        xx(2) = (rndm(4)*2-1)*sradius
        dxmag2 = (x(1)-xx(1))**2+(x(2)-xx(2))**2+sdistance**2
        dxmag = sqrt(dxmag2)
        dxmag3 = dxmag2*dxmag
        Rpart = strength*(cos(k*dxmag)/dxmag3+k*sin(k*dxmag)/dxmag2)
        Ipart = strength*(cos(k*dxmag)/dxmag3-k*cos(k*dxmag)/dxmag2)
        x1 = xx(1)*1e3
        x2 = xx(2)*1e3
        call hfill(idN,x1,x2,1.)
        call hfill(idN+1,x1,x2,Rpart)
        call hfill(idN+2,x1,x2,Rpart*Rpart)
        call hfill(idN+3,x1,x2,Ipart)
        call hfill(idN+4,x1,x2,Ipart*Ipart)
      enddo
      end

      subroutine huygenspro(id2,id1,count)
      integer id2,id1,count
      real x,y
      do i=1,count
        call hrndm2(id2,x,y)
        r=sqrt(x**2+y**2)
        call hfill(id1,r,0.,1/r)
      enddo
      end