Newer
Older
#define ORDER 1
int no_tps = ORDER; // arbitrary TPSA order is defined locally
extern bool freq_map;
#include "tracy_lib.h"
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
#include <sys/time.h>
//***************************************************************************************
//
// MAIN CODE
//
//****************************************************************************************
int main(int argc, char *argv[])
{
const long seed = 1121;
const bool All = TRUE;
iniranf(seed); setrancut(2.0);
// turn on globval.Cavity_on and globval.radiation to get proper synchr radiation damping
// IDs accounted too if: wiggler model and symplectic integrator (method = 1)
globval.H_exact = false;
globval.quad_fringe = false; // quadrupole fringe field
globval.radiation = false; // synchrotron radiation
globval.emittance = false; // emittance
globval.pathlength = false;
// globval.bpm = 0;
// overview, on energy: 25-15
//const double x_max_FMA = 20e-3, y_max_FMA = 10e-3; //const x_max_FMA = 25e-3, y_max_FMA = 15e-3;
//const int n_x = 80, n_y = 80, n_tr = 2048;
// overview, off energy: 10-10
const double x_max_FMA = 10e-3, delta_FMA = 10e-2;
const int n_x = 80, n_dp = 80, n_tr = 2048;
//
// zoom, on energy: 8-2.5
//const double x_max_FMA = 8e-3, y_max_FMA = 2.5e-3;
//const int n_x = 64, n_y = 15, n_tr = 2048;
// zoom, off energy: 7-3
//const double x_max_FMA = 3e-3, delta_FMA = 7e-2;
//const int n_x = 28, n_dp = 56, n_tr = 2048;
double nux=0, nuz=0, ksix=0, ksiz=0;
bool chroma;
double dP = 0.0;
long lastpos = -1L;
char str1[S_SIZE];
/************************************************************************
start read in files and flags
*************************************************************************/
read_script(argv[1], true);
/************************************************************************
end read in files and flags
*************************************************************************/
// if (true)
// // Read_Lattice("/home/nadolski/codes/tracy/maille/soleil/solamor2_tracy3");
// Read_Lattice(argv[1]); //sets some globval params
// else
// rdmfile("flat_file.dat"); //instead of reading lattice file, get data from flat file
//no_sxt(); //turns off sextupoles
// Ring_GetTwiss(true, 0e-2); //gettwiss computes one-turn matrix arg=(w or w/o chromat, dp/p)
//get_matching_params_scl();
//get_alphac2();
//GetEmittance(ElemIndex("cav"), true);
//prt_lat("linlat.out", globval.bpm, true); /* print lattice file for nsrl-ii*/
prtmfile("flat_file.dat"); // writes flat file /* very important file for debug*/
//prt_chrom_lat(); //writes chromatic functions into chromlat.out
// printlatt(); /* print out lattice functions */
/* print lattice file */
// prt_lat("linlatBNL.out", globval.bpm, All); // BNL print for all elements
printlatt(); /* SOLEIL print out lattice functions */
printglob();
// Flag factory
// bool TuneTracFlag = true;
// bool ChromTracFlag = true;
//*************************************************************
//=============================================================
// Chamber factory
if (ChamberFlag == false)
ChamberOff(); // turn off vacuum chamber setting, use the default one
else if (ChamberNoU20Flag == true)
DefineChNoU20(); // using vacuum chamber setting but without undulator U20
else if (ReadChamberFlag == true)
ReadCh(chamber_file); /* read vacuum chamber from a file "chamber_file" , soleil version*/
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
//LoadApers("Apertures.dat", 1.0, 1.0); /* read vacuum chamber definition for bnl */
PrintCh();
// compute tunes by tracking (should be the same as by DA)
if (TuneTracFlag == true) {
GetTuneTrac(1026L, 0.0, &nux, &nuz);
fprintf(stdout,"From tracking: nux = % f nuz = % f \n",nux,nuz);
}
// compute chromaticities by tracking (should be the same as by DA)
if (ChromTracFlag == true){
GetChromTrac(2L, 1026L, 1e-5, &ksix, &ksiz);
fprintf(stdout,"From tracking: ksix= % f ksiz= % f \n",ksix,ksiz);
}
if (FitTuneFlag == true){
fprintf(stdout, "\n Fitting tunes\n");
FitTune(ElemIndex("qp7"),ElemIndex("qp9"), targetnux, targetnuz);
Ring_GetTwiss(chroma=true, 0.0); /* Compute and get Twiss parameters */
printglob(); /* print parameter list */
}
if (FitChromFlag == true){
fprintf(stdout, "\n Fitting chromaticities\n");
FitChrom(ElemIndex("sx9"),ElemIndex("sx10"), targetksix, targetksiz);
Ring_GetTwiss(chroma=true, 0.0); /* Compute and get Twiss parameters */
printglob(); /* print parameter list */
}
//SetKLpar(ElemIndex("QT"), 1, 2L, 0.001026770838382);
// coupling calculation
if (CouplingFlag == true){
Ring_GetTwiss(chroma=true, 0.0); /* Compute and get Twiss parameters */
printlatt(); /* dump linear lattice functions into "linlat.dat" */
Coupling_Edwards_Teng();
printglob(); /* print parameter list */
}
// add coupling by random rotating of the quadrupoles
if (ErrorCouplingFlag == true){
SetErr();
Ring_GetTwiss(chroma=true, 0.0); /* Compute and get Twiss parameters */
printlatt(); /* dump linear lattice functions into "linlat.dat" */
Coupling_Edwards_Teng();
printglob(); /* print parameter list */
}
// WARNING Fit tunes and chromaticities before applying errors !!!!
//set multipoles in all magnets
if (MultipoleFlag == true ){
if (ThinsextFlag ==true){
fprintf(stdout, "\n Setting Multipoles for lattice with thin sextupoles \n");
Multipole_thinsext(); /* for thin sextupoles */
Ring_GetTwiss(chroma=true, 0.0); /* Compute and get Twiss parameters */
printglob();
}
else{
fprintf(stdout, "\n Setting Multipoles for lattice with thick sextupoles \n");
Multipole_thicksext(); /* for thick sextupoles */
Ring_GetTwiss(chroma=true, 0.0); /* Compute and get Twiss parameters */
printglob();
}
}
/* print parameter list */
// PX2 chicane
// if (PX2Flag ==true){
// setPX2chicane();
// //get closed orbit
// getcod (0.0, &lastpos);
// printcod();
// Ring_GetTwiss(chroma=true, 0.0); /* Compute and get Twiss parameters */
// printglob(); /* print parameter list */
// }
// Computes FMA
if (FmapFlag == true){
if (ChamberFlag == true ){
if (ExperimentFMAFlag == true)
fmap(40,12,258,-20e-3,5e-3,0.0,true); // for experimental
if (DetailedFMAFlag == true)
fmap(100,50,1026,20e-3,5e-3,0.0,true);
}
else{
if (ExperimentFMAFlag == true)
fmap(40,12,258,-32e-3,5e-3,0.0,true);
if (DetailedFMAFlag == true)
fmap(200,100,1026,32e-3,7e-3,0.0,true);
}
}
if (CodeComparaisonFlag){
// SOLEIL
fmap(100,50,1026,32e-3,7e-3,0.0,true);
//fmap(200,100,1026,-32e-3,7e-3,0.0,true);
}
if (MomentumAccFlag == true){
//MomentumAcceptance(10L, 28L, 0.01, 0.05, 4L, -0.01, -0.05, 4L);
MomentumAcceptance(1L, 28L, 0.01, 0.05, 40L, -0.01, -0.05, 40L);
// MomentumAcceptance(1L, 108L, 0.01, 0.05, 100L, -0.01, -0.05, 100L);
}
// computes Tuneshift with amplitudes
if (TuneShiftFlag == true){
if (ChamberFlag == true ){
NuDx(31L,21L,516L,0.025,0.005,dP);
//NuDp(31L,516L,0.06);
//NuDp(31L,1026L,0.06);
}
else{
NuDx(50L,30L,516L,0.035,0.02,dP);
NuDp(31L,1026L,0.06);
}
}
// if (SigmaFlag == true){printsigma();
// }
//
// induced amplitude
if (InducedAmplitudeFlag == true){
InducedAmplitude(193L);
}
if (EtaFlag == true){
// compute cod and twiss parameters for different energy offsets
for (int ii=0; ii<=40; ii++) {
dP = -0.02+ 0.001*ii;
Ring_GetTwiss(chroma=false, dP); /* Compute and get Twiss parameters */
printlatt(); /* dump linear lattice functions into "linlat.dat" */
getcod (dP, lastpos);
// printcod();
prt_cod("cod.out", globval.bpm, true);
//system("mv linlat.out linlat_ooo.out");
sprintf(str1, "mv cod.out cod_%02d.out", ii);
system(str1);
sprintf(str1, "mv linlat.out linlat_%02d.out", ii);
system(str1);
}
}
if (PhaseSpaceFlag == true){
int t_sec,t_usec,t_msec;
struct timeval tv,tv1;
struct timezone tz,tz1 ;
struct tm *tm,*tm1;
gettimeofday(&tv, &tz);
tm=localtime(&tv.tv_sec);
t_sec = - tm->tm_sec;
printf(" %d:%02d:%02d %d \n", tm->tm_hour, tm->tm_min,
tm->tm_sec, tv.tv_usec);
Phase(0.001,0,0.001,0,0.001,0, 1);
gettimeofday(&tv1, &tz1);
tm1=localtime(&tv1.tv_sec);
t_sec += tm->tm_sec;
t_usec = tv1.tv_usec - tv.tv_usec;
t_msec = t_sec*1000 + t_usec/1000;
printf(" %d:%02d:%02d %d \n", tm1->tm_hour, tm1->tm_min,
tm1->tm_sec, tv1.tv_usec);
printf("%d %d %d \n",t_sec*1000,t_usec/1000, t_msec);
printf("the simulation time for phase space in tracy 3 is %ld milliseconds\n",t_msec);
//
// IBS & TOUSCHEK
//
int k, n_turns;
double sigma_s, sigma_delta, tau, alpha_z, beta_z, gamma_z, eps[3];
FILE *outf;
double sum_delta[globval.Cell_nLoc+1][2];
double sum2_delta[globval.Cell_nLoc+1][2];
GetEmittance(ElemIndex("cav"), true);
// initialize momentum aperture arrays
for(k = 0; k <= globval.Cell_nLoc; k++){
sum_delta[k][0] = 0.0; sum_delta[k][1] = 0.0;
sum2_delta[k][0] = 0.0; sum2_delta[k][1] = 0.0;
}
//sigma_delta = 7.70e-04; //410:7.70e-4, 411:9.57e-4, 412:9.12e-4
//globval.eps[X_] = 0.326e-9; //410:3.26e-10, 411:2.63e-10, 412:2.01e-10
globval.eps[Y_] = 0.008e-9;
//sigma_s = 9.73e-3; //410:9.73e-3, 411:12.39e-3, 412:12.50e-3/10.33e-3
//globval.eps[Z_] = sigma_delta*sigma_s;
//globval.delta_RF given by cav voltage in lattice file
//globval.delta_RF = 6.20e-2; //410:6.196e-2, 411:5.285e-2, 412:4.046e-2/5.786e-2
// n_turns = 490; //410:490(735), 411:503(755), 412:439(659)/529(794)
n_turns = 40; //410:490(735), 411:503(755), 412:439(659)/529(794)
-globval.Ascr[ct_][ct_]*globval.Ascr[delta_][ct_]
-globval.Ascr[ct_][delta_]*globval.Ascr[delta_][delta_];
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
beta_z = sqr(globval.Ascr[ct_][ct_]) + sqr(globval.Ascr[ct_][delta_]);
gamma_z = (1+sqr(alpha_z))/beta_z;
sigma_delta = sqrt(gamma_z*globval.eps[Z_]);
sigma_s = sqrt(beta_z*globval.eps[Z_]);//50e-3
beta_z = sqr(sigma_s)/globval.eps[Z_];
alpha_z = sqrt(beta_z*gamma_z-1);
// INCLUDE LC (LC changes sigma_s and eps_z, but has no influence on sigma_delta)
if (false) {
double newLength, bunchLengthening;
newLength = 50e-3;
bunchLengthening = newLength/sigma_s;
sigma_s = newLength;
globval.eps[Z_] = globval.eps[Z_]*bunchLengthening;
beta_z = beta_z*bunchLengthening;
gamma_z = gamma_z/bunchLengthening;
alpha_z = sqrt(beta_z*gamma_z-1); // this doesn't change
}
//globval.eps[X_] = 0.362e-9;
//sigma_delta = 1.04e-3;
//sigma_s = 14.8e-3;
Touschek(Qb, globval.delta_RF, globval.eps[X_], globval.eps[Y_],
sigma_delta, sigma_s);
// initialize eps_IBS with eps_SR
for(k = 0; k < 3; k++)
eps[k] = globval.eps[k];
for(k = 0; k < 20; k++) //prototype (looping because IBS routine doesn't check convergence)
IBS(Qb, globval.eps, eps, alpha_z, beta_z);
}
// TOUSCHEK TRACKING
// Calculate Touschek lifetime
// with the momentum acceptance which is determined by
// the RF acceptance delta_RF and the momentum aperture
// at each element location which is tracked over n turns,
//the vacuum chamber is read from the file "chamber_file"
// finally, write the momentum acceptance to file "mom_aper.out".
if (TousTrackFlag == true) {
ReadCh(chamber_file);
// LoadApers("/home/simon/projects/src/lattice/Apertures.dat", 1, 1);
tau = Touschek(Qb, globval.delta_RF, false,
globval.eps[X_], globval.eps[Y_],
sigma_delta, sigma_s,
n_turns, true, sum_delta, sum2_delta); //the TRUE flag requires apertures loaded
printf("Touschek lifetime = %10.3e hrs\n", tau/3600.0);
outf = file_write("mom_aper.out");
for(k = 0; k <= globval.Cell_nLoc; k++)
fprintf(outf, "%4d %7.2f %5.3f %6.3f\n",
k, Cell[k].S, 1e2*sum_delta[k][0], 1e2*sum_delta[k][1]);
fclose(outf);
}
}
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
/*********************************************************************************
Delicated for max4 lattice. To load alignment error files and do correction
----------------------------------------------------------------------------------
*********************************************************************************/
if (false) {
//prt_cod("cod.out", globval.bpm, true); //prints a specific closed orbit with corrector strengths
globval.bpm = ElemIndex("bpm_m"); //definition for max4 lattice, bpm
// globval.bpm = ElemIndex("bpm");
globval.hcorr = ElemIndex("corr_h"); globval.vcorr = ElemIndex("corr_v"); //definition for max4 lattice, correctors
// globval.hcorr = ElemIndex("ch"); globval.vcorr = ElemIndex("cv");
globval.gs = ElemIndex("GS"); globval.ge = ElemIndex("GE"); //definition for max4 lattice, girder maker
//compute response matrix (needed for OCO)
gcmat(globval.bpm, globval.hcorr, 1); gcmat(globval.bpm, globval.vcorr, 2);
//print response matrix (routine in lsoc.cc)
//prt_gcmat(globval.bpm, globval.hcorr, 1); prt_gcmat(globval.bpm, globval.vcorr, 2);
//gets response matrix, does svd, evaluates correction for N seed orbits
//get_cod_rms_scl(dx, dy, dr, n_seed)
//get_cod_rms_scl(100e-6, 100e-6, 1.0e-3, 100); //trim coils aren't reset when finished
//for alignments errors check LoadAlignTol (in nsls_ii_lib.cc) and AlignErr.dat
//CorrectCOD_N("/home/simon/projects/src/lattice/AlignErr.dat", 3, 1, 1);
//for field errors check LoadFieldErr(in nsls_ii_lib.cc) and FieldErr.dat
//LoadFieldErr("/home/simon/projects/src/lattice/FieldErr.dat", true, 1.0, true);
//for alignments errors check LoadAlignTol (in nsls_ii_lib.cc) and AlignErr.dat
//LoadAlignTol("/home/simon/projects/src/lattice/AlignErr.dat", true, 1.0, true, 1);
//LoadAlignTol("/home/simon/projects/out/20091126/AlignErr.dat", true, 1.0, true, 1);
//prt_cod("cod_err.out", globval.bpm, true); //prints a specific closed orbit with corrector strengths
// delicated for max4 lattice
//load alignment errors and field errors, correct orbit, repeat N times, and get statistics
get_cod_rms_scl_new(1); //trim coils aren't reset when finished
//for aperture limitations use LoadApers (in nsls_ii_lib.cc) and Apertures.dat
//globval.Aperture_on = true;
//LoadApers("/home/simon/projects/src/lattice/Apertures.dat", 1, 1);
}
/*******************************************************************************
Call nsls-ii_lib.cc
-----------------------------------------------------------------------------
*******************************************************************************/
//
// tune shift with amplitude
double delta = 0;
if (false) {
cout << endl;
cout << "computing tune shifts" << endl;
get_ksi2(delta); // this gets the chromas and writes them into chrom2.out
// get_ksi2(5.0e-2); // this gets the chromas and writes them into chrom2.out
}
if (false) {
//fmap(n_x, n_y, n_tr, x_max_FMA, y_max_FMA, 0.0, true, false);
// fmapdp(n_x, n_dp, n_tr, x_max_FMA, -delta_FMA, 1e-3, true, false); // always use -delta_FMA (+delta_FMA appears broken)
fmapdp(n_x, n_dp, n_tr, x_max_FMA, -delta_FMA, 1e-3, true); // always use -delta_FMA (+delta_FMA appears broken)
} else {
globval.Cavity_on = true; // this gives longitudinal motion
globval.radiation = false; // this adds ripple around long. ellipse (needs many turns to resolve damp.)
//globval.Aperture_on = true;
//LoadApers("/home/simon/projects/src/lattice/Apertures.dat", 1, 1);
//get_dynap_scl(delta, 512);
}