source: MML/trunk/at/atphysics/findelemraddiffmatrix.m @ 5

Last change on this file since 5 was 4, checked in by zhangj, 11 years ago

Initial import--MML version from SOLEIL@2013
File size: 11.3 KB
Line 
1function [B, M, O] = findelemraddifmat(ELEM,orbit,varargin)
2%FINDELEMRADDIFMAT calculates element 'radiation diffusion matrix' B
3% [B, M, ORBITOUT] = FINDELEMRADDIFMAT(ELEM, ORBITIN);
4% Ohmi, Kirata, Oide 'From the beam-envelope matrix to synchrotron
5%    radiation integrals', Phys.Rev.E  Vol.49 p.751 (1994)
6
7
8% Fourth order-symplectic integrator constants
9
10DRIFT1   = 0.6756035959798286638
11DRIFT2   = -0.1756035959798286639
12KICK1    =  1.351207191959657328
13KICK2    = -1.702414383919314656
14
15% Physical constants used in calculations
16
17TWOPI       = 6.28318530717959
18CGAMMA      = 8.846056192e-05                   % [m]/[GeV^3] Ref[1] (4.1) 
19M0C2        = 5.10999060e5              % Electron rest mass [eV]
20LAMBDABAR   = 3.86159323e-13                    % Compton wavelength/2pi [m]
21CER         = 2.81794092e-15            % Classical electron radius [m]
22CU          = 1.323094366892892         % 55/(24*sqrt(3))
23
24
25function b2 = B2perp(B, irho, r6)
26% Calculates sqr(|e x B|) , where e is a unit vector in the direction of
27% velocity. Components of the  velocity vector:
28% ex = xpr;
29% ey = ypr;
30% ez = (1+x*irho);
31
32{       E = [r(2)/(1+r(5));r(4)/(1+r(5));1+r(1)*irho];
33    b2 = sum(cross(E/norm(E),B).^2);
34    b2 = dot(c,c);
35}
36 
37
38function rout = thinkickrad(rin, PolynomA, PolynomB, L, irho, E0, max_order)
39% Propagate particle through a thin multipole with radiation
40% Calculate field from polynomial coefficients
41P = i*PolynomA(1:max_order+1)+PolynomB(1:max_order+1);
42Z = cumprod([1, (rin(1)+i*rin(3))*ones(1,max_order)]);
43S = sum(P.*Z);
44Bx = real(S); By = imag(S);
45
46B2P = B2perp([Bx By +irho 0], irho, r);
47CRAD = CGAMMA*ELEM.Energy^3/(TWOPI*1e27);
48
49% Propagate particle
50rout = rin;
51
52% Loss of energy (dp/p) due to radiation
53rout(5) = rin(5) - CRAD*(1+rin(5))^2*B2P*...
54    (1+rin(1)*irho + (rin(1)^2+rin(3)^2)/2/(1+rin(5))^2)*L;
55
56% Change in transverse momentum due to radiation
57%   Angle does not change but dp/p changes due to radiation
58%   and therefore transverse canonical momentum changes
59%   px = x'*(1+dp/p)
60%   py = y'*(1+dp/p)
61rout(2 4]) = rin([2 4])*(1+rout(5))/(1+rin(5));
62
63% transverse kick due to magnetic field
64rout(2) = rout(2) - L*(Bx-(rin(5)-rin(1)*irho)*irho);
65rout(4) = rout(4) + L*By;
66
67% pathlength
68rout(6) = rout(6) + L*irho*rin(1);
69
70
71
72function M = thinkickM(rin, PolynomA, PolynomB, L, irho, max_order)
73%     Calculate the symplectic (no radiation) transfer matrix of a
74%    thin multipole kick near the entrance point orbit_in
75%    For elements with straight coordinate system irho = 0
76%    For curved elements the B polynomial (PolynomB in MATLAB)
77%    MUST NOT include the guide field  By0 = irho * E0 /(c*e)
78
79{   P = i*PolynomA(2:max_order+1)+PolynomB(2:max_order+1);
80    Z = cumprod([1, (rin(1)+i*rin(3))*ones(1,max_order-1)]);
81    dB = sum(P.*(1:max_order).*Z);
82
83    M = eye(6);
84
85
86
87        M(2,1)   = -L*real(dB);
88        M(2,3)   =  L*imag(dB);
89        M(4,1)   =  L*imag(dB);
90        M(4,3)   =  L*real(dB);
91        M(2,5)   =  L*irho;
92        M(2,1)   =  M(2,1) - L*irho*irho;
93        M(6,1)   =  L*irho;
94
95}
96
97
98
99function B66 = thinkickB(rin, PolynomA, PolynomB, L, irho, E0, max_order)
100%    Calculate Ohmi's diffusion matrix of a thin multipole  element
101%    For elements with straight coordinate system irho = 0
102%    For curved elements the B polynomial (PolynomB in MATLAB)
103%    MUST NOT include the guide field  By0 = irho * E0 /(c*e)
104%    The result is stored in a preallocated 1-dimentional array B66
105%    of 36 elements of matrix B arranged column-by-column
106
107P = i*PolynomA(1:max_order+1)+PolynomB(1:max_order+1);
108Z = cumprod([1, (rin(1)+i*rin(3))*ones(1,max_order)]);
109S = sum(P.*Z);
110Bx = real(S); By = imag(S);
111
112B2P = B2perp([Bx By +irho 0], irho, r);
113B3P = B2P^(3/2);
114
115p_norm = 1/(1+rin(5));
116p_norm2 = p_norm^2;
117
118BB = CU * CER * LAMBDABAR *  pow(E0/M0C2,5) * L * B3P * (1+rin(5))^4*
119                                (1+rin(1)*irho + (rin(2)^2+rin(4)^2)*p_norm2/2);
120
121
122B66 = zeros(6);
123B66(2,2)    = BB*rin(2)^2*p_norm2;
124B66(2,4)    = BB*rin(2)*rin(4)*p_norm2;
125B66(4,2)    = B66(2,4);
126B66(4,4)    = BB*rin(4)^2*p_norm2;
127B66(5,2)    = BB*rin(2)*p_norm;
128B66(2,5)    = B66(5,2);
129B66(5,4)    = BB*rin(4)*p_norm;
130B66(4,5)    = B66(5,4);
131B66(5,5)    = BB;
132
133
134function = mvoid drift_propagateB(double *orb_in, double L,  double *B)
135{       /* Propagate cumulative Ohmi's diffusion matrix B through a drift
136           B is a (*double) pointer to 1-dimentional array
137           containing 36 elements of matrix elements arranged column-by-column
138           as in MATLAB representation
139
140           The relationship between indexes when a 6-by-6 matrix is
141           represented in MATLAB as one-dimentional array containing
142           36 elements arranged column-by-column is
143           [i][j] <---> [i+6*j]
144        */
145               
146        int m;
147               
148        double *DRIFTMAT = (double*)mxCalloc(36,sizeof(double));
149        for(m=0;m<36;m++)
150                DRIFTMAT[m] = 0;
151        /* Set diagonal elements to 1   */
152        for(m=0;m<6;m++)
153                DRIFTMAT[m*7] = 1;
154
155        DRIFTMAT[6]  =  L/(1+orb_in[4]);
156        DRIFTMAT[20] =  DRIFTMAT[6];
157        DRIFTMAT[24] = -L*orb_in[1]/SQR(1+orb_in[4]);
158        DRIFTMAT[26] = -L*orb_in[3]/SQR(1+orb_in[4]);
159        DRIFTMAT[11] =  L*orb_in[1]/SQR(1+orb_in[4]);
160        DRIFTMAT[23] =  L*orb_in[3]/SQR(1+orb_in[4]);   
161        DRIFTMAT[29] = -L*(SQR(orb_in[1])+SQR(orb_in[3]))/((1+orb_in[4])*SQR(1+orb_in[4]));
162
163        ATsandwichmmt(DRIFTMAT,B);
164        mxFree(DRIFTMAT);
165       
166}
167
168
169void edge_propagateB(double inv_rho, double angle, double *B)
170
171{       /* Propagate  Ohmi's diffusion matrix B
172           through a focusing edge  B -> E*B*E'
173            where  E is a linear map of an edge
174        */
175        int m;
176        double psi = inv_rho*tan(angle);
177       
178        for(m=0;m<6;m++)
179                {       B[1+6*m] += psi*B[6*m];
180                        B[3+6*m] -= psi*B[2+6*m];
181                }
182        for(m=0;m<6;m++)
183                {       B[m+6*1] += psi*B[m+6*0];
184                        B[m+6*3] -= psi*B[m+6*2];
185                }
186}
187
188void FindElemB(double *orbit_in, double le, double irho, double *A, double *B,
189                                        double *pt1, double* pt2,double *PR1, double *PR2,
190                                        double entrance_angle,  double exit_angle,     
191                                        int max_order, int num_int_steps,
192                                        double E0, double *BDIFF)
193
194{       /* Find Ohmi's diffusion matrix BDIFF of a thick multipole
195           BDIFF - cumulative Ohmi's diffusion is initialized to 0
196           BDIFF is preallocated 1 dimensional array to store 6-by-6 matrix
197           columnwise
198        */
199       
200        int m; 
201        double  *MKICK, *BKICK;
202
203        /* 4-th order symplectic integrator constants */
204        double SL, L1, L2, K1, K2;
205        SL = le/num_int_steps;
206        L1 = SL*DRIFT1;
207        L2 = SL*DRIFT2;
208        K1 = SL*KICK1;
209        K2 = SL*KICK2;
210       
211       
212        /* Allocate memory for thin kick matrix MKICK
213           and a diffusion matrix BKICK
214        */
215        MKICK = (double*)mxCalloc(36,sizeof(double));
216        BKICK = (double*)mxCalloc(36,sizeof(double));
217        for(m=0; m < 6; m++)
218                {       MKICK[m] = 0;
219                        BKICK[m] = 0;
220                }
221       
222        /* Transform orbit to a local coordinate system of an element */
223       
224        ATaddvv(orbit_in,pt1); 
225        ATmultmv(orbit_in,PR1);
226
227        /* This coordinate transformation does not affect
228           the cumulative diffusion matrix BDIFF
229           E*BDIFF*E' :   BDIFF stays zero     
230
231        */     
232        smpledge(orbit_in, irho, entrance_angle);       /* change in the input orbit
233                                                                                                   from edge focusing
234                                                                                                */
235       
236        edge_propagateB(irho,entrance_angle,BDIFF);             /* propagate the initial
237                                                                                                           MRAD and BDIFF through
238                                                                                                           the entrance edge
239                                                                                                        */
240
241        /* Propagate orbit_in and BDIFF through a 4-th orderintegrator */
242
243        for(m=0; m < num_int_steps; m++) /* Loop over slices    */                     
244                {               drift_propagateB(orbit_in,L1, BDIFF);
245                                ATdrift6(orbit_in,L1);
246                               
247                                thinkickM(orbit_in, A,B, K1, irho, max_order, MKICK);
248                                thinkickB(orbit_in, A,B, K1, irho, max_order, E0, BKICK);
249                                ATsandwichmmt(MKICK,BDIFF);
250                                ATaddmm(BKICK,BDIFF);
251                                thinkickrad(orbit_in, A, B, K1, irho, E0, max_order);
252               
253                                drift_propagateB(orbit_in,L2, BDIFF);
254                                ATdrift6(orbit_in,L2);
255                               
256                                thinkickM(orbit_in, A,B, K2, irho, max_order, MKICK);
257                                thinkickB(orbit_in, A,B, K2, irho, max_order, E0, BKICK);
258                                ATsandwichmmt(MKICK,BDIFF);
259                                ATaddmm(BKICK,BDIFF);
260                                thinkickrad(orbit_in, A, B, K2, irho, E0, max_order);
261       
262                                drift_propagateB(orbit_in,L2, BDIFF);
263                                ATdrift6(orbit_in,L2);
264                               
265                                thinkickM(orbit_in, A,B, K1, irho, max_order, MKICK);
266                                thinkickB(orbit_in, A,B, K1, irho, max_order, E0, BKICK);
267                                ATsandwichmmt(MKICK,BDIFF);
268                                ATaddmm(BKICK,BDIFF);
269                                thinkickrad(orbit_in, A, B,  K1, irho, E0, max_order);
270
271                                drift_propagateB(orbit_in,L1, BDIFF);
272                                ATdrift6(orbit_in,L1);
273                } 
274                smpledge(orbit_in, irho, exit_angle);
275                edge_propagateB(irho,exit_angle,BDIFF);
276                               
277                ATsandwichmmt(PR2,BDIFF);
278                                                                                               
279                mxFree(MKICK);
280                mxFree(BKICK);
281}
282
283
284void mexFunction(       int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
285/* The calling syntax of this mex-function from MATLAB is
286   FindMPoleRadDiffMatrix(ELEMENT, ORBIT)
287   ELEMENT is the element structure with field names consistent with
288           a multipole transverse field model.
289   ORBIT is a 6-by-1 vector of the closed orbit at the entrance (calculated elsewhere)
290*/
291{       int m,n;
292        double le, ba, *A, *B; 
293        double irho;
294        const mxArray * globvalptr;
295        mxArray *E0ptr;
296        double E0;              /* Design energy [eV] to be obtained from MATLAB global workspace */
297        int max_order, num_int_steps;
298        double entrance_angle, exit_angle ;
299        double *BDIFF;
300        mxArray  *mxtemp;
301
302        double *orb, *orb0;
303        double *pt1, *pt2, *PR1, *PR2;
304
305
306        m = mxGetM(prhs[1]);
307        n = mxGetN(prhs[1]);
308        if(!(m==6 && n==1))
309                mexErrMsgTxt("Second argument must be a 6-by-1 column vector");
310   
311        /* ALLOCATE memory for the output array */
312        plhs[0] = mxCreateDoubleMatrix(6,6,mxREAL);
313        BDIFF = mxGetPr(plhs[0]);
314
315
316        /* If the ELEMENT sructure does not have fields PolynomA and PolynomB
317           return zero matrix and  exit
318        */
319        if(mxGetField(prhs[0],0,"PolynomA") == NULL ||  mxGetField(prhs[0],0,"PolynomB") == NULL)
320                return;
321       
322
323        /* retrieve the value of design Energy [GeV]
324           contained in MATLAB global variable GLOBVAL.
325           GLOBVAL is a MATLAB structure
326           GLOBVAL.E0 contains the design energy of the ring [eV]
327    */
328
329        globvalptr=mexGetArrayPtr("GLOBVAL","global");
330        if(globvalptr != NULL)
331                {       E0ptr = mxGetField(globvalptr,0,"E0");
332                        if(E0ptr !=NULL)
333                                E0 = mxGetScalar(E0ptr);
334                        else
335                                mexErrMsgTxt("Global variable GLOBVAL does not have a field 'E0'");
336                }
337        else
338                mexErrMsgTxt("Global variable GLOBVAL does not exist");
339
340        orb0 = mxGetPr(prhs[1]);
341        /* make local copy of the input closed orbit vector */
342        orb = (double*)mxCalloc(6,sizeof(double));
343        for(m=0;m<6;m++)
344                orb[m] = orb0[m];
345   
346        /* Retrieve element information */
347       
348        le = mxGetScalar(mxGetField(prhs[0],0,"Length"));
349       
350        /* If ELEMENT has a zero length, return zeros matrix end exit */
351        if(le == 0)
352                return;
353       
354        A = mxGetPr(mxGetField(prhs[0],0,"PolynomA"));
355        B = mxGetPr(mxGetField(prhs[0],0,"PolynomB"));
356
357       
358
359               
360        mxtemp = mxGetField(prhs[0],0,"NumIntSteps");
361   if(mxtemp != NULL)
362                num_int_steps = (int)mxGetScalar(mxtemp);
363        else
364                mexErrMsgTxt("Field 'NumIntSteps' not found in the ELEMENT structure");
365
366        mxtemp = mxGetField(prhs[0],0,"MaxOrder");
367   if(mxtemp != NULL)
368                max_order = (int)mxGetScalar(mxtemp);
369        else
370                mexErrMsgTxt("Field 'MaxOrder' not found in the ELEMENT structure");
371
372
373        mxtemp = mxGetField(prhs[0],0,"BendingAngle");
374   if(mxtemp != NULL)
375                {       ba = mxGetScalar(mxtemp);
376                        irho = ba/le;
377                }
378        else
379                {       ba = 0;
380                        irho = 0;
381                }
382               
383        mxtemp = mxGetField(prhs[0],0,"EntranceAngle");
384        if(mxtemp != NULL)
385                entrance_angle = mxGetScalar(mxtemp);
386        else
387                        entrance_angle =0;
388
389        mxtemp = mxGetField(prhs[0],0,"ExitAngle");
390        if(mxtemp != NULL)
391                exit_angle = mxGetScalar(mxtemp);
392        else
393                        exit_angle =0;
394
395        pt1 = mxGetPr(mxGetField(prhs[0],0,"T1"));
396        pt2 = mxGetPr(mxGetField(prhs[0],0,"T2"));
397        PR1 = mxGetPr(mxGetField(prhs[0],0,"R1"));
398        PR2 = mxGetPr(mxGetField(prhs[0],0,"R2"));
399       
400
401        FindElemB(orb, le, irho, A, B,
402                                        pt1, pt2, PR1, PR2,
403                                        entrance_angle,         exit_angle,
404                                        max_order, num_int_steps, E0, BDIFF);
405}
406
407
Note: See TracBrowser for help on using the repository browser.