source: MML/trunk/at/simulator/element/user/StrCorrMPoleSymplectic4Pass.c @ 4

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

Initial import--MML version from SOLEIL@2013
File size: 12.7 KB
Line 
1#include "mex.h"
2#include<math.h>
3#include "elempass.h"
4#include "../atlalib.c"
5
6
7#define DRIFT1    0.6756035959798286638
8#define DRIFT2   -0.1756035959798286639
9#define KICK1     1.351207191959657328
10#define KICK2    -1.702414383919314656
11
12
13
14void strthinkick(double* r, double* A, double* B, double L, int max_order)
15/*****************************************************************************
16Calculate and apply a multipole kick to a 6-dimentional
17phase space vector in a straight element ( quadrupole)
18
19IMPORTANT !!!
20The reference coordinate system is straight but the field expansion may still
21contain dipole terms: PolynomA(1), PolynomB(1) - in MATLAB notation,
22A[0], B[0] - C,C++ notation
23
24
25   Note: in the US convention the transverse multipole field is written as:
26
27                         max_order+1
28                           ----
29                           \                       n-1
30           (B + iB  )/ B rho  =  >   (ia  + b ) (x + iy)
31         y    x            /       n    n
32                                ----
33                          n=1
34        is a polynomial in (x,y) with the highest order = MaxOrder
35       
36
37        Using different index notation
38   
39                         max_order
40                           ----
41                           \                       n
42           (B + iB  )/ B rho  =  >   (iA  + B ) (x + iy)
43         y    x            /       n    n
44                               ----
45                          n=0
46
47        A,B: i=0 ... max_order
48   [0] - dipole, [1] - quadrupole, [2] - sextupole ...
49   units for A,B[i] = 1/[m]^(i+1)
50        Coeficients are stroed in the PolynomA, PolynomB field of the element
51        structure in MATLAB
52
53        A[i] (C++,C) =  PolynomA(i+1) (MATLAB)
54        B[i] (C++,C) =  PolynomB(i+1) (MATLAB)
55        i = 0 .. MaxOrder
56
57******************************************************************************/
58{  int i;
59        double ReSum = B[max_order];
60        double ImSum = A[max_order];
61        double ReSumTemp;
62        for(i=max_order-1;i>=0;i--)
63        {   ReSumTemp = ReSum*r[0] - ImSum*r[2] + B[i];
64            ImSum = ImSum*r[0] +  ReSum*r[2] + A[i];
65            ReSum = ReSumTemp;
66        }
67
68    r[1] -=  L*ReSum;
69    r[3] +=  L*ImSum;
70}
71
72void fastdrift(double* r, double NormL)
73
74/*   NormL=(Physical Length)/(1+delta)  is computed externally to speed up calculations
75     in the loop if momentum deviation (delta) does not change
76     such as in 4-th order symplectic integrator w/o radiation
77*/
78
79{   double dx = NormL*r[1];
80    double dy = NormL*r[3];
81    r[0]+= dx;
82    r[2]+= dy;
83    r[5]+= NormL*(r[1]*r[1]+r[3]*r[3])/(2*(1+r[4]));
84}
85
86void StrMPoleSymplectic4Pass(double *r, double le, double *A, double *B,
87                                        int max_order, int num_int_steps,
88                                        double *T1, double *T2, 
89                                        double *R1, double *R2, int num_particles)
90{       int c,m;
91        double norm, NormL1, NormL2;   
92        double *r6;
93        bool useT1, useT2, useR1, useR2;
94        double SL, L1, L2, K1, K2;
95        SL = le/num_int_steps;
96        L1 = SL*DRIFT1;
97        L2 = SL*DRIFT2;
98        K1 = SL*KICK1;
99        K2 = SL*KICK2;
100       
101        if(T1==NULL)
102            useT1=false;
103        else 
104            useT1=true; 
105           
106    if(T2==NULL)
107            useT2=false; 
108        else 
109            useT2=true; 
110       
111        if(R1==NULL)
112            useR1=false; 
113        else 
114            useR1=true; 
115           
116    if(R2==NULL)
117            useR2=false;
118        else 
119            useR2=true;
120        for(c = 0;c<num_particles;c++)  /*Loop over particles  */
121                        {       r6 = r+c*6;     
122                            if(!mxIsNaN(r6[0]))
123                            {   
124                                        /*  misalignment at entrance  */
125                                        if(useT1)
126                                    ATaddvv(r6,T1);
127                                if(useR1)
128                                    ATmultmv(r6,R1);
129                                        /*  integrator  */
130                                        for(m=0; m < num_int_steps; m++) /*  Loop over slices */
131                                                {       r6 = r+c*6;     
132                                                        norm = 1/(1+r6[4]);
133                                                        NormL1 = L1*norm;
134                                                        NormL2 = L2*norm;
135                                                        fastdrift(r6, NormL1);
136                                                strthinkick(r6, A, B,  K1, max_order);
137                                                fastdrift(r6, NormL2);
138                                                strthinkick(r6, A, B, K2, max_order);
139                                                fastdrift(r6, NormL2);
140                                                strthinkick(r6, A, B,  K1, max_order);
141                                                fastdrift(r6, NormL1); 
142                                                } 
143                                       
144                                        /* Misalignment at exit */     
145                                if(useR2)
146                                    ATmultmv(r6,R2);
147                            if(useT2)   
148                                    ATaddvv(r6,T2);
149                            }
150                        }
151}
152
153
154
155ExportMode int* passFunction(const mxArray *ElemData, int *FieldNumbers,
156                                                                double *r_in, int num_particles, int mode)
157#define NUM_FIELDS_2_REMEMBER 10
158{       int fnum;
159        double *A , *B;
160        double  *pr1, *pr2, *pt1, *pt2, *ka;   
161        int max_order, num_int_steps;
162        double le;
163        int *returnptr;
164        int *NewFieldNumbers;
165
166        switch(mode)
167                {       case NO_LOCAL_COPY:     /* NOT used in AT1.3 Get fields by names from MATLAB workspace   */
168                                {       
169                               
170                                }       break; 
171       
172                        case MAKE_LOCAL_COPY:   /* Find field numbers first
173                                                   Save a list of field number in an array
174                                                   and make returnptr point to that array
175                                                 */
176                                {       
177                                        /* Allocate memory for integer array of
178                                           field numbers for faster futurereference
179                                        */
180                                        NewFieldNumbers = (int*)mxCalloc(NUM_FIELDS_2_REMEMBER,sizeof(int));
181
182                                        /*  Populate */
183                                       
184                                        fnum = mxGetFieldNumber(ElemData,"PolynomA");
185                                        if(fnum<0) 
186                                            mexErrMsgTxt("Required field 'PolynomA' was not found in the element data structure"); 
187                                        NewFieldNumbers[0] = fnum;
188                                        A = mxGetPr(mxGetFieldByNumber(ElemData,0,fnum));
189                                       
190                                       
191                                        fnum = mxGetFieldNumber(ElemData,"PolynomB");
192                                        if(fnum<0) 
193                                            mexErrMsgTxt("Required field 'PolynomB' was not found in the element data structure"); 
194                                        NewFieldNumbers[1] = fnum;
195                                        B = mxGetPr(mxGetFieldByNumber(ElemData,0,fnum));
196                                       
197                                       
198                                       
199                                        fnum = mxGetFieldNumber(ElemData,"MaxOrder");
200                                        if(fnum<0) 
201                                            mexErrMsgTxt("Required field 'MaxOrder' was not found in the element data structure"); 
202                                        NewFieldNumbers[2] = fnum;
203                                        max_order = (int)mxGetScalar(mxGetFieldByNumber(ElemData,0,fnum));
204                                       
205                                       
206                                        fnum = mxGetFieldNumber(ElemData,"NumIntSteps");
207                                        if(fnum<0) 
208                                            mexErrMsgTxt("Required field 'NumIntSteps' was not found in the element data structure"); 
209                                        NewFieldNumbers[3] = fnum;
210                                        num_int_steps = (int)mxGetScalar(mxGetFieldByNumber(ElemData,0,fnum));
211                                       
212                                       
213                                        fnum = mxGetFieldNumber(ElemData,"Length");
214                                        if(fnum<0) 
215                                            mexErrMsgTxt("Required field 'Length' was not found in the element data structure"); 
216                                        NewFieldNumbers[4] = fnum;
217                                        le = mxGetScalar(mxGetFieldByNumber(ElemData,0,fnum));
218                                       
219                                        fnum = mxGetFieldNumber(ElemData,"R1");
220                                        NewFieldNumbers[5] = fnum;
221                                        if(fnum<0)
222                                            pr1 = NULL;
223                                        else
224                                            pr1 = mxGetPr(mxGetFieldByNumber(ElemData,0,fnum));
225                                       
226
227                                        fnum = mxGetFieldNumber(ElemData,"R2");
228                                        NewFieldNumbers[6] = fnum;
229                                        if(fnum<0)
230                                            pr2 = NULL;
231                                        else
232                                            pr2 = mxGetPr(mxGetFieldByNumber(ElemData,0,fnum));
233                                       
234                                        fnum = mxGetFieldNumber(ElemData,"T1");
235                                            NewFieldNumbers[7] = fnum;
236                                        if(fnum<0)
237                                            pt1 = NULL;
238                                        else
239                                            pt1 = mxGetPr(mxGetFieldByNumber(ElemData,0,fnum));
240                                       
241                                        fnum = mxGetFieldNumber(ElemData,"T2");
242                                            NewFieldNumbers[8] = fnum;
243                                        if(fnum<0)
244                                            pt2 = NULL;
245                                        else
246                                            pt2 = mxGetPr(mxGetFieldByNumber(ElemData,0,fnum));
247                                                               
248                                        /* Optional: Kick angles, see section below for explanation */
249                                        /* Kicks from multipole elements can be specified as angles. This handles the
250                                           case where corrector coils are used in sextupoles and used for orbit
251                                           correction. */
252                                        fnum = mxGetFieldNumber(ElemData,"KickAngle");
253                                        NewFieldNumbers[9] = fnum;
254                                        if(fnum<0)
255                                            ka = NULL;
256                                        else
257                                            ka = mxGetPr(mxGetFieldByNumber(ElemData,0,fnum));
258                                                               
259                                        returnptr = NewFieldNumbers;
260
261                                }       break;
262
263                        case    USE_LOCAL_COPY: /* Get fields from MATLAB using field numbers
264                                                                            The second argument ponter to the array of field
265                                                                            numbers is previously created with
266                                                                            QuadLinPass( ..., MAKE_LOCAL_COPY)
267                                                                                       
268                                                                        */
269                                {       A = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[0]));
270                                        B = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[1]));
271                                        max_order = (int)mxGetScalar(mxGetFieldByNumber(ElemData,0,FieldNumbers[2]));
272                                        num_int_steps = (int)mxGetScalar(mxGetFieldByNumber(ElemData,0,FieldNumbers[3]));
273                                        le = mxGetScalar(mxGetFieldByNumber(ElemData,0,FieldNumbers[4]));
274                                       
275                                        /* Optional fields */
276                                        if(FieldNumbers[5]<0)
277                                            pr1 = NULL;
278                                        else
279                                            pr1 = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[5]));
280                                       
281                                        if(FieldNumbers[6]<0)
282                                            pr2 = NULL;
283                                        else
284                                            pr2 = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[6]));
285                                           
286                                        if(FieldNumbers[7]<0)
287                                            pt1 = NULL;
288                                        else   
289                                            pt1 = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[7]));
290                                           
291                                        if(FieldNumbers[8]<0)
292                                            pt2 = NULL;
293                                        else 
294                                            pt2 = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[8]));
295                   
296                                        if(FieldNumbers[9]<0)
297                                            ka = NULL;
298                                        else
299                                            ka = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[9]));
300
301                                        returnptr = FieldNumbers;
302                                }       break;
303                        default:
304                                {       mexErrMsgTxt("No match for calling mode in function StrMPoleSymplectic4Pass\n");
305                                }
306                }
307
308
309    if(ka!=NULL)
310    {
311        /* Positive angle must correspond to -ve B field since +ve B field corresponds
312           to a bend with the same curvature as the bend magnets.
313        */
314        B[0] -= sin(ka[0])/le;
315        A[0] += sin(ka[1])/le;
316    }
317
318    StrMPoleSymplectic4Pass(r_in, le, A, B, max_order, num_int_steps,pt1, pt2, pr1, pr2, num_particles);
319   
320    if(ka!=NULL)
321    {
322        B[0] += sin(ka[0])/le;
323        A[0] -= sin(ka[1])/le;
324    }
325   
326    return(returnptr);
327}
328
329void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
330{       int m,n;
331        double *r_in;
332        double le, *A, *B, *pr1, *pr2, *pt1, *pt2, *ka; 
333        int max_order, num_int_steps;
334        mxArray *tmpmxptr;
335
336    if(nrhs)
337        {
338        /* ALLOCATE memory for the output array of the same size as the input  */
339        m = mxGetM(prhs[1]);
340        n = mxGetN(prhs[1]);
341        if(m!=6) 
342            mexErrMsgTxt("Second argument must be a 6 x N matrix");
343       
344
345        tmpmxptr =mxGetField(prhs[0],0,"PolynomA");
346        if(tmpmxptr)
347            A = mxGetPr(tmpmxptr);
348        else
349            mexErrMsgTxt("Required field 'PolynomA' was not found in the element data structure"); 
350                                   
351        tmpmxptr =mxGetField(prhs[0],0,"PolynomB");
352        if(tmpmxptr)   
353            B = mxGetPr(tmpmxptr);
354        else
355            mexErrMsgTxt("Required field 'PolynomB' was not found in the element data structure");
356                                           
357        tmpmxptr = mxGetField(prhs[0],0,"MaxOrder");
358        if(tmpmxptr)
359            max_order = (int)mxGetScalar(tmpmxptr);
360        else
361            mexErrMsgTxt("Required field 'MaxOrder' was not found in the element data structure");
362                                       
363        tmpmxptr = mxGetField(prhs[0],0,"NumIntSteps");
364        if(tmpmxptr)   
365            num_int_steps = (int)mxGetScalar(tmpmxptr);
366        else
367            mexErrMsgTxt("Required field 'NumIntSteps' was not found in the element data structure");   
368                                   
369        tmpmxptr = mxGetField(prhs[0],0,"Length");
370        if(tmpmxptr)
371            le = mxGetScalar(tmpmxptr);
372        else
373            mexErrMsgTxt("Required field 'Length' was not found in the element data structure");   
374                                       
375                                   
376        /* Optionnal arguments */ 
377        /* Kicks from multipole elements can be specified as angles. This handles the
378        case where corrector coils are used in sextupoles and used for orbit
379        correction. */
380        tmpmxptr = mxGetField(prhs[0],0,"KickAngle");
381        if(tmpmxptr)
382            ka = mxGetPr(tmpmxptr);
383        else
384            ka = NULL;
385   
386        tmpmxptr = mxGetField(prhs[0],0,"R1");
387        if(tmpmxptr)
388            pr1 = mxGetPr(tmpmxptr);
389        else
390            pr1=NULL; 
391           
392        tmpmxptr = mxGetField(prhs[0],0,"R2");
393        if(tmpmxptr)
394            pr2 = mxGetPr(tmpmxptr);
395        else
396            pr2=NULL; 
397
398        tmpmxptr = mxGetField(prhs[0],0,"T1");
399        if(tmpmxptr)
400            pt1=mxGetPr(tmpmxptr);
401        else
402            pt1=NULL;
403
404        tmpmxptr = mxGetField(prhs[0],0,"T2");
405        if(tmpmxptr)
406            pt2=mxGetPr(tmpmxptr);
407        else
408            pt2=NULL;
409       
410
411        if(ka!=NULL)
412        {
413            /* Positive angle must correspond to -ve B field since +ve B field corresponds
414               to a bend with the same curvature as the bend magnets.
415            */
416            B[0] -= sin(ka[0])/le;
417            A[0] += sin(ka[1])/le;
418        }
419
420
421        plhs[0] = mxDuplicateArray(prhs[1]);
422        r_in = mxGetPr(plhs[0]);
423        StrMPoleSymplectic4Pass(r_in, le, A, B, max_order, num_int_steps,pt1, pt2, pr1, pr2, n);
424   
425        if(ka!=NULL)
426        {
427            B[0] += sin(ka[0])/le;
428            A[0] -= sin(ka[1])/le;
429        }
430   
431
432        }
433        else
434        {   /* return list of required fields */
435            plhs[0] = mxCreateCellMatrix(5,1);
436            mxSetCell(plhs[0],0,mxCreateString("Length"));
437            mxSetCell(plhs[0],1,mxCreateString("PolynomA"));
438            mxSetCell(plhs[0],2,mxCreateString("PolynomB"));
439            mxSetCell(plhs[0],3,mxCreateString("MaxOrder"));
440            mxSetCell(plhs[0],4,mxCreateString("NumIntSteps"));     
441           
442        if(nlhs>1) /* Required and optional fields */ 
443            {   plhs[1] = mxCreateCellMatrix(4,1);
444                mxSetCell(plhs[1],0,mxCreateString("T1"));
445                mxSetCell(plhs[1],1,mxCreateString("T2"));
446                mxSetCell(plhs[1],2,mxCreateString("R1"));
447                mxSetCell(plhs[1],3,mxCreateString("R2"));
448                mxSetCell(plhs[1],4,mxCreateString("KickAngle"));
449            }
450        }
451}
452
453
454
Note: See TracBrowser for help on using the repository browser.