source: MML/trunk/at/simulator/element/IdTablePass.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: 16.0 KB
Line 
1/* IdTablePass.c
2   Accelerator Toolbox
3   Created: 13/11/08
4   Z.Marti­ zeus@cells.es
5
6   Based in the matlab routine:
7        WigTablePass.m - The tracking table is described in
8        P. Elleaume, "A new approach to the electron beam dynamics in undulators
9        and wigglers", EPAC92.
10 
11*/
12
13#include <math.h>
14#include "mex.h"
15#include "elempass.h"
16#include "atlalib.c"
17#include "interpolate.c"
18#include "matrix.h"
19
20
21#define SQR(X) X*X
22
23
24double *GLOBAL_x,*GLOBAL_y,*GLOBAL_xkick1,*GLOBAL_ykick1,*GLOBAL_xkick,*GLOBAL_ykick,*GLOBAL_xkick2,*GLOBAL_ykick2; 
25int GLOBAL_m,GLOBAL_n;
26
27/*Definition of the interpolated functions*/
28double Map_x(double x,double y)
29{
30    double f;
31    /*cubic interpolation*/
32    /*splin2(GLOBAL_y,GLOBAL_x,GLOBAL_xkick,GLOBAL_xkick2,GLOBAL_n,GLOBAL_m,y,x,&f);*/
33   
34    /*biliniar interpolation*/
35    linint(GLOBAL_y,GLOBAL_x,GLOBAL_xkick,GLOBAL_m,GLOBAL_n,y,x,&f);
36    return f;
37}
38
39double Map_y(double x,double y)
40{
41    double f;
42     /*cubic interpolation*/
43    /*splin2(GLOBAL_y,GLOBAL_x,GLOBAL_ykick,GLOBAL_ykick2,GLOBAL_m,GLOBAL_n,y,x,&f);*/
44       
45    /*biliniar interpolation*/
46    linint(GLOBAL_y,GLOBAL_x,GLOBAL_ykick,GLOBAL_m,GLOBAL_n,y,x,&f);
47    return f;
48}
49
50double Map1_x(double x,double y)
51{
52    double f;
53    /*cubic interpolation*/
54    /*splin2(GLOBAL_y,GLOBAL_x,GLOBAL_xkick1,GLOBAL_xkick2,GLOBAL_n,GLOBAL_m,y,x,&f);*/
55   
56    /*biliniar interpolation*/
57    linint(GLOBAL_y,GLOBAL_x,GLOBAL_xkick1,GLOBAL_m,GLOBAL_n,y,x,&f);
58    return f;
59}
60
61double Map1_y(double x,double y)
62{
63    double f;
64     /*cubic interpolation*/
65    /*splin2(GLOBAL_y,GLOBAL_x,GLOBAL_ykick1,GLOBAL_ykick2,GLOBAL_m,GLOBAL_n,y,x,&f);*/
66       
67    /*biliniar interpolation*/
68    linint(GLOBAL_y,GLOBAL_x,GLOBAL_ykick1,GLOBAL_m,GLOBAL_n,y,x,&f);
69    return f;
70}
71
72
73void markaslost(double *r6)
74{       int i;
75   
76    r6[0] = mxGetNaN();
77       
78    for(i=1;i<6;i++)
79                r6[i] =0;
80}
81
82
83/* Set T1, T2, R1, R2 to NULL pointers to ignore misalignmets*/
84void IdKickMapModelPass(double *r, double le, double *xkick1, double *ykick1, double *xkick, double *ykick, double *x, double *y,int n,int m, int Nslice, double *T1, double *T2, double *R1, double *R2, int num_particles)
85{       double *r6,f,L1,deltaxp,deltayp,deltaxp1,deltayp1,*limitsptr;   
86        int c;
87    bool useT1, useT2, useR1, useR2;
88   
89    /*Act as AperturePass*/
90    limitsptr=(double*)mxCalloc(4,sizeof(double));
91    limitsptr[0]=x[0];
92    limitsptr[1]=x[n-1];
93    limitsptr[2]=y[0];
94    limitsptr[3]=y[m-1];
95
96    /*globalize*/
97   
98    /* For cubic interpolation only*/
99   
100    /*GLOBAL_xkick2=(double*)mxCalloc(n*m,sizeof(double));
101        GLOBAL_ykick2=(double*)mxCalloc(n*m,sizeof(double));
102    splie2(y,x,xkick,m,n,GLOBAL_xkick2);
103    splie2(y,x,ykick,m,n,GLOBAL_ykick2); */
104
105    GLOBAL_x=x;
106    GLOBAL_y=y;
107    GLOBAL_xkick1=xkick1;
108    GLOBAL_ykick1=ykick1;
109    GLOBAL_xkick=xkick;
110    GLOBAL_ykick=ykick;
111    GLOBAL_m=m; /* y used as colums*/
112    GLOBAL_n=n; /* x used as rows*/
113     
114        if(T1==NULL)
115            useT1=false;
116        else 
117            useT1=true; 
118           
119    if(T2==NULL)
120            useT2=false; 
121        else 
122            useT2=true; 
123       
124        if(R1==NULL)
125            useR1=false; 
126        else 
127            useR1=true; 
128           
129    if(R2==NULL)
130            useR2=false;
131        else 
132            useR2=true;
133             
134
135         
136    L1=le/(2*Nslice);
137        for(c = 0;c<num_particles;c++)
138                {       
139                    r6 = r+c*6;
140           
141                        if(!mxIsNaN(r6[0]) & mxIsFinite(r6[4]))
142                        /*
143                       function bend6 internally calculates the square root
144                           of the energy deviation of the particle
145                           To protect against DOMAIN and OVERFLOW error, check if the
146                           fifth component of the phase spacevector r6[4] is finite
147                        */
148                        {                       
149                if(r6[0]<limitsptr[0] || r6[0]>limitsptr[1] || r6[2]<limitsptr[2] || r6[2]>limitsptr[3])
150                {
151                    markaslost(r6);
152                }
153                else
154                {
155                    /* Misalignment at entrance */
156                    if(useT1) ATaddvv(r6,T1);
157                    if(useR1) ATmultmv(r6,R1);
158                    /*Tracking in the main body*/
159                    for(m=0; m < Nslice; m++) /* Loop over slices*/                     
160                    {             
161                        r6 = r+c*6;               
162                        ATdrift6(r6,L1);       
163                        if (!mxIsNaN(r6[0])&&!mxIsNaN(r6[2])) 
164                        {     /*The kick from IDs varies quadratically, not linearly, with energy.   */
165                            deltaxp = (1.0/Nslice)*Map_x(r6[0],r6[2])/(1.0+r6[4]);       
166                            deltayp = (1.0/Nslice)*Map_y(r6[0],r6[2])/(1.0+r6[4]); 
167                            deltaxp1 = (1.0/Nslice)*Map1_x(r6[0],r6[2]);       
168                            deltayp1 = (1.0/Nslice)*Map1_y(r6[0],r6[2]); 
169                            r6[1] = r6[1] + deltaxp + deltaxp1; 
170                            r6[3] = r6[3] + deltayp + deltayp1;
171                        }
172                        ATdrift6(r6,L1);       
173                    } 
174                    /* Misalignment at exit */ 
175                    if(useR2) ATmultmv(r6,R2);
176                    if(useT2) ATaddvv(r6,T2);                   
177                }
178                }
179                }       
180}
181
182
183
184/*Just for debugg!! Subtitute routine for a drift*/
185void DriftPass(double *r_in, double le, int num_particles)
186/* le - physical length
187   r_in - 6-by-N matrix of initial conditions reshaped into
188   1-d array of 6*N elements
189*/
190{       int c, c6;
191        double p_norm, NormL;
192    mexPrintf("Length: %f\n",le);
193    mexPrintf("Initial %e %e %e %e %e %e\n",r_in[0],r_in[1],r_in[2],r_in[3],r_in[4],r_in[5]);
194        for(c = 0;c<num_particles;c++)
195                {       c6 = c*6;
196                    if(!mxIsNaN(r_in[c6]))
197                        {    p_norm = 1/(1+r_in[c6+4]); 
198                            NormL  = le*p_norm;
199                            r_in[c6+0]+= NormL*r_in[c6+1];
200                            r_in[c6+2]+= NormL*r_in[c6+3];
201                            r_in[c6+5]+= NormL*p_norm*(r_in[c6+1]*r_in[c6+1]+r_in[c6+3]*r_in[c6+3])/2;
202                         }
203                       
204                }
205    mexPrintf("Final %e %e %e %e %e %e\n",r_in[0],r_in[1],r_in[2],r_in[3],r_in[4],r_in[5]);
206}
207
208
209
210
211
212
213ExportMode int* passFunction(const mxArray *ElemData,int *FieldNumbers,
214                                                                double *r_in,int num_particles,int mode)
215
216#define NUM_FIELDS_2_REMEMBER 12
217
218{       int *returnptr,n,m;
219        int *NewFieldNumbers,fnum,Nslice;
220        double le;
221        double  *pr1, *pr2, *pt1, *pt2, *xkick, *ykick, *xkick1, *ykick1, *x, *y; 
222
223       
224        switch(mode)
225                {       case NO_LOCAL_COPY:     /* Get fields by names from MATLAB workspace   */
226                                {       /* NO_LOCAL_COPY - obsolete since AT1.3 */
227                                    returnptr = NULL;
228                                }       break; 
229
230                        case MAKE_LOCAL_COPY:   /* Find field numbers first
231                                                                           Save a list of field number in an array
232                                                                           and make returnptr point to that array
233                                                                        */
234                                {       
235                                        /* Populate */                                 
236                                        /*mexPrintf("Make copy\n");*/
237                    NewFieldNumbers = (int*)mxCalloc(NUM_FIELDS_2_REMEMBER,sizeof(int));
238                                       
239                                       
240                                       
241                                        fnum = mxGetFieldNumber(ElemData,"Length");
242                                        if(fnum<0) 
243                                            mexErrMsgTxt("Required field 'Length' was not found in the element data structure"); 
244                                        NewFieldNumbers[0] = fnum;
245                                        le = mxGetScalar(mxGetFieldByNumber(ElemData,0,fnum));
246
247                                        fnum = mxGetFieldNumber(ElemData,"xkick1");
248                                        if(fnum<0) 
249                                            mexErrMsgTxt("Required field 'xkick1' was not found in the element data structure");                                       
250                                        NewFieldNumbers[1] = fnum;
251                    xkick1 = mxGetPr(mxGetFieldByNumber(ElemData,0,fnum));
252                    n = mxGetN(mxGetFieldByNumber(ElemData,0,fnum));
253                    m = mxGetM(mxGetFieldByNumber(ElemData,0,fnum));
254                   
255                                        fnum = mxGetFieldNumber(ElemData,"ykick1");
256                                        if(fnum<0) 
257                                            mexErrMsgTxt("Required field 'ykick1' was not found in the element data structure");                                       
258                                        NewFieldNumbers[2] = fnum;
259                     ykick1 = mxGetPr(mxGetFieldByNumber(ElemData,0,fnum));
260                     
261                   
262                                        fnum = mxGetFieldNumber(ElemData,"xkick");
263                                        if(fnum<0) 
264                                            mexErrMsgTxt("Required field 'xkick' was not found in the element data structure");                                         
265                                        NewFieldNumbers[3] = fnum;
266                    xkick = mxGetPr(mxGetFieldByNumber(ElemData,0,fnum));
267                   
268                                        fnum = mxGetFieldNumber(ElemData,"ykick");
269                                        if(fnum<0) 
270                                            mexErrMsgTxt("Required field 'ykick' was not found in the element data structure");                                         
271                                        NewFieldNumbers[4] = fnum;
272                     ykick = mxGetPr(mxGetFieldByNumber(ElemData,0,fnum));
273                   
274                                        fnum = mxGetFieldNumber(ElemData,"xtable");
275                                        if(fnum<0) 
276                                            mexErrMsgTxt("Required field 'x' was not found in the element data structure");                                     
277                                        NewFieldNumbers[5] = fnum;
278                    x = mxGetPr(mxGetFieldByNumber(ElemData,0,fnum));
279                   
280                                        fnum = mxGetFieldNumber(ElemData,"ytable");
281                                        if(fnum<0) 
282                                            mexErrMsgTxt("Required field 'y' was not found in the element data structure");                                     
283                                        NewFieldNumbers[6] = fnum;
284                    y = mxGetPr(mxGetFieldByNumber(ElemData,0,fnum));
285                   
286                   
287                   
288                                                       
289                   
290                   
291                   
292                   
293                   
294                    fnum = mxGetFieldNumber(ElemData,"Nslice");
295                                        if(fnum<0) 
296                                            mexErrMsgTxt("Required field 'Nslice' was not found in the element data structure");                                       
297                                        NewFieldNumbers[7] = fnum;
298                    Nslice = (int)mxGetScalar(mxGetFieldByNumber(ElemData,0,fnum));
299                   
300                                       
301                                        /* Optional fields */
302
303                               
304                    fnum = mxGetFieldNumber(ElemData,"R1");
305                                        NewFieldNumbers[8] = fnum;
306                                        if(fnum<0)
307                                            pr1 = NULL;
308                                        else
309                                            pr1 = mxGetPr(mxGetFieldByNumber(ElemData,0,fnum));
310                                       
311
312                    fnum = mxGetFieldNumber(ElemData,"R2");
313                                        NewFieldNumbers[9] = fnum;
314                                        if(fnum<0)
315                                            pr2 = NULL;
316                                        else
317                                            pr2 = mxGetPr(mxGetFieldByNumber(ElemData,0,fnum));
318                                       
319                    fnum = mxGetFieldNumber(ElemData,"T1");
320                        NewFieldNumbers[10] = fnum;
321                                        if(fnum<0)
322                                            pt1 = NULL;
323                                        else
324                                            pt1 = mxGetPr(mxGetFieldByNumber(ElemData,0,fnum));
325                                       
326                        fnum = mxGetFieldNumber(ElemData,"T2");
327                        NewFieldNumbers[11] = fnum;
328                                        if(fnum<0)
329                                            pt2 = NULL;
330                                        else
331                                            pt2 = mxGetPr(mxGetFieldByNumber(ElemData,0,fnum));
332
333                                        returnptr = NewFieldNumbers;
334
335                                }       break;
336
337                        case    USE_LOCAL_COPY: /* Get fields from MATLAB using field numbers
338                                                                            The second argument ponter to the array of field
339                                                                                numbers is previously created with
340                                                                                BendLinearPass( ..., MAKE_LOCAL_COPY)
341                                                                        */
342                                                                                       
343                                {       
344                    /*mexPrintf("Aqui hi vas?\n");*/
345                    le = mxGetScalar(mxGetFieldByNumber(ElemData,0,FieldNumbers[0]));
346                    xkick1 = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[1]));
347                    ykick1 = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[2]));
348                    xkick = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[3]));
349                    ykick = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[4]));
350                    x = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[5]));
351                    y = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[6]));
352                    Nslice = (int)mxGetScalar(mxGetFieldByNumber(ElemData,0,FieldNumbers[7]));       
353                    n = mxGetN(mxGetFieldByNumber(ElemData,0,FieldNumbers[1]));
354                    m = mxGetM(mxGetFieldByNumber(ElemData,0,FieldNumbers[1]));
355                   
356                                        /* Optional fields */
357                                        if(FieldNumbers[8]<0)
358                                            pr1 = NULL;
359                                        else
360                                            pr1 = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[8]));
361                                       
362                                        if(FieldNumbers[9]<0)
363                                            pr2 = NULL;
364                                        else
365                                            pr2 = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[9]));
366                                       
367                                           
368                                        if(FieldNumbers[10]<0)
369                                            pt1 = NULL;
370                                        else   
371                                            pt1 = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[10]));
372                                           
373                                        if(FieldNumbers[11]<0)
374                                            pt2 = NULL;
375                                        else 
376                                            pt2 = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[11]));
377
378
379
380                                        returnptr = FieldNumbers;
381
382                                }       break;
383
384                }     
385   
386   
387    /*DriftPass(r_in,le,num_particles);*/
388    IdKickMapModelPass(r_in, le,xkick1,ykick1,xkick,ykick,x,y,n,m,Nslice, 
389                                                        pt1, pt2, pr1, pr2, num_particles);
390       
391    return(returnptr);
392}
393
394
395
396
397
398
399
400void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
401{       int m,n,m_map,n_map,Nslice;
402        double *r_in,le,*xkick,*ykick,*xkick1,*ykick1,*x,*y; 
403        double  *pr1, *pr2, *pt1, *pt2;
404        mxArray *tmpmxptr;
405 
406        if(nrhs)
407        {
408        /* ALLOCATE memory for the output array of the same size as the input */
409            m = mxGetM(prhs[1]);
410            n = mxGetN(prhs[1]);
411            if(m!=6) 
412            {mexErrMsgTxt("Second argument must be a 6 x N matrix");}   
413           
414            tmpmxptr = mxGetField(prhs[0],0,"Length");
415            if(tmpmxptr)
416                le = mxGetScalar(tmpmxptr);
417            else
418                mexErrMsgTxt("Required field 'Length' was not found in the element data structure"); 
419           
420        tmpmxptr = mxGetField(prhs[0],0,"xkick1");
421            if(tmpmxptr)
422        {
423            xkick1 = mxGetPr(tmpmxptr);
424            n_map = mxGetN(tmpmxptr);       
425            m_map = mxGetM(tmpmxptr);
426        }
427            else
428                mexErrMsgTxt("Required field 'xkick1' was not found in the element data structure");               
429       
430            tmpmxptr = mxGetField(prhs[0],0,"ykick1");
431            if(tmpmxptr)
432                ykick1 = mxGetPr(tmpmxptr);
433            else
434                mexErrMsgTxt("Required field 'ykick1' was not found in the element data structure");   
435       
436       
437            tmpmxptr = mxGetField(prhs[0],0,"xkick");
438            if(tmpmxptr)
439        {
440            xkick = mxGetPr(tmpmxptr);
441        }
442            else
443                mexErrMsgTxt("Required field 'xkick' was not found in the element data structure");         
444       
445            tmpmxptr = mxGetField(prhs[0],0,"ykick");
446            if(tmpmxptr)
447                ykick = mxGetPr(tmpmxptr);
448            else
449                mexErrMsgTxt("Required field 'ykick' was not found in the element data structure");           
450       
451            tmpmxptr = mxGetField(prhs[0],0,"xtable");
452            if(tmpmxptr)
453                x = mxGetPr(tmpmxptr);
454            else
455                mexErrMsgTxt("Required field 'x' was not found in the element data structure");           
456       
457            tmpmxptr = mxGetField(prhs[0],0,"ytable");
458            if(tmpmxptr)
459                y = mxGetPr(tmpmxptr);
460            else
461                mexErrMsgTxt("Required field 'y' was not found in the element data structure");   
462       
463            tmpmxptr = mxGetField(prhs[0],0,"Nslice");
464            if(tmpmxptr)
465                Nslice = (int)mxGetScalar(tmpmxptr);
466            else
467                mexErrMsgTxt("Required field 'Nslice' was not found in the element data structure");       
468       
469       
470       
471        /*Optional fields*/ 
472       
473       
474            tmpmxptr = mxGetField(prhs[0],0,"R1");
475            if(tmpmxptr)
476                pr1 = mxGetPr(tmpmxptr);
477            else
478                pr1=NULL; 
479           
480            tmpmxptr = mxGetField(prhs[0],0,"R2");
481            if(tmpmxptr)
482                pr2 = mxGetPr(tmpmxptr);
483            else
484                pr2=NULL; 
485           
486            tmpmxptr = mxGetField(prhs[0],0,"T1");
487            if(tmpmxptr)
488                pt1=mxGetPr(tmpmxptr);
489            else
490                pt1=NULL;
491           
492            tmpmxptr = mxGetField(prhs[0],0,"T2");
493            if(tmpmxptr)
494                pt2=mxGetPr(tmpmxptr);
495            else
496                pt2=NULL; 
497       
498
499
500        plhs[0] = mxDuplicateArray(prhs[1]);
501            r_in = mxGetPr(plhs[0]);
502
503                IdKickMapModelPass(r_in, le, xkick1,ykick1, xkick,ykick,x,y,n_map,m_map,Nslice, 
504                                                        pt1, pt2, pr1, pr2, n);
505        }
506        else                             
507        {   /* return list of required fields */
508            plhs[0] = mxCreateCellMatrix(8,1);
509            mxSetCell(plhs[0],0,mxCreateString("Length"));
510            mxSetCell(plhs[0],1,mxCreateString("xkick1"));
511            mxSetCell(plhs[0],2,mxCreateString("ykick1"));
512            mxSetCell(plhs[0],3,mxCreateString("xkick"));
513            mxSetCell(plhs[0],4,mxCreateString("ykick"));
514            mxSetCell(plhs[0],5,mxCreateString("x"));
515        mxSetCell(plhs[0],6,mxCreateString("y"));
516        mxSetCell(plhs[0],7,mxCreateString("Nslice"));
517                       
518            if(nlhs>1) /* Required and optional fields */ 
519            {   plhs[1] = mxCreateCellMatrix(4,1);
520                mxSetCell(plhs[1],0,mxCreateString("T1"));
521                mxSetCell(plhs[1],1,mxCreateString("T2"));
522                mxSetCell(plhs[1],2,mxCreateString("R1"));
523                mxSetCell(plhs[1],3,mxCreateString("R2"));
524            }
525        }
526 }
Note: See TracBrowser for help on using the repository browser.