source: MML/trunk/at/simulator/element/user/ThinMPoleBendPass.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: 7.0 KB
Line 
1/* ThinMPolePassBend
2   Accelerator Toolbox
3   
4   Revision 3/10/04
5   A.Terebilo terebilo@ssrl.slac.stanford.edu
6     
7   Created 07/29/2003
8   C. Steier CSteier@lbl.gov
9   
10   interprets dipole term in PolynomB as bending angle instead of
11   corrector kick (reference system changes direction).
12*/
13
14
15
16#include "mex.h"
17#include "elempass.h"
18#include "../atlalib.c"
19
20
21#define DRIFT1    0.6756035959798286638
22#define DRIFT2   -0.1756035959798286639
23#define KICK1     1.351207191959657328
24#define KICK2    -1.702414383919314656
25
26
27
28void strthinkick(double* r, double* A, double* B,  int max_order)
29/******************************************************************************
30Calculate and apply a multipole kick to a 6-dimentional
31phase space vector in a straight element ( quadrupole)
32
33IMPORTANT !!!
34The reference coordinate system is straight but the field expansion may still
35contain dipole terms: PolynomA(1), PolynomB(1) - in MATLAB notation,
36A[0], B[0] - C,C++ notation
37
38
39   Note: in the US convention the transverse multipole field is written as:
40
41                         max_order+1
42                           ----
43                           \                       n-1
44           (B + iB  )/ B rho  =  >   (ia  + b ) (x + iy)
45         y    x            /       n    n
46                                ----
47                          n=1
48        is a polynomial in (x,y) with the highest order = MaxOrder
49       
50
51        Using different index notation
52   
53                         max_order
54                           ----
55                           \                       n
56           (B + iB  )/ B rho  =  >   (iA  + B ) (x + iy)
57         y    x            /       n    n
58                               ----
59                          n=0
60
61        A,B: i=0 ... max_order
62   [0] - dipole, [1] - quadrupole, [2] - sextupole ...
63   units for A,B[i] = 1/[m]^(i+1)
64        Coeficients are stroed in the PolynomA, PolynomB field of the element
65        structure in MATLAB
66
67        A[i] (C++,C) =  PolynomA(i+1) (MATLAB)
68        B[i] (C++,C) =  PolynomB(i+1) (MATLAB)
69        i = 0 .. MaxOrder
70
71******************************************************************************/
72{  int i;
73        double ReSum = B[max_order];
74        double ImSum = A[max_order];
75        double B0temp,A0temp;
76        double ReSumTemp;
77       
78        B0temp=B[0];
79        B[0]=0.0;
80        A0temp=A[0];
81        A[0]=0.0;
82       
83        /* thin multipole kick: quadrupole and higher */
84       
85        for(i=max_order-1;i>=0;i--)
86        {   ReSumTemp = ReSum*r[0] - ImSum*r[2] + B[i];
87            ImSum = ImSum*r[0] +  ReSum*r[2] + A[i];
88            ReSum = ReSumTemp;
89        }
90
91    r[1] -=  ReSum;
92    r[3] += ImSum;
93   
94    B[0]=B0temp;
95    A[0]=A0temp;
96   
97    /* dispersion prime from energy dependend dipole kick */
98    r[1] += B0temp*r[4];
99    r[3] -= A0temp*r[4];
100   
101    /* path lengthening */
102    r[5] -= B0temp*r[0]-A0temp*r[2];
103}
104
105
106
107
108
109void ThinMPolePass(double *r, double *A, double *B, int max_order, int num_particles)
110
111
112{       int c;
113        double *r6;   
114       
115        for(c = 0;c<num_particles;c++)  /* Loop over particles */
116                        {               r6 = r+c*6;     
117                                    strthinkick(r6, A, B, max_order);
118                        }
119}
120
121
122
123ExportMode int* passFunction(const mxArray *ElemData, int *FieldNumbers,
124                                                                double *r_in, int num_particles, int mode)
125
126#define NUM_FIELDS_2_REMEMBER 3
127
128
129{       double *A , *B;
130        int max_order;
131        int *returnptr;
132        int *NewFieldNumbers,fnum;
133        mxArray *tmpmxptr;
134
135       
136        switch(mode)
137                {       case NO_LOCAL_COPY:     /* Get fields by names from MATLAB workspace  */
138                                {       tmpmxptr =mxGetField(ElemData,0,"PolynomA");
139                                    if(tmpmxptr)
140                                        A = mxGetPr(tmpmxptr);
141                                    else
142                                        mexErrMsgTxt("Required field 'PolynomA' was not found in the element data structure"); 
143                                   
144                                        tmpmxptr =mxGetField(ElemData,0,"PolynomB");
145                                    if(tmpmxptr)   
146                                            B = mxGetPr(tmpmxptr);
147                                        else
148                                            mexErrMsgTxt("Required field 'PolynomB' was not found in the element data structure");
149                                           
150                                    tmpmxptr = mxGetField(ElemData,0,"MaxOrder");
151                                    if(tmpmxptr)
152                                        max_order = (int)mxGetScalar(tmpmxptr);
153                                        else
154                                            mexErrMsgTxt("Required field 'MaxOrder' was not found in the element data structure");
155                                        returnptr = NULL;
156                               
157                                }       break; 
158       
159                        case MAKE_LOCAL_COPY:   /* Find field numbers first
160                                                                           Save a list of field number in an array
161                                                                           and make returnptr point to that array
162                                                                        */
163                                {       
164                                        /* Allocate memory for integer array of
165                                          field numbers for faster futurereference
166                                        */
167               
168                                        NewFieldNumbers = (int*)mxCalloc(NUM_FIELDS_2_REMEMBER,sizeof(int));
169
170                                        /* Populate */
171                                        fnum = mxGetFieldNumber(ElemData,"PolynomA");
172                                        if(fnum<0) 
173                                            mexErrMsgTxt("Required field 'PolynomA' was not found in the element data structure"); 
174                                        NewFieldNumbers[0] = fnum;
175                                       
176                                        fnum = mxGetFieldNumber(ElemData,"PolynomB");
177                                        if(fnum<0) 
178                                            mexErrMsgTxt("Required field 'PolynomB' was not found in the element data structure"); 
179                                        NewFieldNumbers[1] = fnum;
180                                       
181                                       
182                                       
183                                        fnum = mxGetFieldNumber(ElemData,"MaxOrder");
184                                        if(fnum<0) 
185                                            mexErrMsgTxt("Required field 'MaxOrder' was not found in the element data structure"); 
186                                        NewFieldNumbers[2] = fnum;
187                                                                               
188                                        A = mxGetPr(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[0]));
189                                        B = mxGetPr(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[1]));
190                                        max_order = (int)mxGetScalar(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[2]));
191                                                                       
192                                        returnptr = NewFieldNumbers;
193
194                                }       break;
195
196                        case    USE_LOCAL_COPY: /* Get fields from MATLAB using field numbers
197                                                                           The second argument ponter to the array of field
198                                                                           numbers is previously created with
199                                                                           QuadLinPass( ..., MAKE_LOCAL_COPY)
200                                                                        */
201                                                                                       
202                                {       A = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[0]));
203                                        B = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[1]));
204                                        max_order = (int)mxGetScalar(mxGetFieldByNumber(ElemData,0,FieldNumbers[2]));
205                               
206                                        returnptr = FieldNumbers;
207                                }       break;
208                        default:
209                                {       mexErrMsgTxt("No match for calling mode in function ThinMPolePass\n");
210                                }
211                }
212
213        ThinMPolePass(r_in,A, B, max_order,num_particles);
214       
215        return(returnptr);
216
217}
218
219 
220
221
222
223
224
225
226void mexFunction(       int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
227{       int m,n;
228        double *r_in;
229        double *A, *B ; 
230        int max_order;
231        mxArray *tmpmxptr;
232
233
234/* ALLOCATE memory for the output array of the same size as the input */
235        m = mxGetM(prhs[1]);
236        n = mxGetN(prhs[1]);
237        if(m!=6) 
238                mexErrMsgTxt("Second argument must be a 6 x N matrix");
239       
240
241        tmpmxptr =mxGetField(prhs[0],0,"PolynomA");
242        if(tmpmxptr)
243            A = mxGetPr(tmpmxptr);
244        else
245                mexErrMsgTxt("Required field 'PolynomA' was not found in the element data structure"); 
246                                   
247        tmpmxptr =mxGetField(prhs[0],0,"PolynomB");
248        if(tmpmxptr)   
249            B = mxGetPr(tmpmxptr);
250        else
251            mexErrMsgTxt("Required field 'PolynomB' was not found in the element data structure");
252                                           
253        tmpmxptr = mxGetField(prhs[0],0,"MaxOrder");
254        if(tmpmxptr)
255            max_order = (int)mxGetScalar(tmpmxptr);
256        else
257                mexErrMsgTxt("Required field 'MaxOrder' was not found in the element data structure");
258               
259    plhs[0] = mxDuplicateArray(prhs[1]);
260        r_in = mxGetPr(plhs[0]);
261
262        ThinMPolePass(r_in,A,B,max_order, n);
263       
264
265}
266
267
268
Note: See TracBrowser for help on using the repository browser.