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

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

Initial import--MML version from SOLEIL@2013
File size: 16.7 KB
Line 
1/* TaylormapPass.c
2   Accelerator Toolbox
3   Created on 09/26/2003
4   W. Wan wwan@lbl.gov */
5
6#include "mex.h"
7#include<stdio.h>
8#include<math.h>
9#include "../atlalib.c"
10#include "elempass.h"
11
12#define PI        3.141592653589793
13
14void TaylormapPass(double *r_in, double le,
15                        double *x, int nx, double *px, int npx,
16                        double *y, int ny, double *py, int npy,
17                        double *delta, int ndelta, double *t, int nt,
18                        double *T1, double *T2,
19                        double *R1, double *R2, int num_particles)
20/*
21   r_in -- 6-by-N matrix of initial conditions reshaped into
22           1-d array of 6*N elements
23   le -- physical length
24   x -- Taylor map of final x
25   nx -- number of monomials in the map x
26   px -- Taylor map of final px
27   npx -- number of monomials in the map px
28   y -- Taylor map of final y
29   ny -- number of monomials in the map y
30   py -- Taylor map of final py
31   npy -- number of monomials in the map py
32   delta -- Taylor map of final delta p / p
33   ndelta -- number of monomials in the map delta
34   t -- Taylor map of final t
35   nt -- number of monomials in the map t
36   T1 -- 6 x 1 translation at entrance
37   T2 -- 6 x 1 translation at exit
38   R1 -- 6 x 6 rotation matrix at the entrance
39   R2 -- 6 x 6 rotation matrix at the exit
40   num_particles -- number of particles
41*/
42#define NUM_COLUMNS 9
43{
44 int c,i,j;
45 double *r6; 
46 double *rtmp;
47 double tmp;
48
49 rtmp = (double*)mxCalloc(6,sizeof(double));
50
51 for(c = 0;c<num_particles;c++) {
52  r6 = r_in+c*6;
53  for(i=0;i<6;i++) {
54   rtmp[i] = r6[i];
55   r6[i] = 0;
56  }
57
58  /* linearized misalignment transformation at the entrance */
59  ATaddvv(rtmp,T1);
60  ATmultmv(rtmp,R1);
61
62  /* Taylor map */
63
64  /* x final */
65  for(i=0;i<nx;i++) {
66   tmp = x[i+1*nx];
67   for(j=0;j<6;j++) {
68    tmp *= pow(rtmp[j], (int)x[i+(j+3)*nx]);
69   }
70   r6[0] += tmp;
71  } 
72  /* px final */
73  for(i=0;i<npx;i++) {
74   tmp = px[i+1*npx];
75   for(j=0;j<6;j++) {
76    tmp *= pow(rtmp[j], (int)px[i+(j+3)*npx]);
77   }
78   r6[1] += tmp;
79  }
80  /* y final */
81  for(i=0;i<ny;i++) {
82   tmp = y[i+1*ny];
83   for(j=0;j<6;j++) {
84    tmp *= pow(rtmp[j], (int)y[i+(j+3)*ny]);
85   }
86   r6[2] += tmp;
87  }
88  /* py final */
89  for(i=0;i<npy;i++) {
90   tmp = py[i+1*npy];
91   for(j=0;j<6;j++) {
92    tmp *= pow(rtmp[j], (int)py[i+(j+3)*npy]);
93   }
94   r6[3] += tmp;
95  }
96  /* delta E final */
97  for(i=0;i<ndelta;i++) {
98   tmp = delta[i+1*ndelta];
99   for(j=0;j<6;j++) {
100    tmp *= pow(rtmp[j], (int)delta[i+(j+3)*ndelta]);
101   }
102   r6[4] += tmp;
103  }
104  /* t final */
105  for(i=0;i<nt;i++) {
106   tmp = t[i+1*nt];
107   for(j=0;j<6;j++) {
108    tmp *= pow(rtmp[j], (int)t[i+(j+3)*nt]);
109   }
110   r6[5] += tmp;
111  }
112
113  /* linearized misalignment transformation at the exit */
114  ATmultmv(r6,R2);
115  ATaddvv(r6,T2);
116 }
117 mxFree(rtmp);
118}
119
120ExportMode int* passFunction(const mxArray *ElemData,int *FieldNumbers,
121                             double *r_in, int num_particles, int mode)
122
123#define NUM_FIELDS_2_REMEMBER 17
124{
125 int *returnptr;
126 int fnum, *NewFieldNumbers;
127 mxArray *tmpmxptr;
128 double le;
129 double *pr1, *pr2, *pt1, *pt2;
130 double *x, *px, *y, *py, *t, *delta;
131 int nx, npx, ny, npy, nt, ndelta;
132
133 switch(mode) {
134  case NO_LOCAL_COPY: { /* Get fields by names from MATLAB workspace  */
135   tmpmxptr=mxGetField(ElemData,0,"Length");
136   if(tmpmxptr) {
137    le = mxGetScalar(tmpmxptr);
138    returnptr = NULL;
139   }
140   else
141    mexErrMsgTxt("Required field 'Length' was not found in the element data structure"); 
142   tmpmxptr=mxGetField(ElemData,0,"MapX");
143   if(tmpmxptr) {
144    x = mxGetPr(tmpmxptr);
145    returnptr = NULL;
146   }
147   else
148    mexErrMsgTxt("Required field 'MapX' was not found in the element data structure");
149   tmpmxptr=mxGetField(ElemData,0,"NumOfMonoX");
150   if(tmpmxptr) {
151    nx = (int)mxGetScalar(tmpmxptr);
152    returnptr = NULL;
153   }
154   else
155    mexErrMsgTxt("Required field 'NumOfMonoX' was not found in the element data structure");
156   tmpmxptr=mxGetField(ElemData,0,"MapPx");
157   if(tmpmxptr) {
158    px = mxGetPr(tmpmxptr);
159    returnptr = NULL;
160   }
161   else
162    mexErrMsgTxt("Required field 'MapPx' was not found in the element data structure");
163   tmpmxptr=mxGetField(ElemData,0,"NumOfMonoPx");
164   if(tmpmxptr) {
165    npx = (int)mxGetScalar(tmpmxptr);
166    returnptr = NULL;
167   }
168   else
169    mexErrMsgTxt("Required field 'NumOfMonoPx' was not found in the element data structure");
170   tmpmxptr=mxGetField(ElemData,0,"MapY");
171   if(tmpmxptr) {
172    y = mxGetPr(tmpmxptr);
173    returnptr = NULL;
174   }
175   else
176    mexErrMsgTxt("Required field 'MapY' was not found in the element data structure");
177   tmpmxptr=mxGetField(ElemData,0,"NumOfMonoY");
178   if(tmpmxptr) {
179    ny = (int)mxGetScalar(tmpmxptr);
180    returnptr = NULL;
181   }
182   else
183    mexErrMsgTxt("Required field 'NumOfMonoY' was not found in the element data structure");
184   tmpmxptr=mxGetField(ElemData,0,"MapPy");
185   if(tmpmxptr) {
186    py = mxGetPr(tmpmxptr);
187    returnptr = NULL;
188   }
189   else
190    mexErrMsgTxt("Required field 'MapPy' was not found in the element data structure");
191   tmpmxptr=mxGetField(ElemData,0,"NumOfMonoPy");
192   if(tmpmxptr) {
193    npy = (int)mxGetScalar(tmpmxptr);
194    returnptr = NULL;
195   }
196   else
197    mexErrMsgTxt("Required field 'NumOfMonoPy' was not found in the element data structure");
198   tmpmxptr=mxGetField(ElemData,0,"MapDelta");
199   if(tmpmxptr) {
200    delta = mxGetPr(tmpmxptr);
201    returnptr = NULL;
202   }
203   else
204    mexErrMsgTxt("Required field 'MapDelta' was not found in the element data structure");
205   tmpmxptr=mxGetField(ElemData,0,"NumOfMonoDelta");
206   if(tmpmxptr) {
207    ndelta = (int)mxGetScalar(tmpmxptr);
208    returnptr = NULL;
209   }
210   else
211    mexErrMsgTxt("Required field 'NumOfMonoDelta' was not found in the element data structure");
212   tmpmxptr=mxGetField(ElemData,0,"MapT");
213   if(tmpmxptr) {
214    t = mxGetPr(tmpmxptr);
215    returnptr = NULL;
216   }
217   else
218    mexErrMsgTxt("Required field 'MapT' was not found in the element data structure");
219   tmpmxptr=mxGetField(ElemData,0,"NumOfMonoT");
220   if(tmpmxptr) {
221    nt = (int)mxGetScalar(tmpmxptr);
222    returnptr = NULL;
223   }
224   else
225    mexErrMsgTxt("Required field 'NumOfMonoT' was not found in the element data structure");
226   tmpmxptr=mxGetField(ElemData,0,"T1");
227   if(tmpmxptr) {
228    pt1 = mxGetPr(tmpmxptr);
229    returnptr = NULL;
230   }
231   else
232    mexErrMsgTxt("Required field 'T1' was not found in the element data structure");
233   tmpmxptr=mxGetField(ElemData,0,"T2");
234   if(tmpmxptr) {
235    pt2 = mxGetPr(tmpmxptr);
236    returnptr = NULL;
237   }
238   else
239    mexErrMsgTxt("Required field 'T2' was not found in the element data structure");
240   tmpmxptr=mxGetField(ElemData,0,"R1");
241   if(tmpmxptr) {
242    pr1 = mxGetPr(tmpmxptr);
243    returnptr = NULL;
244   }
245   else
246    mexErrMsgTxt("Required field 'R1' was not found in the element data structure");
247   tmpmxptr=mxGetField(ElemData,0,"R2");
248   if(tmpmxptr) {
249    pr2 = mxGetPr(tmpmxptr);
250    returnptr = NULL;
251   }
252   else
253    mexErrMsgTxt("Required field 'R2' was not found in the element data structure");
254  }
255  break;
256  case MAKE_LOCAL_COPY: { /* Find field numbers first
257                             Save a list of field number in an array
258                             and make returnptr point to that array */
259   NewFieldNumbers = (int*)mxCalloc(NUM_FIELDS_2_REMEMBER,sizeof(int));
260   fnum = mxGetFieldNumber(ElemData,"Length");
261   if(fnum<0)
262    mexErrMsgTxt("Required field 'Length' was not found in the element data structure"); 
263   NewFieldNumbers[0] = fnum;
264   fnum = mxGetFieldNumber(ElemData,"MapX");
265   if(fnum<0)
266    mexErrMsgTxt("Required field 'MapX' was not found in the element data structure");
267   NewFieldNumbers[1] = fnum;
268   fnum = mxGetFieldNumber(ElemData,"NumOfMonoX");
269   if(fnum<0)
270    mexErrMsgTxt("Required field 'NumOfMonoX' was not found in the element data structure");
271   NewFieldNumbers[2] = fnum;
272   fnum = mxGetFieldNumber(ElemData,"MapPx");
273   if(fnum<0)
274    mexErrMsgTxt("Required field 'MapPx' was not found in the element data structure");
275   NewFieldNumbers[3] = fnum;
276   fnum = mxGetFieldNumber(ElemData,"NumOfMonoPx");
277   if(fnum<0)
278    mexErrMsgTxt("Required field 'NumOfMonoPx' was not found in the element data structure");
279   NewFieldNumbers[4] = fnum;
280   fnum = mxGetFieldNumber(ElemData,"MapY");
281   if(fnum<0)
282    mexErrMsgTxt("Required field 'MapY' was not found in the element data structure");
283   NewFieldNumbers[5] = fnum;
284   fnum = mxGetFieldNumber(ElemData,"NumOfMonoY");
285   if(fnum<0)
286    mexErrMsgTxt("Required field 'NumOfMonoY' was not found in the element data structure");
287   NewFieldNumbers[6] = fnum;
288   fnum = mxGetFieldNumber(ElemData,"MapPy");
289   if(fnum<0)
290    mexErrMsgTxt("Required field 'MapPy' was not found in the element data structure");
291   NewFieldNumbers[7] = fnum;
292   fnum = mxGetFieldNumber(ElemData,"NumOfMonoPy");
293   if(fnum<0)
294    mexErrMsgTxt("Required field 'NumOfMonoPy' was not found in the element data structure");
295   NewFieldNumbers[8] = fnum;
296   fnum = mxGetFieldNumber(ElemData,"MapDelta");
297   if(fnum<0)
298    mexErrMsgTxt("Required field 'MapDelta' was not found in the element data structure");
299   NewFieldNumbers[9] = fnum;
300   fnum = mxGetFieldNumber(ElemData,"NumOfMonoDelta");
301   if(fnum<0)
302    mexErrMsgTxt("Required field 'NumOfMonoDelta' was not found in the element data structure");
303   NewFieldNumbers[10] = fnum;
304   fnum = mxGetFieldNumber(ElemData,"MapT");
305   if(fnum<0)
306    mexErrMsgTxt("Required field 'MapT' was not found in the element data structure");
307   NewFieldNumbers[11] = fnum;
308   fnum = mxGetFieldNumber(ElemData,"NumOfMonoT");
309   if(fnum<0)
310    mexErrMsgTxt("Required field 'NumOfMonoT' was not found in the element data structure");
311   NewFieldNumbers[12] = fnum;
312   fnum = mxGetFieldNumber(ElemData,"T1");
313   if(fnum<0)
314    mexErrMsgTxt("Required field 'T1' was not found in the element data structure");
315   NewFieldNumbers[13] = fnum;
316   fnum = mxGetFieldNumber(ElemData,"T2");
317   if(fnum<0)
318    mexErrMsgTxt("Required field 'T2' was not found in the element data structure");
319   NewFieldNumbers[14] = fnum;
320   fnum = mxGetFieldNumber(ElemData,"R1");
321   if(fnum<0)
322    mexErrMsgTxt("Required field 'R1' was not found in the element data structure");
323   NewFieldNumbers[15] = fnum;
324   fnum = mxGetFieldNumber(ElemData,"R2");
325   if(fnum<0)
326    mexErrMsgTxt("Required field 'R2' was not found in the element data structure");
327   NewFieldNumbers[16] = fnum;
328
329   le = mxGetScalar(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[0]));
330   x = mxGetPr(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[1]));
331   nx = (int)mxGetScalar(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[2]));
332   px = mxGetPr(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[3]));
333   npx = (int)mxGetScalar(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[4]));
334   y = mxGetPr(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[5]));
335   ny = (int)mxGetScalar(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[6]));
336   py = mxGetPr(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[7]));
337   npy = (int)mxGetScalar(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[8]));
338   delta = mxGetPr(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[9]));
339   ndelta = (int)mxGetScalar(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[10]));
340   t = mxGetPr(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[11]));
341   nt = (int)mxGetScalar(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[12]));
342   pt1 = mxGetPr(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[13]));
343   pt2 = mxGetPr(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[14]));
344   pr1 = mxGetPr(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[15]));
345   pr2 = mxGetPr(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[16]));
346
347   returnptr = NewFieldNumbers;
348  }
349  break;
350  case USE_LOCAL_COPY: { /* Get fields from MATLAB using field numbers
351                            The second argument ponter to the array of field
352                            numbers is previously created with
353                            QuadLinPass( ..., MAKE_LOCAL_COPY) */
354   le = mxGetScalar(mxGetFieldByNumber(ElemData,0,FieldNumbers[0]));
355   x = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[1]));
356   nx = (int)mxGetScalar(mxGetFieldByNumber(ElemData,0,FieldNumbers[2]));
357   px = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[3]));
358   npx = (int)mxGetScalar(mxGetFieldByNumber(ElemData,0,FieldNumbers[4]));
359   y = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[5]));
360   ny = (int)mxGetScalar(mxGetFieldByNumber(ElemData,0,FieldNumbers[6]));
361   py = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[7]));
362   npy = (int)mxGetScalar(mxGetFieldByNumber(ElemData,0,FieldNumbers[8]));
363   delta = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[9]));
364   ndelta = (int)mxGetScalar(mxGetFieldByNumber(ElemData,0,FieldNumbers[10]));
365   t = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[11]));
366   nt = (int)mxGetScalar(mxGetFieldByNumber(ElemData,0,FieldNumbers[12]));
367   pt1 = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[13]));
368   pt2 = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[14]));
369   pr1 = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[15]));
370   pr2 = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[16]));
371   returnptr = FieldNumbers;
372  }
373  break;
374  default: {
375   mexErrMsgTxt("No match for calling mode in function TaylormapPass\n");
376  }
377 }
378 TaylormapPass(r_in, le, x, nx, px, npx, y, ny, py, npy, delta, ndelta, t, nt,
379                    pt1, pt2, pr1, pr2, num_particles);
380 return(returnptr);
381}
382
383void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
384{
385 int m,n; double *r_in;
386 mxArray *tmpmxptr;
387 double le;
388 double *pr1, *pr2, *pt1, *pt2;
389 double *x, *px, *y, *py, *t, *delta;
390 int nx, npx, ny, npy, nt, ndelta;
391
392 /* ALLOCATE memory for the output array of the same size as the input  */
393 m = mxGetM(prhs[1]);
394 n = mxGetN(prhs[1]);
395 if(m!=6) 
396  mexErrMsgTxt("Second argument must be a 6 x N matrix");
397 tmpmxptr=mxGetField(prhs[0],0,"Length");
398 if(tmpmxptr)
399   le = mxGetScalar(tmpmxptr);
400 else
401   mexErrMsgTxt("Required field 'Length' was not found in the element data structure"); 
402 tmpmxptr=mxGetField(prhs[0],0,"MapX");
403 if(tmpmxptr)
404   x = mxGetPr(tmpmxptr);
405 else
406   mexErrMsgTxt("Required field 'MapX' was not found in the element data structure");
407 tmpmxptr=mxGetField(prhs[0],0,"NumOfMonoX");
408 if(tmpmxptr)
409   nx = (int)mxGetScalar(tmpmxptr);
410 else
411   mexErrMsgTxt("Required field 'NumOfMonoX' was not found in the element data structure");
412 tmpmxptr=mxGetField(prhs[0],0,"MapPx");
413 if(tmpmxptr)
414   px = mxGetPr(tmpmxptr);
415 else
416   mexErrMsgTxt("Required field 'MapPx' was not found in the element data structure");
417 tmpmxptr=mxGetField(prhs[0],0,"NumOfMonoPx");
418 if(tmpmxptr)
419   npx = (int)mxGetScalar(tmpmxptr);
420 else
421   mexErrMsgTxt("Required field 'NumOfMonoPx' was not found in the element data structure");
422 tmpmxptr=mxGetField(prhs[0],0,"MapY");
423 if(tmpmxptr)
424   y = mxGetPr(tmpmxptr);
425 else
426   mexErrMsgTxt("Required field 'MapY' was not found in the element data structure");
427 tmpmxptr=mxGetField(prhs[0],0,"NumOfMonoY");
428 if(tmpmxptr)
429   ny = (int)mxGetScalar(tmpmxptr);
430 else
431   mexErrMsgTxt("Required field 'NumOfMonoY' was not found in the element data structure");
432 tmpmxptr=mxGetField(prhs[0],0,"MapPy");
433 if(tmpmxptr)
434   py = mxGetPr(tmpmxptr);
435 else
436   mexErrMsgTxt("Required field 'MapPy' was not found in the element data structure");
437 tmpmxptr=mxGetField(prhs[0],0,"NumOfMonoPy");
438 if(tmpmxptr)
439   npy = (int)mxGetScalar(tmpmxptr);
440 else
441   mexErrMsgTxt("Required field 'NumOfMonoPy' was not found in the element data structure");
442 tmpmxptr=mxGetField(prhs[0],0,"MapDelta");
443 if(tmpmxptr)
444   delta = mxGetPr(tmpmxptr);
445 else
446   mexErrMsgTxt("Required field 'MapDelta' was not found in the element data structure");
447 tmpmxptr=mxGetField(prhs[0],0,"NumOfMonoDelta");
448 if(tmpmxptr)
449   ndelta = (int)mxGetScalar(tmpmxptr);
450 else
451   mexErrMsgTxt("Required field 'NumOfMonoDelta' was not found in the element data structure");
452 tmpmxptr=mxGetField(prhs[0],0,"MapT");
453 if(tmpmxptr)
454   t = mxGetPr(tmpmxptr);
455 else
456   mexErrMsgTxt("Required field 'MapT' was not found in the element data structure");
457 tmpmxptr=mxGetField(prhs[0],0,"NumOfMonoT");
458 if(tmpmxptr)
459   nt = (int)mxGetScalar(tmpmxptr);
460 else
461   mexErrMsgTxt("Required field 'NumOfMonoT' was not found in the element data structure");
462 tmpmxptr=mxGetField(prhs[0],0,"R1");
463 if(tmpmxptr)
464   pr1 = mxGetPr(tmpmxptr);
465 else
466   mexErrMsgTxt("Required field 'R1' was not found in the element data structure");
467 tmpmxptr=mxGetField(prhs[0],0,"R2");
468 if(tmpmxptr)
469   pr2 = mxGetPr(tmpmxptr);
470 else
471   mexErrMsgTxt("Required field 'R2' was not found in the element data structure");
472 tmpmxptr=mxGetField(prhs[0],0,"T1");
473 if(tmpmxptr)
474   pt1 = mxGetPr(tmpmxptr);
475 else
476   mexErrMsgTxt("Required field 'T1' was not found in the element data structure");
477 tmpmxptr=mxGetField(prhs[0],0,"T2");
478 if(tmpmxptr)
479   pt2 = mxGetPr(tmpmxptr);
480 else
481   mexErrMsgTxt("Required field 'T2' was not found in the element data structure");
482 plhs[0] = mxDuplicateArray(prhs[1]);
483 r_in = mxGetPr(plhs[0]);
484 TaylormapPass(r_in, le, x, nx, px, npx, y, ny, py, npy, delta, ndelta, t, nt,
485                    pt1, pt2, pr1, pr2, n);
486}
Note: See TracBrowser for help on using the repository browser.