Skip to content
Snippets Groups Projects
Commit 2dabfe0b authored by zhang's avatar zhang
Browse files

02/11/2010 

1) Call functions to fit tunes for soleil ring.
2) Print the quadrupole and sextupoles values before and after the fitting to the screen.
3) Add the routines to calculate Touschek lifetime.
parent a30f25f4
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 189 additions and 6 deletions
tracy/tools/soltracy.cc
+ 189
− 6
View file @ 2dabfe0b
...@@ -39,7 +39,8 @@ extern bool freq_map; ...@@ -39,7 +39,8 @@ extern bool freq_map;
// for time handling // for time handling
uint32_t start, stop; uint32_t start, stop;
/************************************************************************
/************************************************************************
start read in files and flags start read in files and flags
*************************************************************************/ *************************************************************************/
if (argc > 1){ if (argc > 1){
...@@ -103,16 +104,84 @@ extern bool freq_map; ...@@ -103,16 +104,84 @@ extern bool freq_map;
} }
if (FitTuneFlag == true){ //generic function, to fit tunes using 1 family of quadrupoles
fprintf(stdout, "\n Fitting tunes: %s %s, targetnux = %f, targetnuz = %f \n",qm1,qm2,targetnux,targetnuz); if (FitTuneFlag == true){
FitTune(ElemIndex(qm1),ElemIndex(qm2), targetnux, targetnuz); double qf_bn = 0.0,qd_bn = 0.0;
Ring_GetTwiss(chroma=true, 0.0); /* Compute and get Twiss parameters */ 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);
printglob(); /* print parameter list */ printglob(); /* print parameter list */
} }
/* fit the chromaticities*/
if (FitChromFlag == true){ if (FitChromFlag == true){
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); fprintf(stdout, "\n Fitting chromaticities: %s %s, targetksix = %f, targetksiz = %f\n",sxm1,sxm2,targetksix,targetksiz);
FitChrom(ElemIndex(sxm1),ElemIndex(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 */ Ring_GetTwiss(chroma=true, 0.0); /* Compute and get Twiss parameters */
printglob(); /* print parameter list */ printglob(); /* print parameter list */
} }
...@@ -217,10 +286,34 @@ extern bool freq_map; ...@@ -217,10 +286,34 @@ extern bool freq_map;
// MOMENTUM ACCEPTANCE // MOMENTUM ACCEPTANCE
if (MomentumAccFlag == true){ if (MomentumAccFlag == 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( MomentumAcceptance(
_MomentumAccFlag_istart, _MomentumAccFlag_istop, _MomentumAccFlag_istart, _MomentumAccFlag_istop,
_MomentumAccFlag_deltaminp, _MomentumAccFlag_deltamaxp, _MomentumAccFlag_nstepp, _MomentumAccFlag_deltaminp, _MomentumAccFlag_deltamaxp, _MomentumAccFlag_nstepp,
_MomentumAccFlag_deltaminn, _MomentumAccFlag_deltamaxn, _MomentumAccFlag_nstepn); _MomentumAccFlag_deltaminn, _MomentumAccFlag_deltamaxn, _MomentumAccFlag_nstepn);
/* reset back to the initial values*/
globval.Cavity_on = cavityflag;
globval.radiation = radiationflag;
} }
...@@ -253,7 +346,97 @@ extern bool freq_map; ...@@ -253,7 +346,97 @@ extern bool freq_map;
printf("the simulation time for phase space in tracy 3 is \n"); printf("the simulation time for phase space in tracy 3 is \n");
stop = stampstop(start); stop = stampstop(start);
} }
// ????????????? 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);
}
}
return 0; return 0;
} }
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

The place to be for your source code