Newer
Older
/* Current revision $Revision$
On branch $Name$
Latest change $Date$ by $Author$
*/
#define ORDER 1
int no_tps = ORDER; // arbitrary TPSA order is defined locally
extern bool freq_map;
#include "tracy_lib.h"
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
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
void compare_cod(void)
{
const double dPmin = 1e-10, dPmax = 1e-6;
int k;
const int kmax = 20;
double dp, dpstep = dPmax/kmax;
FILE * outf;
const char fic[] = "compare_cod.out";
long lastpos = 0L;
Vector vector_findcod;
/* Opening file */
if ((outf = fopen(fic, "w")) == NULL) {
fprintf(stdout, "compare_cod: error while opening file %s\n", fic);
exit_(1);
}
fprintf(stdout,"\n");
for (k = 0; k <= kmax; k++){
dp = dPmin + k*dpstep;
getcod(dp, lastpos); // get cod for printout
CopyVec(6, globval.CODvect, vector_findcod);
findcod(dp);
// fprintf(stdout,"%14.6e %14.6e %14.6e %14.6e %14.6e %14.6e %14.6e %14.6e %14.6e\n",
// dp, globval.CODvect[0], vector_findcod[0], globval.CODvect[1],vector_findcod[1],
// globval.CODvect[2], vector_findcod[2], globval.CODvect[3], vector_findcod[3]);
fprintf(stdout,"%14.6e %14.6e %14.6e %14.6e %14.6e %14.6e %14.6e %14.6e %14.6e\n",
dp,
globval.CODvect[0], (globval.CODvect[0]-vector_findcod[0])/globval.CODvect[0],
globval.CODvect[1], (globval.CODvect[1]-vector_findcod[1])/globval.CODvect[1],
globval.CODvect[2], (globval.CODvect[2]-vector_findcod[2])/globval.CODvect[2],
globval.CODvect[3], (globval.CODvect[3]-vector_findcod[3])/globval.CODvect[3]);
}
fclose(outf);
}
#define LOG10 log(10.0)
void Getchrom2(double dP)
{
long lastpos = 0L;
double dPlocal = 0.0, expo = 0.0, TEMP = 0.0, TEMPX = 0.0, TEMPZ = 0.0;
Vector2 alpha={0.0,0.0}, beta={0.0,0.0}, gamma={0.0,0.0}, nu={0.0,0.0},
nu0={0.0,0.0}, Chrom={0.0,0.0};
trace = false;
if (dP != 0.0) {
fprintf(stderr,"Ring_Getchrom: Warning this is NOT the CHROMA, dP=%e\n",dP);
}
/* initialization */
globval.Chrom[0] = 1e38;
globval.Chrom[1] = 1e38;
expo = log(globval.dPcommon) / LOG10;
do
{
Chrom[0] = globval.Chrom[0];
Chrom[1] = globval.Chrom[1];
dPlocal = exp(expo*LOG10);
/* Get cod for energy dP - dPlocal*/
GetCOD(globval.CODimax, globval.CODeps, dP - dPlocal*0.5,
lastpos);
if (!status.codflag)
{ /* if no cod */
fprintf(stderr,"Ring_Getchrom: Lattice is unstable for dP - dPlocal=% .5e\n",
dP - dPlocal*0.5);
return;
}
/* get tunes for energy dP - dPlocal/2 from oneturn matrix */
Cell_GetABGN(globval.OneTurnMat, alpha, beta, gamma, nu0);
/* Get cod for energy dP + dPlocal*/
GetCOD(globval.CODimax, globval.CODeps, dP + dPlocal*0.5,
lastpos);
if (!status.codflag) { /* if no cod */
fprintf(stderr,"Ring_Getchrom Lattice is unstable for dP + dPlocal=% .5e \n",
dP + dPlocal*0.5);
return;
}
/* get tunes for energy dP + dPlocal/2 from oneturn matrix */
Cell_GetABGN(globval.OneTurnMat, alpha, beta, gamma, nu);
if (!globval.stable) {
printf("Ring_Getchrom: Lattice is unstable\n");
}
/* Get chromaticities by numerical differentiation*/
globval.Chrom[0] = (nu[0] - nu0[0]) / dPlocal;
globval.Chrom[1] = (nu[1] - nu0[1]) / dPlocal;
TEMP=sqrt(
fabs(globval.Chrom[0]*globval.Chrom[0]- Chrom[0]*Chrom[0])
+fabs(globval.Chrom[1]*globval.Chrom[1]- Chrom[1]*Chrom[1])
);
TEMPX = sqrt(fabs(globval.Chrom[0]*globval.Chrom[0]- Chrom[0]*Chrom[0]));
TEMPZ = sqrt(fabs(globval.Chrom[1]*globval.Chrom[1]- Chrom[1]*Chrom[1]));
// TEST CHROMA convergence
if (trace) {
fprintf(stdout,"\nexpo % e xix = % e xiz = % e TEMP = %e TEMPX %+e TEMPZ %+e\n",
expo,Chrom[0],Chrom[1],TEMP, TEMPX, TEMPZ);
fprintf(stdout,"expo % e nux = % e nuz = % e dPlocal= %+e\n",
expo,nu0[0],nu0[1],dP-0.5*dPlocal);
fprintf(stdout,"expo % e nux = % e nuz = % e dPlocal= %+e\n",
expo,nu[0],nu[1],dP+0.5*dPlocal);
}
fprintf(stdout, "%+e %+.12e %+.12e %+.12e %+.12e\n", dPlocal,
globval.Chrom[0], fabs(globval.Chrom[0]-Chrom[0])/Chrom[0],
globval.Chrom[1], fabs(globval.Chrom[1]-Chrom[1])/Chrom[1]);
expo += 0.1;
} while (expo<-2);
status.chromflag = true;
}
//***************************************************************************************
//
// 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];
/************************************************************************
start read in files and flags
*************************************************************************/
if (argc > 1){
read_script(argv[1], true);}
else{
fprintf(stdout, "Not enough parameters\nSyntax is program parameterfile\n");
return 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();
double V_RF;
V_RF = get_RFVoltage(ElemIndex("cav"));
set_RFVoltage(ElemIndex("cav"), 3e6); // 3 MV
V_RF = get_RFVoltage(ElemIndex("cav"));
// 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);
}
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 */
}
// 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(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();
}
}
/******************************************************************************************/
// COMPUTATION PART after setting the model
/******************************************************************************************/
//first print the full lattice with error as a flat file
prtmfile("flat_file_error.dat"); // writes flat file /* very important file for debug*/
// computes TuneShift with amplitudes
if (AmplitudeTuneShiftFlag == true){
if (ChamberFlag == true ){
_AmplitudeTuneShift_nypoint, _AmplitudeTuneShift_nturn,
_AmplitudeTuneShift_xmax, _AmplitudeTuneShift_ymax,
_AmplitudeTuneShift_delta);
//NuDx(31L,21L,516L,0.025,0.005,dP);
}
else{ // Utility ?
}
}
if (EnergyTuneShiftFlag == true){
if (ChamberFlag == true ){
_EnergyTuneShift_nturn, _EnergyTuneShift_deltamax);
//NuDp(31L,1026L,0.06);
}
else{ // utility ?
// Computes FMA
if (FmapFlag == true){
if (ChamberFlag == true ){
if (ExperimentFMAFlag == true)
fmap( _FmapFlag_nxpoint, _FmapFlag_nypoint,
_FmapFlag_nturn, _FmapFlag_xmax, _FmapFlag_ymax,
_FmapFlag_delta, _FmapFlag_diffusion);
//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);
}
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){
MomentumAcceptance(
_MomentumAccFlag_istart, _MomentumAccFlag_istop,
_MomentumAccFlag_deltaminp, _MomentumAccFlag_deltamaxp, _MomentumAccFlag_nstepp,
_MomentumAccFlag_deltaminn, _MomentumAccFlag_deltamaxn, _MomentumAccFlag_nstepn);
// MomentumAcceptance(1L, 108L, 0.01, 0.05, 100L, -0.01, -0.05, 100L);
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
}
// 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);
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
/*********************************************************************************
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);
}