source: PSPA/madxPSPA/src/mad_elemerr.c @ 430

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

import madx-5.01.00

File size: 26.1 KB
Line 
1#include "madx.h"
2
3static int
4find_index_in_table(char *cols[], const char *name )
5{
6  for (int i = 0; strcmp(cols[i], " "); i++)
7    if (string_icmp(cols[i], name) == 0)
8      return i;
9
10  return -1; // not found
11}
12
13static void
14pro_error_make_efield_table(void)
15{
16  struct table *ttb = efield_table;
17  struct node *nanf;
18  struct node *nend;
19  int         j;
20  struct sequence* mysequ = current_sequ;
21
22  setbuf(stdout,(char *)0);
23  nanf = mysequ->ex_start;
24  nend = mysequ->ex_end;
25
26      while (nanf != nend) {
27        if(nanf->sel_err == 1) {
28           string_to_table_curr("efield","name",nanf->name);
29/* */
30                  /*
31           printf("=> %s %e %e %e\n",nanf->name,nanf->p_fd_err,nanf->p_al_err);
32                  */
33           if(nanf->p_fd_err != NULL) {
34              int from_col = find_index_in_table(efield_table_cols, "k0l");
35              int to_col = find_index_in_table(efield_table_cols, "k20sl");
36              int ncols = to_col - from_col + 1;
37              for (j=0; j < ncols; j++) {
38                 ttb->d_cols[from_col+j][ttb->curr] = nanf->p_fd_err->a[j];
39                  /*
40                 printf("Field: %d %e\n",j,ttb->d_cols[j][ttb->curr]);
41                  */
42              }
43           }
44           if(nanf->p_al_err != NULL) {
45              int from_col = find_index_in_table(efield_table_cols, "dx");
46              int to_col = find_index_in_table(efield_table_cols, "mscaly");
47              int ncols = to_col - from_col + 1;
48              for (j=0; j < ncols; j++) {
49                 ttb->d_cols[from_col+j][ttb->curr] = nanf->p_al_err->a[j];
50                  /*
51                 printf("Align: %d %e\n",j,ttb->d_cols[j][ttb->curr]);
52                  */
53              }
54           }
55           /* AL: RF-errors */
56           if(nanf->p_ph_err != NULL) {
57              ttb->d_cols[find_index_in_table(efield_table_cols, "rfm_freq")][ttb->curr] = nanf->rfm_freq;
58              ttb->d_cols[find_index_in_table(efield_table_cols, "rfm_harmon")][ttb->curr] = nanf->rfm_harmon;
59              ttb->d_cols[find_index_in_table(efield_table_cols, "rfm_lag")][ttb->curr] = nanf->rfm_lag;
60              int from_col = find_index_in_table(efield_table_cols, "p0l");
61              int to_col = find_index_in_table(efield_table_cols, "p20sl");
62              int ncols = to_col - from_col + 1;
63              for (j=0; j < ncols; j++) {
64                 ttb->d_cols[from_col+j][ttb->curr] = nanf->p_ph_err->a[j];
65                  /*
66                 printf("Align: %d %e\n",j,ttb->d_cols[j][ttb->curr]);
67                  */
68              }
69           }
70/* */
71           augment_count("efield");
72        }
73        nanf = nanf->next;
74      }
75}
76
77static void
78error_seterr(struct in_cmd* cmd)
79{
80
81/* read the errors from a named table  and stores
82   them in the nodes of the sequence.       
83   Subsequent Twiss will use them correctly.
84   ===> Must be preceded by a call to "read_table"
85   ===> (unless table exists in memory !)
86*/
87
88  int from_col, to_col, row, col, i;
89
90  struct node *node, *node_end;
91
92  char     name[NAME_L];
93  char   slname[NAME_L];
94
95  char     nname[NAME_L];
96  char   slnname[NAME_L];
97
98  char    *namtab, namtab_buf[NAME_L];
99  int      t1;
100
101  struct   table *err;
102
103/* set up pointers to current sequence for later use */
104  struct sequence* mysequ = current_sequ;
105  node     = mysequ->ex_start;
106  node_end = mysequ->ex_end;
107
108  if ((namtab = command_par_string("table",cmd->clone)) != NULL) {
109    printf("Want to use named table: %s\n",namtab);
110    if ((t1 = name_list_pos(namtab, table_register->names)) > -1)
111      printf("The table ==> %s <=== was found \n",namtab);
112    else {
113      warning("No such error table in memory:", namtab);
114      exit(-77);
115    }
116  }
117  else {
118    if (get_option("debug")) {
119      printf("No table name requested\n");
120      printf("Use default name\n");
121    }
122
123    strcpy(namtab=namtab_buf,"error");
124    if ((t1 = name_list_pos(namtab, table_register->names)) > -1)
125      printf("The default table ==> %s <=== was found \n",namtab);
126    else {
127      warning("No default error table in memory:", namtab);
128      exit(-77);
129    }
130  }
131 
132  err = table_register->tables[t1];
133
134  for (row = 1; row <= err->curr; row++) {
135    if (string_from_table_row(namtab, "name", &row, name)) break;
136
137    // probably useless...
138    stolower(name);
139    strcpy(slname,strip(name));
140    supp_tb(slname);
141
142    for (node = mysequ->ex_start; node != node_end; node = node->next) {
143    // probably useless...
144      strcpy(nname,node->name);
145      stolower(nname);
146      strcpy(slnname,strip(nname));
147      supp_tb(slnname);
148
149      if(strcmp(slname, slnname) == 0) break;
150    }
151
152    /* We have now the input and the node, generate array and selection flag */
153    if (!strcmp(slname, slnname)) {
154      node->sel_err = 1;
155      node->p_fd_err = new_double_array(FIELD_MAX); // zero initialized
156      node->p_fd_err->curr = FIELD_MAX;
157      node->p_al_err = new_double_array(ALIGN_MAX); // zero initialized
158      node->p_al_err->curr = ALIGN_MAX;
159      node->p_ph_err = new_double_array(RFPHASE_MAX); // zero initialized
160      node->p_ph_err->curr = RFPHASE_MAX;
161
162      from_col = find_index_in_table(efield_table_cols, "k0l");
163      to_col   = find_index_in_table(efield_table_cols, "k20sl");
164      if (from_col > 0 && to_col > 0)
165        for (i=0, col=from_col; col <= to_col; col++, i++)
166          node->p_fd_err->a[i] = err->d_cols[col][row-1];
167
168      from_col = find_index_in_table(efield_table_cols, "dx");
169      to_col   = find_index_in_table(efield_table_cols, "mscaly");
170      if (from_col > 0 && to_col > 0)
171        for (i=0, col=from_col; col <= to_col; col++, i++)
172          node->p_al_err->a[i] = err->d_cols[col][row-1];
173
174      col = find_index_in_table(efield_table_cols, "rfm_freq");
175      node->rfm_freq   = col < 0 ? 0 : err->d_cols[col][row-1];
176
177      col = find_index_in_table(efield_table_cols, "rfm_harmon");
178      node->rfm_harmon = col < 0 ? 0 : err->d_cols[col][row-1];
179
180      col = find_index_in_table(efield_table_cols, "rfm_lag");
181      node->rfm_lag    = col < 0 ? 0 : err->d_cols[col][row-1];
182
183      from_col = find_index_in_table(efield_table_cols, "p0l");
184      to_col   = find_index_in_table(efield_table_cols, "p20sl");
185      if (from_col > 0 && to_col > 0)
186        for (i=0, col=from_col; col <= to_col; col++, i++)
187          node->p_ph_err->a[i] = err->d_cols[col][row-1];
188    }
189  }
190}
191
192static void
193error_esave(struct in_cmd* cmd)
194{
195    char *ef_table_file;
196/*  if(efield_table == NULL) { */
197       efield_table = make_table("efield", "efield", efield_table_cols,
198                               efield_table_types, 10000);
199       add_to_table_list(efield_table, table_register);
200       pro_error_make_efield_table();
201/*  }                          */
202    ef_table_file = command_par_string("file",cmd->clone);
203    out_table("efield",efield_table,ef_table_file);
204}
205
206static void
207error_ealign(struct in_cmd* cmd)
208{
209  struct node *ndexe;
210  struct node *nextnode;
211  int i;
212  int chcount[3] = {0,0,0};
213  double val[ALIGN_MAX] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0};
214  static char att[ALIGN_MAX][7] = {"dx","dy","ds","dphi","dtheta","dpsi","mrex","mrey","mredx","mredy","arex","arey","mscalx","mscaly"};
215
216  struct command_parameter_list* pl = current_error->par;
217  struct sequence* mysequ = current_sequ;
218
219  nextnode = mysequ->ex_start;
220  ndexe = mysequ->ex_end;
221
222  while (nextnode != ndexe) {
223
224      if(nextnode->sel_err == 1) {
225        if(nextnode->p_al_err == NULL) {
226          chcount[0]++;
227          nextnode->p_al_err = new_double_array(ALIGN_MAX);
228          nextnode->p_al_err->curr = ALIGN_MAX;
229        } else {
230          if(add_error_opt == 1) {
231              chcount[2]++;
232          } else {
233              chcount[1]++;
234          }
235        }
236          for(i=0;i<ALIGN_MAX;i++){
237               val[i] = pl->parameters[i]->double_value;
238               if(pl->parameters[i]->expr != NULL) {
239                  if(pl->parameters[i]->expr->status != 1) {
240                      val[i] = command_par_value(att[i], cmd->clone);
241                      pl->parameters[i]->expr->status = 0;
242                  }
243               }
244               if(add_error_opt == 1) {
245                  nextnode->p_al_err->a[i] += val[i];
246               } else {
247                  nextnode->p_al_err->a[i] = val[i];
248               }
249          }
250      }  /* end of treatment of selected node */
251      nextnode = nextnode->next;
252  }  /* end of loop over all nodes */
253  if(chcount[0] != 0)
254  fprintf(prt_file, "Assigned alignment errors to %d elements\n",chcount[0]);
255  if(chcount[1] != 0)
256  fprintf(prt_file, "Replaced alignment errors for %d elements\n",chcount[1]);
257  if(chcount[2] != 0)
258  fprintf(prt_file, "Added alignment errors to %d elements\n",chcount[2]);
259
260}
261
262static void
263error_eprint(struct in_cmd* cmd)
264{
265  struct node *ndexe;
266  struct node *nextnode;
267  static char pln_alig[ALIGN_MAX][7] = {"dx","dy","ds","dphi","dtheta","dpsi","mrex","mrey","mredx","mredy","arex","arey","mscalx","mscaly"};
268  static float alig_fact[ALIGN_MAX] = {1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1.0,1.0};
269
270  int i;
271  struct sequence* mysequ = current_sequ;
272  int mycount;
273  int fll;
274
275  ndexe = mysequ->ex_end;
276  nextnode = mysequ->ex_start;
277
278  mycount = 0;
279 
280  fll = command_par_value("full", cmd->clone);
281
282  while (nextnode != ndexe) {
283
284    if((nextnode->sel_err == 1) || (fll > 0) ){
285       if(nextnode->p_al_err != NULL) {
286         fprintf(prt_file, "\n\nAlignment errors for element %s \n",nextnode->name);
287         fprintf(prt_file,"\nDisplacements in [mm], rotations in [mrad] \n");
288         fprintf(prt_file," %6s %10s %12s %10s %10s %10s %10s %10s %10s %10s %10s %10s %10s %10s\n",
289                  pln_alig[0], pln_alig[1],pln_alig[2],pln_alig[3],
290                  pln_alig[4], pln_alig[5],pln_alig[6],pln_alig[7],
291                  pln_alig[8], pln_alig[9],pln_alig[10],pln_alig[11],
292                  pln_alig[12],pln_alig[13]);
293         for(i=0;i<nextnode->p_al_err->curr;i++) { 
294            fprintf(prt_file, "%10.6f ",alig_fact[i]*nextnode->p_al_err->a[i]);
295         }
296         fprintf(prt_file, "\n");
297       } else {
298         fprintf(prt_file, "\n\nAlignment errors for element %s \n",nextnode->name);
299         fprintf(prt_file,"\nDisplacements in [mm], rotations in [mrad] \n");
300         fprintf(prt_file," %6s %10s %12s %10s %10s %10s %10s %10s %10s %10s %10s %10s %10s %10s\n",
301                  pln_alig[0], pln_alig[1],pln_alig[2],pln_alig[3],
302                  pln_alig[4], pln_alig[5],pln_alig[6],pln_alig[7],
303                  pln_alig[8], pln_alig[9],pln_alig[10],pln_alig[11],
304                  pln_alig[12],pln_alig[13]);
305         fprintf(prt_file, "\n");
306       }
307
308       if(nextnode->p_fd_err != NULL) {
309         mycount++;
310         if(mycount <= 50000) {
311           /* fprintf(prt_file,"%s %d\n",nextnode->name,(int)nextnode->p_fd_err); */
312           fprintf(prt_file, "\n\nField errors for element %s \n",nextnode->name);
313           fprintf(prt_file, "Multipole order:     Normal:           Skew: \n");
314           for(i=0;i<EFIELD_TAB;i++) {
315              fprintf(prt_file, "%8d          %8e      %8e\n",i/2,
316                     nextnode->p_fd_err->a[i],
317                     nextnode->p_fd_err->a[i+1]);
318              i++;
319           }
320           fprintf(prt_file, "\n");
321         }
322       } else {
323         mycount++;
324         if(mycount <= 50000) {
325           fprintf(prt_file, "\n\nField errors for element %s \n",nextnode->name);
326           fprintf(prt_file, "Multipole order:     Normal:           Skew: \n");
327         }
328       }
329    }
330
331      nextnode = nextnode->next;
332  }
333
334}
335
336static void
337error_efcomp(struct in_cmd* cmd)
338{
339  //  struct name_list* nl;
340  //  int pos;
341 
342  //  struct node *ndexe;
343  //  struct node *nextnode;
344  //  int    lvec;
345  int    hyst = 0;
346  //  int    flgmgt = 0;
347  int chcount[3] = {0,0,0};
348  char rout_name[] = "error_efcomp";
349  //  double norfac; /* factor for normalization at reference radius */
350  int    n = -1;     /* order of reference multipole */
351  double rr = 0.0;   /* reference radius for multipole error */
352  double rrr = 0.0;  /* reference radius for multipole error */
353  double freq = 0.0; /* frequency for RF-Multiples */
354  int harmon = 0;    /* harmonic number for RF-Multipoles */
355  double lag = 0.0;  /* lag for RF-Multipoles */
356  //  struct double_array *ptr;
357  //  struct double_array *pcoef;
358  double h_co_n[FIELD_MAX/2][4];
359  double h_co_s[FIELD_MAX/2][4];
360  //  double *hco_n;
361  //  double *hco_s;
362  //  double *nvec;
363  //  double deer; 
364  //  double ref_str;
365  //  double ref_strn;
366  //  double ref_len;
367  //  double nlength;
368  //  double nvec0, nvec1, nvec2, nvec3;
369  //  double val[4] = {0, 0, 0, 0};
370  static const char *atts[6] = {"order","radius","hyster","rfm_freq", "rfm_harmon", "rfm_lag"};
371  static const char *attv[6] = {"dkn","dks","dknr","dksr","dpn","dps"};
372  const size_t attv_len = sizeof attv/sizeof *attv;
373  const size_t atts_len = sizeof atts/sizeof *atts;
374  int iattv[attv_len];
375  struct sequence* mysequ = current_sequ;
376 
377  // double *nvec = (double *)mycalloc("error_efcomp",1000, sizeof(double));
378
379  struct node *ndexe = mysequ->ex_end;
380  struct node *nextnode = mysequ->ex_start;
381 
382  struct name_list* nl = cmd->clone->par_names;
383
384  const int opt_debug = get_option("debug");
385 
386  /* here comes a kludge, check which of the assignment vectors is there */
387  /*
388    i = 0;
389    while((cmd->tok_list->p[i]) != NULL) {
390    for(k=0;k<4;k++) {
391    if(strcmp(cmd->tok_list->p[i],attv[k]) == 0) {
392    iattv[k] = 1;
393    }
394    }
395    i++;
396    }
397  */
398  for(unsigned int k=0; k<attv_len; k++) {
399    iattv[k] = 0;
400    int pos = name_list_pos(attv[k],nl);
401    if(nl->inform[pos] > 0) {
402      if (opt_debug)
403        fprintf(prt_file, "set iattv %d for %s to 1\n",iattv[k],attv[k]);
404      iattv[k] = 1;
405    }
406  }
407 
408  for(int i=0;i<FIELD_MAX/2;i++) {
409    for(int j=0;j<4;j++) {
410      h_co_n[i][j] = 0.0;
411      h_co_s[i][j] = 0.0;
412    }
413  }
414 
415  while (nextnode != ndexe) { /*loop over elements and get strengths in vector*/
416    current_node = nextnode;
417    int flgmgt = node_value("magnet");
418    if((nextnode->sel_err == 1) && (flgmgt == 1))  {
419      if(nextnode->p_fd_err == NULL) {
420        chcount[0]++;
421        nextnode->p_fd_err = new_double_array(FIELD_MAX);
422        nextnode->p_fd_err->curr = FIELD_MAX;
423      } else {
424        if(add_error_opt == 1) {
425          chcount[2]++;
426        } else {
427          chcount[1]++;
428        }
429      }
430     
431      if (opt_debug)
432        fprintf(prt_file, "field for %s %s %d\n",
433                nextnode->name,nextnode->base_name,nextnode->sel_err);
434     
435      /* now get order (n), radius (rr) and hyster flag (hyst) from command, if any */
436      /* AL: added 'freq' option for RF-Multipoles */
437      for(unsigned int i=0;i<atts_len;i++){
438        double val = command_par_value(atts[i],cmd->clone);
439        if(i==0) {
440          n = val;
441          /* debug printout */
442          if (opt_debug)
443            fprintf(prt_file, "order  is %d\n",n);
444        } else if (i==1) {
445          rrr = val;
446          rr  = fabs(rrr);
447          /* debug printout */
448          if (opt_debug)
449            fprintf(prt_file, "radius is %f\n",val);
450        } else if (i==2) {
451          hyst = val;
452          /* debug printout */
453          if (opt_debug)
454            fprintf(prt_file, "hyster flag is %d\n",(int)val);
455        } else if (i==3) {
456          freq = val;
457          nextnode->rfm_freq = freq;
458          /* debug printout */
459          if (opt_debug)
460            fprintf(prt_file, "freq flag is %d\n",(int)val);
461        } else if (i==4) {
462          harmon = (int)val;
463          nextnode->rfm_harmon = harmon;
464          /* debug printout */
465          if (opt_debug)
466            fprintf(prt_file, "harmon flag is %d\n",(int)val);
467        } else if (i==5) {
468          lag = lag;
469          nextnode->rfm_lag = lag;
470          /* debug printout */
471          if (opt_debug)
472            fprintf(prt_file, "lag flag is %d\n",(int)val);
473        }
474      }
475     
476      /* now get coefficients for time memory effects in magnets */
477      {
478        struct double_array *pcoef = command_par_array("hcoeffn",cmd->clone);
479        if(pcoef != NULL) {
480          double *hco_n = &h_co_n[0][0];
481          for(int j=0;j<pcoef->curr;j++) {
482            *hco_n = pcoef->a[j];
483            hco_n++;
484          }
485          if (opt_debug) {
486            for(int j=0;j<FIELD_MAX/2;j++) {
487              printf("COEFF: %d %e %e %e %e\n",j,h_co_n[j][0],h_co_n[j][1],h_co_n[j][2],h_co_n[j][3]);
488            }
489          }
490        }
491      }
492      {
493        struct double_array *pcoef = command_par_array("hcoeffs",cmd->clone);
494        if(pcoef != NULL) {
495          double *hco_s = &h_co_s[0][0];
496          for(int j=0;j<pcoef->curr;j++) {
497            *hco_s = pcoef->a[j];
498            hco_s++;
499          }
500          if (opt_debug) {
501            for(int j=0;j<FIELD_MAX/2;j++) {
502              printf("COEFF: %d %e %e %e %e\n",j,h_co_s[j][0],h_co_s[j][1],h_co_s[j][2],h_co_s[j][3]);
503            }
504          }
505        }
506      }
507      /* get length of node and check if magnet */
508      double ref_str  = 0.0;
509      double ref_strn = 0.0;
510      double nlength = node_value("l");
511      // double ref_len = nlength; // never used
512      if (opt_debug)
513        fprintf(prt_file, "original length is %f\n",nlength);
514     
515      if(strcmp(nextnode->base_name,"multipole") == 0) {
516        double *nvec;
517        if ((nvec = (double *)mycalloc("error_efcomp",1000, sizeof(double)))) {
518          int lvec;
519          if(rrr > 0 ) {
520            get_node_vector("knl",&lvec,nvec);
521          } else {
522            get_node_vector("ksl",&lvec,nvec);
523          }
524          if (opt_debug) {
525            for(int i=0;i<4;i++) {
526              fprintf(prt_file, "original field = %d is %f\n",i,nvec[i]);
527            }
528          }
529          if (opt_debug)
530            fprintf(prt_file, "====n====>>> %d %f %f \n\n",n,nvec[n],nlength);
531          ref_str = nvec[n];
532          ref_strn = fabs(ref_str);
533          myfree(rout_name,nvec);
534        }
535      } else if (strcmp(nextnode->base_name,"sbend") == 0) {
536        double nvec0 = node_value("k0");
537        if (opt_debug) {
538          fprintf(prt_file, "original field0 is %f\n",nvec0);
539          fprintf(prt_file, "====0====>>> %d %f %f \n\n",n,nvec0,nlength);
540        }
541        ref_str = nvec0*nlength;
542        ref_strn = fabs(nvec0);
543      } else if (strcmp(nextnode->base_name,"rbend") == 0) {
544        double nvec0 = node_value("k0");
545        if (opt_debug) {
546          fprintf(prt_file, "original field0 is %f\n",nvec0);
547          fprintf(prt_file, "====0====>>> %d %f %f \n\n",n,nvec0,nlength);
548        }
549        ref_str = nvec0*nlength;
550        ref_strn = fabs(nvec0);
551      } else if (strcmp(nextnode->base_name,"quadrupole") == 0) {
552        double nvec1 = node_value("k1");
553        if (opt_debug) {
554          fprintf(prt_file, "original field1 is %f\n",nvec1);
555          fprintf(prt_file, "====1====>>> %d %f %f \n\n",n,nvec1,nlength);
556        }
557        ref_str = nvec1*nlength;
558        ref_strn = fabs(nvec1);
559      } else if (strcmp(nextnode->base_name,"sextupole") == 0) {
560        double nvec2 = node_value("k2");
561        if (opt_debug) {
562          fprintf(prt_file, "original field2 is %f\n",nvec2);
563          fprintf(prt_file, "====2====>>> %d %f %f \n\n",n,nvec2,nlength);
564        }
565        ref_str = nvec2*nlength;
566        ref_strn = fabs(nvec2);
567      } else if (strcmp(nextnode->base_name,"octupole") == 0) {
568        double nvec3 = node_value("k3");
569        if (opt_debug) {
570          fprintf(prt_file, "original field3 is %f\n",nvec3);
571          fprintf(prt_file, "====3====>>> %d %f %f \n\n",n,nvec3,nlength);
572        }
573        ref_str = nvec3*nlength;
574        ref_strn = fabs(nvec3);
575      }
576     
577      /*  edbug print out field components , not done for production version
578       */
579     
580      /* normal components -> 2j, skew components 2j+1 */
581      /* AL: the following 'if' is not necessary */
582      if(flgmgt == 1) {
583       
584        for(unsigned i=0;i<attv_len;i++)  {   /* loop over possible commands */
585         
586          if (opt_debug) fprintf(prt_file, "%s %d\n",attv[i],iattv[i]);
587         
588          struct double_array *ptr = command_par_array(attv[i],cmd->clone);
589          if((ptr != NULL) && (iattv[i] == 1))  { /* command [i] found ? */
590           
591            for(int j=0;j<ptr->curr;j++)  { /* loop over all parameters */
592             
593              /* start field error assignment */
594              /* NORMAL COMPONENTS, ABSOLUTE ERRORS */
595              if(i==0)  {
596                if(add_error_opt == 1) {
597                  nextnode->p_fd_err->a[2*j]   += ptr->a[j];
598                } else {
599                  nextnode->p_fd_err->a[2*j]    = ptr->a[j];
600                }
601               
602                /* SKEW COMPONENTS, ABSOLUTE ERRORS */
603              } else if(i==1) {
604                if(add_error_opt == 1) {
605                  nextnode->p_fd_err->a[2*j+1] += ptr->a[j];
606                } else {
607                  nextnode->p_fd_err->a[2*j+1]  = ptr->a[j];
608                }
609               
610                /* NORMAL COMPONENTS, RELATIVE ERRORS, MAY BE CORRECTED FOR MEMORY EFFECTS */
611              } else if(i==2) {
612                if(fabs(rr) < 1.0E-6) {
613                  printf("++++++ error: trying to assign relative field errors \n");
614                  printf("       with no or zero reference radius specified\n");
615                  exit(-1);
616                }
617                double norfac = pow(rr,(n-j)) * (fact(j)/fact(n));
618               
619                /* if flag for hysteresis correction is set, use coefficients for correction */
620                double deer = 0.0;
621                if(hyst == 1) {
622                  deer = h_co_n[j][3]*pow(ref_strn,3) + h_co_n[j][2]*pow(ref_strn,2) + 
623                    h_co_n[j][1]*pow(ref_strn,1) + h_co_n[j][0];
624                  if (opt_debug) 
625                    printf("after correction (n): %d %e %e %e %e\n",
626                           j,ref_strn,ptr->a[j],deer,(ptr->a[j] + deer));
627                }
628                /*
629                  if (opt_debug)
630                  fprintf(prt_file, "norm(n): %d %d %f %f\n",n,j,rr,norfac);
631                */
632                if(add_error_opt == 1) {
633                  nextnode->p_fd_err->a[2*j]   += (ptr->a[j]+deer)*ref_str*norfac;
634                } else {
635                  nextnode->p_fd_err->a[2*j]    = (ptr->a[j]+deer)*ref_str*norfac;
636                }
637               
638                /* SKEW COMPONENTS, RELATIVE ERRORS, MAY BE CORRECTED FOR MEMORY EFFECTS */
639              } else if(i==3) {
640                if(fabs(rr) < 1.0E-6) {
641                  printf("++++++ error: trying to assign relative field errors \n");
642                  printf("       with no or zero reference radius specified\n");
643                  exit(-1);
644                }
645                double norfac = pow(rr,(n-j)) * (fact(j)/fact(n));
646               
647                /* if flag for hysteresis correction is set, use coefficients for correction */
648                double deer = 0.0;
649                if(hyst == 1) {
650                  deer = h_co_s[j][3]*pow(ref_strn,3) + h_co_s[j][2]*pow(ref_strn,2) + 
651                    h_co_s[j][1]*pow(ref_strn,1) + h_co_s[j][0];
652                  if (opt_debug) 
653                    printf("after correction (s): %d %e %e %e %e\n",
654                           j,ref_strn,ptr->a[j],deer,(ptr->a[j] + deer));
655                }
656                /*
657                  if (opt_debug)
658                  fprintf(prt_file, "norm(s): %d %d %f %f\n",n,j,rr,norfac);
659                */
660                if(add_error_opt == 1) {
661                  nextnode->p_fd_err->a[2*j+1] += (ptr->a[j]+deer)*ref_str*norfac;
662                } else {
663                  nextnode->p_fd_err->a[2*j+1]  = (ptr->a[j]+deer)*ref_str*norfac;
664                }
665
666                /* RF-PHASE OF NORMAL COMPONENTS */
667              } else if(i==4) {
668                nextnode->p_ph_err->a[2*j] = ptr->a[j];     
669
670                /* RF-PHASE OF SKEW COMPONENTS */
671              } else if(i==5) {
672                nextnode->p_ph_err->a[2*j+1] = ptr->a[j];     
673
674              } /*  end  of field error assignment */ 
675            }
676          }
677        }
678      }
679    }  /* end of treatment of selected node */
680    nextnode = nextnode->next;
681  }      /* end of loop over all nodes */
682  if(chcount[0] != 0)
683    fprintf(prt_file, "Assigned field errors to %d elements\n",chcount[0]);
684  if(chcount[1] != 0)
685    fprintf(prt_file, "Replaced field errors for %d elements\n",chcount[1]);
686  if(chcount[2] != 0)
687    fprintf(prt_file, "Added field errors to %d elements\n",chcount[2]);
688  //  myfree(rout_name,nvec);
689}
690
691static void
692error_efield(struct in_cmd* cmd)
693{
694  int i;
695  /*
696  struct node *ndexs, *ndexe;
697  struct node *nextnode;
698  struct name_list* nl = current_error->par_names;
699  struct command_parameter_list* pl = current_error->par;
700  struct sequence* mysequ = current_sequ;
701  */
702
703  fprintf(prt_file, "in efield routine\n");
704  fprintf(prt_file, "efield command not yet implemented\n");
705
706  if (get_option("debug"))
707  {
708   for(i=0;i<cmd->tok_list->curr;i++)
709     fprintf(prt_file, "command(s): %s\n",cmd->tok_list->p[i]);
710  }
711}
712
713static void
714error_eoption(struct in_cmd* cmd)
715{
716  struct name_list* nl = cmd->clone->par_names;
717  int i, debug;
718  int val, pos, seed;
719  int is, ia;
720  static  int  ia_seen = 0;
721
722  is = 0; ia = 0;
723 
724  i = 0;
725  while(cmd->tok_list->p[i] != NULL) {
726     if(strcmp("add",cmd->tok_list->p[i]) == 0) {
727          ia = 1;
728     }
729     if(strcmp("seed",cmd->tok_list->p[i]) == 0) {
730          is = 1;
731     }
732     i++;
733  }
734  if ((debug=get_option("debug"))) printf("FOUND: %d %d \n",ia,is);
735
736  if ((debug=get_option("debug")))  {
737     fprintf(prt_file, "in eoption routine\n");
738     for(i=0;i<cmd->tok_list->curr;i++) {
739        fprintf(prt_file, "command(s): %s\n",cmd->tok_list->p[i]);
740     }
741  }
742
743  if ((pos = name_list_pos("seed", nl)) > -1)
744    {
745     if (nl->inform[pos])
746       {
747        seed = command_par_value("seed", cmd->clone);
748        init55(seed);
749       }
750    }
751
752  /* change only if present in command or not yet set */
753  if ((ia == 1) || (ia_seen != 1))
754  {
755       val = command_par_value("add", cmd->clone);
756       if(val == 0) {
757         if (debug)  fprintf(prt_file, "add option not set\n");
758         add_error_opt = 0;
759       }else {
760         if (debug)  fprintf(prt_file, "add option set\n");
761         add_error_opt = 1;
762       }
763  }
764
765
766  if(ia == 1) ia_seen = 1;
767
768  if ((debug=get_option("debug"))) printf("err_add eoption: %d seen: %d\n",add_error_opt,ia_seen);
769
770}
771
772// public interface
773
774int
775node_al_errors(double* errors)
776  /* returns the alignment errors of a node */
777{
778  if (current_node->p_al_err == NULL) return 0;
779  else
780  {
781    copy_double(current_node->p_al_err->a, errors,
782                current_node->p_al_err->curr);
783    return current_node->p_al_err->curr;
784  }
785}
786
787int
788node_fd_errors(double* errors)
789  /* returns the field errors of a node */
790{
791  if (current_node->p_fd_err == NULL) return 0;
792  else
793  {
794    copy_double(current_node->p_fd_err->a, errors,
795                current_node->p_fd_err->curr);
796    return current_node->p_fd_err->curr;
797  }
798}
799
800int
801node_rf_errors(double* errors, double *freq, double *harmon, double *lag )
802  /* AL: returns the phase errors of a node */
803{
804  if (current_node->p_ph_err == NULL) return 0;
805  else
806  {
807    *freq = current_node->rfm_freq;
808    *harmon = current_node->rfm_harmon;
809    *lag = current_node->rfm_lag;
810    copy_double(current_node->p_ph_err->a, errors,
811                current_node->p_ph_err->curr);
812    return current_node->p_ph_err->curr;
813  }
814}
815
816void
817pro_error(struct in_cmd* cmd)
818{
819
820 if (strcmp(cmd->tok_list->p[0], "eoption") == 0)
821     {
822      error_eoption(cmd);
823      cmd->clone_flag = 1; /* do not drop */
824      current_eopt = cmd->clone;
825      return;
826     }
827  if (get_option("debug")) fprintf(prt_file, "enter ERROR module\n");
828  if (current_sequ == NULL || current_sequ->ex_start == NULL)
829    {
830     warning("ERROR, but no active sequence:", "ignored");
831     return;
832    }
833        setbuf(stdout,(char *)0);
834
835 if (error_select->curr > 0) set_selected_errors();
836
837 if (strcmp(cmd->tok_list->p[0], "ealign") == 0)
838          {
839           error_ealign(cmd);
840          }
841 else if (strcmp(cmd->tok_list->p[0], "efield") == 0)
842          {
843           error_efield(cmd);
844          }
845 else if (strcmp(cmd->tok_list->p[0], "efcomp") == 0)
846          {
847           error_efcomp(cmd);
848          }
849 else if (strcmp(cmd->tok_list->p[0], "eprint") == 0)
850          {
851           error_eprint(cmd);
852          }
853 else if (strcmp(cmd->tok_list->p[0], "seterr") == 0)
854          {
855           error_seterr(cmd);
856          }
857 else if (strcmp(cmd->tok_list->p[0], "esave") == 0)
858          {
859           error_esave(cmd);
860          }
861}
862
Note: See TracBrowser for help on using the repository browser.