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