source: PSPA/madxPSPA/src/mad_orbit.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: 100.8 KB
Line 
1#include "madx.h"
2
3// private types
4
5struct val_mic {
6  double before[2];
7  double after[2];
8};
9
10struct id_mic {
11  int   id_ttb;
12  int   enable;
13  struct val_mic val;
14  struct node* p_node;
15  struct id_mic *next;
16  struct id_mic *previous;
17};
18
19struct id_mic2 {
20  int   id_ttb[2];
21  int   enable;
22  struct val_mic val;
23  struct node* p_node;
24  struct node* p_node_s1;
25  struct node* p_node_s2;
26  struct id_mic2 *next;
27  struct id_mic2 *previous;
28};
29
30struct orb_cor {
31  double qx0;
32  double qy0;
33  double units;
34  struct id_mic *cor_table;
35  struct id_mic *mon_table;
36};
37
38struct orb_cor2 {
39  double qx0;
40  double qy0;
41  double units;
42  struct id_mic2 *cor_table;
43  struct id_mic2 *mon_table;
44};
45
46// forward declarations
47
48static void    correct_correct(struct in_cmd*);
49static void    correct_usemonitor(struct in_cmd*);
50static void    correct_usekick(struct in_cmd*);
51static void    correct_putorbit(struct in_cmd*);
52static void    correct_getorbit(struct in_cmd*); /* empty */
53static void    correct_option(struct in_cmd*);
54static void    correct_readcorr(struct in_cmd*);
55static void    correct_setcorr(struct in_cmd*);
56
57static void    correct_correct1(struct in_cmd*);
58static int     pro_correct_getactive(int ip, int *nm, int *nx, int *nc, double *corvec, double *monvec,char *conm);
59static void    pro_correct_write_results(double *monvec, double *resvec, double *corvec, int *nx, int *nc, int *nm, int imon, int icor, int ip);
60static void    pro_correct_fill_mon_table(int ip ,char *name, double old, double new_);
61static void    pro_correct_fill_corr_table(int ip ,char *name, double old, double new_);
62static void    pro_correct_make_mon_table(void);
63static void    pro_correct_make_corr_table(void);
64static double* pro_correct_response_line(int ip, int nc, int nm);
65static double* pro_correct_response_ring(int ip, int nc, int nm);
66static int     pro_correct_filter(int iplane, double sigcut);
67static void    pro_correct_write_cocu_table(void);
68static void    pro_correct_prtwiss(void);
69static int     pro_correct_getcorrs(struct in_cmd*);
70static int     pro_correct_getorbit_ext(struct in_cmd*);
71static int     pro_correct_getorbit(struct in_cmd*);
72static int     pro_correct_gettables(int iplane, struct in_cmd*);
73static int     pro_correct_getcommands(struct in_cmd*);
74// static void    pro_correct_option(struct in_cmd*);
75
76static void    correct_correct2(struct in_cmd*);
77static void    pro_correct2_fill_mon_table(int ip ,char *name, double old, double new_);
78static void    pro_correct2_fill_corr_table(int b, int ip ,char *name, double old, double new_);
79static void    pro_correct2_make_mon_table(void);
80static void    pro_correct2_make_corr_table(void);
81static double* pro_correct2_response_ring(int ip, int nc, int nm);
82static int     pro_correct2_getactive(int ip, int *nm, int *nx, int *nc, double *corvec, double *monvec,char *conm);
83static int     pro_correct2_getcorrs(struct in_cmd*);
84static int     pro_correct2_getorbit(struct in_cmd*);
85static int     pro_correct2_gettables(int iplane, struct in_cmd*);
86static void    pro_correct2_write_results(double *monvec, double *resvec, double *corvec, int *nx, int *nc, int *nm, int imon, int icor, int ip);
87
88static void    fill_orbit_table(struct table* t_out, struct table* t_in);
89
90// private interface
91
92static double
93crms(double *r, int m)
94{
95  double xave = {0.0};
96  double xrms = {0.0};
97  int    i;
98
99  for(i=0; i<m; i++) {
100    xave = xave + r[i];
101  }
102  xave = xave/m;
103  for(i=0; i<m; i++) {
104    xrms = xrms + (xave - r[i])*(xave - r[i]);
105  }
106  xrms = sqrt(xrms/m);
107
108  return(xrms);
109}
110
111static double
112cprp(double *r, int m)
113{
114  double xhi  = {-9999.};
115  double xlo  = {9999.};
116  double xptp = {0.0};
117  int    i;
118
119  for(i=0; i<m; i++) {
120     if(r[i] < xlo) xlo = r[i];
121     if(r[i] > xhi) xhi = r[i];
122  }
123     xptp = xhi - xlo;
124
125  return(xptp);
126}
127
128static double
129copk(double *r, int m)
130{
131  double xpk  = {-9999.};
132  int    i;
133
134  for(i=0; i<m; i++) {
135     if(fabs(r[i]) > xpk) xpk = fabs(r[i]);
136  }
137
138  return(xpk);
139}
140
141static int
142c_micit(double *dmat,char *conm, double *monvec,double *corvec,double *resvec,int *nx,float rms,int imon,int icor,int niter)
143{
144  char rout_name[] = "c_micit";
145  int *ny;
146  int ifail;
147  float *ax,*cinx,*xinx,*resx;
148  float *rho,*ptop,*rmss,*xrms,*xptp,*xiter;
149
150  /* allocate auxiliary vectors used by correction algorithms */
151  ny    = mycalloc("c_micit_ny",icor,sizeof(int));
152  ax    = mycalloc("c_micit_ax",imon*icor,sizeof(float));
153  cinx  = mycalloc("c_micit_cinx",icor,sizeof(float));
154  xinx  = mycalloc("c_micit_xinx",imon,sizeof(float));
155  resx  = mycalloc("c_micit_resx",imon,sizeof(float));
156  rho   = mycalloc("c_micit_rho",3*icor,sizeof(float));
157  ptop  = mycalloc("c_micit_ptop",icor,sizeof(float));
158  rmss  = mycalloc("c_micit_rmss",icor,sizeof(float));
159  xrms  = mycalloc("c_micit_xrms",icor,sizeof(float));
160  xptp  = mycalloc("c_micit_xptp",icor,sizeof(float));
161  xiter = mycalloc("c_micit_xiter",icor,sizeof(float));
162
163  micit_(dmat, conm, monvec, corvec, resvec, nx, &rms, &imon, &icor, &niter, ny, ax, cinx, xinx, resx,rho,ptop,rmss,xrms,xptp,xiter,&ifail);
164
165  myfree(rout_name,ny); myfree(rout_name,ax); myfree(rout_name,cinx);
166  myfree(rout_name,xinx); myfree(rout_name,resx); myfree(rout_name,rho);
167  myfree(rout_name,ptop); myfree(rout_name,rmss); myfree(rout_name,xrms);
168  myfree(rout_name,xptp); myfree(rout_name,xiter);
169
170  return(ifail);
171}
172
173static void
174c_haveit(double *dmat,double *monvec,double *corvec,double *resvec,int *nx,int imon,int icor)
175{
176  char rout_name[] = "c_haveit";
177  double *cb,*xmeas,*xres,*y,*z,*xd;
178
179  cb   = mycalloc("c_haveit_cb",icor,sizeof(double));
180  xmeas= mycalloc("c_haveit_xmeas",imon,sizeof(double));
181  xres = mycalloc("c_haveit_xres",imon,sizeof(double));
182  y    = mycalloc("c_haveit_y",icor*imon,sizeof(double));
183  z    = mycalloc("c_haveit_z",icor*icor,sizeof(double));
184  xd   = mycalloc("c_haveit_xd",icor,sizeof(double));
185
186  haveit_(dmat,monvec,corvec,resvec,nx,&imon,&icor,cb,xmeas,xres,y,z,xd);
187
188  myfree(rout_name,cb); myfree(rout_name,xmeas); myfree(rout_name,xres);
189  myfree(rout_name,y);  myfree(rout_name,z); myfree(rout_name,xd);
190
191  return;
192}
193
194static int
195c_svddec(double *dmat, int imon, int icor, int *sing, double *sngcut, double *sngval)
196
197{
198  char rout_name[] = "c_svddev";
199  int    flag;
200  int    dbg;
201
202  double *s, *u, *v, *w, *ut, *vt, *wt;
203  double *ws, *wv;
204  int    *sw;
205
206  s   = mycalloc("c_svddec_s",icor*imon,sizeof(double));
207  u   = mycalloc("c_svddec_u",icor*imon,sizeof(double));
208  v   = mycalloc("c_svddec_v",icor*imon,sizeof(double));
209  w   = mycalloc("c_svddec_w",icor*imon,sizeof(double));
210  ut  = mycalloc("c_svddec_ut",icor*imon,sizeof(double));
211  vt  = mycalloc("c_svddec_vt",icor*imon,sizeof(double));
212  wt  = mycalloc("c_svddec_wt",icor*imon,sizeof(double));
213  ws  = mycalloc("c_svddec_ws",icor,sizeof(double));
214  wv  = mycalloc("c_svddec_wv",icor,sizeof(double));
215  sw  = mycalloc("c_svddec_sw",icor,sizeof(int));
216
217  dbg = debug_correct_opt;
218
219  if(imon >= icor ) {
220      svddec_m_(dmat,s,u,v,w,ut,vt,wt,ws,wv,sw,sngcut,sngval,&imon,&icor,&flag,sing,&dbg);
221  } else {
222      svddec_c_(dmat,s,u,v,w,ut,vt,wt,ws,wv,sw,sngcut,sngval,&imon,&icor,&flag,sing,&dbg);
223  }
224  myfree(rout_name,s); myfree(rout_name,u);
225  myfree(rout_name,v); myfree(rout_name,w);
226  myfree(rout_name,ut); myfree(rout_name,vt);
227  myfree(rout_name,wt); myfree(rout_name,ws);
228  myfree(rout_name,wv); myfree(rout_name,sw);
229
230  return(flag);
231}
232
233static int 
234c_svdcorr(double *dmat, double *xin, double *cor, double *res, int *nx, int imon, int icor)
235{
236  char rout_name[] = "c_svdcorr";
237  int    flag;
238  int    dbg;
239
240  double *s, *u, *v, *w, *ut, *vt, *wt;
241  double *xa, *xb, *xp, *wv, *ws;
242  int    *sw;
243
244  s   = mycalloc("c_svdcorr_s",icor*imon,sizeof(double));
245  u   = mycalloc("c_svdcorr_u",icor*imon,sizeof(double));
246  v   = mycalloc("c_svdcorr_v",icor*imon,sizeof(double));
247  w   = mycalloc("c_svdcorr_w",icor*imon,sizeof(double));
248  ut  = mycalloc("c_svdcorr_ut",icor*imon,sizeof(double));
249  vt  = mycalloc("c_svdcorr_vt",icor*imon,sizeof(double));
250  wt  = mycalloc("c_svdcorr_wt",icor*imon,sizeof(double));
251
252  xa  = mycalloc("c_svdcorr_xa",imon,sizeof(double));
253  xb  = mycalloc("c_svdcorr_xb",imon,sizeof(double));
254  xp  = mycalloc("c_svdcorr_xp",imon,sizeof(double));
255  ws  = mycalloc("c_svdcorr_xp",icor,sizeof(double));
256  wv  = mycalloc("c_svdcorr_xp",icor,sizeof(double));
257
258  sw  = mycalloc("c_svdcorr_sw",icor,sizeof(int));
259
260  dbg = debug_correct_opt;
261
262  if(imon >= icor ) {
263      svdcorr_m_(dmat,s,u,v,w,ut,vt,wt,xin,cor,res,
264                 xa,xb,xp,ws,wv,sw,
265                 nx,&imon,&icor,&flag,&dbg);
266  } else {
267      svdcorr_c_(dmat,s,u,v,w,ut,vt,wt,xin,cor,res,
268                 xa,xb,xp,ws,wv,sw,
269                 nx,&imon,&icor,&flag,&dbg);
270  }
271  myfree(rout_name,s); myfree(rout_name,u); myfree(rout_name,v);
272  myfree(rout_name,w); myfree(rout_name,ut); myfree(rout_name,vt);
273  myfree(rout_name,wt); myfree(rout_name,sw); myfree(rout_name,xa);
274  myfree(rout_name,xb); myfree(rout_name,xp); myfree(rout_name,ws);
275  myfree(rout_name,wv);
276
277  return(flag);
278}
279
280static void
281fill_orbit_table(struct table* t_out, struct table* t_in)
282  /* fills a table with orbit values at monitor positions */
283{
284  int i, j, pos;
285  t_out->curr = 0;
286  for (i = 0; i < t_in->curr; i++)
287  {
288    if (strstr(t_in->s_cols[1][i], "monitor"))
289    {
290      for (j = 0; j < t_out->num_cols; j++)
291      {
292        if ((pos = name_list_pos(t_out->columns->names[j],
293                                 t_in->columns)) > -1)
294        {
295          if (t_out->columns->inform[j] < 3)
296            t_out->d_cols[j][t_out->curr] = t_in->d_cols[pos][i];
297          else t_out->s_cols[j][t_out->curr]
298                 = tmpbuff(t_in->s_cols[pos][i]);
299        }
300        else
301        {
302          if (t_out->columns->inform[j] < 3)
303            t_out->d_cols[j][t_out->curr] = zero;
304          else t_out->s_cols[j][t_out->curr] = tmpbuff(blank);
305        }
306      }
307      t_out->curr++;
308    }
309  }
310}
311
312/* revert to old version after Thys Risselada's fix of Micado */ 
313
314static void
315correct_setcorr(struct in_cmd* cmd)
316{
317
318/* read the correctors from named table and stores
319   them in the nodes of the sequence at
320   "->chkick" and "->cvkick". Subsequent Twiss will
321   use them correctly.
322   ===> Must be preceded by a call to "read_table"
323   ===> (unless table exists in memory !)
324   ===> Watch out, does not yet take care of existing corrector
325   ===> settings already present in sequence
326   ===> Uses table with name specified in parameter: table=
327*/
328
329  int i, ix;
330
331  struct node *ndexe;
332  struct node *nextnode;
333
334  char     name[NAME_L];
335  char   slname[NAME_L];
336
337  char     nname[NAME_L];
338  char   slnname[NAME_L];
339
340  char*    namtab;
341//  int      t1;
342
343  double   xnew, ynew;
344
345/* set up pointers to current sequence for later use */
346  struct sequence* mysequ = current_sequ;
347  nextnode = mysequ->ex_start;
348  ndexe = mysequ->ex_end;
349
350/* printf("Pointers: %d %d %d\n",mysequ,nextnode,ndexe); */
351
352  if ((namtab = command_par_string("table",cmd->clone)) != NULL) {
353       printf("Want to use named table: %s\n",namtab);
354       if (name_list_pos(namtab, table_register->names) > -1) { // (t1 = not used
355          printf("The table ==> %s <=== was found \n",namtab);
356       } else {
357          /* fatal_error("Corrector table requested, but not existing:",namtab); */
358          /* exit(-77); */ 
359          printf("No such corrector table in memory: %s\n",namtab);
360       }
361
362  } else {
363       if (get_option("debug")) {
364         printf("No table name requested\n");
365         printf("Use default name\n");
366       }
367       strcpy(namtab,"corr");
368  }
369
370
371  i = 1; ix=0;
372  while(ix == 0) {
373      ix = string_from_table_row(namtab, "name", &i, name);
374      ix = double_from_table_row(namtab, "px.correction", &i, &xnew);
375      ix = double_from_table_row(namtab, "py.correction", &i, &ynew);
376      if(ix == 0) {
377             stolower(name);
378             strcpy(slname,strip(name));
379             supp_tb(slname);
380
381          /* printf("corrs: %s %d %e %e %e %e\n",name,ix,xold,yold,xnew,ynew); */   
382          nextnode = mysequ->ex_start;
383          while (nextnode != ndexe) {
384             stolower(name);
385             strcpy(slname,strip(name));
386             supp_tb(slname);
387           
388             strcpy(nname,nextnode->name);
389             stolower(nname);
390             strcpy(slnname,strip(nname));
391             supp_tb(slnname);
392         
393             /* printf("seq and input (0): %s %d %s %d\n", nname,strlen(nname),  name,strlen(name));
394             printf("seq d in (2): %s %d %s %d\n",slnname,strlen(slnname),slname,strlen(slname)); */
395       
396             if(strcmp(slname,slnname) == 0) {
397                /*
398                printf("Corrector selection found: %s, %s %d\n",lname,nextnode->name,nextnode->sel_err);
399                printf("corrs: %s %d %e %e %e %e\n",name,ix,xold,yold,xnew,ynew);   
400                printf("corrs in sequence: %s %e %e\n",nextnode->name,nextnode->chkick,nextnode->cvkick);
401                */
402                nextnode->chkick += xnew;
403                nextnode->cvkick += ynew;
404                /*
405                printf("corrs in sequence: %s %e %e\n",nextnode->name,nextnode->chkick,nextnode->cvkick);
406                */
407                nextnode = ndexe;
408             } else {
409                nextnode = nextnode->next;
410             }
411          }
412   
413      } 
414        i++;
415  }
416  return;
417}
418
419static void
420correct_readcorr(struct in_cmd* cmd)
421{
422
423/* read the correctors from table "corr" and stores
424   them in the nodes of the sequence at
425   "->chkick" and "->cvkick". Subsequent Twiss will
426   use them correctly.
427   ===> Must be preceded by a call to "read_table"
428   ===> Watch out, does not yet take care of existing corrector
429   ===> settings already present in sequence
430   ===> Always uses table with name "corr", will change ...
431*/
432
433  int i, ix;
434
435  struct node *ndexe;
436  struct node *nextnode;
437
438  char     name[NAME_L];
439  char    lname[NAME_L];
440  char   slname[NAME_L];
441  char* uslname;
442
443  char     nname[NAME_L];
444  char    lnname[NAME_L];
445  char   slnname[NAME_L];
446  char* uslnname;
447
448  double   xnew, ynew;
449
450/* set up pointers to current sequence for later use */
451  struct sequence* mysequ = current_sequ;
452  nextnode = mysequ->ex_start;
453  ndexe = mysequ->ex_end;
454
455/* printf("Pointers: %d %d %d\n",mysequ,nextnode,ndexe); */
456
457  (void)cmd;
458  i = 1; ix=0;
459  while(ix == 0) {
460      ix = string_from_table_row("corr", "name", &i, name);
461      ix = double_from_table_row("corr", "px.correction", &i, &xnew);
462      ix = double_from_table_row("corr", "py.correction", &i, &ynew);
463      if(ix == 0) {
464          /* printf("corrs: %s %d %e %e %e %e\n",name,ix,xold,yold,xnew,ynew); */   
465          nextnode = mysequ->ex_start;
466          while (nextnode != ndexe) {
467             strcpy(lname,name);
468             stolower(lname);
469             strcpy(slname,strip(lname));
470             uslname = supp_tb(slname);
471           
472             strcpy(nname,nextnode->name);
473             strcpy(lnname,nname);
474             stolower(lnname);
475             strcpy(slnname,strip(lnname));
476             uslnname = supp_tb(slnname);
477         
478             /* printf("seq and input (0): %s %d %s %d\n", nname,strlen(nname),  name,strlen(name));
479             printf("seq d in (1): %s %d %s %d\n",lnname,strlen(lnname),lname,strlen(lname));
480             printf("seq d in (2): %s %d %s %d\n",slnname,strlen(slnname),slname,strlen(slname));
481             printf("seq d in (3): %s %d %s %d\n",uslnname,strlen(uslnname),uslname,strlen(uslname));
482             printf("compare: %s %d %s %d \n",uslname,strlen(uslname),uslnname,strlen(uslnname)); */
483       
484             if(strcmp(uslname,uslnname) == 0) {
485                /*
486                printf("Corrector selection found: %s, %s %d\n",lname,nextnode->name,nextnode->sel_err);
487                printf("corrs: %s %d %e %e %e %e\n",name,ix,xold,yold,xnew,ynew);   
488                printf("corrs in sequence: %s %e %e\n",nextnode->name,nextnode->chkick,nextnode->cvkick);
489                */
490                nextnode->chkick += xnew;
491                nextnode->cvkick += ynew;
492                /*
493                printf("corrs in sequence: %s %e %e\n",nextnode->name,nextnode->chkick,nextnode->cvkick);
494                */
495                nextnode = ndexe;
496             } else {
497                nextnode = nextnode->next;
498             }
499          }
500   
501      } 
502      i++;
503  }
504}
505
506static void
507correct_correct2(struct in_cmd* cmd)
508/* Steering routine for orbit corrections of two beams */
509{
510  char rout_name[] = "correct_correct2";
511
512/*
513  struct name_list* spos = sequences->list;
514  struct table *twb1;
515  struct table *twb2;
516  int idrop;
517  int pos;
518  struct timeb tp;
519  int sflag, svdflg;
520  double  sigcut;         
521*/
522
523  int ix, im, ip; // , it; not used
524  int i,j,nnnseq; // ,err // not used
525  int imon, icor;
526  int ncorr, nmon;
527  int niter;
528//  int resout; not used
529  int twism;
530  int ifail;
531  float  rms;
532  double rrms;
533  double tmp1, tmp2, tmp3, tmp4;
534  char    *clist, *mlist;    /* file names for monitor and corrector output */
535  char    clist1[100], clist2[100];   /* file names for corrector output ring 1 and ring 2 */
536  double  *dmat = {NULL};    /* response matrix, double precision */
537  double  *corvec, *monvec;  /* vectors to hold measured orbit and correctors */
538  double  *resvec;           /* vector to hold corrected orbit */
539  char    *conm;             /* vector to hold corrector names (for MICADO) */
540//  int     *sing;  // not used /* array to store pointer to singular correctors */
541  static int     *nm, *nx, *nc;
542  struct id_mic2  *c;
543  /*
544  struct id_mic2   *m;
545  */
546
547  strcpy(clist1,"\0");
548  strcpy(clist2,"\0");
549
550  printf("for two beams orbit corrections ...\n");
551  ip = pro_correct_getcommands(cmd);
552  im = pro_correct2_gettables(ip,cmd);
553  ncorr = im%10000; nmon  = im/10000;
554  printf("%d monitors and %d correctors found in input\n",nmon,ncorr);
555  if(nmon == 0) {
556    printf("No monitor found in input, no correction done\n");
557    return;
558  }
559  if(ncorr == 0) {
560    printf("No corrector found in input, no correction done\n");
561    return;
562  }
563
564
565
566  /* For debugging set output buffer to zero */
567  if (get_option("debug"))  setbuf(stdout,NULL);
568
569
570  /* Prepare file descriptors for the output */
571  if(command_par_value("resout",cmd->clone) > 0) { // (resout = not used
572     if(fddata == NULL) {
573        if((fddata = fopen("corr.out","w")) == NULL)
574           exit(99);
575     }
576     if(fcdata == NULL) {
577        if((fcdata = fopen("stren.out","w")) == NULL)
578           exit(99);
579     }
580     if(fgdata == NULL) {
581        if((fgdata = fopen("plot.orb","w")) == NULL)
582           exit(99);
583     }
584  }
585
586
587  /* If only Twiss summary is required prepare and write it */
588  if((twism = command_par_value("twissum",cmd->clone)) > 0) {
589     if(ftdata == NULL) {
590        if((ftdata = fopen("twiss.summ","w")) == NULL)
591           exit(99);
592     }
593     j = 1;
594     if((nnnseq = get_variable("n")) == 0) {
595         nnnseq = twism;
596     }
597     double_from_table_row("summ", "xcomax",&j,&tmp1); // err = not used
598     double_from_table_row("summ", "xcorms",&j,&tmp2); // err = not used
599     double_from_table_row("summ", "ycomax",&j,&tmp3); // err = not used
600     double_from_table_row("summ", "ycorms",&j,&tmp4); // err = not used
601     fprintf(ftdata," T: %d %e %e %e %e\n",nnnseq,tmp1,tmp2,tmp3,tmp4);
602     return;
603  }
604
605  /* allocate vectors used by correction algorithms */
606  nx     = mycalloc("correct_correct2_nx",ncorr,sizeof(int));
607  nc     = mycalloc("correct_correct2_nc",ncorr,sizeof(int));
608  nm     = mycalloc("correct_correct2_nm",nmon,sizeof(int));
609  // sing   = mycalloc("correct_correct2_sing",ncorr*2,sizeof(int)); // not used
610  corvec = mycalloc("correct_correct2_corvec",ncorr,sizeof(double));
611  monvec = mycalloc("correct_correct2_monvec",nmon,sizeof(double));
612  resvec = mycalloc("correct_correct2_resvec",nmon,sizeof(double));
613  conm   = mycalloc("correct_correct2_conm",ncorr*16,sizeof(char));
614
615  /* get original settings of correctors from input Twiss-table */
616  pro_correct2_getcorrs(cmd); // it = not used
617  /* get input orbit, default is from input Twiss-table */
618  pro_correct2_getorbit(cmd);  // it = not used
619
620
621  /* find and prepare enabled correctors and monitors, may be repeated */
622  ix = pro_correct2_getactive(ip, nm, nx, nc, corvec, monvec, conm);
623  icor = ix%10000; imon  = ix/10000;
624  printf("%d monitors and %d correctors enabled\n",imon,icor);
625
626
627    if (get_option("debug")) {
628  for (i=0;i<icor;i++) {
629    printf("C: %d %d \n",nx[i],nc[i]);
630  }
631  for (i=0;i<imon;i++) {
632    printf("M: %d %e \n",nm[i],monvec[i]);
633  }
634  }
635
636  if(strcmp("ring",command_par_string("flag",cmd->clone)) == 0) {
637    if(dmat != NULL) myfree(rout_name,dmat);
638    /* icor and imon used to set up correct matrix size !! */
639    dmat =  pro_correct2_response_ring(ip,icor,imon);
640  }
641  else { printf("INVALID MACHINE TYPE\n"); exit(-1);
642  }
643
644  /* MICADO correction, get desired number of correctors from command */
645  corrl = command_par_value("corrlim",cmd->clone);
646  set_variable("corrlim",&corrl);
647  if(strcmp("micado",command_par_string("mode",cmd->clone)) == 0) {
648    printf("enter MICADO correction ...\n");
649    if((niter = command_par_value("ncorr",cmd->clone)) == 0) {
650          printf("Requested %d correctors (\?\?\?) set to %d\n",niter,icor);
651          niter = icor;
652    }
653    else if((niter = command_par_value("ncorr",cmd->clone)) < 0) {
654          printf("Requested %d correctors (\?\?\?) set to 0\n",niter);
655          niter = 0;
656    }
657    else if((niter = command_par_value("ncorr",cmd->clone)) > icor) {
658          printf("Fewer correctors available than requested by ncorr\n");
659          printf("you want %d,  you get %d\n",niter,icor);
660          printf("ncorr reset to %d\n",icor);
661          niter = icor;
662    }
663    rms  = 1000.0*command_par_value("error",cmd->clone);
664    /*frs       micit_(dmat,monvec,corvec,resvec,nx,&rms,&imon,&icor,&niter); */
665    /* printf("Time before micado:  %-6.3f\n",fextim());  */
666    ifail = c_micit(dmat,conm,monvec,corvec,resvec,nx,rms,imon,icor,niter);
667    printf("Back from micado %d\n",ifail);
668    if(ifail != 0) {
669       printf("MICADO correction completed with error code %d\n\n",ifail);
670       warning("MICADO back with error",", no correction done");
671    }
672    rrms = crms(monvec,imon);
673    printf("RMS before %e\n",rrms);
674    rrms = crms(resvec,imon);
675    printf("RMS after  %e\n",rrms);
676    if(fgdata != NULL) {
677    for (i=0; i<nmon; i++) {
678       fprintf(fgdata,"%e %e \n",monvec[i],resvec[i]);
679    }
680    }
681/*
682    for (i=0; i<nmon; i++) {
683       printf("monvec: %d %e \n",i,monvec[i]);
684    }
685    printf("\n");
686    for (i=0; i<nmon; i++) {
687       printf("resvec: %d %e \n",i,resvec[i]);
688    }
689    m = correct_orbit12->mon_table;
690    for (i=0; i<nmon; i++) {
691       printf("resvec: %s %e \n",m[nm[i]].p_node->name,resvec[i]);
692    }
693    printf("\n");
694    for (i=0; i<ncorr; i++) {
695       printf("corvec: %d %e \n",i,corvec[i]);
696    }
697    printf("\n");
698*/
699
700    c = correct_orbit12->cor_table;
701    for(i=0;i<icor;i++) {
702       printf("%s %e\n",c[nc[i]].p_node->name,corvec[nx[i]-1]);
703    }
704    printf("\n");
705 
706    /* printf("Time after micado:  %-6.3f\n",fextim());   */
707    if(ifail != 0) {
708       printf("MICADO correction completed with error code %d\n\n",ifail);
709       warning("MICADO back with error",", no correction done");
710    }
711    if(ifail == 0) {
712       pro_correct2_write_results(monvec, resvec, corvec, nx, nc, nm, imon, icor, ip);
713    }
714  }
715 
716  /* write corrector output to tfs table */
717  if ((clist = command_par_string("clist",cmd->clone)) != NULL) {
718    strcat(clist1,clist);
719    strcat(clist1,"_1");
720    strcat(clist2,clist);
721    strcat(clist2,"_2");
722    out_table("corr1",corr_table1,clist1);
723    out_table("corr2",corr_table2,clist2);
724  }
725
726  /* write monitor output to tfs table */
727  if ((mlist = command_par_string("mlist",cmd->clone)) != NULL) {
728    out_table("mon",mon_table,mlist);
729  }
730
731
732  /* Clean up at the end of the module */
733   
734  myfree(rout_name,nm);myfree(rout_name,dmat);myfree(rout_name,nx);
735  myfree(rout_name,nc);myfree(rout_name,corvec);
736  myfree(rout_name,monvec);myfree(rout_name,resvec); myfree(rout_name,conm);
737}
738
739static int
740pro_correct2_gettables(int iplane, struct in_cmd* cmd)
741{
742
743  char rout_name[] = "pro_correct2_gettables";
744
745  struct id_mic2 *cor_l1,  *cor_l2;
746  struct id_mic2 *mon_l1,  *mon_l2;
747  struct id_mic2 *cor_l12, *mon_l12;
748  struct id_mic2 *prt;
749
750  // struct table *ttb; // not used
751
752  struct table *b1 = NULL;
753  struct table *b2 = NULL;
754
755  char* orbtab1;
756  char* orbtab2;
757
758  int t1, t2;   
759
760  int ebl1, ebl2;
761
762  int j,k;
763//  int set0; // not used
764  int cntm1 = {0};
765  int cntc1 = {0};
766  int cntm2 = {0};
767  int cntc2 = {0};
768  int cntm12 = {0};
769  int cntc12 = {0};
770 
771  double ounits;
772
773/*
774  static char atm[6][4] = {"hmon","vmon","moni","hkic","vkic","kick"};
775*/
776
777/* Get access to tables, for orbit and model the default is twiss_table */
778
779  if ((orbtab1 = command_par_string("beam1tab",cmd->clone)) != NULL) {
780      printf("Want to use orbit from: %s\n",orbtab1);
781    if ((t1 = name_list_pos(orbtab1, table_register->names)) > -1) {
782       b1 = table_register->tables[t1];
783    } else {
784       fatal_error("Beam 1 ORBIT table requested, but not provided:",orbtab1);
785    }
786  } else {
787
788  }
789  if ((orbtab2 = command_par_string("beam2tab",cmd->clone)) != NULL) {
790      printf("Want to use orbit from: %s\n",orbtab2);
791    if ((t2 = name_list_pos(orbtab2, table_register->names)) > -1) {
792       b2 = table_register->tables[t2];
793    } else {
794       fatal_error("Beam 2 ORBIT table requested, but not provided:",orbtab2);
795    }
796  } else {
797
798  }
799
800/* store as globals for later use */
801  if((b1 != NULL) && (b2 != NULL)) {
802     twiss_table_beam1 = b1;
803     twiss_table_beam2 = b2;
804  } else {
805       fatal_error("Beam 1 and 2 orbit tables not found:",orbtab1);
806  }
807 
808/* reserve space for orbit correction structures */
809  if(correct_orbit12 == NULL) {
810    correct_orbit12 = (struct orb_cor2*)mycalloc("pro_correct2_gettables",1, 
811                     sizeof(struct orb_cor2));
812  }
813
814
815
816  if(correct_orbit12->cor_table != NULL) myfree(rout_name,correct_orbit12->cor_table);
817  if(correct_orbit12->mon_table != NULL) myfree(rout_name,correct_orbit12->mon_table);
818
819  correct_orbit12->cor_table = (struct id_mic2 *)mycalloc("pro_correct2_gettables_cor",5200, sizeof(struct id_mic2));
820  correct_orbit12->mon_table = (struct id_mic2 *)mycalloc("pro_correct2_gettables_mon",5200, sizeof(struct id_mic2));
821
822/* orbit table available, get units, if defined */
823   if((ounits = command_par_value("units",cmd->clone)) > 0) {
824          correct_orbit12->units=ounits;                     
825   } else {
826          correct_orbit12->units=1.0;     
827   }
828
829  // ttb = model_table; // not used
830/* no more need, we have b1 and b2 as pointers .. */
831
832  correct_orbit12->mon_table->previous = NULL;
833  correct_orbit12->mon_table->next = NULL;
834  correct_orbit12->cor_table->previous = NULL;
835  correct_orbit12->cor_table->next = NULL;
836
837  mon_l1 = correct_orbit12->mon_table;
838  cor_l1 = correct_orbit12->cor_table;
839 
840  for (j=0; j < b1->curr; j++) {
841   if((strncmp(atm[iplane-1],b1->p_nodes[j]->base_name,4) == 0) ||
842      (strncmp(atm[2],      b1->p_nodes[j]->base_name,4) == 0))  {
843/*    printf("1m: %s %ld\n", b1->p_nodes[j]->name, strstr(".b2", b1->p_nodes[j]->name)); */
844      if(strstr(b1->p_nodes[j]->name,".b1") != NULL) {
845        mon_l1->id_ttb[0] = j;
846        mon_l1->id_ttb[1] = -1;
847        mon_l1->enable = b1->p_nodes[j]->enable;
848        mon_l1->p_node = b1->p_nodes[j];
849        mon_l1->next = mon_l1;
850        mon_l1->next++; mon_l1++;
851        cntm1++;
852      } else {
853/*      printf("Removed: %s\n",b1->p_nodes[j]->name); */
854      }
855    }
856    if((strncmp(atc[iplane-1],b1->p_nodes[j]->base_name,4) == 0) ||
857       (strncmp(atc[2],       b1->p_nodes[j]->base_name,4) == 0))  {
858/*    printf("1c: %s %ld\n", b1->p_nodes[j]->name, b1->p_nodes[j]->name); */
859      if(strstr(b1->p_nodes[j]->name,".b1") != NULL) {
860        cor_l1->id_ttb[0] = j;
861        cor_l1->id_ttb[1] = -1;
862        cor_l1->enable = b1->p_nodes[j]->enable;
863        cor_l1->p_node = b1->p_nodes[j];
864        cor_l1->p_node_s1 = b1->p_nodes[j];
865        cor_l1->p_node_s2 = NULL;           
866        if(command_par_value("corzero",cmd->clone) > 0) { // (set0 = not used
867          if(iplane == 1) cor_l1->p_node_s1->chkick = 0.0;
868          if(iplane == 2) cor_l1->p_node_s1->cvkick = 0.0;
869        }
870        cor_l1->next = cor_l1;
871        cor_l1->next++; cor_l1++;
872        cntc1++;
873      } else {
874/*      printf("Removed: %s\n",b1->p_nodes[j]->name); */
875      }
876    }
877  }
878
879  mon_l2 = mon_l1;
880  cor_l2 = cor_l1;
881  for (j=0; j < b2->curr; j++) {
882   if((strncmp(atm[iplane-1],b2->p_nodes[j]->base_name,4) == 0) ||
883      (strncmp(atm[2],      b2->p_nodes[j]->base_name,4) == 0))  {
884/*    printf("2m: %s %ld\n", b2->p_nodes[j]->name, b2->p_nodes[j]->name); */
885      if(strstr(b2->p_nodes[j]->name,".b2") != NULL) {
886        mon_l2->id_ttb[0] = -1;
887        mon_l2->id_ttb[1] = j;
888        mon_l2->enable = b2->p_nodes[j]->enable;
889        mon_l2->p_node = b2->p_nodes[j];
890        mon_l2->next = mon_l2;
891        mon_l2->next++; mon_l2++;
892        cntm2++;
893      } else {
894/*      printf("Removed: %s\n",b2->p_nodes[j]->name); */
895      }
896    }
897    if((strncmp(atc[iplane-1],b2->p_nodes[j]->base_name,4) == 0) ||
898       (strncmp(atc[2],       b2->p_nodes[j]->base_name,4) == 0))  {
899/*    printf("2c: %s %ld\n", b2->p_nodes[j]->name, b2->p_nodes[j]->name); */
900      if(strstr(b2->p_nodes[j]->name,".b2") != NULL) {
901        cor_l2->id_ttb[0] = -1;
902        cor_l2->id_ttb[1] = j;
903        cor_l2->enable = b2->p_nodes[j]->enable;
904        cor_l2->p_node = b2->p_nodes[j];
905        cor_l2->p_node_s2 = b2->p_nodes[j];
906        cor_l2->p_node_s1 = NULL;           
907        if(command_par_value("corzero",cmd->clone) > 0) { // (set0 = not used
908          if(iplane == 1) cor_l2->p_node_s2->chkick = 0.0;
909          if(iplane == 2) cor_l2->p_node_s2->cvkick = 0.0;
910        }
911        cor_l2->next = cor_l2;
912        cor_l2->next++; cor_l2++;
913        cntc2++;
914      } else {
915/*      printf("Removed: %s\n",b2->p_nodes[j]->name); */
916      }
917    }
918  }
919
920  mon_l12 = mon_l2;
921  cor_l12 = cor_l2;
922  for (j=0; j < b1->curr; j++) {
923   if((strncmp(atm[iplane-1],b1->p_nodes[j]->base_name,4) == 0) ||
924      (strncmp(atm[2],      b1->p_nodes[j]->base_name,4) == 0))  {
925/*    printf("12m: %s \n", b1->p_nodes[j]->name); */
926      if((strstr(b1->p_nodes[j]->name,".b1") == NULL)  &&
927         (strstr(b1->p_nodes[j]->name,".b2") == NULL)) {
928        mon_l12->id_ttb[0] = j;
929         for (k=0; k < b2->curr; k++) {
930           if(strcmp(b2->p_nodes[k]->name,b1->p_nodes[j]->name) == 0) {
931            mon_l12->id_ttb[1] = k;
932           }
933         }
934        mon_l12->enable = b1->p_nodes[j]->enable;
935        mon_l12->p_node = b1->p_nodes[j];
936        mon_l12->next = mon_l12;
937        mon_l12->next++; mon_l12++;
938        cntm12++;
939      } else {
940/*      printf("Removed: %s\n",b1->p_nodes[j]->name); */
941      }
942    }
943    if((strncmp(atc[iplane-1],b1->p_nodes[j]->base_name,4) == 0) ||
944       (strncmp(atc[2],       b1->p_nodes[j]->base_name,4) == 0))  {
945/*    printf("12c: %s \n", b1->p_nodes[j]->name);     */
946      if((strstr(b1->p_nodes[j]->name,".b1") == NULL)  &&
947         (strstr(b1->p_nodes[j]->name,".b2") == NULL)) {
948         cor_l12->id_ttb[0] = j;
949         for (k=0; k < b2->curr; k++) {
950           if(strcmp(b2->p_nodes[k]->name,b1->p_nodes[j]->name) == 0) {
951            cor_l12->id_ttb[1] = k;
952           }
953         }
954        cor_l12->p_node = b1->p_nodes[j];
955        cor_l12->p_node_s1 = b1->p_nodes[cor_l12->id_ttb[0]];
956        cor_l12->p_node_s2 = b2->p_nodes[cor_l12->id_ttb[1]];
957        ebl1 = b1->p_nodes[cor_l12->id_ttb[0]]->enable;
958        ebl2 = b2->p_nodes[cor_l12->id_ttb[1]]->enable;
959        cor_l12->enable = ebl1*ebl2;                 
960        if(command_par_value("corzero",cmd->clone) > 0) { // (set0 = not used
961          if(iplane == 1) cor_l12->p_node_s1->chkick = 0.0;
962          if(iplane == 2) cor_l12->p_node_s1->cvkick = 0.0;
963          if(iplane == 1) cor_l12->p_node_s2->chkick = 0.0;
964          if(iplane == 2) cor_l12->p_node_s2->cvkick = 0.0;
965        }
966        cor_l12->next = cor_l12;
967        cor_l12->next++; cor_l12++;
968        cntc12++;
969      } else {
970/*      printf("Removed: %s\n",b1->p_nodes[j]->name); */
971      }
972    }
973  }
974  /* terminate linked list   */
975  mon_l12--; mon_l12->next = NULL;
976  cor_l12--; cor_l12->next = NULL;
977
978  printf("mons and corrs (beam 1)   : %ld %ld\n",(long int)cntm1, (long int)cntc1);
979  printf("mons and corrs (beam 2)   : %ld %ld\n",(long int)cntm2, (long int)cntc2);
980  printf("mons and corrs (beam 1+2) : %ld %ld\n",(long int)cntm12, (long int)cntc12);
981
982    if (get_option("debug")) {
983     prt = correct_orbit12->mon_table;
984     while(prt != NULL) {
985       printf("Monitors beam12: %s %ld %ld\n",prt->p_node->name,(long int)prt->id_ttb[0],(long int)prt->id_ttb[1]);
986       prt = prt->next;
987     }
988
989     prt = correct_orbit12->cor_table;
990     while(prt != NULL) {
991       printf("Correctors beam12: %s %ld %ld\n",prt->p_node->name,(long int)prt->id_ttb[0],(long int)prt->id_ttb[1]);
992       prt = prt->next;
993     }
994    }
995
996/*
997     prt = correct_orbit12->cor_table;
998     while(prt != NULL) {
999       printf("Correctors beam12: %s %ld %ld\n",prt->p_node->name,prt->id_ttb[0],prt->id_ttb[1]);
1000       for (j=0; j < b2->curr; j++) {
1001         if(strcmp(b2->p_nodes[j]->name,prt->p_node->name) == 0) {
1002            prt->id_ttb[1] = j;
1003            printf("matched correctors beam12: %s %ld %ld\n",prt->p_node->name,prt->id_ttb[0],prt->id_ttb[1]);
1004         }
1005       }
1006       prt = prt->next;
1007       
1008     }
1009*/
1010
1011  if(corr_table1 == NULL) {
1012    corr_table1 = make_table("corr1", "corr1", corr_table_cols,
1013            corr_table_types, 15000);
1014    add_to_table_list(corr_table1, table_register);
1015  }
1016  if(corr_table2 == NULL) {
1017    corr_table2 = make_table("corr2", "corr2", corr_table_cols,
1018            corr_table_types, 15000);
1019    add_to_table_list(corr_table2, table_register);
1020  }
1021  pro_correct2_make_corr_table();
1022
1023  if(mon_table == NULL) {
1024    mon_table = make_table("mon", "mon", mon_table_cols,
1025           mon_table_types, 15000);
1026    add_to_table_list(mon_table, table_register);
1027    pro_correct2_make_mon_table();
1028  }
1029
1030  return(10000*(cntm1+ cntm2+ cntm12) + (cntc1 + cntc2 + cntc12));
1031}
1032
1033static int
1034pro_correct2_getorbit(struct in_cmd* cmd)
1035{
1036  struct name_list* nl;
1037  int i;
1038  int pos;
1039/*
1040  int pps, ppt;
1041*/
1042  struct id_mic2 *m;  /* access to tables for monitors and correctors */
1043  double **da1;
1044  double **da2;
1045  double xlimit;
1046
1047  char   strx[40];
1048  char   stry[40];
1049
1050  int    posx, posy, pospx, pospy;
1051
1052  da1 = twiss_table_beam1->d_cols;
1053  da2 = twiss_table_beam2->d_cols;
1054
1055  nl = cmd->clone->par_names;
1056
1057  m = correct_orbit12->mon_table;
1058
1059  strcpy(strx,"x");
1060  strcpy(stry,"y");
1061
1062  if((posx = name_list_pos(strx,twiss_table_beam1->columns)) < 0) { 
1063      fatal_error("orbit x not found in input table",", MAD-X terminates ");
1064  }
1065  if((posy = name_list_pos(stry,twiss_table_beam1->columns)) < 0) { 
1066      fatal_error("orbit y not found in input table",", MAD-X terminates ");
1067  }
1068  if (get_option("debug")) {
1069    if((pospx = name_list_pos("px",twiss_table_beam1->columns)) < 0) { 
1070        warning("orbit px not found in input table",", MAD-X continues ");
1071    }
1072    if((pospy = name_list_pos("py",twiss_table_beam1->columns)) < 0) { 
1073        warning("orbit py not found in input table",", MAD-X continues ");
1074    }
1075    printf("====c1===>  %d %d %d %d \n",posx,posy,pospx,pospy);
1076  }
1077
1078
1079  while(m) {
1080
1081/* If correction to target orbit, subtract the wanted orbit ... */
1082    if(m->id_ttb[0] > 0) {
1083    m->val.before[0] = m->p_node->other_bv*da1[ 9][m->id_ttb[0]];
1084    m->val.before[1] = m->p_node->other_bv*da1[11][m->id_ttb[0]];
1085    m->val.before[0] = m->p_node->other_bv*da1[ 9][m->id_ttb[0]]*1000.;
1086    m->val.before[1] = m->p_node->other_bv*da1[11][m->id_ttb[0]]*1000.;
1087    } else if (m->id_ttb[1] > 0) {
1088    m->val.before[0] = m->p_node->other_bv*da2[ 9][m->id_ttb[1]];
1089    m->val.before[1] = m->p_node->other_bv*da2[11][m->id_ttb[1]];
1090    m->val.before[0] = m->p_node->other_bv*da2[ 9][m->id_ttb[1]]*1000.;
1091    m->val.before[1] = m->p_node->other_bv*da2[11][m->id_ttb[1]]*1000.;
1092    } else {
1093      printf("BIG SHIT .... \n");
1094      exit(-10);
1095    }
1096
1097    pos = name_list_pos("monon", nl);
1098       if(nl->inform[pos] > 0) {
1099          xlimit = command_par_value("monon",cmd->clone);
1100          if(frndm() > xlimit) {
1101             m->enable = 0;
1102             printf("Monitor %s disabled\n",m->p_node->name);
1103          }
1104       }
1105    if (get_option("debug")) {
1106        printf("m-list: %d %d %s %s\n",m->id_ttb[0],m->id_ttb[1],m->p_node->name,m->p_node->base_name);
1107        printf("initial reading: %e %e\n\n",m->val.before[0],m->val.before[1]);
1108    }
1109       /*
1110       */
1111    m = m->next;
1112  };
1113  i = 0;
1114  return(i);
1115}
1116
1117static int
1118pro_correct2_getcorrs(struct in_cmd* cmd)
1119{
1120  int i;
1121  struct id_mic2 *c;  /* access to tables for monitors and correctors */
1122  // double **da1; // not used
1123  // double **da2; // not used
1124  (void)cmd;
1125/*
1126  double xlimit;
1127*/
1128  // da1 = twiss_table_beam1->d_cols; // not used
1129  // da2 = twiss_table_beam2->d_cols; // not used
1130
1131  c = correct_orbit12->cor_table;
1132  while(c) {
1133    if(c->id_ttb[0] > 0) {
1134/*     c->val.before[0] = da1[59][c->id_ttb[0]]*1000.; */
1135/*     c->val.before[1] = da1[60][c->id_ttb[0]]*1000.; */
1136       c->val.before[0] = c->p_node_s1->chkick*1000.;
1137       c->val.before[1] = c->p_node_s1->cvkick*1000.;
1138    } else if(c->id_ttb[1] > 0) {
1139/*     c->val.before[0] = da2[59][c->id_ttb[1]]*1000.; */
1140/*     c->val.before[1] = da2[60][c->id_ttb[1]]*1000.; */
1141       c->val.before[0] = c->p_node_s2->chkick*1000.;
1142       c->val.before[1] = c->p_node_s2->cvkick*1000.;
1143    }
1144    if (get_option("debug")) {
1145      printf("c-list: %d %d %s %s\n",c->id_ttb[0],c->id_ttb[1],c->p_node->name,c->p_node->base_name);
1146      printf("initial strengths: %e %e\n",c->val.before[0],c->val.before[1]);
1147    }
1148    /*
1149    */
1150
1151    c = c->next;
1152  };
1153  i = 0;
1154  return(i);
1155}
1156
1157static int
1158pro_correct2_getactive(int ip, int *nm, int *nx, int *nc, double *corvec, double *monvec,char *conm)
1159{
1160  int  imon, icor;
1161  int  imona, icora;
1162  struct id_mic2 *m, *c;
1163
1164  m = correct_orbit12->mon_table;
1165  imon = 0;
1166  imona = 0;
1167  while(m) {
1168    if (get_option("debug")) {
1169      printf("from list: %d %d %s %s\n",m->id_ttb[0],m->id_ttb[1],m->p_node->name,m->p_node->base_name);
1170      printf("orbit readings: %d %f %f\n",ip,m->val.before[0],m->val.before[1]);
1171    }
1172    if(m->enable == 1) {
1173      monvec[imon] = m->val.before[ip-1];
1174      nm[imon] = imona;
1175      imon++;
1176    }
1177    imona++;
1178    m = m->next;
1179  };
1180
1181  c = correct_orbit12->cor_table;
1182  icor = 0;
1183  icora = 0;
1184  while(c) {
1185    if (get_option("debug")) {
1186      printf("from list: %d %d %d %s %s\n",c->enable,c->id_ttb[0],c->id_ttb[1],c->p_node->name,c->p_node->base_name);
1187      printf("kicker readings: %f %f\n",c->val.before[0],c->val.before[1]);
1188    }
1189    if(c->enable == 1) {
1190      corvec[icor] = c->val.before[ip-1];
1191      nx[icor] = icora;
1192      nc[icor] = icora;
1193      strcpy(conm,c->p_node->name);
1194      conm+=16;
1195      /*          printf("nc: %d %d \n",icor,nc[icor]); */
1196      icor++;
1197    }
1198    icora++;
1199    c = c->next;
1200  };
1201  return(10000*imon + icor);
1202}
1203
1204static double*
1205pro_correct2_response_ring(int ip, int nc, int nm)
1206{
1207  int    ic, im;
1208  struct id_mic2 *m, *c;  /* access to tables for monitors and correctors */
1209
1210  double **da1;
1211  double **da2;
1212  double bx_c,by_c,pix_c,piy_c;
1213  double bx_m,by_m,pix_m,piy_m;
1214  double qx0, qy0;
1215  double respx1, respy1;
1216  double respx, respy;
1217  double *dmat;
1218  int  *imat;
1219  int    mp;
1220  int    i_zero, i_one;
1221  int icb;
1222  int    i, j;
1223
1224  setbuf(stdout,(char *)0);
1225
1226  ic = 0; im = 0;
1227  i_zero = 0; i_one = 1;
1228
1229  da1 = twiss_table_beam1->d_cols;
1230  da2 = twiss_table_beam2->d_cols;
1231
1232  dmat =  mycalloc("pro_correct2_response_ring",nc*nm,sizeof(double));
1233  imat =  mycalloc("pro_correct2_response_ring",nc*nm,sizeof(int));
1234
1235/* initialize imat: */
1236  for (i=0; i<nc; i++) {
1237     for (j=0; j<nm; j++) {
1238         setupi_(&i_zero, imat, &j, &i, &nm, &nc);
1239     }
1240  }
1241
1242  c = correct_orbit12->cor_table;
1243  ic = 0;
1244
1245  while(c) {
1246    if (get_option("debug")) {
1247    printf("corrector flag: %d\n",c->enable);
1248    }
1249    if(c->enable == 1) {
1250
1251      for(icb = 0; icb < 2; icb++) {
1252      if(c->id_ttb[icb] > 0) {
1253
1254         if(icb == 0) {
1255           correct_orbit12->qx0 = da1[5][twiss_table_beam1->curr-1];
1256           correct_orbit12->qy0 = da1[8][twiss_table_beam1->curr-1];
1257           qx0 = correct_orbit12->qx0;
1258           qy0 = correct_orbit12->qy0;
1259           if(c->id_ttb[icb] > 0) {
1260              bx_c = da1[3][c->id_ttb[icb]];
1261              by_c = da1[6][c->id_ttb[icb]];
1262              pix_c = da1[5][c->id_ttb[icb]];
1263              piy_c = da1[8][c->id_ttb[icb]];
1264           } else {
1265              bx_c = 0.0;                   
1266              by_c = 0.0;                     
1267              pix_c = 0.0;                     
1268              piy_c = 0.0;                         
1269           }
1270         } else {
1271           correct_orbit12->qx0 = da2[5][twiss_table_beam2->curr-1];
1272           correct_orbit12->qy0 = da2[8][twiss_table_beam2->curr-1];
1273           qx0 = correct_orbit12->qx0;
1274           qy0 = correct_orbit12->qy0;
1275           if(c->id_ttb[icb] > 0) {
1276              bx_c = da2[3][c->id_ttb[icb]];
1277              by_c = da2[6][c->id_ttb[icb]];
1278              pix_c = da2[5][c->id_ttb[icb]];
1279              piy_c = da2[8][c->id_ttb[icb]];
1280           } else {
1281              bx_c = 0.0;                   
1282              by_c = 0.0;                     
1283              pix_c = 0.0;                     
1284              piy_c = 0.0;                         
1285           }
1286         }
1287
1288         m = correct_orbit12->mon_table;
1289         im = 0;
1290         while(m) {
1291              if (get_option("debug")) {
1292            printf("monitor flag: %d\n",m->enable);
1293              }
1294            if(m->enable == 1) {
1295             if((m->id_ttb[icb] > 0) && (c->id_ttb[icb] > 0)) { 
1296              if(m->id_ttb[icb] > 0) { 
1297                if(icb == 0) {
1298                   mp = m->id_ttb[icb];
1299                   bx_m = da1[3][mp];
1300                   by_m = da1[6][mp];
1301                   pix_m = da1[5][mp];
1302                   piy_m = da1[8][mp];
1303                } else {
1304                   mp = m->id_ttb[icb];
1305                   bx_m = da2[3][mp];
1306                   by_m = da2[6][mp];
1307                   pix_m = da2[5][mp];
1308                   piy_m = da2[8][mp];
1309                }
1310              } else {
1311                   bx_m = 0.0;       
1312                   by_m = 0.0;       
1313                   pix_m = 0.0;       
1314                   piy_m = 0.0;         
1315              }
1316     
1317              respx = 0.0;
1318              respy = 0.0;
1319
1320/*  print Twiss parameters ... */
1321              if (get_option("debug"))  {
1322                printf("%s %d %e %e %e %e -- %s %e %e %e %e\n",
1323                c->p_node->name,icb,bx_c,by_c,pix_c,piy_c,
1324                m->p_node->name,bx_m,by_m,pix_m,piy_m);
1325              }
1326
1327              if(ip == 1) {
1328                  respx1 = cos((fabs(pix_m - pix_c)*twopi) - qx0*pi);
1329                  respx = respx1*sqrt(bx_m*bx_c)/(2.0*sin(pi*qx0));
1330                  if(icb != 0) { respx = respx; }
1331                  setup_(&respx, dmat, &im, &ic, &nm, &nc);
1332              }  else if (ip == 2) {
1333                  respy1 = cos((fabs(piy_m - piy_c)*twopi) - qy0*pi);
1334                  respy = respy1*sqrt(by_m*by_c)/(2.0*sin(pi*qy0));
1335                  if(icb != 0) { respy = respy; }
1336                  setup_(&respy, dmat, &im, &ic, &nm, &nc);
1337              }
1338              if((fabs(respy) > 0.000006) || (fabs(respx) > 0.000006)) { 
1339                 if (get_option("debug"))  {
1340                    printf("true %d %d",ic,im); 
1341                 }
1342                    setupi_(&i_one, imat, &im, &ic, &nm, &nc);
1343              } else {
1344                 if (get_option("debug"))  {
1345                    printf("false "); 
1346                 }
1347                    setupi_(&i_zero, imat, &im, &ic, &nm, &nc);
1348              }
1349              if (get_option("debug"))  {
1350                printf("Response:  %d %d %e %e %e \n",ic,im,respx, respy,fabs(respy));
1351              }
1352         }
1353              im++;
1354            }
1355            m = m->next;
1356         };
1357      }
1358      }
1359      ic++;
1360    }
1361    c = c->next;
1362  };
1363  if (get_option("debug"))  {
1364     primat_(imat,&nm, &nc);
1365     prdmat_(dmat,&nm, &nc);
1366     printf("\n");
1367     printf("\n");
1368  }
1369  return(dmat);
1370}
1371
1372static void
1373pro_correct2_write_results(double *monvec, double *resvec, double *corvec, int *nx, int *nc, int *nm, int imon, int icor, int ip)
1374{
1375/*                                              */
1376/* Writes a summary of the correction           */
1377/* Writes correctors strengths into sequences   */
1378/* Fills TFS tables for correctors and monitors */
1379/* Fills the 'stren.out' output                 */
1380/* Makes various prints on request              */
1381/*                                              */
1382  int i;
1383  int rst;
1384  double corrm;
1385  struct id_mic2 *m, *c;      /* access to tables for monitors and correctors */
1386
1387  m = correct_orbit12->mon_table;
1388  c = correct_orbit12->cor_table;
1389
1390  if(fddata != NULL) {
1391     rst = get_variable("n");
1392     fprintf(fddata,"%d %d %e %e %e %e %e %e\n",ip,rst,cprp(monvec,imon),cprp(resvec,imon),crms(monvec,imon),crms(resvec,imon),copk(monvec,imon),copk(resvec,imon));
1393  }
1394
1395  if(print_correct_opt > 0) {
1396     printf("CORRECTION SUMMARY:   \n\n");
1397     printf("rms before correction: %f mm\nrms after correction:  %f mm\n\n",crms(monvec,imon),crms(resvec,imon));
1398     printf("ptp before correction: %f mm\nptp after correction:  %f mm\n\n",cprp(monvec,imon),cprp(resvec,imon));
1399  }
1400
1401  if(print_correct_opt > 1) {
1402  printf("Monitor:  Before:     After:    Difference:\n");
1403  printf("           (mm)        (mm)         (mm)   \n");
1404  }
1405
1406  for(i=0;i<imon;i++) {
1407    if(print_correct_opt > 1) {
1408      printf("%s   %-4.3f     %-4.3f     %-4.3f\n",m[nm[i]].p_node->name,monvec[i],resvec[i],resvec[i]-monvec[i]);
1409    }
1410    m[nm[i]].val.after[ip-1] = resvec[i];
1411    pro_correct2_fill_mon_table(ip,m[nm[i]].p_node->name,monvec[i],resvec[i]);
1412  }
1413
1414  corrm = copk(corvec,icor);
1415
1416  printf("Max strength: %e should be less than %e\n",corrm,corrl);
1417  if(corrm > corrl) {
1418     printf("++++++ warning: maximum corrector strength larger than limit\n");
1419  }
1420  set_variable("corrmax",&corrm);
1421  if(print_correct_opt > 1) {
1422     printf("Max strength: %e\n",copk(corvec,icor));
1423     printf("Corrector:  Before:     After:    Difference:\n");
1424     printf("             (mrad)     (mrad)       (mrad)  \n");
1425  }
1426
1427  for(i=0;i<icor;i++) {   /* loop over all correctors */
1428
1429    c[nc[i]].val.after[ip-1] = corvec[nx[i]-1];
1430    if(print_correct_opt > 1) {
1431      printf("%s %-3.6f %-3.6f %-3.6f\n",c[nc[i]].p_node->name,
1432                                         c[nc[i]].val.before[ip-1],
1433                                         corvec[nx[i]-1]+c[nc[i]].val.before[ip-1], 
1434                                         corvec[nx[i]-1]);
1435    }
1436
1437    if(ip == 1) {
1438      /* Fill horizontal corrections for beam 1  */
1439      if(c[nc[i]].id_ttb[0] > 0) {
1440         c[nc[i]].p_node_s1->chkick += c[nc[i]].p_node_s1->other_bv*0.001*corvec[nx[i]-1];
1441         pro_correct2_fill_corr_table(0,
1442                                      ip,
1443                                      c[nc[i]].p_node->name,
1444                                      c[nc[i]].val.before[ip-1]*0.001,
1445                                      c[nc[i]].p_node_s1->chkick);
1446/*                                    c[nc[i]].p_node_s1->other_bv*0.001*corvec[nx[i]-1]); */
1447         if(fcdata != NULL) {
1448           fprintf(fcdata,"[1] %s = %e;\n",c[nc[i]].p_node->name,c[nc[i]].p_node_s1->other_bv*0.001*corvec[nx[i]-1]);
1449         }
1450      }
1451      /* Fill horizontal corrections for beam 2  */
1452      if(c[nc[i]].id_ttb[1] > 0) {
1453         c[nc[i]].p_node_s2->chkick += 0.001*corvec[nx[i]-1];
1454         pro_correct2_fill_corr_table(1,
1455                                      ip,
1456                                      c[nc[i]].p_node->name,
1457                                      c[nc[i]].val.before[ip-1]*0.001,
1458                                      c[nc[i]].p_node_s2->chkick);
1459/*                                    c[nc[i]].p_node_s2->other_bv*0.001*corvec[nx[i]-1]); */
1460         if(fcdata != NULL) {
1461           fprintf(fcdata,"[2] %s = %e;\n",c[nc[i]].p_node->name,0.001*corvec[nx[i]-1]);
1462         }
1463      }
1464
1465    } else if (ip == 2) {
1466      /* Fill vertical corrections for beam 1  */
1467      if(c[nc[i]].id_ttb[0] > 0) {
1468         c[nc[i]].p_node_s1->cvkick += c[nc[i]].p_node_s1->other_bv*0.001*corvec[nx[i]-1];
1469         pro_correct2_fill_corr_table(0,
1470                                      ip,
1471                                      c[nc[i]].p_node->name,
1472                                      c[nc[i]].val.before[ip-1]*0.001,
1473                                      c[nc[i]].p_node_s1->cvkick);
1474/*                                    c[nc[i]].p_node_s1->other_bv*0.001*corvec[nx[i]-1]); */
1475         if(fcdata != NULL) {
1476           fprintf(fcdata,"[1] %s = %e;\n",c[nc[i]].p_node->name,c[nc[i]].p_node_s1->other_bv*0.001*corvec[nx[i]-1]);
1477         }
1478      }
1479      if(c[nc[i]].id_ttb[1] > 0) {
1480      /* Fill vertical corrections for beam 2  */
1481         c[nc[i]].p_node_s2->cvkick += 0.001*corvec[nx[i]-1];
1482         pro_correct2_fill_corr_table(1,
1483                                      ip,
1484                                      c[nc[i]].p_node->name,
1485                                      c[nc[i]].val.before[ip-1]*0.001,
1486                                      c[nc[i]].p_node_s2->cvkick);
1487/*                                    c[nc[i]].p_node_s2->other_bv*0.001*corvec[nx[i]-1]); */
1488         if(fcdata != NULL) {
1489           fprintf(fcdata,"[2] %s = %e;\n",c[nc[i]].p_node->name,0.001*corvec[nx[i]-1]);
1490         }
1491      }
1492    }
1493
1494  }  /* end loop over correctors */
1495}
1496
1497static void
1498correct_correct1(struct in_cmd* cmd)
1499/* Steering routine for orbit corrections of one beam */
1500{
1501  char rout_name[] = "correct_correct";
1502  int ix, im, ip, idrop; // , it not used
1503  int j,nnnseq; // ,err not used
1504  int imon, icor;
1505  int ncorr, nmon;
1506  int niter;
1507//  int resout; // not used
1508  int twism;
1509  int dbg = 0;
1510  int ifail, sflag; // , svdflg; // not used
1511  float  rms;
1512  double sngcut, sngval;
1513  double tmp1, tmp2, tmp3, tmp4;
1514  double  sigcut;            /* number of sigmas (normalized) for filter cut */
1515  char    *clist, *mlist;    /* file names for monitor and corrector output */
1516  double  *dmat = {NULL};    /* response matrix, double precision */
1517  double  *corvec, *monvec;  /* vectors to hold measured orbit and correctors */
1518  double  *resvec;           /* vector to hold corrected orbit */
1519  char    *conm;             /* vector to hold corrector names (for MICADO) */
1520  int     *sing;             /* array to store pointer to singular correctors */
1521  static int     *nm, *nx, *nc;
1522  struct id_mic  *corl;
1523
1524  ip = pro_correct_getcommands(cmd);
1525  im = pro_correct_gettables(ip, cmd);
1526  ncorr = im%10000; nmon  = im/10000;
1527  printf("%d monitors and %d correctors found in input\n",nmon,ncorr);
1528  if(nmon == 0) {
1529    printf("No monitor found in input, no correction done\n");
1530    return;
1531  }
1532  if(ncorr == 0) {
1533    printf("No corrector found in input, no correction done\n");
1534    return;
1535  }
1536
1537
1538  /* For debugging set output buffer to zero */
1539  if (get_option("debug"))  setbuf(stdout,NULL);
1540
1541
1542  /* Prepare file descriptors for the output */
1543  if(command_par_value("resout",cmd->clone) > 0) { // (resout = not used
1544     if(fddata == NULL) {
1545        if((fddata = fopen("corr.out","w")) == NULL)
1546           exit(99);
1547     }
1548     if(fcdata == NULL) {
1549        if((fcdata = fopen("stren.out","w")) == NULL)
1550           exit(99);
1551     }
1552  }
1553
1554
1555  /* If only Twiss summary is required prepare and write it */
1556  if((twism = command_par_value("twissum",cmd->clone)) > 0) {
1557     if(ftdata == NULL) {
1558        if((ftdata = fopen("twiss.summ","w")) == NULL)
1559           exit(99);
1560     }
1561     j = 1;
1562     if((nnnseq = get_variable("n")) == 0) {
1563         nnnseq = twism;
1564     }
1565     double_from_table_row("summ", "xcomax",&j,&tmp1); // err = not used
1566     double_from_table_row("summ", "xcorms",&j,&tmp2); // err = not used
1567     double_from_table_row("summ", "ycomax",&j,&tmp3); // err = not used
1568     double_from_table_row("summ", "ycorms",&j,&tmp4); // err = not used
1569     fprintf(ftdata," T: %d %e %e %e %e\n",nnnseq,tmp1,tmp2,tmp3,tmp4);
1570     return;
1571  }
1572
1573
1574  /* allocate vectors used by correction algorithms */
1575  nx  =  mycalloc("correct_correct_nx",ncorr,sizeof(int));
1576  nc  =  mycalloc("correct_correct_nc",ncorr,sizeof(int));
1577  nm  =  mycalloc("correct_correct_nm",nmon,sizeof(int));
1578  sing =  mycalloc("correct_correct_sing",ncorr*2,sizeof(int));
1579  corvec  = mycalloc("correct_correct_corvec",ncorr,sizeof(double));
1580  monvec  = mycalloc("correct_correct_monvec",nmon,sizeof(double));
1581  resvec  = mycalloc("correct_correct_resvec",nmon,sizeof(double));
1582  conm = mycalloc("correct_correct_conm",ncorr*16,sizeof(char));
1583
1584
1585  /* get original settings of correctors from input Twiss-table */
1586  pro_correct_getcorrs(cmd); // it = not used
1587
1588  /* get input orbit, default is from input Twiss-table */
1589  /* if flag "extern" is true: can be from external table */
1590  if(command_par_value("extern",cmd->clone))
1591     pro_correct_getorbit_ext(cmd); // it = not used
1592  else
1593     pro_correct_getorbit(cmd); // it = not used
1594
1595
1596  /* find and prepare enabled correctors and monitors, may be repeated */
1597  ix = pro_correct_getactive(ip, nm, nx, nc, corvec, monvec, conm);
1598  icor = ix%10000; imon  = ix/10000;
1599  printf("%d monitors and %d correctors enabled\n",imon,icor);
1600
1601
1602  /* normalized cut on beam position, if requested */
1603  if((sigcut = command_par_value("moncut",cmd->clone)) > 0) {
1604     idrop = pro_correct_filter(ip,sigcut);
1605     printf("Disabled %d monitors with %-2.2f sigma cut\n",idrop,sigcut);
1606     ix = pro_correct_getactive(ip, nm, nx, nc, corvec, monvec, conm);
1607     icor = ix%10000; imon  = ix/10000;
1608     printf("After filter of %-2.2f sigma:\n",sigcut);
1609     printf("%d monitors and %d correctors enabled\n",imon,icor);
1610  }
1611
1612
1613  /* set up response matrix for ring or line */
1614  corl = correct_orbit->cor_table;
1615  if(strcmp("ring",command_par_string("flag",cmd->clone)) == 0) {
1616    if(dmat != NULL) myfree(rout_name,dmat);
1617    /* icor and imon used to set up correct matrix size !! */
1618    dmat  = pro_correct_response_ring(ip,icor,imon);
1619    if(command_par_value("cond",cmd->clone) == 1) { // (svdflg = not used
1620       sngcut = command_par_value("sngcut", cmd->clone);
1621       sngval = command_par_value("sngval", cmd->clone);
1622       printf("SVD conditioning requested ...\n");
1623       if (get_option("debug")) printf("Conditioning parameters: %e %e\n",sngcut, sngval);
1624
1625       /* printf("Time before svd-comd:  %-6.3f\n",fextim());    */
1626       sflag=c_svddec(dmat,imon,icor,sing,&sngcut,&sngval);
1627       printf("Initially found %d singular values\n",sflag);
1628       /* printf("Time after svd-cond:  %-6.3f\n",fextim());     */
1629       /* printf("sflag: %d\n",sflag); */
1630        printf("sflag: %d\n",sflag); 
1631       for(ix=0;ix<sflag;ix++) {
1632         corl[nx[sing[2*ix+0]]].enable = 0;
1633           printf("Removed:   %d %s\n",nx[sing[2*ix+0]],
1634                 corl[nx[sing[2*ix+0]]].p_node->name);
1635         if(dbg == 1) {
1636           printf("Removed:   %d %s\n",nx[sing[2*ix+0]],
1637                 corl[nx[sing[2*ix+0]]].p_node->name);
1638         }
1639       }
1640       ix = pro_correct_getactive(ip, nm, nx, nc, corvec, monvec, conm);
1641       icor = ix%10000; imon  = ix/10000;
1642       printf("After SVD conditioning:             \n");
1643       printf("%d monitors and %d correctors enabled\n\n",imon,icor);
1644       if(dmat != NULL) myfree(rout_name,dmat);
1645       /* icor and imon used to set up correct matrix size !! */
1646       dmat  = pro_correct_response_ring(ip,icor,imon);
1647       sflag=c_svddec(dmat,imon,icor,sing,&sngcut,&sngval);
1648       printf("Finally found %d singular values\n",sflag);
1649    }
1650  }
1651  else if(strcmp("line",command_par_string("flag",cmd->clone)) == 0) {
1652    if(dmat != NULL) myfree(rout_name,dmat);
1653          printf("make response for line\n");
1654    dmat  = pro_correct_response_line(ip,icor,imon); 
1655
1656    if(command_par_value("cond",cmd->clone) == 1) { // (svdflg = not used
1657       sngcut = command_par_value("sngcut", cmd->clone);
1658       sngval = command_par_value("sngval", cmd->clone);
1659       printf("SVD conditioning requested ...\n");
1660       if (get_option("debug")) printf("Conditioning parameters: %e %e\n",sngcut, sngval);
1661
1662       /* printf("Time before svd-comd:  %-6.3f\n",fextim());    */
1663       sflag=c_svddec(dmat,imon,icor,sing,&sngcut,&sngval);
1664       printf("Initially found %d singular values\n",sflag);
1665       /* printf("Time after svd-cond:  %-6.3f\n",fextim());     */
1666       /* printf("sflag: %d\n",sflag); */
1667       for(ix=0;ix<sflag;ix++) {
1668         corl[nx[sing[2*ix+0]]].enable = 0;
1669         if(dbg == 1) {
1670           printf("Removed:   %d %s\n",nx[sing[2*ix+0]],
1671                 corl[nx[sing[2*ix+0]]].p_node->name);
1672         }
1673       }
1674       ix = pro_correct_getactive(ip, nm, nx, nc, corvec, monvec, conm);
1675       icor = ix%10000; imon  = ix/10000;
1676       printf("After SVD conditioning:             \n");
1677       printf("%d monitors and %d correctors enabled\n\n",imon,icor);
1678       if(dmat != NULL) myfree(rout_name,dmat);
1679       /* icor and imon used to set up correct matrix size !! */
1680       dmat  = pro_correct_response_ring(ip,icor,imon);
1681       sflag=c_svddec(dmat,imon,icor,sing,&sngcut,&sngval);
1682       printf("Finally found %d singular values\n",sflag);
1683    }
1684  }
1685 
1686  else { printf("INVALID MACHINE TYPE\n"); exit(-1);
1687  }
1688
1689  if (get_option("debug")) {
1690    pro_correct_prtwiss();
1691    pro_correct_write_cocu_table();
1692  }
1693
1694
1695  /* LSQ correction, use all available correctors */
1696  if(strcmp("lsq",command_par_string("mode",cmd->clone)) == 0) {
1697    /*frs haveit_(dmat,monvec,corvec,resvec,nx,&imon,&icor); */
1698    /* printf("Time before lsq:  %-6.3f\n",fextim());   */
1699    c_haveit(dmat,monvec,corvec,resvec,nx,imon,icor);
1700    /* printf("Time after lsq:  %-6.3f\n",fextim());    */
1701    pro_correct_write_results(monvec, resvec, corvec, nx, nc, nm, imon, icor, ip);
1702  }
1703
1704
1705  /* SVD correction, use all available correctors */
1706  if(strcmp("svd",command_par_string("mode",cmd->clone)) == 0) {
1707    /*frs haveit_(dmat,monvec,corvec,resvec,nx,&imon,&icor); */
1708    /* printf("Time before svd-corr:  %-6.3f\n",fextim());   */
1709       sflag=c_svdcorr(dmat,monvec,corvec,resvec,nx,imon,icor);
1710    /* printf("Time after svd-corr:  %-6.3f\n",fextim());    */
1711    pro_correct_write_results(monvec, resvec, corvec, nx, nc, nm, imon, icor, ip);
1712  }
1713
1714
1715  /* MICADO correction, get desired number of correctors from command */
1716  corrl = command_par_value("corrlim",cmd->clone);
1717  set_variable("corrlim",&corrl);
1718  if(strcmp("micado",command_par_string("mode",cmd->clone)) == 0) {
1719    printf("enter MICADO correction ...\n");
1720    if((niter = command_par_value("ncorr",cmd->clone)) == 0) {
1721          printf("Requested %d correctors (\?\?\?) set to %d\n",niter,icor);
1722          niter = icor;
1723    }
1724    else if((niter = command_par_value("ncorr",cmd->clone)) < 0) {
1725          printf("Requested %d correctors (\?\?\?) set to 0\n",niter);
1726          niter = 0;
1727    }
1728    else if((niter = command_par_value("ncorr",cmd->clone)) > icor) {
1729          printf("Fewer correctors available than requested by ncorr\n");
1730          printf("you want %d,  you get %d\n",niter,icor);
1731          printf("ncorr reset to %d\n",icor);
1732          niter = icor;
1733    }
1734    rms  = 1000.0*command_par_value("error",cmd->clone);
1735    /*frs       micit_(dmat,monvec,corvec,resvec,nx,&rms,&imon,&icor,&niter); */
1736    /* printf("Time before micado:  %-6.3f\n",fextim());  */
1737    ifail = c_micit(dmat,conm,monvec,corvec,resvec,nx,rms,imon,icor,niter);
1738    /* printf("Time after micado:  %-6.3f\n",fextim());   */
1739
1740    if(ifail != 0) {
1741       printf("MICADO correction completed with error code %d\n\n",ifail);
1742       warning("MICADO back with error",", no correction done");
1743    }
1744if(ifail == 0) {
1745    pro_correct_write_results(monvec, resvec, corvec, nx, nc, nm, imon, icor, ip);
1746 }
1747  }
1748
1749
1750  /* write corrector output to tfs table */
1751  if ((clist = command_par_string("clist",cmd->clone)) != NULL) {
1752    out_table("corr",corr_table,clist);
1753  }
1754  /* write monitor output to tfs table */
1755  if ((mlist = command_par_string("mlist",cmd->clone)) != NULL) {
1756    out_table("mon",mon_table,mlist);
1757  }
1758
1759
1760  /* Clean up at the end of the module */
1761  myfree(rout_name,nm);myfree(rout_name,dmat);myfree(rout_name,nx);
1762  myfree(rout_name,nc);myfree(rout_name,corvec);
1763  myfree(rout_name,monvec);myfree(rout_name,resvec); myfree(rout_name,conm);
1764  return;
1765}
1766
1767static void
1768correct_correct(struct in_cmd* cmd)
1769/* Steering routine for orbit corrections */
1770{
1771/*
1772  char rout_name[] = "correct_correct";
1773*/
1774  char  *orbtab1, *orbtab2;
1775 
1776
1777  /* Call for one or two ring orbit correction */
1778  if(command_par_value("tworing",cmd->clone)) {
1779           if ((orbtab1 = command_par_string("beam1tab",cmd->clone)) == NULL) {
1780              fatal_error("Two beam correction requested but no table supplied for beam 1",orbtab1);
1781           }
1782           if ((orbtab2 = command_par_string("beam2tab",cmd->clone)) == NULL) {
1783              fatal_error("Two beam correction requested but no table supplied for beam 2",orbtab2);
1784           }
1785           printf("Want to use orbits from: %s and : %s\n",orbtab1,orbtab2);
1786           correct_correct2(cmd);
1787  }  else {
1788           printf("Want to correct orbit of a single ring\n");                   
1789           if ((orbtab1 = command_par_string("beam1tab",cmd->clone)) != NULL) {
1790              warning(" "," ");
1791              warning("Single beam correction requested but beam 1 table supplied:",orbtab1);
1792              warning("Requested table ignored:",orbtab1);
1793              warning(" "," ");
1794           }
1795           if ((orbtab2 = command_par_string("beam2tab",cmd->clone)) != NULL) {
1796              warning(" "," ");
1797              warning("Single beam correction requested but beam 2 table supplied:",orbtab2);
1798              warning("Requested table ignored:",orbtab2);
1799              warning(" "," ");
1800           }
1801           correct_correct1(cmd);
1802  }
1803}
1804
1805#if 0 // not used...
1806static void
1807pro_correct_option(struct in_cmd* cmd)
1808{
1809  struct name_list* nl = cmd->clone->par_names;
1810  int i, debug;
1811  int val, pos, seed;
1812
1813  if ((debug=get_option("debug")))  {
1814     fprintf(prt_file, "in coption routine\n");
1815     for(i=0;i<cmd->tok_list->curr;i++) {
1816        fprintf(prt_file, "command(s): %s\n",cmd->tok_list->p[i]);
1817     }
1818  }
1819  if ((pos = name_list_pos("seed", nl)) > -1)
1820    {
1821     if (nl->inform[pos])
1822       {
1823        seed = command_par_value("seed", cmd->clone);
1824        init55(seed);
1825       }
1826    }
1827 val = command_par_value("print", cmd->clone);
1828 if(val == 0) {
1829      if (debug)  fprintf(prt_file, "print option not set\n");
1830      print_correct_opt = 0;
1831 }else {
1832      if (debug)  fprintf(prt_file, "print option set\n");
1833      print_correct_opt = val;
1834 }
1835}
1836#endif
1837
1838static int
1839pro_correct_getcommands(struct in_cmd* cmd)
1840{
1841
1842  static char att[10][8] = {"iterate","plane","ncorr","error","clist",
1843                            "mlist", "flag","mode","",""};
1844
1845  static int iplane = 1;
1846  char plane[20];
1847
1848  if (get_option("debug")) printf("enter CORRECT module\n");
1849  if (current_sequ == NULL || current_sequ->ex_start == NULL)
1850    {
1851      warning("CORRECT, but no active sequence:", "ignored");
1852      return(-1);
1853    }
1854
1855  strcpy(plane,command_par_string(att[1],cmd->clone));
1856  if(strcmp("x",plane) == 0) {
1857    iplane = 1;
1858  } else if (strcmp("y",plane) == 0) {
1859    iplane = 2;
1860  } else if (strcmp("h",plane) == 0) {
1861    iplane = 1;
1862  } else if (strcmp("v",plane) == 0) {
1863    iplane = 2;
1864  } else {
1865    printf("No valid plane specified, x plane used \n");
1866    iplane = 1;
1867  }
1868
1869  return(iplane);
1870}
1871
1872static int
1873pro_correct_gettables(int iplane, struct in_cmd* cmd)
1874{
1875
1876  char rout_name[] = "pro_correct_gettables";
1877
1878  struct id_mic *cor_l;
1879  struct id_mic *mon_l;
1880
1881  struct table *ttb;
1882
1883  char* orbtab;
1884  char* tartab;
1885  char* modtab;
1886
1887  int j;
1888  int pps, ppt;
1889//  int set0; // not used
1890  int cntm = {0};
1891  int cntc = {0};
1892
1893  double ounits;
1894
1895/*
1896  static char atm[6][4] = {"hmon","vmon","moni","hkic","vkic","kick"};
1897*/
1898
1899/* Get access to tables, for orbit and model the default is twiss_table */
1900
1901  if ((orbtab = command_par_string("orbit",cmd->clone)) != NULL) {
1902      printf("Want to use orbit from: %s\n",orbtab);
1903    if ((pps = name_list_pos(orbtab, table_register->names)) > -1) {
1904       orbin_table = table_register->tables[pps];
1905    } else {
1906       fatal_error("ORBIT table for correction requested, but not provided:",orbtab);
1907    }
1908  } else {
1909       if((orbin_table = twiss_table) == NULL) {
1910         printf("FATAL ERROR:\n");
1911         printf("You request the ORBIT from a non-existing TWISS table\n");
1912         printf("You MUST run TWISS before trying to correct the orbit\n");
1913         printf("MAD-X stops\n");
1914         exit(81);
1915       } else {
1916         if (get_option("debug")) {
1917            printf("TWISS table: %p\n", (void*)twiss_table);
1918         }
1919       }
1920       pps = -1;
1921  }
1922
1923
1924  if ((tartab = command_par_string("target",cmd->clone)) != NULL) {
1925       printf("Want to use target orbit from: %s\n",tartab);
1926    if ((ppt = name_list_pos(tartab, table_register->names)) > -1) {
1927       target_table = table_register->tables[ppt];
1928    } else {
1929       fatal_error("TARGET table for correction requested, but not provided:",tartab);
1930    }
1931  } else {
1932       if (get_option("debug")) {
1933         printf("No target orbit requested\n");
1934       }
1935       ppt = -1;
1936  }
1937
1938  if ((modtab = command_par_string("model",cmd->clone)) != NULL) {
1939       printf("Want to use model orbit from: %s\n",modtab);
1940    if ((ppt = name_list_pos(modtab, table_register->names)) > -1) {
1941       model_table = table_register->tables[ppt];
1942    } else {
1943       fatal_error("MODEL table for correction requested, but not provided:",modtab);
1944    }
1945  } else {
1946       if((model_table = twiss_table) == NULL) {
1947         printf("FATAL ERROR:\n");
1948         printf("You request the MODEL from a non-existing TWISS table\n");
1949         printf("You MUST run TWISS before trying to correct the orbit\n");
1950         printf("MAD-X stops\n");
1951         exit(81);
1952       } else {
1953         if (get_option("debug")) {
1954            printf("TWISS table: %p\n", (void*)twiss_table);
1955         }
1956       }
1957       ppt = -1;
1958  }
1959
1960
1961
1962       if (get_option("debug")) {
1963            printf("The tables are: %p %p %p %p\n",
1964                    (void*)orbin_table, (void*)twiss_table,
1965                    (void*)target_table, (void*)model_table);
1966       }
1967       if (get_option("debug")) {
1968       }
1969
1970  if(correct_orbit == NULL) {
1971    correct_orbit = mycalloc("pro_correct_gettables",1, sizeof(struct orb_cor));
1972  }
1973       if (get_option("debug")) {
1974           printf("-0-\n");
1975       }
1976
1977/*    if(corr_table == NULL) {   */
1978    corr_table = make_table("corr", "corr", corr_table_cols,
1979            corr_table_types, 5000);
1980       if (get_option("debug")) {
1981           printf("-01-\n");
1982       }
1983    add_to_table_list(corr_table, table_register);
1984       if (get_option("debug")) {
1985           printf("-02-\n");
1986       }
1987    pro_correct_make_corr_table();
1988       if (get_option("debug")) {
1989           printf("-03-\n");
1990       }
1991/*    }                          */
1992       if (get_option("debug")) {
1993           printf("-1-\n");
1994       }
1995
1996/*    if(mon_table == NULL) { */
1997    mon_table = make_table("mon", "mon", mon_table_cols,
1998           mon_table_types, 5000);
1999       if (get_option("debug")) {
2000           printf("-11-\n");
2001       }
2002    add_to_table_list(mon_table, table_register);
2003       if (get_option("debug")) {
2004           printf("-12-\n");
2005       }
2006    pro_correct_make_mon_table();
2007       if (get_option("debug")) {
2008           printf("-13-\n");
2009       }
2010/*    }                       */
2011       if (get_option("debug")) {
2012           printf("-2-\n");
2013       }
2014
2015
2016  if(correct_orbit->cor_table != NULL) myfree(rout_name,correct_orbit->cor_table);
2017  if(correct_orbit->mon_table != NULL) myfree(rout_name,correct_orbit->mon_table);
2018  correct_orbit->cor_table = (struct id_mic *)mycalloc("pro_correct_gettables_cor",5200, sizeof(struct id_mic));
2019  correct_orbit->mon_table = (struct id_mic *)mycalloc("pro_correct_gettables_mon",5200, sizeof(struct id_mic));
2020
2021/* orbit table available, get units, if defined */
2022   if((ounits = command_par_value("units",cmd->clone)) > 0) {
2023          correct_orbit->units=ounits;                     
2024   } else {
2025          correct_orbit->units=1.0;     
2026   }
2027
2028  ttb = model_table;
2029  correct_orbit->mon_table->previous = NULL;
2030  correct_orbit->mon_table->next = NULL;
2031  correct_orbit->cor_table->previous = NULL;
2032  correct_orbit->cor_table->next = NULL;
2033       if (get_option("debug")) {
2034           printf("-3-\n");
2035       }
2036
2037  mon_l = correct_orbit->mon_table;
2038  cor_l = correct_orbit->cor_table;
2039  for (j=0; j < ttb->curr; j++) {
2040   if((strncmp(atm[iplane-1],ttb->p_nodes[j]->base_name,4) == 0) ||
2041      (strncmp(atm[2],      ttb->p_nodes[j]->base_name,4) == 0))  {
2042      mon_l->id_ttb = j;
2043      mon_l->enable = ttb->p_nodes[j]->enable;
2044      mon_l->p_node = ttb->p_nodes[j];
2045      mon_l->next = mon_l;
2046      mon_l->next++; mon_l++;
2047      cntm++;
2048    }
2049    if((strncmp(atc[iplane-1],ttb->p_nodes[j]->base_name,4) == 0) ||
2050       (strncmp(atc[2],       ttb->p_nodes[j]->base_name,4) == 0))  {
2051      cor_l->id_ttb = j;
2052      cor_l->enable = ttb->p_nodes[j]->enable;
2053      cor_l->p_node = ttb->p_nodes[j];
2054      if(command_par_value("corzero",cmd->clone) > 0) { // (set0 = not used
2055        if(iplane == 1) cor_l->p_node->chkick = 0.0;
2056        if(iplane == 2) cor_l->p_node->cvkick = 0.0;
2057      }
2058      cor_l->next = cor_l;
2059      cor_l->next++; cor_l++;
2060      cntc++;
2061    }
2062  }
2063       if (get_option("debug")) {
2064           printf("-4-\n");
2065       }
2066  mon_l--; mon_l->next = NULL;
2067  cor_l--; cor_l->next = NULL;
2068       if (get_option("debug")) {
2069         printf("done: %d %d\n",cntm,cntc);
2070       }
2071
2072  return(10000*cntm + cntc);
2073}
2074
2075static int
2076pro_correct_getorbit(struct in_cmd* cmd)
2077{
2078  struct name_list* nl;
2079  int i;
2080//  char *tartab; // not used
2081  int pos;
2082  struct id_mic *m;  /* access to tables for monitors and correctors */
2083  struct table *ttb;
2084  struct table *tar = NULL;
2085  double **da1;
2086  double **da2 = NULL;
2087  double xlimit;
2088  double rx, ry, dpsi;
2089
2090  char   strx[40];
2091  char   stry[40];
2092
2093  int    posx, posy, pospx, pospy;
2094  int    tosx = -1;
2095  int    tosy = -1;
2096  int    tospx, tospy;
2097
2098  ttb = orbin_table;
2099  da1 = ttb->d_cols;
2100
2101  if(target_table != NULL) {
2102     tar = target_table;
2103     da2 = tar->d_cols;
2104  }
2105
2106  nl = cmd->clone->par_names;
2107
2108  m = correct_orbit->mon_table;
2109
2110  strcpy(strx,"x");
2111  strcpy(stry,"y");
2112                                                                                                         
2113  if((posx = name_list_pos(strx,ttb->columns)) < 0) { 
2114      fatal_error("orbit x not found in input table",", MAD-X terminates ");
2115  }
2116  if((posy = name_list_pos(stry,ttb->columns)) < 0) { 
2117      fatal_error("orbit y not found in input table",", MAD-X terminates ");
2118  }
2119  if (get_option("debug")) {
2120    if((pospx = name_list_pos("px",ttb->columns)) < 0) { 
2121        fatal_error("orbit px not found in input table",", MAD-X terminates ");
2122    }
2123    if((pospy = name_list_pos("py",ttb->columns)) < 0) { 
2124        fatal_error("orbit py not found in input table",", MAD-X terminates ");
2125    }
2126    printf("====c1===>  %d %d %d %d \n",posx,posy,pospx,pospy);
2127  }
2128
2129  if (command_par_string("target",cmd->clone) != NULL) { // (tartab = not used
2130    if((tosx = name_list_pos("x",tar->columns)) < 0) { 
2131        fatal_error("target orbit x not found in table",", MAD-X terminates ");
2132    }
2133    if((tosy = name_list_pos("y",tar->columns)) < 0) { 
2134        fatal_error("target orbit y not found in table",", MAD-X terminates ");
2135    }
2136    if (get_option("debug")) {
2137      if((tospx = name_list_pos("px",tar->columns)) < 0) { 
2138          fatal_error("target orbit px not found in table",", MAD-X terminates ");
2139      }
2140      if((tospy = name_list_pos("py",tar->columns)) < 0) { 
2141          fatal_error("target orbit px not found in table",", MAD-X terminates ");
2142      }
2143      printf("====c1===>  %d %d %d %d \n",tosx,tosy,tospx,tospy);
2144    }
2145  }
2146
2147
2148  while(m) {
2149
2150/* If correction to target orbit, subtract the wanted orbit ... */
2151  if (command_par_string("target",cmd->clone) != NULL) { // (tartab = not used
2152    m->val.before[0] =  da1[posx][m->id_ttb] - da2[tosx][m->id_ttb];
2153    m->val.before[1] =  da1[posy][m->id_ttb] - da2[tosy][m->id_ttb];
2154    m->val.before[0] = (da1[posx][m->id_ttb] - da2[tosx][m->id_ttb])*1000.*correct_orbit->units;
2155    m->val.before[1] = (da1[posy][m->id_ttb] - da2[tosy][m->id_ttb])*1000.*correct_orbit->units;
2156  } else {
2157    m->val.before[0] = da1[posx][m->id_ttb];
2158    m->val.before[1] = da1[posy][m->id_ttb];
2159    m->val.before[0] = da1[posx][m->id_ttb]*1000.*correct_orbit->units;
2160    m->val.before[1] = da1[posy][m->id_ttb]*1000.*correct_orbit->units;
2161  }
2162
2163    pos = name_list_pos("monon", nl);
2164    if(nl->inform[pos] > 0) {
2165          xlimit = command_par_value("monon",cmd->clone);
2166          if(frndm() > xlimit) {
2167             m->enable = 0;
2168             printf("Monitor %s disabled\n",m->p_node->name);
2169          }
2170    }
2171
2172    /* scaling error should come first, monitor alignment not scaled ... */
2173    pos = name_list_pos("monscale", nl);
2174    if(nl->inform[pos] > 0) {
2175      if((command_par_value("monscale",cmd->clone)) == 1) {
2176        if(m->p_node->p_al_err != NULL) {
2177          if (get_option("debug")) {
2178           printf("m-list: %d %s %s\n",m->id_ttb,m->p_node->name,m->p_node->base_name);
2179           printf("scales: %e %e\n",m->p_node->p_al_err->a[12],m->p_node->p_al_err->a[13]);
2180          }
2181          m->val.before[0] = m->val.before[0]*(1.0 + m->p_node->p_al_err->a[12]);
2182          m->val.before[1] = m->val.before[1]*(1.0 + m->p_node->p_al_err->a[13]);
2183        }
2184      }
2185    }
2186
2187    /* monitor misalignment after all other reading manipulations ! */
2188    pos = name_list_pos("monerror", nl);
2189    if(nl->inform[pos] > 0) {
2190      if((command_par_value("monerror",cmd->clone)) == 1) {
2191        if(m->p_node->p_al_err != NULL) {
2192          if (get_option("debug")) {
2193            printf("m-list: %d %s %s\n",m->id_ttb,m->p_node->name,m->p_node->base_name);
2194            printf("errors: %e %e \n",m->p_node->p_al_err->a[6], m->p_node->p_al_err->a[7]);
2195          }
2196          dpsi = m->p_node->p_al_err->a[5];
2197          rx = m->val.before[0];
2198          ry = m->val.before[1];
2199          printf("\nA: %e %e %e\n",m->val.before[0],m->val.before[1],dpsi);
2200          m->val.before[0] =  rx * cos(dpsi) +  ry * sin(dpsi);
2201          m->val.before[1] = -rx * sin(dpsi) +  ry * cos(dpsi);
2202          printf("B: %e %e %e\n",m->val.before[0],m->val.before[1],dpsi);
2203          m->val.before[0] += m->p_node->p_al_err->a[6]*1000.;
2204          m->val.before[1] += m->p_node->p_al_err->a[7]*1000.;
2205          printf("C: %e %e %e\n",m->val.before[0],m->val.before[1],dpsi);
2206        }
2207      }
2208    }
2209    m = m->next;
2210  };
2211  i = 0;
2212  return(i);
2213}
2214
2215static int
2216pro_correct_getorbit_ext(struct in_cmd* cmd)
2217{
2218  struct name_list* nl;
2219  int i;
2220  int j;
2221//  char *tartab; // not used
2222  int pos;
2223  struct id_mic *m;  /* access to tables for monitors and correctors */
2224  struct table *ttb;
2225  struct table *tar = NULL;
2226  double **da1;
2227  double **da2 = NULL;
2228  double xlimit;
2229  char     name[NAME_L];
2230  char   l1name[NAME_L];
2231  char   l2name[NAME_L];
2232  char   l3name[NAME_L];
2233  char   l4name[NAME_L];
2234  double rx, ry, dpsi;
2235
2236  char   *nam_col;
2237  char   *x_col;
2238  char   *y_col;
2239
2240  char   strx[40];
2241  char   stry[40];
2242  char   strn[40];
2243
2244  int    posx, posy, pospx, pospy;
2245  int    tosx = -1;
2246  int    tosy = -1;
2247  int    tospx, tospy;
2248
2249  int dbk = 0;
2250  int yok;
2251
2252  int    jjx, jjy, jj;
2253
2254  ttb = orbin_table;
2255  da1 = ttb->d_cols;
2256
2257  if(target_table != NULL) {
2258     tar = target_table;
2259     da2 = tar->d_cols;
2260  }
2261
2262  nl = cmd->clone->par_names;
2263
2264  m = correct_orbit->mon_table;
2265
2266  if ((x_col = command_par_string("x_col",cmd->clone)) != NULL) {
2267         printf("X orbit in column: %s\n",x_col);
2268         strcpy(strx,x_col);
2269  } else {
2270         strcpy(strx,"x");
2271  }
2272  if ((y_col = command_par_string("y_col",cmd->clone)) != NULL) {
2273         printf("y orbit in column: %s\n",y_col);
2274         strcpy(stry,y_col);
2275  } else {
2276         strcpy(stry,"y");
2277  }
2278  if ((nam_col = command_par_string("name_col",cmd->clone)) != NULL) {
2279         printf("names in column: %s\n",nam_col);
2280         strcpy(strn,"name");
2281  } else {
2282         strcpy(strn,"name");
2283  }
2284 
2285  if((posx = name_list_pos(strx,ttb->columns)) < 0) { 
2286      fatal_error("orbit x not found in input table",", MAD-X terminates ");
2287  }
2288  if((posy = name_list_pos(stry,ttb->columns)) < 0) { 
2289      fatal_error("orbit y not found in input table",", MAD-X terminates ");
2290  }
2291  if (get_option("debug")) {
2292    if((pospx = name_list_pos("px",ttb->columns)) < 0) { 
2293        warning("orbit px not found in input table",", MAD-X continues ");
2294    }
2295    if((pospy = name_list_pos("py",ttb->columns)) < 0) { 
2296        warning("orbit py not found in input table",", MAD-X continues ");
2297    }
2298    printf("====c1===>  %d %d %d %d \n",posx,posy,pospx,pospy);
2299  }
2300
2301  if (command_par_string("target",cmd->clone) != NULL) { // (tartab = not used
2302    if((tosx = name_list_pos("x",tar->columns)) < 0) { 
2303        fatal_error("target orbit x not found in table",", MAD-X terminates ");
2304    }
2305    if((tosy = name_list_pos("y",tar->columns)) < 0) { 
2306        fatal_error("target orbit y not found in table",", MAD-X terminates ");
2307    }
2308    if (get_option("debug")) {
2309      if((tospx = name_list_pos("px",tar->columns)) < 0) { 
2310          warning("target orbit px not found in table",", MAD-X continues ");
2311      }
2312      if((tospy = name_list_pos("py",tar->columns)) < 0) { 
2313          warning("target orbit px not found in table",", MAD-X continues ");
2314      }
2315      printf("====c1===>  %d %d %d %d \n",tosx,tosy,tospx,tospy);
2316    }
2317  }
2318
2319
2320  if (get_option("debug")) {
2321    printf("Number in table: %d\n",ttb->curr);
2322    for (j=1; j < (ttb->curr)+1; j++) {
2323      i = string_from_table_row(ttb->name, "name", &j, name);
2324    }
2325  }
2326
2327  jj = 0;
2328  while(m) {
2329
2330
2331    strcpy(l1name,m->p_node->name);
2332    stolower(l1name);
2333    strcpy(l2name,strip(l1name));
2334    supp_tb(l2name);
2335
2336     if (dbk ==1) printf("monitor name: %s\n",l2name);
2337
2338    jjx = -1;
2339    jjy = -1;
2340    jj++;
2341    yok = 0;
2342
2343    for (j=1; j < (ttb->curr)+1; j++) {
2344      i = string_from_table_row(ttb->name, "name", &j, name);
2345      strcpy(l3name,name);
2346      stolower(l3name);
2347      strcpy(l4name,strip(l3name));
2348      supp_tb(l4name);
2349      if(strlen(l4name) == strlen(l2name)) {
2350         if(strncmp(l4name,l2name,strlen(l2name)) == 0) {
2351              jjx = j-1;
2352              jjy = jj-1;
2353              yok = 1;
2354           if(dbk ==1) printf("monitor names found: %s %s %d\n",l2name,l4name,yok);
2355             }
2356      }
2357    }
2358     if(dbk ==1) printf("jjx,jjy %d %d\n",jjx,jjy);
2359
2360/* If correction to target orbit, subtract the wanted orbit ... */
2361
2362    if((jjy >= 0) && (yok == 1)) { 
2363/*  if(jjx >= 0)  {  */
2364     if (command_par_string("target",cmd->clone) != NULL) { // (tartab = not used
2365  if(dbk == 1) {
2366       printf("x ==> %d %d %e %e\n",jjx,m->id_ttb,da1[posx][jjx],da2[tosx][jjy]);
2367       printf("y ==> %e %e\n",da1[posy][jjx],da2[tosy][jjy]);
2368  }
2369       m->val.before[0] =  da1[posx][jjx] - da2[tosx][jjy];
2370       m->val.before[1] =  da1[posy][jjx] - da2[tosy][jjy];
2371       m->val.before[0] = (da1[posx][jjx] - da2[tosx][jjy])*1000.*correct_orbit->units;
2372       m->val.before[1] = (da1[posy][jjx] - da2[tosy][jjy])*1000.*correct_orbit->units;
2373  if(dbk == 1) {
2374       printf("bxy ==> %s %d %e %e\n",m->p_node->name,jjx,m->val.before[0],m->val.before[1]);
2375  }
2376     } else {
2377  if(dbk == 1) {
2378       printf("x ==> %e %e\n",da1[posx][jjx],da2[tosx][jjx]);
2379       printf("y ==> %e %e\n",da1[posy][jjx],da2[tosy][jjx]);
2380  }
2381       m->val.before[0] = da1[posx][jjx];
2382       m->val.before[1] = da1[posy][jjx];
2383       m->val.before[0] = da1[posx][jjx]*1000.*correct_orbit->units;
2384       m->val.before[1] = da1[posy][jjx]*1000.*correct_orbit->units;
2385if(dbk == 1) {
2386       printf("bxy ==> %s %d %e %e\n",m->p_node->name,jjx,m->val.before[0],m->val.before[1]);
2387}
2388     }
2389
2390    pos = name_list_pos("monon", nl);
2391    if(nl->inform[pos] > 0) {
2392          xlimit = command_par_value("monon",cmd->clone);
2393          if(frndm() > xlimit) {
2394             m->enable = 0;
2395             printf("Monitor %s disabled\n",m->p_node->name);
2396          }
2397    }
2398
2399    /* scaling error should come first, monitor alignment not scaled ... */
2400    pos = name_list_pos("monscale", nl);
2401    if(nl->inform[pos] > 0) {
2402      if((command_par_value("monscale",cmd->clone)) == 1) {
2403        if(m->p_node->p_al_err != NULL) {
2404          if (get_option("debug")) {
2405           printf("m-list: %d %s %s\n",m->id_ttb,m->p_node->name,m->p_node->base_name);
2406           printf("scales: %e %e\n",m->p_node->p_al_err->a[12],m->p_node->p_al_err->a[13]);
2407          }
2408          m->val.before[0] = m->val.before[0]*(1.0 + m->p_node->p_al_err->a[12]);
2409          m->val.before[1] = m->val.before[1]*(1.0 + m->p_node->p_al_err->a[13]);
2410        }
2411      }
2412    }
2413
2414    /* monitor misalignment after all other reading manipulations ! */
2415    pos = name_list_pos("monerror", nl);
2416     if(nl->inform[pos] > 0) {
2417      if((command_par_value("monerror",cmd->clone)) == 1) {
2418        if(m->p_node->p_al_err != NULL) {
2419          if (get_option("debug")) {
2420            printf("m-list: %d %s %s\n",m->id_ttb,m->p_node->name,m->p_node->base_name);
2421            printf("errors: %e %e \n",m->p_node->p_al_err->a[6], m->p_node->p_al_err->a[7]);
2422          }
2423          dpsi = m->p_node->p_al_err->a[5];
2424          rx = m->val.before[0];
2425          ry = m->val.before[1];
2426          printf("\nA: %e %e %e\n",m->val.before[0],m->val.before[1],dpsi);
2427          m->val.before[0] =  rx * cos(dpsi) +  ry * sin(dpsi);
2428          m->val.before[1] = -rx * sin(dpsi) +  ry * cos(dpsi);
2429          printf("B: %e %e %e\n",m->val.before[0],m->val.before[1],dpsi);
2430          m->val.before[0] += m->p_node->p_al_err->a[6]*1000.;
2431          m->val.before[1] += m->p_node->p_al_err->a[7]*1000.;
2432          printf("C: %e %e %e\n",m->val.before[0],m->val.before[1],dpsi);
2433        }
2434      }
2435     }
2436    } else { m->enable = 0;} /* Only enable monitors found in input */
2437    m = m->next;
2438  };
2439  i = 0;
2440  return(i);
2441}
2442
2443static int
2444pro_correct_getcorrs(struct in_cmd* cmd)
2445{
2446  int i;
2447  struct id_mic *c;  /* access to tables for monitors and correctors */
2448//  struct table *ttb; // not used
2449//  double **da1; // not used
2450
2451  (void)cmd;
2452//  ttb = model_table; // not used
2453
2454  // da1 = ttb->d_cols; // not used
2455
2456  c = correct_orbit->cor_table;
2457  while(c) {
2458/*  c->val.before[0] = da1[59][c->id_ttb]*1000.;   */
2459/*  c->val.before[1] = da1[60][c->id_ttb]*1000.;   */
2460    c->val.before[0] = c->p_node->chkick*1000.;
2461    c->val.before[1] = c->p_node->cvkick*1000.;
2462
2463    /*
2464    printf("c-list: %d %s %s\n",c->id_ttb,c->p_node->name,c->p_node->base_name);
2465    printf("initial strengths: %e %e\n",c->val.before[0],c->val.before[1]);
2466    */
2467
2468    c = c->next;
2469  };
2470  i = 0;
2471  return(i);
2472}
2473
2474static void
2475pro_correct_prtwiss(void)
2476{
2477  int i,j;
2478  int pr_cols;
2479  struct table *ttb;
2480  double **da1;
2481
2482  ttb = model_table;
2483
2484  printf(" %d %d\n",ttb->curr,ttb->num_cols);
2485  for (i=0; i<ttb->curr; i++) {
2486
2487    printf(" %s %s\n",ttb->s_cols[0][i],ttb->s_cols[1][i]);
2488  }
2489
2490  da1 = ttb->d_cols;
2491  for (j=0; j < ttb->curr; j++) {
2492    printf("\n\n");
2493    printf("from table: %s \n",ttb->node_nm->p[j]);
2494    printf("from node:  %s \n",ttb->p_nodes[j]->name);
2495    printf(" %s %s\n",ttb->s_cols[0][j],ttb->s_cols[1][j]);
2496
2497    pr_cols = ttb->num_cols;
2498    pr_cols = 19; /* print only for 20 columns */
2499    for (i=0; i<pr_cols; i++) {
2500      if(&da1[i][0] != NULL) {
2501      printf("%-8s %f\n",twiss_table_cols[i],da1[i][j]);
2502      }
2503    }
2504  }
2505  return;
2506}
2507
2508static void
2509pro_correct_write_cocu_table(void)
2510{
2511  int i,j;
2512  int pr_cols;
2513  int cp[13] = {1,0,2,9,11,3,6,4,7,5,8,15,17};
2514  struct table *ttb;
2515  double **da1;
2516  FILE *fp1;
2517
2518  fp1 = fopen ("cocu_in.opt","w");
2519  ttb = model_table;
2520
2521  pr_cols = ttb->num_cols;
2522  pr_cols = 13; /* print only for 19 columns */
2523  fprintf(fp1,"*");
2524  for (i=0; i<pr_cols; i++) {
2525    fprintf(fp1,"%-8s ",twiss_table_cols[cp[i]]);
2526  }
2527
2528  da1 = ttb->d_cols;
2529  for (j=0; j < ttb->curr; j++) {
2530
2531    fprintf(fp1,"\n%s %s ",ttb->s_cols[1][j],ttb->s_cols[0][j]);
2532    for (i=2; i<pr_cols; i++) {
2533      if(&da1[cp[i]][0] != NULL) {
2534      fprintf(fp1," %f",da1[cp[i]][j]);
2535      }
2536    }
2537  }
2538  return;
2539}
2540
2541static int
2542pro_correct_filter(int iplane, double sigcut)
2543{
2544  int    im, ip, icnt; // ic, no used
2545  struct id_mic *m;  /* access to tables for monitors */
2546
2547  struct table *ttb;
2548  static char  pl[2] = "xy";
2549  double **da1;
2550  double bx_m = -9999.;
2551  double xsig;
2552  double xmea; //, ymea; not used
2553  double xn;
2554
2555  ttb = model_table;
2556  da1 = ttb->d_cols;
2557  im = 0; icnt = 0; // ic = 0; // not used
2558  ip = iplane - 1;
2559
2560  printf("A (normalized) cut of %-2.2f is requested\n",sigcut);
2561
2562      m = correct_orbit->mon_table;
2563      xmea= 0.0; // ymea = 0.0; // not used
2564      while(m) {
2565        if (get_option("debug")) {
2566      printf("monitor flag: %d\n",m->enable);
2567        }
2568      if(m->enable == 1) {
2569          if(ip == 0) {
2570          bx_m = da1[3][m->id_ttb];
2571          } else if(ip == 1) {
2572          bx_m = da1[6][m->id_ttb];
2573          }
2574          xn = m->val.before[ip]/sqrt(bx_m);
2575          xmea += xn;
2576        if (get_option("debug")) {
2577          printf("==> %s %-4.3f %-4.3f \n",m->p_node->name,bx_m,m->val.before[ip]);
2578          printf("==> %-4.3f \n",xn);
2579        }
2580        im++;
2581      }
2582      m = m->next;
2583      };
2584          xmea = xmea/im;
2585        if (get_option("debug")) {
2586          printf("Mean values: %-4.3f \n",xmea);
2587        }
2588      m = correct_orbit->mon_table;
2589      im = 0;
2590      xsig= 0.0;
2591      while(m) {
2592      if(m->enable == 1) {
2593          if(ip == 0) {
2594          bx_m = da1[3][m->id_ttb];
2595          } else if(ip == 1) {
2596          bx_m = da1[6][m->id_ttb];
2597          }
2598          xn = m->val.before[ip]/sqrt(bx_m);
2599          xsig += (xmea - xn)*(xmea - xn);
2600        im++;
2601      }
2602      m = m->next;
2603      };
2604          xsig = sqrt(xsig/im);
2605        if (get_option("debug")) {
2606          printf("Sigma values: %-4.3f \n",xsig);
2607        }
2608
2609      m = correct_orbit->mon_table;
2610      while(m) {
2611      if(m->enable == 1) {
2612          if(ip == 0) {
2613          bx_m = da1[3][m->id_ttb];
2614          } else if(ip == 1) {
2615          bx_m = da1[6][m->id_ttb];
2616          }
2617          xn = (m->val.before[ip]/sqrt(bx_m)) - xmea;
2618          if(fabs(xn) > (sigcut*xsig)) {
2619             printf("disabled %s %c = %-4.3f (%-4.3f), limit is %-2.2f*%-4.3f\n",
2620                     m->p_node->name,pl[ip],xn,m->val.before[ip],sigcut,xsig);
2621             m->enable = 0;
2622             icnt++;
2623          }
2624      }
2625      m = m->next;
2626      };
2627
2628  return(icnt);
2629}
2630
2631static double*
2632pro_correct_response_ring(int ip, int nc, int nm)
2633{
2634  int    ic, im;
2635  struct id_mic *m, *c;  /* access to tables for monitors and correctors */
2636
2637  struct table *ttb;
2638  double **da1;
2639  double bx_c,by_c,pix_c,piy_c;
2640  double bx_m,by_m,pix_m,piy_m;
2641  double qx0, qy0;
2642  double respx1, respy1;
2643  double respx, respy;
2644  double *dmat;
2645
2646  ttb = model_table;
2647  da1 = ttb->d_cols;
2648  ic = 0; im = 0;
2649
2650  dmat = mycalloc("pro_correct_response_ring",nc*nm,sizeof(double));
2651
2652  correct_orbit->qx0 = da1[5][ttb->curr-1];
2653  correct_orbit->qy0 = da1[8][ttb->curr-1];
2654  qx0 = correct_orbit->qx0;
2655  qy0 = correct_orbit->qy0;
2656
2657  c = correct_orbit->cor_table;
2658  ic = 0;
2659  while(c) {
2660    if (get_option("debug")) {
2661    printf("corrector flag: %d\n",c->enable);
2662    }
2663    if(c->enable == 1) {
2664      bx_c = da1[3][c->id_ttb];
2665      by_c = da1[6][c->id_ttb];
2666      pix_c = da1[5][c->id_ttb];
2667      piy_c = da1[8][c->id_ttb];
2668      m = correct_orbit->mon_table;
2669      im = 0;
2670      while(m) {
2671        if (get_option("debug")) {
2672      printf("monitor flag: %d\n",m->enable);
2673        }
2674      if(m->enable == 1) {
2675        bx_m = da1[3][m->id_ttb];
2676        by_m = da1[6][m->id_ttb];
2677        pix_m = da1[5][m->id_ttb];
2678        piy_m = da1[8][m->id_ttb];
2679        respx = 0.0;
2680        respy = 0.0;
2681        if(ip == 1) {
2682            respx1 = cos((fabs(pix_m - pix_c)*twopi) - qx0*pi);
2683            respx = respx1*sqrt(bx_m*bx_c)/(2.0*sin(pi*qx0));
2684            setup_(&respx, dmat, &im, &ic, &nm, &nc);
2685        }  else if (ip == 2) {
2686            respy1 = cos((fabs(piy_m - piy_c)*twopi) - qy0*pi);
2687            respy = respy1*sqrt(by_m*by_c)/(2.0*sin(pi*qy0));
2688            setup_(&respy, dmat, &im, &ic, &nm, &nc);
2689        }
2690        im++;
2691      }
2692      m = m->next;
2693      };
2694      ic++;
2695    }
2696    c = c->next;
2697  };
2698  return(dmat);
2699}
2700
2701static double*
2702pro_correct_response_line(int ip, int nc, int nm)
2703{
2704  int    ic, im;
2705  struct id_mic *m, *c;  /* access to tables for monitors and correctors */
2706
2707  struct table *ttb;
2708  double **da1;
2709  double bx_c,by_c,pix_c,piy_c;
2710  double bx_m,by_m,pix_m,piy_m;
2711//  double qx0, qy0; // not used
2712  double respx1, respy1;
2713  double respx, respy;
2714  double *dmat;
2715
2716  ttb = model_table;
2717  da1 = ttb->d_cols;
2718  ic = 0; im = 0;
2719
2720  dmat =  mycalloc("pro_correct_response_ring",nc*nm,sizeof(double));
2721
2722  correct_orbit->qx0 = da1[5][ttb->curr-1];
2723  correct_orbit->qy0 = da1[8][ttb->curr-1];
2724//  qx0 = correct_orbit->qx0; // not used
2725//  qy0 = correct_orbit->qy0; // not used
2726
2727  c = correct_orbit->cor_table;
2728  ic = 0;
2729  while(c) {
2730    if(c->enable ==1) {
2731      bx_c = da1[3][c->id_ttb];
2732      by_c = da1[6][c->id_ttb];
2733      pix_c = da1[5][c->id_ttb];
2734      piy_c = da1[8][c->id_ttb];
2735      m = correct_orbit->mon_table;
2736      im = 0;
2737      while(m) {
2738      if(m->enable ==1) {
2739        bx_m = da1[3][m->id_ttb];
2740        by_m = da1[6][m->id_ttb];
2741        pix_m = da1[5][m->id_ttb];
2742        piy_m = da1[8][m->id_ttb];
2743        respx = 0.0;
2744        respy = 0.0;
2745        if(ip == 1) {
2746            if(pix_m > pix_c) {
2747              respx1 = sin((pix_m - pix_c)*twopi);
2748              respx = respx1*sqrt(bx_m*bx_c);
2749            } else {
2750              respx = 0.0;
2751            }
2752            setup_(&respx, dmat, &im, &ic, &nm, &nc);
2753        }  else if (ip == 2) {
2754            if(piy_m > piy_c) {
2755              respy1 = sin((piy_m - piy_c)*twopi);
2756              respy = respy1*sqrt(by_m*by_c);
2757            } else {
2758              respy = 0.0;
2759            }
2760            setup_(&respy, dmat, &im, &ic, &nm, &nc);
2761        }
2762        im++;
2763      }
2764      m = m->next;
2765      };
2766      ic++;
2767    }
2768    c = c->next;
2769  };
2770  return(dmat);
2771}
2772
2773static void
2774pro_correct_make_corr_table(void)
2775{
2776  struct table *ttb;
2777  int j;
2778
2779/*
2780  static char atm[5][4] = {"hmon","vmon","hkic","vkic","kick"};
2781*/
2782
2783/*
2784  ttb = orbin_table;
2785*/
2786  ttb = model_table;
2787
2788  for (j=0; j < ttb->curr; j++) {
2789    if((strncmp(atc[0],ttb->p_nodes[j]->base_name,4) == 0) ||
2790       (strncmp(atc[1],ttb->p_nodes[j]->base_name,4) == 0) ||
2791       (strncmp(atc[2],ttb->p_nodes[j]->base_name,4) == 0))  {
2792      string_to_table_curr("corr","name",ttb->p_nodes[j]->name);
2793      augment_count("corr");
2794    }
2795  }
2796}
2797
2798static void
2799pro_correct2_make_corr_table(void)
2800{
2801  struct id_mic2 *ttb;
2802
2803/*
2804  static char atm[5][4] = {"hmon","vmon","hkic","vkic","kick"};
2805*/
2806
2807  ttb = correct_orbit12->cor_table;
2808
2809  while(ttb != NULL) {
2810    if((strncmp(atc[0],ttb->p_node->base_name,4) == 0) ||
2811       (strncmp(atc[1],ttb->p_node->base_name,4) == 0) ||
2812       (strncmp(atc[2],ttb->p_node->base_name,4) == 0))  {
2813       if(ttb->id_ttb[0] > 0) {
2814         string_to_table_curr("corr1","name",ttb->p_node->name);
2815         augment_count("corr1");
2816       }
2817
2818       if(ttb->id_ttb[1] > 0) {
2819         string_to_table_curr("corr2","name",ttb->p_node->name);
2820         augment_count("corr2");
2821       }
2822    }
2823    ttb = ttb->next;
2824  }
2825}
2826
2827static void
2828pro_correct_make_mon_table(void)
2829{
2830  struct table *ttb;
2831  int j;
2832
2833/*
2834  static char atm[3][4] = {"hmon","vmon","moni"};
2835*/
2836
2837  ttb = model_table;
2838
2839  for (j=0; j < ttb->curr; j++) {
2840    if((strncmp(atm[0],ttb->p_nodes[j]->base_name,4) == 0) ||
2841       (strncmp(atm[1],ttb->p_nodes[j]->base_name,4) == 0) ||
2842       (strncmp(atm[2],ttb->p_nodes[j]->base_name,4) == 0))  {
2843      string_to_table_curr("mon","name",ttb->p_nodes[j]->name);
2844      augment_count("mon");
2845    }
2846  }
2847}
2848
2849static void
2850pro_correct2_make_mon_table(void)
2851{
2852  struct id_mic2 *ttb;
2853/*
2854  static char atm[3][4] = {"hmon","vmon","moni"};
2855*/
2856
2857  ttb = correct_orbit12->mon_table;
2858
2859  while(ttb != NULL) {
2860    if((strncmp(atm[0],ttb->p_node->base_name,4) == 0) ||
2861       (strncmp(atm[1],ttb->p_node->base_name,4) == 0) ||
2862       (strncmp(atm[2],ttb->p_node->base_name,4) == 0))  {
2863      string_to_table_curr("mon","name",ttb->p_node->name);
2864      augment_count("mon");
2865    }
2866    ttb = ttb->next;
2867  }
2868}
2869
2870static void
2871pro_correct_fill_corr_table(int ip ,char *name, double old, double new)
2872{
2873  struct table *cor;
2874
2875  int j;
2876
2877  cor = corr_table;
2878
2879  for (j=0; j < cor->curr; j++) {
2880    if(strcmp(name,cor->s_cols[0][j]) == 0) {
2881      cor->d_cols[ip][j] = old;
2882      cor->d_cols[ip+2][j] = new;
2883    }
2884  }
2885}
2886
2887static void
2888pro_correct2_fill_corr_table(int b, int ip ,char *name, double old, double new)
2889{
2890  struct table *cor = NULL;
2891
2892  int j;
2893 
2894  if((b != 1) && (b != 0)) {
2895      char buf[64];
2896      sprintf(buf, "%d", b);
2897      fatal_error("Invalid beam requested:", buf);
2898  }
2899
2900  if(b == 0) cor =  corr_table1;
2901  if(b == 1) cor =  corr_table2;
2902
2903  for (j=0; j < cor->curr; j++) {
2904    if(strcmp(name,cor->s_cols[0][j]) == 0) {
2905      cor->d_cols[ip][j] = old;
2906      cor->d_cols[ip+2][j] = new;
2907    }
2908  }
2909}
2910
2911static void
2912pro_correct_fill_mon_table(int ip ,char *name, double old, double new)
2913{
2914  struct table *mon;
2915
2916  int j;
2917
2918  mon =  mon_table;
2919
2920  for (j=0; j < mon->curr; j++) {
2921    if(strcmp(name,mon->s_cols[0][j]) == 0) {
2922      mon->d_cols[ip][j] = old*0.001;
2923      mon->d_cols[ip+2][j] = new*0.001;
2924    }
2925  }
2926}
2927
2928static void
2929pro_correct2_fill_mon_table(int ip ,char *name, double old, double new)
2930{
2931  struct table *mon;
2932
2933  int j;
2934
2935  mon =  mon_table;
2936
2937  for (j=0; j < mon->curr; j++) {
2938    if(strcmp(name,mon->s_cols[0][j]) == 0) {
2939      mon->d_cols[ip][j] = old*0.001;
2940      mon->d_cols[ip+2][j] = new*0.001;
2941    }
2942  }
2943}
2944
2945static void
2946pro_correct_write_results(double *monvec, double *resvec, double *corvec, int *nx, int *nc, int *nm, int imon, int icor, int ip)
2947/*                                              */
2948/* Writes a summary of the correction           */
2949/* Writes correctors strengths into sequences   */
2950/* Fills TFS tables for correctors and monitors */
2951/* Fills the 'stren.out' output                 */
2952/* Makes various prints on request              */
2953/*                                              */
2954{
2955  int i;
2956  int rst;
2957  double corrm;
2958  struct id_mic *m, *c;      /* access to tables for monitors and correctors */
2959
2960  m = correct_orbit->mon_table;
2961  c = correct_orbit->cor_table;
2962
2963  if(fddata != NULL) {
2964     rst = get_variable("n");
2965     fprintf(fddata,"%d %d %e %e %e %e %e %e\n",ip,rst,cprp(monvec,imon),cprp(resvec,imon),crms(monvec,imon),crms(resvec,imon),copk(monvec,imon),copk(resvec,imon));
2966  }
2967
2968  if(print_correct_opt > 0) {
2969     printf("CORRECTION SUMMARY:   \n\n");
2970     printf("rms before correction: %f mm\nrms after correction:  %f mm\n\n",crms(monvec,imon),crms(resvec,imon));
2971     printf("ptp before correction: %f mm\nptp after correction:  %f mm\n\n",cprp(monvec,imon),cprp(resvec,imon));
2972  }
2973
2974  if(print_correct_opt > 1) {
2975  printf("Monitor:  Before:     After:    Difference:\n");
2976  printf("           (mm)        (mm)         (mm)   \n");
2977  }
2978
2979  for(i=0;i<imon;i++) {
2980    if(print_correct_opt > 1) {
2981      printf("%s   %-4.3f     %-4.3f     %-4.3f\n",m[nm[i]].p_node->name,monvec[i],resvec[i],resvec[i]-monvec[i]);
2982    }
2983    m[nm[i]].val.after[ip-1] = resvec[i];
2984    pro_correct_fill_mon_table(ip,m[nm[i]].p_node->name,monvec[i],resvec[i]);
2985  }
2986
2987  corrm = copk(corvec,icor);
2988
2989  printf("Max strength: %e should be less than %e\n",corrm,corrl);
2990  if(corrm > corrl) {
2991     printf("++++++ warning: maximum corrector strength larger than limit\n");
2992  }
2993  set_variable("corrmax",&corrm);
2994  if(print_correct_opt > 1) {
2995     printf("Max strength: %e\n",copk(corvec,icor));
2996     printf("Corrector:  Before:     After:    Difference:\n");
2997     printf("             (mrad)     (mrad)       (mrad)  \n");
2998  }
2999
3000  for(i=0;i<icor;i++) {  /* loop over all correctors */
3001
3002    if(print_correct_opt > 1) {
3003      printf("%s %-3.6f %-3.6f %-3.6f\n",c[nc[i]].p_node->name,
3004                                         c[nc[i]].val.before[ip-1],
3005                                         corvec[nx[i]-1]+c[nc[i]].val.before[ip-1], 
3006                                         corvec[nx[i]-1]);
3007    }
3008
3009    c[nc[i]].val.after[ip-1] = corvec[nx[i]-1];
3010    if(ip == 1) {
3011      c[nc[i]].p_node->chkick += c[nc[i]].p_node->other_bv*0.001*corvec[nx[i]-1];
3012      pro_correct_fill_corr_table(ip,
3013                                  c[nc[i]].p_node->name,
3014                                  c[nc[i]].val.before[ip-1]*0.001,
3015                                  c[nc[i]].p_node->chkick);
3016/*                                c[nc[i]].p_node->other_bv*0.001*corvec[nx[i]-1]); */
3017      if(fcdata != NULL) {
3018         fprintf(fcdata,"%s = %e;\n",c[nc[i]].p_node->name,c[nc[i]].p_node->other_bv*0.001*corvec[nx[i]-1]);
3019      }
3020    } else if (ip == 2) {
3021      c[nc[i]].p_node->cvkick += c[nc[i]].p_node->other_bv*0.001*corvec[nx[i]-1];
3022      pro_correct_fill_corr_table(ip,
3023                                  c[nc[i]].p_node->name,
3024                                  c[nc[i]].val.before[ip-1]*0.001,
3025                                  c[nc[i]].p_node->cvkick);
3026/*                                c[nc[i]].p_node->other_bv*0.001*corvec[nx[i]-1]); */
3027      if(fcdata != NULL) {
3028         fprintf(fcdata,"%s = %e;\n",c[nc[i]].p_node->name,c[nc[i]].p_node->other_bv*0.001*corvec[nx[i]-1]);
3029      }
3030    }
3031  } /* end of loop ove correctors */
3032}
3033
3034static int
3035pro_correct_getactive(int ip, int *nm, int *nx, int *nc, double *corvec, double *monvec,char *conm)
3036{
3037  int  imon, icor;
3038  int  imona, icora;
3039  struct id_mic *m, *c;
3040
3041  m = correct_orbit->mon_table;
3042  imon = 0;
3043  imona = 0;
3044  while(m) {
3045    if (get_option("debug")) {
3046      printf("from list: %d %s %s\n",m->id_ttb,m->p_node->name,m->p_node->base_name);
3047      printf("orbit readings: %d %f %f\n",ip,m->val.before[0],m->val.before[1]);
3048    }
3049    if(m->enable == 1) {
3050      monvec[imon] = m->val.before[ip-1];
3051      nm[imon] = imona;
3052      imon++;
3053    }
3054    imona++;
3055    m = m->next;
3056  };
3057
3058  c = correct_orbit->cor_table;
3059  icor = 0;
3060  icora = 0;
3061  while(c) {
3062    if (get_option("debug")) {
3063      printf("from list: %d %d %s %s\n",c->enable,c->id_ttb,c->p_node->name,c->p_node->base_name);
3064      printf("kicker readings: %f %f\n",c->val.before[0],c->val.before[1]);
3065    }
3066    if(c->enable == 1) {
3067      corvec[icor] = c->val.before[ip-1];
3068      nx[icor] = icora;
3069      nc[icor] = icora;
3070      strcpy(conm,c->p_node->name);
3071      conm+=16;
3072      /*          printf("nc: %d %d \n",icor,nc[icor]); */
3073      icor++;
3074    }
3075    icora++;
3076    c = c->next;
3077  };
3078  return(10000*imon + icor);
3079}
3080
3081static void
3082correct_option(struct in_cmd* cmd)
3083{
3084  struct name_list* nl = cmd->clone->par_names;
3085  int i, debug;
3086  int val, pos, seed;
3087
3088  if ((debug=get_option("debug")))  {
3089     fprintf(prt_file, "in coption routine\n");
3090     for(i=0;i<cmd->tok_list->curr;i++) {
3091        fprintf(prt_file, "command(s): %s\n",cmd->tok_list->p[i]);
3092     }
3093  }
3094  if ((pos = name_list_pos("seed", nl)) > -1)
3095    {
3096     if (nl->inform[pos])
3097       {
3098        seed = command_par_value("seed", cmd->clone);
3099        init55(seed);
3100       }
3101    }
3102 val = command_par_value("print", cmd->clone);
3103 if(val == 0) {
3104      if (debug)  fprintf(prt_file, "print option not set\n");
3105      print_correct_opt = 0;
3106 }else {
3107      if (debug)  fprintf(prt_file, "print option set\n");
3108      print_correct_opt = val;
3109 }
3110 val = command_par_value("debug", cmd->clone);
3111 if(val == 0) {
3112      if (debug)  fprintf(prt_file, "debug option not set\n");
3113      debug_correct_opt = 0;
3114 }else {
3115      if (debug)  fprintf(prt_file, "debug option set\n");
3116      debug_correct_opt = val;
3117 }
3118
3119}
3120
3121static void
3122correct_getorbit(struct in_cmd* cmd)
3123{
3124  (void)cmd;
3125}
3126
3127static void
3128correct_putorbit(struct in_cmd* cmd)
3129{
3130  int i;
3131  struct name_list* nl;
3132  char* filename = command_par_string("file", cmd->clone);
3133  char* table_name;
3134  current_twiss = clone_command(find_command("twiss", defined_commands));
3135  nl = current_twiss->par_names;
3136  for (i = 0; i < nl->curr; i++) nl->inform[i] = 0;
3137  pro_twiss();
3138  table_name = permbuff("orbit");
3139  orbit_table = make_table(table_name, "orbit", orbit_table_cols,
3140            orbit_table_types, current_sequ->n_nodes);
3141  add_to_table_list(orbit_table, table_register);
3142  fill_orbit_table(orbit_table, orbin_table);
3143  out_table("orbit", orbit_table, filename);
3144  current_twiss = delete_command(current_twiss);
3145}
3146
3147static void
3148correct_usekick(struct in_cmd* cmd)
3149{
3150  char temp[12];
3151  int count = set_enable("kicker", cmd);
3152  sprintf(temp, "%d", count);
3153  put_info(temp, "corrector(s) affected");
3154}
3155
3156static void
3157correct_usemonitor(struct in_cmd* cmd)
3158{
3159  char temp[12];
3160  int count = set_enable("monitor", cmd);
3161  sprintf(temp, "%d", count);
3162  put_info(temp, "monitor(s) affected");
3163}
3164
3165// public interface
3166
3167void
3168store_orbit(struct command* comm, double* orbit)
3169{
3170  struct name_list* nl = comm->par_names;
3171  if (nl->inform[name_list_pos("x", nl)])
3172    orbit[0] = command_par_value("x",comm);
3173  if (nl->inform[name_list_pos("px", nl)])
3174    orbit[1] = command_par_value("px",comm);
3175  if (nl->inform[name_list_pos("y", nl)])
3176    orbit[2] = command_par_value("y",comm);
3177  if (nl->inform[name_list_pos("py", nl)])
3178    orbit[3] = command_par_value("py",comm);
3179  if (nl->inform[name_list_pos("t", nl)])
3180    orbit[4] = command_par_value("t",comm);
3181  if (nl->inform[name_list_pos("pt", nl)])
3182    orbit[5] = command_par_value("pt",comm);
3183}
3184
3185void
3186pro_correct(struct in_cmd* cmd)
3187{
3188  if (strcmp(cmd->tok_list->p[0], "correct") == 0)
3189    {
3190     correct_correct(cmd);
3191    }
3192  else if (strcmp(cmd->tok_list->p[0], "usekick") == 0)
3193    {
3194     correct_usekick(cmd);
3195    }
3196  else if (strcmp(cmd->tok_list->p[0], "usemonitor") == 0)
3197    {
3198     correct_usemonitor(cmd);
3199    }
3200  else if (strcmp(cmd->tok_list->p[0], "getorbit") == 0)
3201    {
3202     correct_getorbit(cmd);
3203    }
3204  else if (strcmp(cmd->tok_list->p[0], "putorbit") == 0)
3205    {
3206     correct_putorbit(cmd);
3207    }
3208  else if (strcmp(cmd->tok_list->p[0], "readmytable") == 0)
3209    {
3210     read_my_table(cmd);
3211    }
3212  else if (strcmp(cmd->tok_list->p[0], "readcorr") == 0)
3213    {
3214     correct_readcorr(cmd);
3215    }
3216  else if (strcmp(cmd->tok_list->p[0], "setcorr") == 0)
3217    {
3218     correct_setcorr(cmd);
3219    }
3220  else if (strcmp(cmd->tok_list->p[0], "coption") == 0)
3221    {
3222     correct_option(cmd);
3223    }
3224}
3225
3226int
3227locf_(char *iadr)   
3228#define NADUPW 4   /* Number of ADdress Units Per Word */
3229#define LADUPW 2   /* Logarithm base 2 of ADdress Units Per Word */
3230{
3231        return (uintptr_t) iadr >> LADUPW;
3232}
3233
3234void
3235f_ctof(int *j, char *string, int *nel)
3236{
3237  long i, flg = 0;
3238
3239  for(i = 0; i < *nel; i++) {
3240    if(flg == 1) {
3241      string[i] = ' ';
3242      continue;
3243    }
3244
3245    if(string[i] == '\0') {
3246      string[i] = ' ';
3247      flg = 1;
3248      continue;
3249    }
3250  }
3251  *j = i;
3252}
3253
Note: See TracBrowser for help on using the repository browser.