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++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*