source: MML/trunk/at/simulator/element/user/TaylormapPass3.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: 21.1 KB
Line 
1/* TaylormapPass3.c
2   Accelerator Toolbox
3   Created on 05/06/2004
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 TaylormapPass3(double *r_in, double le, int no, int nv,
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   no -- order of the Taylor map
25   nv -- number of variables in the input
26   x -- Taylor map of final x
27   nx -- number of monomials in the map x
28   px -- Taylor map of final px
29   npx -- number of monomials in the map px
30   y -- Taylor map of final y
31   ny -- number of monomials in the map y
32   py -- Taylor map of final py
33   npy -- number of monomials in the map py
34   delta -- Taylor map of final delta p / p
35   ndelta -- number of monomials in the map delta
36   t -- Taylor map of final t
37   nt -- number of monomials in the map t
38   T1 -- 6 x 1 translation at entrance
39   T2 -- 6 x 1 translation at exit
40   R1 -- 6 x 6 rotation matrix at the entrance
41   R2 -- 6 x 6 rotation matrix at the exit
42   num_particles -- number of particles
43*/
44#define NUM_COLUMNS 9
45{
46 int c,i,j,itmp;
47 int *ptmpx;
48 int *ptmppx;
49 int *ptmpy;
50 int *ptmppy;
51 int *ptmpdelta;
52 int *ptmpt;
53 double *r6; 
54 double *rtmp;
55 double tmp;
56 double *product;
57
58 rtmp = (double*)mxCalloc(6,sizeof(double));
59 ptmpx = (int*)mxCalloc(nx,sizeof(int));
60 ptmppx = (int*)mxCalloc(npx,sizeof(int));
61 ptmpy = (int*)mxCalloc(ny,sizeof(int));
62 ptmppy = (int*)mxCalloc(npy,sizeof(int));
63 ptmpdelta = (int*)mxCalloc(ndelta,sizeof(int));
64 ptmpt = (int*)mxCalloc(nt,sizeof(int));
65 product = (double*)mxCalloc( (int)pow(no+1,nv),sizeof(double));
66
67 for(c = 0;c<num_particles;c++) {
68  r6 = r_in+c*6;
69  for(i=0;i<6;i++) {
70   rtmp[i] = r6[i];
71   r6[i] = 0;
72  }
73
74  /* linearized misalignment transformation at the entrance */
75  ATaddvv(rtmp,T1);
76  ATmultmv(rtmp,R1);
77
78  /* Taylor map */
79
80  for(i=0;i<((int)pow(no+1,nv));i++) {
81   product[i] = 0;
82  }
83
84  product[0] = 1;
85  itmp = 1;
86  for(i=0;i<nv;i++) {
87   tmp = 1;
88   for(j=0;j<no;j++) {
89    tmp *= rtmp[i];
90    product[(j+1)*itmp] = tmp;
91   }
92   itmp *= no+1;
93  }
94
95  /* x */
96  for(i=0;i<nx;i++) {
97   itmp = 1;
98   ptmpx[i] = 0;
99   for(j=0;j<nv;j++) {
100    ptmpx[i] += (int)x[i+(j+3)*nx]*itmp;
101    itmp *= no+1;
102   }
103   if (product[ptmpx[i]]==0) {
104    itmp = 1;
105    tmp = 1;
106    for(j=0;j<nv;j++) {
107     if ((int)x[i+(j+3)*nx]>0) {
108      tmp *= product[(int)x[i+(j+3)*nx]*itmp];
109     }
110     itmp *= no+1;
111    }
112    product[ptmpx[i]] = tmp;
113   }
114  }
115  /* px */
116  for(i=0;i<npx;i++) {
117   itmp = 1;
118   ptmppx[i] = 0;
119   for(j=0;j<nv;j++) {
120    ptmppx[i] += (int)px[i+(j+3)*npx]*itmp;
121    itmp *= no+1;
122   }
123   if (product[ptmppx[i]]==0) {
124    itmp = 1;
125    tmp = 1;
126    for(j=0;j<nv;j++) {
127     if ((int)px[i+(j+3)*npx]>0) {
128      tmp *= product[(int)px[i+(j+3)*npx]*itmp];
129     }
130     itmp *= no+1;
131    }
132    product[ptmppx[i]] = tmp;
133   }
134  }
135  /* y */
136  for(i=0;i<ny;i++) {
137   itmp = 1;
138   ptmpy[i] = 0;
139   for(j=0;j<nv;j++) {
140    ptmpy[i] += (int)y[i+(j+3)*ny]*itmp;
141    itmp *= no+1;
142   }
143   if (product[ptmpy[i]]==0) {
144    itmp = 1;
145    tmp = 1;
146    for(j=0;j<nv;j++) {
147     if ((int)y[i+(j+3)*ny]>0) {
148      tmp *= product[(int)y[i+(j+3)*ny]*itmp];
149     }
150     itmp *= no+1;
151    }
152    product[ptmpy[i]] = tmp;
153   }
154  }
155  /* py */
156  for(i=0;i<npy;i++) {
157   itmp = 1;
158   ptmppy[i] = 0;
159   for(j=0;j<nv;j++) {
160    ptmppy[i] += (int)py[i+(j+3)*npy]*itmp;
161    itmp *= no+1;
162   }
163   if (product[ptmppy[i]]==0) {
164    itmp = 1;
165    tmp = 1;
166    for(j=0;j<nv;j++) {
167     if ((int)py[i+(j+3)*npy]>0) {
168      tmp *= product[(int)py[i+(j+3)*npy]*itmp];
169     }
170     itmp *= no+1;
171    }
172    product[ptmppy[i]] = tmp;
173   }
174  }
175  /* delta */
176  for(i=0;i<ndelta;i++) {
177   itmp = 1;
178   ptmpdelta[i] = 0;
179   for(j=0;j<nv;j++) {
180    ptmpdelta[i] += (int)delta[i+(j+3)*ndelta]*itmp;
181    itmp *= no+1;
182   }
183   if (product[ptmpdelta[i]]==0) {
184    itmp = 1;
185    tmp = 1;
186    for(j=0;j<nv;j++) {
187     if ((int)delta[i+(j+3)*ndelta]>0) {
188      tmp *= product[(int)delta[i+(j+3)*ndelta]*itmp];
189     }
190     itmp *= no+1;
191    }
192    product[ptmpdelta[i]] = tmp;
193   }
194  }
195  /* t */
196  for(i=0;i<nt;i++) {
197   itmp = 1;
198   ptmpt[i] = 0;
199   for(j=0;j<nv;j++) {
200    ptmpt[i] += (int)t[i+(j+3)*nt]*itmp;
201    itmp *= no+1;
202   }
203   if (product[ptmpt[i]]==0) {
204    itmp = 1;
205    tmp = 1;
206    for(j=0;j<nv;j++) {
207     if ((int)t[i+(j+3)*nt]>0) {
208      tmp *= product[(int)t[i+(j+3)*nt]*itmp];
209     }
210     itmp *= no+1;
211    }
212    product[ptmpt[i]] = tmp;
213   }
214  }
215
216
217  /* x final */
218  for(i=0;i<nx;i++) {
219   tmp = x[i+1*nx]*product[ptmpx[i]];
220   r6[0] += tmp;
221  }
222
223  /* px final */
224  for(i=0;i<npx;i++) {
225   tmp = px[i+1*npx]*product[ptmppx[i]];
226   r6[1] += tmp;
227  }
228
229  /* y final */
230  for(i=0;i<ny;i++) {
231   tmp = y[i+1*ny]*product[ptmpy[i]];
232   r6[2] += tmp;
233  }
234
235  /* py final */
236  for(i=0;i<npy;i++) {
237   tmp = py[i+1*npy]*product[ptmppy[i]];
238   r6[3] += tmp;
239  }
240
241  /* delta final */
242  for(i=0;i<ndelta;i++) {
243   tmp = delta[i+1*ndelta]*product[ptmpdelta[i]];
244   r6[4] += tmp;
245  }
246
247  /* t final */
248  for(i=0;i<nt;i++) {
249   tmp = t[i+1*nt]*product[ptmpt[i]];
250   r6[5] += tmp;
251  }
252
253  /* linearized misalignment transformation at the exit */
254  ATmultmv(r6,R2);
255  ATaddvv(r6,T2);
256 }
257 mxFree(rtmp);
258 mxFree(ptmpx);
259 mxFree(ptmppx);
260 mxFree(ptmpy);
261 mxFree(ptmppy);
262 mxFree(ptmpdelta);
263 mxFree(ptmpt);
264 mxFree(product);
265}
266
267ExportMode int* passFunction(const mxArray *ElemData,int *FieldNumbers,
268                             double *r_in, int num_particles, int mode)
269
270#define NUM_FIELDS_2_REMEMBER 19
271{
272 int *returnptr;
273 int fnum, *NewFieldNumbers;
274 mxArray *tmpmxptr;
275 double le;
276 int no, nv;
277 double *pr1, *pr2, *pt1, *pt2;
278 double *x, *px, *y, *py, *t, *delta;
279 int nx, npx, ny, npy, nt, ndelta;
280
281 switch(mode) {
282  case NO_LOCAL_COPY: { /* Get fields by names from MATLAB workspace  */
283   tmpmxptr=mxGetField(ElemData,0,"Length");
284   if(tmpmxptr) {
285    le = mxGetScalar(tmpmxptr);
286    returnptr = NULL;
287   }
288   else
289    mexErrMsgTxt("Required field 'Length' was not found in the element data structure"); 
290   tmpmxptr=mxGetField(ElemData,0,"Order");
291   if(tmpmxptr) {
292    no = (int)mxGetScalar(tmpmxptr);
293    returnptr = NULL;
294   }
295   else
296    mexErrMsgTxt("Required field 'Order' was not found in the element data structure");
297   tmpmxptr=mxGetField(ElemData,0,"NumOfVar");
298   if(tmpmxptr) {
299    nv = (int)mxGetScalar(tmpmxptr);
300    returnptr = NULL;
301   }
302   else
303    mexErrMsgTxt("Required field 'NumOfVar' was not found in the element data structure");
304   tmpmxptr=mxGetField(ElemData,0,"MapX");
305   if(tmpmxptr) {
306    x = mxGetPr(tmpmxptr);
307    returnptr = NULL;
308   }
309   else
310    mexErrMsgTxt("Required field 'MapX' was not found in the element data structure");
311   tmpmxptr=mxGetField(ElemData,0,"NumOfMonoX");
312   if(tmpmxptr) {
313    nx = (int)mxGetScalar(tmpmxptr);
314    returnptr = NULL;
315   }
316   else
317    mexErrMsgTxt("Required field 'NumOfMonoX' was not found in the element data structure");
318   tmpmxptr=mxGetField(ElemData,0,"MapPx");
319   if(tmpmxptr) {
320    px = mxGetPr(tmpmxptr);
321    returnptr = NULL;
322   }
323   else
324    mexErrMsgTxt("Required field 'MapPx' was not found in the element data structure");
325   tmpmxptr=mxGetField(ElemData,0,"NumOfMonoPx");
326   if(tmpmxptr) {
327    npx = (int)mxGetScalar(tmpmxptr);
328    returnptr = NULL;
329   }
330   else
331    mexErrMsgTxt("Required field 'NumOfMonoPx' was not found in the element data structure");
332   tmpmxptr=mxGetField(ElemData,0,"MapY");
333   if(tmpmxptr) {
334    y = mxGetPr(tmpmxptr);
335    returnptr = NULL;
336   }
337   else
338    mexErrMsgTxt("Required field 'MapY' was not found in the element data structure");
339   tmpmxptr=mxGetField(ElemData,0,"NumOfMonoY");
340   if(tmpmxptr) {
341    ny = (int)mxGetScalar(tmpmxptr);
342    returnptr = NULL;
343   }
344   else
345    mexErrMsgTxt("Required field 'NumOfMonoY' was not found in the element data structure");
346   tmpmxptr=mxGetField(ElemData,0,"MapPy");
347   if(tmpmxptr) {
348    py = mxGetPr(tmpmxptr);
349    returnptr = NULL;
350   }
351   else
352    mexErrMsgTxt("Required field 'MapPy' was not found in the element data structure");
353   tmpmxptr=mxGetField(ElemData,0,"NumOfMonoPy");
354   if(tmpmxptr) {
355    npy = (int)mxGetScalar(tmpmxptr);
356    returnptr = NULL;
357   }
358   else
359    mexErrMsgTxt("Required field 'NumOfMonoPy' was not found in the element data structure");
360   tmpmxptr=mxGetField(ElemData,0,"MapDelta");
361   if(tmpmxptr) {
362    delta = mxGetPr(tmpmxptr);
363    returnptr = NULL;
364   }
365   else
366    mexErrMsgTxt("Required field 'MapDelta' was not found in the element data structure");
367   tmpmxptr=mxGetField(ElemData,0,"NumOfMonoDelta");
368   if(tmpmxptr) {
369    ndelta = (int)mxGetScalar(tmpmxptr);
370    returnptr = NULL;
371   }
372   else
373    mexErrMsgTxt("Required field 'NumOfMonoDelta' was not found in the element data structure");
374   tmpmxptr=mxGetField(ElemData,0,"MapT");
375   if(tmpmxptr) {
376    t = mxGetPr(tmpmxptr);
377    returnptr = NULL;
378   }
379   else
380    mexErrMsgTxt("Required field 'MapT' was not found in the element data structure");
381   tmpmxptr=mxGetField(ElemData,0,"NumOfMonoT");
382   if(tmpmxptr) {
383    nt = (int)mxGetScalar(tmpmxptr);
384    returnptr = NULL;
385   }
386   else
387    mexErrMsgTxt("Required field 'NumOfMonoT' was not found in the element data structure");
388   tmpmxptr=mxGetField(ElemData,0,"T1");
389   if(tmpmxptr) {
390    pt1 = mxGetPr(tmpmxptr);
391    returnptr = NULL;
392   }
393   else
394    mexErrMsgTxt("Required field 'T1' was not found in the element data structure");
395   tmpmxptr=mxGetField(ElemData,0,"T2");
396   if(tmpmxptr) {
397    pt2 = mxGetPr(tmpmxptr);
398    returnptr = NULL;
399   }
400   else
401    mexErrMsgTxt("Required field 'T2' was not found in the element data structure");
402   tmpmxptr=mxGetField(ElemData,0,"R1");
403   if(tmpmxptr) {
404    pr1 = mxGetPr(tmpmxptr);
405    returnptr = NULL;
406   }
407   else
408    mexErrMsgTxt("Required field 'R1' was not found in the element data structure");
409   tmpmxptr=mxGetField(ElemData,0,"R2");
410   if(tmpmxptr) {
411    pr2 = mxGetPr(tmpmxptr);
412    returnptr = NULL;
413   }
414   else
415    mexErrMsgTxt("Required field 'R2' was not found in the element data structure");
416  }
417  break;
418  case MAKE_LOCAL_COPY: { /* Find field numbers first
419                             Save a list of field number in an array
420                             and make returnptr point to that array */
421   NewFieldNumbers = (int*)mxCalloc(NUM_FIELDS_2_REMEMBER,sizeof(int));
422   fnum = mxGetFieldNumber(ElemData,"Length");
423   if(fnum<0)
424    mexErrMsgTxt("Required field 'Length' was not found in the element data structure"); 
425   NewFieldNumbers[0] = fnum;
426   fnum = mxGetFieldNumber(ElemData,"Order");
427   if(fnum<0)
428    mexErrMsgTxt("Required field 'Order' was not found in the element data structure");
429   NewFieldNumbers[1] = fnum;
430   fnum = mxGetFieldNumber(ElemData,"NumOfVar");
431   if(fnum<0)
432    mexErrMsgTxt("Required field 'NumOfVar' was not found in the element data structure");
433   NewFieldNumbers[2] = fnum;
434   fnum = mxGetFieldNumber(ElemData,"MapX");
435   if(fnum<0)
436    mexErrMsgTxt("Required field 'MapX' was not found in the element data structure");
437   NewFieldNumbers[3] = fnum;
438   fnum = mxGetFieldNumber(ElemData,"NumOfMonoX");
439   if(fnum<0)
440    mexErrMsgTxt("Required field 'NumOfMonoX' was not found in the element data structure");
441   NewFieldNumbers[4] = fnum;
442   fnum = mxGetFieldNumber(ElemData,"MapPx");
443   if(fnum<0)
444    mexErrMsgTxt("Required field 'MapPx' was not found in the element data structure");
445   NewFieldNumbers[5] = fnum;
446   fnum = mxGetFieldNumber(ElemData,"NumOfMonoPx");
447   if(fnum<0)
448    mexErrMsgTxt("Required field 'NumOfMonoPx' was not found in the element data structure");
449   NewFieldNumbers[6] = fnum;
450   fnum = mxGetFieldNumber(ElemData,"MapY");
451   if(fnum<0)
452    mexErrMsgTxt("Required field 'MapY' was not found in the element data structure");
453   NewFieldNumbers[7] = fnum;
454   fnum = mxGetFieldNumber(ElemData,"NumOfMonoY");
455   if(fnum<0)
456    mexErrMsgTxt("Required field 'NumOfMonoY' was not found in the element data structure");
457   NewFieldNumbers[8] = fnum;
458   fnum = mxGetFieldNumber(ElemData,"MapPy");
459   if(fnum<0)
460    mexErrMsgTxt("Required field 'MapPy' was not found in the element data structure");
461   NewFieldNumbers[9] = fnum;
462   fnum = mxGetFieldNumber(ElemData,"NumOfMonoPy");
463   if(fnum<0)
464    mexErrMsgTxt("Required field 'NumOfMonoPy' was not found in the element data structure");
465   NewFieldNumbers[10] = fnum;
466   fnum = mxGetFieldNumber(ElemData,"MapDelta");
467   if(fnum<0)
468    mexErrMsgTxt("Required field 'MapDelta' was not found in the element data structure");
469   NewFieldNumbers[11] = fnum;
470   fnum = mxGetFieldNumber(ElemData,"NumOfMonoDelta");
471   if(fnum<0)
472    mexErrMsgTxt("Required field 'NumOfMonoDelta' was not found in the element data structure");
473   NewFieldNumbers[12] = fnum;
474   fnum = mxGetFieldNumber(ElemData,"MapT");
475   if(fnum<0)
476    mexErrMsgTxt("Required field 'MapT' was not found in the element data structure");
477   NewFieldNumbers[13] = fnum;
478   fnum = mxGetFieldNumber(ElemData,"NumOfMonoT");
479   if(fnum<0)
480    mexErrMsgTxt("Required field 'NumOfMonoT' was not found in the element data structure");
481   NewFieldNumbers[14] = fnum;
482   fnum = mxGetFieldNumber(ElemData,"T1");
483   if(fnum<0)
484    mexErrMsgTxt("Required field 'T1' was not found in the element data structure");
485   NewFieldNumbers[15] = fnum;
486   fnum = mxGetFieldNumber(ElemData,"T2");
487   if(fnum<0)
488    mexErrMsgTxt("Required field 'T2' was not found in the element data structure");
489   NewFieldNumbers[16] = fnum;
490   fnum = mxGetFieldNumber(ElemData,"R1");
491   if(fnum<0)
492    mexErrMsgTxt("Required field 'R1' was not found in the element data structure");
493   NewFieldNumbers[17] = fnum;
494   fnum = mxGetFieldNumber(ElemData,"R2");
495   if(fnum<0)
496    mexErrMsgTxt("Required field 'R2' was not found in the element data structure");
497   NewFieldNumbers[18] = fnum;
498
499   le = mxGetScalar(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[0]));
500   no = (int)mxGetScalar(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[1]));
501   nv = (int)mxGetScalar(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[2]));
502   x = mxGetPr(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[3]));
503   nx = (int)mxGetScalar(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[4]));
504   px = mxGetPr(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[5]));
505   npx = (int)mxGetScalar(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[6]));
506   y = mxGetPr(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[7]));
507   ny = (int)mxGetScalar(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[8]));
508   py = mxGetPr(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[9]));
509   npy = (int)mxGetScalar(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[10]));
510   delta = mxGetPr(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[11]));
511   ndelta = (int)mxGetScalar(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[12]));
512   t = mxGetPr(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[13]));
513   nt = (int)mxGetScalar(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[14]));
514   pt1 = mxGetPr(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[15]));
515   pt2 = mxGetPr(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[16]));
516   pr1 = mxGetPr(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[17]));
517   pr2 = mxGetPr(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[18]));
518
519   returnptr = NewFieldNumbers;
520  }
521  break;
522  case USE_LOCAL_COPY: { /* Get fields from MATLAB using field numbers
523                            The second argument ponter to the array of field
524                            numbers is previously created with
525                            QuadLinPass( ..., MAKE_LOCAL_COPY) */
526   le = mxGetScalar(mxGetFieldByNumber(ElemData,0,FieldNumbers[0]));
527   no = (int)mxGetScalar(mxGetFieldByNumber(ElemData,0,FieldNumbers[1]));
528   nv = (int)mxGetScalar(mxGetFieldByNumber(ElemData,0,FieldNumbers[2]));
529   x = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[3]));
530   nx = (int)mxGetScalar(mxGetFieldByNumber(ElemData,0,FieldNumbers[4]));
531   px = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[5]));
532   npx = (int)mxGetScalar(mxGetFieldByNumber(ElemData,0,FieldNumbers[6]));
533   y = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[7]));
534   ny = (int)mxGetScalar(mxGetFieldByNumber(ElemData,0,FieldNumbers[8]));
535   py = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[9]));
536   npy = (int)mxGetScalar(mxGetFieldByNumber(ElemData,0,FieldNumbers[10]));
537   delta = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[11]));
538   ndelta = (int)mxGetScalar(mxGetFieldByNumber(ElemData,0,FieldNumbers[12]));
539   t = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[13]));
540   nt = (int)mxGetScalar(mxGetFieldByNumber(ElemData,0,FieldNumbers[14]));
541   pt1 = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[15]));
542   pt2 = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[16]));
543   pr1 = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[17]));
544   pr2 = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[18]));
545   returnptr = FieldNumbers;
546  }
547  break;
548  default: {
549   mexErrMsgTxt("No match for calling mode in function TaylormapPass3\n");
550  }
551 }
552 TaylormapPass3(r_in, le, no, nv, x, nx, px, npx, y, ny, py, npy, delta, ndelta, t, nt,
553                    pt1, pt2, pr1, pr2, num_particles);
554 return(returnptr);
555}
556
557void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
558{
559 int m,n; double *r_in;
560 mxArray *tmpmxptr;
561 double le;
562 int no, nv;
563 double *pr1, *pr2, *pt1, *pt2;
564 double *x, *px, *y, *py, *t, *delta;
565 int nx, npx, ny, npy, nt, ndelta;
566
567 /* ALLOCATE memory for the output array of the same size as the input  */
568 m = mxGetM(prhs[1]);
569 n = mxGetN(prhs[1]);
570 if(m!=6) 
571  mexErrMsgTxt("Second argument must be a 6 x N matrix");
572 tmpmxptr=mxGetField(prhs[0],0,"Length");
573 if(tmpmxptr)
574   le = mxGetScalar(tmpmxptr);
575 else
576   mexErrMsgTxt("Required field 'Length' was not found in the element data structure"); 
577 tmpmxptr=mxGetField(prhs[0],0,"Order");
578 if(tmpmxptr)
579   no = (int)mxGetScalar(tmpmxptr);
580 else
581   mexErrMsgTxt("Required field 'Order' was not found in the element data structure");
582 tmpmxptr=mxGetField(prhs[0],0,"NumOfVariables");
583 if(tmpmxptr)
584   nv = (int)mxGetScalar(tmpmxptr);
585 else
586   mexErrMsgTxt("Required field 'NumOfVariables' was not found in the element data structure");
587 tmpmxptr=mxGetField(prhs[0],0,"MapX");
588 if(tmpmxptr)
589   x = mxGetPr(tmpmxptr);
590 else
591   mexErrMsgTxt("Required field 'MapX' was not found in the element data structure");
592 tmpmxptr=mxGetField(prhs[0],0,"NumOfMonoX");
593 if(tmpmxptr)
594   nx = (int)mxGetScalar(tmpmxptr);
595 else
596   mexErrMsgTxt("Required field 'NumOfMonoX' was not found in the element data structure");
597 tmpmxptr=mxGetField(prhs[0],0,"MapPx");
598 if(tmpmxptr)
599   px = mxGetPr(tmpmxptr);
600 else
601   mexErrMsgTxt("Required field 'MapPx' was not found in the element data structure");
602 tmpmxptr=mxGetField(prhs[0],0,"NumOfMonoPx");
603 if(tmpmxptr)
604   npx = (int)mxGetScalar(tmpmxptr);
605 else
606   mexErrMsgTxt("Required field 'NumOfMonoPx' was not found in the element data structure");
607 tmpmxptr=mxGetField(prhs[0],0,"MapY");
608 if(tmpmxptr)
609   y = mxGetPr(tmpmxptr);
610 else
611   mexErrMsgTxt("Required field 'MapY' was not found in the element data structure");
612 tmpmxptr=mxGetField(prhs[0],0,"NumOfMonoY");
613 if(tmpmxptr)
614   ny = (int)mxGetScalar(tmpmxptr);
615 else
616   mexErrMsgTxt("Required field 'NumOfMonoY' was not found in the element data structure");
617 tmpmxptr=mxGetField(prhs[0],0,"MapPy");
618 if(tmpmxptr)
619   py = mxGetPr(tmpmxptr);
620 else
621   mexErrMsgTxt("Required field 'MapPy' was not found in the element data structure");
622 tmpmxptr=mxGetField(prhs[0],0,"NumOfMonoPy");
623 if(tmpmxptr)
624   npy = (int)mxGetScalar(tmpmxptr);
625 else
626   mexErrMsgTxt("Required field 'NumOfMonoPy' was not found in the element data structure");
627 tmpmxptr=mxGetField(prhs[0],0,"MapDelta");
628 if(tmpmxptr)
629   delta = mxGetPr(tmpmxptr);
630 else
631   mexErrMsgTxt("Required field 'MapDelta' was not found in the element data structure");
632 tmpmxptr=mxGetField(prhs[0],0,"NumOfMonoDelta");
633 if(tmpmxptr)
634   ndelta = (int)mxGetScalar(tmpmxptr);
635 else
636   mexErrMsgTxt("Required field 'NumOfMonoDelta' was not found in the element data structure");
637 tmpmxptr=mxGetField(prhs[0],0,"MapT");
638 if(tmpmxptr)
639   t = mxGetPr(tmpmxptr);
640 else
641   mexErrMsgTxt("Required field 'MapT' was not found in the element data structure");
642 tmpmxptr=mxGetField(prhs[0],0,"NumOfMonoT");
643 if(tmpmxptr)
644   nt = (int)mxGetScalar(tmpmxptr);
645 else
646   mexErrMsgTxt("Required field 'NumOfMonoT' was not found in the element data structure");
647 tmpmxptr=mxGetField(prhs[0],0,"R1");
648 if(tmpmxptr)
649   pr1 = mxGetPr(tmpmxptr);
650 else
651   mexErrMsgTxt("Required field 'R1' was not found in the element data structure");
652 tmpmxptr=mxGetField(prhs[0],0,"R2");
653 if(tmpmxptr)
654   pr2 = mxGetPr(tmpmxptr);
655 else
656   mexErrMsgTxt("Required field 'R2' was not found in the element data structure");
657 tmpmxptr=mxGetField(prhs[0],0,"T1");
658 if(tmpmxptr)
659   pt1 = mxGetPr(tmpmxptr);
660 else
661   mexErrMsgTxt("Required field 'T1' was not found in the element data structure");
662 tmpmxptr=mxGetField(prhs[0],0,"T2");
663 if(tmpmxptr)
664   pt2 = mxGetPr(tmpmxptr);
665 else
666   mexErrMsgTxt("Required field 'T2' was not found in the element data structure");
667 plhs[0] = mxDuplicateArray(prhs[1]);
668 r_in = mxGetPr(plhs[0]);
669 TaylormapPass3(r_in, le, no, nv, x, nx, px, npx, y, ny, py, npy, delta, ndelta, t, nt,
670                    pt1, pt2, pr1, pr2, n);
671}
Note: See TracBrowser for help on using the repository browser.