diff --git a/tracy/tracy/src/soleillib.cc b/tracy/tracy/src/soleillib.cc
index 37d2ec604ff049bd13074d61994043e9f2de4ade..836417e7a9c0d5db8aab467f3bbb351a572665b6 100644
--- a/tracy/tracy/src/soleillib.cc
+++ b/tracy/tracy/src/soleillib.cc
@@ -273,7 +273,7 @@ void Hcofonction(long pos, double dP)
/****************************************************************************/
-/* void SetErr(void)
+/* void SetErr(long seed,double fac)
Purpose:
Set error
@@ -282,8 +282,8 @@ void Hcofonction(long pos, double dP)
Distribution gaussienne d'ecart type fac et coupe a normcut*sigma
Input:
- none
-
+ seed: random seed number
+ fac : RMS value of the rotation angle of the quadrupole
Output:
none
@@ -305,29 +305,22 @@ void Hcofonction(long pos, double dP)
Test if normal quad sin(theta) = 0. Do not work if tilt error
****************************************************************************/
-void SetErr(void)
+void SetErr(long seed,double fac)
{
- double fac = 0.0, normcut = 0.0;
- long seed = 0L;
+ double normcut = 0.0;
long i = 0L;
// CellType Cell;
double theta = 0.0;
int pair = 0;
-
- /* coupling */
-// fac = 0.00064; //normal faux
-// fac = 0.00075; //normal faux
-// fac = 0.00155; //normal faux
-// fac = 0.00076; // mais coupe en deux
-// fac = 0.001335/2.0;
-
-// Normal
-// fac = 0.00155/2.0;
-// setrancut(normcut=2L);
-// iniranf(seed);
-
- fac = 0.0014/2.0;
-// fac = 0.00115/2.0;
+ bool prt=false;
+
+
+ if(!prt){
+ printf("\n");
+ printf(" Setting random rotation errors to quadrupole:\n");
+ printf(" random seed number is: %ld, rms value of the rotation error is: %lf \n",seed,fac);
+ }
+
setrancut(normcut=2L);
iniranf(seed);
@@ -564,11 +557,11 @@ void DefineChNoU20(void)
1) line begin with "#" is comment line
2) first line Name1: Start
first line Name2: All
- 3) Name1 and Name2 should be unique in the lattice
- 4) Name1 is defined before Name2 in the lattice
+ 3) MK1 and MK2 should be unique in the lattice
+ 4) MK1 is defined before MK2 in the lattice
5)
- Name1: start element of the section for the aperture
- Name2: end element of the section for the aperture
+ MK1: marker before the start element of the section for the aperture
+ Mk2: marker after the end element of the section for the aperture
dxmin: minimum x value of vacuum chamber
dxmax: maxmum x value of vacuum chamber
dymin: minimum y value of vacuum chamber
@@ -602,13 +595,12 @@ void ReadCh(const char *AperFile)
int i=0, LineNum=0;
double dxmin=0.0, dxmax=0.0, dymin=0.0, dymax=0.0; // min and max x and apertures
FILE *fp;
-
bool prt = false;
fp = file_read(AperFile);
printf("\n");
- printf("...Load and Set Apertures to lattice elements\n");
+ printf("Loading and setting vacuum apertures to lattice elements...\n");
while (fgets(line, max_str, fp) != NULL) {
LineNum++; /* count the line number that has been read*/
@@ -657,6 +649,7 @@ void ReadCh(const char *AperFile)
k1 = Elem_GetPos(Fnum1, 1);
k2 = Elem_GetPos(Fnum2, 1);
/* set the vacuum chamber*/
+ //read the marker before the first element, and the markder after the last elment
for(i=1; i<globval.Cell_nLoc; i++){
if(i>=k1 && i<k2){
Cell[i].maxampl[X_][0] = dxmin;
@@ -675,8 +668,8 @@ void ReadCh(const char *AperFile)
}
}
- else /* print the comment line */
- printf("%s", line);
+ // else /* print the comment line */
+ // printf("%s", line);
}
fclose(fp);
}
@@ -1804,7 +1797,7 @@ void Multipole_thicksext(char const *fic_hcorr, char const *fic_vcorr, char cons
double b2 = 0.0, b3 = 0.0;
double dBoB2 = 0.0, dBoB3 = 0.0, dBoB4 = 0.0, dBoB5 = 0.0, dBoB6 = 0.0,
dBoB7 = 0.0, dBoB9 = 0.0, dBoB11 = 0.0, dBoB15 = 0.0, dBoB21 = 0.0,
- dBoB27;
+ dBoB27 = 0.0;
double dBoB6C = 0.0, dBoB6L = 0.0, dBoB10C = 0.0, dBoB10L = 0.0,
dBoB14C = 0.0, dBoB14L = 0.0, dBoB3C = 0.0, dBoB3L = 0.0,
dBoB4C = 0.0, dBoB4L = 0.0;
@@ -1991,7 +1984,7 @@ void Multipole_thicksext(char const *fic_hcorr, char const *fic_vcorr, char cons
x06i = x03i*x03i;
x07i = x06i*x0i;
- dBoB2 = 2.2e-4*0; /* gradient, used for curve trajectory simulation */
+ dBoB2 = 2.2e-4*1; /* gradient, used for curve trajectory simulation */
dBoB3 = -3.0e-4*1; /* hexapole */
dBoB4 = 2.0e-5*1; /* octupole */
dBoB5 = -1.0e-4*1; /* decapole */
@@ -2145,8 +2138,8 @@ void Multipole_thicksext(char const *fic_hcorr, char const *fic_vcorr, char cons
/* multipoles from dipolar unallowed component */
dBoB5 = 5.4e-4*1;
dBoB7 = 3.3e-4*1;
- dBoB5rms = 4.7e-4*1;
- dBoB7rms = 2.1e-4*1;
+ dBoB5rms = 4.7e-4*1; // for test
+ dBoB7rms = 2.1e-4*1; // for test
/* allowed multipoles */
dBoB9 = -4.7e-4*1;
@@ -2367,7 +2360,7 @@ void Multipole_thicksext(char const *fic_hcorr, char const *fic_vcorr, char cons
x0i = 1.0/35e-3; /* 1/radius */
x02i = x0i*x0i;
- dBoB4 = -0.680*0; /* Octupole */
+ dBoB4 = -0.680*1; /* Octupole */
/* open skew quaI (A) *
310
@@ -5033,3 +5026,558 @@ void Enveloppe2(double x, double px, double y, double py, double dp, double ntur
}
}
+/****************************************************************************/
+/* double get_RFVoltage(const int Fnum)
+
+ Purpose:
+ Get RF voltage of family Fnum
+
+ Input:
+ Fnum: family name
+
+ Output:
+ none
+
+ Return:
+ RF voltage
+
+ Global variables:
+ none
+
+ Specific functions:
+ none
+
+ Comments:
+ 10/2010 by L.Nadolski
+****************************************************************************/
+double get_RFVoltage(const int Fnum){
+
+ double V_RF = 0.0;
+ bool prt = false;
+
+ V_RF = Cell[Elem_GetPos(Fnum, 1)].Elem.C->Pvolt; //RF voltage in Volts
+ if (prt) fprintf(stdout, "RF voltage of cavity %s is %f MV \n",
+ Cell[Elem_GetPos(Fnum, 1)].Elem.PName, V_RF/1e6);
+ return V_RF;
+}
+
+/****************************************************************************/
+/* void set_RFVoltage(const int Fnum, const double V_RF)
+
+ Purpose:
+ Set RF voltage to the first kid in the family Fnum
+
+ Input:
+ Fnum: family name
+
+ Output:
+ none
+
+ Return:
+ RF voltage
+
+ Global variables:
+ none
+
+ Specific functions:
+ none
+
+ Comments:
+ 10/2010 by L.Nadolski
+****************************************************************************/
+void set_RFVoltage(const int Fnum, const double V_RF){
+
+ int k, n = 0;
+
+
+ n = GetnKid(Fnum);
+ bool prt = false;
+
+ for (k=1; k <=n; k++){
+ Cell[Elem_GetPos(Fnum, k)].Elem.C->Pvolt = V_RF; // in Volts
+ }
+ if(prt)
+ fprintf(stdout, "Setting cavity %s to %f MV \n",
+ Cell[Elem_GetPos(Fnum, 1)].Elem.PName, V_RF/1e6);
+}
+
+
+/****************************************************************************************************/
+/* void ReadFieldErr(const char *FieldErrorFile)
+
+ Purpose:
+ Read multipole errors from a file
+ The input format of the file is:
+
+ keyWords raduis when the field error is meausred "r0", field error order "n",
+ field error component "Bn", field error component "An"; "n","Bn,""An",...
+
+ Input:
+
+
+ Output:
+ none
+
+ Return:
+
+
+ Global variables:
+ none
+
+ Specific functions:
+ none
+
+ Comments:
+ 10/2010 Written by Jianfeng Zhang
+*****************************************************************************************************/
+void ReadFieldErr(const char *FieldErrorFile)
+{
+ char line[max_str], name[max_str], *prm;
+ int n = 0; /* field error order*/
+ double Bn = 0.0, An = 0.0, r0 = 0.0; /* field error components and radius when the field error is measured */
+ /* conversion number from A to T.m for soleil*/
+ double _convHcorr = 8.14e-4,_convVcorr = 4.642e-4, _convQt = 93.83e-4;
+ FILE *inf;
+
+ const bool prt = false;
+
+ inf = file_read(FieldErrorFile);
+
+ printf("\n");
+ /* read lines*/
+ while (fgets(line, max_str, inf) != NULL) {
+ /* check the line is whether comment line or null line*/
+ if (strstr(line, "#") == NULL && strcmp(line,"\n") != 0) {
+ /*read and assign the key words and measure radius*/
+ sscanf(line, "%s", name);
+ sscanf(line, " %*s %lf", &r0);
+ printf("\nsetting multipole error to: %-5s r0 = %7.1le\n",name,r0);
+ // skip first two parameters
+ strtok(line, " \t");
+ strtok(NULL, " \t");
+
+ while (((prm = strtok(NULL, " \t")) != NULL) &&
+ (strcmp(prm, "\n") != 0)) {
+ /* read and assign n Bn An*/
+ sscanf(prm, "%d", &n);
+ prm = get_prm(); /*move the pointer to the next block of the line, delimiter is table key */
+ sscanf(prm, "%lf", &Bn);
+ prm = get_prm();
+ sscanf(prm, "%lf", &An);
+
+ if (!prt)
+ printf(" n = %2d, Bn = %9.1e, An = %9.1e\n", n, Bn, An);
+
+ /* set multipole error to horizontal correctors of soleil ring*/
+ if(strcmp("hcorr", name) == 0)
+ SetCorrQtErr_fam(fic_hcorr,globval.hcorr,_convHcorr,r0,n,Bn,An);
+ /* set multipole error to vertical correctors of soleil ring*/
+ else if(strcmp("vcorr", name) == 0)
+ SetCorrQtErr_fam(fic_vcorr,globval.vcorr,_convVcorr,r0,n,Bn,An);
+ /* set multipole error to skew quadrupoles of soleil ring*/
+ else if(strcmp("qt", name) == 0)
+ SetCorrQtErr_fam(fic_skew,globval.qt,_convQt,r0,n,Bn,An);
+ else
+ /* set errors for other multipole*/
+ AddFieldErrors(name, r0, n, Bn, An) ;
+ }
+ } else
+ continue;
+ // printf("%s", line);
+ }
+
+ fclose(inf);
+}
+
+/***********************************************************************
+void AddFieldErrors(const char *name, const double r0,
+ const int n, const double Bn, const double An)
+
+ Purpose:
+ Add field error of the elements with the same type or single element,
+ with the previous value, and then the summation value replaces
+ the previous value.
+
+ Input:
+ name type name or element name
+ r0 radius at which error is measured, error field is relative
+ to the design field strength when r0 !=0
+ n order of the error
+ Bn relative B component for the n-th error
+ An relative A component for the n-th error
+
+
+ Output:
+ None
+
+ Return:
+ None
+
+ Global variables
+ None
+
+ Specific functions:
+ None
+
+ Comments:
+ 10/2010 Written by Jianfeng Zhang
+**********************************************************************/
+void AddFieldErrors(const char *name, const double r0,
+ const int n, const double Bn, const double An)
+{
+ int Fnum = 0;
+
+ if (strcmp("all", name) == 0) {
+ printf("all: not yet implemented\n");
+ } else if (strcmp("dip", name) == 0) {
+ AddFieldValues_type(Dip, r0, n, Bn, An);
+ } else if (strcmp("quad", name) == 0) {
+ AddFieldValues_type(Quad, r0, n, Bn, An);
+ } else if (strcmp("sext", name) == 0) {
+ AddFieldValues_type(Sext, r0, n, Bn, An);
+ } else {
+ Fnum = ElemIndex(name);
+ if(Fnum > 0)
+ AddFieldValues_fam(Fnum, r0, n, Bn, An);
+ else
+ printf("SetFieldErrors: undefined element %s\n", name);
+ }
+}
+
+
+/***********************************************************************
+void SetFieldValues_type(const int N, const double r0,
+ const int n, const double Bn, const double An)
+
+ Purpose:
+ Add the field error of the upright multipole with the design order "type"
+ with the previous value, and then the summation value replaces the previous value.
+ Input:
+ N type name
+ r0 radius at which error is measured, error field is relative
+ to the design field strength when r0 != 0
+ n order of the error
+ Bn relative B component of n-th error
+ An relative A component of n-th error
+
+
+
+ Output:
+ None
+
+ Return:
+ None
+
+ Global variables
+ None
+
+ Specific functions:
+ None
+
+ Comments:
+ 14/10/2010 Written by Jianfeng Zhang
+**********************************************************************/
+void AddFieldValues_type(const int N, const double r0,
+ const int n, const double Bn, const double An)
+{
+ double bnL = 0.0, anL = 0.0, KLN = 0.0;
+ int k = 0;
+
+
+ if (r0 == 0.0) {
+ // input is: (b_n*L), (a_n*L) ????
+ set_bnL_sys_type(N, n, Bn, An);
+ } else {
+ // find the strength for multipole
+ for(k = 1; k <= globval.Cell_nLoc; k++)
+ {
+ //only set upright multipole, NOT set skew multipole(skew quadrupole,etc)
+ if ((Cell[k].Elem.Pkind == Mpole) && Cell[k].Elem.M->n_design == N && Cell[k].Elem.M->PdTpar == 0)
+ {
+ //find the integrated design field strength
+ if(N == 1)
+ KLN = Cell[k].Elem.PL*Cell[k].Elem.M->Pirho; /*dipole angle*/
+ else
+ KLN = GetKLpar(Cell[k].Fnum, Cell[k].Knum, N);/*other multipoles*/
+
+
+ //absolute integrated multipole error strength
+ bnL = Bn/pow(r0, n-N)*KLN;
+ anL = An/pow(r0, n-N)*KLN;
+
+
+ //NOT add the multipole errors of short quadrupole to long quadrupole qp2 & qp7 of soleil ring
+ if(N == 2 && strncmp(Cell[k].Elem.PName,"qp2",3)==0)
+ add_bnL_sys_elem(Cell[k].Fnum, Cell[k].Knum, n, 0, 0);
+ else if(N == 2 && strncmp(Cell[k].Elem.PName,"qp7",3)==0)
+ add_bnL_sys_elem(Cell[k].Fnum, Cell[k].Knum, n, 0, 0);
+ else
+ //add errors to multipoles except qp2, qp7
+ add_bnL_sys_elem(Cell[k].Fnum, Cell[k].Knum, n, bnL,anL);
+
+ }
+ }
+ }
+
+}
+/***********************************************************************
+void AddFieldValues_fam(const int Fnum, const double r0,
+ const int n, const double Bn, const double An)
+
+ Purpose:
+ add field error of all the kids in a family, with the previous value,
+ and then the summation value replaces the previous value.
+
+ Input:
+ Fnum family name
+ r0 radius at which error is measured
+ n order of the error
+ Bn relative B component for the n-th error
+ An relative A component for the n-th error
+ new_rnd bool flag to set new random number
+
+
+ Output:
+ None
+
+ Return:
+ None
+
+ Global variables
+ None
+
+ Specific functions:
+ None
+
+ Comments:
+ 10/2010 Written by Jianfeng Zhang
+**********************************************************************/
+void AddFieldValues_fam(const int Fnum, const double r0,
+ const int n, const double Bn, const double An)
+{
+ int loc = 0, ElemNum = 0, N = 0, k = 0;
+ double bnL = 0.0, anL = 0.0, KLN = 0.0;
+
+ loc = Elem_GetPos(Fnum, 1); /*element index of first kid*/
+ N = Cell[loc].Elem.M->n_design;/*design field order*/
+
+ if (r0 == 0.0) {
+ if(Cell[loc].Elem.PL != 0){
+ bnL = Cell[loc].Elem.PL*Bn;
+ anL = Cell[loc].Elem.PL*An;
+ }
+ else{
+ bnL = Bn;
+ anL = An;
+ }
+ // input is: (b_n*L), (a_n*L)
+ set_bnL_sys_fam(Fnum, n, bnL, anL);
+ } else {
+ // find the integrated design field strength for multipole
+ if (Cell[loc].Elem.M->n_design == 1)
+ KLN = Cell[loc].Elem.PL*Cell[loc].Elem.M->Pirho; /* dipole angle */
+ else
+ KLN = GetKLpar(Cell[loc].Fnum, Cell[loc].Knum, N);/* other multipole*/
+
+ /* absolute integrated field strength*/
+ bnL = Bn/pow(r0, n-N)*KLN;
+ anL = An/pow(r0, n-N)*KLN;
+
+ //add absolute multipole field error for the family
+ for(k = 1; k <= GetnKid(Fnum); k++){
+ ElemNum = Elem_GetPos(Fnum, k); /*get the element index*/
+ add_bnL_sys_elem(Cell[ElemNum].Fnum, Cell[ElemNum].Knum, n, bnL, anL);
+ }
+ }
+}
+
+
+/***********************************************************************
+void add_bnL_sys_elem(const int Fnum, const int Knum,
+ const int n, const double bnL, const double anL)
+
+ Purpose:
+ Add the field error with the previous value, then
+ the summmation value replace the previous value,
+ in the PBsys definition of multipole.
+
+ Input:
+ Fnum family index
+ Knum kids index
+ n order of the error
+ bnL absolute integrated B component for the n-th error
+ anL absolute integrated A component for the n-th error
+
+
+
+ Output:
+ None
+
+ Return:
+ None
+
+ Global variables
+ None
+
+ Specific functions:
+ None
+
+ Comments:
+ 10/2010 Written Jianfeng Zhang
+**********************************************************************/
+void add_bnL_sys_elem(const int Fnum, const int Knum,
+ const int n, const double bnL, const double anL)
+{
+ elemtype elem;
+ const bool prt = false;
+
+ elem = Cell[Elem_GetPos(Fnum, Knum)].Elem;
+
+ if (elem.PL != 0.0) {
+ elem.M->PBsys[HOMmax+n] += bnL/elem.PL;
+ elem.M->PBsys[HOMmax-n] += anL/elem.PL;
+ } else {
+ // thin kick
+ elem.M->PBsys[HOMmax+n] += bnL; elem.M->PBsys[HOMmax-n] += anL;
+ }
+
+ Mpole_SetPB(Fnum, Knum, n); //set for Bn component
+ Mpole_SetPB(Fnum, Knum, -n); //set for An component
+
+ if (prt)
+ printf("add the n=%d component of %s, bnL = %e, %e, anL = %e, %e\n",
+ n, Cell[Elem_GetPos(Fnum, Knum)].Elem.PName,
+ bnL, elem.M->PB[HOMmax+n],anL, elem.M->PB[HOMmax-n]);
+}
+
+/***********************************************************************
+void SetCorrQtErr_fam(char const *fic, const int Fnum, const double conv, const double r0,
+ const int n, const double Bn, const double An)
+
+ Purpose:
+ Set multipole field error to the thick sextupole which also functions as
+ skew quadrupoles, horizontal and vertical correctors which are used for
+ orbit correction.
+
+ Input:
+ fic file name with measured corrector value or qt values
+ Fnum family index of horizontal or vertical corrector or skew quadrupole
+ conv conversion from A to T.m for soleil
+ r0 radius at which error is measured
+ n order of the error
+ Bn integrated B component for the n-th error
+ An integrated A component for the n-th error
+
+ Output:
+ None
+
+ Return:
+ None
+
+ Global variables
+ None
+
+ Specific functions:
+ None
+
+ Comments:
+
+ a.) Measured corrector value is read from a file "fic"
+ b.) correctors are at the same location of some sextupoles,
+ so their multipole errors are added to the thick sextupoles
+ which also functions as these correctors.
+
+ 10/2010 Written by Jianfeng Zhang
+**********************************************************************/
+void SetCorrQtErr_fam(char const *fic, const int Fnum, const double conv, const double r0,
+ const int n, const double Bn, const double An)
+{
+ int i = 0, N = 0, corr_index = 0;
+ double bnL = 0.0, anL = 0.0;
+ double brho = 0.0, conv_strength = 0.0;
+ double corr; /* skew quadrupole horizontal or vertical corrector error, read from a file*/
+ int corrlistThick[120]; /* index of associated sextupole*/
+
+ FILE *fi;
+
+
+
+ brho = globval.Energy/0.299792458; /* magnetic rigidity */
+
+ // assign the design order
+ if(Cell[Elem_GetPos(Fnum,1)].Elem.M->n_design == 2 )
+ N = 2; /* skew quadrupole*/
+ else
+ N = 1; /* correctors, they act like dipoles, so N =1, but in the lattice reading, their n_design = 0!!!!*/
+
+
+ /* Open file with multipole errors*/
+ if ((fi = fopen(fic,"r")) == NULL)
+ {
+ fprintf(stderr, "Error while opening file %s \n",fic);
+ exit_(1);
+ }
+
+
+ /* find index of sextupole associated with the corrector */
+ // solution 1: find by names
+ // solution 2: use a predfined list
+ // solution 3: smothing smart ???
+ for (i=0; i< GetnKid(Fnum); i++){
+ if (trace) fprintf(stdout, "%d\n", i);
+
+ corr_index = Elem_GetPos(Fnum, i+1);
+
+ if (Cell[corr_index-1].Elem.PName[0] == 's' && Cell[corr_index-1].Elem.PName[1] == 'x')
+ corrlistThick[i] = corr_index-1;
+ else{
+
+ if (Cell[corr_index+1].Elem.PName[0] == 's' && Cell[corr_index+1].Elem.PName[1] == 'x')
+ corrlistThick[i] = corr_index+1;
+ else{
+
+ if (Cell[corr_index+2].Elem.PName[0] == 's' && Cell[corr_index+2].Elem.PName[1] == 'x')
+ corrlistThick[i] = corr_index+2;
+ else{
+
+ if (Cell[corr_index-2].Elem.PName[0] == 's' && Cell[corr_index-2].Elem.PName[1] == 'x')
+ corrlistThick[i] = corr_index-2;
+ else{
+
+ if (Cell[corr_index+3].Elem.PName[0] == 's' && Cell[corr_index+3].Elem.PName[1] == 'x')
+ corrlistThick[i] = corr_index+3;
+ else{
+
+ if (Cell[corr_index-3].Elem.PName[0] == 's' && Cell[corr_index-3].Elem.PName[1] == 'x')
+ corrlistThick[i] = corr_index-3;
+ else fprintf(stdout, "Warning: Sextupole not found associated with corrector or skew quadrupole! \n");
+ }
+ }
+ }
+ }
+ }
+ }
+
+
+ // add the multipole errors to the associated sextupole
+ for (i = 0; i < GetnKid(Fnum); i++)
+ {
+ fscanf(fi,"%le \n", &corr); /* read the corrector values from a file */
+
+ if (r0 == 0.0) {
+ // input is: (b_n*L), (a_n*L) ???
+ add_bnL_sys_elem(Cell[corrlistThick[i]].Fnum,Cell[corrlistThick[i]].Knum, n, Bn, An);
+ } else {
+ conv_strength = corr*conv/brho;
+ // absolute integrated error field strength
+ bnL = Bn/pow(r0, n-N)*conv_strength;
+ anL = An/pow(r0, n-N)*conv_strength;
+
+ add_bnL_sys_elem(Cell[corrlistThick[i]].Fnum,Cell[corrlistThick[i]].Knum, n, bnL, anL);
+
+ }
+ }
+ fclose(fi); /* close corrector file */
+}
+
+