Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

StrCorrMPoleSymplectic4RadPass.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: 14.3 KB

Line
1	#include "mex.h"
2	#include<math.h>
3	#include "../atlalib.c"
4	#include "elempass.h"
5
6
7	#define DRIFT1 0.6756035959798286638
8	#define DRIFT2 -0.1756035959798286639
9	#define KICK1 1.351207191959657328
10	#define KICK2 -1.702414383919314656
11
12
13
14
15	#define SQR(X) ((X)*(X))
16
17
18
19	double StrB2perp(double bx, double by,
20	double x, double xpr, double y, double ypr)
21	/* Calculates sqr(\|B x e\|) , where e is a unit vector in the direction of velocity */
22
23	{ double v_norm2;
24	v_norm2 = 1/(1 + SQR(xpr) + SQR(ypr));
25
26	/* components of the normalized velocity vector
27	double ex, ey, ez;
28	ex = xpr;
29	ey = ypr;
30	ez = 1;
31	*/
32
33	return((SQR(by) + SQR(bx) + SQR(bxypr - byxpr) )*v_norm2) ;
34
35	}
36
37
38
39
40	void strthinkickrad(double* r, double* A, double* B, double L, double E0, int max_order)
41
42	/*****************************************************************************
43	Calculate and apply
44	(a) multipole kick
45	(b) momentum kick due to classical radiation
46	to a 6-dimentional phase space vector in a straight element ( quadrupole)
47
48	IMPORTANT !!!
49	The reference coordinate system is straight but the field expansion may still
50	contain dipole terms: PolynomA(1), PolynomB(1) - in MATLAB notation,
51	A[0], B[0] - C,C++ notation
52
53
54	Note: According to US convention the transverse multipole field is written as:
55
56	max_order+1
57	----
58	\ n-1
59	(B + iB )/ B rho = > (ia + b ) (x + iy)
60	y x / n n
61	----
62	n=1
63	is a polynomial in (x,y) with the highest order = MaxOrder
64
65
66	Using different index notation
67
68	max_order
69	----
70	\ n
71	(B + iB )/ B rho = > (iA + B ) (x + iy)
72	y x / n n
73	----
74	n=0
75
76	A,B: i=0 ... max_order
77	[0] - dipole, [1] - quadrupole, [2] - sextupole ...
78	units for A,B[i] = 1/[m]^(i+1)
79	Coeficients are stroed in the PolynomA, PolynomB field of the element
80	structure in MATLAB
81
82	A[i] (C++,C) = PolynomA(i+1) (MATLAB)
83	B[i] (C++,C) = PolynomB(i+1) (MATLAB)
84	i = 0 .. MaxOrder
85
86	******************************************************************************/
87	{ int i;
88	double ReSumTemp;
89	double ReSum = B[max_order];
90	double ImSum = A[max_order];
91	double x ,xpr, y, ypr, p_norm, B2P;
92
93	#define TWOPI 6.28318530717959
94	#define CGAMMA 8.846056192e-05
95
96
97	double CRAD = CGAMMAE0E0E0/(TWOPI1e27); /* [m]/[GeV^3] M.Sands (4.1) */
98
99
100
101
102	/* recursively calculate the local transvrese magnetic field
103	Bx = ImSum, By = ReSum
104	*/
105	for(i=max_order-1;i>=0;i--)
106	{ ReSumTemp = ReSumr[0] - ImSumr[2] + B[i];
107	ImSum = ImSumr[0] + ReSumr[2] + A[i];
108	ReSum = ReSumTemp;
109	}
110
111
112	/* calculate angles from momentums */
113	p_norm = 1/(1+r[4]);
114	x = r[0];
115	xpr = r[1]*p_norm;
116	y = r[2];
117	ypr = r[3]*p_norm;
118
119	/* For instantaneous rate of energy loss due to classical radiation
120	need to calculate \|n x B\|^2, n unit vector in the direction of velocity
121	*/
122	B2P = StrB2perp(ImSum, ReSum , x , xpr, y ,ypr);
123
124
125	r[4] = r[4] - CRAD(1+r[4])(1+r[4])B2P(1 + (SQR(xpr)+SQR(ypr))/2 )*L;
126
127	/* recalculate momentums from angles after losing energy for radiation */
128	p_norm = 1/(1+r[4]);
129	r[1] = xpr/p_norm;
130	r[3] = ypr/p_norm;
131
132
133	r[1] -= L*ReSum;
134	r[3] += L*ImSum;
135
136	}
137
138
139
140	void StrMPoleSymplectic4RadPass(double r, double le, double A, double *B,
141	int max_order, int num_int_steps,
142	double T1, double T2,
143	double R1, double R2,
144	double E0,
145	int num_particles)
146
147
148	{ int c,m;
149	double *r6;
150	double SL, L1, L2, K1, K2;
151	bool useT1, useT2, useR1, useR2;
152	SL = le/num_int_steps;
153	L1 = SL*DRIFT1;
154	L2 = SL*DRIFT2;
155	K1 = SL*KICK1;
156	K2 = SL*KICK2;
157
158
159	if(T1==NULL)
160	useT1=false;
161	else
162	useT1=true;
163
164	if(T2==NULL)
165	useT2=false;
166	else
167	useT2=true;
168
169	if(R1==NULL)
170	useR1=false;
171	else
172	useR1=true;
173
174	if(R2==NULL)
175	useR2=false;
176	else
177	useR2=true;
178	for(c = 0;c<num_particles;c++) /* Loop over particles */
179	{ r6 = r+c*6;
180	if(!mxIsNaN(r6[0]))
181	{
182
183	/* misalignment at entrance */
184	if(useT1)
185	ATaddvv(r6,T1);
186	if(useR1)
187	ATmultmv(r6,R1);
188
189	/* integrator */
190	for(m=0; m < num_int_steps; m++) /* Loop over slices */
191	{ r6 = r+c*6;
192
193	ATdrift6(r6,L1);
194	strthinkickrad(r6, A, B, K1, E0, max_order);
195	ATdrift6(r6,L2);
196	strthinkickrad(r6, A, B, K2, E0, max_order);
197	ATdrift6(r6,L2);
198	strthinkickrad(r6, A, B, K1, E0, max_order);
199	ATdrift6(r6,L1);
200	}
201	/* Misalignment at exit */
202	if(useR2)
203	ATmultmv(r6,R2);
204	if(useT2)
205	ATaddvv(r6,T2);
206
207	}
208	}
209	}
210
211
212
213
214
215
216	ExportMode int* passFunction(const mxArray ElemData, int FieldNumbers,
217	double *r_in, int num_particles, int mode)
218
219	#define NUM_FIELDS_2_REMEMBER 11
220
221
222	{ double A , B;
223	double pr1, pr2, pt1, pt2, *ka;
224	double E0; /* Design energy [eV] */
225
226
227
228	int max_order, num_int_steps;
229	double le;
230	int *returnptr;
231	int *NewFieldNumbers, fnum;
232
233	switch(mode)
234	{ case NO_LOCAL_COPY: /* Obsolete in AT1.3 */
235	{
236
237	} break;
238
239	case MAKE_LOCAL_COPY: /* Find field numbers first
240	Save a list of field number in an array
241	and make returnptr point to that array
242	*/
243	{ /* Allocate memory for integer array of
244	field numbers for faster futurereference
245	*/
246	NewFieldNumbers = (int*)mxCalloc(NUM_FIELDS_2_REMEMBER,sizeof(int));
247
248	/* Populate */
249
250	fnum = mxGetFieldNumber(ElemData,"PolynomA");
251	if(fnum<0)
252	mexErrMsgTxt("Required field 'PolynomA' was not found in the element data structure");
253	NewFieldNumbers[0] = fnum;
254	A = mxGetPr(mxGetFieldByNumber(ElemData,0,fnum));
255
256
257	fnum = mxGetFieldNumber(ElemData,"PolynomB");
258	if(fnum<0)
259	mexErrMsgTxt("Required field 'PolynomB' was not found in the element data structure");
260	NewFieldNumbers[1] = fnum;
261	B = mxGetPr(mxGetFieldByNumber(ElemData,0,fnum));
262
263
264
265	fnum = mxGetFieldNumber(ElemData,"MaxOrder");
266	if(fnum<0)
267	mexErrMsgTxt("Required field 'MaxOrder' was not found in the element data structure");
268	NewFieldNumbers[2] = fnum;
269	max_order = (int)mxGetScalar(mxGetFieldByNumber(ElemData,0,fnum));
270
271
272	fnum = mxGetFieldNumber(ElemData,"NumIntSteps");
273	if(fnum<0)
274	mexErrMsgTxt("Required field 'NumIntSteps' was not found in the element data structure");
275	NewFieldNumbers[3] = fnum;
276	num_int_steps = (int)mxGetScalar(mxGetFieldByNumber(ElemData,0,fnum));
277
278
279	fnum = mxGetFieldNumber(ElemData,"Length");
280	if(fnum<0)
281	mexErrMsgTxt("Required field 'Length' was not found in the element data structure");
282	NewFieldNumbers[4] = fnum;
283	le = mxGetScalar(mxGetFieldByNumber(ElemData,0,fnum));
284
285	fnum = mxGetFieldNumber(ElemData,"Energy");
286	if(fnum<0)
287	mexErrMsgTxt("Required field 'Energy' was not found in the element data structure");
288	NewFieldNumbers[5] = fnum;
289	E0 = mxGetScalar(mxGetFieldByNumber(ElemData,0,fnum));
290
291
292	fnum = mxGetFieldNumber(ElemData,"R1");
293	NewFieldNumbers[6] = fnum;
294	if(fnum<0)
295	pr1 = NULL;
296	else
297	pr1 = mxGetPr(mxGetFieldByNumber(ElemData,0,fnum));
298
299
300	fnum = mxGetFieldNumber(ElemData,"R2");
301	NewFieldNumbers[7] = fnum;
302	if(fnum<0)
303	pr2 = NULL;
304	else
305	pr2 = mxGetPr(mxGetFieldByNumber(ElemData,0,fnum));
306
307
308	fnum = mxGetFieldNumber(ElemData,"T1");
309	NewFieldNumbers[8] = fnum;
310	if(fnum<0)
311	pt1 = NULL;
312	else
313	pt1 = mxGetPr(mxGetFieldByNumber(ElemData,0,fnum));
314
315
316	fnum = mxGetFieldNumber(ElemData,"T2");
317	NewFieldNumbers[9] = fnum;
318	if(fnum<0)
319	pt2 = NULL;
320	else
321	pt2 = mxGetPr(mxGetFieldByNumber(ElemData,0,fnum));
322
323
324	/* Optional: Kick angles, see section below for explanation */
325	/* Kicks from multipole elements can be specified as angles. This handles the
326	case where corrector coils are used in sextupoles and used for orbit
327	correction. */
328	fnum = mxGetFieldNumber(ElemData,"KickAngle");
329	NewFieldNumbers[10] = fnum;
330	if(fnum<0)
331	ka = NULL;
332	else
333	ka = mxGetPr(mxGetFieldByNumber(ElemData,0,fnum));
334
335
336	returnptr = NewFieldNumbers;
337
338	} break;
339
340	case USE_LOCAL_COPY: /* Get fields from MATLAB using field numbers
341	The second argument ponter to the array of field
342	numbers is previously created with
343	QuadLinPass( ..., MAKE_LOCAL_COPY)
344	*/
345
346	{ A = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[0]));
347	B = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[1]));
348	max_order = (int)mxGetScalar(mxGetFieldByNumber(ElemData,0,FieldNumbers[2]));
349	num_int_steps = (int)mxGetScalar(mxGetFieldByNumber(ElemData,0,FieldNumbers[3]));
350	le = mxGetScalar(mxGetFieldByNumber(ElemData,0,FieldNumbers[4]));
351	E0 = mxGetScalar(mxGetFieldByNumber(ElemData,0,FieldNumbers[5]));
352
353
354	/* Optional fields */
355	if(FieldNumbers[6]<0)
356	pr1 = NULL;
357	else
358	pr1 = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[6]));
359
360	if(FieldNumbers[7]<0)
361	pr2 = NULL;
362	else
363	pr2 = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[7]));
364
365
366	if(FieldNumbers[8]<0)
367	pt1 = NULL;
368	else
369	pt1 = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[8]));
370
371	if(FieldNumbers[9]<0)
372	pt2 = NULL;
373	else
374	pt2 = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[9]));
375
376	if(FieldNumbers[10]<0)
377	ka = NULL;
378	else
379	ka = mxGetPr(mxGetFieldByNumber(ElemData,0,FieldNumbers[9]));
380
381	returnptr = FieldNumbers;
382
383	} break;
384	default:
385	{ mexErrMsgTxt("No match for calling mode in function StrMPoleSymplectic4RadPass\n");
386	}
387	}
388
389
390	if(ka!=NULL)
391	{
392	/* Positive angle must correspond to -ve B field since +ve B field corresponds
393	to a bend with the same curvature as the bend magnets.
394	*/
395	B[0] -= sin(ka[0])/le;
396	A[0] += sin(ka[1])/le;
397	}
398
399	StrMPoleSymplectic4RadPass(r_in, le, A, B, max_order, num_int_steps,pt1, pt2, pr1, pr2, E0, num_particles);
400
401	if(ka!=NULL)
402	{
403	B[0] += sin(ka[0])/le;
404	A[0] -= sin(ka[1])/le;
405	}
406
407	return(returnptr);
408
409	}
410
411
412
413
414
415
416
417
418
419	void mexFunction( int nlhs, mxArray plhs[], int nrhs, const mxArray prhs[])
420	{ int m,n;
421	double *r_in;
422	double le, A, B;
423	int max_order, num_int_steps;
424	double pr1, pr2, pt1, pt2, *ka;
425	mxArray *tmpmxptr;
426	double E0; /* Design energy [eV] */
427
428
429	if(nrhs)
430	{
431	/* ALLOCATE memory for the output array of the same size as the input */
432	m = mxGetM(prhs[1]);
433	n = mxGetN(prhs[1]);
434	if(m!=6)
435	mexErrMsgTxt("Second argument must be a 6 x N matrix");
436
437
438
439	tmpmxptr =mxGetField(prhs[0],0,"PolynomA");
440	if(tmpmxptr)
441	A = mxGetPr(tmpmxptr);
442	else
443	mexErrMsgTxt("Required field 'PolynomA' was not found in the element data structure");
444
445	tmpmxptr =mxGetField(prhs[0],0,"PolynomB");
446	if(tmpmxptr)
447	B = mxGetPr(tmpmxptr);
448	else
449	mexErrMsgTxt("Required field 'PolynomB' was not found in the element data structure");
450
451	tmpmxptr = mxGetField(prhs[0],0,"MaxOrder");
452	if(tmpmxptr)
453	max_order = (int)mxGetScalar(tmpmxptr);
454	else
455	mexErrMsgTxt("Required field 'MaxOrder' was not found in the element data structure");
456
457	tmpmxptr = mxGetField(prhs[0],0,"NumIntSteps");
458	if(tmpmxptr)
459	num_int_steps = (int)mxGetScalar(tmpmxptr);
460	else
461	mexErrMsgTxt("Required field 'NumIntSteps' was not found in the element data structure");
462
463	tmpmxptr = mxGetField(prhs[0],0,"Length");
464	if(tmpmxptr)
465	le = mxGetScalar(tmpmxptr);
466	else
467	mexErrMsgTxt("Required field 'Length' was not found in the element data structure");
468
469	/* Kicks from multipole elements can be specified as angles. This handles the
470	case where corrector coils are used in sextupoles and used for orbit
471	correction. */
472	tmpmxptr = mxGetField(prhs[0],0,"KickAngle");
473	if(tmpmxptr)
474	ka = mxGetPr(tmpmxptr);
475	else
476	ka = NULL;
477
478	tmpmxptr = mxGetField(prhs[0],0,"Energy");
479	if(tmpmxptr)
480	E0 = mxGetScalar(tmpmxptr);
481	else
482	mexErrMsgTxt("Required field 'Energy' was not found in the element data structure");
483
484	tmpmxptr=mxGetField(prhs[0],0,"R1");
485	if(tmpmxptr)
486	pr1 = mxGetPr(tmpmxptr);
487	else
488	pr1 = NULL;
489
490
491	tmpmxptr=mxGetField(prhs[0],0,"R2");
492	if(tmpmxptr)
493	pr2 = mxGetPr(tmpmxptr);
494	else
495	pr2 = NULL;
496
497	tmpmxptr=mxGetField(prhs[0],0,"T1");
498	if(tmpmxptr)
499	pt1 = mxGetPr(tmpmxptr);
500	else
501	pt1 = NULL;
502
503	tmpmxptr=mxGetField(prhs[0],0,"T2");
504	if(tmpmxptr)
505	pt2 = mxGetPr(tmpmxptr);
506	else
507	pt2 = NULL;
508
509
510	if(ka!=NULL)
511	{
512	/* Positive angle must correspond to -ve B field since +ve B field corresponds
513	to a bend with the same curvature as the bend magnets.
514	*/
515	B[0] -= sin(ka[0])/le;
516	A[0] += sin(ka[1])/le;
517	}
518
519	plhs[0] = mxDuplicateArray(prhs[1]);
520	r_in = mxGetPr(plhs[0]);
521
522	StrMPoleSymplectic4RadPass(r_in, le, A, B, max_order, num_int_steps,
523	pt1, pt2, pr1, pr2, E0, n);
524
525	if(ka!=NULL)
526	{
527	B[0] -= sin(ka[0])/le;
528	A[0] -= sin(ka[1])/le;
529	}
530
531
532	}
533	else
534	{ /* return list of required fields */
535	plhs[0] = mxCreateCellMatrix(6,1);
536	mxSetCell(plhs[0],0,mxCreateString("Length"));
537	mxSetCell(plhs[0],1,mxCreateString("PolynomA"));
538	mxSetCell(plhs[0],2,mxCreateString("PolynomB"));
539	mxSetCell(plhs[0],3,mxCreateString("MaxOrder"));
540	mxSetCell(plhs[0],4,mxCreateString("NumIntSteps"));
541	mxSetCell(plhs[0],5,mxCreateString("Energy"));
542
543	if(nlhs>1) /* Required and optional fields */
544	{ plhs[1] = mxCreateCellMatrix(4,1);
545	mxSetCell(plhs[1],0,mxCreateString("T1"));
546	mxSetCell(plhs[1],1,mxCreateString("T2"));
547	mxSetCell(plhs[1],2,mxCreateString("R1"));
548	mxSetCell(plhs[1],3,mxCreateString("R2"));
549	mxSetCell(plhs[1],4,mxCreateString("KickAngle"));
550	}
551	}
552
553
554
555	}
556
557
558

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: MML/trunk/at/simulator/element/user/StrCorrMPoleSymplectic4RadPass.c @ 4

Download in other formats: