source: PSPA/madxPSPA/src/mad_eval.c @ 441

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

import madx-5.01.00
File size: 21.1 KB
Line 
1#include "madx.h"
2
3static double
4act_value(int pos, struct name_list* chunks)
5  /* returns the actual value of a variable, element, or command parameter */
6{
7  char* name = chunks->names[pos];
8  char comm[NAME_L];
9  char par[NAME_L];
10  char *p, *n = name, *q = comm;
11  double val = zero;
12  struct element* el;
13  struct command* cmd = NULL;
14  if ((p = strstr(name, "->")) == NULL) /* variable */
15  {
16    if ((current_variable = find_variable(name, variable_list)) == NULL)
17    {
18      if (get_option("verify"))
19        warning("undefined variable set to zero:", name);
20      current_variable = new_variable(name, zero, 1, 1, NULL, NULL);
21      val = zero;
22      add_to_var_list(current_variable, variable_list, 0);
23    }
24    else val = variable_value(current_variable);
25  }
26  else /* element or command parameter */
27  {
28    while (n < p)  *(q++) = *(n++);
29    *q = '\0';
30    q = par; n++; n++;
31    while (*n != '\0')  *(q++) = *(n++);
32    *q = '\0';
33    if (strncmp(comm, "beam", 4) == 0)
34    {
35      cmd = current_beam = find_command("default_beam", beam_list);
36      if ((p = strchr(comm, '%')) != NULL)
37      {
38        if ((current_beam = find_command(++p, beam_list)) == NULL)
39          current_beam = cmd;
40      }
41      val = command_par_value(par, current_beam);
42    }
43    else if ((el = find_element(comm, element_list)) != NULL)
44      val = el_par_value(par, el);
45    else if ((cmd = find_command(comm, stored_commands)) != NULL)
46      val = command_par_value(par, cmd);
47    else if ((cmd = find_command(comm, beta0_list)) != NULL)
48      val = command_par_value(par, cmd);
49    else if ((cmd = find_command(comm, defined_commands)) != NULL)
50      val = command_par_value(par, cmd);
51  }
52  return val;
53}
54
55static int
56simple_logic_expr(int nit, char* toks[])
57  /* evaluates a logical expression of type expr1 op expr2
58     where op is one of ==, <>, <, >, <=, >=
59     returns -1 if invalid
60     0 if false
61     1 if true */
62{
63  int i, t1, t2, l1_start, l1_end, l2_start, l2_end,
64    logex = -1, brack = 0;
65  double val1, val2;
66  char c;
67
68  for (i = 0; i < nit; i++)
69  {
70    c = *toks[i];
71    switch (c)
72    {
73      case '(':
74        brack++;
75        break;
76      case ')':
77        brack--;
78        break;
79      case '<':
80      case '>':
81      case '=':
82        if (brack == 0)  goto found;
83    }
84  }
85  return -1;
86  found:
87  l1_start = 0;
88  l2_start = (*toks[i+1] == '=' || *toks[i+1] == '>') ? i + 2 : i + 1;
89  if ((t1 = loc_expr(toks, i, l1_start, &l1_end)) == 0) return -1;
90  if (t1 == 2)
91  {
92    if (polish_expr(l1_end + 1 - l1_start, &toks[l1_start]) == 0)
93      val1 = polish_value(deco, join(&toks[l1_start], l1_end + 1 - l1_start));
94    else return -1;
95  }
96  else  val1 = simple_double(toks, l1_start, l1_end);
97  if ((t2 = loc_expr(toks, nit, l2_start, &l2_end)) == 0) return -1;
98  if (t2 == 2)
99  {
100    if (polish_expr(l2_end + 1 - l2_start, &toks[l2_start]) == 0)
101      val2 = polish_value(deco, join(&toks[l2_start], l2_end + 1 - l2_start));
102    else return -1;
103  }
104  else  val2 = simple_double(toks, l2_start, l2_end);
105  if (c == '<')
106  {
107    if (*toks[i+1] == '=')     logex = val1 <= val2 ? 1 : 0;
108    else if(*toks[i+1] == '>') logex = val1 != val2 ? 1 : 0;
109    else                       logex = val1 <  val2 ? 1 : 0;
110  }
111  else if (c == '>')
112  {
113    if (*toks[i+1] == '=')  logex = val1 >= val2 ? 1 : 0;
114    else                    logex = val1 >  val2 ? 1 : 0;
115  }
116  else if (c == '=' && *toks[i+1] == '=')  logex = val1 == val2 ? 1 : 0;
117  return logex;
118}
119
120static int
121logic_expr(int nit, char* toks[])
122{
123  /* evaluates composite logical expression of type simple1 op simple2,
124     where simple1 and simple2 are simple logical expressions, and
125     op is either || or &&;
126     returns -1 for invalid
127     0 for false
128     1 for true */
129  int i = 0, k;
130  char c = ' ';
131
132  for (i = 0; i < nit; i++)
133  {
134    if (*toks[i] == '|' || *toks[i] == '&')
135    {
136      c = *toks[i]; break;
137    }
138  }
139  if (i == nit) return simple_logic_expr(nit, toks);
140  if (++i == nit || *toks[i] != c) return -1;
141  k = simple_logic_expr(i-1, toks);
142  if (c == '|')
143  {
144    if (k != 0) return 1;
145    else if (++i == nit)  return -1;
146    else return simple_logic_expr(nit-i, &toks[i]);
147  }
148  else
149  {
150    if (k == 0) return 0;
151    else if (++i == nit)  return -1;
152    else return simple_logic_expr(nit-i, &toks[i]);
153  }
154}
155
156// public interface
157
158void
159deco_init(void)
160  /* initializes Polish decoding */
161{
162  expr_chunks = new_name_list("expr_chunks", 2000);
163  cat = new_int_array(MAX_ITEM);
164  deco = new_int_array(MAX_ITEM);
165  d_var = new_int_array(MAX_ITEM);
166  oper = new_int_array(MAX_ITEM);
167  func = new_int_array(MAX_ITEM);
168  cat_doubles = new_double_array(MAX_ITEM);
169  doubles = new_double_array(MAX_D_ITEM);
170  twiss_deltas = new_double_array(MAX_ITEM);
171}
172
173int
174act_special(int type, char* statement)
175  /* acts on special commands (IF{..} etc.) */
176{
177  struct char_array* loc_buff = NULL;
178  struct char_array* loc_w = NULL;
179  int cnt_1, start_2, rs, re, level = pro->curr, ls = strlen(statement);
180  int ret_val = 0;
181  struct char_p_array* logic;
182  int logex = 0;
183  char *cp = statement;
184  if (ls < IN_BUFF_SIZE) ls = IN_BUFF_SIZE;
185  if (level == pro->max) grow_in_buff_list(pro);
186  if (pro->buffers[level] == NULL)
187    pro->buffers[level] = new_in_buffer(ls);
188  else
189  {
190    while(pro->buffers[level]->c_a->max < ls)
191      grow_char_array(pro->buffers[level]->c_a);
192  }
193  if (type == 5) /* macro */ return make_macro(statement);
194  else if (type == 6) /* line */ return make_line(statement);
195  logic = new_char_p_array(1000);
196  loc_buff = new_char_array(ls);
197  loc_w = new_char_array(ls);
198  get_bracket_range(statement, '{', '}', &rs, &re);
199  if (re < 0) fatal_error("missing '{' or '}' in statement:",statement);
200  cnt_1 = rs; start_2 = rs + 1;
201  mystrcpy(loc_buff, statement); loc_buff->c[re] =  '\0';
202  while(aux_buff->max < cnt_1) grow_char_array(aux_buff);
203  strncpy(aux_buff->c, statement, cnt_1); aux_buff->c[cnt_1] = '\0';
204  switch (type)
205  {
206    case 1:  /* if */
207      pro->buffers[level]->flag = 0;
208    case 3:  /* else if */
209      if (pro->buffers[level]->flag < 0)
210      {
211        ret_val = -1;
212        break;
213      }
214      if (pro->buffers[level]->flag == 0)
215      {
216        pre_split(aux_buff->c, loc_w, 0);
217        mysplit(loc_w->c, tmp_l_array);
218        get_bracket_t_range(tmp_l_array->p, '(', ')', 0, tmp_l_array->curr,
219                            &rs, &re);
220        rs++;
221        if ((logex = logic_expr(re-rs, &tmp_l_array->p[rs])) > 0)
222        {
223          pro->buffers[level]->flag = 1;
224          pro->curr++;
225          /* now loop over statements inside {...} */
226          pro_input(&loc_buff->c[start_2]);
227          pro->curr--;
228        }
229        else if (logex < 0) warning("illegal if construct set false:", cp);
230      }
231      break;
232    case 2: /* else */
233      if (pro->buffers[level]->flag < 0)
234      {
235        ret_val = -1;
236        break;
237      }
238      if (pro->buffers[level]->flag == 0)
239      {
240        pro->curr++;
241        /* now loop over statements inside {...} */
242        pro_input(&loc_buff->c[start_2]);
243        pro->curr--;
244        pro->buffers[level]->flag = -1;
245      }
246      break;
247    case 4: /* while */
248      pre_split(aux_buff->c, loc_w, 0);
249      mysplit(loc_w->c, logic);
250      get_bracket_t_range(logic->p, '(', ')', 0, logic->curr,
251                          &rs, &re);
252      pro->curr++; rs++;
253      while ((logex = logic_expr(re-rs, &logic->p[rs])) > 0)
254      {
255        /* now loop over statements inside {...} */
256        pro_input(&loc_buff->c[start_2]);
257      }
258      pro->curr--;
259      break;
260    default:
261      ret_val = -1;
262  }
263  if (loc_buff != NULL) delete_char_array(loc_buff);
264  if (loc_w != NULL) delete_char_array(loc_w);
265  delete_char_p_array(logic, 0);
266  return ret_val;
267}
268
269void
270pro_input(char* statement)
271  /* processes one special (IF() etc.), or one normal statement after input */
272{
273  int type, code, nnb, ktmp;
274  char* sem;
275  int rs, re, start = 0, l = strlen(statement);
276
277  clearerrorflag(); /*reset global error flag */
278
279  while (start < l)
280  {
281    if ((type = in_spec_list(&statement[start])))
282    {
283      if (type == 6)
284      {
285        get_bracket_range(&statement[start], '(', ')', &rs, &re);
286        ktmp = re+1;
287        if (re > rs && strchr(&statement[ktmp], ':')) /* formal arg.s */
288        {
289          get_bracket_range(&statement[ktmp], '(', ')', &rs, &re);
290          rs += ktmp; re += ktmp;
291        }
292      }
293      else get_bracket_range(&statement[start], '{', '}', &rs, &re);
294      if (re > rs)
295      {
296        re += start + 1;
297        if (re < l && next_non_blank(&statement[re]) == ';')
298          re += next_non_blank_pos(&statement[re]) + 1;
299      }
300      if((code = act_special(type, &statement[start])) < 0)
301      {
302        if (get_option("warn"))
303        {
304          switch (code)
305          {
306            case -1:
307              warning("statement illegal in this context,", "skipped");
308              break;
309            case -2:
310              warning("statement label is protected keyword,", "skipped");
311            case -3:
312              warning("statement not recognised:", statement);
313              break;
314            default:
315              fatal_error("illegal return code","from act_special");
316          }
317        }
318      }
319      if (re > rs && re < l)
320      {
321        if (next_non_blank_pos(&statement[re]) < 0) start = l;
322        else start = re;
323      }
324      else start = l;
325    }
326    else
327    {
328      if ((sem = mystrchr(&statement[start], ';')) == NULL) return;
329      if (sem > &statement[start]) /* skip empty ';' */
330      {
331        *sem = '\0';
332        this_cmd = new_in_cmd(400);
333        pre_split(&statement[start], work, 1);
334        check_table(work->c); check_tabstring(work->c);
335        this_cmd->tok_list->curr = mysplit(work->c, this_cmd->tok_list);
336        if ((type = decode_command()) < 0) /* error */
337        {
338          if (get_option("warn"))
339          {
340            switch (type)
341            {
342              case -1:
343                warning("statement illegal in this context,", "skipped");
344                break;
345              case -2:
346                warning("statement label is protected keyword,","skipped");
347              case -3:
348                warning("statement not recognised:", work->c);
349                break;
350              default:
351                fatal_error("illegal return code","from decode_command");
352            }
353          }
354        }
355        else process();
356        if (stop_flag)  return;
357        *sem = ';';
358      }
359      sem++;
360      start = sem - statement;
361      if (start < l)
362      {
363        if ((nnb = next_non_blank_pos(sem)) < 0)  start = l;
364        else start += nnb;
365      }
366    }
367  }
368}
369
370void
371process(void)  /* steering routine: processes one command */
372{
373  int pos;
374  char* name;
375  struct element* el;
376  if (this_cmd != NULL)
377  {
378
379    switch (this_cmd->type)
380    {
381      case 0: /* executable commands */
382        exec_command();
383        if (stop_flag)
384        {
385          if (this_cmd)
386          {
387            if (this_cmd->clone != NULL)
388              this_cmd->clone = delete_command(this_cmd->clone);
389            this_cmd = delete_in_cmd(this_cmd);
390          }
391          return;
392        }
393        break;
394      case 1: /* element definition */
395        enter_element(this_cmd);
396        buffer_in_cmd(this_cmd);
397        break;
398      case 2: /* variable definition */
399        enter_variable(this_cmd);
400        break;
401      case 3: /* sequence start or end */
402        enter_sequence(this_cmd);
403        break;
404      case 4:
405        name = this_cmd->tok_list->p[0];
406        if (sequ_is_on)
407          /* element or sequence reference in sequence */
408        {
409          if ((pos = name_list_pos(name, sequences->list)) < 0)
410            enter_element(this_cmd);
411          else
412          {
413            this_cmd->cmd_def = find_command("sequence", defined_commands);
414            this_cmd->clone = clone_command(this_cmd->cmd_def);
415            strcpy(this_cmd->clone->name, name);
416            scan_in_cmd(this_cmd);
417            enter_sequ_reference(this_cmd, sequences->sequs[pos]);
418          }
419        }
420        else
421          /* element parameter definition */
422        {
423          if ((el = find_element(name, element_list)) == NULL)
424            warning("skipped, command or element unknown:", name);
425          else
426          {
427            this_cmd->cmd_def = el->def;
428            this_cmd->clone = clone_command(this_cmd->cmd_def);
429            strcpy(this_cmd->clone->name, name);
430            scan_in_cmd(this_cmd);
431            update_element(el, this_cmd->clone);
432          }
433        }
434        break;
435      default:
436        warning("unknown command type:",
437                join_b(this_cmd->tok_list->p, this_cmd->tok_list->curr));
438    }
439    if (this_cmd != NULL && (this_cmd->type == 0 || this_cmd->type == 2))
440    {
441      if (this_cmd->clone != NULL)
442      {
443        if (this_cmd->clone_flag == 0)
444          this_cmd->clone = delete_command(this_cmd->clone);
445        else add_to_command_list(this_cmd->clone->name,
446                                 this_cmd->clone, stored_commands, 0);
447      }
448      if (this_cmd->label != NULL) buffer_in_cmd(this_cmd);
449      else this_cmd = delete_in_cmd(this_cmd);
450    }
451  }
452}
453
454int
455polish_expr(int c_item, char** item)   /* split input */
456  /* deco output array containing:
457     expression in Polish notation of length deco->curr,
458     coded as 0-, 1+, 2*, 3/, 4^ (power),
459     6 evaluate function
460     100000000 + n = variable n (refers to vars),
461     200000000 + n = function n (refers to functs),
462     400000000 + n = real n (refers to doubles)
463     -- Example: suppose a, b are variables 0 and 4, exp is function 3:
464     then     3 * a * q[l] * q[k1] / exp((b - 1.57)^2) + 1.57
465     would result in
466     400000000 100000000 2 100000001 2 100000002 2
467     100000003 400000001 0 400000002 3 200000003 3 400000001 1
468     where 3 = real 0, 1.57 = real 1, 2 = real 2
469     a = vars 0, q[l] vars 1, q[k1] vars 2, exp functs 3
470  */
471{
472  int i, j, error, op, id, stack = 0, l_deco, l_double;
473  int up[100][3] = {{-1, -1, -1}};
474
475  l_deco = deco->curr = 0;
476  l_double = doubles->curr;
477  error = scan_expr(c_item, item);
478  if (error) return error;
479  for (i = 0; i < cat->curr; i++)
480  {
481
482    /* categories: 1: variable, 3: floating constant, 4: operator
483       6: left par., 7: right par.     */
484    switch (cat->i[i])
485    {
486      case 1:                              /* variable */
487        if (l_deco == deco->max) grow_int_array(deco);
488        deco->i[l_deco++] = 100000000 + d_var->i[i];
489        break;
490      case 3:                              /* constant */
491        if (l_deco == deco->max) grow_int_array(deco);
492        if (l_double == doubles->max) grow_double_array(doubles);
493        doubles->a[l_double] = cat_doubles->a[i];
494        deco->i[l_deco++] = 400000000 + l_double++;
495        doubles->curr = l_double;
496        break;
497      case 4:
498        if ((op = oper->i[i]) < 5)           /* operator */
499        {
500          id = op / 2;
501          for (j = 2; j >= id; j--)
502          {
503            if (up[stack][j] > -1)
504            {
505              if (l_deco == deco->max) grow_int_array(deco);
506              deco->i[l_deco++] = up[stack][j];
507              up[stack][j] = -1;
508            }
509          }
510          up[stack][id] = op;
511        }
512        else
513        {
514          if (l_deco == deco->max) grow_int_array(deco);
515          deco->i[l_deco++] = 200000000 + func->i[i];  /* function */
516        }
517        break;
518      case 6:      /*  '(' */
519        stack++;
520        for (j = 0; j < 3; j++)  up[stack][j] = -1;
521        break;
522      case 7:      /*  ')' */
523        for (j = 2; j >= 0; j--)
524        {
525          if (up[stack][j] > -1)
526          {
527            if (l_deco == deco->max) grow_int_array(deco);
528            deco->i[l_deco++] = up[stack][j];
529          }
530        }
531        stack--;
532        break;
533      default:
534        return 9;
535    }   /* end switch */
536  }     /* end loop over categories */
537  for (j = 2; j >= 0; j--)   /* clear stack */
538  {
539    if (up[stack][j] > -1)
540    {
541      if (l_deco == deco->max) grow_int_array(deco);
542      deco->i[l_deco++] = up[stack][j];
543    }
544  }
545  deco->curr = l_deco;
546  return 0;
547}
548
549double
550polish_value(struct int_array* deco, char* expr_string)
551  /* coded input (see below) */
552  /* description see polish_expression */
553{
554  int i, k, kc, c_stack = -1;
555  double stack[MAX_ITEM];
556  char tmp[20];
557
558  if (++polish_cnt > MAX_LOOP)
559    fatal_error("circular call in expression", expr_string);
560  stack[0] = 0;
561  for (i = 0; i < deco->curr; i++)   /* decoding loop */
562  {
563    k = deco->i[i];
564    if ( k < 5)     /* operator */
565    {
566      if (c_stack < 0)
567      {
568        fatal_error("stack underflow in expression:", expr_string);
569      }
570      else if (c_stack == 0)
571      {
572        stack[1] = stack[0]; stack[0] = 0;
573      }
574      else  c_stack--;
575
576
577      switch(k)
578      {
579        case 0:
580          stack[c_stack] -= stack[c_stack+1];
581          break;
582        case 1:
583          stack[c_stack] += stack[c_stack+1];
584          break;
585        case 2:
586          stack[c_stack] *= stack[c_stack+1];
587          break;
588        case 3:
589          if (stack[c_stack+1] == 0.0)
590          {
591            warning("division by zero, result set to zero, expr:", expr_string);
592            stack[c_stack] = 0.0;
593            break;
594          }
595          stack[c_stack] /= stack[c_stack+1];
596          break;
597        case 4:
598          stack[c_stack] = pow(stack[c_stack],stack[c_stack+1]);
599          break;
600        default:
601          fatal_error("illegal operator, Polish, expr.:", expr_string);
602      }
603    }
604    else
605    {
606      kc = k / 100000000;  k -= 100000000 * kc;
607      switch(kc)
608      {
609        case 1:            /* variable */
610          stack[++c_stack] = act_value(k, expr_chunks);
611          break;
612        case 4:            /* real constant */
613          stack[++c_stack] = doubles->a[k];
614          break;
615        case 2:            /* function */
616          switch(k-1)      /* the offset is due to dummyfunction */
617          {
618            case 0:
619              stack[c_stack] = fabs(stack[c_stack]);
620              break;
621            case 1:
622              stack[c_stack] = sqrt(stack[c_stack]);
623              break;
624            case 2:
625              stack[c_stack] = exp(stack[c_stack]);
626              break;
627            case 3:
628              stack[c_stack] = log(stack[c_stack]);
629              break;
630            case 4:
631              stack[c_stack] = log10(stack[c_stack]);
632              break;
633            case 5:
634              stack[c_stack] = sin(stack[c_stack]);
635              break;
636            case 6:
637              stack[c_stack] = cos(stack[c_stack]);
638              break;
639            case 7:
640              stack[c_stack] = tan(stack[c_stack]);
641              break;
642            case 8:
643              stack[c_stack] = asin(stack[c_stack]);
644              break;
645            case 9:
646              stack[c_stack] = acos(stack[c_stack]);
647              break;
648            case 10:
649              stack[c_stack] = atan(stack[c_stack]);
650              break;
651            case 11:
652              stack[c_stack] = sinh(stack[c_stack]);
653              break;
654            case 12:
655              stack[c_stack] = cosh(stack[c_stack]);
656              break;
657            case 13:
658              stack[c_stack] = tanh(stack[c_stack]);
659              break;
660            case 14:
661              stack[c_stack] = frndm();
662              break;
663            case 15:
664              stack[c_stack] = grndm();
665              break;
666            case 16:
667              stack[c_stack] = tgrndm(stack[c_stack]);
668              break;
669            case 17:
670              stack[c_stack] = table_value();
671              break;
672            case 18: /* function "exist" */
673              continue; /* value in stack not changed */
674              break;
675            case 19:
676              stack[c_stack] = floor(stack[c_stack]);
677              break;
678            case 20:
679              stack[c_stack] = ceil(stack[c_stack]);
680              break;
681            case 21:
682              stack[c_stack] = rint(stack[c_stack]);
683              break;
684            case 22: {
685              double int_part;
686              stack[c_stack] = modf(stack[c_stack], &int_part);
687              } break;
688            default:
689              fatal_error("polish_value: illegal function in expr:",
690                          expr_string);
691          }
692          break;
693        default:
694          /* if you get here, then most likely someone has created
695             more than 100000000 double precision numbers */
696          sprintf(tmp, "%d", k-1);
697          fatal_error("illegal type in Polish decoding: ", tmp);
698          exit(1);
699      }
700    }
701  }       /* end of decoding loop */
702  polish_cnt--;
703  return stack[0];
704}
705
706void
707print_value(struct in_cmd* cmd)
708{
709  char** toks = &cmd->tok_list->p[cmd->decl_start];
710  int n = cmd->tok_list->curr - cmd->decl_start;
711  int j, s_start = 0, end, nitem; // , type // not used
712  while (s_start < n)
713  {
714    for (j = s_start; j < n; j++) if (*toks[j] == ',') break;
715    if (loc_expr(toks, j, s_start, &end) > 0) // type = // not used
716    {
717      nitem = end + 1 - s_start;
718      if (polish_expr(nitem, &toks[s_start]) == 0)
719        fprintf(prt_file, v_format("%s = %F ;\n"),
720                spec_join(&toks[s_start], nitem), 
721                polish_value(deco, join(&toks[s_start], nitem)));
722      else
723      {
724        warning("invalid expression:", spec_join(&toks[s_start], nitem));
725        return;
726      }
727      s_start = end+1;
728      if (s_start < n-1 && *toks[s_start] == ',') s_start++;
729    }
730    else
731    {
732      warning("invalid expression:", spec_join(&toks[s_start], n));
733      return;
734    }
735  }
736}
737
Note: See TracBrowser for help on using the repository browser.