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