*
* xraymono: calculates the optics transport function for an asymmetric
*           crystal monochromator.  The magnitude of the crystal momentum
*           q for the desired set of crystal planes, and the Darwin width
*           for that q must be calculated by the user.  For help, see url
*
*              http://www.chess.cornell.edu/oldchess/operatns/xrclcdwn.htm
*
* author: Richard Jones, University of Connecticut
* version: March 2, 2007
*
* Model:
*  1. Coordinate system
*     The scattering takes place in the xy plane.  The direction of the
*     incident wave is theta, defined wrt. the x axis by the usual ccw
*     sign convention.  The diffracted wave direction thetaprime is defined
*     in a similar way.
*  2. Selecting the desired planes
*     Only one set of crystal planes is considered in the model.  It is
*     assumed that the geometry of the monochromator excludes significant
*     diffractive yield from any other than the desired set of planes.
*     The desired set of planes is selected by assigning physical values
*     to the variables qmag (keV) and qfwhm (eV).  The qfwhm parameter is
*     the spread in the magnitude qmag that comes about from zero-point
*     motion of the crystal, and is related to the Darwin width by
*         qfwhm = 2 k cos(theta_Bragg) Darwin_width
*     which is independent of k.
*  3. Crystal cut and orientation
*     The angle from the y-hat direction to the outward-pointing unit
*     normal to the crystal surface is denoted beta.  The angle from the
*     y-hat direction to the outward-pointing unit normal to the crystal
*     planes for the selected reflection is denoted alpha.  For symmetric
*     crystals, beta and alpha are the same value.  Their allowed ranges are
*       * alpha: [-pi/2,pi/2] (radians) : crystal planes angle wrt hor.
*       * beta:  [-pi/2,pi/2] (radians) : polished surface angle wrt hor.
 
* Xraymono: set up the monochromator angles to select a particular
*     X-ray energy k (keV) and asymmetry factor b, defined as
*          b = sin(theta_Bragg+delta)/sin(theta_Bragg-delta)
*     where delta = alpha - beta is the angle between q and the surface normal
*     and  sin(theta_Bragg) = q/(2k)
      subroutine Xraymono(k,b)
      real k,b
      double precision pi,alpha,beta,qmag,qfwhm
      common /monoconfig/pi,alpha,beta,qmag,qfwhm
      integer mset,maxset
      real qmagtab,qfwhmtab
      common /msetcom/mset,maxset,qmagtab(12),qfwhmtab(12)
      data mset,maxset/1,12/
      data qmagtab(1),qfwhmtab(1)/3.954352,0.511246/   ! Si (111) planes
      data qmagtab(2),qfwhmtab(2)/11.86306,0.105920/   ! Si (333) planes
      data qmagtab(3),qfwhmtab(3)/6.457430,0.369586/   ! Si (220) planes
      data qmagtab(4),qfwhmtab(4)/9.132186,0.226969/   ! Si (400) planes
      data qmagtab(5),qfwhmtab(5)/9.951568,0.140974/   ! Si (331) planes
      data qmagtab(6),qfwhmtab(6)/7.572008,0.209791/   ! Si (311) planes
      data qmagtab(7),qfwhmtab(7)/6.020636,0.354588/   ! diamond (111) planes
      data qmagtab(8),qfwhmtab(8)/18.06191,0.048427/   ! diamond (333) planes
      data qmagtab(9),qfwhmtab(9)/9.831658,0.190329/   ! diamond (220) planes
      data qmagtab(10),qfwhmtab(10)/13.90406,0.105626/ ! diamond (400) planes
      data qmagtab(11),qfwhmtab(11)/15.15160,0.065005/ ! diamond (331) planes
      data qmagtab(12),qfwhmtab(12)/11.52864,0.101601/ ! diamond (311) planes
      character*8 material
      character*5 planes
      common /msetcomc/material(12),planes(12)
      data material/6*'silicon',6*'diamond'/
      data planes/'(111)','(333)','(220)','(400)','(331)','(311)',
     +            '(111)','(333)','(220)','(400)','(331)','(311)'/
      double precision tandelta
      data pi/0/
      qmag=qmagtab(mset)
      qfwhm=qfwhmtab(mset)
      if (k.lt.q/2) then
         alpha=pi/2
      else
         alpha=asin(qmag/(2*k))
      endif
      tandelta=tan(alpha)*(b-1)/(b+1)
      beta=alpha-atan(tandelta)
      if (pi.eq.0) then
        pi=2*atan2(1.,0.)
        print *, 'X-ray monochromator initialized for ',
     +            material(mset),planes(mset)
        print *, '   * crystal momentum:',qmag,' keV/c'
        print *, '   * X-ray energy:',k,' keV'
        print *, '   * Bragg angle:',alpha,' (',alpha*180/pi,' degrees)'
        print *, '   * surface angle:',beta,' (',beta*180/pi,' degrees)'
        print *, '   * asymmetry parameter:',b
      endif
      end

