source: MML/trunk/at/simulator/element/CavityPass.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: 6.4 KB
Line 
1/* CavityPass.c
2   Accelerator Toolbox
3   Revision 3/10/04
4   A.Terebilo terebilo@ssrl.slac.stanford.edu
5*/
6
7#include "mex.h"
8#include "elempass.h"
9#include <math.h>
10#define TWOPI  6.28318530717959
11#define C0      2.99792458e8
12
13
14void CavityPass(double *r_in, double le, double nv, double freq, double lag, int num_particles)
15/* le - physical length
16   nv - peak voltage (V) normalized to the design enegy (eV)
17   r is a 6-by-N matrix of initial conditions reshaped into
18   1-d array of 6*N elements
19*/
20{       int c, c6;
21        double halflength , p_norm, NormL;     
22        if(le == 0)
23                {       for(c = 0;c<num_particles;c++)
24                        {       c6 = c*6;
25                            if(!mxIsNaN(r_in[c6]))
26                                    r_in[c6+4] += -nv*sin(TWOPI*freq*(r_in[c6+5]-lag)/C0);
27                        }
28                }
29        else
30                {       halflength = le/2;
31                        for(c = 0;c<num_particles;c++)
32                        {       c6 = c*6;
33                                if(!mxIsNaN(r_in[c6])) 
34                                {   p_norm = 1/(1+r_in[c6+4]);                         
35                                    NormL  = halflength*p_norm;
36                                    /* Prropagate through a drift equal to half cavity length */
37                                    r_in[c6+0]+= NormL*r_in[c6+1];
38                                r_in[c6+2]+= NormL*r_in[c6+3];
39                                r_in[c6+5]+= NormL*p_norm*(r_in[c6+1]*r_in[c6+1]+r_in[c6+3]*r_in[c6+3])/2;
40                               
41                                    /* Longitudinal momentum kick */
42                                    r_in[c6+4] += -nv*sin(TWOPI*freq*(r_in[c6+5]-lag)/C0);
43                                    p_norm = 1/(1+r_in[c6+4]);                         
44                                    NormL  = halflength*p_norm;
45                                    /* Prropagate through a drift equal to half cavity length */
46                                    r_in[c6+0]+= NormL*r_in[c6+1];
47                                r_in[c6+2]+= NormL*r_in[c6+3];
48                                r_in[c6+5]+= NormL*p_norm*(r_in[c6+1]*r_in[c6+1]+r_in[c6+3]*r_in[c6+3])/2;
49                }
50                        }
51                }
52
53}
54
55
56
57
58ExportMode int* passFunction(const mxArray *ElemData,int *FieldNumbers,
59                                double *r_in, int num_particles, int mode)
60
61
62#define NUM_FIELDS_2_REMEMBER 5
63
64{       double le, volt, freq, energy, lag;
65        int *returnptr;
66        int *NewFieldNumbers, fnum;
67
68   
69        switch(mode)
70                {       case NO_LOCAL_COPY:     /* Obsolete in AT1.3  */
71                                {       
72                                   
73                                }       break; 
74                       
75                        case MAKE_LOCAL_COPY:   /* Find field numbers first
76                                                                            Save a list of field number in an array
77                                                                                 and make returnptr point to that array
78                                                                    */
79                                {       
80                                        NewFieldNumbers = (int*)mxCalloc(NUM_FIELDS_2_REMEMBER,sizeof(int));
81                                       
82                                        fnum = mxGetFieldNumber(ElemData,"Length");
83                                        if(fnum<0) 
84                                            mexErrMsgTxt("Required field 'Length' was not found in the element data structure"); 
85                                        NewFieldNumbers[0] = fnum;
86                                       
87                                        fnum = mxGetFieldNumber(ElemData,"Voltage");
88                                        if(fnum<0) 
89                                            mexErrMsgTxt("Required field 'Voltage' was not found in the element data structure"); 
90                                        NewFieldNumbers[1] = fnum;
91                                       
92                                        fnum = mxGetFieldNumber(ElemData,"Energy");
93                                        if(fnum<0) 
94                                            mexErrMsgTxt("Required field 'Energy' was not found in the element data structure"); 
95                                        NewFieldNumbers[2] = fnum;
96                                       
97                                       
98                                        fnum = mxGetFieldNumber(ElemData,"Frequency");
99                                        if(fnum<0) 
100                                            mexErrMsgTxt("Required field 'Frequency' was not found in the element data structure"); 
101                                        NewFieldNumbers[3] = fnum;
102                   
103                    /* Optional field TimeLag */
104                    fnum = mxGetFieldNumber(ElemData,"TimeLag");
105                    NewFieldNumbers[4] = fnum;
106                                        if(fnum<0) 
107                                            lag = 0;
108                    else
109                        lag = mxGetScalar(mxGetFieldByNumber(ElemData,0,fnum));
110                                       
111                                       
112                                        le = mxGetScalar(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[0]));
113                                        volt = mxGetScalar(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[1]));
114                                        energy = mxGetScalar(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[2]));
115                                        freq = mxGetScalar(mxGetFieldByNumber(ElemData,0,NewFieldNumbers[3]));
116
117                                        returnptr = NewFieldNumbers;
118                                }       break;
119
120                        case    USE_LOCAL_COPY: /* Get fields from MATLAB using field numbers
121                                                                                 The second argument ponter to the array of field
122                                                                                 numbers is previously created with
123                                                                                 QuadLinPass( ..., MAKE_LOCAL_COPY)
124                                                                    */
125                                                                                       
126                                {       le = mxGetScalar(mxGetFieldByNumber(ElemData,0,FieldNumbers[0]));
127                                        volt = mxGetScalar(mxGetFieldByNumber(ElemData,0,FieldNumbers[1]));
128                                        energy = mxGetScalar(mxGetFieldByNumber(ElemData,0,FieldNumbers[2]));
129                                        freq = mxGetScalar(mxGetFieldByNumber(ElemData,0,FieldNumbers[3]));
130                    /* Optional field TimeLag */
131                    if(FieldNumbers[4]<0) 
132                                            lag = 0;
133                                        else
134                        lag = mxGetScalar(mxGetFieldByNumber(ElemData,0,FieldNumbers[4]));
135                   
136                                        returnptr = FieldNumbers;
137                                }       break;
138
139                        default:
140                                {       mexErrMsgTxt("No match found for calling mode in function CavityPass\n");
141                                }
142        }
143
144       
145       
146        CavityPass(r_in, le, volt/energy, freq, lag, num_particles);
147        return(returnptr);
148}
149
150
151
152
153
154
155
156
157
158void mexFunction(       int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
159{       double volt,freq, energy, lag;
160        int m,n;
161        double *r_in, le;   
162        mxArray *tmpmxptr;
163       
164        if(nrhs)
165        {
166
167        /* ALLOCATE memory for the output array of the same size as the input */
168        m = mxGetM(prhs[1]);
169        n = mxGetN(prhs[1]);
170        if(m!=6) 
171                mexErrMsgTxt("Second argument must be a 6 x N matrix");
172
173       
174        tmpmxptr=mxGetField(prhs[0],0,"Length");
175        if(tmpmxptr)
176                le = mxGetScalar(tmpmxptr);
177        else
178                mexErrMsgTxt("Required field 'Length' was not found in the element data structure"); 
179                               
180                                           
181        tmpmxptr=mxGetField(prhs[0],0,"Voltage");
182        if(tmpmxptr)
183                volt = mxGetScalar(tmpmxptr);
184        else
185                mexErrMsgTxt("Required field 'Voltage' was not found in the element data structure");
186               
187    tmpmxptr=mxGetField(prhs[0],0,"Energy");
188        if(tmpmxptr)
189                energy = mxGetScalar(tmpmxptr);
190        else
191                mexErrMsgTxt("Required field 'Energy' was not found in the element data structure");
192                                           
193        tmpmxptr=mxGetField(prhs[0],0,"Frequency");
194        if(tmpmxptr)
195                freq = mxGetScalar(tmpmxptr);
196        else
197                mexErrMsgTxt("Required field 'Frequency' was not found in the element data structure"); 
198   
199    tmpmxptr=mxGetField(prhs[0],0,"TimeLag");
200        if(tmpmxptr)
201                lag = mxGetScalar(tmpmxptr);
202        else
203                lag = 0; 
204
205
206    plhs[0] = mxDuplicateArray(prhs[1]);
207    r_in = mxGetPr(plhs[0]);
208        CavityPass(r_in, le, volt/energy, freq, lag, n);
209    }
210    else
211    {   /* return list of required fields */
212            plhs[0] = mxCreateCellMatrix(4,1);
213            mxSetCell(plhs[0],0,mxCreateString("Length"));
214            mxSetCell(plhs[0],1,mxCreateString("Voltage"));
215            mxSetCell(plhs[0],2,mxCreateString("Energy"));
216            mxSetCell(plhs[0],3,mxCreateString("Frequency"));
217            if(nlhs>1) /* optional fields */ 
218            {   plhs[1] = mxCreateCellMatrix(1,1); 
219            mxSetCell(plhs[1],0,mxCreateString("TimeLag"));
220            }
221        }
222
223}
Note: See TracBrowser for help on using the repository browser.