source: MML/trunk/at/simulator/element/user/GWigSymplecticPass.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: 14.5 KB
Line 
1/* GWigSymplecticPass.c for
2   Accelerator Toolbox
3*/
4
5/*
6 *---------------------------------------------------------------------------
7 * Modification Log:
8 * -----------------
9 * .02  2003-06-18     J. Li, jing@fel.duke.edu
10 *                              Cleanup the code
11 *
12 * .01  2003-04-20     YK Wu, wu@fel.duke.edu
13 *                              GWiggler interface
14 *
15 *---------------------------------------------------------------------------
16 *  Accelerator Physics Group, Duke FEL Lab, www.fel.duke.edu 
17 */
18
19#include "mex.h"
20#include "elempass.h"
21#include "../atlalib.c"
22#include <stdlib.h>
23#include <math.h>
24#include "gwig.h"
25
26#define GWIG
27
28#include "gwig.c"
29
30/******************************************************************************/
31/* PHYSICS SECTION ************************************************************/
32
33void GWigInit(struct gwig *Wig, double Ltot, double Lw, double Bmax,
34              int Nstep, int Nmeth, int NHharm, int NVharm,
35              double *pBy, double *pBx, double *T1, double *T2, 
36              double *R1, double *R2)
37{
38  double *tmppr;
39  double design_energy;
40  int    i;
41  double kw;
42  mxArray *E0Field;
43  const mxArray *GLOBVALPTR = mexGetVariablePtr("GLOBVAL","global");
44
45  /* Get the design energy for normalization from GLOBVAL.E0 in
46   * global workspace
47   */
48  if(GLOBVALPTR == NULL)
49    mexPrintf("GLOBVALPTR = NULL\n");
50  if(GLOBVALPTR != NULL){
51    E0Field = mxGetField(GLOBVALPTR,0,"E0");
52    if(E0Field != NULL)
53      design_energy = mxGetScalar(E0Field)*1e-9; /* convert to GeV */
54    else
55      mexErrMsgTxt("Field 'E0' is not defined in GLOBVAL structure");
56  }
57  else
58    mexErrMsgTxt("global variable GLOBVAL does not exist");
59
60  Wig->E0 = design_energy;
61  Wig->Pmethod = Nmeth;
62  Wig->PN = Nstep;
63  Wig->Nw = (int)(Ltot / Lw);
64  Wig->NHharm = NHharm;
65  Wig->NVharm = NVharm;
66  Wig->PB0 = Bmax;
67  Wig->Lw  = Lw;
68
69
70  kw = 2.0e0*PI/(Wig->Lw);
71  Wig->Zw = 0.0;
72  Wig->Aw = 0.0;
73  tmppr = pBy;
74  for (i = 0; i < NHharm; i++){
75    tmppr++;
76    Wig->HCw[i] = 0.0;
77    Wig->HCw_raw[i] = *tmppr;
78
79    tmppr++;
80    Wig->Hkx[i]     = (*tmppr) * kw;
81
82    tmppr++;
83    Wig->Hky[i]     = (*tmppr) * kw;
84
85    tmppr++;
86    Wig->Hkz[i]     = (*tmppr) * kw;
87
88    tmppr++;
89    Wig->Htz[i]     =  *tmppr;
90
91    tmppr++;
92  }
93
94  tmppr = pBx;
95  for (i = 0; i < NVharm; i++){
96    tmppr++;
97    Wig->VCw[i] = 0.0;
98    Wig->VCw_raw[i] = *tmppr;
99
100    tmppr++;
101    Wig->Vkx[i]     = (*tmppr) * kw;
102
103    tmppr++;
104    Wig->Vky[i]     = (*tmppr) * kw;
105
106    tmppr++;
107    Wig->Vkz[i]     = (*tmppr) * kw;
108
109    tmppr++;
110    Wig->Vtz[i]     =  *tmppr;
111
112    tmppr++;
113  }
114 
115  for (i = NHharm ; i< WHmax; i++) {
116    Wig->HCw[i] = 0.0;
117    Wig->HCw_raw[i] = 0.0;
118    Wig->Hkx[i] = 0.0;
119    Wig->Hky[i] = 0.0;
120    Wig->Hkz[i] = 0.0;
121    Wig->Htz[i] = 0.0;
122  }
123  for (i = NVharm ; i< WHmax; i++) {
124    Wig->VCw[i] = 0.0;
125    Wig->VCw_raw[i] = 0.0;
126    Wig->Vkx[i] = 0.0;
127    Wig->Vky[i] = 0.0;
128    Wig->Vkz[i] = 0.0;
129    Wig->Vtz[i] = 0.0;
130  }
131}
132
133#define second 2
134#define fourth 4
135void GWigSymplecticPass(double *r, double Ltot, double Lw, double Bmax,
136                        int Nstep, int Nmeth, int NHharm, int NVharm,
137                        double *pBy, double *pBx, double *T1, double *T2, 
138                        double *R1, double *R2, int num_particles)
139{       
140
141  int c, i;
142  double *r6;
143  double *tmppr;
144  struct gwig Wig;
145
146  GWigInit(&Wig, Ltot, Lw, Bmax, Nstep, Nmeth,NHharm,NVharm,pBy,pBx,T1,T2,R1,R2);
147
148  for(c = 0;c<num_particles;c++){       
149    r6 = r+c*6; 
150   
151    if(!mxIsNaN(r6[0]) & mxIsFinite(r6[4]))
152    {
153    switch (Nmeth) {
154      case second :
155        GWigPass_2nd(&Wig, r6);
156        break;
157      case fourth:
158        GWigPass_4th(&Wig, r6);
159        break;
160      default:
161        printf("Invalid method ...\n");
162        break;
163    }
164    } 
165  }
166
167}
168
169/********** END PHYSICS SECTION ***********************************************/
170/******************************************************************************/
171
172ExportMode int* passFunction(const mxArray *ElemData, int *FieldNumbers,
173                             double *r_in, int num_particles, int mode)
174
175#define NUM_FIELDS_2_REMEMBER 11
176
177{       
178  double *pr1, *pr2, *pt1, *pt2;
179  int m, n;
180  double *pBy, *pBx;
181  double Ltot, Lw, Bmax; 
182  int Nstep, Nmeth;
183  int NHharm, NVharm;
184  mxArray *tmpmxptr;
185  int *returnptr,fnum;
186  int *NewFieldNumbers;
187
188  switch(mode)
189    {   
190    case NO_LOCAL_COPY: /* Get fields by names from MATLAB workspace  */
191      { 
192        tmpmxptr = mxGetField(ElemData,0,"Length");
193        if(tmpmxptr)
194          Ltot = mxGetScalar(tmpmxptr);
195        else
196          mexErrMsgTxt("Required field 'Length' was not found in the element data structure"); 
197       
198        tmpmxptr = mxGetField(ElemData,0,"Lw");
199        if(tmpmxptr)
200          Lw = mxGetScalar(tmpmxptr);
201        else
202          mexErrMsgTxt("Required field 'Lw' was not found in the element data structure"); 
203       
204        tmpmxptr = mxGetField(ElemData,0,"Bmax");
205        if(tmpmxptr)
206          Bmax = mxGetScalar(tmpmxptr);
207        else
208          mexErrMsgTxt("Required field 'Bmax' was not found in the element data structure"); 
209       
210        tmpmxptr = mxGetField(ElemData,0,"Nstep");
211        if(tmpmxptr)
212          Nstep = (int)mxGetScalar(tmpmxptr);
213        else
214          mexErrMsgTxt("Required field 'Nstep' was not found in the element data structure"); 
215 
216        tmpmxptr = mxGetField(ElemData,0,"Nmeth");
217        if(tmpmxptr)
218          Nmeth = (int)mxGetScalar(tmpmxptr);
219        else
220          mexErrMsgTxt("Required field 'Nmeth' was not found in the element data structure"); 
221
222        tmpmxptr = mxGetField(ElemData,0,"NHharm");
223        if(tmpmxptr)
224          NHharm = (int)mxGetScalar(tmpmxptr);
225        else
226          mexErrMsgTxt("Required field 'NHharm' was not found in the element data structure"); 
227       
228        tmpmxptr = mxGetField(ElemData,0,"NVharm");
229        if(tmpmxptr)
230          NVharm = (int)mxGetScalar(tmpmxptr);
231        else
232          mexErrMsgTxt("Required field 'NVharm' was not found in the element data structure"); 
233
234        tmpmxptr = mxGetField(ElemData,0,"By");
235        if(tmpmxptr)
236          pBy = mxGetPr(tmpmxptr);
237        else
238          mexErrMsgTxt("Required field 'By' was not found in the element data structure"); 
239
240        tmpmxptr = mxGetField(ElemData,0,"Bx");
241        if(tmpmxptr)
242          pBx = mxGetPr(tmpmxptr);
243        else
244          mexErrMsgTxt("Required field 'Bx' was not found in the element data structure"); 
245 
246        tmpmxptr = mxGetField(ElemData,0,"R1");
247        if(tmpmxptr)
248          pr1 = mxGetPr(tmpmxptr);
249        else
250          mexErrMsgTxt("Required field 'R1' was not found in the element data structure"); 
251       
252        tmpmxptr = mxGetField(ElemData,0,"R2");
253        if(tmpmxptr)
254          pr2 = mxGetPr(tmpmxptr);
255        else
256          mexErrMsgTxt("Required field 'R2' was not found in the element data structure"); 
257       
258        tmpmxptr = mxGetField(ElemData,0,"T1");
259        if(tmpmxptr)
260          pt1 = mxGetPr(tmpmxptr);
261        else
262          mexErrMsgTxt("Required field 'T1' was not found in the element data structure"); 
263       
264        tmpmxptr = mxGetField(ElemData,0,"T2");
265        if(tmpmxptr)
266          pt2 = mxGetPr(tmpmxptr);
267        else
268          mexErrMsgTxt("Required field 'T2' was not found in the element data structure");         
269
270        returnptr = NULL;
271       
272      } break; 
273    case MAKE_LOCAL_COPY:       /* Find field numbers first
274                                 * Save a list of field number in an array
275                                 *  and make returnptr point to that array
276                                 */
277      {                                         
278        NewFieldNumbers = (int*)mxCalloc(NUM_FIELDS_2_REMEMBER,sizeof(int));
279                                       
280        fnum = mxGetFieldNumber(ElemData,"Length");
281        if(fnum<0) 
282          mexErrMsgTxt("Required field 'Length' was not found in the element data structure"); 
283        NewFieldNumbers[0] = fnum;
284                                       
285        fnum = mxGetFieldNumber(ElemData,"Lw");
286        if(fnum<0) 
287          mexErrMsgTxt("Required field 'Lw' was not found in the element data structure"); 
288        NewFieldNumbers[1] = fnum;
289                                       
290
291        fnum = mxGetFieldNumber(ElemData,"Bmax");
292        if(fnum<0) 
293          mexErrMsgTxt("Required field 'Bmax' was not found in the element data structure"); 
294        NewFieldNumbers[2] = fnum;
295                                       
296
297        fnum = mxGetFieldNumber(ElemData,"Nstep");
298        if(fnum<0) 
299          mexErrMsgTxt("Required field 'Nstep' was not found in the element data structure"); 
300        NewFieldNumbers[3] = fnum;
301                                       
302                                       
303        fnum = mxGetFieldNumber(ElemData,"Nmeth");
304        if(fnum<0) 
305          mexErrMsgTxt("Required field 'Nmeth' was not found in the element data structure"); 
306        NewFieldNumbers[4] = fnum;
307                                       
308        fnum = mxGetFieldNumber(ElemData,"NHharm");
309        if(fnum<0) 
310          mexErrMsgTxt("Required field 'NHharm' was not found in the element data structure"); 
311        NewFieldNumbers[5] = fnum;
312
313        fnum = mxGetFieldNumber(ElemData,"NVharm");
314        if(fnum<0) 
315          mexErrMsgTxt("Required field 'NVharm' was not found in the element data structure"); 
316        NewFieldNumbers[6] = fnum;       
317         
318        fnum = mxGetFieldNumber(ElemData,"By");
319        if(fnum<0) 
320          mexErrMsgTxt("Required field 'By' was not found in the element data structure"); 
321        NewFieldNumbers[7] = fnum;
322                       
323        fnum = mxGetFieldNumber(ElemData,"Bx");
324        if(fnum<0) 
325          mexErrMsgTxt("Required field 'Bx' was not found in the element data structure"); 
326        NewFieldNumbers[8] = fnum;
327                       
328        fnum = mxGetFieldNumber(ElemData,"R1");
329        if(fnum<0) 
330          mexErrMsgTxt("Required field 'R1' was not found in the element data structure"); 
331        NewFieldNumbers[9] = fnum;
332       
333       
334        fnum = mxGetFieldNumber(ElemData,"R2");
335        if(fnum<0) 
336          mexErrMsgTxt("Required field 'R2' was not found in the element data structure"); 
337        NewFieldNumbers[10] = fnum;
338       
339       
340        fnum = mxGetFieldNumber(ElemData,"T1");
341        if(fnum<0) 
342          mexErrMsgTxt("Required field 'T1' was not found in the element data structure"); 
343        NewFieldNumbers[11] = fnum;
344       
345       
346        fnum = mxGetFieldNumber(ElemData,"T2");
347        if(fnum<0) 
348          mexErrMsgTxt("Required field 'T2' was not found in the element data structure"); 
349        NewFieldNumbers[12] = fnum;
350       
351        Ltot   = mxGetScalar(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[0]));
352        Lw     = mxGetScalar(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[1]));
353        Bmax   = mxGetScalar(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[2]));
354        Nstep  = mxGetScalar(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[3]));
355        Nmeth  = mxGetScalar(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[4]));
356        NHharm = mxGetScalar(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[5]));
357        NVharm = mxGetScalar(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[6]));
358        pBy    = mxGetPr(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[7]));
359        pBx    = mxGetPr(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[8]));
360        pr1    = mxGetPr(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[9]));
361        pr2    = mxGetPr(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[10]));
362        pt1    = mxGetPr(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[11]));
363        pt2    = mxGetPr(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[12]));
364       
365        returnptr = NewFieldNumbers;
366      } break;
367    case        USE_LOCAL_COPY: /* Get fields from MATLAB using field numbers
368                                   The second argument ponter to the array of
369                                   field numbers is previously created with
370                                   QuadLinPass( ..., MAKE_LOCAL_COPY)
371                                */
372      { 
373        Ltot   = mxGetScalar(mxGetFieldByNumber(ElemData,0,FieldNumbers[0]));
374        Lw     = mxGetScalar(mxGetFieldByNumber(ElemData,0,FieldNumbers[1]));
375        Bmax   = mxGetScalar(mxGetFieldByNumber(ElemData,0,FieldNumbers[2]));
376        Nstep  = mxGetScalar(mxGetFieldByNumber(ElemData,0,FieldNumbers[3]));
377        Nmeth  = mxGetScalar(mxGetFieldByNumber(ElemData,0,FieldNumbers[4]));
378        NHharm = mxGetScalar(mxGetFieldByNumber(ElemData,0,FieldNumbers[5]));
379        NVharm = mxGetScalar(mxGetFieldByNumber(ElemData,0,FieldNumbers[6]));
380        pBy    = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[7]));
381        pBx    = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[8]));
382        pr1    = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[9]));
383        pr2    = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[10]));
384        pt1    = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[11]));
385        pt2    = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[12]));
386       
387        returnptr = FieldNumbers;
388      } break;
389    default:
390      { mexErrMsgTxt("No match for calling mode in function QuadLinPass\n");
391      }
392    }
393
394  GWigSymplecticPass(r_in,Ltot,Lw,Bmax,Nstep,Nmeth,NHharm,NVharm,pBy,pBx,pt1,pt2,pr1,pr2,num_particles);
395  return(returnptr);   
396}
397
398
399/********** END WINDOWS DLL GATEWAY SECTION ***************************************/
400
401/********** MATLAB GATEWAY  ***************************************/
402
403void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
404{       
405  int m, n;
406  double *r_in;
407  double *pBy, *pBx;
408  double *pr1, *pr2, *pt1, *pt2 ;
409  double Ltot, Lw, Bmax; 
410  int Nstep, Nmeth;
411  int NHharm, NVharm;
412  mxArray *tmpmxptr;
413       
414  /* ALLOCATE memory for the output array of the same size as the input */
415  m = mxGetM(prhs[1]);
416  n = mxGetN(prhs[1]);
417  if(m!=6) 
418    {mexErrMsgTxt("Second argument must be a 6 x N matrix");}   
419       
420  tmpmxptr = mxGetField(prhs[0],0,"Length");
421  if(tmpmxptr)
422    Ltot = mxGetScalar(tmpmxptr);
423  else
424    mexErrMsgTxt("Required field 'Length' was not found in the element data structure"); 
425       
426  tmpmxptr = mxGetField(prhs[0],0,"Lw");
427  if(tmpmxptr)
428    Lw = mxGetScalar(tmpmxptr);
429  else
430    mexErrMsgTxt("Required field 'Lw' was not found in the element data structure"); 
431 
432  tmpmxptr = mxGetField(prhs[0],0,"Bmax");
433  if(tmpmxptr)
434    Bmax = mxGetScalar(tmpmxptr);
435  else
436    mexErrMsgTxt("Required field 'Bmax' was not found in the element data structure"); 
437 
438  tmpmxptr = mxGetField(prhs[0],0,"Nstep");
439  if(tmpmxptr)
440    Nstep = (int)mxGetScalar(tmpmxptr);
441  else
442    mexErrMsgTxt("Required field 'Nstep' was not found in the element data structure"); 
443 
444  tmpmxptr = mxGetField(prhs[0],0,"Nmeth");
445  if(tmpmxptr)
446    Nmeth = (int)mxGetScalar(tmpmxptr);
447  else
448    mexErrMsgTxt("Required field 'Nmeth' was not found in the element data structure"); 
449
450  tmpmxptr = mxGetField(prhs[0],0,"NHharm");
451  if(tmpmxptr)
452    NHharm = (int)mxGetScalar(tmpmxptr);
453  else
454    mexErrMsgTxt("Required field 'NHharm' was not found in the element data structure"); 
455 
456  tmpmxptr = mxGetField(prhs[0],0,"NVharm");
457  if(tmpmxptr)
458    NVharm = (int)mxGetScalar(tmpmxptr);
459  else
460    mexErrMsgTxt("Required field 'NVharm' was not found in the element data structure"); 
461 
462  tmpmxptr = mxGetField(prhs[0],0,"By");
463  if(tmpmxptr)
464    pBy = mxGetPr(tmpmxptr);
465  else
466    mexErrMsgTxt("Required field 'By' was not found in the element data structure"); 
467
468  tmpmxptr = mxGetField(prhs[0],0,"Bx");
469  if(tmpmxptr)
470    pBx = mxGetPr(tmpmxptr);
471  else
472    mexErrMsgTxt("Required field 'Bx' was not found in the element data structure"); 
473 
474  tmpmxptr = mxGetField(prhs[0],0,"R1");
475  if(tmpmxptr)
476    pr1 = mxGetPr(tmpmxptr);
477  else
478    mexErrMsgTxt("Required field 'R1' was not found in the element data structure"); 
479 
480  tmpmxptr = mxGetField(prhs[0],0,"R2");
481  if(tmpmxptr)
482    pr2 = mxGetPr(tmpmxptr);
483  else
484    mexErrMsgTxt("Required field 'R2' was not found in the element data structure"); 
485 
486  tmpmxptr = mxGetField(prhs[0],0,"T1");
487  if(tmpmxptr)
488    pt1 = mxGetPr(tmpmxptr);
489  else
490    mexErrMsgTxt("Required field 'T1' was not found in the element data structure"); 
491 
492  tmpmxptr = mxGetField(prhs[0],0,"T2");
493  if(tmpmxptr)
494    pt2 = mxGetPr(tmpmxptr);
495  else
496    mexErrMsgTxt("Required field 'T2' was not found in the element data structure");       
497
498  plhs[0] = mxDuplicateArray(prhs[1]);
499  r_in = mxGetPr(plhs[0]);
500 
501  GWigSymplecticPass(r_in,Ltot,Lw,Bmax,Nstep,Nmeth,NHharm,NVharm,pBy,pBx,pt1,pt2,pr1,pr2,n);
502}
503
504
Note: See TracBrowser for help on using the repository browser.