source: MML/trunk/at/simulator/element/user/ThinEPUPass.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: 13.3 KB
Line 
1/* ThinEPUPass.c
2   Accelerator Toolbox
3*/
4
5#include <stdlib.h>
6#include <stdio.h>
7#include <string.h>
8#include <math.h>
9#include "mex.h"
10#include "elempass.h"
11#include "../atlalib.c"
12
13void fill(double *pout, double *pin, int m, int n)
14{ int i, j;
15  for (i=0; i<m; i++) {
16     for (j=0; j<n; j++) {
17        pout[i+j*m] = pin[i+j*m];
18     }
19  }
20}
21
22void ThinEPUPass(double *r, int nx, int ny, double *XEPU, double *YEPU, double *PXEPU, double *PYEPU,
23                 double *R1, double *R2, double *T1, double *T2, int num_particles)
24
25{ int c;
26  double *r6, *xi, *yi;   
27  char *method;
28  int buflen;
29  double xkick, ykick;
30  int num_in, num_out;
31  mxArray *input_array[6], *output_array[1];
32
33  num_in = 6;
34  num_out = 1;
35  input_array[0] = mxCreateDoubleMatrix(ny,nx,mxREAL);
36  input_array[1] = mxCreateDoubleMatrix(ny,nx,mxREAL);
37  input_array[2] = mxCreateDoubleMatrix(ny,nx,mxREAL);
38  input_array[3] = mxCreateDoubleMatrix(1,1,mxREAL);
39  input_array[4] = mxCreateDoubleMatrix(1,1,mxREAL);
40  xi = (double*)mxCalloc(1,sizeof(double));
41  yi = (double*)mxCalloc(1,sizeof(double));
42  fill(mxGetPr(input_array[0]), XEPU, ny, nx);
43  fill(mxGetPr(input_array[1]), YEPU, ny, nx);
44  mxSetM(input_array[0], ny) ;
45  mxSetN(input_array[0], nx) ;
46  mxSetM(input_array[1], ny) ;
47  mxSetN(input_array[1], nx) ;
48  buflen = 6;
49  method = mxCalloc(buflen,sizeof(char));
50  strncpy(method,"*cubic",buflen);
51  /*
52  mexPrintf(" method = ");
53  for(c = 0;c<buflen;c++) {
54     mexPrintf("%c", method[c]);
55  }
56  mexPrintf("\n");
57  */
58  input_array[5] = mxCreateString(method);
59  for (c = 0;c<num_particles;c++){/* Loop over particles */
60     r6 = r+c*6;
61     /* linearized misalignment transformation at the entrance */
62     ATaddvv(r6,T1);
63     ATmultmv(r6,R1);
64     xi[0] = r6[0]; 
65     yi[0] = r6[2]; 
66     fill(mxGetPr(input_array[3]), xi, 1, 1);
67     fill(mxGetPr(input_array[4]), yi, 1, 1);
68     mxSetM(input_array[3], 1) ;
69     mxSetN(input_array[3], 1) ;
70     mxSetM(input_array[4], 1) ;
71     mxSetN(input_array[4], 1) ;
72     fill(mxGetPr(input_array[2]), PXEPU, ny, nx);
73     mxSetM(input_array[2], ny) ;
74     mxSetN(input_array[2], nx) ;
75     mexCallMATLAB(num_out, output_array, num_in, input_array, "interp2");
76     xkick = mxGetScalar(output_array[0]);
77     fill(mxGetPr(input_array[2]), PYEPU, ny, nx);
78     mxSetM(input_array[2], ny) ;
79     mxSetN(input_array[2], nx) ;
80     mexCallMATLAB(num_out, output_array, num_in, input_array, "interp2");
81     ykick = mxGetScalar(output_array[0]);
82
83     r6[1] += xkick;
84     r6[3] += ykick;
85           
86     /* linearized misalignment transformation at the exit */
87     ATmultmv(r6,R2);
88     ATaddvv(r6,T2);
89  }
90  mxDestroyArray(input_array[5]);
91  mxFree(method);
92  mxFree(yi);
93  mxFree(xi);
94  mxDestroyArray(input_array[4]);
95  mxDestroyArray(input_array[3]);
96  mxDestroyArray(input_array[2]);
97  mxDestroyArray(input_array[1]);
98  mxDestroyArray(input_array[0]);
99}
100
101ExportMode int* passFunction(const mxArray *ElemData, int *FieldNumbers,
102double *r_in, int num_particles, int mode)
103
104#define NUM_FIELDS_2_REMEMBER 11
105{
106  double le;
107  int nx, ny;
108  double *XEPU, *YEPU, *PXEPU, *PYEPU;
109  double *pr1, *pr2, *pt1, *pt2;
110  int *returnptr;
111  int *NewFieldNumbers,fnum;
112  mxArray *tmpmxptr;
113
114  switch(mode) {
115     case NO_LOCAL_COPY: { /* Get fields by names from MATLAB workspace  */
116        tmpmxptr=mxGetField(ElemData,0,"Length");
117        if(tmpmxptr) {
118           le = mxGetScalar(tmpmxptr);
119           returnptr = NULL;
120        }
121        else
122           mexErrMsgTxt("Required field 'Length' was not found in the element data structure");
123        tmpmxptr =mxGetField(ElemData,0,"NumX");
124        if(tmpmxptr) {
125           nx = (int)mxGetScalar(tmpmxptr);
126           returnptr = NULL;
127        }
128        else
129           mexErrMsgTxt("Required field 'NumX' was not found in the element data structure"); 
130        tmpmxptr =mxGetField(ElemData,0,"NumY");
131        if(tmpmxptr) { 
132           ny = (int)mxGetScalar(tmpmxptr);
133           returnptr = NULL;
134        }
135        else
136           mexErrMsgTxt("Required field 'NumY' was not found in the element data structure");
137        tmpmxptr = mxGetField(ElemData,0,"XGrid");
138        if(tmpmxptr) {
139           XEPU = mxGetPr(tmpmxptr);
140           returnptr = NULL;
141        }
142        else
143           mexErrMsgTxt("Required field 'XGrid' was not found in the element data structure");
144        tmpmxptr = mxGetField(ElemData,0,"YGrid");
145        if(tmpmxptr) {
146           YEPU = mxGetPr(tmpmxptr);
147           returnptr = NULL;
148        }
149        else
150           mexErrMsgTxt("Required field 'YGrid' was not found in the element data structure");
151        tmpmxptr = mxGetField(ElemData,0,"PxGrid");
152        if(tmpmxptr) {
153           PXEPU = mxGetPr(tmpmxptr);
154           returnptr = NULL;
155        }
156        else
157           mexErrMsgTxt("Required field 'PxGrid' was not found in the element data structure");
158        tmpmxptr = mxGetField(ElemData,0,"PyGrid");
159        if(tmpmxptr) {
160           PYEPU = mxGetPr(tmpmxptr);
161           returnptr = NULL;
162        }
163        else
164           mexErrMsgTxt("Required field 'PyGrid' was not found in the element data structure");
165        tmpmxptr=mxGetField(ElemData,0,"R1");
166        if(tmpmxptr) {
167           pr1 = mxGetPr(tmpmxptr);
168           returnptr = NULL;
169        }
170        else
171           mexErrMsgTxt("Required field 'R1' was not found in the element data structure");
172        tmpmxptr=mxGetField(ElemData,0,"R2");
173        if(tmpmxptr) {
174           pr2 = mxGetPr(tmpmxptr);
175           returnptr = NULL;
176        }
177        else
178           mexErrMsgTxt("Required field 'R2' was not found in the element data structure");
179        tmpmxptr=mxGetField(ElemData,0,"T1");
180        if(tmpmxptr) {
181           pt1 = mxGetPr(tmpmxptr);
182           returnptr = NULL;
183        }
184        else
185           mexErrMsgTxt("Required field 'T1' was not found in the element data structure");
186        tmpmxptr=mxGetField(ElemData,0,"T2");
187        if(tmpmxptr) {
188           pt2 = mxGetPr(tmpmxptr);
189           returnptr = NULL;
190        }
191        else
192           mexErrMsgTxt("Required field 'T2' was not found in the element data structure");
193     }
194     break;
195     case MAKE_LOCAL_COPY: { /* Find field numbers first
196                                Save a list of field number in an array
197                                and make returnptr point to that array */
198        /* Allocate memory for integer array of
199           field numbers for faster futurereference */
200
201        NewFieldNumbers = (int*)mxCalloc(NUM_FIELDS_2_REMEMBER,sizeof(int));
202
203        /* Populate */
204        fnum = mxGetFieldNumber(ElemData,"Length");
205        if(fnum<0)
206           mexErrMsgTxt("Required field 'Length' was not found in the element data structure");
207        NewFieldNumbers[0] = fnum;
208        fnum = mxGetFieldNumber(ElemData,"NumX");
209        if(fnum<0) 
210           mexErrMsgTxt("Required field 'NumX' was not found in the element data structure"); 
211        NewFieldNumbers[1] = fnum;
212        fnum = mxGetFieldNumber(ElemData,"NumY");
213        if(fnum<0) 
214           mexErrMsgTxt("Required field 'NumY' was not found in the element data structure"); 
215        NewFieldNumbers[2] = fnum;
216        fnum = mxGetFieldNumber(ElemData,"XGrid");
217        if(fnum<0) 
218           mexErrMsgTxt("Required field 'XGrid' was not found in the element data structure"); 
219        NewFieldNumbers[3] = fnum;
220        fnum = mxGetFieldNumber(ElemData,"YGrid");
221        if(fnum<0) 
222           mexErrMsgTxt("Required field 'YGrid' was not found in the element data structure"); 
223        NewFieldNumbers[4] = fnum;
224        fnum = mxGetFieldNumber(ElemData,"PxGrid");
225        if(fnum<0)
226           mexErrMsgTxt("Required field 'PxGrid' was not found in the element data structure");
227        NewFieldNumbers[5] = fnum;
228        fnum = mxGetFieldNumber(ElemData,"PyGrid");
229        if(fnum<0)
230           mexErrMsgTxt("Required field 'PyGrid' was not found in the element data structure");
231        NewFieldNumbers[6] = fnum;
232        fnum = mxGetFieldNumber(ElemData,"R1");
233        if(fnum<0)
234           mexErrMsgTxt("Required field 'R1' was not found in the element data structure");
235        NewFieldNumbers[7] = fnum;
236        fnum = mxGetFieldNumber(ElemData,"R2");
237        if(fnum<0)
238           mexErrMsgTxt("Required field 'R2' was not found in the element data structure");
239        NewFieldNumbers[8] = fnum;
240        fnum = mxGetFieldNumber(ElemData,"T1");
241        if(fnum<0)
242           mexErrMsgTxt("Required field 'T1' was not found in the element data structure");
243        NewFieldNumbers[9] = fnum;
244        fnum = mxGetFieldNumber(ElemData,"T2");
245        if(fnum<0)
246           mexErrMsgTxt("Required field 'T2' was not found in the element data structure");
247        NewFieldNumbers[10] = fnum;
248
249        le = mxGetScalar(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[0]));
250        nx = (int)mxGetScalar(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[1]));
251        ny = (int)mxGetScalar(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[2]));
252        XEPU = mxGetPr(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[3]));
253        YEPU = mxGetPr(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[4]));
254        PXEPU = mxGetPr(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[5]));
255        PYEPU = mxGetPr(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[6]));
256        pr1 = mxGetPr(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[7]));
257        pr2 = mxGetPr(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[8]));
258        pt1 = mxGetPr(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[9]));
259        pt2 = mxGetPr(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[10]));
260
261        returnptr = NewFieldNumbers;
262     }
263     break;
264     case USE_LOCAL_COPY: { /* Get fields from MATLAB using field numbers
265                               The second argument ponter to the array of field
266                               numbers is previously created with
267                               QuadLinPass( ..., MAKE_LOCAL_COPY) */
268
269        le = mxGetScalar(mxGetFieldByNumber(ElemData,0,FieldNumbers[0]));
270        nx = (int)mxGetScalar(mxGetFieldByNumber(ElemData,0,FieldNumbers[1]));
271        ny = (int)mxGetScalar(mxGetFieldByNumber(ElemData,0,FieldNumbers[2]));
272        XEPU = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[3]));
273        YEPU = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[4]));
274        PXEPU = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[5]));
275        PYEPU = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[6]));
276        pr1 = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[7]));
277        pr2 = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[8]));
278        pt1 = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[9]));
279        pt2 = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[10]));
280
281        returnptr = FieldNumbers;
282     }
283     break;
284     default: {
285        mexErrMsgTxt("No match for calling mode in function ThinEPUPass\n");
286     }
287  }
288  ThinEPUPass(r_in, nx, ny, XEPU, YEPU, PXEPU, PYEPU, pr1, pr2, pt1, pt2, num_particles);
289
290  return(returnptr);
291}
292
293void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
294{
295  int m,n;
296  double *r_in;
297  mxArray *tmpmxptr;
298  double le;
299  int nx, ny;
300  double *XEPU, *YEPU, *PXEPU, *PYEPU;
301  double *pr1, *pr2, *pt1, *pt2;
302
303  /* ALLOCATE memory for the output array of the same size as the input */
304  m = mxGetM(prhs[1]);
305  n = mxGetN(prhs[1]);
306  if(m!=6) 
307     mexErrMsgTxt("Second argument must be a 6 x N matrix");
308
309  tmpmxptr = mxGetField(prhs[0],0,"Length");
310  if(tmpmxptr)
311     le = mxGetScalar(tmpmxptr);
312  else
313     mexErrMsgTxt("Required field 'Length' was not found in the element data structure");
314  tmpmxptr = mxGetField(prhs[0],0,"NumX");
315  if(tmpmxptr)
316     nx = (int)mxGetScalar(tmpmxptr);
317  else
318     mexErrMsgTxt("Required field 'NumX' was not found in the element data structure");
319  tmpmxptr = mxGetField(prhs[0],0,"NumY");
320  if(tmpmxptr)
321     ny = (int)mxGetScalar(tmpmxptr);
322  else
323     mexErrMsgTxt("Required field 'NumY' was not found in the element data structure");
324  tmpmxptr=mxGetField(prhs[0],0,"XGrid");
325  if(tmpmxptr)
326     XEPU = mxGetPr(tmpmxptr);
327  else
328     mexErrMsgTxt("Required field 'XGrid' was not found in the element data structure"); 
329  tmpmxptr=mxGetField(prhs[0],0,"YGrid");
330  if(tmpmxptr)
331     YEPU = mxGetPr(tmpmxptr);
332  else
333     mexErrMsgTxt("Required field 'YGrid' was not found in the element data structure");
334  tmpmxptr=mxGetField(prhs[0],0,"PxGrid");
335  if(tmpmxptr)
336     PXEPU = mxGetPr(tmpmxptr);
337  else
338     mexErrMsgTxt("Required field 'PxGrid' was not found in the element data structure"); 
339  tmpmxptr=mxGetField(prhs[0],0,"PyGrid");
340  if(tmpmxptr)
341     PYEPU = mxGetPr(tmpmxptr);
342  else
343     mexErrMsgTxt("Required field 'PyGrid' was not found in the element data structure");
344  tmpmxptr=mxGetField(prhs[0],0,"R1");
345  if(tmpmxptr)
346     pr1 = mxGetPr(tmpmxptr);
347  else
348     mexErrMsgTxt("Required field 'R1' was not found in the element data structure");
349  tmpmxptr=mxGetField(prhs[0],0,"R2");
350  if(tmpmxptr)
351     pr2 = mxGetPr(tmpmxptr);
352  else
353     mexErrMsgTxt("Required field 'R2' was not found in the element data structure");
354  tmpmxptr=mxGetField(prhs[0],0,"T1");
355  if(tmpmxptr)
356     pt1 = mxGetPr(tmpmxptr);
357  else
358     mexErrMsgTxt("Required field 'T1' was not found in the element data structure");
359  tmpmxptr=mxGetField(prhs[0],0,"T2");
360  if(tmpmxptr)
361     pt2 = mxGetPr(tmpmxptr);
362  else
363     mexErrMsgTxt("Required field 'T2' was not found in the element data structure");
364  plhs[0] = mxDuplicateArray(prhs[1]);
365  r_in = mxGetPr(plhs[0]);
366  ThinEPUPass(r_in, nx, ny, XEPU, YEPU, PXEPU, PYEPU, pr1, pr2, pt1, pt2, n);
367}
Note: See TracBrowser for help on using the repository browser.