Newer
Older
#define ORDER 1
int no_tps = ORDER; // arbitrary TPSA order is defined locally
extern bool freq_map;
#include "tracy_lib.h"
//***************************************************************************************
//
// MAIN CODE
//
//****************************************************************************************
int main(int argc, char *argv[])
{
const long seed = 1121;
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
// const double x_max_FMA = 10e-3, delta_FMA = 10e-2;
// const int n_x = 801, n_dp = 80, n_tr = 2048;
bool chroma;
double dP = 0.0;
long lastpos = -1L;
char str1[S_SIZE];
// for time handling
uint32_t start, stop;
/************************************************************************
start read in files and flags
*************************************************************************/
fprintf(stdout, "Not enough parameters\nSyntax is program parameterfile\n");
exit_(1);
/************************************************************************
end read in files and flags
*************************************************************************/
prtmfile("flat_file.dat"); // writes flat file /* very important file for debug*/
printlatt(); /* SOLEIL print out lattice functions */
printglob();
// Flag factory
//set RF voltage
if (RFvoltageFlag == true){
printf("\nSetting RF voltage:\n");
printf(" Old RF voltage is: %lf [MV]\n",get_RFVoltage(ElemIndex("cav"))/1e6);
set_RFVoltage(ElemIndex("cav"), RFvolt);
printf(" New RF voltage is: %lf [MV]\n",get_RFVoltage(ElemIndex("cav"))/1e6);
}
// 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 "Apertures.dat" , soleil version*/
PrintCh(); // print chamber into chamber.out
//get_matching_params_scl(); // get tunes and beta functions at entrance
get_alphac2(); // compute up to 3rd order mcf
//cout << endl << "computing tune shifts" << endl;
//dnu_dA(10e-3, 5e-3, 0.002);
//get_ksi2(0.0); // this gets the chromas and writes them into chrom2.out
// 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){
start = stampstart();
GetChromTrac(2L, 1026L, 1e-5, &ksix, &ksiz);
stop = stampstop(start);
fprintf(stdout,"From tracking: ksix= % f ksiz= % f \n",ksix,ksiz);
}
107
108
109
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
//generic function, to fit tunes using 1 family of quadrupoles
if (FitTuneFlag == true){
double qf_bn = 0.0,qd_bn = 0.0;
double qf_bn_fit = 0.0, qd_bn_fit = 0.0;
/* quadrupole field strength before fitting*/
qf_bn = Cell[Elem_GetPos(ElemIndex(qf), 1)].Elem.M->PB[HOMmax+2];
qd_bn = Cell[Elem_GetPos(ElemIndex(qd), 1)].Elem.M->PB[HOMmax+2];
/* fitting tunes*/
fprintf(stdout, "\n Fitting tunes: %s %s, targetnux = %f, targetnuz = %f \n",qf,qd,targetnux,targetnuz);
FitTune(ElemIndex(qf),ElemIndex(qd),targetnux, targetnuz);
/* integrated field strength after fitting*/
qf_bn_fit = Cell[Elem_GetPos(ElemIndex(qf), 1)].Elem.M->PB[HOMmax+2];
qd_bn_fit = Cell[Elem_GetPos(ElemIndex(qd), 1)].Elem.M->PB[HOMmax+2];
/* print out the quadrupole strengths before and after the fitting*/
printf("Before the fitting, the quadrupole field strength is: \n");
printf(" %s = %f, %s = %f\n",qf,qf_bn,qd,qd_bn);
printf("After the fitting, the quadrupole field strength is: \n");
printf(" %s = %f, %s = %f\n",qf,qf_bn_fit,qd,qd_bn_fit);
/* Compute and get Twiss parameters */
Ring_GetTwiss(chroma=true, 0.0);
printglob(); /* print parameter list */
}
/* sepcific for soleil ring in which the quadrupole is cut into 2 parts*/
if (FitTune4Flag == true){
double qf1_bn = 0.0, qf2_bn = 0.0,qd1_bn = 0.0,qd2_bn = 0.0;
double qf1_bn_fit = 0.0, qf2_bn_fit = 0.0,qd1_bn_fit = 0.0,qd2_bn_fit = 0.0;
/* quadrupole field strength before fitting*/
qf1_bn = Cell[Elem_GetPos(ElemIndex(qf1), 1)].Elem.M->PB[HOMmax+2];
qf2_bn = Cell[Elem_GetPos(ElemIndex(qf2), 1)].Elem.M->PB[HOMmax+2];
qd1_bn = Cell[Elem_GetPos(ElemIndex(qd1), 1)].Elem.M->PB[HOMmax+2];
qd2_bn = Cell[Elem_GetPos(ElemIndex(qd2), 1)].Elem.M->PB[HOMmax+2];
/* fitting tunes*/
fprintf(stdout, "\n Fitting tunes for Soleil ring: %s %s %s %s, targetnux = %f, targetnuz = %f \n",qf1,qf2,qd1,qd2,targetnux,targetnuz);
FitTune4(ElemIndex(qf1),ElemIndex(qf2),ElemIndex(qd1),ElemIndex(qd2),targetnux, targetnuz);
/* integrated field strength after fitting*/
qf1_bn_fit = Cell[Elem_GetPos(ElemIndex(qf1), 1)].Elem.M->PB[HOMmax+2];
qf2_bn_fit = Cell[Elem_GetPos(ElemIndex(qf2), 1)].Elem.M->PB[HOMmax+2];
qd1_bn_fit = Cell[Elem_GetPos(ElemIndex(qd1), 1)].Elem.M->PB[HOMmax+2];
qd2_bn_fit = Cell[Elem_GetPos(ElemIndex(qd2), 1)].Elem.M->PB[HOMmax+2];
/* print out the quadrupole strengths before and after the fitting*/
printf("Before the fitting, the quadrupole field strength is: \n");
printf(" %s = %f, %s = %f, %s = %f, %s = %f\n",qf1,qf1_bn, qf2,qf2_bn,qd1,qd1_bn,qd2,qd2_bn);
printf("After the fitting, the quadrupole field strength is: \n");
printf(" %s = %f, %s = %f, %s = %f, %s = %f\n",qf1,qf1_bn_fit, qf2,qf2_bn_fit,qd1,qd1_bn_fit,qd2,qd2_bn_fit);
/* Compute and get Twiss parameters */
Ring_GetTwiss(chroma=true, 0.0);
double sxm1_bn = 0.0, sxm2_bn = 0.0;
double sxm1_bn_fit = 0.0, sxm2_bn_fit = 0.0;
sxm1_bn = Cell[Elem_GetPos(ElemIndex(sxm1), 1)].Elem.M->PB[HOMmax+3];
sxm2_bn = Cell[Elem_GetPos(ElemIndex(sxm2), 1)].Elem.M->PB[HOMmax+3];
/* fit the chromaticites*/
fprintf(stdout, "\n Fitting chromaticities: %s %s, targetksix = %f, targetksiz = %f\n",sxm1,sxm2,targetksix,targetksiz);
FitChrom(ElemIndex(sxm1),ElemIndex(sxm2), targetksix, targetksiz);
sxm1_bn_fit = Cell[Elem_GetPos(ElemIndex(sxm1), 1)].Elem.M->PB[HOMmax+3];
sxm2_bn_fit = Cell[Elem_GetPos(ElemIndex(sxm2), 1)].Elem.M->PB[HOMmax+3];
/* print out the sextupole strengths before and after the fitting*/
printf("Before the fitting, the sextupole field strength is: \n");
printf(" %s = %f, %s = %f\n",sxm1,sxm1_bn, sxm2,sxm2_bn);
printf("After the fitting, the sextupole field strength is: \n");
printf(" %s = %f, %s = %f\n",sxm1,sxm1_bn_fit, sxm2,sxm2_bn_fit);
Ring_GetTwiss(chroma=true, 0.0); /* Compute and get Twiss parameters */
printglob(); /* print parameter list */
}
// coupling calculation
if (CouplingFlag == true){
Ring_GetTwiss(chroma=true, 0.0); /* Compute and get Twiss parameters */
printlatt(); /* dump linear lattice functions into "linlat.dat" */
GetEmittance(ElemIndex("cav"), true);
Coupling_Edwards_Teng();
printglob(); /* print parameter list */
}
// add coupling by random rotating of the quadrupoles
if (ErrorCouplingFlag == true){
Ring_GetTwiss(chroma=true, 0.0); /* Compute and get Twiss parameters */
printlatt(); /* dump linear lattice functions into "linlat.dat" */
Coupling_Edwards_Teng();
GetEmittance(ElemIndex("cav"), true);
printglob(); /* print parameter list */
}
// WARNING Fit tunes and chromaticities before applying errors !!!!
//set multipoles in all magnets
// read multipole error from a file
if (ReadMultipoleFlag == true){
fprintf(stdout, "\n Read Multipoles file for lattice with thick sextupoles \n");
ReadFieldErr(multipole_file);
}
if (MultipoleFlag == true ){
if (ThinsextFlag ==true){
fprintf(stdout, "\n Setting Multipoles for lattice with thin sextupoles \n");
Multipole_thinsext(fic_hcorr,fic_vcorr,fic_skew); /* 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(fic_hcorr,fic_vcorr,fic_skew); /* for thick sextupoles */
Ring_GetTwiss(chroma=true, 0.0); /* Compute and get Twiss parameters */
printglob();
}
}
//first print the full lattice with error as a flat file
prtmfile("flat_file_error.dat"); // writes flat file /* very important file for debug*/
printlatt(); /* SOLEIL print out lattice functions */
printglob();
/******************************************************************************************/
// COMPUTATION PART after setting the model
/******************************************************************************************/
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
// computes TuneShift with amplitudes
if (AmplitudeTuneShiftFlag == true){
if (ChamberFlag == true ){
NuDx(_AmplitudeTuneShift_nxpoint,
_AmplitudeTuneShift_nypoint, _AmplitudeTuneShift_nturn,
_AmplitudeTuneShift_xmax, _AmplitudeTuneShift_ymax,
_AmplitudeTuneShift_delta);
//NuDx(31L,21L,516L,0.025,0.005,dP);
}
else{ // Utility ?
NuDx(50L,30L,516L,0.035,0.02,dP);
}
}
if (EnergyTuneShiftFlag == true){
if (ChamberFlag == true ){
NuDp(_EnergyTuneShift_npoint,
_EnergyTuneShift_nturn, _EnergyTuneShift_deltamax);
//NuDp(31L,1026L,0.06);
}
else{ // utility ?
NuDp(31L,1026L,0.06);
}
}
fmap( _FmapFlag_nxpoint, _FmapFlag_nypoint,
_FmapFlag_nturn, _FmapFlag_xmax, _FmapFlag_ymax,
_FmapFlag_delta, _FmapFlag_diffusion);
// Compute FMA dp
if (FmapdpFlag == true){
fmapdp( _FmapdpFlag_nxpoint, _FmapdpFlag_nepoint,
_FmapdpFlag_nturn, _FmapdpFlag_xmax, _FmapdpFlag_emax,
_FmapdpFlag_z, _FmapdpFlag_diffusion);
}
if (CodeComparaisonFlag){
fmap(200,100,1026,-32e-3,7e-3,0.0,true);
}
bool cavityflag, radiationflag;
/* record the initial values*/
cavityflag = globval.Cavity_on;
radiationflag = globval.radiation;
if(strncmp("6D",DimTrack,2)==0){
globval.Cavity_on = true;
globval.radiation = true;
}
else if(strncmp("4D",DimTrack,2)==0){
globval.Cavity_on = false;
globval.radiation = false;
}
else{
printf("MomentumAccFlag: Error!!! DimTrack must be '4D' or '6D'\n");
exit_(1);
};
MomentumAcceptance(
_MomentumAccFlag_istart, _MomentumAccFlag_istop,
_MomentumAccFlag_deltaminp, _MomentumAccFlag_deltamaxp, _MomentumAccFlag_nstepp,
_MomentumAccFlag_deltaminn, _MomentumAccFlag_deltamaxn, _MomentumAccFlag_nstepn);
/* reset back to the initial values*/
globval.Cavity_on = cavityflag;
globval.radiation = radiationflag;
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
}
// 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){
start = stampstart();
Phase(0.001,0,0.001,0,0.001,0, 1);
printf("the simulation time for phase space in tracy 3 is \n");
stop = stampstop(start);
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
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
// ????????????? NSRL version, Check with Soleil version "MomentumAcceptance"
// IBS & TOUSCHEK
int k, n_turns;
double sigma_s, sigma_delta, tau, alpha_z, beta_z, gamma_z, eps[3];
FILE *outf;
const double Qb = 5e-9; // e charge in one bunch
if (TouschekFlag == true) {
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;
}
globval.eps[Y_] = 0.008e-9;
n_turns = 40;
// get the twiss parameters
alpha_z =
-globval.Ascr[ct_][ct_]*globval.Ascr[delta_][ct_]
-globval.Ascr[ct_][delta_]*globval.Ascr[delta_][delta_];
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
}
Touschek(Qb, globval.delta_RF, globval.eps[X_], globval.eps[Y_],
sigma_delta, sigma_s);
// Intra Beam Scattering(IBS)
if (IBSFlag == true) {
// 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) {
globval.Aperture_on = 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);
}
}