source: PSPA/parmelaPSPA/tags/v0.1/tank.f @ 447

      subroutine tank(wtp,dwtp,iend,ibeg)
c---tank transformation.
c---treat each cell as drift-gap-drift
c---------------------------------------------------------------------------
      save
c
      include 'param_sz.h'
      include 'constcom.h'
      include 'pcordcom.h'
      include 'syscom.h'
      include 'ucom.h'
c
c--------------------------------------------------------------------------
c*
      data fd/.45/
      ne=rne
      p=wtp
      dpp=dwtp
      pp=p+dpp
      bz=bgz/gamma
      dzt=bz*dcon*dpp
      zz=z-zloc(ne)+el(1,ne)
      e0t=el(4,ne)
      xnct=el(5,ne)
      xphi=el(6,ne)
      phi=el(7,ne)
      cl=el(8,ne)
      hcl=.5*cl
      xnc=aint(zz/cl)
      fz=zz-xnc*cl
c---dtc = dist to center of cell
      dtc=hcl-fz
c---dte = dist to end of cell
      dte=hcl
      if(dtc.lt.0.)dte=hcl+dtc
      if(dtc.le.0.)go to 20
c---drift to cell center, unless too far
   10 if(dtc.gt.dzt)go to 30
      dpp=dpp*dtc/dzt
      call drift(dpp,iend,ibeg)
c---update the phase
      p=p+dpp
      dpp=pp-p
c---gap transformation
      ph=amod(p+phi+xnc*xphi*180.,360.)*radian-.5*pi
      cp=cos(ph)
      sp=sin(ph)
      xp=bgx/bgz
      yp=bgy/bgz
c---drift backward to middle of first half of cell
      x=x-fd*hcl*xp
      y=y-fd*hcl*yp
c---cayrsq is (kr)**2
      cayrsq=el(10,ne)*(x**2+y**2)
      fi0=1.+cayrsq/4.*(1.+cayrsq/16.*(1.+cayrsq/36.))
      fi1okr=.5*(1.+.125*cayrsq*(1.+cayrsq/24.*(1.+cayrsq/48.)))
      g=sqrt(1.+bgz**2)
      fg=e0t*cl*fi0/erest
      dg=fg*(.5*cp+sp/pi)
      gb=g+.5*dg
      g=g+dg
      if (g.lt.1.) go to 12
      if(cayrsq.eq.0.)go to 11
      bsq=1.-1./gb**2
      bbar=sqrt(bsq)
      ff=-pi*e0t*fi1okr/(bbar*erest)
     1   *(cp*(1.+bsq)/pi+.5*sp*(1.-bsq))
      bgx=bgx+ff*x
      bgy=bgy+ff*y
      bgz=sqrt(g**2-1.)
      xp=bgx/bgz
      yp=bgy/bgz
c---drift to middle of second half of cell
   11 x=x+2.*fd*hcl*xp
      y=y+2.*fd*hcl*yp
c---cayrsq is (kr)**2
      cayrsq=el(10,ne)*(x**2+y**2)
      fi0=1.+cayrsq/4.*(1.+cayrsq/16.*(1.+cayrsq/36.))
      fi1okr=.5*(1.+.125*cayrsq*(1.+cayrsq/24.*(1.+cayrsq/48.)))
      fg=e0t*cl*fi0/erest
      dg=fg*(.5*cp-sp/pi)
      gb=g+.5*dg
      g=g+dg
c-----allow particle to go backwards
   12 if (g.lt.1.) then
      g=2-g
      bgz=-sqrt(g**2-1.)
      else
      bgz=sqrt(g**2-1.)
      endif
      if(gb.lt.1.)gb=2-gb
c      if (bgz.le.0) return
      if (cayrsq.le.0.) go to 15
      bsq=1.-1./gb**2
      bbar=sqrt(bsq)
      ff= pi*e0t*fi1okr/(bbar*erest)
     1   *(cp*(1.+bsq)/pi-.5*sp*(1.-bsq))
      bgx=bgx+ff*x
      bgy=bgy+ff*y
      xp=bgx/bgz
      yp=bgy/bgz
      x=x-.5*hcl*xp
      y=y-.5*hcl*yp
   15 gamma=sqrt(1.+bgx**2+bgy**2+bgz**2)
c---calculate new drift dist
      bz=bgz/gamma
      dzt=bz*dcon*dpp
c---drift to end of cell, unless too far
   20 if(dte.gt.dzt)go to 30
      dpp=dpp*dte/dzt
      call drift(dpp,iend,ibeg)
c---update the phase
      p=p+dpp
      dpp=pp-p
      dzt=dzt-dte
      dtc=hcl
      dte=hcl
      xnc=xnc+1.
      if(xnc.lt.xnct)go to 10
c---end of tank
      dwtp=p-wtp
      z=zloc(ne)
      return
c---drift and exit
   30 continue
      call drift(dpp,iend,ibeg)
      p=p+dpp
      dwtp=p-wtp
      return
      end
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*