* rseqset: select a unique random number sequence seq,
*          any 32-bit integer will do.
      integer function rseqset(seq)
      integer seq
      call rluxgo(2,seq,0,0)
      rseqset=seq
      end

* monoset: select the monochromator crystal and set of planes
*          from the tabulated list
      integer function monoset(iset)
      integer iset
      integer mset,maxset
      real qmagtab,qfwhmtab
      common /msetcom/mset,maxset,qmagtab(12),qfwhmtab(12)
      character*8 material
      character*5 planes
      common /msetcomc/material(12),planes(12)
      double precision pi,alpha,beta,qmag,qfwhm
      common /monoconfig/pi,alpha,beta,qmag,qfwhm
      if ((iset.ge.1).and.(iset.le.maxset)) then
        mset=iset
      endif
      print *, 'Monochromator configured for ',
     +          material(mset),planes(mset)
      monoset=mset
      pi=0
      end

* amono: compute thetaprime (radians) from inputs k (keV)
*        and theta (radians), and return the reflectivity
*        of the crystal for the given set of waves.
      real function amono(k,theta,thetaprime)
      real k,theta
      real thetaprime
      double precision pi,alpha,beta,qmag,qfwhm
      common /monoconfig/pi,alpha,beta,qmag,qfwhm
      double precision costhe1,costhe2
      double precision sinthe1,sinthe2
      double precision the1,the2
      double precision dq(2),dq2
      the1=beta-theta
      costhe1=cos(the1)
      sinthe1=sin(the1)
      if (sinthe1.gt.0) then   ! reflect upward
        costhe2=costhe1-(qmag/k)*sin(alpha-beta)
        the2=acos(costhe2)
      else                     ! reflect downward
        costhe2=costhe1+(qmag/k)*sin(alpha-beta)
        the2=-acos(costhe2)
      endif
      sinthe2=sin(the2)
      thetaprime=the2+beta
      if (sinthe1.gt.0) then   ! reflect upward
        dq(1)=k*(costhe2-costhe1)+qmag*sin(alpha-beta)
        dq(2)=k*(sinthe2+sinthe1)-qmag*cos(alpha-beta)
      else                     ! reflect downward
        dq(1)=k*(costhe2-costhe1)-qmag*sin(alpha-beta)
        dq(2)=k*(sinthe2+sinthe1)+qmag*cos(alpha-beta)
      endif
      dq2=(dq(1)**2+dq(2)**2)*1e6
      amono=exp(-0.5*(dq2/(qfwhm/2.3)**2)**3)
      end
      
* theta2: return thetaprime (radians) from inputs k (keV)
*         and theta (radians), with no regard to the size of amono.
      real function theta2(k,theta)
      real k,theta
      real amono
      external amono
      real a,thetaprime
      a=amono(k,theta,thetaprime)
      theta2=thetaprime
      end

* refl: return reflectivity of the monochromator crystal from
*       inputs k (keV) and theta (radians)
      real function refl(k,theta)
      real k,theta
      real amono
      external amono
      real thetaprime
      refl=amono(k,theta,thetaprime)
      end

* MCgen: generate a set of n Monte Carlo rays from a broadband source
*        of rms size sigma (m) a distance D (m) away from the monochromator.
*        A single bounce from the crystal surface is simulated in MCgen.
*        The generated points are weighted by the reflectivity and stored
*        in a ntuple in the current paw working area (memory or disk file).
      integer function MCgen(n,sigma,D)
      integer n
      real sigma,D
      double precision pi,alpha,beta,qmag,qfwhm
      common /monoconfig/pi,alpha,beta,qmag,qfwhm
      integer id
      character*80 ntform
      real a,k,theta,thetaprime
      common /ntcommon/a,k,theta,thetaprime
      data ntform/'a:R,k:R,theta:R,thetaprime:R'/
      integer events,ngen
      real rvec(4)
      real rho,phi
      real k0,kmin,kmax
      real amono
      external amono
      logical hexist
      external hexist
      k0=qmag/(2*sin(alpha))
      kmin=qmag/(2*sin(alpha+20*sigma/D))
      kmax=qmag/(2*sin(alpha-20*sigma/D))
      id=1
      do while (hexist(id))
        id=id+1
      enddo
      call hbnt(id,'X-ray monochromator with one bounce, asym',' ')
      call hbname(id,'xraymono',a,ntform)
      events=0
      ngen=0
      do while (events.lt.n)
        ngen=ngen+1
        call ranlux(rvec,4)
        k=rvec(1)*kmin+(1-rvec(1))*kmax
        phi=rvec(2)*2*pi
        rho=sqrt(-2*log(rvec(3)))
        theta=rho*cos(phi)*(sigma/D)
        a=amono(k,theta,thetaprime)
        if (a.gt.rvec(4)) then
          call hfnt(id)
          events=events+1
        endif
      enddo
      MCgen=ngen
      end

* MCgen2: generate a set of n Monte Carlo rays from a broadband source
*        of rms size sigma (m) a distance D (m) away from the monochromator.
*        A second bounce from a symmetric crystal that returns the central
*        beam ray to the original horizontal direction is simulated in MCgen2.
*        The generated points are weighted by the reflectivity product and
*        stored in a ntuple in the current paw working area (memory or disk).
      integer function MCgen2(n,sigma,D)
      integer n
      real sigma,D
      double precision pi,alpha,beta,qmag,qfwhm
      common /monoconfig/pi,alpha,beta,qmag,qfwhm
      double precision alpha_saved,beta_saved
      integer id
      character*80 ntform
      real a,k,theta,thetaprime
      common /ntcommon/a,k,theta,thetaprime
      data ntform/'a:R,k:R,theta:R,thetaprime:R'/
      real theta1,theta2
      integer events,ngen
      real rvec(4)
      real rho,phi
      real k0,kmin,kmax
      real amono
      external amono
      logical hexist
      external hexist
      alpha_saved=alpha
      beta_saved=beta
      k0=qmag/(2*sin(alpha))
      kmin=qmag/(2*sin(alpha+20*sigma/D))
      kmax=qmag/(2*sin(alpha-20*sigma/D))
      id=1
      do while (hexist(id))
        id=id+1
      enddo
      call hbnt(id,'X-ray monochromator with two bounces, asym+sym',' ')
      call hbname(id,'xraymono',a,ntform)
      events=0
      ngen=0
      do while (events.lt.n)
        ngen=ngen+1
        call ranlux(rvec,4)
        k=rvec(1)*kmin+(1-rvec(1))*kmax
        phi=rvec(2)*2*pi
        rho=sqrt(-2*log(rvec(3)))
        theta=rho*cos(phi)*(sigma/D)
        a=amono(k,theta,theta1)
        beta=alpha
        a=a*amono(k,theta1,thetaprime)
        if (a.gt.rvec(4)) then
          call hfnt(id)
          events=events+1
        endif
        beta=beta_saved
      enddo
      MCgen2=ngen
      end

