Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

BndMPoleSymplectic4RadPass.c @ 8

Last change on this file since 8 was 4, checked in by zhangj, 10 years ago
Initial import--MML version from SOLEIL@2013
File size: 17.9 KB

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

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

Context Navigation

source: MML/trunk/at/simulator/element/BndMPoleSymplectic4RadPass.c @ 8

Download in other formats: