From 2dabfe0b6e24e474b7e262c10a104968933e39c6 Mon Sep 17 00:00:00 2001
From: zhang <zhang@9a6e40ed-f3a0-4838-9b4a-bf418f78e88d>
Date: Tue, 2 Nov 2010 15:36:23 +0000
Subject: [PATCH] 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.
---
tracy/tools/soltracy.cc | 195 ++++++++++++++++++++++++++++++++++++++--
1 file changed, 189 insertions(+), 6 deletions(-)
diff --git a/tracy/tools/soltracy.cc b/tracy/tools/soltracy.cc
index 700c31a..da45e5e 100644
--- a/tracy/tools/soltracy.cc
+++ b/tracy/tools/soltracy.cc
@@ -39,7 +39,8 @@ extern bool freq_map;
// for time handling
uint32_t start, stop;
- /************************************************************************
+
+ /************************************************************************
start read in files and flags
*************************************************************************/
if (argc > 1){
@@ -103,16 +104,84 @@ extern bool freq_map;
}
- if (FitTuneFlag == true){
- fprintf(stdout, "\n Fitting tunes: %s %s, targetnux = %f, targetnuz = %f \n",qm1,qm2,targetnux,targetnuz);
- FitTune(ElemIndex(qm1),ElemIndex(qm2), targetnux, targetnuz);
- Ring_GetTwiss(chroma=true, 0.0); /* Compute and get Twiss parameters */
+ //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);
printglob(); /* print parameter list */
}
+ /* fit the chromaticities*/
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);
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 */
}
@@ -217,10 +286,34 @@ extern bool freq_map;
// MOMENTUM ACCEPTANCE
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(
_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;
}
@@ -253,7 +346,97 @@ extern bool freq_map;
printf("the simulation time for phase space in tracy 3 is \n");
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;
}
-
--
GitLab