source: PSPA/parmelaPSPA/trunk/wiggler.f @ 315

      subroutine wiggler (dwtp,iend)
c--- LANL
c---aply ideal fields for linear wiggler
c---dwtp is the time interval (phase in degrees) for the integration.
c---if particle reaches end of wiggler, iend will be set to 1
c---and dwtp will be changed to the time interval used.
c---the parameters for the wiggler are stored in array el
c---L,APER,IOUT,PERIODS,STEPS,B0,KX,KY,KW
c---------------------------------------------------------------------------
c
      include 'param_sz.h'
      include 'constcom.h'
      include 'pcordcom.h'
      include 'syscom.h'
c
      real lambda
c--------------------------------------------------------------------------
c*
      ne=rne
      iprin=el(10,ne)
      if (iprin.ne.0) then
        write (1,1)
    1   format ('     x     bgx       y     bgy       z     bgz')
        write (1,2) x,bgx,y,bgy,z,bgz
    2   format (6f8.3)
      endif
      zstart=zloc(ne-1)
c      zend=zloc(ne)
      zw=z-zstart
c---zw is the z-coordinate with respect to the wiggler
      dwt=dwtp
      wigl=el(1,ne)
      periods=el(4,ne)
      lambda=wigl/periods
      steps=el(5,ne)
c--steps is the minimum number of integration steps per wiggler period
      dzmax=lambda/steps
      betax=bgx/gamma
      betay=bgy/gamma
      betaz=bgz/gamma
      dz=betaz*dcon*dwt
c---dz is the z-distance that a particle with normalized
c---velocity betaz would travel in a time interval of dwt
      ns=1
      if (dz.gt.dzmax) then
c---need to take smaller step
         ns=dz/dzmax + .99
         dwt=dwt/ns
         dz=betaz*dcon*dwt
      endif
      dx=betax*dcon*dwt
      dy=betay*dcon*dwt
      b0=el(6,ne)
      cayx=el(7,ne)
      cayy=el(8,ne)
      cayw=el(9,ne)
      pn=sqrt(bgx**2+bgy**2+bgz**2)
c---pn is the total normalized momentum
      con1=dwt*wavel/(brhof*360.)
      tdwt=0.
c---tdwt is the total time interval used.
      iend=0
c---start loop on integration steps
      do 10 n=1,ns
c---estimate z at end of step
      z2=zw+dz
      if (z2.gt.wigl) then
c---force stop at end of wiggler
         dz=wigl-zw
         dwt=dz/(dcon*betaz)
         dx=betax*dcon*dwt
         dy=betay*dcon*dwt
         con1=dwt*wavel/(brhof*360.)
         iend=1
      endif
c---do integration by drift-impulse-drift
      z1=zw+.5*dz
      x1=x+.5*dx
      y1=y+.5*dy
c---calculate field components
      zz=amod(z1,lambda)
      ckz=cos(cayw*zz)
      skz=sin(cayw*zz)
      eky=exp(cayy*y1)
      cky=.5*(eky+1./eky)
      sky=.5*(eky-1./eky)
      if (cayx.eq.0.) then
         ckx=1.
         skx=0.
         bx=0.
      else
         ekx=exp(cayx*x1)
         ckx=.5*(ekx+1./ekx)
         skx=.5*(ekx-1./ekx)
         bx=b0*skx*sky*ckz*cayx/cayy
      endif
      by=b0*ckx*cky*ckz
      bz=-b0*ckx*sky*skz*cayw/cayy
c---apply impulse
      bgx=bgx+con1*(betay*bz-betaz*by)
      bgy=bgy+con1*(betaz*bx-betax*bz)
c     bgz=bgz+con1*(betax*by-betay*bx)
       bgz=sqrt(pn**2-bgx**2-bgy**2)
c---adjust so that total momentum is unchanged
c     sqr=sqrt(bgx**2+bgy**2+bgz**2)
c     bgx=bgx*pn/sqr
c     bgy=bgy*pn/sqr
c     bgz=bgz*pn/sqr
      betax=bgx/gamma
      betay=bgy/gamma
      betaz=bgz/gamma
      dx=betax*dcon*dwt
      dy=betay*dcon*dwt
      dz=betaz*dcon*dwt
      x=x1+.5*dx
      y=y1+.5*dy
      zw=z1+.5*dz
      z=zstart+zw
      if (iprin.ne.0) then
        write (1,2) x,bgx,y,bgy,z,bgz
      endif
      tdwt=tdwt+dwt
      if (iend.ne.0) go to 20
   10 continue
      return
c---particle has reached the end of wiggler before
c---end of time interval
   20 dwtp=tdwt
c      z=zend
      return
      end
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*