1 | #include "madx.h" |
---|
2 | |
---|
3 | // types and constants |
---|
4 | |
---|
5 | #define MIN_DOUBLE 1.e-36 |
---|
6 | |
---|
7 | struct table; |
---|
8 | struct aper_node /* aperture limit node */ |
---|
9 | { |
---|
10 | char name[NAME_L]; |
---|
11 | double n1; |
---|
12 | double s; |
---|
13 | char apertype[NAME_L]; |
---|
14 | double aperture[4]; |
---|
15 | double aper_tol[3]; |
---|
16 | double deltap_twiss; |
---|
17 | }; |
---|
18 | |
---|
19 | struct aper_e_d /* element displacement */ |
---|
20 | { |
---|
21 | char name[NAME_L]; /* element name */ |
---|
22 | int curr; /* # of rows */ |
---|
23 | double tab[E_D_MAX][3]; /* the table of read values */ |
---|
24 | }; |
---|
25 | |
---|
26 | // static interface |
---|
27 | |
---|
28 | static struct aper_e_d* true_tab; |
---|
29 | |
---|
30 | /* struct aper_e_d* offs_tab;*/ |
---|
31 | static struct table* offs_tab; |
---|
32 | |
---|
33 | static struct aper_node* |
---|
34 | aperture(char *table, struct node* use_range[], struct table* tw_cp, int *tw_cnt, struct aper_node*); |
---|
35 | |
---|
36 | /* next function replaced 19 june 2007 BJ */ |
---|
37 | /* Improved for potential zero divide 20feb08 BJ */ |
---|
38 | |
---|
39 | static int |
---|
40 | aper_rectellipse(double* ap1, double* ap2, double* ap3, double* ap4, int* quarterlength, double tablex[], double tabley[]) |
---|
41 | { |
---|
42 | double x, y, angle, alfa, theta, dangle; |
---|
43 | int i,napex; |
---|
44 | |
---|
45 | /* printf("++ %10.5f %10.5f %10.5f %10.5f\n",*ap1,*ap2,*ap3,*ap4);*/ |
---|
46 | |
---|
47 | if ( *ap1 < MIN_DOUBLE*1.e10 || *ap2 < MIN_DOUBLE*1.e10) { |
---|
48 | fatal_error("Illegal zero or too small value value for ap1 or ap2"," "); |
---|
49 | } |
---|
50 | if ( *ap3 < MIN_DOUBLE*1.e10 || *ap4 < MIN_DOUBLE*1.e10) { |
---|
51 | fatal_error("Illegal zero or too small value value for ap3 or ap4"," "); |
---|
52 | } |
---|
53 | |
---|
54 | |
---|
55 | /* Produces a table of only the first quadrant coordinates */ |
---|
56 | /* aper_fill_quads() completes the polygon */ |
---|
57 | |
---|
58 | if (*quarterlength) napex=9; |
---|
59 | else napex=19; |
---|
60 | |
---|
61 | /* PIECE OF USELESS CODE NOW THAT CASE RECTANGLE IS CORRECT, BJ 28feb08 |
---|
62 | special case : rectangle |
---|
63 | does not work -- to be reworked |
---|
64 | |
---|
65 | if ( *ap3 < 0.) { |
---|
66 | tablex[0]=(*ap1); |
---|
67 | tabley[0]=(*ap2-0.0001); |
---|
68 | tablex[1]=(*ap1-0.00002); |
---|
69 | tabley[1]=(*ap2-0.00002); |
---|
70 | tablex[2]=(*ap1-0.0001); |
---|
71 | tabley[2]=(*ap2); |
---|
72 | *quarterlength=2; |
---|
73 | return 0; |
---|
74 | } |
---|
75 | END OF PIECE OF USELESS CODE, to be removed after some use |
---|
76 | (I expect no negative feedback ...) |
---|
77 | */ |
---|
78 | |
---|
79 | /* find angles where rectangle and circle crosses */ |
---|
80 | |
---|
81 | if ( (*ap1) >= (*ap3) ) { |
---|
82 | alfa = 0.; |
---|
83 | } |
---|
84 | else { |
---|
85 | y=sqrt((*ap3)*(*ap3)-(*ap1)*(*ap1)); |
---|
86 | alfa=atan2(y,*ap1); |
---|
87 | } |
---|
88 | |
---|
89 | if ( (*ap2) >= (*ap4) ) { |
---|
90 | theta = 0.; |
---|
91 | } |
---|
92 | else { |
---|
93 | x=sqrt(((*ap3)*(*ap3)) * (1 - ((*ap2)*(*ap2)) / ((*ap4)*(*ap4)))); |
---|
94 | y=sqrt((*ap3)*(*ap3)-x*x); |
---|
95 | theta=atan2(x,y); |
---|
96 | } |
---|
97 | |
---|
98 | dangle=(pi/2-(alfa+theta))/napex; |
---|
99 | |
---|
100 | if (!((0 < dangle) && (dangle < pi/2))) |
---|
101 | { |
---|
102 | return -1; |
---|
103 | } |
---|
104 | |
---|
105 | /*write coordinates for first quadrant*/ |
---|
106 | |
---|
107 | for ( i=0 ; i<=napex; i++ ) { |
---|
108 | angle = alfa + i*dangle; |
---|
109 | tablex[i]=(*ap3)*cos(angle); |
---|
110 | tabley[i]=(*ap4)*sin(angle); |
---|
111 | |
---|
112 | if (i >= MAXARRAY/4) |
---|
113 | { |
---|
114 | fatal_error("Memory full in aper_rectellipse", "Number of coordinates exceeds set limit"); } |
---|
115 | } |
---|
116 | |
---|
117 | *quarterlength=i-1; |
---|
118 | |
---|
119 | return 0; |
---|
120 | } |
---|
121 | |
---|
122 | static void |
---|
123 | aper_adj_quad(double angle, double x, double y, double* xquad, double* yquad) |
---|
124 | { |
---|
125 | int quadrant; |
---|
126 | quadrant=angle/(pi/2)+1; |
---|
127 | switch (quadrant) |
---|
128 | { |
---|
129 | case 1: *xquad=x; *yquad=y; break; |
---|
130 | case 2: *xquad=-x; *yquad=y; break; |
---|
131 | case 3: *xquad=-x; *yquad=-y; break; |
---|
132 | case 4: *xquad=x; *yquad=-y; break; |
---|
133 | } |
---|
134 | } |
---|
135 | |
---|
136 | static void |
---|
137 | aper_adj_halo_si(double ex, double ey, double betx, double bety, double bbeat, double halox[], double haloy[], int halolength, double haloxsi[], double haloysi[]) |
---|
138 | { |
---|
139 | int j; |
---|
140 | |
---|
141 | for (j=0;j<=halolength+1;j++) |
---|
142 | { |
---|
143 | haloxsi[j]=halox[j]*bbeat*sqrt(ex*betx); |
---|
144 | haloysi[j]=haloy[j]*bbeat*sqrt(ey*bety); |
---|
145 | } |
---|
146 | } |
---|
147 | |
---|
148 | static double |
---|
149 | aper_online(double xm, double ym, double startx, double starty, |
---|
150 | double endx, double endy, double dist_limit) |
---|
151 | |
---|
152 | /* BJ 13.02.2009. |
---|
153 | - added check |x-e| < dist_limit |
---|
154 | - removed useless calculations of sqrt |
---|
155 | - made consistent use of dist_limit and min_double */ |
---|
156 | /* BJ 13.02.2009 |
---|
157 | check if point x = (xm,ym) is in the segment [s,e] with |
---|
158 | s = (startx,starty) and e = (endx,endy) by computing |
---|
159 | cosfi = (x-s).(x-e) / |x-s||x-e|. cosfi = -1 : x is in |
---|
160 | first check if |x-s| and |x-e| are not too small.If yes for one of them : in |
---|
161 | if OK , the zero divide check must be superfluous. But keep it anyway. |
---|
162 | */ |
---|
163 | |
---|
164 | { |
---|
165 | double cosfi=1 , aaa, sss, eee; |
---|
166 | |
---|
167 | sss = sqrt((xm-startx)*(xm-startx)+ (ym-starty)*(ym-starty)); |
---|
168 | eee = sqrt((xm-endx)*(xm-endx) + (ym-endy)*(ym-endy)); |
---|
169 | |
---|
170 | if ( sss <= dist_limit) |
---|
171 | { |
---|
172 | cosfi=-1; |
---|
173 | } |
---|
174 | else if ( eee <= dist_limit) |
---|
175 | { |
---|
176 | cosfi=-1; |
---|
177 | } |
---|
178 | else |
---|
179 | { |
---|
180 | aaa = sss * eee ; |
---|
181 | if ( aaa < MIN_DOUBLE ) |
---|
182 | { |
---|
183 | fatal_error("Attempt to zero divide ", "In aper_online"); |
---|
184 | } |
---|
185 | |
---|
186 | cosfi= ((xm-startx)*(xm-endx)+(ym-starty)*(ym-endy)) / aaa; |
---|
187 | } |
---|
188 | |
---|
189 | if (cosfi <= -1+dist_limit) |
---|
190 | { |
---|
191 | cosfi=-1; |
---|
192 | } |
---|
193 | return cosfi; |
---|
194 | } |
---|
195 | |
---|
196 | /* NEW VERSION of aper_race, 20feb08 BJ, potential zero-divide issues cleared */ |
---|
197 | |
---|
198 | static void |
---|
199 | aper_race(double xshift, double yshift, double r, double angle, double* x, double* y) |
---|
200 | { |
---|
201 | double angle0, angle1, angle2, alfa, gamma, theta; |
---|
202 | int quadrant; |
---|
203 | |
---|
204 | /* this function calculates the displacement of the beam centre |
---|
205 | due to tolerance uncertainty for every angle */ |
---|
206 | |
---|
207 | quadrant=angle/(pi/2)+1; |
---|
208 | |
---|
209 | if (xshift==0 && yshift==0 && r==0) |
---|
210 | { |
---|
211 | *x=0; *y=0; |
---|
212 | return; |
---|
213 | } |
---|
214 | |
---|
215 | switch (quadrant) /*adjusting angle to first quadrant*/ |
---|
216 | { |
---|
217 | case 1: angle=angle; break; |
---|
218 | case 2: angle=pi-angle; break; |
---|
219 | case 3: angle=angle-pi; break; |
---|
220 | case 4: angle=twopi-angle; break; |
---|
221 | } |
---|
222 | |
---|
223 | if (angle==pi/2) |
---|
224 | { |
---|
225 | *x=0; |
---|
226 | *y=yshift+r; |
---|
227 | } |
---|
228 | else |
---|
229 | { |
---|
230 | /*finding where arc starts and ends*/ |
---|
231 | angle0=atan2( yshift , xshift+r ); |
---|
232 | angle1=atan2( r+yshift , xshift ); |
---|
233 | |
---|
234 | /*different methods is needed, depending on angle*/ |
---|
235 | if (angle <= angle0 + MIN_DOUBLE * 1.e10 ) |
---|
236 | { |
---|
237 | *x=xshift+r; |
---|
238 | *y=tan(angle)*(xshift+r); |
---|
239 | } |
---|
240 | else if (angle<angle1) |
---|
241 | { |
---|
242 | /* if this is a circle, angle2 useless */ |
---|
243 | if (!xshift && !yshift) angle2 = 0; |
---|
244 | else angle2 = atan2( yshift , xshift ); |
---|
245 | |
---|
246 | alfa = fabs(angle-angle2); |
---|
247 | if (alfa < MIN_DOUBLE * 1.e10) |
---|
248 | { |
---|
249 | /* alfa==0 is a simpler case */ |
---|
250 | *x=cos(angle)*(r+sqrt(xshift*xshift+yshift*yshift)); |
---|
251 | *y=sin(angle)*(r+sqrt(xshift*xshift+yshift*yshift)); |
---|
252 | } |
---|
253 | else |
---|
254 | { |
---|
255 | /* solving sine rule w.r.t. gamma */ |
---|
256 | gamma=asin(sqrt(xshift*xshift+yshift*yshift)/r*sin(alfa)); |
---|
257 | /*theta is the last corner in the triangle*/ |
---|
258 | theta=pi-(alfa+gamma); |
---|
259 | *x=cos(angle)*r*sin(theta)/sin(alfa); |
---|
260 | *y=sin(angle)*r*sin(theta)/sin(alfa); |
---|
261 | } |
---|
262 | } |
---|
263 | /* upper flat part */ |
---|
264 | else |
---|
265 | { |
---|
266 | *y=r+yshift; |
---|
267 | *x=(r+yshift)*tan(pi/2-angle); |
---|
268 | } |
---|
269 | } |
---|
270 | } |
---|
271 | |
---|
272 | static int |
---|
273 | aper_chk_inside(double p, double q, double pipex[], double pipey[], double dist_limit, int pipelength) |
---|
274 | { |
---|
275 | int i; |
---|
276 | double n12, salfa, calfa, alfa=0; |
---|
277 | |
---|
278 | /*checks first whether p,q is exact on a pipe coordinate*/ |
---|
279 | for (i=0;i<=pipelength;i++) |
---|
280 | { |
---|
281 | if (-1 == aper_online(p, q, pipex[i], pipey[i], pipex[i+1], pipey[i+1], dist_limit)) |
---|
282 | { |
---|
283 | return 0; |
---|
284 | } |
---|
285 | } |
---|
286 | |
---|
287 | /*calculates and adds up angle from centre between all coordinates*/ |
---|
288 | for (i=0;i<=pipelength;i++) |
---|
289 | { |
---|
290 | n12=sqrt(((pipex[i]-p)*(pipex[i]-p) + (pipey[i]-q)*(pipey[i]-q)) |
---|
291 | * ((pipex[i+1]-p)*(pipex[i+1]-p) + (pipey[i+1]-q)*(pipey[i+1]-q))); |
---|
292 | |
---|
293 | salfa=((pipex[i]-p)*(pipey[i+1]-q) - (pipey[i]-q)*(pipex[i+1]-p))/n12; |
---|
294 | |
---|
295 | calfa=((pipex[i]-p)*(pipex[i+1]-p) + (pipey[i]-q)*(pipey[i+1]-q))/n12; |
---|
296 | |
---|
297 | alfa += atan2(salfa, calfa); |
---|
298 | } |
---|
299 | |
---|
300 | /*returns yes to main if total angle is at least twopi*/ |
---|
301 | if (sqrt(alfa*alfa)>=(twopi-dist_limit)) |
---|
302 | { |
---|
303 | return 1; |
---|
304 | } |
---|
305 | |
---|
306 | return 0; |
---|
307 | } |
---|
308 | |
---|
309 | static void |
---|
310 | aper_intersect(double a1, double b1, double a2, double b2, double x1, double y1, double x2, double y2, |
---|
311 | int ver1, int ver2, double *xm, double *ym) |
---|
312 | { |
---|
313 | (void)y1; |
---|
314 | |
---|
315 | if (ver1 && ver2 && x1==x2) |
---|
316 | { |
---|
317 | *xm=x2; |
---|
318 | *ym=y2; |
---|
319 | } |
---|
320 | else if (ver1) |
---|
321 | { |
---|
322 | *xm=x1; |
---|
323 | *ym=a2*x1+b2; |
---|
324 | } |
---|
325 | else if (ver2) |
---|
326 | { |
---|
327 | *xm=x2; |
---|
328 | *ym=a1*x2+b1; |
---|
329 | } |
---|
330 | else |
---|
331 | { |
---|
332 | *xm=(b1-b2)/(a2-a1); |
---|
333 | *ym=a1*(*xm)+b1; |
---|
334 | } |
---|
335 | } |
---|
336 | |
---|
337 | static int |
---|
338 | aper_linepar(double x1,double y1,double x2,double y2,double *a,double *b) |
---|
339 | { |
---|
340 | int vertical=0; |
---|
341 | |
---|
342 | if ( fabs(x1-x2)< MIN_DOUBLE) |
---|
343 | { |
---|
344 | *a=0; |
---|
345 | *b=y1-(*a)*x1; |
---|
346 | vertical=1; |
---|
347 | } else { |
---|
348 | *a=(y1-y2)/(x1-x2); |
---|
349 | *b=y1-(*a)*x1; |
---|
350 | } |
---|
351 | |
---|
352 | return vertical; |
---|
353 | } |
---|
354 | |
---|
355 | static void |
---|
356 | aper_fill_quads(double polyx[], double polyy[], int quarterlength, int* halolength) |
---|
357 | { |
---|
358 | int i=quarterlength+1, j; |
---|
359 | |
---|
360 | /* receives two tables with coordinates for the first quadrant */ |
---|
361 | /* and mirrors them across x and y axes */ |
---|
362 | |
---|
363 | /*copying first quadrant coordinates to second quadrant*/ |
---|
364 | for (j=quarterlength;j>=0;j--) |
---|
365 | { |
---|
366 | polyx[i]=polyx[j]; |
---|
367 | polyy[i]=polyy[j]; |
---|
368 | aper_adj_quad(pi/2, polyx[i], polyy[i], &polyx[i], &polyy[i]); |
---|
369 | i++; |
---|
370 | } |
---|
371 | |
---|
372 | /*copying first quadrant coordinates to third quadrant*/ |
---|
373 | for (j=0;j<=quarterlength;j++) |
---|
374 | { |
---|
375 | polyx[i]=polyx[j]; |
---|
376 | polyy[i]=polyy[j]; |
---|
377 | aper_adj_quad(pi, polyx[i], polyy[i], &polyx[i], &polyy[i]); |
---|
378 | i++; |
---|
379 | } |
---|
380 | |
---|
381 | /*copying first quadrant coordinates to fourth quadrant*/ |
---|
382 | for (j=quarterlength;j>=0;j--) |
---|
383 | { |
---|
384 | polyx[i]=polyx[j]; |
---|
385 | polyy[i]=polyy[j]; |
---|
386 | aper_adj_quad(pi*3/2, polyx[i], polyy[i], &polyx[i], &polyy[i]); |
---|
387 | i++; |
---|
388 | } |
---|
389 | |
---|
390 | /*sets the last point equal to the first, to complete the shape. |
---|
391 | Necessary for compatibility with aper_calc function*/ |
---|
392 | polyx[i]=polyx[0]; |
---|
393 | polyy[i]=polyy[0]; |
---|
394 | |
---|
395 | *halolength=i-1; |
---|
396 | } |
---|
397 | |
---|
398 | static void |
---|
399 | aper_read_twiss(char* table, int* jslice, double* s, double* x, double* y, |
---|
400 | double* betx, double* bety, double* dx, double* dy) |
---|
401 | { |
---|
402 | double_from_table_row(table, "s", jslice, s); |
---|
403 | double_from_table_row(table, "x", jslice, x); |
---|
404 | double_from_table_row(table, "y", jslice, y); |
---|
405 | double_from_table_row(table, "betx", jslice, betx); |
---|
406 | double_from_table_row(table, "bety", jslice, bety); |
---|
407 | double_from_table_row(table, "dx", jslice, dx); |
---|
408 | double_from_table_row(table, "dy", jslice, dy); |
---|
409 | } |
---|
410 | |
---|
411 | static int |
---|
412 | aper_external_file(char *file, double tablex[], double tabley[]) |
---|
413 | { |
---|
414 | /* receives the name of file containing coordinates. Puts coordinates into tables. */ |
---|
415 | int i=0; |
---|
416 | FILE *filept; |
---|
417 | |
---|
418 | if (file != NULL) |
---|
419 | { |
---|
420 | if ((filept=fopen(file, "r")) == NULL) |
---|
421 | { |
---|
422 | warning("Can not find file: ", file); |
---|
423 | return -1; |
---|
424 | } |
---|
425 | |
---|
426 | /*start making table*/ |
---|
427 | while (2==fscanf(filept, "%lf %lf", &tablex[i], &tabley[i])) |
---|
428 | { |
---|
429 | i++; |
---|
430 | if (i >= MAXARRAY) |
---|
431 | { |
---|
432 | fatal_error("Memory full. ", "Number of coordinates exceeds set limit"); |
---|
433 | } |
---|
434 | } |
---|
435 | |
---|
436 | tablex[i]=tablex[0]; |
---|
437 | tabley[i]=tabley[0]; |
---|
438 | fclose(filept); |
---|
439 | } |
---|
440 | return i-1; |
---|
441 | } |
---|
442 | |
---|
443 | static int |
---|
444 | aper_bs(char* apertype, double* ap1, double* ap2, double* ap3, double* ap4, int* pipelength, double pipex[], double pipey[]) |
---|
445 | { |
---|
446 | int i, err, quarterlength=0; |
---|
447 | |
---|
448 | /* "var1 .. 4" represents values in the aperture array of each element */ |
---|
449 | /* After they are read: */ |
---|
450 | /* *ap1 = half width rectangle */ |
---|
451 | /* *ap2 = half height rectangle */ |
---|
452 | /* *ap3 = half horizontal axis ellipse */ |
---|
453 | /* *ap4 = half vertical axis ellipse */ |
---|
454 | /* returns 1 on success, 0 on failure */ |
---|
455 | |
---|
456 | (*ap1)=(*ap2)=(*ap3)=(*ap4)=0; |
---|
457 | |
---|
458 | /* printf("-- #%s#\n",apertype); */ |
---|
459 | |
---|
460 | if (!strcmp(apertype,"circle")) |
---|
461 | { |
---|
462 | *ap3=get_aperture(current_node, "var1"); /*radius circle*/ |
---|
463 | |
---|
464 | *ap1 = *ap2 = *ap4 = *ap3; |
---|
465 | |
---|
466 | if (*ap3) /* check if r = 0, skip calc if r = 0 */ |
---|
467 | { |
---|
468 | err=aper_rectellipse(ap1, ap2, ap3, ap4, &quarterlength, pipex, pipey); |
---|
469 | if (!err) aper_fill_quads(pipex, pipey, quarterlength, pipelength); |
---|
470 | } |
---|
471 | else err = -1; |
---|
472 | } |
---|
473 | |
---|
474 | else if (!strcmp(apertype,"ellipse")) |
---|
475 | { |
---|
476 | *ap3 = get_aperture(current_node, "var1"); /*half hor axis ellipse*/ |
---|
477 | *ap4 = get_aperture(current_node, "var2"); /*half ver axis ellipse*/ |
---|
478 | |
---|
479 | *ap1 = *ap3; *ap2 = *ap4; |
---|
480 | |
---|
481 | err=aper_rectellipse(ap1, ap2, ap3, ap4, &quarterlength, pipex, pipey); |
---|
482 | if (!err) aper_fill_quads(pipex, pipey, quarterlength, pipelength); |
---|
483 | } |
---|
484 | |
---|
485 | else if (!strcmp(apertype,"rectangle")) |
---|
486 | { |
---|
487 | |
---|
488 | *ap1 = get_aperture(current_node, "var1"); /*half width rect*/ |
---|
489 | *ap2 = get_aperture(current_node, "var2"); /*half height rect*/ |
---|
490 | /* next changed 28 feb 2008 BJ */ |
---|
491 | *ap3 = *ap4 = sqrt((*ap1) * (*ap1) + ((*ap2) * (*ap2))) - 1.e-15; |
---|
492 | /* printf("-- %10.5f %10.5f %10.5f %10.5f\n",*ap1,*ap2,*ap3,*ap4); */ |
---|
493 | |
---|
494 | err=aper_rectellipse(ap1, ap2, ap3, ap4, &quarterlength, pipex, pipey); |
---|
495 | if (!err) aper_fill_quads(pipex, pipey, quarterlength, pipelength); |
---|
496 | } |
---|
497 | |
---|
498 | else if (!strcmp(apertype,"lhcscreen")) |
---|
499 | { |
---|
500 | *ap1=get_aperture(current_node, "var1"); /*half width rect*/ |
---|
501 | *ap2=get_aperture(current_node, "var2"); /*half height rect*/ |
---|
502 | *ap3=get_aperture(current_node, "var3"); /*radius circle*/ |
---|
503 | |
---|
504 | (*ap4) = (*ap3); |
---|
505 | |
---|
506 | err=aper_rectellipse(ap1, ap2, ap3, ap4, &quarterlength, pipex, pipey); |
---|
507 | if (!err) aper_fill_quads(pipex, pipey, quarterlength, pipelength); |
---|
508 | } |
---|
509 | |
---|
510 | else if (!strcmp(apertype,"marguerite")) |
---|
511 | { |
---|
512 | printf("\nApertype %s not yet supported.", apertype); |
---|
513 | err=-1; |
---|
514 | } |
---|
515 | |
---|
516 | else if (!strcmp(apertype,"rectellipse")) |
---|
517 | { |
---|
518 | *ap1=get_aperture(current_node, "var1"); /*half width rect*/ |
---|
519 | *ap2=get_aperture(current_node, "var2"); /*half height rect*/ |
---|
520 | *ap3=get_aperture(current_node, "var3"); /*half hor axis ellipse*/ |
---|
521 | *ap4=get_aperture(current_node, "var4"); /*half ver axis ellipse*/ |
---|
522 | |
---|
523 | if (*ap1==0) /*this will not be 0 in the future*/ |
---|
524 | { |
---|
525 | *ap1=*ap3; |
---|
526 | } |
---|
527 | if (*ap2==0) /*this will not be 0 in the future*/ |
---|
528 | { |
---|
529 | *ap2=*ap4; |
---|
530 | } |
---|
531 | |
---|
532 | err=aper_rectellipse(ap1, ap2, ap3, ap4, &quarterlength, pipex, pipey); |
---|
533 | if (!err) aper_fill_quads(pipex, pipey, quarterlength, pipelength); |
---|
534 | } |
---|
535 | |
---|
536 | else if (!strcmp(apertype,"racetrack")) |
---|
537 | { |
---|
538 | *ap1=get_aperture(current_node, "var1"); /*half width rect*/ |
---|
539 | *ap2=get_aperture(current_node, "var2"); /*half height rect*/ |
---|
540 | *ap3=get_aperture(current_node, "var3"); /*radius circle*/ |
---|
541 | |
---|
542 | *ap4 = *ap3; |
---|
543 | |
---|
544 | err=aper_rectellipse(ap3, ap3, ap3, ap4, &quarterlength, pipex, pipey); |
---|
545 | |
---|
546 | if (!err) |
---|
547 | { |
---|
548 | /*displaces the quartercircle*/ |
---|
549 | for (i=0;i<=quarterlength;i++) |
---|
550 | { |
---|
551 | pipex[i] += (*ap1); |
---|
552 | pipey[i] += (*ap2); |
---|
553 | } |
---|
554 | |
---|
555 | aper_fill_quads(pipex, pipey, quarterlength, pipelength); |
---|
556 | } |
---|
557 | } |
---|
558 | |
---|
559 | else if (strlen(apertype)) |
---|
560 | { |
---|
561 | *pipelength = aper_external_file(apertype, pipex, pipey); |
---|
562 | *ap1 = *ap2 = *ap3 = *ap4 = 0; |
---|
563 | if (*pipelength > -1) err=0; else err=-1; |
---|
564 | } |
---|
565 | |
---|
566 | else |
---|
567 | { |
---|
568 | *pipelength = -1; |
---|
569 | err=-1; |
---|
570 | } |
---|
571 | |
---|
572 | return err+1; |
---|
573 | } |
---|
574 | |
---|
575 | static int |
---|
576 | aper_tab_search(int cnt, struct aper_e_d* tab, char* name, int* pos) |
---|
577 | { |
---|
578 | /* looks for node *name in tab[], returns 1 if found, and its pos */ |
---|
579 | int i=-1, found=0; |
---|
580 | |
---|
581 | while (i < cnt && found == 0) |
---|
582 | { |
---|
583 | i++; |
---|
584 | if (strcmp(name,tab[i].name) == 0) found=1; |
---|
585 | } |
---|
586 | *pos=i; |
---|
587 | |
---|
588 | return found; |
---|
589 | } |
---|
590 | |
---|
591 | static int |
---|
592 | aper_tab_search_tfs(struct table* t, char* name, double* row) |
---|
593 | { |
---|
594 | /* looks for node *name in tab[], returns 1 if found, and its pos |
---|
595 | NAME S_IP X_OFF DX_OFF DDX_OFF Y_OFF DY_OFF DDY_OFF |
---|
596 | function value return: |
---|
597 | 1 ok |
---|
598 | 0 not found |
---|
599 | -1 column does not exist |
---|
600 | */ |
---|
601 | |
---|
602 | int i=0, found=0; |
---|
603 | int name_pos, s_ip_pos, x_off_pos, dx_off_pos, ddx_off_pos, y_off_pos, dy_off_pos, ddy_off_pos; |
---|
604 | |
---|
605 | /* get column positions */ |
---|
606 | mycpy(c_dum->c, "name"); |
---|
607 | if ((name_pos = name_list_pos(c_dum->c, t->columns)) < 0) return -1; |
---|
608 | mycpy(c_dum->c, "s_ip"); |
---|
609 | if ((s_ip_pos = name_list_pos(c_dum->c, t->columns)) < 0) return -1; |
---|
610 | mycpy(c_dum->c, "x_off"); |
---|
611 | if ((x_off_pos = name_list_pos(c_dum->c, t->columns)) < 0) return -1; |
---|
612 | mycpy(c_dum->c, "dx_off"); |
---|
613 | if ((dx_off_pos = name_list_pos(c_dum->c, t->columns)) < 0) return -1; |
---|
614 | mycpy(c_dum->c, "ddx_off"); |
---|
615 | if ((ddx_off_pos = name_list_pos(c_dum->c, t->columns)) < 0) return -1; |
---|
616 | mycpy(c_dum->c, "y_off"); |
---|
617 | if ((y_off_pos = name_list_pos(c_dum->c, t->columns)) < 0) return -1; |
---|
618 | mycpy(c_dum->c, "dy_off"); |
---|
619 | if ((dy_off_pos = name_list_pos(c_dum->c, t->columns)) < 0) return -1; |
---|
620 | mycpy(c_dum->c, "ddy_off"); |
---|
621 | if ((ddy_off_pos = name_list_pos(c_dum->c, t->columns)) < 0) return -1; |
---|
622 | |
---|
623 | /* printf("\n column names ok %d\n",t->curr);*/ |
---|
624 | |
---|
625 | while (i < t->curr && found == 0) |
---|
626 | { |
---|
627 | i++; |
---|
628 | if( !strcmp(t->s_cols[name_pos][i-1],name)) { |
---|
629 | row[1] = t->d_cols[s_ip_pos][i-1]; |
---|
630 | row[2] = t->d_cols[x_off_pos][i-1]; |
---|
631 | row[3] = t->d_cols[dx_off_pos][i-1]; |
---|
632 | row[4] = t->d_cols[ddx_off_pos][i-1]; |
---|
633 | row[5] = t->d_cols[y_off_pos][i-1]; |
---|
634 | row[6] = t->d_cols[dy_off_pos][i-1]; |
---|
635 | row[7] = t->d_cols[ddy_off_pos][i-1]; |
---|
636 | found = 1; |
---|
637 | } |
---|
638 | } |
---|
639 | |
---|
640 | return found; |
---|
641 | |
---|
642 | } |
---|
643 | |
---|
644 | static int |
---|
645 | aper_e_d_read(char* e_d_name, struct aper_e_d** e_d_tabp, int* cnt, char* refnode) |
---|
646 | { |
---|
647 | /* Reads data for special displacements of some magnets */ |
---|
648 | int i=1, j, k, l, e_d_flag=0, curr_e_d_max = E_D_LIST_CHUNK, new_e_d_max; |
---|
649 | char comment[100]="empty"; |
---|
650 | char *strpt; |
---|
651 | FILE *e_d_pt; |
---|
652 | struct aper_e_d* e_d_tab_loc; |
---|
653 | struct aper_e_d* e_d_tab = *e_d_tabp; |
---|
654 | |
---|
655 | |
---|
656 | if (e_d_name != NULL) |
---|
657 | { |
---|
658 | if((e_d_pt = fopen(e_d_name,"r")) == NULL) |
---|
659 | { |
---|
660 | printf("\nFile does not exist: %s\n",e_d_name); |
---|
661 | } |
---|
662 | else |
---|
663 | { |
---|
664 | /* part for reading reference node */ |
---|
665 | while (strncmp(comment,"reference:",10) && i != EOF) |
---|
666 | { |
---|
667 | /*fgets(buf, 100, e_d_pt);*/ |
---|
668 | i = fscanf(e_d_pt, "%s", comment); |
---|
669 | stolower(comment); |
---|
670 | } |
---|
671 | |
---|
672 | if (i == EOF) rewind(e_d_pt); |
---|
673 | else |
---|
674 | { |
---|
675 | if (strlen(comment) != 10) |
---|
676 | { |
---|
677 | strpt=strchr(comment,':'); |
---|
678 | strpt++; |
---|
679 | strcpy(refnode, strpt); |
---|
680 | } |
---|
681 | else i = fscanf(e_d_pt, "%s", refnode); |
---|
682 | |
---|
683 | stolower(refnode); |
---|
684 | strcat(refnode, ":1"); |
---|
685 | } |
---|
686 | /* printf("\nReference node: %s\n",refnode);*/ |
---|
687 | /* end reading reference node */ |
---|
688 | |
---|
689 | i=0; |
---|
690 | |
---|
691 | while (i != EOF) |
---|
692 | { |
---|
693 | i=fscanf(e_d_pt, "%s", e_d_tab[*cnt].name); |
---|
694 | /*next while-loop treats comments*/ |
---|
695 | while ( e_d_tab[*cnt].name[0] == '!' && i != EOF) |
---|
696 | { |
---|
697 | fgets(comment, 100, e_d_pt); |
---|
698 | i=fscanf(e_d_pt, "%s", e_d_tab[*cnt].name); |
---|
699 | } |
---|
700 | |
---|
701 | stolower(e_d_tab[*cnt].name); |
---|
702 | |
---|
703 | if (i != EOF) |
---|
704 | { |
---|
705 | strcat(e_d_tab[*cnt].name, ":1"); |
---|
706 | |
---|
707 | k=0; j=3; |
---|
708 | while (j == 3 && k < E_D_MAX) |
---|
709 | { |
---|
710 | j=fscanf(e_d_pt, "%lf %lf %lf", &e_d_tab[*cnt].tab[k][0], |
---|
711 | &e_d_tab[*cnt].tab[k][1], |
---|
712 | &e_d_tab[*cnt].tab[k][2]); |
---|
713 | k++; |
---|
714 | |
---|
715 | if (e_d_tab[*cnt].curr == E_D_MAX) printf("\nToo many points of x,y displacement...\n"); |
---|
716 | } |
---|
717 | |
---|
718 | e_d_tab[*cnt].curr=k-2; |
---|
719 | |
---|
720 | (*cnt)++; |
---|
721 | |
---|
722 | if (*cnt == curr_e_d_max) /* grow e_d array */ |
---|
723 | { |
---|
724 | /* printf("\nToo many special elements...(less than %d expected)\n", E_D_MAX); */ |
---|
725 | new_e_d_max = curr_e_d_max + E_D_LIST_CHUNK; |
---|
726 | printf("\ngrowin e_d_max array to %d\n", new_e_d_max); |
---|
727 | |
---|
728 | e_d_tab_loc = (struct aper_e_d*) mycalloc("Aperture",new_e_d_max,sizeof(struct aper_e_d) ); |
---|
729 | |
---|
730 | for( l=0 ; l < curr_e_d_max; l++) |
---|
731 | { |
---|
732 | e_d_tab_loc[l] = e_d_tab[l]; |
---|
733 | } |
---|
734 | |
---|
735 | |
---|
736 | myfree("Aperture",e_d_tab); |
---|
737 | |
---|
738 | e_d_tab = e_d_tab_loc; |
---|
739 | |
---|
740 | curr_e_d_max = new_e_d_max; |
---|
741 | |
---|
742 | } |
---|
743 | |
---|
744 | i=j; |
---|
745 | } |
---|
746 | } /* while !EOF */ |
---|
747 | |
---|
748 | printf("\nUsing extra displacements from file \"%s\"\n",e_d_name); |
---|
749 | e_d_flag=1; fclose(e_d_pt); |
---|
750 | (*cnt)--; |
---|
751 | } |
---|
752 | } |
---|
753 | |
---|
754 | *e_d_tabp = e_d_tab; |
---|
755 | |
---|
756 | return e_d_flag; |
---|
757 | } |
---|
758 | |
---|
759 | static struct table* |
---|
760 | aper_e_d_read_tfs(char* e_d_name, int* cnt, char* refnode) |
---|
761 | { |
---|
762 | /* Reads displacement data in tfs format */ |
---|
763 | struct table* t = NULL; |
---|
764 | struct char_p_array* tcpa = NULL; |
---|
765 | struct name_list* tnl = NULL; |
---|
766 | int i, k, error = 0; |
---|
767 | short sk; |
---|
768 | char *cc, *tmp, *name; |
---|
769 | int tempcount; |
---|
770 | tempcount = 0; |
---|
771 | |
---|
772 | (void)cnt; |
---|
773 | |
---|
774 | if (e_d_name == NULL) return NULL; |
---|
775 | |
---|
776 | printf("\n Reading offsets from tfs \"%s\"\n",e_d_name); |
---|
777 | |
---|
778 | |
---|
779 | if ((tab_file = fopen(e_d_name, "r")) == NULL) |
---|
780 | { |
---|
781 | warning("cannot open file:", e_d_name); return NULL; |
---|
782 | } |
---|
783 | |
---|
784 | while (fgets(aux_buff->c, aux_buff->max, tab_file)) |
---|
785 | { |
---|
786 | tempcount++; |
---|
787 | |
---|
788 | cc = strtok(aux_buff->c, " \"\n"); |
---|
789 | |
---|
790 | if (*cc == '@') |
---|
791 | { |
---|
792 | if ((tmp = strtok(NULL, " \"\n")) != NULL |
---|
793 | && strcmp(tmp, "REFERENCE") == 0) /* search for reference node */ |
---|
794 | { |
---|
795 | if ((name = strtok(NULL, " \"\n")) != NULL) /* skip format */ |
---|
796 | { |
---|
797 | if ((name = strtok(NULL, " \"\n")) != NULL) { |
---|
798 | strcpy(refnode, name); |
---|
799 | stolower(refnode); |
---|
800 | strcat(refnode, ":1"); |
---|
801 | } |
---|
802 | } |
---|
803 | } |
---|
804 | /* printf("\n+++++ Reference node: %s\n",refnode); */ |
---|
805 | } |
---|
806 | |
---|
807 | |
---|
808 | else if (*cc == '*' && tnl == NULL) |
---|
809 | { |
---|
810 | tnl = new_name_list("table_names", 20); |
---|
811 | while ((tmp = strtok(NULL, " \"\n")) != NULL) |
---|
812 | add_to_name_list(permbuff(stolower(tmp)), 0, tnl); |
---|
813 | } |
---|
814 | else if (*cc == '$' && tcpa == NULL) |
---|
815 | { |
---|
816 | |
---|
817 | if (tnl == NULL) |
---|
818 | { |
---|
819 | warning("formats before names","skipped"); return NULL; |
---|
820 | } |
---|
821 | tcpa = new_char_p_array(20); |
---|
822 | while ((tmp = strtok(NULL, " \"\n")) != NULL) |
---|
823 | { |
---|
824 | if (tcpa->curr == tcpa->max) grow_char_p_array(tcpa); |
---|
825 | if (strcmp(tmp, "%s") == 0) tnl->inform[tcpa->curr] = 3; |
---|
826 | else if (strcmp(tmp, "%hd") == 0) tnl->inform[tcpa->curr] = 1; |
---|
827 | else if (strcmp(tmp, "%d") == 0) tnl->inform[tcpa->curr] = 1; |
---|
828 | else tnl->inform[tcpa->curr] = 2; |
---|
829 | tcpa->p[tcpa->curr++] = permbuff(tmp); |
---|
830 | } |
---|
831 | } |
---|
832 | else |
---|
833 | { |
---|
834 | if(t == NULL) |
---|
835 | { |
---|
836 | if (tcpa == NULL) |
---|
837 | { |
---|
838 | warning("TFS table without formats,","skipped"); error = 1; |
---|
839 | } |
---|
840 | else if (tnl == NULL) |
---|
841 | { |
---|
842 | warning("TFS table without column names,","skipped"); error = 1; |
---|
843 | } |
---|
844 | else if (tnl->curr == 0) |
---|
845 | { |
---|
846 | warning("TFS table: empty column name list,","skipped"); |
---|
847 | error = 1; |
---|
848 | } |
---|
849 | else if (tnl->curr != tcpa->curr) |
---|
850 | { |
---|
851 | warning("TFS table: number of names and formats differ,", |
---|
852 | "skipped"); |
---|
853 | error = 1; |
---|
854 | } |
---|
855 | if (error) |
---|
856 | { |
---|
857 | delete_name_list(tnl); return NULL; |
---|
858 | } |
---|
859 | if(e_d_name != NULL) { |
---|
860 | t = new_table(e_d_name, "OFFSETS", 500, tnl); |
---|
861 | } else { |
---|
862 | t = new_table(e_d_name, "OFFSETS", 500, tnl); |
---|
863 | } |
---|
864 | t->curr = 0; |
---|
865 | } |
---|
866 | |
---|
867 | for (i = 0; i < tnl->curr; i++) |
---|
868 | { |
---|
869 | if (t->curr == t->max) grow_table(t); |
---|
870 | tmp = tcpa->p[i]; |
---|
871 | if (strcmp(tmp,"%s") == 0) { |
---|
872 | t->s_cols[i][t->curr] = stolower(tmpbuff(cc)); |
---|
873 | strcat(t->s_cols[i][t->curr], ":1"); |
---|
874 | } |
---|
875 | else if (strcmp(tmp,"%d") == 0 ) |
---|
876 | { |
---|
877 | sscanf(cc, tmp, &k); t->d_cols[i][t->curr] = k; |
---|
878 | } |
---|
879 | else if (strcmp(tmp,"%hd") == 0 ) |
---|
880 | { |
---|
881 | sscanf(cc, tmp, &sk); t->d_cols[i][t->curr] = sk; |
---|
882 | } |
---|
883 | else sscanf(cc, tmp, &t->d_cols[i][t->curr]); |
---|
884 | if (i+1 < tnl->curr) |
---|
885 | { |
---|
886 | if ((cc =strtok(NULL, " \"\n")) == NULL) |
---|
887 | { |
---|
888 | warning("incomplete table line starting with:", aux_buff->c); |
---|
889 | return NULL; |
---|
890 | } |
---|
891 | } |
---|
892 | } |
---|
893 | t->curr++; |
---|
894 | } |
---|
895 | } |
---|
896 | |
---|
897 | fclose(tab_file); |
---|
898 | t->origin = 1; |
---|
899 | /* next line commented : avoid memory error at 2nd APERTURE command */ |
---|
900 | /* when the offset file has the same same, BJ 8apr2009 */ |
---|
901 | /* add_to_table_list(t, table_register);*/ |
---|
902 | return t; |
---|
903 | } |
---|
904 | |
---|
905 | static void |
---|
906 | aper_header(struct table* aper_t, struct aper_node *lim) |
---|
907 | /* puts beam and aperture parameters at start of the aperture table */ |
---|
908 | { |
---|
909 | int i, h_length = 25; // not used, nint=1; |
---|
910 | double dtmp, vtmp[4]; // not used, deltap_twiss; |
---|
911 | char tmp[NAME_L], name[NAME_L], *stmp; |
---|
912 | |
---|
913 | strncpy(name, lim->name, sizeof name); |
---|
914 | |
---|
915 | /* =================================================================*/ |
---|
916 | /* ATTENTION: if you add header lines, augment h_length accordingly */ |
---|
917 | /* =================================================================*/ |
---|
918 | |
---|
919 | |
---|
920 | /* many modif to make the header being standard; BJ 25feb2008 */ |
---|
921 | |
---|
922 | if (aper_t == NULL) return; |
---|
923 | stmp = command_par_string("pipefile", this_cmd->clone); |
---|
924 | if (stmp) h_length++; |
---|
925 | stmp = command_par_string("halofile", this_cmd->clone); |
---|
926 | if (stmp) h_length += 1; else h_length += 4; |
---|
927 | |
---|
928 | // printf("\nheader %d \n",h_length); |
---|
929 | |
---|
930 | /* beam properties */ |
---|
931 | if (aper_t->header == NULL) aper_t->header = new_char_p_array(h_length); |
---|
932 | strncpy(tmp, current_sequ->name, sizeof tmp); |
---|
933 | sprintf(c_dum->c, v_format("@ SEQUENCE %%%02ds \"%s\""),strlen(tmp),stoupper(tmp)); |
---|
934 | aper_t->header->p[aper_t->header->curr++] = tmpbuff(c_dum->c); |
---|
935 | i = get_string("beam", "particle", tmp); |
---|
936 | sprintf(c_dum->c, v_format("@ PARTICLE %%%02ds \"%s\""),i,stoupper(tmp)); |
---|
937 | aper_t->header->p[aper_t->header->curr++] = tmpbuff(c_dum->c); |
---|
938 | dtmp = get_value("beam", "mass"); |
---|
939 | sprintf(c_dum->c, v_format("@ MASS %%le %F"), dtmp); |
---|
940 | aper_t->header->p[aper_t->header->curr++] = tmpbuff(c_dum->c); |
---|
941 | dtmp = get_value("beam", "energy"); |
---|
942 | sprintf(c_dum->c, v_format("@ ENERGY %%le %F"), dtmp); |
---|
943 | aper_t->header->p[aper_t->header->curr++] = tmpbuff(c_dum->c); |
---|
944 | dtmp = get_value("beam", "pc"); |
---|
945 | sprintf(c_dum->c, v_format("@ PC %%le %F"), dtmp); |
---|
946 | aper_t->header->p[aper_t->header->curr++] = tmpbuff(c_dum->c); |
---|
947 | dtmp = get_value("beam", "gamma"); |
---|
948 | sprintf(c_dum->c, v_format("@ GAMMA %%le %F"), dtmp); |
---|
949 | aper_t->header->p[aper_t->header->curr++] = tmpbuff(c_dum->c); |
---|
950 | |
---|
951 | |
---|
952 | /* aperture command properties */ |
---|
953 | |
---|
954 | dtmp = command_par_value("exn", this_cmd->clone); |
---|
955 | sprintf(c_dum->c, v_format("@ EXN %%le %F"), dtmp); |
---|
956 | aper_t->header->p[aper_t->header->curr++] = tmpbuff(c_dum->c); |
---|
957 | dtmp = command_par_value("eyn", this_cmd->clone); |
---|
958 | sprintf(c_dum->c, v_format("@ EYN %%le %F"), dtmp); |
---|
959 | aper_t->header->p[aper_t->header->curr++] = tmpbuff(c_dum->c); |
---|
960 | dtmp = command_par_value("dqf", this_cmd->clone); |
---|
961 | sprintf(c_dum->c, v_format("@ DQF %%le %F"), dtmp); |
---|
962 | aper_t->header->p[aper_t->header->curr++] = tmpbuff(c_dum->c); |
---|
963 | dtmp = command_par_value("betaqfx", this_cmd->clone); |
---|
964 | sprintf(c_dum->c, v_format("@ BETAQFX %%le %F"), dtmp); |
---|
965 | aper_t->header->p[aper_t->header->curr++] = tmpbuff(c_dum->c); |
---|
966 | dtmp = command_par_value("dparx", this_cmd->clone); |
---|
967 | sprintf(c_dum->c, v_format("@ PARAS_DX %%le %g"), dtmp); |
---|
968 | aper_t->header->p[aper_t->header->curr++] = tmpbuff(c_dum->c); |
---|
969 | dtmp = command_par_value("dpary", this_cmd->clone); |
---|
970 | sprintf(c_dum->c, v_format("@ PARAS_DY %%le %g"), dtmp); |
---|
971 | aper_t->header->p[aper_t->header->curr++] = tmpbuff(c_dum->c); |
---|
972 | dtmp = command_par_value("dp", this_cmd->clone); |
---|
973 | sprintf(c_dum->c, v_format("@ DP_BUCKET_SIZE %%le %F"), dtmp); |
---|
974 | aper_t->header->p[aper_t->header->curr++] = tmpbuff(c_dum->c); |
---|
975 | |
---|
976 | // LD: summary table is corrupted by embedded twiss, use saved value |
---|
977 | // double_from_table_row("summ","deltap",&nint,&deltap_twiss); |
---|
978 | sprintf(c_dum->c, v_format("@ TWISS_DELTAP %%le %F"), lim->deltap_twiss); |
---|
979 | aper_t->header->p[aper_t->header->curr++] = tmpbuff(c_dum->c); |
---|
980 | |
---|
981 | dtmp = command_par_value("cor", this_cmd->clone); |
---|
982 | sprintf(c_dum->c, v_format("@ CO_RADIUS %%le %F"), dtmp); |
---|
983 | aper_t->header->p[aper_t->header->curr++] = tmpbuff(c_dum->c); |
---|
984 | dtmp = command_par_value("bbeat", this_cmd->clone); |
---|
985 | sprintf(c_dum->c, v_format("@ BETA_BEATING %%le %F"), dtmp); |
---|
986 | aper_t->header->p[aper_t->header->curr++] = tmpbuff(c_dum->c); |
---|
987 | dtmp = command_par_value("nco", this_cmd->clone); |
---|
988 | sprintf(c_dum->c, v_format("@ NB_OF_ANGLES %%d %g"), dtmp*4); |
---|
989 | aper_t->header->p[aper_t->header->curr++] = tmpbuff(c_dum->c); |
---|
990 | |
---|
991 | /* if a filename with halo coordinates is given, need not show halo */ |
---|
992 | stmp = command_par_string("halofile", this_cmd->clone); |
---|
993 | if (stmp) |
---|
994 | { |
---|
995 | strncpy(tmp, stmp, sizeof tmp); |
---|
996 | sprintf(c_dum->c, v_format("@ HALOFILE %%%02ds \"%s\""),strlen(tmp),stoupper(tmp)); |
---|
997 | aper_t->header->p[aper_t->header->curr++] = tmpbuff(c_dum->c); |
---|
998 | } |
---|
999 | else |
---|
1000 | { |
---|
1001 | i = command_par_vector("halo", this_cmd->clone, vtmp); |
---|
1002 | sprintf(c_dum->c, v_format("@ HALO_PRIM %%le %g"),vtmp[0]); |
---|
1003 | aper_t->header->p[aper_t->header->curr++] = tmpbuff(c_dum->c); |
---|
1004 | sprintf(c_dum->c, v_format("@ HALO_R %%le %g"),vtmp[1]); |
---|
1005 | aper_t->header->p[aper_t->header->curr++] = tmpbuff(c_dum->c); |
---|
1006 | sprintf(c_dum->c, v_format("@ HALO_H %%le %g"),vtmp[2]); |
---|
1007 | aper_t->header->p[aper_t->header->curr++] = tmpbuff(c_dum->c); |
---|
1008 | sprintf(c_dum->c, v_format("@ HALO_V %%le %g"),vtmp[3]); |
---|
1009 | aper_t->header->p[aper_t->header->curr++] = tmpbuff(c_dum->c); |
---|
1010 | } |
---|
1011 | /* show filename with pipe coordinates if given */ |
---|
1012 | stmp = command_par_string("pipefile", this_cmd->clone); |
---|
1013 | if (stmp) |
---|
1014 | { |
---|
1015 | strncpy(tmp, stmp, sizeof tmp); |
---|
1016 | sprintf(c_dum->c, v_format("@ PIPEFILE %%%02ds \"%s\""), strlen(tmp), stoupper(tmp) ); |
---|
1017 | aper_t->header->p[aper_t->header->curr++] = tmpbuff(c_dum->c); |
---|
1018 | } |
---|
1019 | |
---|
1020 | sprintf(c_dum->c, v_format("@ n1min %%le %g"), lim->n1); |
---|
1021 | aper_t->header->p[aper_t->header->curr++] = tmpbuff(c_dum->c); |
---|
1022 | set_value("beam","n1min",&lim->n1); |
---|
1023 | |
---|
1024 | sprintf(c_dum->c, v_format("@ at_element %%%02ds \"%s\""), strlen(name), stoupper(name) ); |
---|
1025 | aper_t->header->p[aper_t->header->curr++] = tmpbuff(c_dum->c); |
---|
1026 | } |
---|
1027 | |
---|
1028 | static void |
---|
1029 | aper_surv(double init[], int nint) |
---|
1030 | { |
---|
1031 | struct in_cmd* aper_survey; |
---|
1032 | struct name_list* asnl; |
---|
1033 | int aspos; |
---|
1034 | |
---|
1035 | /* Constructs artificial survey command, the result is the */ |
---|
1036 | /* table 'survey' which can be accessed from all functions. */ |
---|
1037 | /* init[0] = x0, init[1] = y0, init[2] = z0, */ |
---|
1038 | /* init[3] = theta0, init[4] = phi0, init[5] = psi0 */ |
---|
1039 | |
---|
1040 | aper_survey = new_in_cmd(10); |
---|
1041 | aper_survey->type = 0; |
---|
1042 | aper_survey->clone = aper_survey->cmd_def = clone_command(find_command("survey",defined_commands)); |
---|
1043 | asnl = aper_survey->cmd_def->par_names; |
---|
1044 | aspos = name_list_pos("table", asnl); |
---|
1045 | aper_survey->cmd_def->par->parameters[aspos]->string = "survey"; |
---|
1046 | aper_survey->cmd_def->par_names->inform[aspos] = 1; |
---|
1047 | |
---|
1048 | aspos = name_list_pos("x0", asnl); |
---|
1049 | aper_survey->cmd_def->par->parameters[aspos]->double_value = init[0]; |
---|
1050 | aper_survey->cmd_def->par_names->inform[aspos] = 1; |
---|
1051 | |
---|
1052 | aspos = name_list_pos("y0", asnl); |
---|
1053 | aper_survey->cmd_def->par->parameters[aspos]->double_value = init[1]; |
---|
1054 | aper_survey->cmd_def->par_names->inform[aspos] = 1; |
---|
1055 | |
---|
1056 | aspos = name_list_pos("z0", asnl); |
---|
1057 | aper_survey->cmd_def->par->parameters[aspos]->double_value = init[2]; |
---|
1058 | aper_survey->cmd_def->par_names->inform[aspos] = 1; |
---|
1059 | |
---|
1060 | aspos = name_list_pos("theta0", asnl); |
---|
1061 | aper_survey->cmd_def->par->parameters[aspos]->double_value = init[3]; |
---|
1062 | aper_survey->cmd_def->par_names->inform[aspos] = 1; |
---|
1063 | |
---|
1064 | aspos = name_list_pos("phi0", asnl); |
---|
1065 | aper_survey->cmd_def->par->parameters[aspos]->double_value = init[4]; |
---|
1066 | aper_survey->cmd_def->par_names->inform[aspos] = 1; |
---|
1067 | |
---|
1068 | aspos = name_list_pos("psi0", asnl); |
---|
1069 | aper_survey->cmd_def->par->parameters[aspos]->double_value = init[5]; |
---|
1070 | aper_survey->cmd_def->par_names->inform[aspos] = 1; |
---|
1071 | |
---|
1072 | /* frs: suppressing the survey file created by the internal survey command */ |
---|
1073 | aspos = name_list_pos("file", asnl); |
---|
1074 | aper_survey->cmd_def->par->parameters[aspos]->string = NULL; |
---|
1075 | aper_survey->cmd_def->par_names->inform[aspos] = 0; |
---|
1076 | |
---|
1077 | current_survey=aper_survey->clone; |
---|
1078 | pro_survey(aper_survey); |
---|
1079 | |
---|
1080 | double_from_table_row("survey","x",&nint, &init[0]); |
---|
1081 | double_from_table_row("survey","y",&nint, &init[1]); |
---|
1082 | double_from_table_row("survey","z",&nint, &init[2]); |
---|
1083 | double_from_table_row("survey","theta",&nint, &init[3]); |
---|
1084 | double_from_table_row("survey","phi",&nint, &init[4]); |
---|
1085 | double_from_table_row("survey","psi",&nint, &init[5]); |
---|
1086 | } |
---|
1087 | |
---|
1088 | static void |
---|
1089 | aper_trim_ws(char* string, int len) |
---|
1090 | { |
---|
1091 | int c=0; |
---|
1092 | |
---|
1093 | /* Replaces the first ws or : in a string with a '\0', */ |
---|
1094 | /* thus translating a FORTRAN-like attribute string to */ |
---|
1095 | /* C compatibility, or washes the ':1' from node names */ |
---|
1096 | |
---|
1097 | while (string[c] && string[c]!=' ' && c<len-1) c++; |
---|
1098 | |
---|
1099 | string[c]=0; |
---|
1100 | if (c<len-1) string[c+1]=' '; /*adds a ws to avoid two \0 in a row*/ |
---|
1101 | } |
---|
1102 | |
---|
1103 | static void |
---|
1104 | aper_write_table(char* name, double* n1, double* n1x_m, double* n1y_m, |
---|
1105 | double* rtol, double* xtol, double* ytol, |
---|
1106 | char* apertype,double* ap1,double* ap2,double* ap3,double* ap4, |
---|
1107 | double* on_ap, double* on_elem, double* spec,double* s, |
---|
1108 | double* x, double* y, double* betx, double* bety,double* dx, double* dy, |
---|
1109 | char *table) |
---|
1110 | { |
---|
1111 | string_to_table_curr(table, "name", name); |
---|
1112 | double_to_table_curr(table, "n1", n1); |
---|
1113 | double_to_table_curr(table, "n1x_m", n1x_m); |
---|
1114 | double_to_table_curr(table, "n1y_m", n1y_m); |
---|
1115 | double_to_table_curr(table, "rtol", rtol); |
---|
1116 | double_to_table_curr(table, "xtol", xtol); |
---|
1117 | double_to_table_curr(table, "ytol", ytol); |
---|
1118 | string_to_table_curr(table, "apertype", apertype); |
---|
1119 | double_to_table_curr(table, "aper_1", ap1); |
---|
1120 | double_to_table_curr(table, "aper_2", ap2); |
---|
1121 | double_to_table_curr(table, "aper_3", ap3); |
---|
1122 | double_to_table_curr(table, "aper_4", ap4); |
---|
1123 | double_to_table_curr(table, "on_ap", on_ap); |
---|
1124 | double_to_table_curr(table, "on_elem", on_elem); |
---|
1125 | double_to_table_curr(table, "spec", spec); |
---|
1126 | double_to_table_curr(table, "s", s); |
---|
1127 | double_to_table_curr(table, "x", x); |
---|
1128 | double_to_table_curr(table, "y", y); |
---|
1129 | double_to_table_curr(table, "betx", betx); |
---|
1130 | double_to_table_curr(table, "bety", bety); |
---|
1131 | double_to_table_curr(table, "dx", dx); |
---|
1132 | double_to_table_curr(table, "dy", dy); |
---|
1133 | |
---|
1134 | augment_count(table); |
---|
1135 | } |
---|
1136 | |
---|
1137 | static double |
---|
1138 | aper_calc(double p, double q, double* minhl, double halox[], double haloy[], |
---|
1139 | int halolength,double haloxadj[],double haloyadj[], |
---|
1140 | double newhalox[], double newhaloy[], double pipex[], double pipey[], |
---|
1141 | int pipelength, double notsimple) |
---|
1142 | { |
---|
1143 | int i=0, j=0, c=0, ver1, ver2; |
---|
1144 | double dist_limit=0.0000000001; |
---|
1145 | double a1, b1, a2, b2, xm, ym, h, l; |
---|
1146 | |
---|
1147 | for (c=0;c<=halolength+1;c++) |
---|
1148 | { |
---|
1149 | haloxadj[c]=halox[c]+p; |
---|
1150 | haloyadj[c]=haloy[c]+q; |
---|
1151 | } |
---|
1152 | |
---|
1153 | c=0; |
---|
1154 | |
---|
1155 | /*if halo centre is inside beam pipe, calculate smallest H/L ratio*/ |
---|
1156 | if (aper_chk_inside(p, q, pipex, pipey, dist_limit, pipelength)) |
---|
1157 | { |
---|
1158 | if (notsimple) |
---|
1159 | { |
---|
1160 | /*Adds extra apexes first:*/ |
---|
1161 | for (j=0;j<=halolength;j++) |
---|
1162 | { |
---|
1163 | newhalox[c]=haloxadj[j]; |
---|
1164 | newhaloy[c]=haloyadj[j]; |
---|
1165 | c++; |
---|
1166 | |
---|
1167 | for (i=0;i<=pipelength;i++) |
---|
1168 | { |
---|
1169 | /*Find a and b parameters for line*/ |
---|
1170 | ver1=aper_linepar(p, q, pipex[i], pipey[i], &a1, &b1); |
---|
1171 | ver2=aper_linepar(haloxadj[j], haloyadj[j], |
---|
1172 | haloxadj[j+1], haloyadj[j+1], &a2, &b2); |
---|
1173 | |
---|
1174 | /*find meeting coordinates for infinitely long lines*/ |
---|
1175 | aper_intersect(a1, b1, a2, b2, pipex[i], pipey[i], |
---|
1176 | haloxadj[j], haloyadj[j], ver1, ver2, &xm, &ym); |
---|
1177 | |
---|
1178 | /*eliminate intersection points not between line limits*/ |
---|
1179 | if (-1 == aper_online(xm, ym, haloxadj[j], haloyadj[j], |
---|
1180 | haloxadj[j+1], haloyadj[j+1], dist_limit)) /*halo line*/ |
---|
1181 | { |
---|
1182 | if (-1 != aper_online(p, q, pipex[i], pipey[i], xm, ym, |
---|
1183 | dist_limit)) /*test line*/ |
---|
1184 | { |
---|
1185 | newhalox[c]=xm; |
---|
1186 | newhaloy[c]=ym; |
---|
1187 | c++; |
---|
1188 | } |
---|
1189 | } |
---|
1190 | } |
---|
1191 | } |
---|
1192 | |
---|
1193 | halolength=c-1; |
---|
1194 | for (j=0;j<=halolength;j++) |
---|
1195 | { |
---|
1196 | haloxadj[j]=newhalox[j]; |
---|
1197 | haloyadj[j]=newhaloy[j]; |
---|
1198 | } |
---|
1199 | |
---|
1200 | } |
---|
1201 | |
---|
1202 | /*Calculates smallest ratio:*/ |
---|
1203 | for (i=0;i<=pipelength;i++) |
---|
1204 | { |
---|
1205 | for (j=0;j<=halolength;j++) |
---|
1206 | { |
---|
1207 | /*Find a and b parameters for line*/ |
---|
1208 | ver1=aper_linepar(p, q, haloxadj[j], haloyadj[j], &a1, &b1); |
---|
1209 | ver2=aper_linepar(pipex[i], pipey[i], pipex[i+1], pipey[i+1], &a2, &b2); |
---|
1210 | |
---|
1211 | /*find meeting coordinates for infinitely long lines*/ |
---|
1212 | aper_intersect(a1, b1, a2, b2, haloxadj[j], haloyadj[j], |
---|
1213 | pipex[i], pipey[i], ver1, ver2, &xm, &ym); |
---|
1214 | |
---|
1215 | /*eliminate intersection points not between line limits*/ |
---|
1216 | if (-1 == aper_online(xm, ym, pipex[i], pipey[i], pipex[i+1], pipey[i+1], |
---|
1217 | dist_limit)) /*pipe line*/ |
---|
1218 | { |
---|
1219 | if (-1 != aper_online(p, q, haloxadj[j], haloyadj[j], xm, ym, |
---|
1220 | dist_limit)) /*test line*/ |
---|
1221 | { |
---|
1222 | h=sqrt((xm-p)*(xm-p)+(ym-q)*(ym-q)); |
---|
1223 | l=sqrt((haloxadj[j]-p)*(haloxadj[j]-p) |
---|
1224 | + (haloyadj[j]-q)*(haloyadj[j]-q)); |
---|
1225 | if (h/l < *minhl) |
---|
1226 | { |
---|
1227 | *minhl=h/l; |
---|
1228 | } |
---|
1229 | } |
---|
1230 | } |
---|
1231 | } |
---|
1232 | } |
---|
1233 | } |
---|
1234 | else /*if halo centre is outside of beam pipe*/ |
---|
1235 | { |
---|
1236 | *minhl=0; |
---|
1237 | return -1; |
---|
1238 | } |
---|
1239 | |
---|
1240 | return 0; |
---|
1241 | } |
---|
1242 | |
---|
1243 | // public interface |
---|
1244 | |
---|
1245 | double |
---|
1246 | get_apertol(struct node* node, char* par) |
---|
1247 | /* returns aper_tol parameter 'i' where i is integer at the end of par; |
---|
1248 | e.g. aptol_1 gives i = 1 etc. (count starts at 1) */ |
---|
1249 | { |
---|
1250 | int i, k, n = strlen(par); |
---|
1251 | double val = zero, vec[100]; |
---|
1252 | for (i = 0; i < n; i++) if(isdigit(par[i])) break; |
---|
1253 | if (i == n) return val; |
---|
1254 | sscanf(&par[i], "%d", &k); k--; |
---|
1255 | if ((n = element_vector(node->p_elem, "aper_tol", vec)) > k) val = vec[k]; |
---|
1256 | return val; |
---|
1257 | } |
---|
1258 | |
---|
1259 | double |
---|
1260 | get_aperture(struct node* node, char* par) |
---|
1261 | /* returns aperture parameter 'i' where i is integer at the end of par; |
---|
1262 | e.g. aper_1 gives i = 1 etc. (count starts at 1) */ |
---|
1263 | { |
---|
1264 | int i, k, n = strlen(par); |
---|
1265 | double val = zero, vec[100]; |
---|
1266 | for (i = 0; i < n; i++) if(isdigit(par[i])) break; |
---|
1267 | if (i == n) return val; |
---|
1268 | sscanf(&par[i], "%d", &k); k--; |
---|
1269 | if ((n = element_vector(node->p_elem, "aperture", vec)) > k) val = vec[k]; |
---|
1270 | return val; |
---|
1271 | } |
---|
1272 | |
---|
1273 | void |
---|
1274 | pro_aperture(struct in_cmd* cmd) |
---|
1275 | { |
---|
1276 | struct node *use_range[2]; |
---|
1277 | struct table* tw_cp; |
---|
1278 | char *file, *range, tw_name[NAME_L], table_ap[NAME_L]="aperture", *table=table_ap; |
---|
1279 | int tw_cnt, rows; |
---|
1280 | double interval; |
---|
1281 | |
---|
1282 | embedded_twiss_cmd = cmd; |
---|
1283 | |
---|
1284 | /* check for valid sequence, beam and Twiss table */ |
---|
1285 | if (current_sequ != NULL && current_sequ->ex_start != NULL && sequence_length(current_sequ) != 0) |
---|
1286 | { |
---|
1287 | if (attach_beam(current_sequ) == 0) |
---|
1288 | fatal_error("Aperture module - sequence without beam:", current_sequ->name); |
---|
1289 | } |
---|
1290 | else fatal_error("Aperture module - no active sequence", ""); |
---|
1291 | |
---|
1292 | if (!sequ_check_valid_twiss(current_sequ)) |
---|
1293 | { |
---|
1294 | warning("Aperture module, no valid TWISS table present", "Aperture command ignored"); |
---|
1295 | return; |
---|
1296 | } |
---|
1297 | |
---|
1298 | range = command_par_string("range", this_cmd->clone); |
---|
1299 | if (get_ex_range(range, current_sequ, use_range) == 0) |
---|
1300 | { |
---|
1301 | warning("Illegal range.","Aperture command ignored"); |
---|
1302 | return; |
---|
1303 | } |
---|
1304 | |
---|
1305 | current_node = use_range[0]; |
---|
1306 | |
---|
1307 | /* navigate to starting point in Twiss table */ |
---|
1308 | tw_cp=current_sequ->tw_table; |
---|
1309 | tw_cnt=1; /* table starts at 1 */ |
---|
1310 | |
---|
1311 | while (1) { |
---|
1312 | string_from_table_row(tw_cp->name, "name", &tw_cnt, tw_name); |
---|
1313 | aper_trim_ws(tw_name, NAME_L); |
---|
1314 | |
---|
1315 | if (!strcmp(tw_name,current_node->name)) break; |
---|
1316 | |
---|
1317 | if (++tw_cnt > tw_cp->curr) { |
---|
1318 | warning("Aperture command ignored: could not find range start in Twiss table:", current_node->name); |
---|
1319 | return; |
---|
1320 | } |
---|
1321 | } |
---|
1322 | |
---|
1323 | /* approximate # of needed rows in aperture table */ |
---|
1324 | interval = command_par_value("interval", this_cmd->clone); |
---|
1325 | rows = current_sequ->n_nodes + 2 * (sequence_length(current_sequ)/interval); |
---|
1326 | |
---|
1327 | /* make empty aperture table */ |
---|
1328 | aperture_table=make_table(table, table, ap_table_cols, ap_table_types, rows); |
---|
1329 | aperture_table->dynamic=1; |
---|
1330 | add_to_table_list(aperture_table, table_register); |
---|
1331 | |
---|
1332 | /* calculate apertures and fill table */ |
---|
1333 | struct aper_node limit_node = { "none", -1, -1, "none", {-1,-1,-1,-1}, {-1,-1,-1}, 0.0 }; |
---|
1334 | struct aper_node *limit_pt = &limit_node; |
---|
1335 | |
---|
1336 | limit_pt = aperture(table, use_range, tw_cp, &tw_cnt, limit_pt); |
---|
1337 | |
---|
1338 | if (limit_pt->n1 != -1) { |
---|
1339 | // printf("\n\nAPERTURE LIMIT: %s, n1: %g, at: %g\n\n", limit_pt->name, limit_pt->n1, limit_pt->s); |
---|
1340 | |
---|
1341 | aper_header(aperture_table, limit_pt); |
---|
1342 | |
---|
1343 | file = command_par_string("file", this_cmd->clone); |
---|
1344 | if (file != NULL) { |
---|
1345 | out_table(table, aperture_table, file); |
---|
1346 | } |
---|
1347 | if (strcmp(aptwfile,"dummy")) out_table(tw_cp->name, tw_cp, aptwfile); |
---|
1348 | } |
---|
1349 | else warning("Could not run aperture command.","Aperture command ignored"); |
---|
1350 | |
---|
1351 | /* set pointer to updated Twiss table */ |
---|
1352 | current_sequ->tw_table=tw_cp; |
---|
1353 | } |
---|
1354 | |
---|
1355 | static struct aper_node* |
---|
1356 | aperture(char *table, struct node* use_range[], struct table* tw_cp, int *tw_cnt, struct aper_node *lim_pt) |
---|
1357 | { |
---|
1358 | int stop=0, nint=1, jslice=1, first, ap=1; // , err not used |
---|
1359 | int true_flag, true_node=0, offs_node=0, do_survey=0; |
---|
1360 | int truepos=0, true_cnt=0, offs_cnt=0; |
---|
1361 | int halo_q_length=1, halolength, pipelength, namelen=NAME_L, ntol; // nhalopar, not used |
---|
1362 | double surv_init[6]={0, 0, 0, 0, 0, 0}; |
---|
1363 | double surv_x=zero, surv_y=zero, elem_x=0, elem_y=0; |
---|
1364 | double xa=0, xb=0, xc=0, ya=0, yb=0, yc=0; |
---|
1365 | double on_ap=1, on_elem=0; |
---|
1366 | double mass, energy, exn, eyn, dqf, betaqfx, dp, dparx, dpary; |
---|
1367 | double cor, bbeat, nco, halo[4], interval, spec, ex, ey, notsimple; |
---|
1368 | double s=0, x=0, y=0, betx=0, bety=0, dx=0, dy=0, ratio, n1, nr, length; |
---|
1369 | double xeff=0,yeff=0; |
---|
1370 | double n1x_m, n1y_m; |
---|
1371 | double s_start, s_curr, s_end; |
---|
1372 | double node_s=-1, node_n1=-1; |
---|
1373 | double aper_tol[3], ap1, ap2, ap3, ap4; |
---|
1374 | double dispx, dispy, tolx, toly; |
---|
1375 | double dispxadj=0, dispyadj=0, coxadj, coyadj, tolxadj=0, tolyadj=0; |
---|
1376 | double angle, dangle, deltax, deltay; |
---|
1377 | double xshift, yshift, r; |
---|
1378 | double halox[MAXARRAY], haloy[MAXARRAY], haloxsi[MAXARRAY], haloysi[MAXARRAY]; |
---|
1379 | double haloxadj[MAXARRAY], haloyadj[MAXARRAY], newhalox[MAXARRAY], newhaloy[MAXARRAY]; |
---|
1380 | double pipex[MAXARRAY], pipey[MAXARRAY]; |
---|
1381 | double parxd,paryd; |
---|
1382 | char *halofile, *truefile, *offsfile; |
---|
1383 | char refnode[NAME_L]; |
---|
1384 | char *cmd_refnode; |
---|
1385 | char apertype[NAME_L]; |
---|
1386 | char name[NAME_L]; |
---|
1387 | char tol_err_mess[80] = ""; |
---|
1388 | |
---|
1389 | struct node* rng_glob[2]; |
---|
1390 | // struct aper_node limit_node = {"none", -1, -1, "none", {-1,-1,-1,-1}, {-1,-1,-1}}; |
---|
1391 | // struct aper_node* lim_pt = &limit_node; |
---|
1392 | |
---|
1393 | int is_zero_len; |
---|
1394 | |
---|
1395 | true_tab = mycalloc("Aperture",E_D_LIST_CHUNK,sizeof(struct aper_e_d) ); |
---|
1396 | /* offs_tab = (struct aper_e_d*) mycalloc("Aperture",E_D_LIST_CHUNK,sizeof(struct aper_e_d));*/ |
---|
1397 | |
---|
1398 | printf("\nProcessing apertures from %s to %s...\n",use_range[0]->name,use_range[1]->name); |
---|
1399 | |
---|
1400 | /* read command parameters */ |
---|
1401 | halofile = command_par_string("halofile", this_cmd->clone); |
---|
1402 | /* removed IW 240205 */ |
---|
1403 | /* pipefile = command_par_string("pipefile", this_cmd->clone); */ |
---|
1404 | exn = command_par_value("exn", this_cmd->clone); |
---|
1405 | eyn = command_par_value("eyn", this_cmd->clone); |
---|
1406 | dqf = command_par_value("dqf", this_cmd->clone); |
---|
1407 | betaqfx = command_par_value("betaqfx", this_cmd->clone); |
---|
1408 | dp = command_par_value("dp", this_cmd->clone); |
---|
1409 | dparx = command_par_value("dparx", this_cmd->clone); |
---|
1410 | dpary = command_par_value("dpary", this_cmd->clone); |
---|
1411 | cor = command_par_value("cor", this_cmd->clone); |
---|
1412 | bbeat = command_par_value("bbeat", this_cmd->clone); |
---|
1413 | nco = command_par_value("nco", this_cmd->clone); |
---|
1414 | command_par_vector("halo", this_cmd->clone, halo); // nhalopar = // not used |
---|
1415 | interval = command_par_value("interval", this_cmd->clone); |
---|
1416 | spec = command_par_value("spec", this_cmd->clone); |
---|
1417 | notsimple = command_par_value("notsimple", this_cmd->clone); |
---|
1418 | truefile = command_par_string("trueprofile", this_cmd->clone); |
---|
1419 | offsfile = command_par_string("offsetelem", this_cmd->clone); |
---|
1420 | |
---|
1421 | cmd_refnode = command_par_string("refnode", this_cmd->clone); |
---|
1422 | |
---|
1423 | mass = get_value("beam", "mass"); |
---|
1424 | energy = get_value("beam", "energy"); |
---|
1425 | |
---|
1426 | /* fetch deltap as set by user in the former TWISS command */ |
---|
1427 | /* will be used below for displacement associated to parasitic dipersion */ |
---|
1428 | |
---|
1429 | double_from_table_row("summ","deltap",&nint,&lim_pt->deltap_twiss); |
---|
1430 | printf ("+++++++ deltap from TWISS %12.6g\n",lim_pt->deltap_twiss); |
---|
1431 | |
---|
1432 | /* calculate emittance and delta angle */ |
---|
1433 | ex=mass*exn/energy; ey=mass*eyn/energy; |
---|
1434 | dangle=twopi/(nco*4); |
---|
1435 | |
---|
1436 | /* check if trueprofile and offsetelem files exist */ |
---|
1437 | true_flag = aper_e_d_read(truefile, &true_tab, &true_cnt, refnode); |
---|
1438 | /* offs_flag = aper_e_d_read(offsfile, &offs_tab, &offs_cnt, refnode);*/ |
---|
1439 | offs_tab = aper_e_d_read_tfs(offsfile, &offs_cnt, refnode); |
---|
1440 | |
---|
1441 | if (cmd_refnode != NULL) { |
---|
1442 | strcpy(refnode, cmd_refnode); |
---|
1443 | strcat(refnode, ":1"); |
---|
1444 | } |
---|
1445 | |
---|
1446 | printf("\nreference node: %s\n",refnode); |
---|
1447 | |
---|
1448 | /* build halo polygon based on input ratio values or coordinates */ |
---|
1449 | if ((halolength = aper_external_file(halofile, halox, haloy)) > -1) |
---|
1450 | ; |
---|
1451 | else if (aper_rectellipse(&halo[2], &halo[3], &halo[1], &halo[1], &halo_q_length, halox, haloy)) { |
---|
1452 | warning("Not valid parameters for halo. ", "Unable to make polygon."); |
---|
1453 | |
---|
1454 | /* IA */ |
---|
1455 | myfree("Aperture",true_tab); |
---|
1456 | myfree("Aperture",offs_tab); |
---|
1457 | |
---|
1458 | return lim_pt; |
---|
1459 | } |
---|
1460 | else aper_fill_quads(halox, haloy, halo_q_length, &halolength); |
---|
1461 | |
---|
1462 | /* check for externally given pipe polygon */ |
---|
1463 | /* changed this recently, IW 240205 */ |
---|
1464 | /* pipelength = aper_external_file(pipefile, pipex, pipey); |
---|
1465 | if ( pipelength > -1) ext_pipe=1; */ |
---|
1466 | |
---|
1467 | /* get initial twiss parameters, from start of first element in range */ |
---|
1468 | aper_read_twiss(tw_cp->name, tw_cnt, &s_end, &x, &y, &betx, &bety, &dx, &dy); |
---|
1469 | // LD: shift further results by one step (?) and finish outside the table |
---|
1470 | // (*tw_cnt)++; |
---|
1471 | aper_adj_halo_si(ex, ey, betx, bety, bbeat, halox, haloy, halolength, haloxsi, haloysi); |
---|
1472 | |
---|
1473 | /* calculate initial normal+parasitic disp. */ |
---|
1474 | /* modified 27feb08 BJ */ |
---|
1475 | parxd = dparx*sqrt(betx/betaqfx)*dqf; |
---|
1476 | paryd = dpary*sqrt(bety/betaqfx)*dqf; |
---|
1477 | |
---|
1478 | /* Initialize n1 limit value */ |
---|
1479 | lim_pt->n1=999999; |
---|
1480 | |
---|
1481 | int n = 0; |
---|
1482 | |
---|
1483 | while (!stop) { |
---|
1484 | ++n; |
---|
1485 | |
---|
1486 | strcpy(name,current_node->name); |
---|
1487 | aper_trim_ws(name, NAME_L); |
---|
1488 | |
---|
1489 | is_zero_len = 0; |
---|
1490 | |
---|
1491 | /* the first node in a sequence can not be sliced, hence: */ |
---|
1492 | if (current_sequ->range_start == current_node) first=1; else first=0; |
---|
1493 | |
---|
1494 | length=node_value("l"); |
---|
1495 | double_from_table_row(tw_cp->name, "s", tw_cnt, &s_end); |
---|
1496 | s_start=s_end-length; |
---|
1497 | s_curr=s_start; |
---|
1498 | |
---|
1499 | node_string("apertype", apertype, &namelen); |
---|
1500 | aper_trim_ws(apertype, NAME_L); |
---|
1501 | |
---|
1502 | if (!strncmp("drift",name,5)) on_elem=-999999; |
---|
1503 | else on_elem=1; |
---|
1504 | |
---|
1505 | if ( (offs_tab != NULL) && (strcmp(refnode, name) == 0)) do_survey=1; |
---|
1506 | /* printf("\nname: %s, ref: %s, do_survey:: %d\n",name,refnode, do_survey);*/ |
---|
1507 | |
---|
1508 | /* read data for tol displacement of halo */ |
---|
1509 | get_node_vector("aper_tol",&ntol,aper_tol); |
---|
1510 | if (ntol == 3) { |
---|
1511 | r = aper_tol[0]; |
---|
1512 | xshift = aper_tol[1]; |
---|
1513 | yshift = aper_tol[2]; |
---|
1514 | } |
---|
1515 | else r=xshift=yshift=0; |
---|
1516 | |
---|
1517 | /*read aperture data and make polygon tables for beam pipe*/ |
---|
1518 | /* IW 250205 */ |
---|
1519 | /* if (ext_pipe == 0) */ |
---|
1520 | ap=aper_bs(apertype, &ap1, &ap2, &ap3, &ap4, &pipelength, pipex, pipey); |
---|
1521 | |
---|
1522 | if (ap == 0 || first == 1) { |
---|
1523 | /* if no pipe can be built, the n1 is set to inf and Twiss parms read for reference*/ |
---|
1524 | n1=999999; n1x_m=999999; n1y_m=999999; on_ap=-999999; nint=1; |
---|
1525 | |
---|
1526 | aper_read_twiss(tw_cp->name, tw_cnt, &s_end, &x, &y, &betx, &bety, &dx, &dy); |
---|
1527 | |
---|
1528 | aper_write_table(name, &n1, &n1x_m, &n1y_m, &r, &xshift, &yshift, apertype, |
---|
1529 | &ap1, &ap2, &ap3, &ap4, &on_ap, &on_elem, &spec, |
---|
1530 | &s_end, &x, &y, &betx, &bety, &dx, &dy, table); |
---|
1531 | on_ap=1; |
---|
1532 | |
---|
1533 | double_to_table_row(tw_cp->name, "n1", tw_cnt, &n1); |
---|
1534 | (*tw_cnt)++; |
---|
1535 | |
---|
1536 | /* calc disp and adj halo to have ready for next node */ |
---|
1537 | /* modified 27feb08 BJ */ |
---|
1538 | parxd = dparx*sqrt(betx/betaqfx)*dqf; |
---|
1539 | paryd = dpary*sqrt(bety/betaqfx)*dqf; |
---|
1540 | |
---|
1541 | aper_adj_halo_si(ex, ey, betx, bety, bbeat, halox, haloy, halolength, haloxsi, haloysi); |
---|
1542 | |
---|
1543 | /*do survey to have ready init for next node */ |
---|
1544 | if (do_survey) { |
---|
1545 | rng_glob[0] = current_sequ->range_start; |
---|
1546 | rng_glob[1] = current_sequ->range_end; |
---|
1547 | current_sequ->range_start = current_sequ->range_end = current_node; |
---|
1548 | aper_surv(surv_init, nint); |
---|
1549 | double_from_table_row("survey","x",&nint, &surv_x); |
---|
1550 | double_from_table_row("survey","y",&nint, &surv_y); |
---|
1551 | current_sequ->range_start = rng_glob[0]; |
---|
1552 | current_sequ->range_end = rng_glob[1]; |
---|
1553 | } |
---|
1554 | } /* end loop 'if no pipe ' */ |
---|
1555 | else { |
---|
1556 | node_n1=999999; |
---|
1557 | true_node=0; |
---|
1558 | offs_node=0; |
---|
1559 | |
---|
1560 | /* calculate the number of slices per node */ |
---|
1561 | if (true_flag == 0) |
---|
1562 | nint=length/interval; |
---|
1563 | else { |
---|
1564 | true_node=aper_tab_search(true_cnt, true_tab, name, &truepos); |
---|
1565 | |
---|
1566 | if (true_node) nint=true_tab[truepos].curr; |
---|
1567 | else nint=length/interval; |
---|
1568 | } |
---|
1569 | /* printf("\nname: %s, nint: %d",name,nint); */ |
---|
1570 | |
---|
1571 | if (!nint) nint=1; |
---|
1572 | |
---|
1573 | /* don't interpolate 0-length elements*/ |
---|
1574 | |
---|
1575 | if (fabs(length) < MIN_DOUBLE ) is_zero_len = 1; |
---|
1576 | |
---|
1577 | /* slice the node, call survey if necessary, make twiss for slices*/ |
---|
1578 | interpolate_node(&nint); |
---|
1579 | |
---|
1580 | /* do survey */ |
---|
1581 | if (do_survey) { |
---|
1582 | double offs_row[8] = { 0 }; |
---|
1583 | |
---|
1584 | aper_surv(surv_init, nint); |
---|
1585 | |
---|
1586 | offs_node = aper_tab_search_tfs(offs_tab, name, offs_row); |
---|
1587 | if (offs_node) { |
---|
1588 | /* printf("\nusing offset");*/ |
---|
1589 | xa=offs_row[4]; |
---|
1590 | xb=offs_row[3]; |
---|
1591 | xc=offs_row[2]; |
---|
1592 | ya=offs_row[7]; |
---|
1593 | yb=offs_row[6]; |
---|
1594 | yc=offs_row[5]; |
---|
1595 | } |
---|
1596 | } |
---|
1597 | |
---|
1598 | embedded_twiss(); |
---|
1599 | |
---|
1600 | /* Treat each slice, for all angles */ |
---|
1601 | for (jslice=0; jslice <= nint; jslice++) { |
---|
1602 | ratio=999999; |
---|
1603 | if (jslice != 0) { |
---|
1604 | aper_read_twiss("embedded_twiss_table", &jslice, &s, &x, &y, &betx, &bety, &dx, &dy); |
---|
1605 | |
---|
1606 | s_curr=s_start+s; |
---|
1607 | aper_adj_halo_si(ex, ey, betx, bety, bbeat, halox, haloy, halolength, haloxsi, haloysi); |
---|
1608 | |
---|
1609 | /* calculate normal+parasitic disp.*/ |
---|
1610 | /* modified 27feb08 BJ */ |
---|
1611 | parxd = dparx*sqrt(betx/betaqfx)*dqf; |
---|
1612 | paryd = dpary*sqrt(bety/betaqfx)*dqf; |
---|
1613 | |
---|
1614 | if (do_survey) { |
---|
1615 | double_from_table_row("survey","x",&jslice, &surv_x); |
---|
1616 | double_from_table_row("survey","y",&jslice, &surv_y); |
---|
1617 | } |
---|
1618 | } |
---|
1619 | else { |
---|
1620 | /// jslice==0, parameters from previous node will be used |
---|
1621 | s_curr+=1.e-12; /*to get correct plot at start of elements*/ |
---|
1622 | s=0; /*used to calc elem_x elem_y) */ |
---|
1623 | } |
---|
1624 | |
---|
1625 | xeff = x; |
---|
1626 | yeff = y; |
---|
1627 | |
---|
1628 | /* survey adjustments */ |
---|
1629 | /* BJ 3 APR 2009 : introduced xeff and yeff in order to avoid |
---|
1630 | interferences with x and y as used for the first slice |
---|
1631 | (re-use from end of former node) */ |
---|
1632 | if (offs_node) { |
---|
1633 | elem_x=xa*s*s+xb*s+xc; |
---|
1634 | elem_y=ya*s*s+yb*s+yc; |
---|
1635 | xeff=x+(surv_x-elem_x); |
---|
1636 | yeff=y+(surv_y-elem_y); |
---|
1637 | } |
---|
1638 | |
---|
1639 | /* discrete adjustments */ |
---|
1640 | if (true_node) { |
---|
1641 | xeff+=true_tab[truepos].tab[jslice][1]; |
---|
1642 | yeff+=true_tab[truepos].tab[jslice][2]; |
---|
1643 | } |
---|
1644 | |
---|
1645 | for (angle=0; angle < twopi; angle += dangle) { |
---|
1646 | /* new 27feb08 BJ */ |
---|
1647 | dispx = bbeat*(fabs(dx)*dp + parxd*(fabs(lim_pt->deltap_twiss)+dp) ); |
---|
1648 | dispy = bbeat*(fabs(dy)*dp + paryd*(fabs(lim_pt->deltap_twiss)+dp) ); |
---|
1649 | |
---|
1650 | /*adjust dispersion to worst-case for quadrant*/ |
---|
1651 | aper_adj_quad(angle, dispx, dispy, &dispxadj, &dispyadj); |
---|
1652 | |
---|
1653 | /*calculate displacement co+tol for each angle*/ |
---|
1654 | coxadj=cor*cos(angle); coyadj=cor*sin(angle); |
---|
1655 | |
---|
1656 | /* Error check added 20feb08 BJ */ |
---|
1657 | if ( xshift<0 || yshift < 0 || r<0 ) { |
---|
1658 | sprintf(tol_err_mess,"In element : %s\n",name); |
---|
1659 | fatal_error("Illegal negative tolerance",tol_err_mess); |
---|
1660 | } |
---|
1661 | |
---|
1662 | aper_race(xshift,yshift,r,angle,&tolx,&toly); |
---|
1663 | aper_adj_quad(angle, tolx, toly, &tolxadj, &tolyadj); |
---|
1664 | |
---|
1665 | /* add all displacements */ |
---|
1666 | deltax = coxadj + tolxadj + xeff + dispxadj; |
---|
1667 | deltay = coyadj + tolyadj + yeff + dispyadj; |
---|
1668 | |
---|
1669 | /* send beta adjusted halo and its displacement to aperture calculation */ |
---|
1670 | aper_calc(deltax,deltay,&ratio,haloxsi,haloysi, |
---|
1671 | halolength,haloxadj,haloyadj,newhalox,newhaloy, |
---|
1672 | pipex,pipey,pipelength,notsimple); |
---|
1673 | } |
---|
1674 | |
---|
1675 | nr=ratio*halo[1]; |
---|
1676 | n1=nr/(halo[1]/halo[0]); /* ratio r/n = 1.4 */ |
---|
1677 | |
---|
1678 | n1x_m=n1*bbeat*sqrt(betx*ex); |
---|
1679 | n1y_m=n1*bbeat*sqrt(bety*ey); |
---|
1680 | |
---|
1681 | /* Change below, BJ 23oct2008 */ |
---|
1682 | /* test block 'if (n1 < node_n1)' included in test block */ |
---|
1683 | |
---|
1684 | if (is_zero_len == 0 || jslice == 1) { |
---|
1685 | aper_write_table(name, &n1, &n1x_m, &n1y_m, &r, &xshift, &yshift, apertype, |
---|
1686 | &ap1, &ap2, &ap3, &ap4, &on_ap, &on_elem, &spec, &s_curr, |
---|
1687 | &xeff, &yeff, &betx, &bety, &dx, &dy, table); |
---|
1688 | |
---|
1689 | /* save node minimum n1 */ |
---|
1690 | |
---|
1691 | if (n1 < node_n1) { |
---|
1692 | node_n1=n1; |
---|
1693 | node_s=s_curr; |
---|
1694 | } |
---|
1695 | } // if is_zero_len |
---|
1696 | } // for jslice |
---|
1697 | |
---|
1698 | reset_interpolation(&nint); |
---|
1699 | |
---|
1700 | /* insert minimum node value into Twiss table */ |
---|
1701 | double_to_table_row(tw_cp->name, "n1", tw_cnt, &node_n1); |
---|
1702 | (*tw_cnt)++; |
---|
1703 | |
---|
1704 | /* save range minimum n1 */ |
---|
1705 | if (node_n1 < lim_pt->n1) { |
---|
1706 | strcpy(lim_pt->name,name); |
---|
1707 | lim_pt->n1=node_n1; |
---|
1708 | lim_pt->s=node_s; |
---|
1709 | strcpy(lim_pt->apertype,apertype); |
---|
1710 | lim_pt->aperture[0]=ap1; |
---|
1711 | lim_pt->aperture[1]=ap2; |
---|
1712 | lim_pt->aperture[2]=ap3; |
---|
1713 | lim_pt->aperture[3]=ap4; |
---|
1714 | lim_pt->aper_tol[0]=r; |
---|
1715 | lim_pt->aper_tol[1]=xshift; |
---|
1716 | lim_pt->aper_tol[2]=yshift; |
---|
1717 | } |
---|
1718 | } // if !(ap == 0 || first == 1) |
---|
1719 | |
---|
1720 | if (!strcmp(current_node->name,use_range[1]->name)) stop=1; |
---|
1721 | if (!advance_node()) stop=1; |
---|
1722 | } // while !stop |
---|
1723 | |
---|
1724 | myfree("Aperture",true_tab); |
---|
1725 | if (offs_tab != NULL) myfree("Aperture",offs_tab); |
---|
1726 | |
---|
1727 | return lim_pt; |
---|
1728 | } |
---|
1729 | |
---|