1 | subroutine bend(dwtp,iend) |
---|
2 | c---second order bending magnet transformation for uniform |
---|
3 | c field magnet with curved pole faces. the edge transformations |
---|
4 | c are form pages 71 to 75 of slac report No. 75 by Karl K |
---|
5 | c brown may 1969. T216 and T436 are not used because del=0 |
---|
6 | c rho is calculated for each particle. |
---|
7 | c the transformation in the magnet is exact. |
---|
8 | c---------------------------------------------------------------------- |
---|
9 | save |
---|
10 | c |
---|
11 | include 'param_sz.h' |
---|
12 | include 'constcom.h' |
---|
13 | include 'pcordcom.h' |
---|
14 | include 'syscom.h' |
---|
15 | include 'ucom.h' |
---|
16 | c |
---|
17 | c-------------------------------------------------------------------------- |
---|
18 | c* |
---|
19 | ne=rne |
---|
20 | if(x**2+y**2.gt.el(2,ne)**2) go to 49 |
---|
21 | gr=el(4,ne)/erest |
---|
22 | bg = sqrt(gr*(2.+gr)) |
---|
23 | rho = el(1,ne)/el(5,ne) |
---|
24 | bz=bgz/gamma |
---|
25 | xp=bgx/bgz |
---|
26 | theta1=atan(xp) |
---|
27 | yp=bgy/bgz |
---|
28 | rhoa=rho*sqrt(bgz**2+bgx**2)/bg |
---|
29 | c---check for particle being at entrance edge |
---|
30 | if (ne.gt.1) go to 20 |
---|
31 | if (z.gt.0.) go to 30 |
---|
32 | 10 continue |
---|
33 | con=1.+xp**2+yp**2 |
---|
34 | b1=el(6,ne) |
---|
35 | b2=b1-el(8,ne) |
---|
36 | r1=el(10,ne) |
---|
37 | h=1/rhoa |
---|
38 | if(r1.eq.0)then |
---|
39 | hr1=0 |
---|
40 | else |
---|
41 | hr1=1./r1 |
---|
42 | endif |
---|
43 | xt=x+.5*h*(y*y/cos(b1)**2-x*x*tan(b1)**2) |
---|
44 | c r11 t133 t111 |
---|
45 | xpt=xp+x*h*tan(b1)+x*x*.5*h*hr1/cos(b1)**3+x*xp*h*tan(b1)**2 |
---|
46 | c r22 r21 t211 t212 |
---|
47 | xpt=xpt+y*y*(h*h*(.5+tan(b1)**2)*tan(b1)-.5*h*hr1/cos(b1)**3) |
---|
48 | c t233 |
---|
49 | xpt=xpt-y*yp*h*tan(b1)**2 |
---|
50 | c t234 |
---|
51 | yt=y+x*y*h*tan(b1)**2 |
---|
52 | c r33 t313 |
---|
53 | ypt=yp-h*y*tan(b2)-x*y*h*hr1/cos(b1)**3 |
---|
54 | c r44 r43 t413 |
---|
55 | yp=ypt-x*yp*h*tan(b1)**2-xp*y*h/cos(b1)**2 |
---|
56 | c t414 t423 |
---|
57 | x=xt |
---|
58 | y=yt |
---|
59 | xp=xpt |
---|
60 | bgz=bgz*sqrt(con/(1.+xp**2+yp**2)) |
---|
61 | bz=bgz/gamma |
---|
62 | bgx=xp*bgz |
---|
63 | rhoa=rho*sqrt(bgz**2+bgx**2)/bg |
---|
64 | theta1=atan(xp) |
---|
65 | go to 30 |
---|
66 | 20 continue |
---|
67 | if (z.eq.zloc(ne-1)) go to 10 |
---|
68 | 30 continue |
---|
69 | dz=bz*dcon*dwtp |
---|
70 | ang=dz/(rho+x) |
---|
71 | dz=ang*rho |
---|
72 | zn=z+dz |
---|
73 | if(zn.lt.zloc(ne))go to 35 |
---|
74 | zn=zloc(ne) |
---|
75 | dwtp=dwtp*(zn-z)/dz |
---|
76 | ang=ang*(zn-z)/dz |
---|
77 | 35 continue |
---|
78 | xx=x |
---|
79 | con=1.+xp**2+yp**2 |
---|
80 | arg=sin(theta1+ang)-(xx+rho)*sin(ang)/rhoa |
---|
81 | if(abs(arg).ge.1)go to 50 |
---|
82 | theta2=asin(arg) |
---|
83 | x=-rho+cos(ang)*(xx+rho)+rhoa*(cos(theta2)-cos(theta1+ang)) |
---|
84 | y=y+yp*rhoa*(ang+theta1-theta2)*cos(theta1) |
---|
85 | yp=yp*cos(theta1)/cos(theta2) |
---|
86 | xp=tan(theta2) |
---|
87 | bgz=bgz*sqrt(con/(1+xp**2+yp**2)) |
---|
88 | z = zn |
---|
89 | c---check for particle being at exit edge |
---|
90 | if (z.ne.zloc(ne))go to 40 |
---|
91 | b2 = el(7,ne) |
---|
92 | con=1.+xp**2+yp**2 |
---|
93 | b3 = b2 - el(9,ne) |
---|
94 | r2=el(11,ne) |
---|
95 | h=1/rhoa |
---|
96 | if(r2.eq.0.)then |
---|
97 | hr2=0 |
---|
98 | else |
---|
99 | hr2=1./r2 |
---|
100 | endif |
---|
101 | xt=x+x*x*.5*h*tan(b2)**2-y*y*.5*h/cos(b2)**2 |
---|
102 | c r11 t111 t133 |
---|
103 | xpt=xp+x*h*tan(b2)+x*x*(.5*h*hr2/cos(b2)**3-h*h*.5*tan(b2)**3) |
---|
104 | c r22 r21 t211 |
---|
105 | xpt=xpt-x*xp*h*tan(b2)**2-y*y*(h*h*.5*tan(b2)**3+.5*h*hr2/cos(b2) |
---|
106 | 1 **3) |
---|
107 | c t212 t233 |
---|
108 | xpt=xpt+y*yp*h*tan(b2)**2 |
---|
109 | c t234 |
---|
110 | yt=y-x*y*h*tan(b2)**2 |
---|
111 | c r33 t313 |
---|
112 | ypt=yp-y*h*tan(b3)-x*y*(h*hr2/cos(b2)**3-h*h/cos(b2)**2*tan(b2)) |
---|
113 | c r44 r43 t413 |
---|
114 | yp=ypt+x*yp*h*tan(b2)**2+xp*y*h/cos(b2)**2 |
---|
115 | c t414 t423 |
---|
116 | x=xt |
---|
117 | y=yt |
---|
118 | xp=xpt |
---|
119 | bgz=bgz*sqrt(con/(1.+yp**2+xp**2)) |
---|
120 | iend=1 |
---|
121 | 40 continue |
---|
122 | bgx = xp*bgz |
---|
123 | bgy = yp*bgz |
---|
124 | 49 continue |
---|
125 | return |
---|
126 | 50 continue |
---|
127 | if(x**2+y**2.lt.el(2,ne)**2)x=-el(2,ne) |
---|
128 | return |
---|
129 | end |
---|