1 | subroutine solnoid(dwtp,iend,ibeg,bfield) |
---|
2 | c---first order solenoid transformation |
---|
3 | c--------------------------------------------------------------------------- |
---|
4 | save |
---|
5 | c |
---|
6 | include 'param_sz.h' |
---|
7 | include 'constcom.h' |
---|
8 | include 'pcordcom.h' |
---|
9 | include 'syscom.h' |
---|
10 | include 'ucom.h' |
---|
11 | c-------------------------------------------------------------------------- |
---|
12 | ne=rne |
---|
13 | if(bfield.ne.0.)then |
---|
14 | bfd=bfield |
---|
15 | else |
---|
16 | bfd=el(4,ne) |
---|
17 | endif |
---|
18 | bz=bgz/gamma |
---|
19 | dz=bz*dcon*dwtp |
---|
20 | beta=sqrt(1.-1./gamma**2) |
---|
21 | cay=bfd/(brhof*sqrt(gamma**2-1.)) |
---|
22 | bx=bgx/gamma |
---|
23 | by=bgy/gamma |
---|
24 | dx=bx*dcon*dwtp |
---|
25 | dy=by*dcon*dwtp |
---|
26 | dl=sqrt(dx**2+dy**2+dz**2) |
---|
27 | xp=dx/dl |
---|
28 | yp=dy/dl |
---|
29 | zp=dz/dl |
---|
30 | c---allow particles to go backwards and bgz may go to zero need to |
---|
31 | c reformulate for the general case |
---|
32 | c xp=bgx/bgz |
---|
33 | c yp=bgy/bgz |
---|
34 | c---check for particle being at entrance of solenoid |
---|
35 | if (ne.ge.1) go to 20 |
---|
36 | if(z.gt.0.)go to 30 |
---|
37 | c---particle at entrance. apply fringe transformation. |
---|
38 | 10 continue |
---|
39 | cayz=(bfd-bfld(z,sqrt(x**2+y**2),brfld))/(brhof*sqrt(gamma**2-1.)) |
---|
40 | xp=xp+.5*cayz*y |
---|
41 | yp=yp-.5*cayz*x |
---|
42 | c zp=zp-.5*cayz*(xp*y-yp*x) |
---|
43 | c bgz=bgz*sqrt(con/(1.+xp**2+yp**2)) |
---|
44 | c bz=bgz/gamma |
---|
45 | c fac=1./sqrt(xp**2+yp**2+zp**2) |
---|
46 | c xp=xp*fac |
---|
47 | c yp=yp*fac |
---|
48 | c zp=zp*fac |
---|
49 | if(1.-xp**2-yp**2.ge.0.)then |
---|
50 | zp=sign(sqrt(1.-xp**2-yp**2),zp) |
---|
51 | bz=zp*beta |
---|
52 | dz=bz*dcon*dwtp |
---|
53 | else |
---|
54 | c--- particle has turned around |
---|
55 | xp=xp-.5*cayz*y |
---|
56 | yp=yp+.5*cayz*x |
---|
57 | zp=-zp |
---|
58 | bz=zp*beta |
---|
59 | dz=bz*dcon*dwtp |
---|
60 | endif |
---|
61 | go to 30 |
---|
62 | 20 if((z.eq.zloc(ne-1) .and. bz.gt.0.) .or. |
---|
63 | * (z.eq.zloc(ne) .and. bz.lt.0.)) go to 10 |
---|
64 | 30 if(z+dz.ge.zloc(ne) .and. ne.le.nel)then |
---|
65 | c---particle will cross into next element during this step |
---|
66 | iend=1 |
---|
67 | dz=zloc(ne)-z |
---|
68 | dwtp=dz/(bz*dcon) |
---|
69 | dl=beta*dcon*dwtp |
---|
70 | endif |
---|
71 | if(ne.ge.1.and.dz.lt.0.and.z+dz.le.zloc(ne-1))then |
---|
72 | c---particle will cross into previous element during this step |
---|
73 | ibeg=1 |
---|
74 | dz=zloc(ne-1)-z |
---|
75 | dwtp=dz/(bz*dcon) |
---|
76 | dl=beta*dcon*dwtp |
---|
77 | endif |
---|
78 | cayl=cay*dl |
---|
79 | s=sin(cayl) |
---|
80 | c=cos(cayl) |
---|
81 | x = x + (s*xp + (1.-c)*yp)/cay |
---|
82 | xxp= c*xp + s*yp |
---|
83 | y = y + (s*yp - (1.-c)*xp)/cay |
---|
84 | yp = c*yp - s*xp |
---|
85 | xp = xxp |
---|
86 | z=z+dz |
---|
87 | c---check for particle being at exit of solenoid |
---|
88 | if (iend.eq.0.and.ibeg.eq.0)go to 40 |
---|
89 | c---apply exit fringe transformation |
---|
90 | cayz=(bfld(z,sqrt(x**2+y**2),brfld)-bfd)/(brhof*sqrt(gamma**2-1.)) |
---|
91 | xp = xp + .5*cayz*y |
---|
92 | yp = yp - .5*cayz*x |
---|
93 | c zp = zp - .5*cayz*(xp*y-yp*x) |
---|
94 | c fac=1./sqrt(xp**2+yp**2+zp**2) |
---|
95 | c xp=xp*fac |
---|
96 | c yp=yp*fac |
---|
97 | c zp=zp*fac |
---|
98 | fac=1.-xp**2-yp**2 |
---|
99 | if(fac.ge.0)then |
---|
100 | zp=sign(sqrt(fac),zp) |
---|
101 | bz=zp*beta |
---|
102 | else |
---|
103 | c--- particle has turned around |
---|
104 | xp=xp-.5*cayz*y |
---|
105 | yp=yp+.5*cayz*x |
---|
106 | zp=-zp |
---|
107 | bz=zp*beta |
---|
108 | endif |
---|
109 | 40 bgx=xp*beta*gamma |
---|
110 | bgy=yp*beta*gamma |
---|
111 | bgz=zp*beta*gamma |
---|
112 | if(iend.eq.1)z=zloc(ne) |
---|
113 | if(ibeg.eq.1)z=zloc(ne-1) |
---|
114 | return |
---|
115 | end |
---|
116 | c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++* |
---|