| 1 | subroutine drift(dwtp,iend,ibeg) |
|---|
| 2 | c---------------------------------------------------------------------- |
|---|
| 3 | c |
|---|
| 4 | include 'param_sz.h' |
|---|
| 5 | include 'bfieldcom.h' |
|---|
| 6 | include 'constcom.h' |
|---|
| 7 | include 'flagcom.h' |
|---|
| 8 | include 'misccom.h' |
|---|
| 9 | include 'pcordcom.h' |
|---|
| 10 | include 'pfieldcom.h' |
|---|
| 11 | include 'syscom.h' |
|---|
| 12 | include 'ucom.h' |
|---|
| 13 | c |
|---|
| 14 | c-------------------------------------------------------------------------- |
|---|
| 15 | c* |
|---|
| 16 | ne=rne |
|---|
| 17 | bz=bgz/gamma |
|---|
| 18 | dz=bz*dcon*dwtp |
|---|
| 19 | if(z+dz.ge.zloc(ne).and.ne.le.nel)then |
|---|
| 20 | c----particle will cross into next element during this step. |
|---|
| 21 | iend=1 |
|---|
| 22 | dummy=0. |
|---|
| 23 | if((ip.eq.999).and.(nupar.eq.1)) |
|---|
| 24 | * write(nimp,*) z,dummy,dummy,dummy |
|---|
| 25 | dz=zloc(ne)-z |
|---|
| 26 | dwtp=dz/(bz*dcon) |
|---|
| 27 | endif |
|---|
| 28 | c-----allow particle to go backwards |
|---|
| 29 | if(ne.ge.1.and.dz.lt.0.and.z+dz.le.zloc(ne-1))then |
|---|
| 30 | c----particle will cross into previous element during this step |
|---|
| 31 | ibeg=1 |
|---|
| 32 | dz=zloc(ne-1)-z |
|---|
| 33 | dwtp=dz/(bz*dcon) |
|---|
| 34 | endif |
|---|
| 35 | if(poiflag .and. z.ge.pzmin .and. z.le.pzmax .and. ne.le.nel |
|---|
| 36 | * .and. z.ge.0.) go to 10 |
|---|
| 37 | if(ifld.eq.0 .or. z+dz.lt.zmin .or. z.gt.zmax .or. ne.gt.nel |
|---|
| 38 | * .or. z.lt.0.)then |
|---|
| 39 | x=x+dz*bgx/bgz |
|---|
| 40 | y=y+dz*bgy/bgz |
|---|
| 41 | z=z+dz |
|---|
| 42 | if(iend.eq.1)z=zloc(ne) |
|---|
| 43 | if(ibeg.eq.1)z=zloc(ne-1) |
|---|
| 44 | go to 2 |
|---|
| 45 | endif |
|---|
| 46 | 10 continue |
|---|
| 47 | if(dwtp.eq.0.) go to 2 |
|---|
| 48 | beta=sqrt(1.-1./gamma**2) |
|---|
| 49 | bx=bgx/gamma |
|---|
| 50 | by=bgy/gamma |
|---|
| 51 | dx=bx*dcon*dwtp |
|---|
| 52 | dy=by*dcon*dwtp |
|---|
| 53 | dl=sqrt(dx*dx+dy*dy+dz*dz) |
|---|
| 54 | xp=dx/dl |
|---|
| 55 | yp=dy/dl |
|---|
| 56 | zp=dz/dl |
|---|
| 57 | r=sqrt(x**2+y**2) |
|---|
| 58 | temp=brhof*sqrt(gamma**2-1.) |
|---|
| 59 | if(temp.ne.0.) then |
|---|
| 60 | if(ifoclal.eq.1) then ! foclal |
|---|
| 61 | cayz=bfldlal(z)/temp |
|---|
| 62 | cayzn=bfldlal(z+dz)/temp |
|---|
| 63 | cayr=0. |
|---|
| 64 | else ! Poisson and Helmholz |
|---|
| 65 | cayz=bfld(z,r,brfld)/temp |
|---|
| 66 | cayzn=bfld(z+dz,r,brfldn)/temp |
|---|
| 67 | if(r.gt.1.e-10) then |
|---|
| 68 | cayr=-(brfld+brfldn)*dl/(r*temp) |
|---|
| 69 | else |
|---|
| 70 | cayr=0. |
|---|
| 71 | endif |
|---|
| 72 | endif |
|---|
| 73 | else |
|---|
| 74 | cayz=0. |
|---|
| 75 | cayzn=0. |
|---|
| 76 | endif |
|---|
| 77 | c----drift half way. |
|---|
| 78 | dl=dl*.5 |
|---|
| 79 | if(cayz.eq.0.)then |
|---|
| 80 | x=x+dl*xp |
|---|
| 81 | y=y+dl*yp |
|---|
| 82 | else |
|---|
| 83 | cayl=cayz*dl |
|---|
| 84 | s=sin(cayl) |
|---|
| 85 | c=cos(cayl) |
|---|
| 86 | x = x + (s*xp + (1.-c)*yp)/cayz |
|---|
| 87 | y = y + (s*yp - (1.-c)*xp)/cayz |
|---|
| 88 | xxp= c*xp + s*yp |
|---|
| 89 | yp = c*yp - s*xp |
|---|
| 90 | xp=xxp |
|---|
| 91 | endif |
|---|
| 92 | z=z+dl*zp |
|---|
| 93 | c----apply fringe transformation |
|---|
| 94 | if(ifoclal.eq.1) then |
|---|
| 95 | xp=xp+.5*(cayzn-cayz)*y*zp |
|---|
| 96 | yp=yp-.5*(cayzn-cayz)*x*zp |
|---|
| 97 | zp=zp-.5*(cayzn-cayz)*(xp*y-yp*x) |
|---|
| 98 | c------allow particle to go backwards |
|---|
| 99 | fac=1./sqrt(xp**2+yp**2+zp**2) |
|---|
| 100 | xp=xp*fac |
|---|
| 101 | yp=yp*fac |
|---|
| 102 | zp=zp*fac |
|---|
| 103 | else |
|---|
| 104 | xp=xp+.5*cayr*y |
|---|
| 105 | yp=yp-.5*cayr*x |
|---|
| 106 | fac=1.-(xp**2+yp**2) |
|---|
| 107 | if(fac.ge.0.) then |
|---|
| 108 | zp=sign(sqrt(fac),zp) |
|---|
| 109 | else |
|---|
| 110 | c---------particle has turned around |
|---|
| 111 | xp=xp-.5*cayr*y |
|---|
| 112 | yp=yp+.5*cayr*x |
|---|
| 113 | zp=-zp |
|---|
| 114 | endif |
|---|
| 115 | endif |
|---|
| 116 | bz=zp*beta |
|---|
| 117 | c----drift remaining distance. |
|---|
| 118 | if(iend.eq.1 .or. ibeg.eq.1)then |
|---|
| 119 | c----adjust dwtp |
|---|
| 120 | dwtp=.5*dwtp+dl*zp/(bz*dcon) |
|---|
| 121 | else |
|---|
| 122 | c----adjust dz and check if new dz makes particle cross into next element. |
|---|
| 123 | dz=.5*dwtp*bz*dcon |
|---|
| 124 | if(z+dz.ge.zloc(ne))then |
|---|
| 125 | c----particle will cross into next element with new dz therefor adjust dz |
|---|
| 126 | c----and dwtp. |
|---|
| 127 | iend=1 |
|---|
| 128 | dztemp=zloc(ne)-z |
|---|
| 129 | dwtp=.5*(dwtp+dztemp/dz*dwtp) |
|---|
| 130 | dl=dl*dztemp/dz |
|---|
| 131 | dx=dl*xp |
|---|
| 132 | dy=dl*yp |
|---|
| 133 | dz=dztemp |
|---|
| 134 | endif |
|---|
| 135 | if(z+dz.le.zloc(ne-1))then |
|---|
| 136 | c----particle will cross into previous element with new dz therefor adjust dz |
|---|
| 137 | c----and dwtp. |
|---|
| 138 | ibeg=1 |
|---|
| 139 | dztemp=zloc(ne-1)-z |
|---|
| 140 | dwtp=.5*(dwtp+dztemp/dz*dwtp) |
|---|
| 141 | dl=dl*dztemp/dz |
|---|
| 142 | dx=dl*xp |
|---|
| 143 | dy=dl*yp |
|---|
| 144 | dz=dztemp |
|---|
| 145 | endif |
|---|
| 146 | endif |
|---|
| 147 | if(cayzn.eq.0.)then |
|---|
| 148 | x=x+xp*dl |
|---|
| 149 | y=y+yp*dl |
|---|
| 150 | else |
|---|
| 151 | cayl=cayzn*dl |
|---|
| 152 | s=sin(cayl) |
|---|
| 153 | c=cos(cayl) |
|---|
| 154 | x = x + (s*xp + (1.-c)*yp)/cayzn |
|---|
| 155 | y = y + (s*yp - (1.-c)*xp)/cayzn |
|---|
| 156 | xxp= c*xp + s*yp |
|---|
| 157 | yp = c*yp - s*xp |
|---|
| 158 | xp = xxp |
|---|
| 159 | endif |
|---|
| 160 | bx=xp*beta |
|---|
| 161 | by=yp*beta |
|---|
| 162 | bgx=bx*gamma |
|---|
| 163 | bgy=by*gamma |
|---|
| 164 | bgz=bz*gamma |
|---|
| 165 | 1 z=z+dz |
|---|
| 166 | if(iend.eq.1)z=zloc(ne) |
|---|
| 167 | cbm 07/03/2001 if(ibeg.eq.1)z=zloc(ne-1) |
|---|
| 168 | if((ne.gt.1).and.(ibeg.eq.1))z=zloc(ne-1) |
|---|
| 169 | 2 continue |
|---|
| 170 | return |
|---|
| 171 | end |
|---|
| 172 | c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++* |
|---|
![(please configure the [header_logo] section in trac.ini)](/PSPA/chrome/site/pspa.gif)
