source: PSPA/parmelaPSPA/trunk/solnoid.f @ 488

      subroutine solnoid(dwtp,iend,ibeg,bfield)
c---first order solenoid transformation
c---------------------------------------------------------------------------
      save
c
      include 'param_sz.h'
      include 'constcom.h'
      include 'pcordcom.h'
      include 'syscom.h'
      include 'ucom.h'
c--------------------------------------------------------------------------
      ne=rne
      if(bfield.ne.0.)then
      bfd=bfield
      else
      bfd=el(4,ne)
      endif
      bz=bgz/gamma
      dz=bz*dcon*dwtp
      beta=sqrt(1.-1./gamma**2)
      cay=bfd/(brhof*sqrt(gamma**2-1.))
      bx=bgx/gamma
      by=bgy/gamma
      dx=bx*dcon*dwtp
      dy=by*dcon*dwtp
      dl=sqrt(dx**2+dy**2+dz**2)
      xp=dx/dl
      yp=dy/dl
      zp=dz/dl
c---allow particles to go backwards and bgz may go to zero need to 
c   reformulate for the general case
c      xp=bgx/bgz
c      yp=bgy/bgz
c---check for particle being at entrance of solenoid
      if (ne.ge.1) go to  20
      if(z.gt.0.)go to 30
c---particle at entrance.  apply fringe transformation.
   10 continue
      cayz=(bfd-bfld(z,sqrt(x**2+y**2),brfld))/(brhof*sqrt(gamma**2-1.))
      xp=xp+.5*cayz*y
      yp=yp-.5*cayz*x
c      zp=zp-.5*cayz*(xp*y-yp*x)
c      bgz=bgz*sqrt(con/(1.+xp**2+yp**2))
c      bz=bgz/gamma
c      fac=1./sqrt(xp**2+yp**2+zp**2)
c      xp=xp*fac
c      yp=yp*fac
c      zp=zp*fac
      if(1.-xp**2-yp**2.ge.0.)then      
      zp=sign(sqrt(1.-xp**2-yp**2),zp)
      bz=zp*beta
      dz=bz*dcon*dwtp
      else
c--- particle has turned around
      xp=xp-.5*cayz*y
      yp=yp+.5*cayz*x
      zp=-zp
      bz=zp*beta
      dz=bz*dcon*dwtp
      endif
      go to 30
   20 if((z.eq.zloc(ne-1) .and. bz.gt.0.) .or.
     * (z.eq.zloc(ne) .and. bz.lt.0.)) go to 10
   30 if(z+dz.ge.zloc(ne) .and. ne.le.nel)then
c---particle will cross into next element during this step
      iend=1
      dz=zloc(ne)-z
      dwtp=dz/(bz*dcon)
      dl=beta*dcon*dwtp
      endif
      if(ne.ge.1.and.dz.lt.0.and.z+dz.le.zloc(ne-1))then
c---particle will cross into previous element during this step
      ibeg=1
      dz=zloc(ne-1)-z
      dwtp=dz/(bz*dcon)
      dl=beta*dcon*dwtp
      endif
      cayl=cay*dl
      s=sin(cayl)
      c=cos(cayl)
      x = x + (s*xp + (1.-c)*yp)/cay
      xxp= c*xp + s*yp
      y = y + (s*yp - (1.-c)*xp)/cay
      yp = c*yp - s*xp
      xp = xxp
      z=z+dz
c---check for particle being at exit of solenoid
      if (iend.eq.0.and.ibeg.eq.0)go to 40
c---apply exit fringe transformation
      cayz=(bfld(z,sqrt(x**2+y**2),brfld)-bfd)/(brhof*sqrt(gamma**2-1.))
      xp = xp + .5*cayz*y
      yp = yp - .5*cayz*x
c      zp = zp - .5*cayz*(xp*y-yp*x)
c      fac=1./sqrt(xp**2+yp**2+zp**2)
c      xp=xp*fac
c      yp=yp*fac
c      zp=zp*fac
      fac=1.-xp**2-yp**2
      if(fac.ge.0)then
      zp=sign(sqrt(fac),zp)
      bz=zp*beta
      else
c--- particle has turned around
      xp=xp-.5*cayz*y
      yp=yp+.5*cayz*x
      zp=-zp
      bz=zp*beta
      endif
  40  bgx=xp*beta*gamma
      bgy=yp*beta*gamma
      bgz=zp*beta*gamma
      if(iend.eq.1)z=zloc(ne)
      if(ibeg.eq.1)z=zloc(ne-1)
      return
      end
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*