source: TRACY3/trunk/tracy/tracy/src/rdmfile.cc @ 11

Last change on this file since 11 was 11, checked in by zhangj, 11 years ago
  • Property svn:executable set to *
File size: 12.4 KB
Line 
1/* Tracy-2
2
3   J. Bengtsson, CBP, LBL      1990 - 1994   Pascal version
4                 SLS, PSI      1995 - 1997
5   M. Boege      SLS, PSI      1998          C translation
6   L. Nadolski   SOLEIL        2002          Link to NAFF, Radia field maps
7   J. Bengtsson  NSLS-II, BNL  2004 -       
8
9
10   To generate a lattice flat file.
11
12   Type codes:
13
14     marker     -1
15     drift       0
16     multipole   1
17     cavity      2
18     thin kick   3
19     wiggler     4
20
21   Integration methods:
22
23     linear, matrix style (obsolete)              0
24     2nd order symplectic integrator (obsolete)   2
25     4th order symplectic integrator              4
26
27   Format:
28
29     name, family no, kid no, element no
30     type code, integration method, no of integration steps
31     apertures: xmin, xmax, ymin, ymax
32
33   The following lines follows depending on element type.
34
35     type
36
37     drift:      L
38
39     multipole:  hor., ver. displ., roll angle (design), roll angle (error)
40                 L, 1/rho, entrance angle, exit angle
41                 no of nonzero multipole coeff., n design
42                 n, b , a
43                     n   n
44                     .
45                     .
46                     .
47
48     wiggler:    L [m], lambda [m]
49                 no of harmonics
50                 harm no, kxV [1/m], BoBrhoV [1/m], kxH, BoBrhoH, phi
51                    ...
52
53     cavity:     cavity voltage/beam energy [eV], omega/c, beam energy [eV]
54
55     thin kick:  hor., ver. displacement, roll angle (total)
56                 no of nonzero multipole coeff.
57                 n, b , a
58                     n   n
59                     .
60                     .
61                     .
62
63     kick_map:   scale order <file name>
64
65*/
66
67
68#define snamelen   10
69
70// numerical type codes
71#define marker_   -1
72#define drift_     0
73#define mpole_     1
74#define cavity_    2
75#define thinkick_  3
76#define wiggler_   4
77#define insertion_ 6
78
79
80ifstream  inf;
81
82
83void get_kind(const int kind, elemtype &Elem)
84{
85
86  switch (kind) {
87  case marker_:
88    Elem.Pkind = PartsKind(marker);
89    break;
90  case drift_:
91    Elem.Pkind = PartsKind(drift);
92    Drift_Alloc(&Elem);
93    break;
94  case mpole_:
95    Elem.Pkind = PartsKind(Mpole);
96    Mpole_Alloc(&Elem);
97    Elem.M->Pthick = pthicktype(thick);
98    break;
99  case cavity_:
100    Elem.Pkind = PartsKind(Cavity);
101    Cav_Alloc(&Elem);
102    break;
103  case thinkick_:
104    Elem.Pkind = PartsKind(Mpole);
105    Mpole_Alloc(&Elem);
106    Elem.M->Pthick = pthicktype(thin);
107    break;
108  case wiggler_:
109    Elem.Pkind = PartsKind(Wigl);
110    Wiggler_Alloc(&Elem);
111    break;
112  case insertion_:
113    Elem.Pkind = PartsKind(Insertion);
114    Insertion_Alloc(&Elem);
115    break;
116  default:
117    cout << "get_kind: unknown type " << kind << " " << Elem.PName << endl;
118    exit_(1);
119    break;
120  }
121}
122
123
124void rdmfile(const char *mfile_dat)
125 {
126    char line[max_str], file_name[max_str];
127    int j, nmpole, kind, method, n;
128    long int i;
129    double dTerror;
130
131    bool prt = false;
132
133    cout << endl;
134    cout << "reading machine file: " << mfile_dat << endl;
135
136    file_rd(inf, mfile_dat);
137
138    while (inf.getline(line, max_str) != NULL) {
139        if (prt)
140            printf("%s\n", line);
141        sscanf(line, "%*s %*d %*d %ld", &i);
142
143        Cell[i].dS[X_] = 0.0;
144        Cell[i].dS[Y_] = 0.0;
145        Cell[i].dT[X_] = 1.0;
146        Cell[i].dT[Y_] = 0.0;
147
148        sscanf(line, "%s %ld %ld", Cell[i].Elem.PName, &Cell[i].Fnum,
149                &Cell[i].Knum);
150
151        if (Cell[i].Knum == 1) {
152            strcpy(ElemFam[Cell[i].Fnum - 1].ElemF.PName, Cell[i].Elem.PName);
153            globval.Elem_nFam = max(Cell[i].Fnum, globval.Elem_nFam);
154        }
155
156        if (i > 0) {
157            ElemFam[Cell[i].Fnum - 1].KidList[Cell[i].Knum - 1] = i;
158            ElemFam[Cell[i].Fnum - 1].nKid = max(Cell[i].Knum,
159                    ElemFam[Cell[i].Fnum - 1].nKid);
160        }
161
162        inf.getline(line, max_str);
163        if (prt)
164            printf("%s\n", line);
165        sscanf(line, "%d %d %d", &kind, &method, &n);
166        get_kind(kind, Cell[i].Elem);
167        if (i > 0)
168            ElemFam[Cell[i].Fnum - 1].ElemF.Pkind = Cell[i].Elem.Pkind;
169
170        inf.getline(line, max_str);
171        if (prt)
172            printf("%s\n", line);
173        sscanf(line, "%lf %lf %lf %lf", &Cell[i].maxampl[X_][0],
174                &Cell[i].maxampl[X_][1], &Cell[i].maxampl[Y_][0],
175                &Cell[i].maxampl[Y_][1]);
176
177        Cell[i].Elem.PL = 0.0;
178
179        switch (Cell[i].Elem.Pkind) {
180        case undef:
181            cout << "rdmfile: unknown type " << i << endl;
182            exit_(1);
183            break;
184        case marker:
185            break;
186        case drift:
187            inf.getline(line, max_str);
188            if (prt)
189                printf("%s\n", line);
190            sscanf(line, "%lf", &Cell[i].Elem.PL);
191            break;
192        case Cavity:
193            inf.getline(line, max_str);
194            if (prt)
195                printf("%s\n", line);
196            sscanf(line, "%lf %lf %d %lf", &Cell[i].Elem.C->Pvolt,
197                    &Cell[i].Elem.C->Pfreq, &Cell[i].Elem.C->Ph,
198                    &globval.Energy);
199            globval.Energy *= 1e-9;
200            Cell[i].Elem.C->Pvolt *= globval.Energy * 1e9;
201            Cell[i].Elem.C->Pfreq *= c0 / (2.0 * M_PI);
202            break;
203        case Mpole:
204            Cell[i].Elem.M->Pmethod = method;
205            Cell[i].Elem.M->PN = n;
206
207            if (Cell[i].Elem.M->Pthick == thick) {
208                inf.getline(line, max_str);
209                if (prt)
210                    printf("%s\n", line);
211                sscanf(line, "%lf %lf %lf %lf", &Cell[i].dS[X_],
212                        &Cell[i].dS[Y_], &Cell[i].Elem.M->PdTpar, &dTerror);
213                Cell[i].dT[X_] = cos(dtor(dTerror + Cell[i].Elem.M->PdTpar));
214                Cell[i].dT[Y_] = sin(dtor(dTerror + Cell[i].Elem.M->PdTpar));
215
216                inf.getline(line, max_str);
217                if (prt)
218                    printf("%s\n", line);
219                sscanf(line, "%lf %lf %lf %lf %lf  %lf %lf", &Cell[i].Elem.PL,
220                        &Cell[i].Elem.M->Pirho, &Cell[i].Elem.M->PTx1,
221                        &Cell[i].Elem.M->PTx2, &Cell[i].Elem.M->PH1, &Cell[i].Elem.M->PH2, &Cell[i].Elem.M->Pgap);
222                if (Cell[i].Elem.M->Pirho != 0.0)
223                    Cell[i].Elem.M->Porder = 1;
224            } else {
225                inf.getline(line, max_str);
226                if (prt)
227                    printf("%s\n", line);
228                sscanf(line, "%lf %lf %lf", &Cell[i].dS[X_], &Cell[i].dS[Y_],
229                        &dTerror);
230                Cell[i].dT[X_] = cos(dtor(dTerror));
231                Cell[i].dT[Y_] = sin(dtor(dTerror));
232            }
233
234            Cell[i].Elem.M->Pc0 = sin(Cell[i].Elem.PL * Cell[i].Elem.M->Pirho
235                    / 2.0);
236            Cell[i].Elem.M->Pc1 = cos(dtor(Cell[i].Elem.M->PdTpar))
237                    * Cell[i].Elem.M->Pc0;
238            Cell[i].Elem.M->Ps1 = sin(dtor(Cell[i].Elem.M->PdTpar))
239                    * Cell[i].Elem.M->Pc0;
240
241            inf.getline(line, max_str);
242            if (prt)
243                printf("%s\n", line);
244            sscanf(line, "%d %d", &nmpole, &Cell[i].Elem.M->n_design);
245            for (j = 1; j <= nmpole; j++) {
246                inf.getline(line, max_str);
247                if (prt)
248                    printf("%s\n", line);
249                sscanf(line, "%d", &n);
250                sscanf(line, "%*d %lf %lf", &Cell[i].Elem.M->PB[HOMmax + n],
251                        &Cell[i].Elem.M->PB[HOMmax - n]);
252                Cell[i].Elem.M->PBpar[HOMmax + n] = Cell[i].Elem.M->PB[HOMmax
253                        + n];
254                Cell[i].Elem.M->PBpar[HOMmax - n] = Cell[i].Elem.M->PB[HOMmax
255                        - n];
256                Cell[i].Elem.M->Porder = max(n, Cell[i].Elem.M->Porder);
257            }
258            break;
259        case Wigl:
260            Cell[i].Elem.W->Pmethod = method;
261            Cell[i].Elem.W->PN = n;
262
263            inf.getline(line, max_str);
264            if (prt)
265                printf("%s\n", line);
266            sscanf(line, "%lf %lf", &Cell[i].Elem.PL, &Cell[i].Elem.W->lambda);
267
268            inf.getline(line, max_str);
269            if (prt)
270                printf("%s\n", line);
271            sscanf(line, "%d", &Cell[i].Elem.W->n_harm);
272
273            if (Cell[i].Knum == 1)
274                Wiggler_Alloc(&ElemFam[Cell[i].Fnum - 1].ElemF);
275            for (j = 0; j < Cell[i].Elem.W->n_harm; j++) {
276                inf.getline(line, max_str);
277                if (prt)
278                    printf("%s\n", line);
279                sscanf(line, "%d %lf %lf %lf %lf %lf",
280                        &Cell[i].Elem.W->harm[j], &Cell[i].Elem.W->kxV[j],
281                        &Cell[i].Elem.W->BoBrhoV[j], &Cell[i].Elem.W->kxH[j],
282                        &Cell[i].Elem.W->BoBrhoH[j], &Cell[i].Elem.W->phi[j]);
283                ElemFam[Cell[i].Fnum - 1].ElemF.W->BoBrhoV[j]
284                        = Cell[i].Elem.W->BoBrhoV[j];
285                ElemFam[Cell[i].Fnum - 1].ElemF.W->BoBrhoH[j]
286                        = Cell[i].Elem.W->BoBrhoH[j];
287            }
288            break;
289        case Insertion:
290            Cell[i].Elem.ID->Pmethod = method;
291            Cell[i].Elem.ID->PN = n;
292
293            inf.getline(line, max_str);
294            if (prt)
295                printf("%s\n", line);
296            sscanf(line, "%lf %d %s", &Cell[i].Elem.ID->scaling1, &n, file_name);
297
298            if (n == 1) {
299                Cell[i].Elem.ID->firstorder = true;
300                Cell[i].Elem.ID->secondorder = false;
301
302                strcpy(Cell[i].Elem.ID->fname1, file_name);
303                Read_IDfile(Cell[i].Elem.ID->fname1, &Cell[i].Elem.PL,
304                        &Cell[i].Elem.ID->nx1, &Cell[i].Elem.ID->nz1,
305                        Cell[i].Elem.ID->tabx1, Cell[i].Elem.ID->tabz1,
306                        Cell[i].Elem.ID->thetax1, Cell[i].Elem.ID->thetaz1);
307            } else if (n == 2) {
308                Cell[i].Elem.ID->firstorder = false;
309                Cell[i].Elem.ID->secondorder = true;
310
311                strcpy(Cell[i].Elem.ID->fname2, file_name);
312                Read_IDfile(Cell[i].Elem.ID->fname2, &Cell[i].Elem.PL,
313                        &Cell[i].Elem.ID->nx2, &Cell[i].Elem.ID->nz2,
314                        Cell[i].Elem.ID->tabx2, Cell[i].Elem.ID->tabz2,
315                        Cell[i].Elem.ID->thetax2, Cell[i].Elem.ID->thetaz2);
316            } else {
317                cout << "rdmfile: undef order " << n << endl;
318                exit_(1);
319            }
320
321            if (Cell[i].Elem.ID->Pmethod == 1)
322                Cell[i].Elem.ID->linear = true;
323            else
324                Cell[i].Elem.ID->linear = false;
325
326            if (!Cell[i].Elem.ID->linear) {
327                Cell[i].Elem.ID->tx1 = dmatrix(1, Cell[i].Elem.ID->nz1, 1,
328                        Cell[i].Elem.ID->nx1);
329                Cell[i].Elem.ID->tz1 = dmatrix(1, Cell[i].Elem.ID->nz1, 1,
330                        Cell[i].Elem.ID->nx1);
331                Cell[i].Elem.ID->tx2 = dmatrix(1, Cell[i].Elem.ID->nz2, 1,
332                        Cell[i].Elem.ID->nx2);
333                Cell[i].Elem.ID->tz2 = dmatrix(1, Cell[i].Elem.ID->nz2, 1,
334                        Cell[i].Elem.ID->nx2);
335                Cell[i].Elem.ID->TabxOrd1 = (double *) malloc(
336                        (Cell[i].Elem.ID->nx1) * sizeof(double));
337                Cell[i].Elem.ID->TabzOrd1 = (double *) malloc(
338                        (Cell[i].Elem.ID->nz1) * sizeof(double));
339                Cell[i].Elem.ID->TabxOrd2 = (double *) malloc(
340                        (Cell[i].Elem.ID->nx2) * sizeof(double));
341                Cell[i].Elem.ID->TabzOrd2 = (double *) malloc(
342                        (Cell[i].Elem.ID->nz2) * sizeof(double));
343                Cell[i].Elem.ID->f2x1 = dmatrix(1, Cell[i].Elem.ID->nz1, 1,
344                        Cell[i].Elem.ID->nx1);
345                Cell[i].Elem.ID->f2z1 = dmatrix(1, Cell[i].Elem.ID->nz1, 1,
346                        Cell[i].Elem.ID->nx1);
347                Cell[i].Elem.ID->f2x2 = dmatrix(1, Cell[i].Elem.ID->nz2, 1,
348                        Cell[i].Elem.ID->nx2);
349                Cell[i].Elem.ID->f2z2 = dmatrix(1, Cell[i].Elem.ID->nz2, 1,
350                        Cell[i].Elem.ID->nx2);
351                Matrices4Spline(Cell[i].Elem.ID, 1);
352                Matrices4Spline(Cell[i].Elem.ID, 2);
353            }
354
355            break;
356        case FieldMap:
357            break;
358        default:
359            cout << "rdmfile: unknown type" << endl;
360            exit_(1);
361            break;
362        }
363
364        if (i == 0)
365            Cell[i].S = 0.0;
366        else
367            Cell[i].S = Cell[i - 1].S + Cell[i].Elem.PL;
368    }
369
370    globval.Cell_nLoc = i;
371
372    globval.dPcommon = 1e-8;
373    globval.CODeps = 1e-14;
374    globval.CODimax = 40;
375
376    SI_init();
377
378    cout << endl;
379    cout << "rdmfile: read " << globval.Cell_nLoc << " elements" << endl;
380
381    inf.close();
382}
Note: See TracBrowser for help on using the repository browser.