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;
/*The following values are set in the Read_lattice( ) in soleilcommon.cc*/
// globval.quad_fringe = false; // quadrupole fringe field
// globval.radiation = false; // synchrotron radiation
// globval.emittance = false; // emittance
// globval.pathlength = false;
// 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;
/************************************************************************
*************************************************************************/
fprintf(stdout, "Not enough parameters\nSyntax is program parameterfile\n");
exit_(1);
/************************************************************************
writes flat file with all the design values of the lattice, very important file for debug
*************************************************************************/
prtmfile("flat_file.dat");
// 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);
}
103
104
105
106
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
//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);
//first print the full lattice with error as a flat file
prtmfile("flat_file_error.dat"); // writes flat file /* very important file for debug*/
Ring_GetTwiss(chroma=true, 0.0); /* Compute and get Twiss parameters */
printglob();
//Old version to multipole errors in soleil lattice, is obsoleted
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 */
prtmfile("flat_file_error.dat"); // writes flat file /* very important file for debug*/
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 */
prtmfile("flat_file_error.dat"); // writes flat file /* very important file for debug*/
Ring_GetTwiss(chroma=true, 0.0); /* Compute and get Twiss parameters */
printglob();
}
}
/******************************************************************************************/
// COMPUTATION PART after setting the model
/******************************************************************************************/
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
// 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;
/* set the dimension for the momentum tracking*/
if(strncmp("6D",TrackDim,2)==0){
else if(strncmp("4D",TrackDim,2)==0){
globval.Cavity_on = false;
globval.radiation = false;
}
else{
printf("MomentumAccFlag: Error!!! DimTrack must be '4D' or '6D'\n");
exit_(1);
};
/* calculate momentum accepetance*/
MomentumAcceptance(
_MomentumAccFlag_istart, _MomentumAccFlag_istop,
_MomentumAccFlag_deltaminp, _MomentumAccFlag_deltamaxp, _MomentumAccFlag_nstepp,
_MomentumAccFlag_deltaminn, _MomentumAccFlag_deltamaxn, _MomentumAccFlag_nstepn);
globval.Cavity_on = cavityflag;
globval.radiation = radiationflag;
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
344
345
}
// 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);
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
440
441
// ????????????? 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);
}
}