source: MML/trunk/at/atphysics/findelemaddiffmatrix.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: 15.5 KB
Line 
1/* findmpoleraddifmatrix.c
2
3   mex-function to calculate radiation diffusion matrix B defined in [2]
4   for multipole elements in MATLAB Accelerator Toolbox
5   A.Terebilo 8/14/00
6
7   References
8   [1] M.Sands 'The Physics of Electron Storage Rings
9   [2] Ohmi, Kirata, Oide 'From the beam-envelope matrix to synchrotron
10   radiation integrals', Phys.Rev.E  Vol.49 p.751 (1994)
11*/
12
13#include "mex.h"
14#include "atlalib.c"
15#include <math.h>
16
17
18/* Fourth order-symplectic integrator constants */
19
20#define DRIFT1    0.6756035959798286638
21#define DRIFT2   -0.1756035959798286639
22#define KICK1     1.351207191959657328
23#define KICK2    -1.702414383919314656
24
25/* Physical constants used in the calculations */
26
27#define TWOPI           6.28318530717959
28#define CGAMMA  8.846056192e-05                         /* [m]/[GeV^3] Ref[1] (4.1)      */
29#define M0C2      5.10999060e5                          /* Electron rest mass [eV]       */
30#define LAMBDABAR 3.86159323e-13                        /* Compton wavelength/2pi [m]    */
31#define CER             2.81794092e-15                  /* Classical electron radius [m] */
32#define CU        1.323094366892892                     /* 55/(24*sqrt(3)) factor        */
33
34
35
36#define SQR(X) ((X)*(X))
37
38
39
40void smpledge(double* r, double inv_rho, double angle)
41{       double psi = inv_rho*tan(angle);
42        r[1]+=r[0]*psi;
43        r[3]-=r[2]*psi;
44}
45
46
47double B2perp(double bx, double by, double irho,
48                            double x, double xpr, double y, double ypr)
49/* Calculates sqr(|e x B|) , where e is a unit vector in the direction of velocity   */
50   
51{       double v_norm2;
52        v_norm2 = 1/(SQR(1+x*irho)+ SQR(xpr) + SQR(ypr));
53
54        /* components of the  velocity vector
55           double ex, ey, ez;
56           ex = xpr;
57           ey = ypr;
58           ez = (1+x*irho);
59        */
60       
61        return((SQR(by*(1+x*irho)) + SQR(bx*(1+x*irho)) + SQR(bx*ypr - by*xpr) )*v_norm2) ;
62
63
64
65}
66 
67
68void thinkickrad(double* r, double* A, double* B, double L, double irho, double E0, int max_order)
69
70/*****************************************************************************
71Calculate and apply a multipole kick to a phase space vector *r in a multipole element.
72The reference coordinate system  may have the curvature given by the inverse
73(design) radius irho. irho = 0 for straight elements
74
75IMPORTANT !!!
76The desighn magnetic field Byo that provides this curvature By0 = irho * E0 /(c*e)
77MUST NOT be included in the dipole term PolynomB(1)(MATLAB notation)(B[0] C notation)
78of the By field expansion
79HOWEVER!!! to calculate the effect of classical radiation the full field must be
80used in the square of the |v x B|.
81When calling B2perp(Bx, By, ...), use the By = ReSum + irho, where ReSum is the
82normalized vertical field - sum of the polynomial terms in PolynomB.
83
84The kick is given by
85
86           e L      L delta      L x
87theta  = - --- B  + -------  -  -----  ,
88     x     p    y     rho           2
89            0                    rho
90
91         e L
92theta  = --- B
93     y    p   x
94           0
95
96Note: in the US convention the field is written as:
97
98                        max_order+1
99                          ----
100                          \                       n-1
101           (B + iB  ) = B rho  >   (ia  + b ) (x + iy)
102        y    x           /       n    n
103                              ----
104                          n=1
105
106Use different index notation
107
108                         max_order
109                           ----
110                           \                       n
111           (B + iB  )/ B rho  =  >   (iA  + B ) (x + iy)
112        y    x             /       n    n
113                                ----
114                           n=0
115
116A,B: i=0 ... i=max_order
117[0] - dipole, [1] - quadrupole, [2] - sextupole ...
118units for A,B[i] = 1/[m]^(i+1)
119Coeficients are stored in the PolynomA, PolynomB field of the element
120structure in MATLAB
121
122
123******************************************************************************/
124{  int i;
125        double ImSum = A[max_order];
126        double ReSum = B[max_order];   
127        double x ,xpr, y, ypr, p_norm,dp_0, B2P;
128        double ReSumTemp;
129        double CRAD = CGAMMA*E0*E0*E0/(TWOPI*1e27);
130       
131        /* recursively calculate the local transvrese magnetic field
132           Bx = ReSum, By = ImSum
133        */
134        for(i=max_order-1;i>=0;i--)
135                {       ReSumTemp = ReSum*r[0] - ImSum*r[2] + B[i];
136                        ImSum = ImSum*r[0] +  ReSum*r[2] + A[i];
137                        ReSum = ReSumTemp;
138                }
139       
140
141        /* calculate angles from momentas */   
142        p_norm = 1/(1+r[4]);
143        x   = r[0];
144        xpr = r[1]*p_norm;
145        y   = r[2];
146        ypr = r[3]*p_norm;
147
148
149        B2P = B2perp(ImSum, ReSum +irho, irho, x , xpr, y ,ypr);
150       
151        dp_0 = r[4]; /* save a copy of the initial value of dp/p */
152
153        r[4] = r[4] - CRAD*SQR(1+r[4])*B2P*(1 + x*irho + (SQR(xpr)+SQR(ypr))/2 )*L;
154
155        /* recalculate momentums from angles after losing energy to radiation   */
156        p_norm = 1/(1+r[4]);
157        r[1] = xpr/p_norm;
158        r[3] = ypr/p_norm;
159
160       
161        r[1] -=  L*(ReSum-(dp_0-r[0]*irho)*irho);
162        r[3] +=  L*ImSum;
163        r[5] +=  L*irho*r[0]; /* pathlength */
164
165
166}
167
168void thinkickM(double* orbit_in, double* A, double* B, double L,
169                                                        double irho, int max_order, double *M66)
170/* Calculate the symplectic (no radiation) transfer matrix of a
171   thin multipole kick near the entrance point orbit_in
172   For elements with straight coordinate system irho = 0
173   For curved elements the B polynomial (PolynomB in MATLAB)
174   MUST NOT include the guide field  By0 = irho * E0 /(c*e)
175
176   M is a (*double) pointer to a preallocated 1-dimentional array
177   of 36 elements of matrix M arranged column-by-column
178*/
179{  int m,n;
180   
181        double ReSumNTemp;
182   double ImSumN = max_order*A[max_order];
183   double ReSumN = max_order*B[max_order];     
184
185   /* Recursively calculate the derivatives
186           ReSumN = (irho/B0)*Re(d(By + iBx)/dx)
187           ImSumN = (irho/B0)*Im(d(By + iBx)/dy)
188    */
189        for(n=max_order-1;n>0;n--)
190                {       ReSumNTemp = (ReSumN*orbit_in[0] - ImSumN*orbit_in[2]) + n*B[n];
191                        ImSumN = ImSumN*orbit_in[0] +  ReSumN*orbit_in[2] + n*A[n];
192                        ReSumN = ReSumNTemp;
193                }
194
195        /* Initialize M66 to a 6-by-6 identity matrix */
196        for(m=0;m<36;m++)
197                M66[m]= 0;
198        /* Set diagonal elements to 1 */
199        for(m=0;m<6;m++)
200                M66[m*7] = 1;
201       
202        /* The relationship between indexes when a 6-by-6 matrix is
203           represented in MATLAB as one-dimentional array containing
204           36 elements arranged column-by-column is
205           [i][j] <---> [i+6*j]
206        */
207
208        M66[1]   = -L*ReSumN;                           /* [1][0] */
209        M66[13]  =  L*ImSumN;                           /* [1][2] */
210        M66[3]   =  L*ImSumN;                           /* [3][0] */
211        M66[15]  =  L*ReSumN;                           /* [3][2] */
212        M66[25]  =  L*irho;                                     /* [1][4] */
213        M66[1]  += -L*irho*irho;                        /* [1][0] */
214        M66[5]   =  L*irho;                                     /* [5][0] */
215
216}
217
218
219
220void thinkickB(double* orbit_in, double* A, double* B, double L,
221                                                        double irho, int max_order, double E0, double *B66)
222
223/* Calculate Ohmi's diffusion matrix of a thin multipole  element
224   For elements with straight coordinate system irho = 0
225   For curved elements the B polynomial (PolynomB in MATLAB)
226   MUST NOT include the guide field  By0 = irho * E0 /(c*e)
227   The result is stored in a preallocated 1-dimentional array B66
228   of 36 elements of matrix B arranged column-by-column
229*/
230
231{       double BB,B2P,B3P;
232        int i;
233        double p_norm = 1/(1+orbit_in[4]);
234        double p_norm2 = SQR(p_norm);
235        double ImSum = A[max_order];
236        double ReSum = B[max_order];   
237        double ReSumTemp;
238       
239        /* recursively calculate the local transvrese magnetic field
240           ReSum = irho*By/B0
241           ImSum = irho*Bx/B0
242        */
243
244        for(i=max_order-1;i>=0;i--)
245                {       ReSumTemp = ReSum*orbit_in[0] - ImSum*orbit_in[2] + B[i];
246                        ImSum = ImSum*orbit_in[0] +  ReSum*orbit_in[2] + A[i];
247                        ReSum = ReSumTemp;
248                }
249       
250       
251        /* calculate |B x n|^3 - the third power of the B field component
252           orthogonal to the normalized velocity vector n
253    */
254        B2P = B2perp(ImSum, ReSum +irho, irho, orbit_in[0] , orbit_in[1]*p_norm ,
255                                                                orbit_in[2] , orbit_in[3]*p_norm );
256        B3P = B2P*sqrt(B2P);
257
258        BB = CU * CER * LAMBDABAR *  pow(E0/M0C2,5) * L * B3P * SQR(SQR(1+orbit_in[4]))*
259                                (1+orbit_in[0]*irho + (SQR(orbit_in[1])+SQR(orbit_in[3]))*p_norm2/2);
260
261       
262        /* When a 6-by-6 matrix is represented in MATLAB as one-dimentional
263           array containing 36 elements arranged column-by-column,
264           the relationship between indexes  is
265           [i][j] <---> [i+6*j]
266
267        */
268       
269        /* initialize B66 to 0 */
270        for(i=0;i<36;i++)
271                B66[i] = 0;
272       
273        /* Populate B66 */
274        B66[7]     = BB*SQR(orbit_in[1])*p_norm2;
275        B66[19]    = BB*orbit_in[1]*orbit_in[3]*p_norm2;
276        B66[9]     = BB*B66[19];
277        B66[21]    = BB*SQR(orbit_in[3])*p_norm2;
278        B66[10] = BB*orbit_in[1]*p_norm;
279        B66[25] = BB*B66[10];
280        B66[22] = BB*orbit_in[3]*p_norm;
281        B66[27] = BB* B66[22];
282        B66[28] = BB;
283}
284
285
286
287
288
289void drift_propagateB(double *orb_in, double L,  double *B)
290{       /* Propagate cumulative Ohmi's diffusion matrix B through a drift
291           B is a (*double) pointer to 1-dimentional array
292           containing 36 elements of matrix elements arranged column-by-column
293           as in MATLAB representation
294
295           The relationship between indexes when a 6-by-6 matrix is
296           represented in MATLAB as one-dimentional array containing
297           36 elements arranged column-by-column is
298           [i][j] <---> [i+6*j]
299        */
300               
301        int m;
302               
303        double *DRIFTMAT = (double*)mxCalloc(36,sizeof(double));
304        for(m=0;m<36;m++)
305                DRIFTMAT[m] = 0;
306        /* Set diagonal elements to 1   */
307        for(m=0;m<6;m++)
308                DRIFTMAT[m*7] = 1;
309
310        DRIFTMAT[6]  =  L/(1+orb_in[4]);
311        DRIFTMAT[20] =  DRIFTMAT[6];
312        DRIFTMAT[24] = -L*orb_in[1]/SQR(1+orb_in[4]);
313        DRIFTMAT[26] = -L*orb_in[3]/SQR(1+orb_in[4]);
314        DRIFTMAT[11] =  L*orb_in[1]/SQR(1+orb_in[4]);
315        DRIFTMAT[23] =  L*orb_in[3]/SQR(1+orb_in[4]);   
316        DRIFTMAT[29] = -L*(SQR(orb_in[1])+SQR(orb_in[3]))/((1+orb_in[4])*SQR(1+orb_in[4]));
317
318        ATsandwichmmt(DRIFTMAT,B);
319        mxFree(DRIFTMAT);
320       
321}
322
323
324void edge_propagateB(double inv_rho, double angle, double *B)
325
326{       /* Propagate  Ohmi's diffusion matrix B
327           through a focusing edge  B -> E*B*E'
328            where  E is a linear map of an edge
329        */
330        int m;
331        double psi = inv_rho*tan(angle);
332       
333        for(m=0;m<6;m++)
334                {       B[1+6*m] += psi*B[6*m];
335                        B[3+6*m] -= psi*B[2+6*m];
336                }
337        for(m=0;m<6;m++)
338                {       B[m+6*1] += psi*B[m+6*0];
339                        B[m+6*3] -= psi*B[m+6*2];
340                }
341}
342
343void FindElemB(double *orbit_in, double le, double irho, double *A, double *B,
344                                        double *pt1, double* pt2,double *PR1, double *PR2,
345                                        double entrance_angle,  double exit_angle,     
346                                        int max_order, int num_int_steps,
347                                        double E0, double *BDIFF)
348
349{       /* Find Ohmi's diffusion matrix BDIFF of a thick multipole
350           BDIFF - cumulative Ohmi's diffusion is initialized to 0
351           BDIFF is preallocated 1 dimensional array to store 6-by-6 matrix
352           columnwise
353        */
354       
355        int m; 
356        double  *MKICK, *BKICK;
357
358        /* 4-th order symplectic integrator constants */
359        double SL, L1, L2, K1, K2;
360        SL = le/num_int_steps;
361        L1 = SL*DRIFT1;
362        L2 = SL*DRIFT2;
363        K1 = SL*KICK1;
364        K2 = SL*KICK2;
365       
366       
367        /* Allocate memory for thin kick matrix MKICK
368           and a diffusion matrix BKICK
369        */
370        MKICK = (double*)mxCalloc(36,sizeof(double));
371        BKICK = (double*)mxCalloc(36,sizeof(double));
372        for(m=0; m < 6; m++)
373                {       MKICK[m] = 0;
374                        BKICK[m] = 0;
375                }
376       
377        /* Transform orbit to a local coordinate system of an element */
378       
379        ATaddvv(orbit_in,pt1); 
380        ATmultmv(orbit_in,PR1);
381
382        /* This coordinate transformation does not affect
383           the cumulative diffusion matrix BDIFF
384           E*BDIFF*E' :   BDIFF stays zero     
385
386        */     
387        smpledge(orbit_in, irho, entrance_angle);       /* change in the input orbit
388                                                                                                   from edge focusing
389                                                                                                */
390       
391        edge_propagateB(irho,entrance_angle,BDIFF);             /* propagate the initial
392                                                                                                           MRAD and BDIFF through
393                                                                                                           the entrance edge
394                                                                                                        */
395
396        /* Propagate orbit_in and BDIFF through a 4-th orderintegrator */
397
398        for(m=0; m < num_int_steps; m++) /* Loop over slices    */                     
399                {               drift_propagateB(orbit_in,L1, BDIFF);
400                                ATdrift6(orbit_in,L1);
401                               
402                                thinkickM(orbit_in, A,B, K1, irho, max_order, MKICK);
403                                thinkickB(orbit_in, A,B, K1, irho, max_order, E0, BKICK);
404                                ATsandwichmmt(MKICK,BDIFF);
405                                ATaddmm(BKICK,BDIFF);
406                                thinkickrad(orbit_in, A, B, K1, irho, E0, max_order);
407               
408                                drift_propagateB(orbit_in,L2, BDIFF);
409                                ATdrift6(orbit_in,L2);
410                               
411                                thinkickM(orbit_in, A,B, K2, irho, max_order, MKICK);
412                                thinkickB(orbit_in, A,B, K2, irho, max_order, E0, BKICK);
413                                ATsandwichmmt(MKICK,BDIFF);
414                                ATaddmm(BKICK,BDIFF);
415                                thinkickrad(orbit_in, A, B, K2, irho, E0, max_order);
416       
417                                drift_propagateB(orbit_in,L2, BDIFF);
418                                ATdrift6(orbit_in,L2);
419                               
420                                thinkickM(orbit_in, A,B, K1, irho, max_order, MKICK);
421                                thinkickB(orbit_in, A,B, K1, irho, max_order, E0, BKICK);
422                                ATsandwichmmt(MKICK,BDIFF);
423                                ATaddmm(BKICK,BDIFF);
424                                thinkickrad(orbit_in, A, B,  K1, irho, E0, max_order);
425
426                                drift_propagateB(orbit_in,L1, BDIFF);
427                                ATdrift6(orbit_in,L1);
428                } 
429                smpledge(orbit_in, irho, exit_angle);
430                edge_propagateB(irho,exit_angle,BDIFF);
431                               
432                ATsandwichmmt(PR2,BDIFF);
433                                                                                               
434                mxFree(MKICK);
435                mxFree(BKICK);
436}
437
438
439void mexFunction(       int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
440/* The calling syntax of this mex-function from MATLAB is
441   FindMPoleRadDiffMatrix(ELEMENT, ORBIT)
442   ELEMENT is the element structure with field names consistent with
443           a multipole transverse field model.
444   ORBIT is a 6-by-1 vector of the closed orbit at the entrance (calculated elsewhere)
445*/
446{       int m,n;
447        double le, ba, *A, *B; 
448        double irho;
449        const mxArray * globvalptr;
450        mxArray *E0ptr;
451        double E0;              /* Design energy [eV] to be obtained from MATLAB global workspace */
452        int max_order, num_int_steps;
453        double entrance_angle, exit_angle ;
454        double *BDIFF;
455        mxArray  *mxtemp;
456
457        double *orb, *orb0;
458        double *pt1, *pt2, *PR1, *PR2;
459
460
461        m = mxGetM(prhs[1]);
462        n = mxGetN(prhs[1]);
463        if(!(m==6 && n==1))
464                mexErrMsgTxt("Second argument must be a 6-by-1 column vector");
465   
466        /* ALLOCATE memory for the output array */
467        plhs[0] = mxCreateDoubleMatrix(6,6,mxREAL);
468        BDIFF = mxGetPr(plhs[0]);
469
470
471        /* If the ELEMENT sructure does not have fields PolynomA and PolynomB
472           return zero matrix and  exit
473        */
474        if(mxGetField(prhs[0],0,"PolynomA") == NULL ||  mxGetField(prhs[0],0,"PolynomB") == NULL)
475                return;
476       
477
478        /* retrieve the value of design Energy [GeV]
479           contained in MATLAB global variable GLOBVAL.
480           GLOBVAL is a MATLAB structure
481           GLOBVAL.E0 contains the design energy of the ring [eV]
482    */
483
484        globvalptr=mexGetArrayPtr("GLOBVAL","global");
485        if(globvalptr != NULL)
486                {       E0ptr = mxGetField(globvalptr,0,"E0");
487                        if(E0ptr !=NULL)
488                                E0 = mxGetScalar(E0ptr);
489                        else
490                                mexErrMsgTxt("Global variable GLOBVAL does not have a field 'E0'");
491                }
492        else
493                mexErrMsgTxt("Global variable GLOBVAL does not exist");
494
495        orb0 = mxGetPr(prhs[1]);
496        /* make local copy of the input closed orbit vector */
497        orb = (double*)mxCalloc(6,sizeof(double));
498        for(m=0;m<6;m++)
499                orb[m] = orb0[m];
500   
501        /* Retrieve element information */
502       
503        le = mxGetScalar(mxGetField(prhs[0],0,"Length"));
504       
505        /* If ELEMENT has a zero length, return zeros matrix end exit */
506        if(le == 0)
507                return;
508       
509        A = mxGetPr(mxGetField(prhs[0],0,"PolynomA"));
510        B = mxGetPr(mxGetField(prhs[0],0,"PolynomB"));
511
512       
513
514               
515        mxtemp = mxGetField(prhs[0],0,"NumIntSteps");
516   if(mxtemp != NULL)
517                num_int_steps = (int)mxGetScalar(mxtemp);
518        else
519                mexErrMsgTxt("Field 'NumIntSteps' not found in the ELEMENT structure");
520
521        mxtemp = mxGetField(prhs[0],0,"MaxOrder");
522   if(mxtemp != NULL)
523                max_order = (int)mxGetScalar(mxtemp);
524        else
525                mexErrMsgTxt("Field 'MaxOrder' not found in the ELEMENT structure");
526
527
528        mxtemp = mxGetField(prhs[0],0,"BendingAngle");
529   if(mxtemp != NULL)
530                {       ba = mxGetScalar(mxtemp);
531                        irho = ba/le;
532                }
533        else
534                {       ba = 0;
535                        irho = 0;
536                }
537               
538        mxtemp = mxGetField(prhs[0],0,"EntranceAngle");
539        if(mxtemp != NULL)
540                entrance_angle = mxGetScalar(mxtemp);
541        else
542                        entrance_angle =0;
543
544        mxtemp = mxGetField(prhs[0],0,"ExitAngle");
545        if(mxtemp != NULL)
546                exit_angle = mxGetScalar(mxtemp);
547        else
548                        exit_angle =0;
549
550        pt1 = mxGetPr(mxGetField(prhs[0],0,"T1"));
551        pt2 = mxGetPr(mxGetField(prhs[0],0,"T2"));
552        PR1 = mxGetPr(mxGetField(prhs[0],0,"R1"));
553        PR2 = mxGetPr(mxGetField(prhs[0],0,"R2"));
554       
555
556        FindElemB(orb, le, irho, A, B,
557                                        pt1, pt2, PR1, PR2,
558                                        entrance_angle,         exit_angle,
559                                        max_order, num_int_steps, E0, BDIFF);
560}
561
562
Note: See TracBrowser for help on using the repository browser.