From 2a935dc2f88a167629abec0cbab8e63d21cb1c9c Mon Sep 17 00:00:00 2001
From: nadolski <nadolski@9a6e40ed-f3a0-4838-9b4a-bf418f78e88d>
Date: Sun, 23 Jun 2013 21:22:19 +0000
Subject: [PATCH] add new number for file
---
tracy/tracy/src/soleillib.cc | 781 ++++++++++++++++-------------------
1 file changed, 367 insertions(+), 414 deletions(-)
diff --git a/tracy/tracy/src/soleillib.cc b/tracy/tracy/src/soleillib.cc
index e14af69..411b086 100644
--- a/tracy/tracy/src/soleillib.cc
+++ b/tracy/tracy/src/soleillib.cc
@@ -1,6 +1,6 @@
/* Tracy-3
- J. Bengtsson, CBP, LBL 1990 - 1994 Pascal version
+ J. Bengtsson, CPB, LBL 1990 - 1994 Pascal version
SLS, PSI 1995 - 1997
M. Boege SLS, PSI 1998 C translation
L. Nadolski SOLEIL 2002 - Link to NAFF, Radia field maps
@@ -399,12 +399,9 @@ void SetErr2(long seed,double fac)
setrancut(normcut=2L);
iniranf(seed);
- for (i = 1L; i <= globval.Cell_nLoc; i++)
- {
- if (Cell[i].Elem.Pkind == Mpole)
- {
- if (Cell[i].Elem.M->n_design == 2L && Cell[i].dT[1] == 0) // exclude skew quads
- {
+ for (i = 1L; i <= globval.Cell_nLoc; i++){
+ if (Cell[i].Elem.Pkind == Mpole){
+ if (Cell[i].Elem.M->n_design == 2L && Cell[i].dT[1] == 0){ // exclude skew quads
if ((pair%2)==0) theta = fac*normranf(); /* random error every 2 elements (quad split into 2) */
pair++;
Cell[i].Elem.M->PBpar[HOMmax-2L] = -Cell[i].Elem.M->PBpar[HOMmax+2L]*sin(2.0*theta);
@@ -424,22 +421,22 @@ void SetErr2(long seed,double fac)
Purpose: read and set the definition of the vacuum chamber
between different sections around the ring from file
- AperFile.dat.
-
- In AperFile.dat,
- 1) line begin with "#" is comment line
- 2) first line Name1: Start
- first line Name2: All
- 3) the numbers of MK1 and MK2 should be the same in the lattice
- 4) MK1 is defined before MK2 in the lattice
- 5)
- 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
- dymax: maxmum y value of vacuum chamber
-
+ AperFile.dat.
+
+ In AperFile.dat,
+ 1) line begin with "#" is comment line
+ 2) first line Name1: Start
+ first line Name2: All
+ 3) the numbers of MK1 and MK2 should be the same in the lattice
+ 4) MK1 is defined before MK2 in the lattice
+ 5)
+ 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
+ dymax: maxmum y value of vacuum chamber
+
Input:
@@ -491,67 +488,67 @@ void ReadCh(const char *AperFile)
/* read the aperture setting */
{
sscanf(line,"%s %s %lf %lf %lf %lf",
- Name1,Name2, &dxmin, &dxmax, &dymin, &dymax);
+ Name1,Name2, &dxmin, &dxmax, &dymin, &dymax);
if (strcmp("Start", Name1)==0 && strcmp("All", Name2)==0) {
- if(prt)
- fprintf(stdout, "setting all apertures to \n"
- " dxmin = %e, dxmax = %e, dymin = %e, dymax = %e\n",
- dxmin, dxmax, dymin, dymax);
- set_aper_type(All, dxmin, dxmax, dymin, dymax);
- // ini_aper(dxmin, dxmax, dymin, dymax);
+ if(prt)
+ fprintf(stdout, "setting all apertures to \n"
+ " dxmin = %e, dxmax = %e, dymin = %e, dymax = %e\n",
+ dxmin, dxmax, dymin, dymax);
+ set_aper_type(All, dxmin, dxmax, dymin, dymax);
+ // ini_aper(dxmin, dxmax, dymin, dymax);
}
else {
/* read the vacuum chamber between section */
- Fnum1 = ElemIndex(Name1);
- Fnum2 = ElemIndex(Name2);
- if(Fnum1>0 && Fnum2>0) {
- /* if element Name1 is defined before element Name2, give error message*/
- if(Fnum1 > Fnum2){
- fprintf(stdout, "\nReadCh(): \n"
- " aperture file, Line %d, Element %s should be defined after Element %s \n",
- LineNum,Name1,Name2);
+ Fnum1 = ElemIndex(Name1);
+ Fnum2 = ElemIndex(Name2);
+ if(Fnum1>0 && Fnum2>0) {
+ /* if element Name1 is defined before element Name2, give error message*/
+ if(Fnum1 > Fnum2){
+ fprintf(stdout, "\nReadCh(): \n"
+ " aperture file, Line %d, Element %s should be defined after Element %s \n",
+ LineNum,Name1,Name2);
exit_(1);
}
- /* if the element is not unique in the lattice, give error message*/
- Kidnum1 = GetnKid(Fnum1);
- Kidnum2 = GetnKid(Fnum2);
- if(Kidnum1 != Kidnum2){
- fprintf(stdout, "\nReadCh(): \n"
- " vacuum aperture file, Line %d, the number of Element %s is not equal to"
- " the number of Element %s in lattice \n", LineNum, Name1, Name2);
+ /* if the element is not unique in the lattice, give error message*/
+ Kidnum1 = GetnKid(Fnum1);
+ Kidnum2 = GetnKid(Fnum2);
+ if(Kidnum1 != Kidnum2){
+ fprintf(stdout, "\nReadCh(): \n"
+ " vacuum aperture file, Line %d, the number of Element %s is not equal to"
+ " the number of Element %s in lattice \n", LineNum, Name1, Name2);
exit_(1);
- }
-
- if(prt)
- fprintf(stdout, "setting apertures to section:\n"
- " %s %s dxmin = %+e, dxmax = %+e, dymin = %+e, dymax = %+e\n",
- Name1, Name2, dxmin, dxmax, dymin, dymax);
-
-
- /* set the vacuum chamber*/
- //read the marker before the first element, and the markder after the last elment
- for(i=0;i<Kidnum1;i++){
- /* find the start and end index of the section*/
- k1 = Elem_GetPos(Fnum1, i+1);
- k2 = Elem_GetPos(Fnum2, i+1);
-
- for(j=1; j<globval.Cell_nLoc; j++){
- if(j>=k1 && j<k2){
- Cell[j].maxampl[X_][0] = dxmin;
- Cell[j].maxampl[X_][1] = dxmax;
- Cell[j].maxampl[Y_][0] = dymin;
- Cell[j].maxampl[Y_][1] = dymax;
- }
- }
- }
+ }
+
+ if(prt)
+ fprintf(stdout, "setting apertures to section:\n"
+ " %s %s dxmin = %+e, dxmax = %+e, dymin = %+e, dymax = %+e\n",
+ Name1, Name2, dxmin, dxmax, dymin, dymax);
+
+
+ /* set the vacuum chamber*/
+ //read the marker before the first element, and the markder after the last elment
+ for(i=0;i<Kidnum1;i++){
+ /* find the start and end index of the section*/
+ k1 = Elem_GetPos(Fnum1, i+1);
+ k2 = Elem_GetPos(Fnum2, i+1);
+
+ for(j=1; j<globval.Cell_nLoc; j++){
+ if(j>=k1 && j<k2){
+ Cell[j].maxampl[X_][0] = dxmin;
+ Cell[j].maxampl[X_][1] = dxmax;
+ Cell[j].maxampl[Y_][0] = dymin;
+ Cell[j].maxampl[Y_][1] = dymax;
+ }
+ }
+ }
}else {
- fprintf(stdout, "\nReadCh(): \n"
- " aperture file, Line %d, lattice does not contain section between element %s and element %s\n",
- LineNum,Name1, Name2);
- exit_(1);
- }
+ fprintf(stdout, "\nReadCh(): \n"
+ " aperture file, Line %d, lattice does not contain section between element %s and element %s\n",
+ LineNum,Name1, Name2);
+ exit_(1);
+ }
}
}
@@ -599,8 +596,8 @@ void ReadCh(const char *AperFile)
****************************************************************************/
void Trac_Tab(double x, double px, double y, double py, double dp,
- long nmax, long pos, long &lastn, long &lastpos, FILE *outf1,
- double Tx[][NTURN])
+ long nmax, long pos, long &lastn, long &lastpos, FILE *outf1,
+ double Tx[][NTURN])
{
bool lostF = true; /* Lost particle Flag */
Vector x1; /* tracking coordinates */
@@ -627,8 +624,8 @@ void Trac_Tab(double x, double px, double y, double py, double dp,
{
Cell_Pass(0,globval.Cell_nLoc, x1, lastpos);
if(trace) fprintf(outf1, "%6ld %+10.5e %+10.5e %+10.5e %+10.5e"
- " %+10.5e %+10.5e \n",
- lastn, x1[x_], x1[px_], x1[y_], x1[py_], x1[delta_], x1[ct_]);
+ " %+10.5e %+10.5e \n",
+ lastn, x1[x_], x1[px_], x1[y_], x1[py_], x1[delta_], x1[ct_]);
i = (lastn)-1;
Tx[x_][i] = x1[x_]; Tx[px_][i] = x1[px_];
Tx[y_][i] = x1[y_]; Tx[py_][i] = x1[py_];
@@ -638,8 +635,8 @@ void Trac_Tab(double x, double px, double y, double py, double dp,
else {
fprintf(stdout, "Trac_Tab: Particle lost \n");
fprintf(stdout, "%6ld %+10.5g %+10.5g %+10.5g"
- " %+10.5g %+10.5g %+10.5g \n",
- lastn, x1[x_], x1[px_], x1[y_], x1[py_], x1[delta_], x1[ct_]);
+ " %+10.5g %+10.5g %+10.5g \n",
+ lastn, x1[x_], x1[px_], x1[y_], x1[py_], x1[delta_], x1[ct_]);
lostF = false;
}
}
@@ -872,7 +869,7 @@ double get_D(const double df_x, const double df_y)
****************************************************************************/
#define NTERM2 10
void fmap(const char *FmapFile, long Nbx, long Nby, long Nbtour, double xmax, double ymax,
- double energy, bool diffusion, bool printloss)
+ double energy, bool diffusion, bool printloss)
{
FILE * outf; //file with the tunes at the grid point
FILE *outf_loss; //file with the information of the lost particle
@@ -919,10 +916,10 @@ void fmap(const char *FmapFile, long Nbx, long Nby, long Nbtour, double xmax, do
fprintf(outf,"# TRACY III -- %s -- %s \n", FmapFile, asctime2(newtime));
fprintf(outf,"# Frequencymap nu = f(x,z) \n");
// fprintf(outf,"# x[mm] z[mm] fx fy"
-// " dfx dfy D=log_10(sqrt(df_x^2+df_y^2))\n");
+// " dfx dfy D=log_10(sqrt(df_x^2+df_y^2))\n");
//
- fprintf(outf,"# xi[m] zi[m] fx fy"
- " dfx dfy\n");
+ fprintf(outf,"# xi[m] zi[m] fx fy"
+ " dfx dfy\n");
//file with lost particle information
@@ -936,10 +933,10 @@ if(printloss){
fprintf(outf_loss,"# Information of the lost particle: "
" lostp = 0 (no particle lost /1(horizontal) /2(vertical) /3(longitudinal) plane \n"
" Starting values & Phase space at turn Nturn and position s\n");
- fprintf(outf_loss,"# xi[m] zi[m] Nturn Plane s[m] x[m]"
- " xp[rad] z[m] zp[rad] delta ctau\n");
+ fprintf(outf_loss,"# xi[m] zi[m] Nturn Plane s[m] x[m]"
+ " xp[rad] z[m] zp[rad] delta ctau\n");
}
-
+
if ((Nbx < 1) || (Nby < 1))
fprintf(stdout,"fmap: Error Nbx=%ld Nby=%ld\n",Nbx,Nby);
@@ -973,9 +970,9 @@ if(printloss){
Get_NAFF(NTERM2, Nbtour, Tab, fx, fy, nb_freq);
Get_freq(fx,fy,&nux1,&nuy1); // gets nux and nuy
if (diffusion) { // diffusion
- // shift data for second round NAFF
+ // shift data for second round NAFF
Get_Tabshift(Tab,Tab0,Nbtour,Nbtour);
- // gets frequency vectors
+ // gets frequency vectors
Get_NAFF(NTERM2, Nbtour, Tab0, fx2, fy2, nb_freq);
Get_freq(fx2,fy2,&nux2,&nuy2); // gets nux and nuy
}
@@ -987,22 +984,22 @@ if(printloss){
// printout value
if (!diffusion){
- fprintf(outf,"%+10.6e %+10.6e %+10.6e %+10.6e\n",
- x, y, nux1, nuy1);
+ fprintf(outf,"%+10.6e %+10.6e %+10.6e %+10.6e\n",
+ x, y, nux1, nuy1);
//fprintf(stdout,"%+14.6e %+14.6e %+14.6e %+14.6e\n",
- // x, y, nux1, nuy1);
+ // x, y, nux1, nuy1);
}
else {
dfx = nux1 - nux2; dfy = nuy1 - nuy2;
fprintf(outf,"%+10.6e %+10.6e %+10.6e %+10.6e %+10.6e %+10.6e\n",
- x, y, nux1, nuy1, dfx, dfy);
+ x, y, nux1, nuy1, dfx, dfy);
//fprintf(stdout,"%+14.6e %+14.6e %+14.6e %+14.6e %+14.6e %+14.6e\n",
- // x, y, nux1, nuy1, dfx, dfy);
+ // x, y, nux1, nuy1, dfx, dfy);
}
//print out the information of the lost particle
if(printloss)
- fprintf(outf_loss, "%+6.3e %+6.3e %8ld %3d %+12.6f %+10.6f %+10.6f %+10.6f %+10.6f %+10.6f %+10.6f \n",
+ fprintf(outf_loss, "%+10.3e %+10.3e %8ld %3d %+12.6f %+10.3e %+10.3e %+10.3e %+10.6e %+10.3e %+10.3e \n",
x, y, lastn,status.lossplane, Cell[lastpos].S, x1[0],x1[1],
x1[2], x1[3], x1[4], x1[5]);
}
@@ -1064,7 +1061,7 @@ if(printloss){
****************************************************************************/
#define NTERM2 10
void fmap_p(const char *FmapFile_p, long Nbx, long Nby, long Nbtour, double xmax,
- double ymax, double energy, bool diffusion, bool printloss, int numprocs, int myid)
+ double ymax, double energy, bool diffusion, bool printloss, int numprocs, int myid)
{
FILE *outf; // output file
FILE *outf_loss; //file with the information of the lost particle
@@ -1082,7 +1079,7 @@ void fmap_p(const char *FmapFile_p, long Nbx, long Nby, long Nbtour, double xmax
struct tm *newtime;
char FmapFile[max_str];
- sprintf(FmapFile,"%d",myid);
+ sprintf(FmapFile,"myid%05d_",myid);
strcat(FmapFile,FmapFile_p);
fprintf(stdout, "%s\n",FmapFile);
@@ -1101,7 +1098,7 @@ void fmap_p(const char *FmapFile_p, long Nbx, long Nby, long Nbtour, double xmax
time(&aclock); /* Get time in seconds */
newtime = localtime(&aclock); /* Convert time to struct */
- if (trace) fprintf(stdout, "Entering fmap ... results in %s\n\n",FmapFile);
+ if (!trace) fprintf(stdout, "Entering fmap_p: ... results in %s\n\n",FmapFile);
/* Opening file */
if ((outf = fopen(FmapFile, "w")) == NULL){
@@ -1125,10 +1122,10 @@ void fmap_p(const char *FmapFile_p, long Nbx, long Nby, long Nbtour, double xmax
if(printloss){
fprintf(outf_loss,"# TRACY III -- %s -- %s \n", FmapLossFile, asctime2(newtime));
fprintf(outf_loss,"# Information of the lost particle: "
- " lostp = 0 (no particle lost /1(horizontal) /2(vertical) /3(longitudinal) plane \n"
- " Starting values & Phase space at turn Nturn and position s\n");
- fprintf(outf_loss,"# xi[m] zi[m] Nturn Plane s[m] x[m]"
- " xp[rad] z[m] zp[rad] delta ctau\n");
+ " lostp = 0 (no particle lost /1(horizontal) /2(vertical) /3(longitudinal) plane \n"
+ " Starting values & Phase space at turn Nturn and position s\n");
+ fprintf(outf_loss,"# xi[m] zi[m] Nturn Plane s[m] x[m]"
+ " xp[rad] z[m] zp[rad] delta ctau\n");
}
}
@@ -1154,7 +1151,7 @@ void fmap_p(const char *FmapFile_p, long Nbx, long Nby, long Nbtour, double xmax
integer=((int)Nbx)/numprocs;
residue=((int)Nbx)-integer*numprocs;
- fprintf(stdout, "fmap_p: myid:%d, integer:%d, residue:%d, numprocs:%d, Nbx:%ld\n\n",myid,integer,residue,
+ if (trace) fprintf(stdout, "fmap_p: myid:%d, integer:%d, residue:%d, numprocs:%d, Nbx:%ld\n\n",myid,integer,residue,
numprocs,Nbx);
//split tracking region (x,z) for each process
@@ -1180,44 +1177,37 @@ void fmap_p(const char *FmapFile_p, long Nbx, long Nby, long Nbtour, double xmax
else
Trac_Simple4DCOD(x,xp,y,zp,energy,ctau,nturn,Tab,&status2);
- if (status2) // if trajectory is stable
- {
+ if (status2){ // if trajectory is stable
// gets frequency vectors
Get_NAFF(NTERM2, Nbtour, Tab, fx, fy, nb_freq);
Get_freq(fx,fy,&nux1,&nuy1); // gets nux and nuy
- if (diffusion)
- { // diffusion
- // shift data for second round NAFF
+ if (diffusion){ // diffusion
+ // shift data for second round NAFF
Get_Tabshift(Tab,Tab0,Nbtour,Nbtour);
- // gets frequency vectors
+ // gets frequency vectors
Get_NAFF(NTERM2, Nbtour, Tab0, fx2, fy2, nb_freq);
Get_freq(fx2,fy2,&nux2,&nuy2); // gets nux and nuy
}
} // unstable trajectory
- else
- { //zeroing output
+ else{ //zeroing output
nux1 = 0.0; nuy1 = 0.0;
nux2 = 0.0; nuy2 = 0.0;
}
// printout value
- if (!diffusion)
- {
+ if (!diffusion){
fprintf(outf,"%+10.6e %+10.6e %+10.6e %+10.6e\n",x, y, nux1, nuy1);
- //fprintf(stdout,"%+14.6e %+14.6e %+14.6e %+14.6e\n",x, y, nux1, nuy1);
}
- else
- {
+ else{
dfx = nux1 - nux2;
dfy = nuy1 - nuy2;
fprintf(outf,"%+10.6e %+10.6e %+10.6e %+10.6e %+10.6e %+10.6e\n", x, y, nux1, nuy1, dfx, dfy);
- //fprintf(stdout,"%+14.6e %+14.6e %14.6e %14.6e %+14.6e %+14.6e\n", x, y, nux1, nuy1, dfx, dfy);
}
//print out the information of the lost particle
if(printloss)
- fprintf(outf_loss, "%+6.3e %+6.3e %8ld %3d %+12.6f %+10.6f %+10.6f %+10.6f %+10.6f %+10.6f %+10.6f \n",
+ fprintf(outf_loss, "%+10.3e %+10.3e %8ld %3d %+12.6f %+10.3e %+10.3e %+10.3e %+10.6e %+10.3e %+10.3e \n",
x, y, lastn,status.lossplane, Cell[lastpos].S, x1[0],x1[1],
x1[2], x1[3], x1[4], x1[5]);
}
@@ -1345,7 +1335,7 @@ if(printloss){
" Starting values & Phase space at turn Nturn and position s\n");
fprintf(outf_loss,"# dpi xi[m] Nturn lostp s[m] x[m]"
- " xp[rad] z[m] zp[rad] delta ctau\n");
+ " xp[rad] z[m] zp[rad] delta ctau\n");
}
if ((Nbx <= 1) || (Nbe <= 1))
fprintf(stdout,"fmapdp: Error Nbx=%ld Nbe=%ld\n",Nbx,Nbe);
@@ -1363,11 +1353,11 @@ if(printloss){
// If 6D Tracking diffusion turn off and x negative for dp negative
if ((globval.Cavity_on == true) && (dp < 0.0)){
x = x0 - sqrt((double)j)*xstep;
- diffusion = false;
+ diffusion = false;
}
else{
// x = x0 + sgn(j)*sqrt((double)abs(j))*xstep;
- x = x0 + sqrt((double)j)*xstep;
+ x = x0 + sqrt((double)j)*xstep;
}
if(printloss)
@@ -1391,7 +1381,7 @@ if(printloss)
// printout value
if (!diffusion){
- fprintf(outf,"%+10.6e %+10.6e %+10.6e %+10.6e\n", dp, x, nux1, nuy1);
+ fprintf(outf,"%+10.6e %+10.6e %+10.6e %+10.6e\n", dp, x, nux1, nuy1);
//fprintf(stdout,"% 10.6e % 10.6e % 10.6e % 10.6e\n", dp, x, nux1, nuy1);
}
else {
@@ -1470,7 +1460,7 @@ void fmapdp_p(const char *FmapdpFile_p, long Nbx, long Nbe, long Nbtour, double
double emax, double y, bool diffusion, bool printloss, int numprocs, int myid)
{
FILE * outf;
-FILE *outf_loss; //file with the information of the lost particle
+ FILE *outf_loss; //file with the information of the lost particle
long i = 0L, j = 0L;
double Tab[DIM][NTURN], Tab0[DIM][NTURN];
double fx[NTERM2], fy[NTERM2], fx2[NTERM2], fy2[NTERM2], dfx, dfy;
@@ -1484,10 +1474,10 @@ FILE *outf_loss; //file with the information of the lost particle
bool status2=true;
struct tm *newtime;
-char FmapdpFile[max_str];
- sprintf(FmapdpFile,"%d",myid);
+ char FmapdpFile[max_str];
+ sprintf(FmapdpFile,"myid%05d_",myid);
strcat(FmapdpFile,FmapdpFile_p);
- fprintf(stdout, "fmap_dp_p: %s\n",FmapdpFile);
+ if (trace) fprintf(stdout, "fmap_dp_p: %s\n",FmapdpFile);
//variables of the returned the tracked particle
long lastn = 0;
@@ -1520,22 +1510,20 @@ if(printloss){
}
}
- if(myid==0)
- {
+ if(myid==0){
fprintf(outf,"# TRACY III -- %s -- %s \n", FmapdpFile_p, asctime2(newtime));
fprintf(outf,"# nu = f(x) \n");
// fprintf(outf,"# dp[%%] x[mm] fx fy dfx dfy\n");
fprintf(outf,"# dp[m] x[m] fx fy dfx dfy\n");
-
if(printloss){
fprintf(outf_loss,"# TRACY III -- %s -- %s \n", FmapdpLossFile, asctime2(newtime));
fprintf(outf_loss,"# Information of the lost particle: "
- " lostp = 0 (no particle lost /1(horizontal) /2(vertical) /3(longitudinal) plane \n"
+ " lostp = 0 (no particle lost /1(horizontal) /2(vertical) /3(longitudinal) plane \n"
" Starting values & Phase space at turn Nturn and position s\n");
fprintf(outf_loss,"# dpi xi[m] Nturn lostp s[m] x[m]"
- " xp[rad] z[m] zp[rad] delta ctau\n");
+ " xp[rad] z[m] zp[rad] delta ctau\n");
}
}
@@ -1548,48 +1536,36 @@ if ((Nbx <= 1) || (Nbe <= 1))
xstep = xmax/sqrt((double)Nbx);
estep = 2.0*emax/Nbe;
-// for (i = 0; i < Nbe; i++)
-// Eace core or process calculate different region of fmapdp according to id number. MSChiu 2011/10/13
+// Each core or process calculate different region of fmapdp according to id number. MSChiu 2011/10/13
int deb,fin;
int integer,residue;
integer=((int)Nbe)/numprocs;
residue=((int)Nbe)-integer*numprocs;
- fprintf(stdout, "fmapdp_p: myid:%d, integer:%d, resideu:%d, numprocs:%d, Nbe:%ld\n\n",myid,integer,residue,numprocs,Nbe);
+ if (trace) fprintf(stdout, "fmapdp_p: myid:%d, integer:%d, resideu:%d, numprocs:%d, Nbe:%ld\n\n",myid,integer,residue,numprocs,Nbe);
//split tracking region for each process
- if(myid<residue)
- {
+ if(myid<residue){
deb=myid*(integer+1);
fin=(myid+1)*(integer+1);
}
- else
- {
+ else{
deb=residue*(integer+1)+(myid-residue)*integer;
fin=residue*(integer+1)+(myid+1-residue)*integer;
}
//begin tracking and FFT, and get tunes for the particle starts from (x,p)
- for (i = deb; i < fin; i++)
- {
+ for (i = deb; i < fin; i++){
dp = -emax + i*estep;
- // for (i = 0; i <= Nbe; i++) {
- // dp = -emax + i*estep;
-// if (!matlab) fprintf(outf,"\n");
- //fprintf(stdout,"\n");
for (j = 0; j<= Nbx; j++) {
-// for (j = -Nbx; j<= Nbx; j++) {
-
// IF 6D Tracking diffusion turn off and x negative for dp negative
if ((globval.Cavity_on == true) && (dp < 0.0)){
- // x = x0 - sgn(j)*sqrt((double)abs(j))*xstep;
- x = x0 - sqrt((double)j)*xstep;
- diffusion = false;
+ x = x0 - sqrt((double)j)*xstep;
+ diffusion = false;
}
else{
- // x = x0 + sgn(j)*sqrt((double)abs(j))*xstep;
- x = x0 + sqrt((double)j)*xstep;
+ x = x0 + sqrt((double)j)*xstep;
}
if(printloss)
@@ -1613,21 +1589,18 @@ if ((Nbx <= 1) || (Nbe <= 1))
// printout value
if (!diffusion){
- fprintf(outf,"% 10.6e % 10.6e % 10.6e % 10.6e\n", dp, x, nux1, nuy1);
- //fprintf(stdout,"% 10.6e % 10.6e % 10.6e % 10.6e\n", dp, x, nux1, nuy1);
+ fprintf(outf,"% 10.6e % 10.6e % 10.6e % 10.6e\n", dp, x, nux1, nuy1);
+ //fprintf(stdout,"% 10.6e % 10.6e % 10.6e % 10.6e\n", dp, x, nux1, nuy1);
}
else {
dfx = nux2 - nux1; dfy = nuy2 - nuy1;
- fprintf(outf,"% 10.6e % 10.6e % 10.6e % 10.6e % 10.6e % 10.6e\n",
+ fprintf(outf,"% 10.6e % 10.6e % 10.6e % 10.6e % 10.6e % 10.6e\n",
dp, x, nux1, nuy2, dfx, dfy);
- //fprintf(stdout,"% 10.6e % 10.6e % 10.6e % 10.6e % 10.6e % 10.6e\n",
- // dp, x, nux1, nuy2, dfx, dfy);
}
if(printloss)
fprintf(outf_loss," %-+8.3f %-+8.3f %-8ld %2d %+12.3f %+10.6f %+10.6f %+10.6f %+10.6f %+10.6f %+10.6f \n",
dp, x, lastn, status.lossplane, Cell[lastpos].S, x1[0],x1[1], x1[2], x1[3], x1[4], x1[5]);
-
}
}
@@ -3472,14 +3445,14 @@ void SetSkewQuad(void)
/****************************************************************************/
/* void MomentumAcceptance(char *MomAccFile, long deb, long fin,
double ep_min, double ep_max, long nstepp, double em_min,
- double em_max, long nstepm, long nturn)
+ double em_max, long nstepm, long nturn)
Purpose:
Compute momemtum acceptance along the ring, track the particle with
- different energy, momentum acceptance is the energy when the particle
- is lost or the last energy if the particle is not lost.
-
+ different energy, momentum acceptance is the energy when the particle
+ is lost or the last energy if the particle is not lost.
+
Based on the version in tracy 2.
-
+
Input:
MomAccFile file to save calculated momentum compact factor
deb first element for momentum acceptance,"debut" is beginning in French
@@ -3520,23 +3493,23 @@ void SetSkewQuad(void)
delta_closed_orbite = dp(cavity)/2
21/10/03 add array for vertical initial conditions using tracking
removed choice of tracking: now this should be done outside
-
+
23/07/10 modify the call variable to the Cell_Pass( ): j-1L --> j (L3435, L3590)
- since the Cell_Pass( ) is tracking from element i0 to i1(tracy 3), and
- the Cell_Pass( ) is tracking from element i0+1L to i1(tracy 2).
- 17/04/11 add number of turn
+ since the Cell_Pass( ) is tracking from element i0 to i1(tracy 3), and
+ the Cell_Pass( ) is tracking from element i0+1L to i1(tracy 2).
+ 17/04/11 add number of turn
27/06/11 fix the bug of the index in the taby and tabpy when calling Trac( );
fix the bug in the vertical closed orbit when calling Trac( ).
19/07/11 add the interface to save calculated momentum compact factor in the
- user defined file.
- add interface for user to define the start vertical amplitude at the
- entrance lattice element which is used to find the 6D closed orbit.
- 22/06/13 Add extra files for information about lost particles
-
+ user defined file.
+ add interface for user to define the start vertical amplitude at the
+ entrance lattice element which is used to find the 6D closed orbit.
+ 22/06/13 Add extra files for information about lost particles
+
****************************************************************************/
void MomentumAcceptance(char *MomAccFile, long deb, long fin,
double ep_min, double ep_max, long nstepp, double em_min,
- double em_max, long nstepm, long nturn, double ymax)
+ double em_max, long nstepm, long nturn, double ymax)
{
double dP = 0.0, dp1 = 0.0, dp2 = 0.0;
long lastpos = 0L,lastn = 0L;
@@ -3787,7 +3760,7 @@ void MomentumAcceptance(char *MomAccFile, long deb, long fin,
set_vectorcod(codvector, dP);
// coordinates around closed orbit specially usefull for 6D
- x0[x_] = codvector[0][x_];
+ x0[x_] = codvector[0][x_];
x0[px_] = codvector[0][px_];
x0[y_] = codvector[0][y_] + ymax;
x0[py_] = codvector[0][py_];
@@ -3913,11 +3886,11 @@ void MomentumAcceptance(char *MomAccFile, long deb, long fin,
Purpose:
Parallel version of MomentumAcceptance( ).
Compute momemtum acceptance along the ring, track the particle with
- different energy, momentum acceptance is the energy when the particle
- is lost or the last energy if the particle is not lost.
-
+ different energy, momentum acceptance is the energy when the particle
+ is lost or the last energy if the particle is not lost.
+
Based on the version in tracy 2.
-
+
Input:
MomAccFile file to save calculated momentum compact factor
deb first element for momentum acceptance,"debut" is beginning in French
@@ -3955,11 +3928,11 @@ Purpose:
Comments:
19/07/11 Add feature to do parallel calculation of momentum compact factor.
- Merged with the version written by Mao-Sen Qiu at Taiwan light source.
+ Merged with the version written by Mao-Sen Qiu at Taiwan light source.
****************************************************************************/
void MomentumAcceptance_p(char *_MomAccFile, long deb, long fin, double ep_min,
double ep_max, long nstepp, double em_min, double em_max,
- long nstepm, long nturn, double ymax, int numprocs,int myid)
+ long nstepm, long nturn, double ymax, int numprocs,int myid)
{
double dP = 0.0, dp1 = 0.0, dp2 = 0.0;
long lastpos = 0L,lastn = 0L;
@@ -3987,20 +3960,19 @@ void MomentumAcceptance_p(char *_MomAccFile, long deb, long fin, double ep_min,
// File opening for writing
char PhaseFile[50];
PhaseFile[0]='\0';
- sprintf(PhaseFile,"%d",myid);
+ sprintf(PhaseFile,"myid%05d_",myid);
strcat(PhaseFile,"phase.out");
//printf("%s\n",PhaseFile);
outf1 = fopen(PhaseFile, "w");
- if(myid==0)
- {
+ if(myid==0) {
fprintf(outf1,"# TRACY III -- %s \n", asctime2(newtime));
fprintf(outf1,"# i x xp y zp dp ctau\n#\n");
}
char MomAccFile[50];
MomAccFile[0]='\0';
- sprintf(MomAccFile,"%d",myid);
+ sprintf(MomAccFile,"myid%05d_",myid);
strcat(MomAccFile,_MomAccFile);
fprintf(stdout, "%s\n",MomAccFile);
@@ -4016,7 +3988,7 @@ void MomentumAcceptance_p(char *_MomAccFile, long deb, long fin, double ep_min,
// pos = deb; // starting position or element index in the ring
/***************************************************************/
- fprintf(stdout,"Computing initial conditions ... \n");
+ fprintf(stdout,"MomentumAcceptance_p: Computing initial conditions ... \n");
/***************************************************************/
// cod search has to be done in 4D since in 6D it is zero
@@ -4029,7 +4001,7 @@ void MomentumAcceptance_p(char *_MomAccFile, long deb, long fin, double ep_min,
taby0 = (double **)malloc((nstepp)*sizeof(double*));
tabpy0 = (double **)malloc((nstepp)*sizeof(double*));
if (taby0 == NULL || tabpy0 == NULL){
- fprintf(stdout,"1 out of memory \n"); return;
+ fprintf(stdout,"MomentumAcceptance_p: 1 out of memory \n"); return;
}
for (i = 1L; i <= nstepp; i++)
@@ -4037,9 +4009,8 @@ void MomentumAcceptance_p(char *_MomAccFile, long deb, long fin, double ep_min,
// Dynamical allocation 0 to nstepp -1
taby0[i-1L] = (double *)malloc((fin+1L)*sizeof(double));
tabpy0[i-1L] = (double *)malloc((fin+1L)*sizeof(double));
- if (taby0[i-1L] == NULL || tabpy0[i-1L] == NULL)
- {
- fprintf(stdout,"2 out of memory \n");
+ if (taby0[i-1L] == NULL || tabpy0[i-1L] == NULL){
+ fprintf(stdout,"MomentumAcceptance_p: 2 out of memory \n");
return;
}
@@ -4065,33 +4036,27 @@ void MomentumAcceptance_p(char *_MomAccFile, long deb, long fin, double ep_min,
// Store vertical initial conditions
// case where deb is not element 1
- if (deb > 1L)
- {
+ if (deb > 1L){
Cell_Pass(1L, deb - 1L, x0, lastpos); // track from 1 to deb-1L element
j = deb -1L;
- if (lastpos != j)
- { // look if stable
+ if (lastpos != j){ // look if stable
taby0 [i- 1L][j] = 1.0;
tabpy0[i- 1L][j] = 1.0;
}
- else
- { // stable case
+ else{ // stable case
taby0 [i - 1L][j] = x0[2] - codvector[deb-1L][2];
tabpy0[i - 1L][j] = x0[3] - codvector[deb-1L][3];
}
}
- else
- { // case where deb is element 1
+ else { // case where deb is element 1
j = deb - 1L;
taby0 [i - 1L][j] = x0[2] - codvector[j][2];
tabpy0[i - 1L][j] = x0[3] - codvector[j][3];
}
- for (j = deb; j < fin; j++)
- { // loop over elements
+ for (j = deb; j < fin; j++){ // loop over elements
Cell_Pass(j, j, x0, lastpos);
- //Cell_Pass(j -1L, j, x0, lastpos);
if (lastpos != j){ // look if stable
taby0 [i - 1L][j] = 1.0;
@@ -4100,7 +4065,6 @@ void MomentumAcceptance_p(char *_MomAccFile, long deb, long fin, double ep_min,
else{ // stable case
taby0 [i - 1L][j] = x0[2] - codvector[j][2];
tabpy0[i - 1L][j] = x0[3] - codvector[j][3];
-// fprintf(stdout,"z0= % e py0= % e\n", taby0 [i - 1L][j], tabpy0 [i - 1L][j]);
}
}
}
@@ -4109,7 +4073,7 @@ void MomentumAcceptance_p(char *_MomAccFile, long deb, long fin, double ep_min,
globval.radiation = radiationflag;
/***************************************************************/
- fprintf(stdout,"Computing positive momentum acceptance ... \n");
+ fprintf(stdout,"MomentumAcceptance_p: Computing positive momentum acceptance ... \n");
/***************************************************************/
//split tracking element region for each process
@@ -4118,35 +4082,32 @@ void MomentumAcceptance_p(char *_MomAccFile, long deb, long fin, double ep_min,
int integer,residue;
//the end element should not less than start element
-if(fin < deb){
- fprintf(stdout, "MomentumAcceptance_p: End element index %ld should be NOT"
+ if(fin < deb){
+ fprintf(stdout, "MomentumAcceptance_p: End element index %ld should be NOT"
"smaller than the start element index %ld\n",fin,deb);
- exit_(1);
-}
+ exit_(1);
+ }
integer=(fin-deb+1)/numprocs;
residue=(fin-deb+1)-integer*numprocs;
- fprintf(stdout, "MomentumAcceptance_p: myid:%d, integer:%d, resideu:%d, numprocs:%d\n",
+ if (trace) fprintf(stdout, "MomentumAcceptance_p: myid:%d, integer:%d, resideu:%d, numprocs:%d\n",
myid,integer,residue,numprocs);
//split tracking element region for each process
//the start element is from deb
- if(myid<residue)
- {
+ if(myid<residue){
debN=myid*(integer+1) +deb;
finN=(myid+1)*(integer+1) +deb;
}
- else
- {
+ else{
debN=residue*(integer+1)+(myid-residue)*integer +deb;
finN=residue*(integer+1)+(myid+1-residue)*integer +deb;
}
//do
- for(pos=debN;pos<finN;pos++)
- {
+ for(pos=debN;pos<finN;pos++) {
getcod(dP=0.0, lastpos); // determine closed orbit
getelem(pos,&Cell);
@@ -4162,15 +4123,14 @@ if(fin < deb){
pos, x, px, y, py, delta, ctau0);
if (trace) fprintf(stdout,"%3ld %6.4g %6.4g %6.4g %6.4g %6.4g %6.4g\n", pos,
- globval.CODvect[0], globval.CODvect[1],globval.CODvect[2],
- globval.CODvect[3], globval.CODvect[4],globval.CODvect[5]);
+ globval.CODvect[x_], globval.CODvect[px_],globval.CODvect[y_],
+ globval.CODvect[py_], globval.CODvect[delta_],globval.CODvect[ct_]);
dp1 = 0.0;
dp2 = 0.0;
i = 0L;
- do // Tracking over nturn
- {
+ do {// Tracking over nturn
i++;
dp1 = dp2;
@@ -4179,7 +4139,6 @@ if(fin < deb){
else
dp2 = ep_max;
- // if (trace) fprintf(stdout, "i=%4ld pos=%4ld dp=%6.4g\n",i,pos,dp2);
if (trace)
fprintf(stdout, "i=%4ld dp=%6.4g pos=%3ld %6.4g %6.4g %6.4g %6.4g %6.4g %6.4g\n",
i, dp2, pos, x, px, y+taby0[i-1L][pos-1L],
@@ -4189,7 +4148,7 @@ if(fin < deb){
pos, taby0[i-1L][pos-1L], tabpy0[i-1L][pos-1L]);
Trac(x, px, y+taby0[i-1L][pos-1L], py+tabpy0[i-1L][pos-1L], dp2+delta ,
- ctau0, nturn, pos, lastn, lastpos, outf1);
+ ctau0, nturn, pos, lastn, lastpos, outf1, x1);
}while (((lastn) == nturn) && (i != nstepp));
@@ -4211,8 +4170,7 @@ if(fin < deb){
}//while(pos != fin);
// free memory
- for (i = 1L; i <= nstepp; i++)
- {
+ for (i = 1L; i <= nstepp; i++){
free(taby0 [i - 1L]);
free(tabpy0[i - 1L]);
}
@@ -4243,7 +4201,7 @@ if(fin < deb){
taby0 = (double **)malloc((nstepm)*sizeof(double*));
tabpy0 = (double **)malloc((nstepm)*sizeof(double*));
if (taby0 == NULL || tabpy0 == NULL){
- fprintf(stdout,"1 out of memory \n"); return;
+ fprintf(stdout,"MomentumAcceptance_p: 1 out of memory \n"); return;
}
for (i = 1L; i <= nstepm; i++){ // loop over energy
@@ -4251,7 +4209,7 @@ if(fin < deb){
taby0[i-1L] = (double *)malloc((fin+1L)*sizeof(double));
tabpy0[i-1L] = (double *)malloc((fin+1L)*sizeof(double));
if (taby0[i-1L] == NULL || tabpy0[i-1L] == NULL){
- fprintf(stdout,"2 out of memory \n"); return;
+ fprintf(stdout,"MomentumAcceptance_p: 2 out of memory \n"); return;
}
// compute dP
@@ -4290,7 +4248,6 @@ if(fin < deb){
j = deb - 1L;
taby0 [i - 1L][j] = x0[2] - codvector[j][2];
tabpy0[i - 1L][j] = x0[3] - codvector[j][3];
-// fprintf(stdout,"z0= % e py0= % e\n", taby0 [i - 1L][j], tabpy0 [i - 1L][j]);
}
for (j = deb; j < fin; j++){ // loop over elements
@@ -4303,8 +4260,6 @@ if(fin < deb){
else{ // stable case
taby0 [i - 1L][j] = x0[2] - codvector[j][2];
tabpy0[i - 1L][j] = x0[3] - codvector[j][3];
-// fprintf(stdout,"dP= % e pos= %ld z0= % e py0= % e\n",
-// dP, j, taby0 [i - 1L][j], tabpy0 [i - 1L][j]);
}
}
}
@@ -4313,30 +4268,28 @@ if(fin < deb){
globval.radiation = radiationflag;
/***************************************************************/
- fprintf(stdout,"Computing negative momentum acceptance ... \n");
+ fprintf(stdout,"MomentumAcceptance_p: Computing negative momentum acceptance ... \n");
/***************************************************************/
// do
- for(pos=debN;pos<finN;pos++)
-{
+ for(pos=debN;pos<finN;pos++){
getcod(dP=0.0, lastpos); // determine closed orbit
getelem(pos,&Cell);
// coordinates around closed orbit which is non zero for 6D tracking
- x = Cell.BeamPos[0];
- px = Cell.BeamPos[1];
- y = Cell.BeamPos[2];
- py = Cell.BeamPos[3];
- delta = Cell.BeamPos[4];
- ctau0 = Cell.BeamPos[5];
- fprintf(stdout,"%3ld %6.4g %6.4g %6.4g %6.4g %6.4g %6.4g\n",
+ x = Cell.BeamPos[x_];
+ px = Cell.BeamPos[px_];
+ y = Cell.BeamPos[y_];
+ py = Cell.BeamPos[py_];
+ delta = Cell.BeamPos[delta_];
+ ctau0 = Cell.BeamPos[ct_];
+ if (trace) fprintf(stdout,"%3ld %6.4g %6.4g %6.4g %6.4g %6.4g %6.4g\n",
pos, x, px, y, py, delta, ctau0);
dp1 = 0.0;
dp2 = 0.0;
i = 0L;
- do // Tracking over nturn
- {
+ do {// Tracking over nturn
i++;
dp1 = dp2;
if (nstepm != 1L) {
@@ -4347,9 +4300,10 @@ if(fin < deb){
}
if (trace)
fprintf(stdout, "i=%4ld pos=%4ld dp=%6.4g %6.4g %6.4g %6.4g %6.4g %6.4g %6.4g\n",
- i,pos,dp2, x, px, y+taby0[i-1L][pos-1L], py+tabpy0[i-1L][pos-1L], dp2+delta, ctau0);
- Trac(x, px, y+taby0[i-1L][pos-1L], py+tabpy0[i-1L][pos-1L],
- dp2+delta , ctau0, nturn, pos, lastn, lastpos, outf1);
+ i,pos,dp2, x, px, y+taby0[i-1L][pos-1L], py+tabpy0[i-1L][pos-1L],
+ dp2+delta, ctau0);
+ Trac(x, px, y+taby0[i-1L][pos-1L], py+tabpy0[i-1L][pos-1L],
+ dp2+delta , ctau0, nturn, pos, lastn, lastpos, outf1, x1);
}
while (((lastn) == nturn) && (i != nstepm));
@@ -4357,10 +4311,10 @@ if(fin < deb){
getelem(lastpos,&Clost);
getelem(pos,&Cell);
- if (!trace) fprintf(stdout, "i=%4ld pos=%4ld dp=%6.4g\n",i,pos,dp2);
- if (!trace) fprintf(stdout, "pos=%4ld z0 =% 10.5f py0 =% 10.5f \n",
+ if (trace) fprintf(stdout, "i=%4ld pos=%4ld dp=%6.4g\n",i,pos,dp2);
+ if (trace) fprintf(stdout, "pos=%4ld z0 =% 10.5f py0 =% 10.5f \n",
pos, taby0[i-1L][pos-1L], tabpy0[i-1L][pos-1L]);
- if (!trace) fprintf(stdout, "%4ld %10.5f %10.5f %10.5f %*s\n",
+ if (trace) fprintf(stdout, "%4ld %10.5f %10.5f %10.5f %*s\n",
pos, Cell.S, dp1, Clost.S, 5, Clost.Elem.PName);
fprintf(outf2,"%4ld %10.3f %-5.5s %+10.2e %5d %10.3f %-5.5s "
"%+10.3e %+10.3e %+10.3e %+10.3e %+10.3e %+10.3e \n",
@@ -4372,11 +4326,11 @@ if(fin < deb){
//while(pos != fin);
// free memory
- for (i = 1L; i <= nstepp; i++)
- {
+ for (i = 1L; i <= nstepp; i++){
free(taby0 [i - 1L]);
free(tabpy0[i - 1L]);
}
+
free(taby0);
free(tabpy0);
@@ -4385,7 +4339,6 @@ if(fin < deb){
fclose(outf2);
}
-
/****************************************************************************/
/* set_vectorcod(double codvector[Cell_nLocMax][6], double dP)
@@ -4608,7 +4561,7 @@ void spectrum(long Nbx, long Nby, long Nbtour, double xmax, double ymax,
****************************************************************************/
void TracCO(double x, double px, double y, double py, double dp, double ctau,
- long nmax, long pos, long &lastn, long &lastpos, FILE *outf1)
+ long nmax, long pos, long &lastn, long &lastpos, FILE *outf1)
{
bool lostF; /* Lost particle Flag */
Vector x1; /* tracking coordinates */
@@ -4643,8 +4596,8 @@ void TracCO(double x, double px, double y, double py, double dp, double ctau,
(lastn)++;
if (!trace) { // print initial conditions
fprintf(outf1, "%6ld %+10.5e %+10.5e %+10.5e %+10.5e"
- " %+10.5e %+10.5e \n",
- lastn, x1[x_], x1[px_], x1[y_], x1[py_], x1[delta_], x1[ct_]);
+ " %+10.5e %+10.5e \n",
+ lastn, x1[x_], x1[px_], x1[y_], x1[py_], x1[delta_], x1[ct_]);
}
// Cell_Pass(pos-1L, globval.Cell_nLoc, x1, lastpos);
@@ -4657,8 +4610,8 @@ void TracCO(double x, double px, double y, double py, double dp, double ctau,
{
fprintf(stdout, "TracCO: Particle lost \n");
fprintf(stdout, "turn=%6ld %+10.5g %+10.5g %+10.5g"
- " %+10.5g %+10.5g %+10.5g \n",
- lastn, x1[x_], x1[px_], x1[y_], x1[py_], x1[delta_], x1[ct_]);
+ " %+10.5g %+10.5g %+10.5g \n",
+ lastn, x1[x_], x1[px_], x1[y_], x1[py_], x1[delta_], x1[ct_]);
}
}
@@ -5210,7 +5163,7 @@ void Coupling_Edwards_Teng(void)
tphi = 0.0, /* tan(phi) */
cphi = 0.0, /* cos(phi) */
sphi = 0.0; /* sin(phi) */
-
+
Matrix M, N, m, n, D, A, B, R, S;
Matrix Rinv, Dinv, nm, MminusN, tS, tn, U2T, dummy, T, U, Sigma;
Vector V;
@@ -5633,9 +5586,9 @@ void PhaseLongitudinalHamiltonien(void)
/* Laskar's integrator is not a good idea here, since the correction factor is
not integrable */
const double D1 = 0.675603595979829E0;
- const double D2 =-0.175603595979829E0;
- const double C2 = 0.135120719195966E1;
- const double C3 =-0.170241438391932E1;
+ const double D2 =-0.175603595979829E0;
+ const double C2 = 0.135120719195966E1;
+ const double C3 =-0.170241438391932E1;
FILE *outf;
const char fic[] = "longitudinal.out";
@@ -5891,7 +5844,7 @@ void Enveloppe2(double x, double px, double y, double py, double dp, double ntur
// }
fprintf(stdout, "xcod=%.5e mm zcod=% .5e mm \n",
- globval.CODvect[0]*1e3, globval.CODvect[2]*1e3);
+ globval.CODvect[0]*1e3, globval.CODvect[2]*1e3);
if ((outf = fopen(fic, "w")) == NULL) {
fprintf(stdout, "Enveloppe: error while opening file %s\n", fic);
@@ -6025,11 +5978,11 @@ void set_RFVoltage(const int Fnum, const double V_RF){
Read multipole errors from a file
The input format of the file is:
-
- seed radom_number ; this set is optional, and only works for the rms error
-
- keyWords sys/rms raduis when the field error is meausred "r0", field error order "n",
- field error component "Bn", field error component "An"; "n","Bn,""An",...
+
+ seed radom_number ; this set is optional, and only works for the rms error
+
+ keyWords sys/rms raduis when the field error is meausred "r0", field error order "n",
+ field error component "Bn", field error component "An"; "n","Bn,""An",...
Input:
@@ -6049,8 +6002,8 @@ void set_RFVoltage(const int Fnum, const double V_RF){
Comments:
10/2010 Written by Jianfeng Zhang
01/2011 Fix the bug for reading the end of line symbol "\n" , "\r",'\r\n'
- at different operation system
- 04/2011 Change the set of 'seed' for rms error in file, now it's mandatory.
+ at different operation system
+ 04/2011 Change the set of 'seed' for rms error in file, now it's mandatory.
*****************************************************************************************************/
void ReadFieldErr(const char *FieldErrorFile)
{
@@ -6086,60 +6039,60 @@ void ReadFieldErr(const char *FieldErrorFile)
/* check the line is whether comment line or null line*/
if (strstr(line, "#") == NULL && strcmp(line,"\n") != 0 &&
strcmp(line,"\r") != 0 &&strcmp(line,"\r\n") != 0) {
-
+
sscanf(line, "%s", name);
-
- if (strcmp("seed", name) == 0) { // the line to set random seed
+
+ if (strcmp("seed", name) == 0) { // the line to set random seed
sscanf(line, "%*s %d", &seed_val);
fprintf(stdout, "ReadFieldErr: setting random seed to %d\n", seed_val);
set_rnd = true;
- iniranf(seed_val);
+ iniranf(seed_val);
} else{//line to set (n Bn An sequence)
-
- /*read and assign the key words and measure radius*/
- sscanf(line, " %*s %s %lf",keywrd, &r0);
- if (prt) fprintf(stdout, "\nsetting <%s> multipole error to: %-5s r0 = %7.1le\n",keywrd,name,r0);
-
- rms = (strcmp("rms", keywrd) == 0)? true : false;
- if (rms && !set_rnd) {
+
+ /*read and assign the key words and measure radius*/
+ sscanf(line, " %*s %s %lf",keywrd, &r0);
+ if (prt) fprintf(stdout, "\nsetting <%s> multipole error to: %-5s r0 = %7.1le\n",keywrd,name,r0);
+
+ rms = (strcmp("rms", keywrd) == 0)? true : false;
+ if (rms && !set_rnd) {
fprintf(stdout, "ReadFieldErr: seed not defined\n");
exit(1);
}
-
- // skip first three parameters
- strtok(line, " \t");
- strtok(NULL, " \t");
- strtok(NULL, " \t");
-
- /* read the end of line symbol '\n','\r','\r\n' at different operation system*/
- while ((prm = strtok(NULL, " \t")) != NULL && strcmp(prm, "\n") != 0 &&
- strcmp(prm, "\r") != 0 && strcmp(prm, "\r\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)
- fprintf(stdout, " n = %2d, Bn = %9.1e, An = %9.1e\n", n, Bn, An);
-
-
- /* set multipole errors to horizontal correctors of soleil ring*/
- if(strcmp("hcorr", name) == 0)
- AddCorrQtErr_fam(fic_hcorr,globval.hcorr,_convHcorr,keywrd,r0,n,Bn,An);
- /* set multipole errors to vertical correctors of soleil ring*/
+
+ // skip first three parameters
+ strtok(line, " \t");
+ strtok(NULL, " \t");
+ strtok(NULL, " \t");
+
+ /* read the end of line symbol '\n','\r','\r\n' at different operation system*/
+ while ((prm = strtok(NULL, " \t")) != NULL && strcmp(prm, "\n") != 0 &&
+ strcmp(prm, "\r") != 0 && strcmp(prm, "\r\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)
+ fprintf(stdout, " n = %2d, Bn = %9.1e, An = %9.1e\n", n, Bn, An);
+
+
+ /* set multipole errors to horizontal correctors of soleil ring*/
+ if(strcmp("hcorr", name) == 0)
+ AddCorrQtErr_fam(fic_hcorr,globval.hcorr,_convHcorr,keywrd,r0,n,Bn,An);
+ /* set multipole errors to vertical correctors of soleil ring*/
else if(strcmp("vcorr", name) == 0)
- AddCorrQtErr_fam(fic_vcorr,globval.vcorr,_convVcorr,keywrd,r0,n,Bn,An);
- /* set multipole errors to skew quadrupoles of soleil ring*/
- else if(strcmp("qt", name) == 0)
- AddCorrQtErr_fam(fic_skew,globval.qt,_convQt,keywrd,r0,n,Bn,An);
- else
- /* set errors for other multipole*/
- AddFieldErrors(name,keywrd, r0, n, Bn, An) ;
- }
+ AddCorrQtErr_fam(fic_vcorr,globval.vcorr,_convVcorr,keywrd,r0,n,Bn,An);
+ /* set multipole errors to skew quadrupoles of soleil ring*/
+ else if(strcmp("qt", name) == 0)
+ AddCorrQtErr_fam(fic_skew,globval.qt,_convQt,keywrd,r0,n,Bn,An);
+ else
+ /* set errors for other multipole*/
+ AddFieldErrors(name,keywrd, r0, n, Bn, An) ;
+ }
}//end of read the (n Bn An) sequence
//end of the line
@@ -6153,8 +6106,8 @@ void ReadFieldErr(const char *FieldErrorFile)
/***********************************************************************
void AddFieldErrors(const char *name, const char *keywrd,const double r0,
- const int n, const double Bn, const double An)
-
+ 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
@@ -6164,7 +6117,7 @@ void AddFieldErrors(const char *name, const char *keywrd,const double r0,
name type name or element name
keyword "rms" or "sys"
"rms": random multipole error
- "sys": systematic multipole error
+ "sys": systematic multipole error
r0 radius at which error is measured, error field is relative
to the design field strength when r0 !=0
n order of the error
@@ -6188,7 +6141,7 @@ void AddFieldErrors(const char *name, const char *keywrd,const double r0,
10/2010 Written by Jianfeng Zhang
**********************************************************************/
void AddFieldErrors(const char *name,const char *keywrd, const double r0,
- const int n, const double Bn, const double An)
+ const int n, const double Bn, const double An)
{
int Fnum = 0;
@@ -6212,8 +6165,8 @@ void AddFieldErrors(const char *name,const char *keywrd, const double r0,
/***********************************************************************
void SetFieldValues_type(const int N, const char *keywrd, const double r0,
- const int n, const double Bn, const double An)
-
+ 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.
@@ -6221,14 +6174,14 @@ void SetFieldValues_type(const int N, const char *keywrd, const double r0,
N type name
keywrd "rms" or "sys"
"rms": random multipole error
- "sys": systematic multipole error
+ "sys": systematic multipole error
r0 radius at which error is measured, error field is relative
to the design field strength when r0 != 0
- if r0 == 0, the Bn and An are absolute value.
+ if r0 == 0, the Bn and An are absolute value.
n order of the error
Bn relative B component of n-th error
An relative A component of n-th error
-
+
Output:
@@ -6251,7 +6204,7 @@ void SetFieldValues_type(const int N, const char *keywrd, const double r0,
short quadrupoles
**********************************************************************/
void AddFieldValues_type(const int N, const char *keywrd, const double r0,
- const int n, const double Bn, const double An)
+ const int n, const double Bn, const double An)
{
double bnL = 0.0, anL = 0.0, KLN = 0.0;
int k = 0;
@@ -6261,46 +6214,46 @@ void AddFieldValues_type(const int N, const char *keywrd, const double r0,
{
//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)
+ {
+ //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
- if (r0 == 0){
- bnL = Bn;
- anL = An;
- }else{
- bnL = Bn/pow(r0, n-N)*KLN;
+ else
+ KLN = GetKLpar(Cell[k].Fnum, Cell[k].Knum, N);/*other multipoles*/
+
+
+ //absolute integrated multipole error strength
+ if (r0 == 0){
+ bnL = Bn;
+ anL = An;
+ }else{
+ 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
- // for the lattice with quadrupoles which are cut into two halves
+
+ //NOT add the multipole errors of short quadrupole to long quadrupole qp2 & qp7 of soleil ring
+ // for the lattice with quadrupoles which are cut into two halves
if(N == 2 && strncmp(Cell[k].Elem.PName,"qp2",3)==0)
- Add_bnL_sys_elem(Cell[k].Fnum, Cell[k].Knum,keywrd, 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, keywrd, n, 0, 0);
- // for the lattice with full quadrupoles
- else if(N == 2 && strncmp(Cell[k].Elem.PName,"q2",2)==0)
- Add_bnL_sys_elem(Cell[k].Fnum, Cell[k].Knum, keywrd, n, 0, 0);
- else if(N == 2 && strncmp(Cell[k].Elem.PName,"q7",2)==0)
- Add_bnL_sys_elem(Cell[k].Fnum, Cell[k].Knum, keywrd, n, 0, 0);
- else
- //add errors to multipoles except qp2, qp7
- Add_bnL_sys_elem(Cell[k].Fnum, Cell[k].Knum, keywrd,n, bnL,anL);
+ Add_bnL_sys_elem(Cell[k].Fnum, Cell[k].Knum,keywrd, 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, keywrd, n, 0, 0);
+ // for the lattice with full quadrupoles
+ else if(N == 2 && strncmp(Cell[k].Elem.PName,"q2",2)==0)
+ Add_bnL_sys_elem(Cell[k].Fnum, Cell[k].Knum, keywrd, n, 0, 0);
+ else if(N == 2 && strncmp(Cell[k].Elem.PName,"q7",2)==0)
+ Add_bnL_sys_elem(Cell[k].Fnum, Cell[k].Knum, keywrd, n, 0, 0);
+ else
+ //add errors to multipoles except qp2, qp7
+ Add_bnL_sys_elem(Cell[k].Fnum, Cell[k].Knum, keywrd,n, bnL,anL);
}
}
}
/***********************************************************************
void AddFieldValues_fam(const int Fnum, const char *keywrd, const double r0,
- const int n, const double Bn, const double An)
-
+ 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.
@@ -6309,14 +6262,14 @@ void AddFieldValues_fam(const int Fnum, const char *keywrd, const double r0,
Fnum family name
keywrd "rms" or "sys"
"rms": random multipole error
- "sys": systematic multipole error
+ "sys": systematic multipole error
r0 radius at which error is measured
for the case of r0 ???????
-
+
n order of the error
Bn relative B component for the n-th error
An relative A component for the n-th error
-
+
Output:
@@ -6335,7 +6288,7 @@ void AddFieldValues_fam(const int Fnum, const char *keywrd, const double r0,
10/2010 Written by Jianfeng Zhang
**********************************************************************/
void AddFieldValues_fam(const int Fnum, const char *keywrd, const double r0,
- const int n, const double Bn, const double An)
+ 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;
@@ -6345,20 +6298,20 @@ void AddFieldValues_fam(const int Fnum, const char *keywrd, const double r0,
// find the integrated design field strength for multipole
- if (Cell[loc].Elem.M->n_design == 1)
+ 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*/
+ KLN = GetKLpar(Cell[loc].Fnum, Cell[loc].Knum, N);/* other multipole*/
- /* absolute integrated field strength*/
- if (r0 == 0){ //?????????
- bnL = Bn;
- anL = An;
- }else{
- bnL = Bn/pow(r0, n-N)*KLN;
+ /* absolute integrated field strength*/
+ if (r0 == 0){ //?????????
+ bnL = Bn;
+ anL = An;
+ }else{
+ bnL = Bn/pow(r0, n-N)*KLN;
anL = An/pow(r0, n-N)*KLN;
- }
- //add absolute multipole field error for the family
+ }
+ //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,keywrd, n, bnL, anL);
@@ -6368,23 +6321,23 @@ void AddFieldValues_fam(const int Fnum, const char *keywrd, const double r0,
/***********************************************************************
void add_bnL_sys_elem(const int Fnum, const int Knum, const char *keywrd,
- const int n, const double bnL, const double anL)
-
+ 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.
+ the summmation value replace the previous value,
+ in the PBsys definition of multipole.
Input:
Fnum family index
Knum kids index
keywrd "rms" or "sys"
"rms": random multipole error
- "sys": systematic multipole error
+ "sys": systematic multipole error
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:
@@ -6405,7 +6358,7 @@ void add_bnL_sys_elem(const int Fnum, const int Knum, const char *keywrd,
Rms error on the quadrupoles: only works for full quadrupole, not for half quadrupole
**********************************************************************/
void Add_bnL_sys_elem(const int Fnum, const int Knum, const char *keywrd,
- const int n, const double bnL, const double anL)
+ const int n, const double bnL, const double anL)
{
elemtype elem;
double *elemMPB; //skew components of the multipole
@@ -6445,19 +6398,19 @@ void Add_bnL_sys_elem(const int Fnum, const int Knum, const char *keywrd,
if (prt)
fprintf(stdout, "add the %s error: n=%d component of %s, bnL = %e, %e, anL = %e, %e\n",
- keywrd,n, Cell[Elem_GetPos(Fnum, Knum)].Elem.PName,
- bnL, elemMPB[HOMmax+n],anL, elemMPB[HOMmax-n]);
+ keywrd,n, Cell[Elem_GetPos(Fnum, Knum)].Elem.PName,
+ bnL, elemMPB[HOMmax+n],anL, elemMPB[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)
-
+ 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
@@ -6465,7 +6418,7 @@ void SetCorrQtErr_fam(char const *fic, const int Fnum, const double conv, const
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
+ An integrated A component for the n-th error
Output:
None
@@ -6483,13 +6436,13 @@ void SetCorrQtErr_fam(char const *fic, const int Fnum, const double conv, const
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.
-
+ so their multipole errors are added to the thick sextupoles
+ which also functions as these correctors.
+
10/2010 Written by Jianfeng Zhang
**********************************************************************/
void AddCorrQtErr_fam(char const *fic, const int Fnum, const double conv, const char *keywrd, const double r0,
- const int n, const double Bn, const double An)
+ const int n, const double Bn, const double An)
{
int i = 0, N = 0, corr_index = 0;
double bnL = 0.0, anL = 0.0;
@@ -6568,12 +6521,12 @@ void AddCorrQtErr_fam(char const *fic, const int Fnum, const double conv, const
Add_bnL_sys_elem(Cell[corrlistThick[i]].Fnum,Cell[corrlistThick[i]].Knum,keywrd, n, Bn, An);
} else {
conv_strength = corr*conv/brho;
- // absolute integrated error field strength
+ // 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,keywrd, n, bnL, anL);
-
+
+ Add_bnL_sys_elem(Cell[corrlistThick[i]].Fnum,Cell[corrlistThick[i]].Knum,keywrd, n, bnL, anL);
+
}
}
fclose(fi); /* close corrector file */
@@ -6647,8 +6600,8 @@ void FitTune4(long qf1,long qf2, long qd1, long qd2, double nux, double nuy)
/**********************************************************************
void PrintTrack(const char *TrackFile, double x, double px, double y,double py,
double delta, double ctau, long int nmax)
-
- Purpose:
+
+ Purpose:
Print the coordinates at each lattice element by tracking around COD
Input:
@@ -6666,7 +6619,7 @@ void PrintTrack(const char *TrackFile, double x, double px, double y,double py,
Comments:
Written by Jianfeng Zhang @ synchro. soleil 05/2011
-**********************************************************************/
+**********************************************************************/
void PrintTrack(const char *TrackFile, double x, double px, double y,double py,
double delta, double ctau, long int nmax)
{
@@ -6685,7 +6638,7 @@ void PrintTrack(const char *TrackFile, double x, double px, double y,double py,
fprintf(outf, "# Tracking using TRACY III-- %s -- %s\n",TrackFile, asctime2(newtime));
fprintf(outf, "#\n");
- // fprintf(outf, "# work tunes: %7.5f %7.5f\n",globval.TotalTune[0], globval.TotalTune[1]);
+ // fprintf(outf, "# work tunes: %7.5f %7.5f\n",globval.TotalTune[0], globval.TotalTune[1]);
fprintf(outf, "# i ElemName S x p_x y p_y");
fprintf(outf, " delta cdt NTurns \n");
fprintf(outf, "# [m] [mm] [mrad] [mm] [mrad]");
@@ -6698,10 +6651,10 @@ void PrintTrack(const char *TrackFile, double x, double px, double y,double py,
x0[y_] = y;
x0[py_] = py;
x0[delta_] = delta;
- x0[ct_] = ctau;
+ x0[ct_] = ctau;
//get the close orbit at each element around the ring
- set_vectorcod(codvector, delta);
- //get the absolute initial coordinates
+ set_vectorcod(codvector, delta);
+ //get the absolute initial coordinates
x2[x_] = x0[x_] + codvector[1][x_];
x2[px_] = x0[px_] + codvector[1][px_];
x2[y_] = x0[y_] + codvector[1][y_];
@@ -6717,41 +6670,41 @@ void PrintTrack(const char *TrackFile, double x, double px, double y,double py,
//print the coordinates at each elements
do {
pos = 1;//track from first element
- (lastn)++;
+ (lastn)++;
for (i = 0; i < nv_; i++) //nv_ = 6
x1[i] = x2[i];
while( pos <= globval.Cell_nLoc){
- Cell_Pass(pos, pos, x2, lastpos);
- //check whether particle is lost
- if (!CheckAmpl(x2, pos)){
- fprintf(stderr,"Error!!! %ld turn, Particle lost at element: %s!",
- lastn, Cell[pos].Elem.PName);
- exit_(1);
- }
- //get the coordinates around the closed orbit
+ Cell_Pass(pos, pos, x2, lastpos);
+ //check whether particle is lost
+ if (!CheckAmpl(x2, pos)){
+ fprintf(stderr,"Error!!! %ld turn, Particle lost at element: %s!",
+ lastn, Cell[pos].Elem.PName);
+ exit_(1);
+ }
+ //get the coordinates around the closed orbit
for (i = x_; i <= py_; i++) //x_=0,px_=1,y_=2,py_=3
xf[i] = x2[i] - codvector[pos][i];
-
+
for (i = delta_; i <= ct_; i++) //delta_=4,ct_=5
if (globval.Cavity_on && (i != ct_))
xf[i] = x2[i] - codvector[pos][i];
else
xf[i] = x2[i];
- if (prt) {
- fprintf(stdout, "%4ld %4ld %s %6.4f %10.4f %10.4f %10.4f %10.4f %10.4f %10.4f %4ld \n",
+ if (prt) {
+ fprintf(stdout, "%4ld %4ld %s %6.4f %10.4f %10.4f %10.4f %10.4f %10.4f %10.4f %4ld \n",
pos, lastpos,Cell[pos].Elem.PName,Cell[pos].S, 1e3*xf[x_], 1e3*xf[px_],
- 1e3*xf[y_], 1e3*xf[py_],1e2*xf[delta_], 1e3*xf[ct_],lastn);
- }
+ 1e3*xf[y_], 1e3*xf[py_],1e2*xf[delta_], 1e3*xf[ct_],lastn);
+ }
fprintf(outf,"%6ld %s %8.4e %12.4e %12.4e %12.4e %12.4e %12.4e %12.4e %4ld \n",
pos,Cell[pos].Elem.PName,Cell[pos].S, 1e3*xf[x_], 1e3*xf[px_],
- 1e3*xf[y_], 1e3*xf[py_],1e2*xf[delta_], 1e3*xf[ct_],lastn);
-
- pos++;
- }//finish of tracking and printing to file
-
+ 1e3*xf[y_], 1e3*xf[py_],1e2*xf[delta_], 1e3*xf[ct_],lastn);
+
+ pos++;
+ }//finish of tracking and printing to file
+
} while ((lastn != nmax) && (lastpos == globval.Cell_nLoc));
// if (globval.MatMeth)
@@ -6777,14 +6730,14 @@ void PrintTrack(const char *TrackFile, double x, double px, double y,double py,
*****************************************************************************************************/
void ReadVirtualSkewQuad(const char *VirtualSkewQuadFile)
{
- int nqtcorr= 0; /* Number of virtual skew quadrupoles */
+ int nqtcorr= 0; /* Number of virtual skew quadrupoles */
int i=0;
long mOrder = 0L;
- double qtcorr[max_str]; /* virtual skew quad value in Amps */
+ double qtcorr[max_str]; /* virtual skew quad value in Amps */
double conv = 0.0, brho = 0.0, corr_strength =0.0;
- brho = globval.Energy/0.299792458; /* magnetic rigidity */
- conv = 93.88e-4; /*conversion des A en T*/
+ brho = globval.Energy/0.299792458; /* magnetic rigidity */
+ conv = 93.88e-4; /*conversion des A en T*/
nqtcorr = GetnKid(ElemIndex("SQ"));
fprintf(stdout, "Number of virtual skew quadrupoles: %d\n",nqtcorr);
@@ -6890,7 +6843,7 @@ void CODCorrect(const char *hcorr_file, const char *vcorr_file,int n_orbit,int n
/* cod = CorrectCOD_N(ae_file, n_orbit, n_scale, k);*/
fprintf(stdout, "\n");
-
+
if (cod){
/* fprintf(stdout, "done with orbit correction, now do coupling",
" correction plus vert. disp\n");*/
@@ -6900,7 +6853,7 @@ void CODCorrect(const char *hcorr_file, const char *vcorr_file,int n_orbit,int n
fprintf(stdout, "!!! Orbit correction failed\n");
exit_(1);
}
- //chk_cod(cod, "iter # %3d error_and_correction");
+ //chk_cod(cod, "iter # %3d error_and_correction");
// for debugging
//print flat lattice
--
GitLab