1 | subroutine wiggler (dwtp,iend) |
---|
2 | c--- LANL |
---|
3 | c---aply ideal fields for linear wiggler |
---|
4 | c---dwtp is the time interval (phase in degrees) for the integration. |
---|
5 | c---if particle reaches end of wiggler, iend will be set to 1 |
---|
6 | c---and dwtp will be changed to the time interval used. |
---|
7 | c---the parameters for the wiggler are stored in array el |
---|
8 | c---L,APER,IOUT,PERIODS,STEPS,B0,KX,KY,KW |
---|
9 | c--------------------------------------------------------------------------- |
---|
10 | c |
---|
11 | include 'param_sz.h' |
---|
12 | include 'constcom.h' |
---|
13 | include 'pcordcom.h' |
---|
14 | include 'syscom.h' |
---|
15 | c |
---|
16 | real lambda |
---|
17 | c-------------------------------------------------------------------------- |
---|
18 | c* |
---|
19 | ne=rne |
---|
20 | iprin=el(10,ne) |
---|
21 | if (iprin.ne.0) then |
---|
22 | write (1,1) |
---|
23 | 1 format (' x bgx y bgy z bgz') |
---|
24 | write (1,2) x,bgx,y,bgy,z,bgz |
---|
25 | 2 format (6f8.3) |
---|
26 | endif |
---|
27 | zstart=zloc(ne-1) |
---|
28 | c zend=zloc(ne) |
---|
29 | zw=z-zstart |
---|
30 | c---zw is the z-coordinate with respect to the wiggler |
---|
31 | dwt=dwtp |
---|
32 | wigl=el(1,ne) |
---|
33 | periods=el(4,ne) |
---|
34 | lambda=wigl/periods |
---|
35 | steps=el(5,ne) |
---|
36 | c--steps is the minimum number of integration steps per wiggler period |
---|
37 | dzmax=lambda/steps |
---|
38 | betax=bgx/gamma |
---|
39 | betay=bgy/gamma |
---|
40 | betaz=bgz/gamma |
---|
41 | dz=betaz*dcon*dwt |
---|
42 | c---dz is the z-distance that a particle with normalized |
---|
43 | c---velocity betaz would travel in a time interval of dwt |
---|
44 | ns=1 |
---|
45 | if (dz.gt.dzmax) then |
---|
46 | c---need to take smaller step |
---|
47 | ns=dz/dzmax + .99 |
---|
48 | dwt=dwt/ns |
---|
49 | dz=betaz*dcon*dwt |
---|
50 | endif |
---|
51 | dx=betax*dcon*dwt |
---|
52 | dy=betay*dcon*dwt |
---|
53 | b0=el(6,ne) |
---|
54 | cayx=el(7,ne) |
---|
55 | cayy=el(8,ne) |
---|
56 | cayw=el(9,ne) |
---|
57 | pn=sqrt(bgx**2+bgy**2+bgz**2) |
---|
58 | c---pn is the total normalized momentum |
---|
59 | con1=dwt*wavel/(brhof*360.) |
---|
60 | tdwt=0. |
---|
61 | c---tdwt is the total time interval used. |
---|
62 | iend=0 |
---|
63 | c---start loop on integration steps |
---|
64 | do 10 n=1,ns |
---|
65 | c---estimate z at end of step |
---|
66 | z2=zw+dz |
---|
67 | if (z2.gt.wigl) then |
---|
68 | c---force stop at end of wiggler |
---|
69 | dz=wigl-zw |
---|
70 | dwt=dz/(dcon*betaz) |
---|
71 | dx=betax*dcon*dwt |
---|
72 | dy=betay*dcon*dwt |
---|
73 | con1=dwt*wavel/(brhof*360.) |
---|
74 | iend=1 |
---|
75 | endif |
---|
76 | c---do integration by drift-impulse-drift |
---|
77 | z1=zw+.5*dz |
---|
78 | x1=x+.5*dx |
---|
79 | y1=y+.5*dy |
---|
80 | c---calculate field components |
---|
81 | zz=amod(z1,lambda) |
---|
82 | ckz=cos(cayw*zz) |
---|
83 | skz=sin(cayw*zz) |
---|
84 | eky=exp(cayy*y1) |
---|
85 | cky=.5*(eky+1./eky) |
---|
86 | sky=.5*(eky-1./eky) |
---|
87 | if (cayx.eq.0.) then |
---|
88 | ckx=1. |
---|
89 | skx=0. |
---|
90 | bx=0. |
---|
91 | else |
---|
92 | ekx=exp(cayx*x1) |
---|
93 | ckx=.5*(ekx+1./ekx) |
---|
94 | skx=.5*(ekx-1./ekx) |
---|
95 | bx=b0*skx*sky*ckz*cayx/cayy |
---|
96 | endif |
---|
97 | by=b0*ckx*cky*ckz |
---|
98 | bz=-b0*ckx*sky*skz*cayw/cayy |
---|
99 | c---apply impulse |
---|
100 | bgx=bgx+con1*(betay*bz-betaz*by) |
---|
101 | bgy=bgy+con1*(betaz*bx-betax*bz) |
---|
102 | c bgz=bgz+con1*(betax*by-betay*bx) |
---|
103 | bgz=sqrt(pn**2-bgx**2-bgy**2) |
---|
104 | c---adjust so that total momentum is unchanged |
---|
105 | c sqr=sqrt(bgx**2+bgy**2+bgz**2) |
---|
106 | c bgx=bgx*pn/sqr |
---|
107 | c bgy=bgy*pn/sqr |
---|
108 | c bgz=bgz*pn/sqr |
---|
109 | betax=bgx/gamma |
---|
110 | betay=bgy/gamma |
---|
111 | betaz=bgz/gamma |
---|
112 | dx=betax*dcon*dwt |
---|
113 | dy=betay*dcon*dwt |
---|
114 | dz=betaz*dcon*dwt |
---|
115 | x=x1+.5*dx |
---|
116 | y=y1+.5*dy |
---|
117 | zw=z1+.5*dz |
---|
118 | z=zstart+zw |
---|
119 | if (iprin.ne.0) then |
---|
120 | write (1,2) x,bgx,y,bgy,z,bgz |
---|
121 | endif |
---|
122 | tdwt=tdwt+dwt |
---|
123 | if (iend.ne.0) go to 20 |
---|
124 | 10 continue |
---|
125 | return |
---|
126 | c---particle has reached the end of wiggler before |
---|
127 | c---end of time interval |
---|
128 | 20 dwtp=tdwt |
---|
129 | c z=zend |
---|
130 | return |
---|
131 | end |
---|
132 | c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++* |
---|