1 | |
---|
2 | #include <nag.h> |
---|
3 | #include <stdio.h> |
---|
4 | #include <nag_stdlib.h> |
---|
5 | #include <math.h> |
---|
6 | #include <nagd01.h> |
---|
7 | #include <nagx01.h> |
---|
8 | #include <stdlib.h> |
---|
9 | #include <string.h> |
---|
10 | #include <nags.h> |
---|
11 | #include <fcntl.h> |
---|
12 | #define NDIM 2 |
---|
13 | |
---|
14 | static double QU1( Integer ndim, double var[]); |
---|
15 | static double QU2( Integer ndim, double var[]); |
---|
16 | static double QU3( Integer ndim, double var[]); |
---|
17 | static double QU4( Integer ndim, double var[]); |
---|
18 | static double GAU(double x); |
---|
19 | |
---|
20 | static double FUN(double x); |
---|
21 | static double FUN_X(double x); |
---|
22 | static double FUN_Y(double y); |
---|
23 | static double LOR1(double t); |
---|
24 | static double LOR2(double t); |
---|
25 | /*static double bessel(double y);*/ |
---|
26 | |
---|
27 | const double pi=3.1415926; |
---|
28 | const double qe=1.602177e-19*3e9; |
---|
29 | double alpha, k, kx, ky, lam_e, om, xv, x0, theta, psi, beta, gam, dum, D; |
---|
30 | double a, b, epsl, sigma_x,sigma_y,sigma_t; |
---|
31 | |
---|
32 | typedef struct {double r, i;} complex; |
---|
33 | typedef struct {complex x, y, z;} vector; |
---|
34 | |
---|
35 | vector Vcreate(double xr, double xi, double yr, double yi, double zr, double zi); |
---|
36 | vector Vadd(vector v1, vector v2); |
---|
37 | vector Vmultiply(vector v, complex c); |
---|
38 | vector Vconjugate(vector v); |
---|
39 | complex Vdot(vector v1, vector v2); |
---|
40 | |
---|
41 | complex Ccreate(double r, double i); |
---|
42 | complex Cadd(complex a, complex b); |
---|
43 | complex Cmultiply(complex a, complex b); |
---|
44 | complex Cdivide(complex a, complex b); |
---|
45 | complex Cconjugate(complex a); |
---|
46 | vector G1,G2, G, nv, nvp, nvv; |
---|
47 | complex G1x, G1y, G1z, G2x, G2y, G2z, Gx, Gy, Gz; |
---|
48 | complex a1, a2, a1t, a2t, f5c; |
---|
49 | |
---|
50 | int main() |
---|
51 | { |
---|
52 | char outfile_name[64]; |
---|
53 | FILE *out_file; |
---|
54 | char line[200], ans,ans1, ans2; |
---|
55 | const double pi=3.1415926; |
---|
56 | const double qe=1.602177e-19*3e9; |
---|
57 | extern double alpha, k, kx, ky, xv, x0, theta, psi, beta, gam, dum, D; |
---|
58 | extern double a,b,sigma_x,sigma_y, sigma_t, lam_e, om; |
---|
59 | int i, ii, ic, j, jj, l, m, n, liw, flag=0; |
---|
60 | int n0, int_thet_1, int_thet_2, jstart, ind, indend; |
---|
61 | Integer ndim=NDIM, minpts, maxpts; |
---|
62 | /* grating parameters below */ |
---|
63 | int grat, Nt; |
---|
64 | double lgr, lgr_e, frac, alpha1, alpha2, hi, el, el1, el2, xx, wid; |
---|
65 | /* beam parameters below */ |
---|
66 | int pulse_shp, ngauss, philmt; |
---|
67 | double y0,y0_av,charlen,nemit, emit, x_rad, n_e, n_e_t, i_av_p ; |
---|
68 | double pulse_len, len50, epsl, sigma_tdum, i_p, x, apb; |
---|
69 | double sigma_t1, sigma_t2, sigma_t3, sigma_t4, sigma_t5, sigma_teff, ssigma_x; |
---|
70 | double mu1, mu2, mu3, mu4, mu5, g1, g2, g3, g4, g5, tot_int; |
---|
71 | double fwhm_x, fwhm_y, fwhm, fwhm_t1, fwhm_t2, fwhm_t3, fwhm_t4, fwhm_t5; |
---|
72 | /* optical system parameters below */ |
---|
73 | double r, r2, domeg_m, theta_u, theta_d, d_theta, d_phi, thet_1, thet_2; |
---|
74 | double phi, phi0, phistep, lam ; |
---|
75 | double aven, aven_p1, aven_p2, aven_single, aven_single_p1, aven_single_p2; |
---|
76 | /* integration parameters below */ |
---|
77 | double ll[2], ul[2]; /* limits for the y,z integral*/ |
---|
78 | double lol, upl, eps, acc, res1, res2, res3, res4 ; |
---|
79 | double epsabs, epsrel, abserr; |
---|
80 | double fmax, tmax, t, f, gtop, gbot; |
---|
81 | double u1, u2, u3, u4, u5; |
---|
82 | double ff, f1, f2, sf3, f_y_c, f_y_s, f_z_c, f_z_s, z, u, s_coh_long; |
---|
83 | double f5, f5_p1, f5_p2; |
---|
84 | double s_coh_xy, s_coh_xy_p1, s_coh_xy_p2, incans, cohans, err; |
---|
85 | double f4[179], en_inc[179]; |
---|
86 | double x0_bin[50], inc_bin[50], coh_bin[50], inc_bin_p1[50], coh_bin_p1[50], inc_bin_p2[50], coh_bin_p2[50]; |
---|
87 | /* general variables */ |
---|
88 | double temp1, temp2, temp3, temp4, ratio, sum1, sum2, sum3, sum4, dumsum, percent; |
---|
89 | double sum_single, sum_single_p1, sum_single_p2 ; |
---|
90 | double sum_th, sum_th_p1, sum_th_p2, sum_th_single, sum_th_single_p1, sum_th_single_p2; |
---|
91 | double sum_tr, sum_tr_p1, sum_tr_p2, sum_tr_single, sum_tr_single_p1, sum_tr_single_p2; |
---|
92 | double sum, sum_p1, sum_p2; |
---|
93 | /* output arrays listed below */ |
---|
94 | double thet[179], en_coh[179], f_p[179]; |
---|
95 | double s_coh_z[179][10]; |
---|
96 | double s_inc_x[179][10], s_inc_x_p1[179][10], s_inc_x_p2[179][10]; |
---|
97 | double s_coh_x[179][10], s_coh_x_p1[179][10], s_coh_x_p2[179][10]; |
---|
98 | double s_coh_tr[179][10], s_coh_tr_p1[179][10], s_coh_tr_p2[179][10]; |
---|
99 | double single_elec[179][15], single_elec_p1[179][15], single_elec_p2[179][15]; |
---|
100 | double bin_single[179][15], bin_single_p1[179][15], bin_single_p2[179][15]; |
---|
101 | double f_b[179][10], f_b_p1[179][10], f_b_p2[179][10]; |
---|
102 | double bin[179][15], bin_p1[179][15], bin_p2[179][15]; |
---|
103 | float outbin[171][8], outbin_single[171][8]; |
---|
104 | double dumm[3], hyperg[500]; |
---|
105 | double bin_1[179][2],dum70[179][2]; |
---|
106 | Nag_QuadProgress QP; |
---|
107 | Nag_QuadSubProgress QPSUB1,QPSUB2; |
---|
108 | Nag_TrigTransform WT_FUNC; |
---|
109 | static NagError fail; |
---|
110 | |
---|
111 | complex G1x, G1y, G1z, G2x, G2y, G2z, Gx, Gy, Gz; |
---|
112 | complex a1, a2, a1t, a2t, f5c; |
---|
113 | |
---|
114 | printf("\n"); printf("\n"); |
---|
115 | printf(" This code calculates the radiated energy, BOTH POLARIZATIONS, (order 1) from a grating\n"); |
---|
116 | printf(" !!!OF FINITE WIDTH!!!\n"); |
---|
117 | printf("\n"); printf("\n"); |
---|
118 | printf(" !!!!!!!!!!!!!! NB1. bunch length is now defined by its FWHM !!!!!!!!!!!!! \n"); |
---|
119 | printf(" !!!!! Gaussians are now defined by their FWHM, rather than SIGMA, as in previous versions!!! \n"); |
---|
120 | printf(" NB2. This version gives the total energy, in both polarisation directions \n"); |
---|
121 | printf("\n"); |
---|
122 | printf(" NB3. scale factor f1 DOES NOT include n_e, which has been moved to S_inc & S_coh \n"); |
---|
123 | printf(" NB4. this version provides data for the single electron yields, diff. & total \n"); |
---|
124 | printf(" NB5. This is the preferred version for the calculation of polarisation \n"); |
---|
125 | |
---|
126 | |
---|
127 | printf("\n"); printf("\n"); |
---|
128 | |
---|
129 | m=1; /*default order 1*/ |
---|
130 | flag=0; |
---|
131 | if (!flag) printf("Default order is 1. Do you want to change it? yes(y) or no(n) \n"); |
---|
132 | line[0]= 'x'; |
---|
133 | do { |
---|
134 | if (line[0] != 'x') printf(" "); |
---|
135 | if (fgets(line, 200, stdin) == NULL) { printf("\nErr: fgets.\n"); exit(1); } |
---|
136 | } |
---|
137 | while (line[0] != 'y' && line[0] != 'n'); |
---|
138 | ans = line[0]; |
---|
139 | if (ans=='y') |
---|
140 | {printf("Enter order required; default is order 1 \n"); |
---|
141 | scanf("%d",&m); |
---|
142 | if (m<=0) m=1; } |
---|
143 | printf(" \n"); |
---|
144 | |
---|
145 | lgr=40; /* grating length*/ |
---|
146 | printf(" Enter grating period in mm\n"); |
---|
147 | scanf("%lf",&el); |
---|
148 | Nt= (int)lgr/el; /*Nt is the number of teeth for this grating*/ |
---|
149 | printf(" A standard grating of 40mm length would have %5d", Nt); printf(" teeth\n"); |
---|
150 | printf(" Enter grating width in mm; if <1, it will be set to the default value of 20mm \n"); |
---|
151 | scanf("%lf",&wid); |
---|
152 | if (wid<=1) wid=20.0; |
---|
153 | |
---|
154 | grat=2; |
---|
155 | /*printf(" Now select type of grating \n"); |
---|
156 | printf(" Type:\n"); |
---|
157 | printf(" -3 for strip grating\n"); |
---|
158 | printf(" -1 for reverse blaze echelle grating \n"); |
---|
159 | printf(" 0 for idealized strip grating \n"); |
---|
160 | printf(" 1 for normal blaze echelle grating \n"); |
---|
161 | printf(" 2 for sawtooth grating \n"); |
---|
162 | printf(" Default is the sawtooth grating \n"); |
---|
163 | scanf("%d",&grat);*/ |
---|
164 | printf("\n"); |
---|
165 | if (grat==-3) |
---|
166 | {printf(" This is a strip grating \n"); |
---|
167 | alpha1=0; |
---|
168 | printf(" Enter height of vertical facet, in mm \n"); |
---|
169 | scanf("%lf",&hi); |
---|
170 | alpha2=0; |
---|
171 | printf(" Enter the fraction of the period taken by the strip \n"); |
---|
172 | scanf("%lf",&frac); |
---|
173 | if (frac<0.0 || frac>1.0) frac=0.5; |
---|
174 | el1=el*frac; |
---|
175 | el2=el-el1;} |
---|
176 | else if (grat==-1) |
---|
177 | {printf("This is a reverse blaze echelle grating \n"); |
---|
178 | frac=0; |
---|
179 | el1=0; |
---|
180 | alpha1=pi/2; |
---|
181 | el2=el; |
---|
182 | printf(" Enter second blaze angle of grating in degrees; MUST BE greater than 90 \n"); |
---|
183 | scanf("%lf",&alpha2); |
---|
184 | alpha2=alpha2*pi/180; |
---|
185 | hi=fabs(el2*tan(alpha2));} |
---|
186 | else if (grat==0) |
---|
187 | {printf(" This is an idealized strip grating \n"); |
---|
188 | alpha1=0; |
---|
189 | printf(" Enter the fraction of the period taken by the strip \n"); |
---|
190 | scanf("%lf",&frac); |
---|
191 | if (frac<0.0 || frac>1.0) frac=0.5; |
---|
192 | el1=el*frac; |
---|
193 | alpha2=0;} |
---|
194 | else if (grat==1) |
---|
195 | {printf("This is an echelle-type grating \n"); |
---|
196 | frac=1.0; |
---|
197 | el1=el; |
---|
198 | el2=0; |
---|
199 | alpha2= pi/2; |
---|
200 | printf("Enter first blaze angle of grating, in degrees \n"); |
---|
201 | scanf("%lf",&alpha1); |
---|
202 | alpha1=alpha1*pi/180; |
---|
203 | hi=fabs(el*tan(alpha1));} |
---|
204 | else |
---|
205 | { |
---|
206 | if (el==1.5) alpha1=30; |
---|
207 | else if (el==1.0) alpha1=35; |
---|
208 | else if(el==0.5) alpha1=40; |
---|
209 | else {printf("Enter first blaze angle of grating in degrees \n"); |
---|
210 | scanf("%lf",&alpha1);} |
---|
211 | alpha1=alpha1*pi/180; |
---|
212 | frac=pow(cos(alpha1),2); |
---|
213 | printf("Fraction of period taken by first blaze is= "); printf("%-4.3f", frac, "\n"); |
---|
214 | el1=el*frac; |
---|
215 | el2=el-el1; |
---|
216 | hi=fabs(el1*tan(alpha1)); /*height of grating indentation*/ |
---|
217 | alpha2= -(0.5*pi- alpha1); /* definition of alpha2 changed, 10/1/10 */ |
---|
218 | /*alpha2=0.5*pi+ alpha1;*/} |
---|
219 | printf("\n"); |
---|
220 | |
---|
221 | r=230; /* distance of Winston cone entrance from beam axis, in mm */ |
---|
222 | r2=21; /* Winston cone entrance diameter, in mm */ |
---|
223 | /* Max. entry angle for Winston cone is 6.3 deg.*/ |
---|
224 | domeg_m=pi*(r2/2)*(r2/2)/pow(r,2); /*max. solid angle */ |
---|
225 | printf("\n"); |
---|
226 | printf(" With these dimensions ,we have: \n"); |
---|
227 | printf("Max. solid angle at 90 is= "); printf("%-7.4f",domeg_m); printf(" sr \n"); |
---|
228 | |
---|
229 | /*beam parameters */ |
---|
230 | |
---|
231 | printf(" Enter value of gamma; if <1 it will be set equal to 55773 \n"); /*default is SLAC */ |
---|
232 | scanf("%lf",&gam); |
---|
233 | if (gam<=1) gam=55773; |
---|
234 | beta=sqrt(pow(gam,2)-1)/gam; |
---|
235 | printf("Type normalized emittance of the beam in mm.mrad \n"); |
---|
236 | scanf("%le",&nemit); |
---|
237 | emit=1e-3*nemit/(beta*gam); /*converted to mm.rad*/ |
---|
238 | esc20:printf("Enter position (in mm) of beam centroid above grating\n"); |
---|
239 | scanf("%lf",&x0); |
---|
240 | printf("\n"); |
---|
241 | |
---|
242 | esc24:printf("Enter FWHM_X(in mm) of Gaussian beam at the centre of grating \n"); |
---|
243 | scanf("%le",&fwhm_x); |
---|
244 | sigma_x= fwhm_x/(2*sqrt(2*log(2))); |
---|
245 | printf(" For this fwhm: \n"); |
---|
246 | printf(" 0.50 of the beam is within a radius of "); printf("%-4.2e",0.6749*sigma_x); printf(" mm \n"); |
---|
247 | printf(" 0.90 of the beam is within a radius of "); printf("%-4.2e",1.6450*sigma_x); printf(" mm \n"); |
---|
248 | printf(" 0.99 of the beam is within a radius of "); printf("%-4.2e",2.5733*sigma_x); printf(" mm \n"); |
---|
249 | y0=1.6450*sigma_x; /*average value of radius over 1/2 grating*/ |
---|
250 | charlen=pow(y0,2)/emit; |
---|
251 | sf3=0; |
---|
252 | j=1; |
---|
253 | for (j=1; j<7; j++) /*first calculate beam profile over grating*/ |
---|
254 | {z=(j-1)*lgr/10; |
---|
255 | x_rad=y0*sqrt(1+pow(z/charlen,2)); /*varying beam half-width*/ |
---|
256 | sf3=sf3+x_rad;} |
---|
257 | y0_av=sf3/6; /*average value of radius over 1/2 grating*/ |
---|
258 | sigma_x= y0_av/1.645; /*average value of sigma_x over 1/2 grating*/ |
---|
259 | printf("Average value of FWHM_X from centre to edge of grating is "); printf("%-4.2e",(2*sigma_x*sqrt(2*log(2)))); printf(" mm\n"); |
---|
260 | printf("Enter FWHM_Y(in mm) of Gaussian beam at waist, at the centre of grating \n"); |
---|
261 | scanf("%le",&fwhm_y); |
---|
262 | sigma_y= fwhm_y/(2*sqrt(2*log(2))); |
---|
263 | printf(" For this fwhm: \n"); |
---|
264 | printf(" 0.50 of the beam is within a radius of "); printf("%-4.2e",0.6749*sigma_y); printf(" mm \n"); |
---|
265 | printf(" 0.90 of the beam is within a radius of "); printf("%-4.2e",1.6450*sigma_y); printf(" mm \n"); |
---|
266 | printf(" 0.99 of the beam is within a radius of "); printf("%-4.2e",2.5733*sigma_y); printf(" mm \n"); |
---|
267 | |
---|
268 | lol=x0-5*sigma_x; /* GD modified 5/6/06*/ |
---|
269 | if (lol<hi) lol=hi; |
---|
270 | upl=x0+5*sigma_x; |
---|
271 | epsabs=1.0e-12; |
---|
272 | epsrel=1.0e-6; |
---|
273 | liw=500; |
---|
274 | d01ajc(GAU,lol,upl,epsabs,epsrel,liw,&res1,&abserr,&QP,&fail); |
---|
275 | res1=res1/(sigma_x*sqrt(2*pi)); |
---|
276 | printf("For a fwhm of "); printf("%-4.2e",fwhm_x); |
---|
277 | printf("mm and x0="); printf("%-4.3f",x0); printf("mm, "); printf("%-4.2f",100*res1); |
---|
278 | printf(" percent of the beam is above the grating \n"); |
---|
279 | printf("\n"); |
---|
280 | NAG_FREE(QP.sub_int_beg_pts); |
---|
281 | NAG_FREE(QP.sub_int_end_pts); |
---|
282 | NAG_FREE(QP.sub_int_result); |
---|
283 | NAG_FREE(QP.sub_int_error); |
---|
284 | |
---|
285 | /* bunch parameters */ |
---|
286 | esc1: printf( "\n"); |
---|
287 | printf("Enter total number of electrons in the bunch \n"); |
---|
288 | scanf("%le",&n_e_t); |
---|
289 | printf("since only "); printf("%-4.1f",100*res1); printf(" percent is above the grating, I am setting N_e="); |
---|
290 | printf("%le",res1*n_e_t); printf("\n"); |
---|
291 | n_e=n_e_t*res1; |
---|
292 | esc25: printf("Enter nominal bunch duration (FWHM), in ps \n"); |
---|
293 | scanf("%lf",&pulse_len); |
---|
294 | i_av_p=n_e_t*1.602177*1e-7/pulse_len; |
---|
295 | printf("The average current in the pulse is= "); printf("%-4.3f",i_av_p); printf(" A \n"); |
---|
296 | |
---|
297 | printf("You now have the following choice for profiles in z-direction \n"); |
---|
298 | printf("Type one of the following: \n"); |
---|
299 | printf(" 0 for a DC beam \n"); |
---|
300 | printf(" 1 for a single Gaussian \n"); printf(" 2 for triangular \n"); printf(" 3 for double-exp. \n"); |
---|
301 | printf(" 4 for parabolic \n"); printf(" 5 for cosine \n"); printf(" 6 NOT IN USE \n"); |
---|
302 | printf(" 7 for Lorentz shape \n"); printf(" 8 for a multi-Gauss shape \n"); |
---|
303 | printf(" which one do you want?? \n"); |
---|
304 | scanf("%d",&pulse_shp); |
---|
305 | if (pulse_shp<0 || pulse_shp>8) |
---|
306 | { printf(" I do not recognize this \n"); |
---|
307 | goto esc1;} |
---|
308 | if (pulse_shp==1) /*Gaussian*/ |
---|
309 | { |
---|
310 | printf("Type asymmetry factor epsl \n"); |
---|
311 | scanf("%lf",&epsl); |
---|
312 | sigma_t= pulse_len/(sqrt(2*log(2))*(1+epsl)); /*length defined by FWHM 15/5/09*/ |
---|
313 | printf(" the sigma of the leading edge is "); printf("%-4.3f",sigma_t); printf(" ps \n"); |
---|
314 | printf(" and sigma of the trailing edge is "); |
---|
315 | printf("%-4.3f",epsl*sigma_t); printf(" ps \n"); |
---|
316 | } |
---|
317 | else if (pulse_shp==2) /*triangular, modified 11/05/09*/ |
---|
318 | { |
---|
319 | printf("Type asymmetry factor epsl \n"); |
---|
320 | scanf("%lf",&epsl); |
---|
321 | sigma_t= pulse_len/((1+epsl)*(1-0.5)); /* now in terms of FWHM, 15/05/09 */ |
---|
322 | printf(" For this bunch shape, the sigma_t is "); printf("%-4.3f",sigma_t); printf(" ps \n"); |
---|
323 | } |
---|
324 | else if (pulse_shp==3) /*double exponential*/ |
---|
325 | { |
---|
326 | printf("Type asymmetry factor epsl \n"); |
---|
327 | scanf("%lf",&epsl); |
---|
328 | sigma_t=pulse_len/((1+epsl)*log(1/0.5)); /*changed 15/05/09, now using FWHM */ |
---|
329 | printf(" the sigma_t is "); printf("%-4.3f",sigma_t); printf(" ps \n"); |
---|
330 | } |
---|
331 | else if (pulse_shp==4) /*parabolic*/ |
---|
332 | { |
---|
333 | printf("Type asymmetry factor epsl \n"); |
---|
334 | scanf("%lf",&epsl); |
---|
335 | sigma_t=pulse_len/((1+epsl)*sqrt(1-0.5)); /*changed 15/05/09, using FWHM*/ |
---|
336 | printf(" For this bunch shape, sigma_t is "); printf("%-4.3f",sigma_t); printf(" ps \n"); |
---|
337 | } |
---|
338 | else if (pulse_shp==5) /*cosine*/ |
---|
339 | { |
---|
340 | epsl=1; |
---|
341 | sigma_t=pulse_len*pi/(4*acos(0.5)); /* changed 15/05/09, using FWHM*/ |
---|
342 | printf(" For this bunch shape, sigma_t is "); printf("%-4.3f",sigma_t); printf(" ps \n"); |
---|
343 | } |
---|
344 | else if (pulse_shp==6) |
---|
345 | { |
---|
346 | goto esc1; |
---|
347 | } |
---|
348 | else if (pulse_shp==7) /*Lorentz*/ |
---|
349 | { |
---|
350 | printf("Enter the asymmetry factor epsl; symmetric shape means eps=1 \n"); |
---|
351 | scanf("%lf",&epsl); |
---|
352 | apb=pulse_len; /* changed 15/5/09*/ |
---|
353 | printf("a+b is = "); printf("%-4.3f",apb); printf(" ps \n"); |
---|
354 | a=apb/(1+epsl); /* no need for sigma_t; shape is defined by a & b */ |
---|
355 | b=a*epsl; |
---|
356 | printf(" a(leading edge)= "); printf("%-4.3f",a); printf("ps \n"); |
---|
357 | printf(" b(trailing edge)= "); printf("%-4.3f",b); printf("ps \n"); |
---|
358 | } |
---|
359 | else if (pulse_shp==8) /*multiple Gaussian, June 08, modified for FWHM 11/05/09*/ |
---|
360 | { |
---|
361 | esc81:fwhm_t1= fwhm_t2= fwhm_t3= fwhm_t4= fwhm_t5=0.0; |
---|
362 | sigma_t1= sigma_t2= sigma_t3= sigma_t4= sigma_t5=0.0; |
---|
363 | g1=g2=g3=g4=g5=0.0; |
---|
364 | mu1=mu2=mu3=mu4=mu5=0.0; |
---|
365 | gtop=gbot=len50=0.0; |
---|
366 | printf(" how many gaussians are you fitting? Min. is 2 and max. is 5 \n"); |
---|
367 | scanf("%1d", &ngauss); |
---|
368 | res1=0; |
---|
369 | printf(" for the 1st peak, enter fwhm(ps), amplitude & position(ps) \n"); |
---|
370 | scanf("%lf %lf %lf", &fwhm_t1, &g1, &mu1); |
---|
371 | sigma_t1= fwhm_t1/(2*sqrt(2*log(2))); |
---|
372 | if (ngauss==2) |
---|
373 | {printf(" for the 2nd peak, enter fwhm(ps), amplitude & position(ps) \n"); |
---|
374 | scanf("%lf %lf %lf", &fwhm_t2, &g2, &mu2); |
---|
375 | sigma_t2= fwhm_t2/(2*sqrt(2*log(2))); |
---|
376 | sigma_t3=sigma_t4=sigma_t5=0; |
---|
377 | g3=g4=g5=0; |
---|
378 | mu3=mu4=mu5=0;} |
---|
379 | else if (ngauss==3) |
---|
380 | {printf(" for the 2nd peak, enter fwhm(ps), amplitude & position(ps) \n"); |
---|
381 | scanf("%lf %lf %lf", &fwhm_t2, &g2, &mu2); |
---|
382 | sigma_t2= fwhm_t2/(2*sqrt(2*log(2))); |
---|
383 | printf(" for the 3rd peak, enter fwhm(ps), amplitude & position(ps) \n"); |
---|
384 | scanf("%lf %lf %lf", &fwhm_t3, &g3, &mu3); |
---|
385 | sigma_t3= fwhm_t3/(2*sqrt(2*log(2))); |
---|
386 | sigma_t4=sigma_t5=0; |
---|
387 | g4=g5=0; |
---|
388 | mu4=mu5=0;} |
---|
389 | else if (ngauss==4) |
---|
390 | {printf(" for the 2nd peak, enter fwhm(ps), amplitude & position(ps) \n"); |
---|
391 | scanf("%lf %lf %lf", &fwhm_t2, &g2, &mu2); |
---|
392 | sigma_t2= fwhm_t2/(2*sqrt(2*log(2))); |
---|
393 | printf(" for the 3rd peak, enter fwhm(ps), amplitude & position(ps) \n"); |
---|
394 | scanf("%lf %lf %lf", &fwhm_t3, &g3, &mu3); |
---|
395 | sigma_t3= fwhm_t3/(2*sqrt(2*log(2))); |
---|
396 | printf(" for the 4th peak, enter fwhm(ps), amplitude & position(ps) \n"); |
---|
397 | scanf("%lf %lf %lf", &fwhm_t4, &g4, &mu4); |
---|
398 | sigma_t4= fwhm_t4/(2*sqrt(2*log(2))); |
---|
399 | sigma_t5=0; |
---|
400 | g5=0; |
---|
401 | mu5=0;} |
---|
402 | else |
---|
403 | {printf(" for the 2nd peak, enter fwhm(ps), amplitude & position(ps) \n"); |
---|
404 | scanf("%lf %lf %lf", &fwhm_t2, &g2, &mu2); |
---|
405 | sigma_t2= fwhm_t2/(2*sqrt(2*log(2))); |
---|
406 | printf(" for the 3rd peak, enter fwhm(ps), amplitude(ps) & position(ps) \n"); |
---|
407 | scanf("%lf %lf %lf", &fwhm_t3, &g3, &mu3); |
---|
408 | sigma_t3= fwhm_t3/(2*sqrt(2*log(2))); |
---|
409 | printf(" for the 4th peak, enter fwhm(ps), amplitude & position(ps) \n"); |
---|
410 | scanf("%lf %lf %lf", &fwhm_t4, &g4, &mu4); |
---|
411 | sigma_t4= fwhm_t4/(2*sqrt(2*log(2))); |
---|
412 | printf(" for the 5th peak, enter fwhm(ps), amplitude & position(ps) \n"); |
---|
413 | scanf("%lf %lf %lf", &fwhm_t5, &g5, &mu5); |
---|
414 | sigma_t5= fwhm_t5/(2*sqrt(2*log(2)));} |
---|
415 | |
---|
416 | /* check total value of this integral, i.e. normalisation factor */ |
---|
417 | tot_int= sqrt(2*pi)*(g1*sigma_t1+g2*sigma_t2+g3*sigma_t3+g4*sigma_t4+g5*sigma_t5); |
---|
418 | printf(" the normalisation factor is "); |
---|
419 | printf("%-4.3e",tot_int); printf(" units \n"); |
---|
420 | sigma_teff= g1*sigma_t1+ g2*sigma_t2+g3*sigma_t3+ g4*sigma_t4+g5*sigma_t5; |
---|
421 | /* find position & value of maximum of normalised distribution; total scan width=15ps */ |
---|
422 | fmax=0; |
---|
423 | tmax=0; |
---|
424 | i=0; |
---|
425 | for (i==0; i<1500; i++) |
---|
426 | { |
---|
427 | t= 1*sigma_teff- i*0.01; /* start well into the positive region*/ |
---|
428 | f=(g1*exp(-(t-mu1)*(t-mu1)/(2*sigma_t1*sigma_t1))+g2*exp(-(t-mu2)*(t-mu2)/(2*sigma_t2*sigma_t2))+ |
---|
429 | g3*exp(-(t-mu3)*(t-mu3)/(2*sigma_t3*sigma_t3))+ |
---|
430 | g4*exp(-(t-mu4)*(t-mu4)/(2*sigma_t4*sigma_t4))+ |
---|
431 | g5*exp(-(t-mu5)*(t-mu5)/(2*sigma_t5*sigma_t5)))/tot_int; |
---|
432 | if (f>= fmax) {fmax=f; tmax=t;} |
---|
433 | } |
---|
434 | printf(" maximum for the (normalised) shape is at t="); |
---|
435 | printf("%-4.3f", tmax); printf(" \n"); |
---|
436 | printf(" and is equal to "); |
---|
437 | printf("%-4.3f", fmax); printf(" \n"); |
---|
438 | /* finished with this; now find points where value of function is about 0.5 */ |
---|
439 | /* and use these as the limits for the evaluation of the normalised integral */ |
---|
440 | /* start on positive side...*/ |
---|
441 | i=0; |
---|
442 | for (i==0; i<1500; i++) |
---|
443 | { |
---|
444 | t= tmax + i*0.01; |
---|
445 | f=(g1*exp(-(t-mu1)*(t-mu1)/(2*sigma_t1*sigma_t1))+g2*exp(-(t-mu2)*(t-mu2)/(2*sigma_t2*sigma_t2))+ |
---|
446 | g3*exp(-(t-mu3)*(t-mu3)/(2*sigma_t3*sigma_t3))+ |
---|
447 | g4*exp(-(t-mu4)*(t-mu4)/(2*sigma_t4*sigma_t4))+ |
---|
448 | g5*exp(-(t-mu5)*(t-mu5)/(2*sigma_t5*sigma_t5)))/tot_int; |
---|
449 | if (f<= 0.5*fmax) {gtop=t; goto esc82;} |
---|
450 | } |
---|
451 | esc82: printf("the upper 50 percent point of the distribution is "); printf("%-4.3f",gtop); |
---|
452 | /* continue on negative...*/ |
---|
453 | i=0; |
---|
454 | for (i==0; i<1500; i++) |
---|
455 | { |
---|
456 | t= tmax - i*0.01; |
---|
457 | f=(g1*exp(-(t-mu1)*(t-mu1)/(2*sigma_t1*sigma_t1))+g2*exp(-(t-mu2)*(t-mu2)/(2*sigma_t2*sigma_t2))+ |
---|
458 | g3*exp(-(t-mu3)*(t-mu3)/(2*sigma_t3*sigma_t3))+ |
---|
459 | g4*exp(-(t-mu4)*(t-mu4)/(2*sigma_t4*sigma_t4))+ |
---|
460 | g5*exp(-(t-mu5)*(t-mu5)/(2*sigma_t5*sigma_t5)))/tot_int; |
---|
461 | if (f<= 0.5*fmax) {gbot=t; goto esc83;} |
---|
462 | } |
---|
463 | esc83: printf(" and the lower one "); printf("%-4.3f",gbot); printf("\n"); |
---|
464 | /* done*/ |
---|
465 | esc80: len50= gtop- gbot; /* this is the assumed bunch length */ |
---|
466 | printf(" Hence, the assumed bunch length (FWHM) is "); printf("%-4.3f",len50); printf(" ps \n"); |
---|
467 | flag=0; |
---|
468 | if (!flag) printf(" Do you want to change anything? yes(y) or no(n) \n"); |
---|
469 | line[0]= 'x'; |
---|
470 | do { |
---|
471 | if (line[0] != 'x') printf(" "); |
---|
472 | if (fgets(line, 200, stdin) == NULL) { printf("\nErr: fgets.\n"); exit(1); } |
---|
473 | } |
---|
474 | while (line[0] != 'y' && line[0] != 'n'); |
---|
475 | ans = line[0]; |
---|
476 | if (ans=='y') |
---|
477 | {printf("go back to input parameters \n"); |
---|
478 | goto esc81; |
---|
479 | } |
---|
480 | printf(" \n"); |
---|
481 | } |
---|
482 | |
---|
483 | /*estimate of grating factor R^2 and of the coherent & incoherent integrals*/ |
---|
484 | |
---|
485 | printf("This code can consider up to 7 steps in the azimuthal direction phi \n"); |
---|
486 | printf(" Define step size: min.=0.1deg, max.=10deg, default=1.0deg. \n"); |
---|
487 | scanf("%lf", &phistep); |
---|
488 | if (phistep<0.1 || phistep>10) phistep=1.0; |
---|
489 | printf("OK, phi step is equal to= "); printf("%-3.1f", phistep); printf(" deg.\n"); |
---|
490 | |
---|
491 | /* set the y-limits for the y-z integrations*/ |
---|
492 | ll[0]= 0.0; ul[0]= wid/2; /*y-integral; wid is the grating width */ |
---|
493 | maxpts= 5000; |
---|
494 | eps=1e-4; /* acceptable error*/ |
---|
495 | |
---|
496 | i=(int)(alpha1*180/pi); |
---|
497 | for ( i=(int)(alpha1*180/pi); i<171; i++) /* angles lower than the blaze angle do not matter! */ |
---|
498 | { thet[i]=i; /*scan angles from 10-170 deg. */ |
---|
499 | theta=thet[i]*pi/180; |
---|
500 | lam=el*(1/beta-cos(theta))/m; /*lam is the wavelength of the emitted radiation*/ |
---|
501 | k= 2*pi/lam; |
---|
502 | if (i==10||i==20||i==30|| i==40||i==50||i==60 || i==70 || i==80||i==90 || i==100 || i==110 || i==120 |
---|
503 | ||i==130 || i==140 || i==150 || i==160 || i==170) |
---|
504 | {printf("%-.1f ", (theta*180/pi)); printf("\n");} |
---|
505 | ii=0; |
---|
506 | for (ii=0; ii<8; ii++) |
---|
507 | {phi0=phistep*ii*pi/180; /* azim. angle phi scanned up to 7xphistep deg.*/ |
---|
508 | kx=k*sin(theta)*cos(phi0); |
---|
509 | ky=k*sin(theta)*sin(phi0); |
---|
510 | lam_e=lam*beta*gam/(2*pi*sqrt(1+pow((beta*gam*sin(theta)*sin(phi0)),2))); /*evanescent w/l */ |
---|
511 | ff= 1/(pi*el); |
---|
512 | /* N.B. all distances from the base of the grating teeth. */ |
---|
513 | |
---|
514 | /* calculate for a range of x0 values and store for integration */ |
---|
515 | lol= x0-5*sigma_x; |
---|
516 | if (lol<=hi) lol=1.001*hi; /* electrons must not hit the grating!!*/ |
---|
517 | upl= x0+5*sigma_x; |
---|
518 | indend=20; /* number of steps in x, i.e. in 0.5sigma steps */ |
---|
519 | ind=0; |
---|
520 | for (ind=0; ind<indend+1; ind++) |
---|
521 | {xv= lol+ ind*0.5*sigma_x; /* proceed in steps of 0.5sigma */ |
---|
522 | /* first facet */ |
---|
523 | res1=res2=res3=res4=0.0; |
---|
524 | ll[1]=0.0; ul[1]=el1; /* z limits for 1st facet */ |
---|
525 | alpha=alpha1; |
---|
526 | psi=0; |
---|
527 | dum=0.0; /* dum accounts for diff. expressions for impact param. for the two facets */ |
---|
528 | D= k*(1/beta-cos(theta))-kx*tan(alpha); |
---|
529 | temp1=temp2=temp3=temp4=0.0; |
---|
530 | minpts=1; |
---|
531 | d01fcc(ndim, QU1, ll, ul, &minpts, maxpts, eps, &res1, &acc, &fail); /* real part of the int. for the x&z components */ |
---|
532 | minpts=1; |
---|
533 | d01fcc(ndim, QU2, ll, ul, &minpts, maxpts, eps, &res2, &acc, &fail); /* imag. part */ |
---|
534 | G1x=Ccreate((ff*res1*tan(alpha1)), (ff*res2*tan(alpha1))); /*changed 25/2/10*/ |
---|
535 | G1z=Ccreate((ff*res1), (ff*res2)); /*changed 25/2/10*/ |
---|
536 | |
---|
537 | temp3=temp4=0.0; |
---|
538 | minpts=1; |
---|
539 | d01fcc(ndim, QU3, ll, ul, &minpts, maxpts, eps, &res3, &acc, &fail); /*real part for the integr. of y component */ |
---|
540 | minpts=1; |
---|
541 | d01fcc(ndim, QU4, ll, ul, &minpts, maxpts, eps, &res4, &acc, &fail); /*imag. part of y-integral */ |
---|
542 | G1y=Ccreate((ff*res3*tan(alpha1)), (-ff*res4*tan(alpha1))); /*changed 25/2/10*/ |
---|
543 | |
---|
544 | /* second facet,including phase correction psi */ |
---|
545 | |
---|
546 | res1=res2=res3=res4=0.0; |
---|
547 | ll[1]=el1; ul[1]=el; /* z limits for 2nd facet */ |
---|
548 | alpha=alpha2; |
---|
549 | psi=2*pi*el1/el; |
---|
550 | dum= el*tan(alpha2); /* dum accounts for diff. expressions for impact param. for the two facets */ |
---|
551 | D= k*(1/beta-cos(theta))-kx*tan(alpha); |
---|
552 | temp1=temp2=temp3=temp4=0.0; |
---|
553 | minpts=1; |
---|
554 | d01fcc(ndim, QU1, ll, ul, &minpts, maxpts, eps, &res1, &acc, &fail); |
---|
555 | minpts=1; |
---|
556 | d01fcc(ndim, QU2, ll, ul, &minpts, maxpts, eps, &res2, &acc, &fail); |
---|
557 | G2x=Ccreate((ff*res1*tan(alpha2)), (ff*res2*tan(alpha2))); /*changed 25/2/10*/ |
---|
558 | G2z=Ccreate((ff*res1), (ff*res2)); /*changed 25/2/10*/ |
---|
559 | |
---|
560 | temp3=temp4=0.0; |
---|
561 | minpts=1; |
---|
562 | d01fcc(ndim, QU3, ll, ul, &minpts, maxpts, eps, &res3, &acc, &fail); |
---|
563 | minpts=1; |
---|
564 | d01fcc(ndim, QU4, ll, ul, &minpts, maxpts, eps, &res4, &acc, &fail); |
---|
565 | G2y=Ccreate((ff*res3*tan(alpha2)), (-ff*res4*tan(alpha2))); /*changed 25/2/10*/ |
---|
566 | /* combine the 2 facets... */ |
---|
567 | Gx= Cadd(G1x, G2x); |
---|
568 | Gy= Cadd(G1y, G2y); |
---|
569 | Gz= Cadd(G1z, G2z); |
---|
570 | /* and create the complex vector G */ |
---|
571 | G= Vcreate( Gx.r, Gx.i, Gy.r, Gy.i, Gz.r, Gz.i); |
---|
572 | nv= Vcreate( (sin(theta)*cos(phi0)), 0, (sin(theta)*sin(phi0)), 0, cos(theta), 0); /*dir. unit vector*/ |
---|
573 | nvp= Vcreate( (cos(theta)*cos(phi0)), 0, (cos(theta)*sin(phi0)), 0, -sin(theta), 0); /*1st pol. unit vector*/ |
---|
574 | nvv= Vcreate( -sin(phi0), 0, cos(phi0), 0, 0, 0); /*2nd pol. unit vector*/ |
---|
575 | |
---|
576 | /* first polarization e1, i.e vector nvp, in plane of n,z. Need G.nvp */ |
---|
577 | a1=Vdot(G,nvp); |
---|
578 | a1t= Cmultiply(a1, Cconjugate(a1)); /*square of magn. in first polarisation*/ |
---|
579 | /* second polarization e2, i.e. vector nvv, vertical to n and z */ |
---|
580 | a2=Vdot(G,nvv); |
---|
581 | a2t= Cmultiply(a2, Cconjugate(a2)); /*square of magn. in second polarisation*/ |
---|
582 | f5c=Cadd(a1t,a2t); /*total */ |
---|
583 | f5= f5c.r; /* output is R^2, in both polarisations, at this specific value of x*/ |
---|
584 | f5_p1= a1t.r; /* polar. p1 */ |
---|
585 | f5_p2= a2t.r; /* .. and p2... */ |
---|
586 | if (ind==10) |
---|
587 | {single_elec[i][ii]= f5; |
---|
588 | single_elec_p1[i][ii]= f5_p1; |
---|
589 | single_elec_p2[i][ii]= f5_p2;} /*store the single-electron R^2 */ |
---|
590 | x0_bin[ind]= xv; |
---|
591 | inc_bin[ind]= f5*exp(-(xv-x0)*(xv-x0)/(2*sigma_x*sigma_x)); |
---|
592 | inc_bin_p1[ind]= f5_p1*exp(-(xv-x0)*(xv-x0)/(2*sigma_x*sigma_x)); |
---|
593 | inc_bin_p2[ind]= f5_p2*exp(-(xv-x0)*(xv-x0)/(2*sigma_x*sigma_x)); |
---|
594 | coh_bin[ind]= sqrt(f5)*exp(-(xv-x0)*(xv-x0)/(2*sigma_x*sigma_x)); |
---|
595 | coh_bin_p1[ind]= sqrt(f5_p1)*exp(-(xv-x0)*(xv-x0)/(2*sigma_x*sigma_x)); |
---|
596 | coh_bin_p2[ind]= sqrt(f5_p2)*exp(-(xv-x0)*(xv-x0)/(2*sigma_x*sigma_x)); |
---|
597 | |
---|
598 | } /*end of x0 scan */ |
---|
599 | |
---|
600 | /* numerical integration over x0 follows, to find the coherent and incoherent integrals*/ |
---|
601 | /* the two arrays below contain the incoherent and coherent parts of the x-integrals */ |
---|
602 | /* arranged so that in col 0 data for phi=0, in 1 for phi=1(10), in 7 for phi=7(70)*/ |
---|
603 | d01gac(indend, x0_bin, inc_bin, &incans, &err, &fail); |
---|
604 | s_inc_x[i][ii]= incans/(sqrt(2*pi)*sigma_x); |
---|
605 | d01gac(indend, x0_bin, coh_bin, &cohans, &err, &fail); |
---|
606 | s_coh_x[i][ii]= cohans*cohans/(2*pi*sigma_x*sigma_x); |
---|
607 | |
---|
608 | d01gac(indend, x0_bin, inc_bin_p1, &incans, &err, &fail); |
---|
609 | s_inc_x_p1[i][ii]= incans/(sqrt(2*pi)*sigma_x); |
---|
610 | d01gac(indend, x0_bin, coh_bin_p1, &cohans, &err, &fail); |
---|
611 | s_coh_x_p1[i][ii]= cohans*cohans/(2*pi*sigma_x*sigma_x); |
---|
612 | |
---|
613 | d01gac(indend, x0_bin, inc_bin_p2, &incans, &err, &fail); |
---|
614 | s_inc_x_p2[i][ii]= incans/(sqrt(2*pi)*sigma_x); |
---|
615 | d01gac(indend, x0_bin, coh_bin_p2, &cohans, &err, &fail); |
---|
616 | s_coh_x_p2[i][ii]= cohans*cohans/(2*pi*sigma_x*sigma_x); |
---|
617 | |
---|
618 | } /* end of ii (phi) scan */ |
---|
619 | |
---|
620 | } /* end of i (theta) scan */ |
---|
621 | |
---|
622 | |
---|
623 | printf(" The following is a list of the values of R^2, total, P1, P2 averaged over x-dimension of bunch\n"); |
---|
624 | printf(" NB!! The dimensionless quantity R^2 must be multiplied by the grating period in order to obtain\n"); |
---|
625 | printf(" the true measure of the efficiencies of two different gratings.\n"); |
---|
626 | printf(" See equation 8a of the Theory note for justification\n"); |
---|
627 | printf(" Each column below is for a phi step of "); printf("%-3.1f",phistep); printf(" deg. \n"); |
---|
628 | printf(" first column is the wavelength, in microns \n"); |
---|
629 | printf(" last column is the average over the 3 values of phi on either side of zero, i.e.\n"); |
---|
630 | printf(" between -%-3.1f & +%-3.1f deg. only, in this case \n",phistep*3, phistep*3); |
---|
631 | i=(int)(alpha1*180/pi); |
---|
632 | for ( i=(int)(alpha1*180/pi); i<171; i++) |
---|
633 | { printf("\n"); |
---|
634 | lam= 1e3*el*(1-cos(i*pi/180))/m; |
---|
635 | printf("%6.1f ", lam); /* in microns */ |
---|
636 | ii = 1; |
---|
637 | sum = sum_p1 = sum_p2 =0; |
---|
638 | for (ii = 1; ii<4; ii++) |
---|
639 | { sum = sum + s_coh_x[i][ii]; |
---|
640 | sum_p1 = sum_p1 + s_coh_x_p1[i][ii]; |
---|
641 | sum_p2 = sum_p2 + s_coh_x_p2[i][ii];} |
---|
642 | s_coh_x[i][8] = (s_coh_x[i][0] + 2*sum) /7; /* average over phi, from -3 to +3 deg.*/ |
---|
643 | s_coh_x_p1[i][8] = (s_coh_x_p1[i][0] + 2*sum_p1)/7; /* average over phi, from -3 to +3 deg.*/ |
---|
644 | s_coh_x_p2[i][8] = (s_coh_x_p2[i][0] + 2*sum_p2)/7; /* average over phi, from -3 to +3 deg.*/ |
---|
645 | for (ii = 0; ii<9; ii++) |
---|
646 | { printf( "%-4.3e ", s_coh_x[i][ii]);} |
---|
647 | printf(" \n"); printf(" "); |
---|
648 | for (ii = 0; ii<9; ii++) |
---|
649 | { printf( "%-4.3e ", s_coh_x_p1[i][ii]);} |
---|
650 | printf("\n"); printf(" "); |
---|
651 | for (ii = 0; ii<9; ii++) |
---|
652 | { printf( "%-4.3e ", s_coh_x_p2[i][ii]);} |
---|
653 | } |
---|
654 | printf("\n"); |
---|
655 | printf(" Each column above is for a phi step of "); printf("%-3.1f",phistep); printf(" deg. \n"); |
---|
656 | printf(" last column is the average over the 3 values of phi on either side of zero, i.e.\n"); |
---|
657 | printf(" between -%-3.1f & +%-3.1f deg. only, in this case \n",phistep*3, phistep*3); |
---|
658 | printf(" Finished with the calculation of R^2; see above for interpretation \n"); |
---|
659 | printf("\n"); |
---|
660 | /* calculation of coherence factors & of radiated energy */ |
---|
661 | |
---|
662 | printf("\nCalculating energy per bunch, for this specific grating (4cm long) \n"); |
---|
663 | f1=2*pi*qe*qe*1e-7; /* in Joules*/ |
---|
664 | printf("The scale factor f1="); printf("%le",f1); printf(" J.cm \n"); |
---|
665 | f2=100/pow(el,2); /* geometry factor f2 in cm-2; grating period in mm */ |
---|
666 | |
---|
667 | i=(int)(alpha1*180/pi); |
---|
668 | for ( i=(int)(alpha1*180/pi); i<171; i++) |
---|
669 | { thet[i]=i; /*scan angles from 10-170 deg. */ |
---|
670 | theta=thet[i]*pi/180; |
---|
671 | lam=el*(1/beta-cos(theta))/m; |
---|
672 | k= 2*pi/lam; |
---|
673 | ii=0; |
---|
674 | for (ii=0; ii<8; ii++) |
---|
675 | {phi0=phistep*ii*pi/180; /* phi scanned in 7 steps, each of size=phistep*/ |
---|
676 | kx=k*sin(theta)*cos(phi0); |
---|
677 | ky=k*sin(theta)*sin(phi0); |
---|
678 | lam_e=lam*beta*gam/(2*pi*sqrt(1+pow((beta*gam*sin(theta)*sin(phi0)),2))); |
---|
679 | /* find the angular distribution factor f4 */ |
---|
680 | f4[i]=pow(m,2)*pow((beta/(1-beta*cos(theta))),3); |
---|
681 | om=2*pi*0.3/lam; /*c=0.3 mm/ps, hence om is in ps^-1*/ |
---|
682 | /*calculate longitudinal coherence effects*/ |
---|
683 | if (pulse_shp==0) |
---|
684 | { |
---|
685 | u=0; /*not necessary, but no coh. for DC beam*/ |
---|
686 | f_z_c=0; |
---|
687 | f_z_s=0; |
---|
688 | } |
---|
689 | else if (pulse_shp==1) |
---|
690 | {u1= om*sigma_t; |
---|
691 | u2= om*sigma_t*epsl; |
---|
692 | temp1=0; temp2=0; |
---|
693 | f_z_c=(exp(-(u1*u1)/2)+ epsl*exp(-(u2*u2)/2))/(1+epsl); /* yes, it is correct!*/ |
---|
694 | /* degenerate hypergeometric function below is taken from G&R page 480 */ |
---|
695 | jj=1; |
---|
696 | hyperg[0]=1; |
---|
697 | for (jj==1; jj<201; jj++) |
---|
698 | {ratio= (u1*u1/2)*(2*jj-1)/(jj*(2*jj+1)); |
---|
699 | hyperg[jj]= hyperg[jj-1]*ratio; |
---|
700 | } |
---|
701 | jj=1; |
---|
702 | sum=1; |
---|
703 | for (jj==1; jj<201; jj++) |
---|
704 | { dumsum= sum+ hyperg[jj]; |
---|
705 | sum=dumsum; |
---|
706 | } |
---|
707 | temp1= fabs(sum*u1*sigma_t*exp(-u1*u1/2)); |
---|
708 | /* now for the trailing part of the bunch */ |
---|
709 | jj=1; |
---|
710 | hyperg[0]=1; |
---|
711 | for (jj==1; jj<201; jj++) |
---|
712 | { ratio= (u2*u2/2)*(2*jj-1)/(jj*(2*jj+1)); |
---|
713 | hyperg[jj]= hyperg[jj-1]*ratio; |
---|
714 | } |
---|
715 | jj=1; |
---|
716 | sum=1; |
---|
717 | for (jj==1; jj<201; jj++) |
---|
718 | { dumsum= sum+ hyperg[jj]; |
---|
719 | sum=dumsum; |
---|
720 | } |
---|
721 | temp2= fabs(sum*u2*sigma_t*epsl*exp(-u2*u2/2)); |
---|
722 | f_z_s=2*(temp1-temp2)/(sqrt(2*pi)*sigma_t*(1+epsl)); /* including nor. factor of 2/sqrt(2pi).. */ |
---|
723 | } |
---|
724 | else if (pulse_shp==2) |
---|
725 | { |
---|
726 | u=om*sigma_t; |
---|
727 | f_z_c=2*(1-(epsl*cos(u)+cos(epsl*u))/(1+epsl))/(epsl*u*u); |
---|
728 | f_z_s=2*((sin(epsl*u)-epsl*sin(u))/(1+epsl))/(epsl*u*u); |
---|
729 | } |
---|
730 | /* integrate to edge of pulse*/ |
---|
731 | else if (pulse_shp==3) |
---|
732 | { |
---|
733 | u=om*sigma_t; |
---|
734 | f_z_c=(1+epsl*u*u)/(1+u*u)/(1+pow((epsl*u),2)); |
---|
735 | f_z_s=(epsl*u-u)/(1+u*u)/(1+pow((epsl*u),2)); |
---|
736 | } |
---|
737 | /*integr. to infinity, so only sigma_t appears in result*/ |
---|
738 | else if (pulse_shp==4) |
---|
739 | { |
---|
740 | u=om*sigma_t; |
---|
741 | f_z_c=3*(sin(u)+sin(epsl*u)/(epsl)-u*cos(u)-u*cos(epsl*u)/(epsl))/((1+epsl)*pow(u,3)); |
---|
742 | f_z_s=3*(1/(epsl*epsl)-1+cos(u)-cos(epsl*u)/(epsl*epsl)+u*sin(u)-u*sin(epsl*u)/(epsl))/((1+epsl)*pow(u,3)); |
---|
743 | /* integrate to edge of pulse i.e. sigma_t*/ |
---|
744 | } |
---|
745 | else if (pulse_shp==5) |
---|
746 | { |
---|
747 | u=om*sigma_t; |
---|
748 | f_z_c=pi*pi*cos(u)/(pi*pi-4*u*u); |
---|
749 | f_z_s=0; |
---|
750 | /*must integrate to edge of pulse, i.e. sigma_t*/ |
---|
751 | } |
---|
752 | else if (pulse_shp==6) |
---|
753 | { |
---|
754 | u=om*sigma_t; |
---|
755 | f_z_c=((1-u*u)*cos(u)+2*u*sin(u))/(pow((1+u*u),2)); |
---|
756 | f_z_s=0; |
---|
757 | /*must integrate from -inf. to edge of pulse, i.e. sigma_t*/ |
---|
758 | } |
---|
759 | else if (pulse_shp==7) /*numerical integr. for Lorentz shape*/ |
---|
760 | { |
---|
761 | /* a and b deal with asymmetric Lorentz shape. a is for leading edge*/ |
---|
762 | /* norm. factor is 2/pi/(a+b) */ |
---|
763 | f_z_c=(b*exp(-om*b)+a*exp(-om*a))/(apb); |
---|
764 | WT_FUNC=Nag_Sine; |
---|
765 | d01asc(LOR1,0.0,om,WT_FUNC,50,500,1e-3,&res1,&abserr,&QPSUB1,&fail); |
---|
766 | d01asc(LOR2,0.0,om,WT_FUNC,50,500,1e-3,&res2,&abserr,&QPSUB2,&fail); |
---|
767 | f_z_s=2*(b*b*res2-a*a*res1)/(pi*(apb)); |
---|
768 | NAG_FREE(QPSUB1.interval_error); |
---|
769 | NAG_FREE(QPSUB1.interval_result); |
---|
770 | NAG_FREE(QPSUB1.interval_flag); |
---|
771 | NAG_FREE(QPSUB2.interval_error); |
---|
772 | NAG_FREE(QPSUB2.interval_result); |
---|
773 | NAG_FREE(QPSUB2.interval_flag); |
---|
774 | } |
---|
775 | else if (pulse_shp==8) /* multiple Gaussian, normalised, modified 15/3/08*/ |
---|
776 | {sigma_teff = g1*sigma_t1+ g2*sigma_t2+g3*sigma_t3+ g4*sigma_t4+g5*sigma_t5; |
---|
777 | u1=om*sigma_t1; u2=om*sigma_t2; |
---|
778 | u3=om*sigma_t3; u4=om*sigma_t4; u5=om*sigma_t5; |
---|
779 | f_z_c=sqrt(2*pi)*(g1*sigma_t1*exp(-u1*u1/2)*cos(mu1*om)+ |
---|
780 | g2*sigma_t2*exp(-u2*u2/2)*cos(mu2*om)+ |
---|
781 | g3*sigma_t3*exp(-u3*u3/2)*cos(mu3*om)+ |
---|
782 | g4*sigma_t4*exp(-u4*u4/2)*cos(mu4*om)+ |
---|
783 | g5*sigma_t5*exp(-u5*u5/2)*cos(mu5*om))/tot_int; |
---|
784 | f_z_s=sqrt(2*pi)*(g1*sigma_t1*exp(-u1*u1/2)*sin(mu1*om)+ |
---|
785 | g2*sigma_t2*exp(-u2*u2/2)*sin(mu2*om)+ |
---|
786 | g3*sigma_t3*exp(-u3*u3/2)*sin(mu3*om)+ |
---|
787 | g4*sigma_t4*exp(-u4*u4/2)*sin(mu4*om)+ |
---|
788 | g5*sigma_t5*exp(-u5*u5/2)*sin(mu5*om))/tot_int; |
---|
789 | } |
---|
790 | |
---|
791 | s_coh_long=pow(f_z_c,2)+pow(f_z_s,2); /*finished with long. coherence factor*/ |
---|
792 | s_coh_z[i][ii]= s_coh_long; /*store in separate array for future use*/ |
---|
793 | f_y_c=exp(-pow((ky*sigma_y),2)/2); /* NB integration here is to infinity; does it matter?? */ |
---|
794 | f_y_s=0.0; |
---|
795 | s_coh_xy= (f_y_c*f_y_c+f_y_s*f_y_s)*s_coh_x[i][ii]; |
---|
796 | s_coh_xy_p1= (f_y_c*f_y_c+f_y_s*f_y_s)*s_coh_x_p1[i][ii]; |
---|
797 | s_coh_xy_p2= (f_y_c*f_y_c+f_y_s*f_y_s)*s_coh_x_p2[i][ii]; |
---|
798 | s_coh_tr[i][ii]= s_coh_xy; /*transverse coherence factor, in new array*/ |
---|
799 | s_coh_tr_p1[i][ii]= s_coh_xy_p1; /*transverse coherence factor, pol. 1, in new array*/ |
---|
800 | s_coh_tr_p2[i][ii]= s_coh_xy_p2; /*transverse coherence factor, pol.2, in new array*/ |
---|
801 | f_b[i][ii]= n_e*s_inc_x[i][ii]+ n_e*n_e*s_coh_xy*s_coh_long; /*overall effect of Ne electrons in a bunch*/ |
---|
802 | f_b_p1[i][ii]= n_e*s_inc_x_p1[i][ii]+ n_e*n_e*s_coh_xy_p1*s_coh_long; /*...in pol.1..*/ |
---|
803 | f_b_p2[i][ii]= n_e*s_inc_x_p2[i][ii]+ n_e*n_e*s_coh_xy_p2*s_coh_long; /*... and pol.2..*/ |
---|
804 | bin[i][0]=i; |
---|
805 | bin[i][ii+1]= f1*f2*f4[i]*(n_e*s_inc_x[i][ii]+ n_e*n_e*s_coh_xy*s_coh_long); |
---|
806 | bin_p1[i][ii+1]= f1*f2*f4[i]*(n_e*s_inc_x_p1[i][ii]+ n_e*n_e*s_coh_xy_p1*s_coh_long); |
---|
807 | bin_p2[i][ii+1]= f1*f2*f4[i]*(n_e*s_inc_x_p2[i][ii]+ n_e*n_e*s_coh_xy_p2*s_coh_long); |
---|
808 | bin_single[i][0]= i; /* store for single electron yield */ |
---|
809 | bin_single[i][ii+1]= f1*f2*f4[i]*single_elec[i][ii]; |
---|
810 | bin_single_p1[i][ii+1]= f1*f2*f4[i]*single_elec_p1[i][ii]; |
---|
811 | bin_single_p2[i][ii+1]= f1*f2*f4[i]*single_elec_p2[i][ii]; |
---|
812 | /*NB!! in col.1 data for phi=0, in 2 for phi=1, in 8 for phi=7*/ |
---|
813 | } /* end of ii (phi) scan */ |
---|
814 | } /* end of i (theta) scan */ |
---|
815 | printf("\n"); printf("\n"); |
---|
816 | |
---|
817 | /* printf("For Phi0=2 ONLY \n"); |
---|
818 | printf(" Angle Lambda S_inc S_coh_z S_coh_tr Bunch_factor Energy \n"); |
---|
819 | printf(" (mm) (J/sr) \n"); |
---|
820 | i=(int)(alpha1*180/pi); |
---|
821 | for ( i=(int)(alpha1*180/pi); i<171; i++) |
---|
822 | { ii=2; |
---|
823 | thet[i]=i; /*scan angles up to 170 deg. */ |
---|
824 | /* theta=thet[i]*pi/180; |
---|
825 | lam=el*(1/beta-cos(theta))/m; |
---|
826 | printf("%5.1f",thet[i]); printf(" "); |
---|
827 | printf("%4.3f",lam); |
---|
828 | printf(" "); printf("%-4.3e",s_inc_x[i][ii]); printf(" "); |
---|
829 | printf("%-4.3e",s_coh_z[i][ii]); printf(" "); printf("%-4.3e",s_coh_tr[i][ii]); printf(" "); |
---|
830 | printf("%-4.3e",f_b[i][ii]); printf(" "); printf("%-4.3e",bin[i][ii+1]); |
---|
831 | printf("\n"); |
---|
832 | printf(" "); printf("%-4.3e",s_inc_x_p1[i][ii]); printf(" "); |
---|
833 | printf("%-4.3e",s_coh_z[i][ii]); printf(" "); printf("%-4.3e",s_coh_tr_p1[i][ii]); printf(" "); |
---|
834 | printf("%-4.3e",f_b_p1[i][ii]); printf(" "); printf("%-4.3e",bin_p1[i][ii+1]); |
---|
835 | printf("\n"); |
---|
836 | printf(" "); printf("%-4.3e",s_inc_x_p2[i][ii]); printf(" "); |
---|
837 | printf("%-4.3e",s_coh_z[i][ii]); printf(" "); printf("%-4.3e",s_coh_tr_p2[i][ii]); printf(" "); |
---|
838 | printf("%-4.3e",f_b_p2[i][ii]); printf(" "); printf("%-4.3e",bin_p2[i][ii+1]); |
---|
839 | printf("\n"); |
---|
840 | } |
---|
841 | printf(" Angle Lambda S_inc S_coh_z S_coh_tr Bunch_factor Energy \n"); |
---|
842 | printf(" (mm) (J/sr) \n"); |
---|
843 | printf("The above apply for Phi2=0 ONLY \n"); |
---|
844 | printf("\n"); |
---|
845 | printf("the following is a list of the INCOHERENT & COHERENT parts of the x-integral; INCOHERENT is first"); |
---|
846 | printf("\n"); |
---|
847 | i=(int)(alpha1*180/pi); |
---|
848 | for ( i=(int)(alpha1*180/pi); i<171; i++) |
---|
849 | { printf("\n"); |
---|
850 | printf("%3d ", i); |
---|
851 | ii=0; |
---|
852 | for (ii=0; ii<8; ii++) |
---|
853 | { printf( "%.3e ", s_inc_x[i][ii]); } |
---|
854 | printf("\n"); printf(" "); |
---|
855 | ii=0; |
---|
856 | for (ii=0; ii<8; ii++) |
---|
857 | { printf( "%.3e ", s_coh_x[i][ii]); } |
---|
858 | } |
---|
859 | printf("\n");*/ |
---|
860 | |
---|
861 | |
---|
862 | printf(" List of single-electron yields \n"); |
---|
863 | i=(int)(alpha1*180/pi); |
---|
864 | for ( i=(int)(alpha1*180/pi); i<171; i++) |
---|
865 | { printf("\n"); |
---|
866 | printf("%3d ", i); |
---|
867 | ii=1; |
---|
868 | for (ii=1; ii<9; ii++) |
---|
869 | { printf( "%.3e ", bin_single[i][ii]); } |
---|
870 | |
---|
871 | } |
---|
872 | printf("\n"); |
---|
873 | |
---|
874 | i=0; |
---|
875 | for (i=0; i<151; i++) |
---|
876 | {dum70[i][0]=i+20; |
---|
877 | dum70[i][1]=bin[i+20][1]; /*copy to DUM70 but only from theta=20 */ |
---|
878 | } |
---|
879 | l=151; /*start of sorting section, N=no. of rows*/ |
---|
880 | n=1; /* no. of columns, i.e. 2 here*/ |
---|
881 | ic=1; /* column index for sorting, i.e. 2nd here*/ |
---|
882 | i=0; |
---|
883 | for (i=0; i<l; i++) |
---|
884 | { j=i+1; |
---|
885 | for (j=i+1; j<l+1; j++) |
---|
886 | { if (dum70[i][ic]<dum70[j][ic]) |
---|
887 | { jj=0; |
---|
888 | for (jj=0; jj<n+1; jj++) |
---|
889 | { dumm[jj]=dum70[i][jj]; |
---|
890 | dum70[i][jj]=dum70[j][jj]; |
---|
891 | dum70[j][jj]=dumm[jj]; |
---|
892 | } |
---|
893 | } |
---|
894 | } |
---|
895 | } |
---|
896 | printf("\n"); |
---|
897 | printf(" Maximum for this case =%-4.3e",dum70[0][1]); |
---|
898 | printf(" J/sr/cm at theta =%6.2f\n",dum70[0][0]); |
---|
899 | |
---|
900 | /*rough estimate of full width at half max. for the power distribution*/ |
---|
901 | i=0; |
---|
902 | for (i=0; i<151; i++) |
---|
903 | { |
---|
904 | if (dum70[i][1]>=0.5*dum70[1][1]) |
---|
905 | { |
---|
906 | bin_1[i][0]=dum70[i][0]; |
---|
907 | bin_1[i][1]=dum70[i][1]; |
---|
908 | } |
---|
909 | else |
---|
910 | { |
---|
911 | bin_1[i][0]=0; |
---|
912 | bin_1[i][1]=0; |
---|
913 | } |
---|
914 | } |
---|
915 | l=151; /*start of sorting section, N=no. of rows*/ |
---|
916 | n=1; /* no. of columns, i.e. 2 here*/ |
---|
917 | ic=0; /* column index for sorting, i.e. 1st here*/ |
---|
918 | i=0; |
---|
919 | for (i=0; i<l; i++) |
---|
920 | { j=i+1; |
---|
921 | for (j=i+1; j<l+1; j++) |
---|
922 | { if (bin_1[i][ic]<bin_1[j][ic]) |
---|
923 | { jj=0; |
---|
924 | for (jj=0; jj<n+1; jj++) |
---|
925 | { dumm[jj]=bin_1[i][jj]; |
---|
926 | bin_1[i][jj]=bin_1[j][jj]; |
---|
927 | bin_1[j][jj]=dumm[jj]; |
---|
928 | } |
---|
929 | } |
---|
930 | } |
---|
931 | } |
---|
932 | theta_u=bin_1[0][0]; |
---|
933 | i=0; |
---|
934 | for (i=0; i<151; i++) |
---|
935 | { |
---|
936 | if (bin_1[i][1]<=0) goto esc2; |
---|
937 | theta_d=bin_1[i][0]; |
---|
938 | } |
---|
939 | esc2:printf("\n"); |
---|
940 | printf("Rough estimate of full-width, half maximum \n"); |
---|
941 | printf(" Lower limit at= %6.2f\n",theta_d); |
---|
942 | printf(" Upper limit at= %6.2f\n",theta_u); |
---|
943 | printf(" \n"); |
---|
944 | /* rough estimate of total radiated energy, per single electron & per bunch */ |
---|
945 | printf("I will now average over the azimuthal angle phi \n"); |
---|
946 | printf("The min. number of steps in phi, INCLUDING phi=0, is 1; max is 8,"); |
---|
947 | printf(" which would take you up to phi= %4.1f \n", (7*phistep)); |
---|
948 | printf(" Enter number of steps; default is 4, which would correspond to phi=3 (for steps of 1 deg.) \n"); |
---|
949 | scanf("%1i",&philmt); /*!! NB. different definition of philmt in this version!!*/ |
---|
950 | if (philmt<1 || philmt>8) philmt=4; |
---|
951 | printf(" I will average up to phi= %-3.1f ",((philmt-1)*phistep)); printf(" \n"); |
---|
952 | i=(int)(alpha1*180/pi); |
---|
953 | for ( i=(int)(alpha1*180/pi); i<171; i++) /* angles lower than the blaze angle do not matter! */ |
---|
954 | {sum= sum_p1= sum_p2= sum_single= sum_single_p1= sum_single_p2= 0.0; |
---|
955 | ii=2; |
---|
956 | if (philmt==1) |
---|
957 | {sum= bin[i][1]; sum_p1=bin_p1[i][1]; sum_p2=bin_p2[i][1]; sum_single= bin_single[i][1]; |
---|
958 | sum_single_p1= bin_single_p1[i][1]; |
---|
959 | sum_single_p2= bin_single_p2[i][1];} |
---|
960 | else |
---|
961 | for (ii=2; ii<philmt+1; ii++) |
---|
962 | {sum=sum+ 2*bin[i][ii]; /* start from col. index 2(phi=1xphistep) and double the sum to account for +-phi */ |
---|
963 | sum_p1=sum_p1+ 2*bin_p1[i][ii]; /* ditto for pol. 1 */ |
---|
964 | sum_p2=sum_p2+ 2*bin_p2[i][ii]; /* and for pol.2 */ |
---|
965 | sum_single= sum_single+ 2*bin_single[i][ii]; |
---|
966 | sum_single_p1= sum_single_p1+ 2*bin_single_p1[i][ii]; |
---|
967 | sum_single_p2= sum_single_p2+ 2*bin_single_p2[i][ii];} |
---|
968 | if (philmt==1) |
---|
969 | {bin[i][9]= bin[i][1]; |
---|
970 | bin_p1[i][9]= bin_p1[i][1]; |
---|
971 | bin_p2[i][9]= bin_p2[i][1]; |
---|
972 | bin_single[i][9]= bin_single[i][1]; |
---|
973 | bin_single_p1[i][9]= bin_single_p1[i][1]; |
---|
974 | bin_single_p2[i][9]= bin_single_p2[i][1];} |
---|
975 | else |
---|
976 | {bin[i][9]=(bin[i][1]+ sum)/(2*philmt-1); /* col. 9 is the store for AVERAGE dE over all ACCEPTED phi's; in J/sr/cm */ |
---|
977 | bin_p1[i][9]=(bin_p1[i][1]+ sum_p1)/(2*philmt-1); /* ditto for pol. 1 */ |
---|
978 | bin_p2[i][9]=(bin_p2[i][1]+ sum_p2)/(2*philmt-1); /* and for pol. 2 */ |
---|
979 | bin_single[i][9]= (bin_single[i][1]+ sum_single)/(2*philmt-1); |
---|
980 | bin_single_p1[i][9]= (bin_single_p1[i][1]+ sum_single_p1)/(2*philmt-1); |
---|
981 | bin_single_p2[i][9]= (bin_single_p2[i][1]+ sum_single_p2)/(2*philmt-1);} |
---|
982 | } |
---|
983 | printf("\n"); |
---|
984 | |
---|
985 | /* integrate over phi up to the max. value=7 */ |
---|
986 | d_theta=1.745e-2; /* d_theta step is always 1 degree */ |
---|
987 | d_phi= phistep*1.745e-2; /* d_phi expressed in rad */ |
---|
988 | i=(int)(alpha1*180/pi); |
---|
989 | for (i=(int)(alpha1*180/pi); i<171; i++) /* integrate dE over ALL phi & store in 10*/ |
---|
990 | { sum= sum_p1= sum_p2= sum_single= sum_single_p1= sum_single_p2= 0.0; |
---|
991 | ii=1; |
---|
992 | for (ii=2; ii<9; ii++) |
---|
993 | {sum = sum + sin(i*pi/180)*d_theta*d_phi*bin[i][ii]; |
---|
994 | sum_p1 = sum_p1 + sin(i*pi/180)*d_theta*d_phi*bin_p1[i][ii]; /* pol. 1 only */ |
---|
995 | sum_p2 = sum_p2 + sin(i*pi/180)*d_theta*d_phi*bin_p2[i][ii]; /* and pol. 2*/ |
---|
996 | sum_single = sum_single + sin(i*pi/180)*d_theta*d_phi*bin_single[i][ii]; |
---|
997 | sum_single_p1 = sum_single_p1 + sin(i*pi/180)*d_theta*d_phi*bin_single_p1[i][ii]; |
---|
998 | sum_single_p2 = sum_single_p2 + sin(i*pi/180)*d_theta*d_phi*bin_single_p2[i][ii]; |
---|
999 | } |
---|
1000 | bin[i][10]= sin(i*pi/180)*d_theta*d_phi*bin[i][1]+ 2*sum; |
---|
1001 | /* col. 10 is a temp. store for phi-integration over ALL phi's, i.e. up to phi=7; in J/sr/cm */ |
---|
1002 | /* only used for the estimate of the total radiated energy, see below */ |
---|
1003 | |
---|
1004 | bin_p1[i][10] = sin(i*pi/180)*d_theta*d_phi*bin_p1[i][1]+ 2*sum_p1; /* pol. 1 only */ |
---|
1005 | bin_p2[i][10] = sin(i*pi/180)*d_theta*d_phi*bin_p2[i][1]+ 2*sum_p2; /* and pol. 2 */ |
---|
1006 | bin_single[i][10] = sin(i*pi/180)*d_theta*d_phi*bin_single[i][1]+ 2*sum_single; |
---|
1007 | bin_single_p1[i][10] = sin(i*pi/180)*d_theta*d_phi*bin_single_p1[i][1] + 2*sum_single_p1; |
---|
1008 | bin_single_p2[i][10] = sin(i*pi/180)*d_theta*d_phi*bin_single_p2[i][1] + 2*sum_single_p2; |
---|
1009 | } |
---|
1010 | /* now integrate over theta */ |
---|
1011 | i=(int)(alpha1*180/pi); |
---|
1012 | sum1 = sum2 = sum3 = sum4 = 0.0; |
---|
1013 | for (i=(int)(alpha1*180/pi); i<171; i++) /* now truly in J! The length is in the f2 factor! */ |
---|
1014 | { sum1= sum1 + bin[i][10]; |
---|
1015 | sum3= sum3 + bin_p1[i][10]; |
---|
1016 | sum4= sum4 + bin_p2[i][10]; |
---|
1017 | sum2= sum2 + bin_single[i][10]; |
---|
1018 | } |
---|
1019 | printf("\n"); |
---|
1020 | |
---|
1021 | printf(" For this grating \n"); |
---|
1022 | printf(" Rough estimate of total radiated energy (single electron) over hemisphere = "); |
---|
1023 | printf("%.3e", sum2*0.1*lgr); /* no need for x2, since +-phi has been taken into account*/ |
---|
1024 | printf(" J"); printf("\n"); |
---|
1025 | printf(" Rough estimate of total radiated energy (bunch) over hemisphere = "); |
---|
1026 | printf("%.3e", sum1*0.1*lgr); /* no need for x2, since +-phi has been taken into account*/ |
---|
1027 | printf(" J"); printf("\n"); |
---|
1028 | percent=sum1*0.1*lgr/(n_e_t*1.602177e-19*(gam-1)*0.511*1e6); |
---|
1029 | printf(" This is equal to "); printf("%.3e",percent); printf(" of the energy in the bunch \n"); |
---|
1030 | printf(" Energy in polarization 1= "); printf("%.3e", sum3*0.1*lgr); printf(" J \n"); |
---|
1031 | printf(" Energy in polarization 2= "); printf("%.3e", sum4*0.1*lgr); printf(" J \n"); |
---|
1032 | printf("\n"); printf("\n"); |
---|
1033 | printf("This case had the following parameters:\n"); |
---|
1034 | printf("gamma= %-5.1f \n",gam); |
---|
1035 | printf("beam height(mm)= %-6.3f \n",x0); |
---|
1036 | printf("pulse length (FWHM,ps)= %-6.3f \n",pulse_len); |
---|
1037 | printf("total number of electrons in the bunch= %-6.3e \n",n_e_t); |
---|
1038 | printf("pulse shape code is= %1i \n",pulse_shp); |
---|
1039 | if (pulse_shp==8) printf("no. of Gaussians= %1i \n", ngauss); |
---|
1040 | else |
---|
1041 | printf("asymmetry factor= %-3.1f \n",epsl); |
---|
1042 | if (pulse_shp==7) |
---|
1043 | {apb=pulse_len; a=apb/(1+epsl); |
---|
1044 | printf("for Lorentz shape a+b= %-5.3f ",apb); |
---|
1045 | printf(" ps and leading edge fraction = %-5.3f \n",a);} |
---|
1046 | printf("grating period(mm)= %-5.3f \n",el); |
---|
1047 | printf("emission order= %1i \n",m); |
---|
1048 | printf("PHI step= %-4.2f \n",phistep); |
---|
1049 | printf("limiting value of phi= %-3.1f ",(phistep*(philmt-1))); printf(" deg. \n"); |
---|
1050 | printf("\n"); |
---|
1051 | |
---|
1052 | flag=0; |
---|
1053 | if (!flag) printf("Do you want a .txt file with the dE values, bunch (p1,p2 & total) & single-electron? yes(y) or no(n) "); |
---|
1054 | line[0] = 'x'; |
---|
1055 | do { |
---|
1056 | if (line[0] != 'x') printf(" "); |
---|
1057 | if (fgets(line, 200, stdin) == NULL) { printf("\nErr: fgets.\n"); exit(1); } |
---|
1058 | } |
---|
1059 | while (line[0] != 'y' && line[0] != 'n'); |
---|
1060 | ans2 = line[0]; |
---|
1061 | if (ans2=='n') goto esc3; |
---|
1062 | printf("Name of output file?? "); |
---|
1063 | scanf("%s",&outfile_name[0]); |
---|
1064 | out_file=fopen(outfile_name,"w"); |
---|
1065 | i=(int)(alpha1*180/pi); |
---|
1066 | for (i=(int)(alpha1*180/pi); i<171; i++) |
---|
1067 | { fprintf(out_file,"\n"); |
---|
1068 | ii=0; |
---|
1069 | for (ii=0; ii<11; ii++) |
---|
1070 | {if (ii==0) fprintf(out_file, "%-5.1f ", bin[i][ii]); |
---|
1071 | else fprintf(out_file, " %.3e ", bin[i][ii]); } |
---|
1072 | fprintf(out_file, "\n"); fprintf(out_file, " "); |
---|
1073 | ii=1; |
---|
1074 | for (ii==1; ii<11; ii++) |
---|
1075 | { fprintf(out_file, " %.3e ", bin_p1[i][ii]);} |
---|
1076 | fprintf(out_file, "\n"); fprintf(out_file, " "); |
---|
1077 | ii=1; |
---|
1078 | for (ii==1; ii<11; ii++) |
---|
1079 | { fprintf(out_file, " %.3e ", bin_p2[i][ii]);} |
---|
1080 | } |
---|
1081 | fprintf(out_file,"\n"); fprintf(out_file,"\n"); |
---|
1082 | fprintf(out_file, "single-electron data below, also in J \n"); |
---|
1083 | fprintf(out_file,"\n"); |
---|
1084 | |
---|
1085 | i=(int)(alpha1*180/pi); |
---|
1086 | for (i=(int)(alpha1*180/pi); i<171; i++) |
---|
1087 | { fprintf(out_file,"\n"); |
---|
1088 | ii=0; |
---|
1089 | for (ii=0; ii<11; ii++) |
---|
1090 | {if (ii==0) fprintf(out_file, "%-5.1f ", bin_single[i][ii]); |
---|
1091 | else fprintf(out_file, " %.3e ", bin_single[i][ii]); } |
---|
1092 | fprintf(out_file, "\n"); fprintf(out_file, " "); |
---|
1093 | ii=1; |
---|
1094 | for (ii==1; ii<11; ii++) |
---|
1095 | { fprintf(out_file, " %.3e ", bin_single_p1[i][ii]);} |
---|
1096 | fprintf(out_file, "\n"); fprintf(out_file, " "); |
---|
1097 | ii=1; |
---|
1098 | for (ii==1; ii<11; ii++) |
---|
1099 | { fprintf(out_file, " %.3e ", bin_single_p2[i][ii]);} |
---|
1100 | } |
---|
1101 | fprintf(out_file,"\n"); |
---|
1102 | fprintf(out_file,"gamma= %6.1f \n",gam); |
---|
1103 | fprintf(out_file,"x0(mm)= %6.3f \n",x0); |
---|
1104 | fprintf(out_file,"fwhm(ps)= %6.3f \n",pulse_len); |
---|
1105 | fprintf(out_file,"no. of electrons= %-4.3e \n",n_e_t); |
---|
1106 | fprintf(out_file,"pulse shape= %1d \n",pulse_shp); |
---|
1107 | if (pulse_shp==8) fprintf(out_file,"no. of Gaussians= %1i \n", ngauss); |
---|
1108 | else |
---|
1109 | fprintf(out_file,"asymmetry factor= %-3.1f \n",epsl); |
---|
1110 | fprintf(out_file,"grating period(mm)= %-5.3f \n",el); |
---|
1111 | fprintf(out_file,"order= %1d \n",m); |
---|
1112 | fprintf(out_file,"phi step= %3.1f \n",phistep); |
---|
1113 | fprintf(out_file,"no. of steps in phi= %2d \n",philmt); |
---|
1114 | fclose(out_file); |
---|
1115 | esc3: printf("\n"); |
---|
1116 | |
---|
1117 | printf("Angle Theta2 Theta1 energy/bunch energy(p1) energy(p2) polarisation\n"); |
---|
1118 | printf(" J J J\n"); |
---|
1119 | i=(int)(alpha1*180/pi); |
---|
1120 | for (i=(int)(alpha1*180/pi); i<171; i++) |
---|
1121 | { |
---|
1122 | thet[i]=i; |
---|
1123 | theta=thet[i]*pi/180; |
---|
1124 | lgr_e = fabs(21/sin(theta)); /* variable effective grating length, determined basically by |
---|
1125 | the cone entrance diameter of 21mm*/ |
---|
1126 | /*lgr_e = lgr;*/ |
---|
1127 | |
---|
1128 | /*thet_1 and thet_2 below are the max. possible angles for a given observation angle*/ |
---|
1129 | /* NB. !!!!definition changed from previous versions of the code!!!*/ |
---|
1130 | /* now determined by the Winston cone acceptance angle */ |
---|
1131 | thet_1 = theta + 6*pi/180; |
---|
1132 | thet_2 = theta - 6*pi/180; |
---|
1133 | if (thet_1<0) thet_1=pi+thet_1; |
---|
1134 | if (thet_2<0) thet_2=pi+thet_2; |
---|
1135 | int_thet_1=(int)180*thet_1/pi; |
---|
1136 | int_thet_2=(int)180*thet_2/pi; |
---|
1137 | if (180*thet_1/pi-int_thet_1>=0.5) int_thet_1++; |
---|
1138 | if (180*thet_2/pi-int_thet_2>=0.5) int_thet_2++; |
---|
1139 | outbin[i][0]=theta*180/pi; outbin_single[i][0]=theta*180/pi; |
---|
1140 | outbin[i][1]=int_thet_2; outbin_single[i][1]=int_thet_2; |
---|
1141 | outbin[i][2]=int_thet_1; outbin_single[i][2]=int_thet_1; |
---|
1142 | /* n0 is total no. of theta angles accepted, modified May 08; must be 13 values */ |
---|
1143 | n0=int_thet_1-int_thet_2; |
---|
1144 | if (n0<=0) n0=-n0; |
---|
1145 | if (int_thet_1<=int_thet_2) jstart= int_thet_1; |
---|
1146 | else |
---|
1147 | jstart=int_thet_2; |
---|
1148 | |
---|
1149 | /* now average over accepted values of theta the already averaged values over phi */ |
---|
1150 | /* remember, col.9 of bin holds the previous averaging over phi*/ |
---|
1151 | sum_th= sum_th_p1= sum_th_p2= sum_th_single= sum_th_single_p1= sum_th_single_p2= 0.0; |
---|
1152 | j = 1; |
---|
1153 | for (j = 1; j<n0+2; j++) |
---|
1154 | { sum_th = sum_th + bin[jstart+j-1][9]; |
---|
1155 | sum_th_p1 = sum_th_p1 + bin_p1[jstart+j-1][9]; |
---|
1156 | sum_th_p2 = sum_th_p2 + bin_p2[jstart+j-1][9]; |
---|
1157 | sum_th_single = sum_th_single + bin_single[jstart+j-1][9]; |
---|
1158 | sum_th_single_p1 = sum_th_single_p1 + bin_single_p1[jstart+j-1][9]; |
---|
1159 | sum_th_single_p2 = sum_th_single_p2 + bin_single_p2[jstart+j-1][9]; |
---|
1160 | }; |
---|
1161 | |
---|
1162 | aven = sum_th/(n0 + 1); /*have averaged over accepted theta & phi angles; in J/sr/cm*/ |
---|
1163 | aven_p1 = sum_th_p1/(n0 + 1); /*ditto for pol.1*/ |
---|
1164 | aven_p2 = sum_th_p2/(n0 + 1); /*and for pol.2 */ |
---|
1165 | aven_single = sum_th_single/(n0 + 1); |
---|
1166 | aven_single_p1 = sum_th_single_p1/(n0 + 1); |
---|
1167 | aven_single_p2 = sum_th_single_p2/(n0 + 1); |
---|
1168 | |
---|
1169 | sum_tr = aven * 0.1 * lgr_e * domeg_m; /*variable grating length, Sept 12; |
---|
1170 | now in J*/ |
---|
1171 | sum_tr_single = aven_single * 0.1 * lgr_e * domeg_m; |
---|
1172 | |
---|
1173 | sum_tr_p1 = aven_p1 * 0.1 * lgr_e * domeg_m; |
---|
1174 | sum_tr_p2 = aven_p2 * 0.1 * lgr_e * domeg_m; |
---|
1175 | sum_tr_single = aven_single * 0.1 * lgr_e * domeg_m; |
---|
1176 | sum_tr_single_p1 = aven_single_p1 * 0.1 * lgr_e * domeg_m; |
---|
1177 | sum_tr_single_p2 = aven_single_p2 * 0.1 * lgr_e * domeg_m; |
---|
1178 | outbin[i][3] = aven; |
---|
1179 | outbin_single[i][3] = aven_single; |
---|
1180 | outbin[i][4] = sum_tr_p1; |
---|
1181 | outbin_single[i][4] = sum_tr_single_p1; /* store single-electr. numbers but don't print */ |
---|
1182 | outbin[i][5] = sum_tr_p2; |
---|
1183 | outbin_single[i][5] = sum_tr_single_p2; |
---|
1184 | outbin[i][6] = (sum_tr_p1-sum_tr_p2)/(sum_tr_p1+sum_tr_p2); |
---|
1185 | outbin_single[i][6] = (sum_tr_single_p1-sum_tr_single_p2)/(sum_tr_single_p1+sum_tr_single_p2); |
---|
1186 | printf("%-5.1f",theta*180/pi); printf(" "); printf("%5d",int_thet_2); printf(" "); |
---|
1187 | printf("%5d",int_thet_1); printf(" "); |
---|
1188 | printf("%-6.3e ",sum_tr); printf(" "); |
---|
1189 | printf("%-6.3e",sum_tr_p1); printf(" "); printf("%-6.3e",sum_tr_p2); printf(" "); |
---|
1190 | printf("%-6.3e\n", (sum_tr_p1-sum_tr_p2)/(sum_tr_p1+sum_tr_p2)); |
---|
1191 | |
---|
1192 | }; |
---|
1193 | printf("Angle Theta2 Theta1 energy/bunch energy(p1) energy(p2) polarisation\n"); |
---|
1194 | printf(" J J J\n"); |
---|
1195 | printf("\n"); |
---|
1196 | flag=0; |
---|
1197 | if (!flag) printf("Do you want an output .txt file with the collected energy values, bunch & single electron? yes(y) or no(n) "); |
---|
1198 | line[0] = 'x'; |
---|
1199 | do { |
---|
1200 | if (line[0] != 'x') printf(" "); |
---|
1201 | if (fgets(line, 200, stdin) == NULL) { printf("\nErr: fgets.\n"); exit(1); } |
---|
1202 | } |
---|
1203 | while (line[0] != 'y' && line[0] != 'n'); |
---|
1204 | ans2 = line[0]; |
---|
1205 | if (ans2=='y') |
---|
1206 | { printf("Name of output file?? "); |
---|
1207 | scanf("%s",&outfile_name[0]); |
---|
1208 | out_file=fopen(outfile_name,"w"); |
---|
1209 | i=(int)(alpha1*180/pi); |
---|
1210 | for (i=(int)(alpha1*180/pi); i<171; i++) |
---|
1211 | { fprintf(out_file,"\n"); ii=0; |
---|
1212 | for (ii=0; ii<7; ii++) |
---|
1213 | {if (ii<3) fprintf(out_file, " %-5.1f \t", outbin[i][ii]); |
---|
1214 | else if (ii==3) fprintf(out_file, " %-6.3e \t", outbin[i][ii]); |
---|
1215 | else fprintf(out_file, " %-6.3e ", outbin[i][ii]); } |
---|
1216 | } |
---|
1217 | fprintf(out_file,"\n"); fprintf(out_file,"\n"); |
---|
1218 | fprintf(out_file, " single electron data follow below \n"); |
---|
1219 | fprintf(out_file,"\n"); |
---|
1220 | i=(int)(alpha1*180/pi); |
---|
1221 | for (i=(int)(alpha1*180/pi); i<171; i++) |
---|
1222 | { fprintf(out_file,"\n"); ii=0; |
---|
1223 | for (ii=0; ii<7; ii++) |
---|
1224 | {if (ii<3) fprintf(out_file, " %-5.1f \t", outbin_single[i][ii]); |
---|
1225 | else if (ii==3) fprintf(out_file, " %-6.3e \t", outbin_single[i][ii]); |
---|
1226 | else fprintf(out_file, " %-6.3e ", outbin_single[i][ii]); } |
---|
1227 | } |
---|
1228 | fprintf(out_file,"\n"); |
---|
1229 | fprintf(out_file," This case had the following parameters: \n"); |
---|
1230 | fprintf(out_file,"gamma= %6.1f \n",gam); |
---|
1231 | fprintf(out_file,"beam height(mm)= %6.3f \n",x0); |
---|
1232 | fprintf(out_file,"pulse length (FWHM, ps)= %6.3f \n",pulse_len); |
---|
1233 | fprintf(out_file,"no. of electrons= %-4.3e \n",n_e_t); |
---|
1234 | fprintf(out_file,"pulse shape code= %1d \n",pulse_shp); |
---|
1235 | if (pulse_shp==8) fprintf(out_file,"no. of Gaussians= %1d \n",ngauss); |
---|
1236 | else if (pulse_shp==7) |
---|
1237 | {fprintf(out_file,"Lorentz shape with a+b= %-5.3f ",apb); |
---|
1238 | fprintf(out_file, "and leading fraction of bunch= %-4.2f \n",a);} |
---|
1239 | else |
---|
1240 | fprintf(out_file,"asymmetry factor= %-3.1f \n",epsl); |
---|
1241 | fprintf(out_file,"grating period(mm)= %-5.3f \n",el); |
---|
1242 | fprintf(out_file,"order= %1d \n",m); |
---|
1243 | fprintf(out_file,"phi step= %3.1f \n",phistep); |
---|
1244 | fprintf(out_file,"limiting phi= %3.1f ",((philmt-1)*phistep)); fprintf(out_file," deg.\n"); |
---|
1245 | fprintf(out_file,"first column= emission angle theta \n"); |
---|
1246 | fprintf(out_file,"second column= theta2 \n"); |
---|
1247 | fprintf(out_file,"third column= theta1 \n"); |
---|
1248 | fprintf(out_file,"fourth column= energy per bunch J, total \n"); |
---|
1249 | fprintf(out_file,"fifth column= energy per bunch (J) in P1 \n"); |
---|
1250 | fprintf(out_file,"sixth column= energy per bunch (J) in P2\n"); |
---|
1251 | fprintf(out_file,"seventh column= degree of polarisation \n"); |
---|
1252 | fprintf(out_file,"total radiated energy per bunch(J)=%-4.3e \n", sum1); |
---|
1253 | fprintf(out_file,"total radiated energy single electron(J)=%-4.3e \n", sum2); |
---|
1254 | fprintf(out_file, " \n "); |
---|
1255 | fclose(out_file); |
---|
1256 | printf(" Done it, goodbye \n"); |
---|
1257 | } |
---|
1258 | else |
---|
1259 | { printf(" OK, re-entry options? "); |
---|
1260 | printf("\n");} |
---|
1261 | /*re-entry options*/ |
---|
1262 | printf("Type:\n"); |
---|
1263 | printf("1 for new beam position \n"); |
---|
1264 | printf("2 for new sigma_x \n"); |
---|
1265 | printf("3 for new bunch shape\n"); |
---|
1266 | printf("4 for output file for CAUCHY & exit\n"); |
---|
1267 | printf("anything else to exit this programme\n"); |
---|
1268 | scanf("%d",&ans1); |
---|
1269 | if (ans1==1) goto esc20; |
---|
1270 | else if (ans1==2) goto esc24; |
---|
1271 | else if (ans1==3) goto esc1; |
---|
1272 | else if (ans1==4) |
---|
1273 | {printf("Name of output file?? "); |
---|
1274 | scanf("%s",&outfile_name[0]); |
---|
1275 | out_file=fopen(outfile_name,"w"); |
---|
1276 | i=(int)(alpha1*180/pi); |
---|
1277 | for (i=(int)(alpha1*180/pi); i<171; i++) |
---|
1278 | { if ( (outbin[i][0] - ((int)(outbin[i][0] / 10))*10 ) == 0 && outbin[i][0] >= 40 && outbin[i][0] <= 140 ) |
---|
1279 | { fprintf(out_file, " %-5.1f\t%-6.3e\t1\n", outbin[i][0], outbin[i][4]); |
---|
1280 | } |
---|
1281 | } |
---|
1282 | fprintf( out_file, "%-5.3f\n", el ); |
---|
1283 | fclose(out_file);} |
---|
1284 | esc999: printf(" End of programme\n"); |
---|
1285 | |
---|
1286 | return (0); |
---|
1287 | } |
---|
1288 | |
---|
1289 | complex Ccreate(double r, double i) |
---|
1290 | { |
---|
1291 | complex t; |
---|
1292 | t.r = r; |
---|
1293 | t.i = i; |
---|
1294 | return t; |
---|
1295 | } |
---|
1296 | |
---|
1297 | complex Cadd(complex a, complex b) |
---|
1298 | { |
---|
1299 | complex t; |
---|
1300 | t.r = a.r + b.r; |
---|
1301 | t.i = a.i + b.i; |
---|
1302 | return t; |
---|
1303 | } |
---|
1304 | |
---|
1305 | complex Cmultiply(complex a, complex b) |
---|
1306 | { |
---|
1307 | complex t; |
---|
1308 | t.r = a.r * b.r - a.i * b.i; |
---|
1309 | t.i = a.i * b.r + a.r * b.i; |
---|
1310 | return t; |
---|
1311 | } |
---|
1312 | |
---|
1313 | complex Cdivide(complex a, complex b) |
---|
1314 | { |
---|
1315 | complex t; |
---|
1316 | double n; |
---|
1317 | n = b.r * b.r + b.i * b.i; |
---|
1318 | t.r = (a.r * b.r + a.i * b.i) / n; |
---|
1319 | t.i = (a.i * b.r - a.r * b.i) / n; |
---|
1320 | return t; |
---|
1321 | } |
---|
1322 | |
---|
1323 | complex Cconjugate(complex a) |
---|
1324 | { |
---|
1325 | complex t; |
---|
1326 | t.r = a.r; |
---|
1327 | t.i = - a.i; |
---|
1328 | return t; |
---|
1329 | } |
---|
1330 | |
---|
1331 | |
---|
1332 | |
---|
1333 | vector Vcreate(double xr, double xi, double yr, double yi, double zr, double zi) |
---|
1334 | { |
---|
1335 | vector temp; |
---|
1336 | |
---|
1337 | temp.x = Ccreate(xr, xi); |
---|
1338 | temp.y = Ccreate(yr, yi); |
---|
1339 | temp.z = Ccreate(zr, zi); |
---|
1340 | |
---|
1341 | return temp; |
---|
1342 | } |
---|
1343 | |
---|
1344 | vector Vadd(vector v1, vector v2) |
---|
1345 | { |
---|
1346 | vector temp; |
---|
1347 | |
---|
1348 | temp.x = Cadd(v1.x, v2.x); |
---|
1349 | temp.y = Cadd(v1.y, v2.y); |
---|
1350 | temp.z = Cadd(v1.z, v2.z); |
---|
1351 | |
---|
1352 | return temp; |
---|
1353 | } |
---|
1354 | |
---|
1355 | vector Vmultiply(vector v, complex c) |
---|
1356 | { |
---|
1357 | vector temp; |
---|
1358 | |
---|
1359 | temp.x = Cmultiply(v.x, c); |
---|
1360 | temp.y = Cmultiply(v.y, c); |
---|
1361 | temp.z = Cmultiply(v.z, c); |
---|
1362 | |
---|
1363 | return temp; |
---|
1364 | } |
---|
1365 | |
---|
1366 | vector Vconjugate(vector v) |
---|
1367 | { |
---|
1368 | vector temp; |
---|
1369 | |
---|
1370 | temp.x = Cconjugate(v.x); |
---|
1371 | temp.y = Cconjugate(v.y); |
---|
1372 | temp.z = Cconjugate(v.z); |
---|
1373 | |
---|
1374 | return temp; |
---|
1375 | } |
---|
1376 | |
---|
1377 | complex Vdot(vector v1, vector v2) |
---|
1378 | { |
---|
1379 | complex temp; |
---|
1380 | |
---|
1381 | temp = Cadd(Cmultiply(v1.x, v2.x), Cmultiply(v1.y, v2.y)); |
---|
1382 | temp = Cadd(temp, Cmultiply(v1.z, v2.z)); |
---|
1383 | |
---|
1384 | return temp; |
---|
1385 | } |
---|
1386 | |
---|
1387 | |
---|
1388 | static double QU1(Integer ndim, double var[]) /*integrand of real part of Gx, Gz*/ |
---|
1389 | { |
---|
1390 | extern double theta, xv, k, kx, alpha, beta, gam, D, dum, psi; |
---|
1391 | double u, tem1, arg, y, z; |
---|
1392 | static NagError fail; |
---|
1393 | |
---|
1394 | y=var[0]; z=var[1]; |
---|
1395 | u= xv -z*tan(alpha)+ dum; |
---|
1396 | arg= k*sqrt(u*u+y*y)/(beta*gam); |
---|
1397 | tem1= s18adc(arg,&fail); |
---|
1398 | return (2*k*u*cos(D*z-kx*dum+psi)*tem1*cos(ky*y)/(sqrt(u*u+y*y)*beta*gam)); |
---|
1399 | } |
---|
1400 | |
---|
1401 | static double QU2(Integer ndim, double var[]) |
---|
1402 | { |
---|
1403 | extern double theta, xv, k, kx, alpha, beta, gam, D, dum, psi; |
---|
1404 | double u, tem1, arg, y, z; |
---|
1405 | static NagError fail; |
---|
1406 | |
---|
1407 | y=var[0]; z=var[1]; |
---|
1408 | u= xv -z*tan(alpha)+ dum; |
---|
1409 | arg= k*sqrt(u*u+y*y)/(beta*gam); |
---|
1410 | tem1= s18adc(arg,&fail); |
---|
1411 | return (2*k*u*sin(D*z-kx*dum+psi)*tem1*cos(ky*y)/(sqrt(u*u+y*y)*beta*gam)); |
---|
1412 | } |
---|
1413 | |
---|
1414 | static double QU3(Integer ndim, double var[]) |
---|
1415 | { |
---|
1416 | extern double theta, xv, k, kx, alpha, beta, gam, D, dum, psi; |
---|
1417 | double u, tem1, arg, y, z; |
---|
1418 | static NagError fail; |
---|
1419 | |
---|
1420 | y=var[0]; z=var[1]; |
---|
1421 | u= xv -z*tan(alpha)+ dum; |
---|
1422 | arg= k*sqrt(u*u+y*y)/(beta*gam); |
---|
1423 | tem1= s18adc(arg,&fail); |
---|
1424 | return (2*k*y*sin(D*z-kx*dum+psi)*tem1*sin(ky*y)/(sqrt(u*u+y*y)*beta*gam)); |
---|
1425 | } |
---|
1426 | |
---|
1427 | static double QU4(Integer ndim, double var[]) |
---|
1428 | { |
---|
1429 | extern double theta, xv, k, kx, alpha, beta, gam, D, dum, psi; |
---|
1430 | double u, tem1, arg, y, z; |
---|
1431 | static NagError fail; |
---|
1432 | |
---|
1433 | y=var[0]; z=var[1]; |
---|
1434 | u= xv -z*tan(alpha)+ dum; |
---|
1435 | arg= k*sqrt(u*u+y*y)/(beta*gam); |
---|
1436 | tem1= s18adc(arg,&fail); |
---|
1437 | return (2*k*y*cos(D*z-kx*dum+psi)*tem1*sin(ky*y)/(sqrt(u*u+y*y)*beta*gam)); |
---|
1438 | } |
---|
1439 | |
---|
1440 | static double GAU(double x) |
---|
1441 | { |
---|
1442 | extern double x0,sigma_x; |
---|
1443 | return (exp(-(x-x0)*(x-x0)/(2*sigma_x*sigma_x))); |
---|
1444 | } |
---|
1445 | |
---|
1446 | static double FUN(double x) |
---|
1447 | { |
---|
1448 | double kk; |
---|
1449 | extern double x0,sigma_x,lam_e; |
---|
1450 | kk=2/lam_e; |
---|
1451 | return (exp(-(x-x0)*(x-x0)/(2*sigma_x*sigma_x))*exp(-kk*x)); |
---|
1452 | } |
---|
1453 | |
---|
1454 | static double FUN_X(double x) |
---|
1455 | { |
---|
1456 | double kk; |
---|
1457 | extern double x0,sigma_x,lam_e; |
---|
1458 | kk=1/lam_e; |
---|
1459 | return (exp(-(x-x0)*(x-x0)/(2*sigma_x*sigma_x))*exp(-kk*x)); |
---|
1460 | } |
---|
1461 | |
---|
1462 | static double LOR1(double t) |
---|
1463 | { |
---|
1464 | extern double a; |
---|
1465 | return (1/(a*a+t*t)); |
---|
1466 | } |
---|
1467 | |
---|
1468 | static double LOR2(double t) |
---|
1469 | { |
---|
1470 | extern double b; |
---|
1471 | return (1/(b*b+t*t)); |
---|
1472 | } |
---|
1473 | |
---|
1474 | |
---|
1475 | |
---|
1476 | |
---|