* MCgen3: generate a set of n Monte Carlo rays from a broadband source
*        of rms size sigma (m) a distance D (m) away from the monochromator.
*        A second bounce from a symmetric crystal that returns the central
*        beam ray to the original horizontal direction is simulated, followed
*        by a reflection from the diamond (220) plane is simulated in MCgen3.
*        The generated points are weighted by the reflectivity product and
*        stored in a ntuple in the current paw working area (memory or disk).
      integer function MCgen3(n,sigma,D)
      integer n
      real sigma,D
      double precision pi,alpha,beta,qmag,qfwhm
      common /monoconfig/pi,alpha,beta,qmag,qfwhm
      integer mset,maxset
      real qmagtab,qfwhmtab
      common /msetcom/mset,maxset,qmagtab(12),qfwhmtab(12)
      double precision alpha_saved,beta_saved
      integer mset_saved
      integer id
      character*80 ntform
      real a,k,theta,thetaprime
      common /ntcommon/a,k,theta,thetaprime
      data ntform/'a:R,k:R,theta:R,thetaprime:R'/
      real theta1,theta2
      integer events,ngen
      real rvec(5)
      real rho,phi
      real k0,kmin,kmax
      real alfa
      equivalence (alfa,alpha)
      real amono
      external amono
      logical hexist
      external hexist
      alpha_saved=alpha
      beta_saved=beta
      mset_saved=mset
      k0=qmag/(2*sin(alpha))
      kmin=qmag/(2*sin(alpha+20*sigma/D))
      kmax=qmag/(2*sin(alpha-20*sigma/D))
      id=1
      do while (hexist(id))
        id=id+1
      enddo
      call hbnt(id,
     +  'X-ray monochromator with three bounces, asym+sym+diamond',' ')
      call hbname(id,'xraymono',a,ntform)
      call hbname(id,'rocking',alfa,'alpha:R*8')
      events=0
      ngen=0
      do while (events.lt.n)
        ngen=ngen+1
        call ranlux(rvec,5)
        k=rvec(1)*kmin+(1-rvec(1))*kmax
        phi=rvec(2)*2*pi
        rho=sqrt(-2*log(rvec(3)))
        theta=rho*cos(phi)*(sigma/D)
        a=amono(k,theta,theta1)
        beta=alpha
        a=a*amono(k,theta1,theta2)
        mset=9               ! switch to diamond (220)
        call xraymono(k0,1.) ! symmetric geometry
        alpha=alpha+(rvec(4)-0.5)*1e-3
        a=a*amono(k,theta2,thetaprime)
        if (a.gt.rvec(5)) then
          call hfnt(id)
          events=events+1
        endif
        alpha=alpha_saved
        beta=beta_saved
        mset=mset_saved
        qmag=qmagtab(mset)
        qfwhm=qfwhmtab(mset)
      enddo
      MCgen3=ngen
      end

* MCgen4: generate a set of n Monte Carlo rays from a broadband source
*        of rms size sigma (m) a distance D (m) away from the monochromator.
*        A second bounce from a symmetric crystal that returns the central
*        beam ray to the original horizontal direction is simulated, followed
*        by a second double-bounce monochromator configured for silicon (220).
*        The generated points are weighted by the reflectivity product and
*        stored in a ntuple in the current paw working area (memory or disk).
      integer function MCgen4(n,sigma,D)
      integer n
      real sigma,D
      double precision pi,alpha,beta,qmag,qfwhm
      common /monoconfig/pi,alpha,beta,qmag,qfwhm
      integer mset,maxset
      real qmagtab,qfwhmtab
      common /msetcom/mset,maxset,qmagtab(12),qfwhmtab(12)
      double precision alpha_saved,beta_saved
      integer mset_saved
      integer id
      character*80 ntform
      real a,k,theta,thetaprime
      common /ntcommon/a,k,theta,thetaprime
      data ntform/'a:R,k:R,theta:R,thetaprime:R'/
      real theta1,theta2
      integer events,ngen
      real rvec(4)
      real rho,phi
      real k0,kmin,kmax
      real amono
      external amono
      logical hexist
      external hexist
      alpha_saved=alpha
      beta_saved=beta
      mset_saved=mset
      k0=qmag/(2*sin(alpha))
      kmin=qmag/(2*sin(alpha+20*sigma/D))
      kmax=qmag/(2*sin(alpha-20*sigma/D))
      id=1
      do while (hexist(id))
        id=id+1
      enddo
      call hbnt(id,
     + 'X-ray monochromator with four bounces, asym+sym+sym+sym',' ')
      call hbname(id,'xraymono',a,ntform)
      events=0
      ngen=0
      do while (events.lt.n)
        ngen=ngen+1
        call ranlux(rvec,4)
        k=rvec(1)*kmin+(1-rvec(1))*kmax
        phi=rvec(2)*2*pi
        rho=sqrt(-2*log(rvec(3)))
        theta=rho*cos(phi)*(sigma/D)
        a=amono(k,theta,theta1)
        beta=alpha
        a=a*amono(k,theta1,theta2)
        mset=3               ! switch to silicon (220)
        call xraymono(k0,1.) ! symmetric geometry
        a=a*amono(k,theta2,theta1)
        a=a*amono(k,theta1,thetaprime)
        if (a.gt.rvec(4)) then
          call hfnt(id)
          events=events+1
        endif
        alpha=alpha_saved
        beta=beta_saved
        mset=mset_saved
        qmag=qmagtab(mset)
        qfwhm=qfwhmtab(mset)
      enddo
      MCgen4=ngen
      end

