* phi_intens.kumac - plot the intensity of the GlueX photon
*                    as a function of theta*E/m and phi of the
*                    bremsstrahlung emission.  The intensity is
*                    plotted for the primary peak at 9 GeV and
*                    an electron beam energy of 12 GeV.
macro make id=100 Emax=12 Epeak=9.0 emit=1e-8 dist=75 coldiam=3.4e-3 cur=2.8
  call cobrems.f77($rsigma([Emax]),$rsigma([Epeak]),$rsigma([emit]),$rsigma([dist]),$rsigma([coldiam]),0)
  id1=[id]+1
  id2=[id]+2
  xpeak=$rsigma([Epeak]/[Emax])
  thetamax=$rsigma([coldiam]/2/[dist]*[Emax]/0.511e-3)
  xmin=$rsigma([xpeak]/([xpeak]+(1-[xpeak])*([thetamax]**2+1)))
  message xpeak,xmin,thetamax are [xpeak],[xmin],[thetamax]
  I0=$rsigma([cur]/1.6e-13)
  if ($hexist([id1]).eq.0 .or. $hexist([id2]).eq.0) then
    message 'Wait, this can take a minute or two to compute...'
    fun2 [id] dNcdxdp([xpeak]/([xpeak]+(1-[xpeak])*(x**2+1)),y)*2*(1-[xpeak])/[xpeak]*([xpeak]/([xpeak]+(1-[xpeak])*(x**2+1)))**2*[I0] 50 0. [thetamax] 50 0. 3.1416
    * make a copy with a shorter title, otherwise segfault!
    h/copy [id] [id1] 'dNcdxdp(theta,phi)'
    fun2 [id] dNidxdt2(sqrt([xpeak]*[xmin]),x**2)*2*([xpeak]-[xmin])/6.2832*[I0] 100 0 $sigma(2*[thetamax]) 50 0 3.1416
    h/copy [id] [id2] 'dNidxdp(theta,phi)'
  endif
  v/del hin*
  v/create hin1(50,50)
  get/cont [id1] hin1
  v/create hin2(100,50)
  get/cont [id2] hin2
  do r=1,50
    vadd hin1(:,[r]) hin2(:,[r]) hin2(:,[r])
  enddo
  *sigma hin2=hin2*([Emax]/([dist]*1e3*5.11e-4))**2
  put/cont [id] hin2
  opt utit
  min [id] 0
  h/pl [id](:50,:) surf
  atit 'theta (m/E)' 'phi (rad)' 'beam intensity (/mm^2!/s)'
  opt htit
return