source: PSPA/madxPSPA/src/mad_mkthin.c @ 476

Last change on this file since 476 was 430, checked in by touze, 11 years ago

import madx-5.01.00

File size: 48.5 KB
Line 
1#include "madx.h"
2
3typedef unsigned char bool;
4
5/* makethin.c
6 Thick to thin lens converter. Helmut Burkhardt
7 Early versions in 2001, 2002 by Mark Hayes
8 */
9
10/* forward declarations of some routines which are only used within makethin */
11static struct element* create_thin_obj(struct element*, int slice_no);
12static struct sequence* seq_diet(struct sequence*);
13static double at_shift(int, int);
14static double q_shift(int, int);
15
16/* this structure is used to store a lookup table of thick to thin
17 element conversions already done */
18struct thin_lookup
19{
20  struct element *thick_elem;
21  struct element *thin_elem;
22  int slice;
23  struct thin_lookup *next;
24};
25
26static struct thin_lookup *my_list = NULL;
27
28/* this structure is used to store a lookup table of thick to thin
29 sequence conversions already done */
30struct thin_sequ_lookup
31{
32  struct sequence *thick_sequ;
33  struct sequence *thin_sequ;
34  struct thin_sequ_lookup *next;
35};
36
37static struct thin_sequ_lookup *my_sequ_list = NULL;
38
39/* this is used when choosing the style of slicing */
40static char* thin_style = NULL;
41static char collim_style[] = "collim";
42/* this is used to return how to slice the selected elements */
43static struct el_list *thin_select_list = NULL;
44
45// private interface
46
47static void
48force_consistent_slices(void)
49/* hbu 10/2005
50 loop over all elements and check that #slices of child and parent agree
51 if not, use the maximum for both
52 */
53{
54//  struct element* el_i;
55  struct command_parameter *child,*parent;
56  int i,el_i_slice_pos,slices,slices_parent;
57  for(i=0; i< element_list->curr; i++) /* loop over element_list */
58  {
59    struct element* el_i = element_list->elem[i];
60    el_i_slice_pos = name_list_pos("slice",el_i->def->par_names);
61    if(el_i_slice_pos>0 && el_i->parent!=NULL && el_i != el_i->parent )
62    {
63      child=el_i->def->par->parameters[el_i_slice_pos];
64      parent=el_i->parent->def->par->parameters[el_i_slice_pos];
65      slices=child->double_value;
66      slices_parent=parent->double_value;
67      if(slices != slices_parent)
68      {
69        if(slices>slices_parent) slices_parent=slices; else slices=slices_parent;
70        child->double_value=slices;
71        parent->double_value=slices_parent;
72      }
73    }
74  }
75}
76
77#if 0 // not used
78static void
79dump_slices(void)
80/* Loops over all current elements and prints the number of slices. Used for debug and info */
81{
82  struct element* el_i;
83  int i,el_i_slice_pos,slices,slices_parent,n_elem_with_slice=0,n_elem_with_slice_gt_1=0;
84  char* parent_name;
85  printf("++++++ dump_slices");
86  printf("            name #slices  derived from #slices\n");
87  for(i=0; i< element_list->curr; i++) /* loop over element_list */
88  {
89    el_i = element_list->elem[i];
90    el_i_slice_pos = name_list_pos("slice",el_i->def->par_names);
91    if(el_i_slice_pos>0)
92    {
93      n_elem_with_slice++;
94      slices=el_i->def->par->parameters[el_i_slice_pos]->double_value;
95      /* look also at parent if existing */
96      slices_parent=0;
97      parent_name="no parent";
98      if(el_i->parent!=NULL)
99      {
100        slices_parent=el_i->parent->def->par->parameters[el_i_slice_pos]->double_value;
101        parent_name=el_i->parent->name;
102      }
103      if(slices>1) n_elem_with_slice_gt_1++;
104    }
105  }
106  printf("------ end of dump slices. There were %4d elements, %3d with slice numbers and %2d with slice numbers>1\n\n",element_list->curr,n_elem_with_slice,n_elem_with_slice_gt_1);
107}
108#endif
109
110static int
111get_slices_from_elem(struct element* elem)
112{
113  int elem_slice_pos=0,slices=1;
114  elem_slice_pos = name_list_pos("slice",elem->def->par_names);
115  if(elem_slice_pos > 0)
116  {
117    slices=elem->def->par->parameters[elem_slice_pos]->double_value;
118  }
119  if (slices==0) slices = 1; /* must always slice to thin */
120  return slices;
121}
122
123/* Has this element already been dieted? returns NULL for NO.*/
124static struct element*
125get_thin(struct element* thick_elem, int slice)
126{
127  struct thin_lookup *cur;
128  if (my_list)
129  {
130    cur = my_list;
131    while (cur)
132    {
133      if (cur->thick_elem == thick_elem && cur->slice == slice)
134      {
135        return cur->thin_elem;
136      }
137      cur = cur->next;
138    }
139  }
140  return NULL;
141}
142
143/* Enter a newly dieted element */
144static void
145put_thin(struct element* thick_elem, struct element* thin_elem, int slice)
146{
147  struct thin_lookup *p,*cur;
148  char rout_name[] = "makethin:put_thin";
149  p = (struct thin_lookup*) mycalloc(rout_name,1, sizeof(struct thin_lookup));
150  p->thick_elem = thick_elem;
151  p->thin_elem = thin_elem;
152  p->slice = slice;
153  p->next = NULL;
154  if (my_list)
155  {
156    cur = my_list;
157    while (cur->next) cur = cur->next;
158    cur->next = p;
159  }
160  else
161  {
162    my_list = p;
163  }
164  return;
165}
166
167/* Has this sequence already been dieted? returns NULL for NO.*/
168static struct sequence*
169get_thin_sequ(struct sequence* thick_sequ)
170{
171  struct thin_sequ_lookup *cur;
172  if (my_sequ_list)
173  {
174    cur = my_sequ_list;
175    while (cur)
176    {
177      if (cur->thick_sequ == thick_sequ)
178      {
179        return cur->thin_sequ;
180      }
181      cur = cur->next;
182    }
183  }
184  return NULL;
185}
186
187/* Enter a newly dieted element */
188static void
189put_thin_sequ(struct sequence* thick_sequ, struct sequence* thin_sequ)
190{
191  struct thin_sequ_lookup *p,*cur;
192  char rout_name[] = "makethin:put_thin_sequ";
193  p = (struct thin_sequ_lookup*) mycalloc(rout_name,1, sizeof(struct thin_sequ_lookup));
194  p->thick_sequ = thick_sequ;
195  p->thin_sequ = thin_sequ;
196  p->next = NULL;
197  if (my_sequ_list)
198  {
199    cur = my_sequ_list;
200    while (cur->next) cur = cur->next;
201    cur->next = p;
202  }
203  else
204  {
205    my_sequ_list = p;
206  }
207  return;
208}
209
210/* makes node name from element name and slice number*/
211static char*
212make_thin_name(char* e_name, int slice)
213{
214  char c_dummy[128];
215  sprintf(c_dummy,"%s..%d", e_name, slice);
216  return buffer(c_dummy);
217}
218
219/* multiply the k by length and divide by slice */
220static struct command_parameter*
221scale_and_slice(struct command_parameter *kn_param, struct command_parameter *length_param, int slices, int slice_no, int angle_conversion, int kl_flag)
222{
223  int last_non_zero=-1,i;
224  struct expression *kn_i_expr;
225  double kn_i_val;
226  if (kn_param == NULL) return NULL;
227 
228  for (i=0; i<kn_param->expr_list->curr; i++)
229  {
230    kn_i_expr = kn_param->expr_list->list[i];
231    kn_i_val  = kn_param->double_array->a[i];
232    if ((kn_i_expr != NULL && zero_string(kn_i_expr->string)==0)  || kn_i_val!=0)
233    {
234      last_non_zero=i;
235      if (kl_flag == 0 && (angle_conversion==0||i>0)) /*hbu apply the angle_conversion==0 check only to zero order multipole */
236      {
237        if ((length_param->expr) || (kn_i_expr))
238        {
239          kn_i_expr = compound_expr(kn_i_expr,kn_i_val,"*",length_param->expr,length_param->double_value); /* multiply expression with length */
240        }
241        else
242        { /* multiply value with length */
243          kn_i_val *= length_param->double_value;
244        }
245      }
246      if (slices > 1)
247      { /* give the correct weight by slice (multiply with the inverse of the number of slices) */
248        if (kn_i_expr)
249        {
250          kn_i_expr = compound_expr(kn_i_expr,kn_i_val,"*",NULL,q_shift(slices,slice_no));
251        }
252        else
253        {
254          kn_i_val *= q_shift(slices,slice_no);
255        }
256      }
257    }
258    if(kn_i_expr) kn_param->expr_list->list[i] = kn_i_expr;
259    kn_param->double_array->a[i] = kn_i_val;
260  } /* for i ..*/
261  if (last_non_zero==-1)
262  {
263    delete_command_parameter(kn_param); kn_param=NULL;
264  }
265  return kn_param;
266}
267
268/* translate k0,k1,k2,k3 & k0s,k1s,k2s,k3s to kn{} and ks{} */
269/* 26/11/01 - removed tilt param */
270static int
271translate_k(struct command_parameter* *kparam,
272            struct command_parameter* *ksparam,
273            struct command_parameter  *angle_param,
274            struct command_parameter  *kn_param,
275            struct command_parameter  *ks_param)
276{
277  int i,angle_conversion=0;
278  /*    char *zero[1]; */
279  /*    zero[0] = buffer("0"); */
280 
281  if ((kparam == NULL) && (ksparam == NULL))
282    fatal_error("translate_k: no kparams to convert","");
283 
284  /* if we have a angle we ignore any given k0 */
285  if (angle_param)
286  {
287    kparam[0] =  new_command_parameter("k0", 2);
288    angle_conversion=1; /* note we do not divide by length, just to multiply again afterwards */
289    if (angle_param->expr)
290    {
291      kparam[0]->expr =  clone_expression(angle_param->expr);
292    }
293    kparam[0]->double_value = angle_param->double_value;
294  }
295 
296  for (i=0; i<4; i++)
297  {
298    /* zero all the parameters */
299    kn_param->expr_list->list[i] = NULL; kn_param->double_array->a[i] = 0;
300    ks_param->expr_list->list[i] = NULL; ks_param->double_array->a[i] = 0;
301    /* copy across the k's */
302    if (kparam[i])
303    {
304      if (kparam[i]->expr)
305      {
306        kn_param->expr_list->list[i] = clone_expression(kparam[i]->expr);
307      }
308      kn_param->double_array->a[i] = kparam[i]->double_value;
309    }
310    if (ksparam[i])
311    {
312      if (ksparam[i]->expr)
313      {
314        ks_param->expr_list->list[i] = clone_expression(ksparam[i]->expr);
315      }
316      ks_param->double_array->a[i] = ksparam[i]->double_value;
317    }
318    /* update the number of k's in our arrays */
319    kn_param->expr_list->curr++; kn_param->double_array->curr++;
320    ks_param->expr_list->curr++; ks_param->double_array->curr++;
321  }
322 
323  return angle_conversion;
324}
325
326/* adds a node to the end of a sequence */
327static void
328seq_diet_add(struct node* node, struct sequence* sequ)
329{
330  if (sequ->start == NULL)
331  { /* first node in new sequence? */
332    sequ->start = node;
333    sequ->end = node;
334    node->next = NULL;
335    node->previous = NULL;
336  }
337  else
338  { /* no? then add to end */
339    sequ->end->next = node;
340    node->previous  = sequ->end;
341    sequ->end = node;
342  }
343  add_to_node_list(node, 0, sequ->nodes);
344 
345  return;
346}
347
348/* adds a sequence to a sequence */
349static void
350seq_diet_add_sequ(struct node* thick_node, struct sequence* sub_sequ, struct sequence* sequ)
351{
352  struct node* node = new_sequ_node(sub_sequ, thick_node->occ_cnt); /* 1 is the occ_cnt*/
353  node->length = 0;
354  node->at_value = thick_node->at_value;
355  if (node->at_expr)
356  {
357    node->at_expr = clone_expression(thick_node->at_expr);
358  }
359  seq_diet_add(node,sequ);
360  return;
361}
362
363static void
364add_lrad(struct command* cmd,struct command_parameter *length_param,int slices)
365{
366  struct command_parameter *l_par;
367  if(length_param)
368  {
369    add_cmd_parameter_new(cmd,0.,"l",1); /* new parameter l with value of 0 */
370    l_par = cmd->par->parameters[cmd->par->curr] = clone_command_parameter(length_param); /* keep what was l */
371    strcpy(l_par->name,"lrad"); /* but rename to lrad and slice : */
372    if (slices > 1) /* divide numbers or expressions by the number of slices */
373    {
374      if (l_par->expr) l_par->expr = compound_expr(l_par->expr,0.,"/",NULL,slices);
375      else l_par->double_value /= slices;
376    }
377    add_to_name_list("lrad",1,cmd->par_names);
378    cmd->par->curr++;
379  }
380}
381
382/* creates the thin magnetic element - recursively for classes from which dericed (parent) */
383static struct element*
384create_thin_multipole(struct element* thick_elem, int slice_no)
385{
386  struct command_parameter *angle_param, *length_param, *kparam[4], *ksparam[4], *kn_param, *ks_param, *at_param, *fint_param;
387  struct element *thin_elem_parent, *thin_elem;
388  struct command* cmd;
389  char *thin_name;
390  int angle_conversion = 0;
391  int slices, minimizefl;
392  int knl_flag = 0,ksl_flag = 0;
393 
394  /* next is new to handle parent with possibly different slice number than child */
395  slices = get_slices_from_elem(thick_elem);
396  at_param = return_param("at",thick_elem);
397 
398  if (thick_elem == thick_elem->parent) return NULL; /* no further parent to consider */
399  else
400  {
401    thin_elem_parent = create_thin_multipole(thick_elem->parent,slice_no); /* slice also the parent */
402  }
403 
404  minimizefl=get_option("minimizeparents") && !at_param && thick_elem == thick_elem->parent;
405  if(minimizefl)
406  {
407    slice_no=slices=1; /* do not slice this one */
408  }
409  if(slice_no > slices && thick_elem!=thick_elem->parent ) /* check, but not for base classes */
410  {
411    slice_no=1;
412  }
413 
414  /* check to see if we've already done this one */
415  thin_elem = get_thin(thick_elem,slice_no);
416  if (thin_elem) return thin_elem;
417 
418  /* issue a warning in case of element parameter combinations not suitable for slicing */
419  fint_param   = return_param_recurse("fint",thick_elem);
420  if(fint_param)
421  {
422    printf("    *** warning %s is a thick %s with fringe fields. These will be lost in the translation to a multipole. Use dipedge.\n",
423           thick_elem->name,thick_elem->parent->name);
424  }
425 
426  length_param = return_param_recurse("l",thick_elem);
427  angle_param  = return_param_recurse("angle",thick_elem);
428  kparam[0]    = return_param_recurse("k0",thick_elem);
429  kparam[1]    = return_param_recurse("k1",thick_elem);
430  kparam[2]    = return_param_recurse("k2",thick_elem);
431  kparam[3]    = return_param_recurse("k3",thick_elem);
432  ksparam[0]   = return_param_recurse("k0s",thick_elem);
433  ksparam[1]   = return_param_recurse("k1s",thick_elem);
434  ksparam[2]   = return_param_recurse("k2s",thick_elem);
435  ksparam[3]   = return_param_recurse("k3s",thick_elem);
436  kn_param     = return_param_recurse("knl",thick_elem);
437  ks_param     = return_param_recurse("ksl",thick_elem);
438  if (kn_param) {kn_param = clone_command_parameter(kn_param); knl_flag++;}
439  if (ks_param) {ks_param = clone_command_parameter(ks_param); ksl_flag++;}
440 
441  /* translate k0,k1,k2,k3,angle */
442  if ((kparam[0] || kparam[1] || kparam[2] || kparam[3] || angle_param
443       || ksparam[0] || ksparam[1] || ksparam[2] || ksparam[3])
444      && (kn_param==NULL && ks_param==NULL))
445  {
446    kn_param = new_command_parameter("knl", 12);
447    kn_param->expr_list = new_expr_list(10);
448    kn_param->double_array = new_double_array(10);
449    ks_param = new_command_parameter("ksl", 12);
450    ks_param->expr_list = new_expr_list(10);
451    ks_param->double_array = new_double_array(10);
452    angle_conversion = translate_k(kparam,ksparam,angle_param,kn_param,ks_param);
453  }
454 
455  kn_param = scale_and_slice(kn_param,length_param,slices,slice_no,
456                             angle_conversion,knl_flag+ksl_flag);
457  ks_param = scale_and_slice(ks_param,length_param,slices,slice_no,
458                             angle_conversion,knl_flag+ksl_flag);
459  /* set up new multipole command */
460  cmd = new_command(buffer("thin_multipole"), 20, 20, /* max num names, max num param */
461                    buffer("element"), buffer("none"), 0, 8); /* 0 is link, multipole is 8 */
462  add_cmd_parameter_new(cmd,1.,"magnet",0); /* parameter magnet with value of 1 and inf=0 */
463  if(!minimizefl)
464  {
465    add_cmd_parameter_clone(cmd,return_param("at"  ,thick_elem),"at"  ,1);
466    add_cmd_parameter_clone(cmd,return_param("from",thick_elem),"from",1);
467    add_lrad(cmd,length_param,slices);
468    add_cmd_parameter_clone(cmd,kn_param,"knl",1);
469    add_cmd_parameter_clone(cmd,ks_param,"ksl",1);
470  }
471  add_cmd_parameter_clone(cmd,return_param_recurse("apertype",thick_elem),"apertype",1);
472  add_cmd_parameter_clone(cmd,return_param_recurse("aperture",thick_elem),"aperture",1);
473  add_cmd_parameter_clone(cmd,return_param("bv",thick_elem),"bv",1);
474  add_cmd_parameter_clone(cmd,return_param_recurse("tilt",thick_elem),"tilt",1);
475  add_cmd_parameter_clone(cmd,return_param_recurse("kmax",    thick_elem),"kmax",    1);
476  add_cmd_parameter_clone(cmd,return_param_recurse("kmin",    thick_elem),"kmin",    1);
477  add_cmd_parameter_clone(cmd,return_param_recurse("calib",   thick_elem),"calib",   1);
478  add_cmd_parameter_clone(cmd,return_param_recurse("polarity",thick_elem),"polarity",1);
479  add_cmd_parameter_clone(cmd,return_param_recurse("mech_sep",thick_elem),"mech_sep",1);
480  /*== jln 11.11.2010 dealt with the new property v_pos as for mech_sep */
481  add_cmd_parameter_clone(cmd,return_param_recurse("v_pos",   thick_elem),"v_pos",1);
482  /*==*/
483  /* create element with this command */
484  if (slices==1 && slice_no==1) thin_name=buffer(thick_elem->name);
485  else
486  {
487    thin_name = make_thin_name(thick_elem->name,slice_no);
488  }
489 
490  if (thin_elem_parent)
491  {
492    thin_elem = make_element(thin_name,thin_elem_parent->name,cmd,-1);
493  }
494  else
495  {
496    thin_elem = make_element(thin_name,"multipole",cmd,-1);
497  }
498  thin_elem->length = 0;
499  thin_elem->bv = el_par_value("bv",thin_elem);
500  if (thin_elem_parent && thin_elem_parent->bv)
501  {
502    thin_elem->bv = thin_elem_parent->bv;
503  }
504  put_thin(thick_elem,thin_elem,slice_no);
505  return thin_elem;
506}
507
508static struct element*
509create_thin_solenoid(struct element* thick_elem, int slice_no)
510/*hbu create thin solenoid element, similar to create_thin_multipole */
511{
512  struct command_parameter *length_param, *ks_param, *at_param ,*ks_par;
513  struct element *thin_elem_parent, *thin_elem;
514  struct command* cmd;
515  char *thin_name;
516  int slices,minimizefl;
517 
518  if (thick_elem == thick_elem->parent) return NULL;
519  else
520  {
521    thin_elem_parent = create_thin_solenoid(thick_elem->parent,slice_no); /*hbu move up to parent */
522  }
523  thin_elem = get_thin(thick_elem,slice_no);
524  if (thin_elem) return thin_elem; /* is already thin */
525  slices = get_slices_from_elem(thick_elem);
526  /* get parameters from the thick solenoid element */
527  length_param  = return_param_recurse("l",thick_elem);
528  ks_param      = return_param_recurse("ks",thick_elem);
529  at_param      = return_param("at",thick_elem);
530 
531  minimizefl=get_option("minimizeparents") && !at_param && thick_elem == thick_elem->parent;
532  if(minimizefl)
533  {
534    slice_no=slices=1; /* do not slice this one */
535  }
536 
537  /* set up new solenoid command */
538  cmd = new_command(buffer("thin_solenoid"), 20, 20, /* max num names, max num param */
539                    buffer("element"), buffer("none"), 0, 9); /* 0 is link, solenoid is 9 */
540  add_cmd_parameter_new(cmd,1.,"magnet",0); /* parameter magnet with value of 1 and inf=0 */
541 
542 
543  if(!minimizefl)
544  {
545    add_cmd_parameter_clone(cmd,return_param("at"  ,thick_elem),"at"  ,1);
546    add_cmd_parameter_clone(cmd,return_param("from",thick_elem),"from",1);
547    add_lrad(cmd,length_param,slices);
548  }
549  add_cmd_parameter_clone(cmd,ks_param,"ks",1); /* keep ks */
550  if(!minimizefl)
551  {
552    if (length_param && ks_param) /* in addition provide   ksi = ks * l /slices */
553    {
554      ks_par = cmd->par->parameters[cmd->par->curr] = clone_command_parameter(ks_param); /* start from clone of ks */
555      strcpy(ks_par->name,"ksi"); /* change name to ksi */
556      if (length_param->expr || ks_par->expr) /* first step is ks * l calculation, expression or value */
557      {
558        ks_par->expr = compound_expr(ks_par->expr,ks_par->double_value,"*",length_param->expr,length_param->double_value); /* multiply expression with length */
559      }
560      else ks_par->double_value *= length_param->double_value; /* multiply value with length */
561      if (slices > 1) /* 2nd step, divide by slices, expression or number */
562      {
563        if (ks_par->expr) ks_par->expr = compound_expr(ks_par->expr,0.,"/",NULL,slices);
564        else ks_par->double_value /= slices;
565      }
566      add_to_name_list("ksi",1,cmd->par_names);
567      cmd->par->curr++;
568    }
569  }
570  add_cmd_parameter_clone(cmd,return_param_recurse("apertype",thick_elem),"apertype",1);
571  add_cmd_parameter_clone(cmd,return_param_recurse("aperture",thick_elem),"aperture",1);
572  add_cmd_parameter_clone(cmd,return_param("bv",thick_elem),"bv",1);
573  add_cmd_parameter_clone(cmd,return_param("tilt",thick_elem),"tilt",1);
574  /* create element with this command */
575  if (slices==1 && slice_no==1) thin_name=buffer(thick_elem->name);
576  else thin_name = make_thin_name(thick_elem->name,slice_no);
577  if (thin_elem_parent)
578  {
579    thin_elem = make_element(thin_name,thin_elem_parent->name,cmd,-1);
580  }
581  else
582  {
583    thin_elem = make_element(thin_name,"solenoid",cmd,-1);
584  }
585  thin_elem->length = 0;
586  thin_elem->bv = el_par_value("bv",thin_elem);
587  if (thin_elem_parent && thin_elem_parent->bv)
588  {
589    thin_elem->bv = thin_elem_parent->bv;
590  }
591  put_thin(thick_elem,thin_elem,slice_no);
592  return thin_elem;
593}
594
595static struct element*
596create_thin_elseparator(struct element* thick_elem, int slice_no)
597/*hbu create thin elseparator element, similar to and made from create_thin_solenoid */
598{
599  struct command_parameter *length_param, *ex_param, *ey_param, *tilt_param, *at_param ,*ex_par, *ey_par;
600  struct element *thin_elem_parent, *thin_elem;
601  struct command* cmd;
602  char *thin_name;
603  int slices,minimizefl;
604 
605  if (thick_elem == thick_elem->parent) return NULL;
606  else
607  {
608    thin_elem_parent = create_thin_elseparator(thick_elem->parent,slice_no); /*hbu move up to parent */
609  }
610  thin_elem = get_thin(thick_elem,slice_no);
611  if (thin_elem) return thin_elem; /* is already thin */
612  slices = get_slices_from_elem(thick_elem);
613  /* get parameters from the thick elseparator element */
614  length_param  = return_param_recurse("l",thick_elem);
615  ex_param      = return_param_recurse("ex",thick_elem);
616  ey_param      = return_param_recurse("ey",thick_elem);
617  tilt_param    = return_param_recurse("tilt",thick_elem);
618  at_param      = return_param("at",thick_elem);
619 
620  minimizefl=get_option("minimizeparents") && !at_param && thick_elem == thick_elem->parent;
621  if(minimizefl)
622  {
623    slice_no=slices=1; /* do not slice this one */
624  }
625 
626  /* set up new solenoid command */
627  cmd = new_command(buffer("thin_elseparator"), 20, 20, /* max num names, max num param */
628                    buffer("element"), buffer("none"), 0, 11); /* 0 is link, elseparator is 11 */
629  add_cmd_parameter_new(cmd,1.,"magnet",0); /* parameter magnet with value of 1 and inf=0 */
630 
631 
632  if(!minimizefl)
633  {
634    add_cmd_parameter_clone(cmd,return_param("at"  ,thick_elem),"at"  ,1);
635    add_cmd_parameter_clone(cmd,return_param("from",thick_elem),"from",1);
636    add_lrad(cmd,length_param,slices);
637  }
638  add_cmd_parameter_clone(cmd,ex_param,"ex",1); /* keep ex */
639  add_cmd_parameter_clone(cmd,ey_param,"ey",1); /* keep ey */
640  add_cmd_parameter_clone(cmd,tilt_param,"tilt",1); /* keep tilt */
641  if(!minimizefl)
642  {
643    /* create ex_l from ex */
644    if (length_param && ex_param) /* in addition provide   ex_l = ex * l /slices */
645    {
646      ex_par = cmd->par->parameters[cmd->par->curr] = clone_command_parameter(ex_param); /* start from clone of ex */
647      strcpy(ex_par->name,"ex_l"); /* change name to ex_l */
648      if (length_param->expr && ex_par->expr) /* first step is ex * l calculation, expression or value */
649      {
650        ex_par->expr = compound_expr(ex_par->expr,ex_par->double_value,"*",length_param->expr,length_param->double_value); /* multiply expression with length */
651      }
652      else ex_par->double_value *= length_param->double_value; /* multiply value with length */
653      if (slices > 1) /* 2nd step, divide by slices, expression or number */
654      {
655        if (ex_par->expr) ex_par->expr = compound_expr(ex_par->expr,0.,"/",NULL,slices);
656        else ex_par->double_value /= slices;
657      }
658      add_to_name_list("ex_l",1,cmd->par_names);
659      cmd->par->curr++;
660    }
661    /* create ey_l from ey */
662    if (length_param && ey_param) /* in addition provide   ey_l = ey * l /slices */
663    {
664      ey_par = cmd->par->parameters[cmd->par->curr] = clone_command_parameter(ey_param); /* start from clone of ey */
665      strcpy(ey_par->name,"ey_l"); /* change name to ey_l */
666      if (length_param->expr && ey_par->expr) /* first step is ey * l calculation, expression or value */
667      {
668        ey_par->expr = compound_expr(ey_par->expr,ey_par->double_value,"*",length_param->expr,length_param->double_value); /* multiply expression with length */
669      }
670      else ey_par->double_value *= length_param->double_value; /* multiply value with length */
671      if (slices > 1) /* 2nd step, divide by slices, expression or number */
672      {
673        if (ey_par->expr) ey_par->expr = compound_expr(ey_par->expr,0.,"/",NULL,slices);
674        else ey_par->double_value /= slices;
675      }
676      add_to_name_list("ey_l",1,cmd->par_names);
677      cmd->par->curr++;
678    }
679  }
680  add_cmd_parameter_clone(cmd,return_param_recurse("apertype",thick_elem),"apertype",1);
681  add_cmd_parameter_clone(cmd,return_param_recurse("aperture",thick_elem),"aperture",1);
682  add_cmd_parameter_clone(cmd,return_param("bv",thick_elem),"bv",1);
683  add_cmd_parameter_clone(cmd,return_param("tilt",thick_elem),"tilt",1);
684  /* create element with this command */
685  if (slices==1 && slice_no==1) thin_name=buffer(thick_elem->name);
686  else thin_name = make_thin_name(thick_elem->name,slice_no);
687  if (thin_elem_parent)
688  {
689    thin_elem = make_element(thin_name,thin_elem_parent->name,cmd,-1);
690  }
691  else
692  {
693    thin_elem = make_element(thin_name,"elseparator",cmd,-1);
694  }
695  thin_elem->length = 0;
696  thin_elem->bv = el_par_value("bv",thin_elem);
697  if (thin_elem_parent && thin_elem_parent->bv)
698  {
699    thin_elem->bv = thin_elem_parent->bv;
700  }
701  put_thin(thick_elem,thin_elem,slice_no);
702  return thin_elem;
703}
704
705/* put in one of those nice marker kind of things */
706static struct node*
707new_marker(struct node *thick_node, double at, struct expression *at_expr)
708{
709  struct node* node=NULL;
710  struct element* elem=NULL;
711 
712  int pos;
713  struct command* p;
714  struct command* clone;
715 
716  if (thick_node->p_elem)
717  {
718    pos = name_list_pos("marker", defined_commands->list);
719    /* clone = clone_command(defined_commands->commands[pos]); */
720    p = defined_commands->commands[pos];
721    clone = new_command(p->name, 20, 20, p->module, p->group, p->link_type,p->mad8_type);
722    /* 20 is the maximum number of par names and parvalues */
723    add_cmd_parameter_clone(clone,return_param_recurse("at",      thick_node->p_elem),"at",      1);
724    add_cmd_parameter_clone(clone,return_param_recurse("from",    thick_node->p_elem),"from",    1);
725    add_cmd_parameter_clone(clone,return_param_recurse("apertype",thick_node->p_elem),"apertype",1);
726    add_cmd_parameter_clone(clone,return_param_recurse("aperture",thick_node->p_elem),"aperture",1);
727    add_cmd_parameter_clone(clone,return_param_recurse("aper_tol",thick_node->p_elem),"aper_tol",1);
728    add_cmd_parameter_clone(clone,return_param_recurse("aper_offset",thick_node->p_elem),"aper_offset",1);
729    add_cmd_parameter_clone(clone,return_param_recurse("kmax",    thick_node->p_elem),"kmax",    1);
730    add_cmd_parameter_clone(clone,return_param_recurse("kmin",    thick_node->p_elem),"kmin",    1);
731    add_cmd_parameter_clone(clone,return_param_recurse("calib",   thick_node->p_elem),"calib",   1);
732    add_cmd_parameter_clone(clone,return_param_recurse("polarity",thick_node->p_elem),"polarity",1);
733    add_cmd_parameter_clone(clone,return_param_recurse("mech_sep",thick_node->p_elem),"mech_sep",1);
734    /*== dealt with the new property v_pos as for mech_sep */
735    add_cmd_parameter_clone(clone,return_param_recurse("v_pos",thick_node->p_elem),"v_pos",1);
736    /*==*/
737    elem = make_element(thick_node->p_elem->name, "marker", clone,-1);
738    node = new_elem_node(elem, thick_node->occ_cnt);
739    strcpy(node->name, thick_node->name);
740    node->occ_cnt = thick_node->occ_cnt;
741    node->at_value = at;
742    if (at_expr) node->at_expr = clone_expression(at_expr);
743    node->from_name = thick_node->from_name;
744  }
745  else
746  {
747    fatal_error("Oh dear, this is not an element!",thick_node->name);
748  }
749 
750  return node;
751}
752
753/* adds a thin elem in sliced nodes to the end of a sequence */
754static void
755seq_diet_add_elem(struct node* node, struct sequence* to_sequ)
756{
757  struct command_parameter *at_param, *length_param;
758  struct expression *l_expr = NULL, *at_expr = NULL;
759  struct node* thin_node;
760  struct element* elem;
761  double length = 0, at = 0;
762  int i,middle=-1,slices = 1;
763  char* old_thin_style;
764 
765  old_thin_style = NULL;
766  if (strstr(node->base_name,"collimator"))
767  {
768    elem = create_thin_obj(node->p_elem,1);
769    old_thin_style = thin_style;
770    thin_style = collim_style;
771  }
772  else if (strstr(node->base_name,"solenoid"))
773  {
774    elem = create_thin_solenoid(node->p_elem,1);  /* create the first thin solenoid slice */
775  }
776  else if (strstr(node->base_name,"elseparator"))
777  {
778    elem = create_thin_elseparator(node->p_elem,1);  /* create the first thin elseparator slice */
779  }
780  else
781  {
782    elem = create_thin_multipole(node->p_elem,1); /* get info from first slice */
783  }
784  slices = get_slices_from_elem(node->p_elem); /*hbu June 2005 */
785 
786  at_param = return_param_recurse("at",elem);
787  length_param = return_param_recurse("l",node->p_elem); /*get original length*/
788  if (length_param) l_expr  = length_param->expr;
789  if (at_param)     at_expr = at_param->expr;
790 
791  at     = el_par_value_recurse("at", elem);
792  length = el_par_value_recurse("l",node->p_elem);
793 
794  if (node->at_expr) at_expr = node->at_expr;
795  if (node->at_value != zero) at = node->at_value;
796  if (node->length   != zero) length = node->length;
797  /* note that a properly created clone node will contain the length of the element */
798  /* this will override all other definitions and hence the already sliced element length
799   is irrelevant */
800 
801  if (slices>1)
802  { /* sets after which element I should put the marker */
803    middle = abs(slices/2);
804  }
805 
806  for (i=0; i<slices; i++)
807  {
808    if (strstr(node->base_name,"collimator"))
809    {
810      elem = create_thin_obj(node->p_elem,i+1);
811    }
812    else if (strstr(node->base_name,"solenoid"))
813    {
814      elem = create_thin_solenoid(node->p_elem,i+1);
815    }
816    else if (strstr(node->base_name,"elseparator"))
817    {
818      elem = create_thin_elseparator(node->p_elem,i+1);
819    }
820    else
821    {
822      elem = create_thin_multipole(node->p_elem,i+1);
823    }
824    thin_node = new_elem_node(elem, node->occ_cnt);
825    thin_node->length   = 0.0;
826    thin_node->from_name = buffer(node->from_name);
827    if (fabs(at_shift(slices,i+1))>0.0)
828    {
829      if (at_expr || l_expr)
830      {
831        thin_node->at_expr =
832        compound_expr(at_expr,at,"+",scale_expr(l_expr,at_shift(slices,i+1)),
833                      length*at_shift(slices,i+1));
834      }
835    }
836    else
837    {
838      if (at_expr) thin_node->at_expr = clone_expression(at_expr);
839    }
840    thin_node->at_value = at + length*at_shift(slices,i+1);
841    if (i==middle) seq_diet_add(new_marker(node,at,at_expr),to_sequ);
842    seq_diet_add(thin_node,to_sequ);
843  }
844  if (strstr(node->base_name,"collimator")) thin_style=old_thin_style;
845  return;
846}
847
848/* creates the thin non-magnetic element - recursively */
849static struct element*
850create_thin_obj(struct element* thick_elem, int slice_no)
851{
852  struct element *thin_elem_parent = NULL , *thin_elem = NULL;
853  struct command* cmd = NULL;
854  struct command_parameter*  length_param= NULL;
855  int length_i = -1,lrad_i = -1,slices=1;
856  char* thin_name = NULL;
857 
858  if (thick_elem == thick_elem->parent)
859  {
860    return NULL;
861  }
862  else
863  {
864    thin_elem_parent = create_thin_obj(thick_elem->parent,slice_no);
865  }
866 
867  /* check to see if we've already done this one */
868  thin_elem = get_thin(thick_elem,slice_no);
869  if (thin_elem) return thin_elem;
870 
871  /* set up new multipole command */
872  cmd = clone_command(thick_elem->def);
873  length_param = return_param_recurse("l",thick_elem);
874  length_i = name_list_pos("l",thick_elem->def->par_names);
875  lrad_i   = name_list_pos("lrad",thick_elem->def->par_names);
876  if (length_param)
877  {
878    if (lrad_i > -1 && thick_elem->def->par_names->inform[lrad_i]>0)
879    {
880      /* already exists so replace lrad */
881      cmd->par->parameters[lrad_i]->double_value = cmd->par->parameters[length_i]->double_value;
882      if (cmd->par->parameters[length_i]->expr)
883      {
884        if (cmd->par->parameters[lrad_i]->expr)
885          delete_expression(cmd->par->parameters[lrad_i]->expr);
886        cmd->par->parameters[lrad_i]->expr =
887        clone_expression(cmd->par->parameters[length_i]->expr);
888      }
889    }
890    else
891    { /* doesn't exist */
892      if (name_list_pos("lrad",thick_elem->base_type->def->par_names)>-1)
893      {
894        /* add lrad only if allowed by element */
895        if (cmd->par->curr == cmd->par->max) grow_command_parameter_list(cmd->par);
896        if (cmd->par_names->curr == cmd->par_names->max)
897          grow_name_list(cmd->par_names);
898        cmd->par->parameters[cmd->par->curr] = clone_command_parameter(length_param);
899        add_to_name_list("lrad",1,cmd->par_names);
900        cmd->par->parameters[name_list_pos("lrad",cmd->par_names)]->expr =
901        clone_expression(cmd->par->parameters[length_i]->expr);
902        cmd->par->curr++;
903      }
904    }
905  }
906 
907  if (length_i > -1)
908  {
909    cmd->par->parameters[length_i]->double_value = 0;
910    cmd->par->parameters[length_i]->expr = NULL;
911  }
912  if (strstr(thick_elem->base_type->name,"collimator"))
913  {
914    slices = get_slices_from_elem(thick_elem);
915  }
916  if (slices==1 && slice_no==1)
917  {
918    thin_name=buffer(thick_elem->name);
919  }
920  else
921  {
922    thin_name=make_thin_name(thick_elem->name,slice_no);
923  }
924 
925  if (thin_elem_parent)
926  {
927    thin_elem = make_element(thin_name,thin_elem_parent->name,cmd,-1);
928  }
929  else
930  {
931    thin_elem = make_element(thin_name,thick_elem->base_type->name,cmd,-1);
932  }
933  thin_elem->length = 0;
934  thin_elem->bv = el_par_value("bv",thin_elem);
935 
936  put_thin(thick_elem,thin_elem,slice_no);
937  return thin_elem;
938}
939
940/* this copies an element node and sets the length to zero
941 and radiation length to the length
942 to be used for "copying" optically neutral elements */
943static struct node*
944copy_thin(struct node* thick_node)
945{
946  struct node* thin_node = NULL;
947 
948  thin_node = clone_node(thick_node, 0);
949  thin_node->length=0;
950  thin_node->p_elem->length=0;
951  /* if we have a non zero length then an lrad has to be created */
952  if (el_par_value("l",thick_node->p_elem)>zero)
953    thin_node->p_elem = create_thin_obj(thick_node->p_elem,1);
954 
955  return thin_node;
956}
957
958/* this decides how to split an individual node and
959 sends it onto the thin_sequ builder */
960static void
961seq_diet_node(struct node* thick_node, struct sequence* thin_sequ)
962{
963  struct node* thin_node;
964  if (thick_node->p_elem)
965  { /* this is an element to split and add */
966    if (el_par_value("l",thick_node->p_elem)==zero) /* if it's already thin copy it directly*/
967    {
968      seq_diet_add(thin_node = copy_thin(thick_node),thin_sequ);
969    }
970    else if(strcmp(thick_node->base_name,"matrix") == 0)
971    { /*hbu. Take matrix as it is, including any length */
972      seq_diet_add(thick_node,thin_sequ);
973    }
974    else
975    { /* we have to slim it down a bit...*/
976      if (strcmp(thick_node->base_name,"marker") == 0      ||
977          strcmp(thick_node->base_name,"instrument") == 0  ||
978          strcmp(thick_node->base_name,"placeholder") == 0 ||
979          strcmp(thick_node->base_name,"hmonitor") == 0    ||
980          strcmp(thick_node->base_name,"vmonitor") == 0    ||
981          strcmp(thick_node->base_name,"monitor") == 0     ||
982          strcmp(thick_node->base_name,"vkicker") == 0     ||
983          strcmp(thick_node->base_name,"hkicker") == 0     ||
984          strcmp(thick_node->base_name,"kicker") == 0      ||
985          strcmp(thick_node->base_name,"tkicker") == 0     ||
986          strcmp(thick_node->base_name,"rfcavity") == 0    ||
987          strcmp(thick_node->base_name,"crabcavity") == 0
988          )
989      {
990        seq_diet_add(thin_node = copy_thin(thick_node),thin_sequ);
991        /* special cavity list stuff */
992        if (strcmp(thin_node->p_elem->base_type->name, "rfcavity") == 0 &&
993            find_element(thin_node->p_elem->name, thin_sequ->cavities) == NULL)
994          add_to_el_list(&thin_node->p_elem, 0, thin_sequ->cavities, 0);
995        /* special crab cavity list stuff */
996        if (strcmp(thin_node->p_elem->base_type->name, "crabcavity") == 0 &&
997            find_element(thin_node->p_elem->name, thin_sequ->crabcavities) == NULL)
998          add_to_el_list(&thin_node->p_elem, 0, thin_sequ->crabcavities, 0);
999      }
1000      else if (strcmp(thick_node->base_name,"rbend") == 0       ||
1001               strcmp(thick_node->base_name,"sbend") == 0       ||
1002               strcmp(thick_node->base_name,"quadrupole") == 0  ||
1003               strcmp(thick_node->base_name,"sextupole") == 0   ||
1004               strcmp(thick_node->base_name,"octupole") == 0    ||
1005               strcmp(thick_node->base_name,"solenoid") == 0    || /*hbu */
1006               strcmp(thick_node->base_name,"multipole") == 0   || /* special spliting required. */
1007               strcmp(thick_node->base_name,"rcollimator") == 0 ||
1008               strcmp(thick_node->base_name,"ecollimator") == 0 ||
1009               strcmp(thick_node->base_name,"elseparator") == 0
1010               )
1011      {
1012        seq_diet_add_elem(thick_node,thin_sequ);
1013      }
1014      else if (strcmp(thick_node->base_name,"drift") == 0)
1015      {
1016        /* ignore this as it makes no sense to slice */
1017      }
1018      else /* new elements not yet implemented for slicing - write a message and do a reasonable default action */
1019      {
1020        fprintf(prt_file, "Found unknown basename %s, doing copy with length set to zero.\n",thick_node->base_name);
1021        seq_diet_add(copy_thin(thick_node),thin_sequ);
1022      }
1023    }
1024  }
1025  else if (thick_node->p_sequ)
1026  { /* this is a sequence to split and add */
1027    seq_diet_add_sequ(thick_node,seq_diet(thick_node->p_sequ),thin_sequ);
1028  }
1029  else
1030  { /* we have no idea what this is - serious error */
1031    fatal_error("node is not element or sequence",thick_node->base_name);
1032  }
1033}
1034
1035/* slim down this sequence - this is the bit to be called recursively */
1036/* this actually creates the thin sequence */
1037static struct sequence*
1038seq_diet(struct sequence* thick_sequ)
1039{
1040  struct node *thick_node = NULL;
1041  struct sequence* thin_sequ;
1042  char name[128];
1043  int pos;
1044 
1045  /* first check to see if it had been already sliced */
1046  if ((thin_sequ=get_thin_sequ(thick_sequ))) return thin_sequ;
1047 
1048  strcpy(name,thick_sequ->name);
1049  fprintf(prt_file, "makethin: slicing sequence : %s\n",name);
1050  thin_sequ = new_sequence(name, thick_sequ->ref_flag);
1051  thin_sequ->start = NULL;
1052  thin_sequ->share = thick_sequ->share;
1053  thin_sequ->nested = thick_sequ->nested;
1054  thin_sequ->length = sequence_length(thick_sequ);
1055  thin_sequ->refpos = buffer(thick_sequ->refpos);
1056  thin_sequ->ref_flag = thick_sequ->ref_flag;
1057  thin_sequ->beam = thick_sequ->beam;
1058  if (thin_sequ->cavities != NULL)  thin_sequ->cavities->curr = 0;
1059  else thin_sequ->cavities = new_el_list(100);
1060  if (thin_sequ->crabcavities != NULL)  thin_sequ->crabcavities->curr = 0;
1061  else thin_sequ->crabcavities = new_el_list(100);
1062  thick_node = thick_sequ->start;
1063  while(thick_node != NULL)
1064  { /* loop over current sequence */
1065    /* the nodes are added to the sequence in seq_diet_add() */
1066    seq_diet_node(thick_node,thin_sequ);
1067    if (thick_node == thick_sequ->end)  break;
1068    thick_node = thick_node->next;
1069  }
1070  thin_sequ->end->next = thin_sequ->start;
1071  /* now we have to move the pointer in the sequences list
1072   to point to our thin sequence */
1073  if ((pos = name_list_pos(name, sequences->list)) < 0)
1074  {
1075    fatal_error("unknown sequence sliced:", name);
1076  }
1077  else
1078  {
1079    sequences->sequs[pos]= thin_sequ;
1080    /* delete_sequence(thick_sequ) */
1081  }
1082 
1083  /* add to list of sequences sliced */
1084  put_thin_sequ(thick_sequ,thin_sequ);
1085 
1086  return thin_sequ;
1087}
1088
1089/*************************************************************************/
1090/* these are the routines to determine the method of splitting */
1091/* note slice number is counted from 1 NOT 0 */
1092
1093/* return at relative shifts from center of unsliced magnet */
1094static double
1095simple_at_shift(int slices, int slice_no)
1096{
1097  int n = slices;
1098  int i = slice_no;
1099
1100  return n>1 ? (2.0*i-1)/(2.0*n)-0.5 : 0.0;
1101}
1102
1103static double
1104teapot_at_shift(int slices, int slice_no)
1105{
1106  int n = slices;
1107  int i = slice_no;
1108
1109  // Formula taken from HBU presentation at LCU, 2012.09.18
1110  return n>1 ? 0.5*n*(1-2*i+n)/(1.0-n*n) : 0.0;
1111}
1112
1113static double
1114collim_at_shift(int slices,int slice_no)
1115{
1116  int n = slices;
1117  int i = slice_no;
1118
1119  return n>1 ? (i-1.0)/(n-1.0)-0.5 : 0.0;
1120}
1121
1122static double
1123default_at_shift(int slices, int slice_no)
1124{
1125  return slices>4 ? simple_at_shift(slices, slice_no)
1126                  : teapot_at_shift(slices, slice_no);
1127}
1128
1129/* previous Teapot limited to 4 slices.
1130static double
1131teapot_at_shift(int slices,int slice_no)
1132{
1133  double at = 0;
1134  switch (slices)
1135  {
1136    case 1:
1137      at = 0.;
1138      break;
1139    case 2:
1140      if (slice_no == 1) at = -1./3.;
1141      if (slice_no == 2) at = +1./3.;
1142      break;
1143    case 3:
1144      if (slice_no == 1) at = -3./8.;
1145      if (slice_no == 2) at = 0.;
1146      if (slice_no == 3) at = +3./8.;
1147      break;
1148    case 4:
1149      if (slice_no == 1) at = -2./5.;
1150      if (slice_no == 2) at = -2./15.;
1151      if (slice_no == 3) at = +2./15.;
1152      if (slice_no == 4) at = +2./5.;
1153      break;
1154  }
1155  // return the simple style if slices > 4
1156  if (slices > 4) at = simple_at_shift(slices,slice_no);
1157  return at;
1158}
1159*/
1160
1161/* return at relative strength shifts from unsliced magnet */
1162static double
1163simple_q_shift(int slices,int slice_no)
1164{
1165  (void)slice_no;
1166  return 1.0/slices;
1167}
1168
1169static double
1170teapot_q_shift(int slices,int slice_no)
1171{
1172  (void)slice_no;
1173  return 1.0/slices;
1174}
1175
1176static double
1177collim_q_shift(int slices,int slice_no)
1178{ /* pointless actually, but it pleases symmetrically */
1179  (void)slice_no;
1180  return 1.0/slices;
1181}
1182
1183static double
1184default_q_shift(int slices, int slice_no)
1185{
1186  return slices>4 ? simple_q_shift(slices, slice_no)
1187                  : teapot_q_shift(slices, slice_no);
1188}
1189
1190/* return at relative shifts from center of unsliced magnet */
1191static double
1192at_shift(int slices,int slice_no)
1193{
1194  if (!slices || !slice_no) {
1195    fatal_error("makethin: invalid slicing for zero slices",thin_style);
1196  }
1197
1198  if (thin_style == NULL) {
1199    return default_at_shift(slices,slice_no);
1200  }
1201  else if (strcmp(thin_style,"simple")==0) {
1202    return simple_at_shift(slices,slice_no);
1203  }
1204  else if (strcmp(thin_style,"teapot")==0) {
1205    return teapot_at_shift(slices,slice_no);
1206  }
1207  else if (strcmp(thin_style,"collim")==0) {
1208    return collim_at_shift(slices,slice_no);
1209  }
1210  else
1211  {
1212    fatal_error("makethin: Style chosen not known:",thin_style);
1213  }
1214  return 0;
1215}
1216
1217/* return at relative strength shifts from unsliced magnet */
1218static double
1219q_shift(int slices,int slice_no)
1220{
1221  if (!slices || !slice_no) {
1222    fatal_error("makethin: invalid slicing for zero slices",thin_style);
1223  }
1224
1225  if (thin_style == NULL) {
1226    return default_q_shift(slices,slice_no);
1227  }
1228  else if (strcmp(thin_style,"simple")==0) {
1229    return simple_q_shift(slices,slice_no);
1230  }
1231  else if (strcmp(thin_style,"teapot")==0) {
1232    return teapot_q_shift(slices,slice_no);
1233  }
1234  else if (strcmp(thin_style,"collim")==0) {
1235    return collim_q_shift(slices,slice_no);
1236  }
1237  else {
1238    fatal_error("makethin: Style chosen not known:",thin_style);
1239  }
1240  return 0;
1241}
1242
1243static void
1244set_selected_elements(void)
1245{ /* New set_selected_elements */
1246  struct element* el_j;
1247  struct name_list* nl;
1248  struct command_parameter_list* pl;
1249  int i, j, pos_slice, pos_full, pos_range, slice, el_j_slice_pos;
1250  bool full_fl,range_fl,slice_fl;
1251  struct node* c_node;    /* for range check.  current node */
1252  struct node* nodes[2];  /* for range check.  first and last in range */
1253  /* Init curr and list->curr in global el_list structure.  selected_elements is passed to add_to_el_list and used at the end as thin_select_list
1254   selected_elements  is only used in makethin (set here and read in and could be named thin_select_list
1255   */
1256  selected_elements->curr = 0;
1257  selected_elements->list->curr = 0;  /* Reset list->curr in global el_list structure.   selected_elements is passed to add_to_el_list */
1258  if (current_sequ == NULL || current_sequ->ex_start == NULL) /* check that there is an active sequence, otherwise crash in get_ex_range */
1259  {
1260    warning("makethin selection without active sequence,", "ignored");
1261    return;
1262  }
1263  /* default is full sequence from start to end */
1264  nodes[0] = current_sequ->ex_start;
1265  nodes[1] = current_sequ->ex_end;
1266  for (i = 0; i < slice_select->curr; i++) /* loop over "select,flag=makethin" commands */
1267  {
1268    nl = slice_select->commands[i]->par_names;
1269    pl = slice_select->commands[i]->par;
1270    pos_full  = name_list_pos("full", nl);
1271    full_fl   = pos_full  > -1 && nl->inform[pos_full];  /* selection with full */
1272    pos_range = name_list_pos("range", nl);
1273    range_fl  = pos_range > -1 && nl->inform[pos_range]; /* selection with range */
1274    pos_slice = name_list_pos("slice", nl);              /* position of slice parameter in select command list */
1275    slice_fl  = pos_slice > -1 && nl->inform[pos_slice]; /* selection with slice */
1276    if (slice_fl) slice = pl->parameters[pos_slice]->double_value; /* Parameter has been read. Slice number from select command */
1277    else slice = 1;
1278    if(full_fl) /* use full sequence from start to end, the default */
1279    {
1280      nodes[0] = current_sequ->ex_start;
1281      nodes[1] = current_sequ->ex_end;
1282    }
1283    if(range_fl)
1284    {
1285      if (current_sequ == NULL || current_sequ->ex_start == NULL) /* check that there is an active sequence, otherwise crash in get_ex_range */
1286      {
1287        warning("makethin range selection without active sequence,", "ignored");
1288        return;
1289      }
1290      if( get_ex_range(pl->parameters[pos_range]->string, current_sequ, nodes) == 0) /* set start nodes[0] and end notes[1] depending on the range string */
1291      {
1292        printf("    +++ warning, empty range");
1293        continue;
1294      }
1295    }
1296    if(slice_fl) /* Set slice number in elements. Add to list of selected_elements */
1297    {
1298      if(range_fl) /* now elements in the sequence in the range */
1299      {
1300        c_node = nodes[0];
1301        while (c_node != NULL) /* loop over nodes in range,  set slice number in elements */
1302        {
1303          el_j = c_node->p_elem;
1304          el_j_slice_pos = name_list_pos("slice",el_j->def->par_names); /* position of slice parameter in element list */
1305          if (pass_select(el_j->name, slice_select->commands[i]) != 0) /* selection on class and pattern done in pass_select. element el_j selected */
1306          { /* the element el_j passes the selection */
1307            if(el_j_slice_pos > 0) el_j->def->par->parameters[el_j_slice_pos]->double_value=slice; /* Set the element slice number to the number of slices given in the select statement. */
1308            if( name_list_pos(el_j->name, selected_elements->list) < 0) /* el_j not yet in selected_elements */
1309            {
1310              add_to_el_list(&el_j, slice, selected_elements, 0);
1311            } /* new selection */
1312          } /* selection */
1313          if (c_node == nodes[1]) break; /* done with last node */
1314          c_node = c_node->next;
1315        } /* end of while loop over nodes in range */
1316      } /* range_fl */
1317      else /* no range_fl */
1318      {
1319        for(j=0; j< element_list->curr; j++) /* loop over element_list */
1320        {
1321          el_j = element_list->elem[j];
1322          el_j_slice_pos = name_list_pos("slice",el_j->def->par_names);
1323          if (pass_select(el_j->name, slice_select->commands[i]) != 0) /* selection on class and pattern done in pass_select. element el_j selected */
1324          { /* the element el_j passes the selection */
1325            if(el_j_slice_pos > 0) el_j->def->par->parameters[el_j_slice_pos]->double_value=slice; /* Set the element slice number to the number of slices given in the select statement. */
1326            if( name_list_pos(el_j->name, selected_elements->list) < 0) /* el_j not yet in selected_elements */
1327            {
1328              add_to_el_list(&el_j, slice, selected_elements, 0);
1329            } /* new selection */
1330          } /* selection */
1331        } /* loop over element_list */
1332      } /* range_fl */
1333    } /* slice_fl */
1334  } /* end of loop over select slice commands */
1335}
1336
1337/*************************************************************************/
1338
1339// public interface
1340
1341/* This converts the MAD-X command to something I can use
1342 if a file has been specified we send the command to exec_save
1343 which writes the file for us */
1344void
1345makethin(struct in_cmd* cmd)
1346{
1347  struct sequence *thick_sequ = NULL ,*thin_sequ = NULL;
1348  struct name_list* nl = cmd->clone->par_names;
1349  struct command_parameter_list* pl = cmd->clone->par;
1350  char *name = NULL;
1351  int pos,pos2;
1352  int k=0;
1353  /*    time_t start; */
1354 
1355  /*    start = time(NULL); */
1356  pos = name_list_pos("style", nl);
1357  if (nl->inform[pos] && (name = pl->parameters[pos]->string))
1358  {
1359    thin_style = buffer(pl->parameters[pos]->string);
1360    fprintf(prt_file, "makethin: style chosen : %s\n",thin_style);
1361  }
1362 
1363  /* first check makethin parameters which influence the selection */
1364 
1365  pos = name_list_pos("minimizeparents", nl);
1366  /* k = true; */  /* Use this to set minimizeparents to true by default. */
1367  if( pos > -1 && nl->inform[pos])
1368  {
1369    k=pl->parameters[pos]->double_value;
1370  }
1371  set_option("minimizeparents", &k);
1372 
1373  pos = name_list_pos("makeconsistent", nl);
1374  if( pos > -1 && nl->inform[pos])
1375  {
1376    k=pl->parameters[pos]->double_value;
1377    set_option("makeconsistent", &k);
1378  }
1379 
1380  if (slice_select->curr > 0)
1381  {
1382    set_selected_elements(); /* makethin selection */
1383    thin_select_list = selected_elements;
1384  }
1385  if (thin_select_list == NULL)
1386  {
1387    warning("makethin: no selection list,","slicing all to one thin lens.");
1388  }
1389  else if (thin_select_list->curr == 0)
1390  {
1391    warning("makethin selection list empty,","slicing all to one thin lens.");
1392  }
1393  if(get_option("makeconsistent"))
1394  {
1395    force_consistent_slices();
1396  }
1397  pos = name_list_pos("sequence", nl);
1398  if (nl->inform[pos] && (name = pl->parameters[pos]->string))
1399  {
1400    if ((pos2 = name_list_pos(name, sequences->list)) >= 0)
1401    {
1402      thick_sequ = sequences->sequs[pos2];
1403      thin_sequ = seq_diet(thick_sequ);
1404      disable_line(thin_sequ->name, line_list);
1405    }
1406    else warning("unknown sequence ignored:", name);
1407  }
1408  else warning("makethin without sequence:", "ignored");
1409 
1410  /* fprintf(prt_file, "makethin: finished in %f seconds.\n",difftime(time(NULL),start)); */
1411  thin_select_list = NULL;
1412}
1413
Note: See TracBrowser for help on using the repository browser.