c sag - Compute the sag in a silicon plate under its own weight if it is
c       supported at the same height by its two ends.  A crude estimate
c       for the fundamental oscillation frequency of the plate in this
c       geometry is sqrt(g/sag).  The plate is L long along the direction
c       spanning the gap between the supports, W wide, and T thick.
c       The differential equation obeyed by this unloaded plate is
c
c                             12 g rho 
c                  y'''' = - -----------
c                              B T^2
c
c       where g is the gravitational acceleration, B is Young's modulus,
c       and rho is the density of the plate material.  Using the fact that
c       the plate is symmetric about the midpoint, put x=0 there and then
c       solve only on the interval [0,L/2).  The boundary conditions are
c
c       1) y(L/2) = 0
c       2) y'(0) = 0
c       3) y''(L/2) = 0
c       4) y'''(0) = 0
c
c       The general solution is y(x) = c0 + c2 x^2 + c4 x^4, with
c
c               g rho                3 g rho L^2              5 g rho L^4
c       c4 = - --------    ,  c2 =  ------------  ,   c0 = - -------------
c               2 B T^2                4 B T^2                  32 B T^2
c
c
c author:   Richard Jones, University of Connecticut
c           May 3, 2009

      subroutine sag(Wcm,Tcm,Lcm)
      real Wcm,Tcm,Lcm        ! all distances in cm
      real W                  ! width of the plates (m)
      real T                  ! thickness of each plate (m)
      real L                  ! length of the plates (m)
      real coef(5)            ! polynomial coefficients for solution y(x)
      real Darwin             ! Darwin width (FWHM) of relevant vector
      common /sagging/W,T,L,coef,Darwin
      data L,T,W/5.0e-2,1.0e-3,7.0e-2/
      data Darwin/5.0e-6/
      real g
      parameter (g=9.8)     ! gravitational acceleration (m/s^2)
      real rho
      parameter (rho=2330)  ! density of silicon (kg/m^3)
      real B                ! Young's modulus of the plate materials (Pascals)
      parameter (B=130e9)   ! for silicon
      real sagita,bow
      if (Wcm.gt.0) then
        W = Wcm/100
      endif
      if (Tcm.gt.0) then
        T = Tcm/100
      endif
      if (Lcm.gt.0) then
        L = Lcm/100
      endif
      print *, 'Silicon plate sag calculation with plate thickness',
     +          T*1e3,' mm,'
      print *, ' width',W*1e3,' mm, and length',L*1e3,' mm.'
      coef(1) = -5*rho*g*L**4/(32*B*T**2)
      coef(2) = 0
      coef(3) = +3*rho*g*L**2/(4*B*T**2)
      coef(4) = 0
      coef(5) = -rho*g/(2*B*T**2)
      sagita = -coef(1)
      bow = coef(2) + 2*coef(3)*(L/2) +
     +                3*coef(4)*(L/2)**2 +
     +                4*coef(5)*(L/2)**3
      print *, 'sag in middle of silicon plate is ',sagita,' m'
      print *, 'fundamental oscillation frequency is ',
     +          sqrt(g/sagita)/(2*3.1415926),' Hz'
      print *, 'bow angle is +/-',bow,' radians'
      print *, 'Note: silicon 331 Darwin width is ',Darwin,
     +         ' microradians at 15 keV.'
      end

      real function saggy(x)  ! shape of saggy plate (m)
      real x                  ! as a function of distance from center (m)
      real W                  ! width of the plates (m)
      real T                  ! thickness of each plate (m)
      real L                  ! length of the plates (m)
      real coef(5)            ! polynomial coefficients for solution y(x)
      real Darwin             ! Darwin width (FWHM) of relevant vector
      common /sagging/W,T,L,coef,Darwin
      saggy = coef(1)+coef(2)*x+coef(3)*x**2+coef(4)*x**3+coef(5)*x**4
      end

      real function saggy1(x) ! slope of saggy plate (radians)
      real x                  ! as a function of distance from center (m)
      real W                  ! width of the plates (m)
      real T                  ! thickness of each plate (m)
      real L                  ! length of the plates (m)
      real coef(5)            ! polynomial coefficients for solution y(x)
      real Darwin             ! Darwin width (FWHM) of relevant vector
      common /sagging/W,T,L,coef,Darwin
      saggy1 = coef(2)+2*coef(3)*x+3*coef(4)*x**2+4*coef(5)*x**3
      end

      subroutine strained(s1,s2) ! set the strain vectors of the mono
      real s1(2)                 ! " for the first crystal (radians/m)
      real s2(2)                 ! " for the second crystal (radians/m)
      real strain(2,2)
      common /straining/strain
      strain(1,1) = s1(1)
      strain(2,1) = s1(2)
      strain(1,2) = s2(1)
      strain(2,2) = s2(2)
      end

      real function dintens(x,y)  ! diffractive intensity of strained mono
      real x,y                    ! as a function of position on crystal (m)
      real strain(2,2)
      common /straining/strain
      real W
      real T
      real L
      real coef(5)
      real Darwin
      common /sagging/W,T,L,coef,Darwin
      real dir(2,2)
      real s1,s2
      real diff
      s1 = sqrt(strain(1,1)**2+strain(2,1)**2)+1e-20
      s2 = sqrt(strain(1,2)**2+strain(2,2)**2)+1e-20
      dir(1,1) = (strain(1,1)*x+strain(2,1)*y)*strain(1,1)/s1
      dir(2,1) = (strain(1,1)*x+strain(2,1)*y)*strain(2,1)/s1
      dir(1,2) = (strain(1,2)*x+strain(2,2)*y)*strain(1,2)/s2
      dir(2,2) = (strain(1,2)*x+strain(2,2)*y)*strain(2,2)/s2
      dintens = exp(-0.5*
     +              (((dir(1,1)-dir(1,2))**2+(dir(2,1)-dir(2,2))**2)
     +               /(Darwin/2.3)**2)**2)
      end