* MCgen5: generate a set of n Monte Carlo rays from a broadband source
*        of rms size sigma (m) a distance D (m) away from the monochromator.
*        A second bounce from a symmetric crystal that returns the central
*        beam ray to the original horizontal direction is simulated, followed
*        by a second double-bounce monochromator configured for silicon (220)
*        and then followed by a final reflection from the diamond (220) plane.
*        The generated points are weighted by the reflectivity product and
*        stored in a ntuple in the current paw working area (memory or disk).
      integer function MCgen5(n,sigma,D)
      integer n
      real sigma,D
      double precision pi,alpha,beta,qmag,qfwhm
      common /monoconfig/pi,alpha,beta,qmag,qfwhm
      integer mset,maxset
      real qmagtab,qfwhmtab
      common /msetcom/mset,maxset,qmagtab(12),qfwhmtab(12)
      double precision alpha_saved,beta_saved
      integer mset_saved
      integer id
      character*80 ntform
      real a,k,theta,thetaprime
      common /ntcommon/a,k,theta,thetaprime
      data ntform/'a:R,k:R,theta:R,thetaprime:R'/
      real theta1,theta2
      integer events,ngen
      real rvec(5)
      real rho,phi
      real k0,kmin,kmax
      real alfa
      equivalence (alfa,alpha)
      real amono
      external amono
      logical hexist
      external hexist
      alpha_saved=alpha
      beta_saved=beta
      mset_saved=mset
      k0=qmag/(2*sin(alpha))
      kmin=qmag/(2*sin(alpha+20*sigma/D))
      kmax=qmag/(2*sin(alpha-20*sigma/D))
      id=1
      do while (hexist(id))
        id=id+1
      enddo
      call hbnt(id,
     + 'X-ray monochromator with five bounces, asym+sym+sym+sym+diamond'
     + ,' ')
      call hbname(id,'xraymono',a,ntform)
      call hbname(id,'rocking',alfa,'alpha:R*8')
      events=0
      ngen=0
      do while (events.lt.n)
        ngen=ngen+1
        call ranlux(rvec,5)
        k=rvec(1)*kmin+(1-rvec(1))*kmax
        phi=rvec(2)*2*pi
        rho=sqrt(-2*log(rvec(3)))
        theta=rho*cos(phi)*(sigma/D)
        a=amono(k,theta,theta1)
        beta=alpha
        a=a*amono(k,theta1,theta2)
        mset=3               ! switch to silicon (220)
        call xraymono(k0,1.) ! symmetric geometry
        a=a*amono(k,theta2,theta1)
        a=a*amono(k,theta1,theta2)
        mset=9               ! switch to diamond (220)
        call xraymono(k0,1.) ! symmetric geometry
        alpha=alpha+(rvec(4)-0.5)*1e-3
        a=a*amono(k,theta2,thetaprime)
        if (a.gt.rvec(5)) then
          call hfnt(id)
          events=events+1
        endif
        alpha=alpha_saved
        beta=beta_saved
        mset=mset_saved
        qmag=qmagtab(mset)
        qfwhm=qfwhmtab(mset)
      enddo
      MCgen5=ngen
      end