source: PSPA/parmelaPSPA/trunk/bend.f @ 405

      subroutine bend(dwtp,iend)
c---second order bending magnet transformation for uniform
c   field magnet with curved pole faces. the edge transformations
c   are form pages 71 to 75 of slac report No. 75  by Karl K
c   brown may 1969. T216 and T436 are not used because del=0
c   rho is calculated for each particle.
c   the transformation in the magnet is exact.
c----------------------------------------------------------------------
      save
c
      include 'param_sz.h'
      include 'constcom.h'
      include 'pcordcom.h'
      include 'syscom.h'
      include 'ucom.h'
c
c--------------------------------------------------------------------------
c*
      ne=rne
      if(x**2+y**2.gt.el(2,ne)**2) go to 49
      gr=el(4,ne)/erest
      bg = sqrt(gr*(2.+gr))
      rho = el(1,ne)/el(5,ne)
      bz=bgz/gamma
      xp=bgx/bgz
      theta1=atan(xp)
      yp=bgy/bgz
      rhoa=rho*sqrt(bgz**2+bgx**2)/bg
c---check for particle being at entrance edge
      if (ne.gt.1) go to 20
      if (z.gt.0.) go to 30
   10 continue
      con=1.+xp**2+yp**2
      b1=el(6,ne)
      b2=b1-el(8,ne)
      r1=el(10,ne)
      h=1/rhoa
      if(r1.eq.0)then
      hr1=0
      else
      hr1=1./r1
      endif
      xt=x+.5*h*(y*y/cos(b1)**2-x*x*tan(b1)**2)
c      r11         t133               t111
      xpt=xp+x*h*tan(b1)+x*x*.5*h*hr1/cos(b1)**3+x*xp*h*tan(b1)**2
c         r22     r21      t211                    t212
      xpt=xpt+y*y*(h*h*(.5+tan(b1)**2)*tan(b1)-.5*h*hr1/cos(b1)**3)
c                  t233
      xpt=xpt-y*yp*h*tan(b1)**2
c                  t234
      yt=y+x*y*h*tan(b1)**2
c       r33   t313
      ypt=yp-h*y*tan(b2)-x*y*h*hr1/cos(b1)**3
c        r44    r43        t413
      yp=ypt-x*yp*h*tan(b1)**2-xp*y*h/cos(b1)**2
c               t414               t423
      x=xt
      y=yt
      xp=xpt
      bgz=bgz*sqrt(con/(1.+xp**2+yp**2))
      bz=bgz/gamma
      bgx=xp*bgz
      rhoa=rho*sqrt(bgz**2+bgx**2)/bg
      theta1=atan(xp)
      go to 30
   20 continue
      if (z.eq.zloc(ne-1)) go to 10
   30 continue
      dz=bz*dcon*dwtp
      ang=dz/(rho+x)
      dz=ang*rho
      zn=z+dz
      if(zn.lt.zloc(ne))go to 35
      zn=zloc(ne)
      dwtp=dwtp*(zn-z)/dz
      ang=ang*(zn-z)/dz
   35 continue
      xx=x
      con=1.+xp**2+yp**2
      arg=sin(theta1+ang)-(xx+rho)*sin(ang)/rhoa
      if(abs(arg).ge.1)go to 50
      theta2=asin(arg)
      x=-rho+cos(ang)*(xx+rho)+rhoa*(cos(theta2)-cos(theta1+ang))
      y=y+yp*rhoa*(ang+theta1-theta2)*cos(theta1)
      yp=yp*cos(theta1)/cos(theta2)
      xp=tan(theta2)
      bgz=bgz*sqrt(con/(1+xp**2+yp**2))
      z = zn
c---check for particle being at exit edge
      if (z.ne.zloc(ne))go to 40
      b2 = el(7,ne)
      con=1.+xp**2+yp**2
      b3 = b2 - el(9,ne)
      r2=el(11,ne)
      h=1/rhoa
      if(r2.eq.0.)then
      hr2=0
      else
      hr2=1./r2
      endif
      xt=x+x*x*.5*h*tan(b2)**2-y*y*.5*h/cos(b2)**2
c      r11     t111               t133
      xpt=xp+x*h*tan(b2)+x*x*(.5*h*hr2/cos(b2)**3-h*h*.5*tan(b2)**3)
c        r22  r21            t211
      xpt=xpt-x*xp*h*tan(b2)**2-y*y*(h*h*.5*tan(b2)**3+.5*h*hr2/cos(b2)
     1 **3)
c                 t212               t233
      xpt=xpt+y*yp*h*tan(b2)**2
c                 t234
      yt=y-x*y*h*tan(b2)**2
c       r33   t313
      ypt=yp-y*h*tan(b3)-x*y*(h*hr2/cos(b2)**3-h*h/cos(b2)**2*tan(b2))
c        r44      r43         t413
      yp=ypt+x*yp*h*tan(b2)**2+xp*y*h/cos(b2)**2
c                t414             t423
      x=xt
      y=yt
      xp=xpt
      bgz=bgz*sqrt(con/(1.+yp**2+xp**2))
      iend=1
   40 continue
      bgx = xp*bgz
      bgy = yp*bgz
   49 continue
      return
   50 continue
      if(x**2+y**2.lt.el(2,ne)**2)x=-el(2,ne)
      return
      end