From 4f24ef71551651e4c3965540e61871c644f3b5c4 Mon Sep 17 00:00:00 2001
From: nadolski <nadolski@9a6e40ed-f3a0-4838-9b4a-bf418f78e88d>
Date: Sat, 22 Jun 2013 22:05:32 +0000
Subject: [PATCH] Format z->y Momentumcompaction with Loss information Printf
-> fprintf
---
tracy/tracy/inc/naffutils.h | 3 +-
tracy/tracy/inc/physlib.h | 5 +-
tracy/tracy/src/naffutils.cc | 6 +-
tracy/tracy/src/physlib.cc | 122 +++-
tracy/tracy/src/soleillib.cc | 1243 ++++++++++++++++++----------------
5 files changed, 767 insertions(+), 612 deletions(-)
diff --git a/tracy/tracy/inc/naffutils.h b/tracy/tracy/inc/naffutils.h
index 25f0e35..8ebf75f 100644
--- a/tracy/tracy/inc/naffutils.h
+++ b/tracy/tracy/inc/naffutils.h
@@ -27,7 +27,8 @@ void Trac_Simple4DCOD(double x, double px, double y, double py, double dp,
double ctau, long nmax, double Tx[][NTURN], bool *status);
void Trac_Simple4DCOD(double x, double px, double y, double py, double dp,
- double ctau, long nmax, double Tx[][NTURN], long& lastn, long& lastpos,ss_vect<double>& x1,bool *status2);
+ double ctau, long nmax, double Tx[][NTURN], long& lastn, long& lastpos,
+ ss_vect<double>& x1,bool *status2);
void Trac_Simple6DCOD(double x, double px, double y, double py, double dp,
double ctau, long nmax, double Tx[][NTURN], bool *status);
diff --git a/tracy/tracy/inc/physlib.h b/tracy/tracy/inc/physlib.h
index 49d45fb..f561164 100644
--- a/tracy/tracy/inc/physlib.h
+++ b/tracy/tracy/inc/physlib.h
@@ -119,7 +119,7 @@ void dynap(FILE *fp, double r, const double delta,
double get_aper(int n, double x[], double y[]);
void GetTrack(const char *file_name,
- long *n, double *x, double *px, double *y, double *py);
+ long *n, double *x, double *px, double *y, double *py);
void Getj(long n, double *x, double *px, double *y, double *py);
@@ -264,6 +264,9 @@ void GetChromTrac(long Nb, long Nbtour, double emax, double *xix, double *xiz);
void GetTuneTrac(long Nbtour, double emax, double *nux, double *nuz);
void Trac(double x, double px, double y, double py, double dp, double ctau,
long nmax, long pos, long &lastn, long &lastpos, FILE *outf1);
+void Trac(double x, double px, double y, double py, double dp, double ctau,
+ long nmax, long pos, long &lastn, long &lastpos, FILE *outf1,
+ ss_vect<double>& x1);
/* close orbit */
// simple precision
diff --git a/tracy/tracy/src/naffutils.cc b/tracy/tracy/src/naffutils.cc
index dfdac96..34833c0 100644
--- a/tracy/tracy/src/naffutils.cc
+++ b/tracy/tracy/src/naffutils.cc
@@ -8,9 +8,9 @@
*/
/*
- Current revision $Revision: 1.4 $
+ Current revision $Revision: 1.5 $
On branch $Name: not supported by cvs2svn $
- Latest change $Date: 2013-06-22 15:17:21 $ by $Author: nadolski $
+ Latest change $Date: 2013-06-22 22:05:32 $ by $Author: nadolski $
*/
#include "complexeheader.h"
@@ -227,7 +227,7 @@ void Trac_Simple4DCOD(double x, double px, double y, double py, double dp,
fprintf(stdout, "Trac_Simple4DCOD (2): Particle lost \n");
fprintf(stdout, "turn=%5ld plane= %1d %+10.5g %+10.5g %+10.5g %+10.5g %+10.5g %+10.5g \n", lastn-1,
status.lossplane, x1[0], x1[1], x1[2], x1[3], x1[4], x1[5]);
- }
+ }
}
}
diff --git a/tracy/tracy/src/physlib.cc b/tracy/tracy/src/physlib.cc
index c2a1bf0..cc55d28 100644
--- a/tracy/tracy/src/physlib.cc
+++ b/tracy/tracy/src/physlib.cc
@@ -6,7 +6,7 @@
L. Nadolski SOLEIL 2002 Link to NAFF, Radia field maps
J. Bengtsson NSLS-II, BNL 2004 -
*/
-/* Current revision $Revision: 1.21 $
+/* Current revision $Revision: 1.22 $
On branch $Name: not supported by cvs2svn $
Latest change by $Author: nadolski $
*/
@@ -890,6 +890,7 @@ void track_(double r, struct LOC_getdynap *LINK) {
****************************************************************************/
void Trac(double x, double px, double y, double py, double dp, double ctau,
long nmax, long pos, long &lastn, long &lastpos, FILE *outf1) {
+
bool lostF; /* Lost particle Flag */
Vector x1; /* tracking coordinates */
Vector2 aperture;
@@ -899,12 +900,111 @@ void Trac(double x, double px, double y, double py, double dp, double ctau,
aperture[0] = 1e0;
aperture[1] = 1e0;
- x1[0] = x;
- x1[1] = px;
- x1[2] = y;
- x1[3] = py;
- x1[4] = dp;
- x1[5] = ctau;
+ x1[x_] = x; x1[px_] = px;
+ x1[y_] = y; x1[py_] = py;
+ x1[delta_] = dp; x1[ct_] = ctau;
+
+ lastn = 0L;
+ (lastpos) = pos;
+ lostF = true;
+
+ if (trace)
+ fprintf(outf1, "\n");
+ 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_]);
+
+ // Cell_Pass(pos -1L, globval.Cell_nLoc, x1, lastpos);
+ Cell_Pass(pos, globval.Cell_nLoc, x1, lastpos);
+
+ if (lastpos == globval.Cell_nLoc)
+ do {
+ (lastn)++; /* 3) continue tracking for nmax turns*/
+ Cell_Pass(0L, pos - 1L, x1, lastpos); /* 1) check whether particle is stable at element from 0 to pos-1L*/
+
+ 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_]);
+
+ if (lastpos == pos - 1L)
+ Cell_Pass(pos, globval.Cell_nLoc, x1, lastpos); /* 2) check particle is stable at element from pos to end*/
+ // Cell_Pass(pos-1L,globval.Cell_nLoc, x1, lastpos);
+
+ } while (((lastn) < nmax) && ((lastpos) == globval.Cell_nLoc));
+
+ if (lastpos == globval.Cell_nLoc)
+ Cell_Pass(0L, pos, x1, lastpos);
+
+ if (trace && lastpos != pos) {
+ fprintf(stdout, "Trac: Particle lost \n");
+ fprintf(stdout, "turn:%6ld plane: %1d"
+ " %+10.5g %+10.5g %+10.5g %+10.5g %+10.5g %+10.5g \n", lastn,
+ status.lossplane, x1[x_], x1[px_], x1[y_], x1[py_], x1[delta_], x1[ct_]);
+ }
+}
+
+/****************************************************************************/
+/* void Trac(double x, double px, double y, double py, double dp, double ctau,
+ long nmax, long pos, long &lastn, long &lastpos, FILE *outf1)
+
+ Purpose:
+ Single particle tracking w/ respect to the chamber centrum
+ Based on the version of tracy 2 in soleillib.c
+ Input:
+ x, px, y, py 4 transverses coordinates
+ ctau c*tau
+ dp energy offset
+ nmax number of turns
+ pos starting position for tracking
+ aperture global physical aperture
+
+ Output:
+ lastn last n (should be nmax if not lost)
+ lastpos last position in the ring
+ outf1 output file with 6D phase at different pos
+ x1 tracking coordinates at the end
+
+ Return:
+ lastn: last tracking turn that particle is not lost
+ lastpos: last element position that particle is not lost
+
+ Global variables:
+ globval
+
+ specific functions:
+ Cell_Pass
+
+ Comments:
+ Absolute TRACKING with respect to the center of the vacuum vessel
+ BUG: last printout is wrong because not at pos but at the end of
+ the ring
+ 26/04/03 print output for phase space is for position pos now
+ 01/12/03 tracking from 0 to pos -1L instead of 0 to pos
+ (wrong observation point)
+
+ 23/07/10 change the call variable in Cell_Pass( ): pos-1L --> pos (L704, L717).
+ 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).
+
+ ****************************************************************************/
+void Trac(double x, double px, double y, double py, double dp, double ctau,
+ long nmax, long pos, long &lastn, long &lastpos, FILE *outf1, ss_vect<double> &x1) {
+
+ bool lostF; /* Lost particle Flag */
+ //Vector x1; /* tracking coordinates */
+ Vector2 aperture;
+
+ /* Compute closed orbit: useful if insertion devices */
+
+ aperture[0] = 1e0;
+ aperture[1] = 1e0;
+
+ x1[x_] = x; x1[px_] = px;
+ x1[y_] = y; x1[py_] = py;
+ x1[delta_] = dp; x1[ct_] = ctau;
lastn = 0L;
(lastpos) = pos;
@@ -915,7 +1015,7 @@ void Trac(double x, double px, double y, double py, double dp, double ctau,
if (trace)
fprintf(outf1,
"%6ld %+10.5e %+10.5e %+10.5e %+10.5e %+10.5e %+10.5e \n",
- lastn, x1[0], x1[1], x1[2], x1[3], x1[4], x1[5]);
+ lastn, x1[x_], x1[px_], x1[y_], x1[py_], x1[delta_], x1[ct_]);
// Cell_Pass(pos -1L, globval.Cell_nLoc, x1, lastpos);
Cell_Pass(pos, globval.Cell_nLoc, x1, lastpos);
@@ -929,7 +1029,7 @@ void Trac(double x, double px, double y, double py, double dp, double ctau,
fprintf(
outf1,
"%6ld %+10.5e %+10.5e %+10.5e %+10.5e %+10.5e %+10.5e \n",
- lastn, x1[0], x1[1], x1[2], x1[3], x1[4], x1[5]);
+ lastn, x1[x_], x1[px_], x1[y_], x1[py_], x1[delta_], x1[ct_]);
if (lastpos == pos - 1L)
Cell_Pass(pos, globval.Cell_nLoc, x1, lastpos); /* 2) check particle is stable at element from pos to end*/
@@ -940,11 +1040,11 @@ void Trac(double x, double px, double y, double py, double dp, double ctau,
if (lastpos == globval.Cell_nLoc)
Cell_Pass(0L, pos, x1, lastpos);
- if (lastpos != pos) {
+ if (trace && lastpos != pos) {
fprintf(stdout, "Trac: Particle lost \n");
fprintf(stdout, "turn:%6ld plane: %1d"
" %+10.5g %+10.5g %+10.5g %+10.5g %+10.5g %+10.5g \n", lastn,
- status.lossplane, x1[0], x1[1], x1[2], x1[3], x1[4], x1[5]);
+ status.lossplane, x1[x_], x1[px_], x1[y_], x1[py_], x1[delta_], x1[ct_]);
}
}
diff --git a/tracy/tracy/src/soleillib.cc b/tracy/tracy/src/soleillib.cc
index 0f9fa2f..1ddd89b 100644
--- a/tracy/tracy/src/soleillib.cc
+++ b/tracy/tracy/src/soleillib.cc
@@ -3,7 +3,7 @@
J. Bengtsson, CBP, 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
+ L. Nadolski SOLEIL 2002 - Link to NAFF, Radia field maps
J. Bengtsson NSLS-II, BNL 2004 -
J. Zhang SOLEIL 2010 ADD SOLEIL PARTS IN TRACY 2.7
*/
@@ -125,9 +125,9 @@ void InducedAmplitude(long spos)
for (k = 0; k <= imax ; k++) {
dP = -dpmax + 2*dpmax*k/imax;
/* Coordonnees initiales */
- x1[0] = 0.0; x1[1] = 0.0;
- x1[2] = 0.0; x1[3] = 0.0;
- x1[4] = dP ; x1[5] = 0.0;
+ x1[x_] = 0.0; x1[px_] = 0.0;
+ x1[y_] = 0.0; x1[py_] = 0.0;
+ x1[delta_] = dP ; x1[ct_] = 0.0;
/* Computes closed orbit and store it in a vector */
set_vectorcod(codvector, dP) ;
@@ -255,7 +255,7 @@ void Hcofonction(long pos, double dP)
Vector2 H;
getcod(dP, lastpos); /* determine closed orbit */
- if (lastpos != globval.Cell_nLoc) printf("Ring unstable for dp=%+e @ pos=%ld\n", dP, lastpos);
+ if (lastpos != globval.Cell_nLoc) fprintf(stdout, "Ring unstable for dp=%+e @ pos=%ld\n", dP, lastpos);
Ring_GetTwiss(true, dP); /* Compute and get Twiss parameters */
getelem(pos, &Cell); /* Position of the element pos */
@@ -318,9 +318,9 @@ void SetErr(long seed,double fac)
if(!prt){
- printf("\n");
- printf(" Setting random rotation errors to quadrupole magnets:\n");
- printf(" random seed number is: %ld, rms value of the rotation error is: %lf rad\n",seed,fac);
+ fprintf(stdout, "\n");
+ fprintf(stdout, " Setting random rotation errors to quadrupole magnets:\n");
+ fprintf(stdout, " random seed number is: %ld, rms value of the rotation error is: %lf rad\n",seed,fac);
}
setrancut(normcut=2L);
@@ -335,7 +335,7 @@ void SetErr(long seed,double fac)
theta = fac*normranf(); /* random error every 2 elements (quad split into 2) */
Cell[i].Elem.M->PBpar[HOMmax-2L] = -Cell[i].Elem.M->PBpar[HOMmax+2L]*sin(2.0*theta);
Cell[i].Elem.M->PBpar[HOMmax+2L] = Cell[i].Elem.M->PBpar[HOMmax+2L]*cos(2.0*theta);
- if (trace) printf("%6s % .5e % .5e % .5e\n",Cell[i].Elem.PName,
+ if (trace) fprintf(stdout, "%6s % .5e % .5e % .5e\n",Cell[i].Elem.PName,
Cell[i].Elem.M->PBpar[HOMmax-2L], Cell[i].Elem.M->PBpar[HOMmax+2L],theta);
Mpole_SetPB(Cell[i].Fnum, Cell[i].Knum, -2L);
@@ -391,9 +391,9 @@ void SetErr2(long seed,double fac)
if(!prt){
- printf("\n");
- printf(" Setting random rotation errors to quadrupole magnets:\n");
- printf(" random seed number is: %ld, rms value of the rotation error is: %lf rad\n",seed,fac);
+ fprintf(stdout, "\n");
+ fprintf(stdout, " Setting random rotation errors to quadrupole magnets:\n");
+ fprintf(stdout, " random seed number is: %ld, rms value of the rotation error is: %lf rad\n",seed,fac);
}
setrancut(normcut=2L);
@@ -409,7 +409,7 @@ void SetErr2(long seed,double fac)
pair++;
Cell[i].Elem.M->PBpar[HOMmax-2L] = -Cell[i].Elem.M->PBpar[HOMmax+2L]*sin(2.0*theta);
Cell[i].Elem.M->PBpar[HOMmax+2L] = Cell[i].Elem.M->PBpar[HOMmax+2L]*cos(2.0*theta);
- if (trace) printf("%6s % .5e % .5e % .5e\n",Cell[i].Elem.PName,
+ if (trace) fprintf(stdout, "%6s % .5e % .5e % .5e\n",Cell[i].Elem.PName,
Cell[i].Elem.M->PBpar[HOMmax-2L], Cell[i].Elem.M->PBpar[HOMmax+2L],theta);
Mpole_SetPB(Cell[i].Fnum, Cell[i].Knum, -2L);
@@ -473,8 +473,8 @@ void ReadCh(const char *AperFile)
fp = file_read(AperFile);
- printf("\n");
- printf("Loading and setting vacuum apertures to lattice elements...\n");
+ fprintf(stdout, "\n");
+ fprintf(stdout, "Loading and setting vacuum apertures to lattice elements...\n");
while ((line=fgets(in, max_str, fp))) {
/* kill preceding whitespace generated by "table" key
@@ -495,7 +495,7 @@ void ReadCh(const char *AperFile)
if (strcmp("Start", Name1)==0 && strcmp("All", Name2)==0) {
if(prt)
- printf("setting all apertures to \n"
+ 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);
@@ -509,7 +509,7 @@ void ReadCh(const char *AperFile)
if(Fnum1>0 && Fnum2>0) {
/* if element Name1 is defined before element Name2, give error message*/
if(Fnum1 > Fnum2){
- printf("\nReadCh(): \n"
+ fprintf(stdout, "\nReadCh(): \n"
" aperture file, Line %d, Element %s should be defined after Element %s \n",
LineNum,Name1,Name2);
exit_(1);
@@ -518,14 +518,14 @@ void ReadCh(const char *AperFile)
Kidnum1 = GetnKid(Fnum1);
Kidnum2 = GetnKid(Fnum2);
if(Kidnum1 != Kidnum2){
- printf("\nReadCh(): \n"
+ 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)
- printf("setting apertures to section:\n"
+ fprintf(stdout, "setting apertures to section:\n"
" %s %s dxmin = %+e, dxmax = %+e, dymin = %+e, dymax = %+e\n",
Name1, Name2, dxmin, dxmax, dymin, dymax);
@@ -547,7 +547,7 @@ void ReadCh(const char *AperFile)
}
}
}else {
- printf("\nReadCh(): \n"
+ 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);
@@ -556,7 +556,7 @@ void ReadCh(const char *AperFile)
}
// else /* print the comment line */
- // printf("%s", line);
+ // fprintf(stdout, "%s", line);
}
fclose(fp);
// turn on the global flag for CheckAmpl()
@@ -608,9 +608,9 @@ void Trac_Tab(double x, double px, double y, double py, double dp,
Vector2 aperture;
aperture[0] = 1e0; aperture[1] = 1e0;
- x1[0] = x; x1[1] = px;
- x1[2] = y; x1[3] = py;
- x1[4] = dp; x1[5] = 0.0;
+ x1[x_] = x; x1[px_] = px;
+ x1[y_] = y; x1[py_] = py;
+ x1[delta_] = dp; x1[ct_] = 0.0;
lastn = 0;
(lastpos)=pos;
@@ -628,18 +628,18 @@ 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[0], x1[1], x1[2], x1[3], x1[4], x1[5]);
+ lastn, x1[x_], x1[px_], x1[y_], x1[py_], x1[delta_], x1[ct_]);
i = (lastn)-1;
- Tx[0][i] = x1[0]; Tx[1][i] = x1[1];
- Tx[2][i] = x1[2]; Tx[3][i] = x1[3];
- Tx[4][i] = x1[4]; Tx[5][i] = x1[5];
+ Tx[x_][i] = x1[x_]; Tx[px_][i] = x1[px_];
+ Tx[y_][i] = x1[y_]; Tx[py_][i] = x1[py_];
+ Tx[delta_][i] = x1[delta_]; Tx[ct_][i] = x1[ct_];
}
else {
- printf("Trac_Tab: Particle lost \n");
+ 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[0], x1[1], x1[2], x1[3], x1[4], x1[5]);
+ lastn, x1[x_], x1[px_], x1[y_], x1[py_], x1[delta_], x1[ct_]);
lostF = false;
}
}
@@ -648,30 +648,30 @@ void Trac_Tab(double x, double px, double y, double py, double dp,
for (i = 1; i < nmax; i++) {
fprintf(outf1, "%6ld %+10.5e %+10.5e %+10.5e %+10.5e %+10.5e %+10.5e \n", i,
- Tx[0][i], Tx[1][i], Tx[2][i], Tx[3][i], Tx[4][i], Tx[5][i]);
+ Tx[x_][i], Tx[px_][i], Tx[y_][i], Tx[py_][i], Tx[delta_][i], Tx[ct_][i]);
}
}
/****************************************************************************/
-/* void TunesShiftWithAmplitude(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
+/* void TunesShiftWithAmplitude(long Nbx, long Nby, long Nbtour, double xmax, double ymax,
double energy)
Purpose:
- Compute nux, nuz with respect to x : nudx.out if xmax!=0
- with respect to z : nudz.out if zmax!=0
+ Compute nux, nuy with respect to x : nudx.out if xmax!=0
+ with respect to y : nudz.out if ymax!=0
for an energy offset delta=energy
over Nbtour turns of the ring
for x varying within [-xmax, xmax] around the closed orbit
- for z varying within [-zmax, zmax] around the closed orbit
+ for y varying within [-ymax, ymax] around the closed orbit
Input:
Nbx horizontal point number
- Nbz vertical point number
+ Nby vertical point number
Nbtour turn number
xmax maximum horizontal amplitude
- zmax maximum vertical amplitude
+ ymax maximum vertical amplitude
energy enrgy offset
Output:
@@ -693,15 +693,15 @@ void Trac_Tab(double x, double px, double y, double py, double dp,
****************************************************************************/
#define nterm 4
void TunesShiftWithAmplitude(const char *NudxFile, const char *NudzFile, long Nbx,
- long Nbz, long Nbtour, double xmax, double zmax, double energy)
+ long Nby, long Nbtour, double xmax, double ymax, double energy)
{
FILE * outf;
int i = 0;
- double Tab[6][NTURN], fx[nterm], fz[nterm];
- double x = 0.0 , xp = 0.0 , z = 0.0 , zp = 0.0;
+ double Tab[6][NTURN], fx[nterm], fy[nterm];
+ double x = 0.0 , xp = 0.0 , y = 0.0 , zp = 0.0;
double x0 = 1e-6, xp0= 0.0 , z0 = 1e-6, zp0 = 0.0;
double xstep = 0.0, zstep = 0.0;
- double nux = 0.0, nuz = 0.0;
+ double nux = 0.0, nuy = 0.0;
int nb_freq[2] = {0, 0};
bool stable = true;
struct tm *newtime;
@@ -709,7 +709,7 @@ void TunesShiftWithAmplitude(const char *NudxFile, const char *NudzFile, long Nb
/* Get time and date */
newtime = GetTime();
- if (!trace) printf("\n Entering TunesShiftWithAmplitude ... results in nudx.out\n\n");
+ if (!trace) fprintf(stdout, "\n Entering TunesShiftWithAmplitude ... results in nudx.out\n\n");
/////////////
// H tuneshift
@@ -725,10 +725,10 @@ void TunesShiftWithAmplitude(const char *NudxFile, const char *NudzFile, long Nb
fprintf(outf,"# Tracy III -- %s -- %s \n", NudxFile, asctime2(newtime));
fprintf(outf,"# nu = f(x) \n");
- fprintf(outf,"# x[m] z[m] fx fz \n");
+ fprintf(outf,"# x[m] z[m] fx fy \n");
- if ((Nbx <= 1) || (Nbz <= 1))
- fprintf(stdout,"TunesShiftWithAmplitude: Error Nbx=%ld Nbz=%ld\n",Nbx,Nbz);
+ if ((Nbx <= 1) || (Nby <= 1))
+ fprintf(stdout,"TunesShiftWithAmplitude: Error Nbx=%ld Nby=%ld\n",Nbx,Nby);
xstep = xmax/Nbx*2.0;
x0 = 1e-6 - xmax;
@@ -737,21 +737,21 @@ void TunesShiftWithAmplitude(const char *NudxFile, const char *NudzFile, long Nb
for (i = 0; i <= Nbx; i++) {
x = x0 + i*xstep ;
xp = xp0 ;
- z = z0 ;
+ y = z0 ;
zp = zp0 ;
- Trac_Simple4DCOD(x,xp,z,zp,energy,0.0,Nbtour,Tab,&stable); // tracking around closed orbit
+ Trac_Simple4DCOD(x,xp,y,zp,energy,0.0,Nbtour,Tab,&stable); // tracking around closed orbit
if (stable) {
- Get_NAFF(nterm, Nbtour, Tab, fx, fz, nb_freq); // gets frequency vectors
- Get_freq(fx,fz,&nux,&nuz); // gets nux and nuz
+ Get_NAFF(nterm, Nbtour, Tab, fx, fy, nb_freq); // gets frequency vectors
+ Get_freq(fx,fy,&nux,&nuy); // gets nux and nuy
}
else { // unstable
- nux = 0.0; nuz = 0.0;
+ nux = 0.0; nuy = 0.0;
}
fprintf(outf,"% 10.6e % 10.6e % 10.6e % 10.6e\n",
- x, z, nux, nuz);
+ x, y, nux, nuy);
}
fclose(outf);
}
@@ -760,7 +760,7 @@ void TunesShiftWithAmplitude(const char *NudxFile, const char *NudzFile, long Nb
// V tuneshift
/////////////
- if (fabs(zmax) > 0.0)
+ if (fabs(ymax) > 0.0)
{
/* Opening file */
if ((outf = fopen(NudzFile, "w")) == NULL) {
@@ -770,27 +770,27 @@ void TunesShiftWithAmplitude(const char *NudxFile, const char *NudzFile, long Nb
fprintf(outf,"# tracy III -- %s -- %s \n", NudzFile, asctime2(newtime));
fprintf(outf,"# nu = f(z) \n");
- fprintf(outf,"# x[mm] z[mm] fx fz \n");
+ fprintf(outf,"# x[mm] z[mm] fx fy \n");
- zstep = zmax/Nbz*2.0;
+ zstep = ymax/Nby*2.0;
x0 = 1e-3;
- z0 = 1e-6 - zmax;
- for (i = 0; i <= Nbz; i++) {
+ z0 = 1e-6 - ymax;
+ for (i = 0; i <= Nby; i++) {
x = x0 ;
xp = xp0;
- z = z0 + i*zstep;
+ y = z0 + i*zstep;
zp = zp0;
- Trac_Simple4DCOD(x,xp,z,zp,energy,0.0,Nbtour,Tab,&stable);
+ Trac_Simple4DCOD(x,xp,y,zp,energy,0.0,Nbtour,Tab,&stable);
if (stable) {
- Get_NAFF(nterm, Nbtour, Tab, fx, fz, nb_freq);
- Get_freq(fx,fz,&nux,&nuz); // gets nux and nuz
+ Get_NAFF(nterm, Nbtour, Tab, fx, fy, nb_freq);
+ Get_freq(fx,fy,&nux,&nuy); // gets nux and nuy
}
else {
- nux = 0.0; nuz =0.0;
+ nux = 0.0; nuy =0.0;
}
fprintf(outf,"% 10.6e % 10.6e % 10.6e % 10.6e\n",
- x, z, nux, nuz);
+ x, y, nux, nuy);
}
fclose(outf);
@@ -816,17 +816,17 @@ double get_D(const double df_x, const double df_y)
/****************************************************************************/
-/* void fmap(const char *FmapFile, long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
+/* void fmap(const char *FmapFile, long Nbx, long Nby, long Nbtour, double xmax, double ymax,
double energy, bool diffusion, bool loss)
Purpose:
- Compute a frequency map of Nbx x Nbz points
+ Compute a frequency map of Nbx x Nby points
For each set of initial conditions the particle is tracked over
Nbtour for an energy offset dp
- Frequency map is based on fixed beam energy, trace x versus z,
+ Frequency map is based on fixed beam energy, trace x versus y,
or, tracking transverse dynamic aperture for fixed momentum
(usually, on-momentum) particle.
@@ -839,7 +839,7 @@ double get_D(const double df_x, const double df_y)
Nbx horizontal step number
Nby vertical step number
xmax horizontal maximum amplitude
- zmax vertical maximum amplitude
+ ymax vertical maximum amplitude
Nbtour number of turn for tracking
energy particle energy offset
diffusion flag to calculate tune diffusion/ tune
@@ -866,24 +866,24 @@ double get_D(const double df_x, const double df_y)
frequency map analysis
11/06/2013 Modified by Jianfeng Zhang @ LAL
- The same as famp(onst char *FmapFile, long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
+ The same as famp(onst char *FmapFile, long Nbx, long Nby, long Nbtour, double xmax, double ymax,
double energy, bool diffusion); but with a flag to print/not print the final information of
the particle; if the final turn is not the "Nbtour", then the particle is lost.
****************************************************************************/
#define NTERM2 10
-void fmap(const char *FmapFile, long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
+void fmap(const char *FmapFile, long Nbx, long Nby, long Nbtour, double xmax, double ymax,
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
long i = 0L, j = 0L;
double Tab[DIM][NTURN], Tab0[DIM][NTURN];
- double fx[NTERM2], fz[NTERM2], fx2[NTERM2], fz2[NTERM2], dfx, dfz;
- double x = 0.0, xp = 0.0, z = 0.0, zp = 0.0;
+ double fx[NTERM2], fy[NTERM2], fx2[NTERM2], fy2[NTERM2], dfx, dfy;
+ double x = 0.0, xp = 0.0, y = 0.0, zp = 0.0;
double x0 = 1e-6, xp0 = 0.0, z0 = 1e-6, zp0 = 0.0;
const double ctau = 0.0;
double xstep = 0.0, zstep = 0.0;
- double nux1 = 0.0, nuz1 = 0.0, nux2 = 0.0, nuz2 = 0.0;
+ double nux1 = 0.0, nuy1 = 0.0, nux2 = 0.0, nuy2 = 0.0;
int nb_freq[2] = {0, 0};
long nturn = Nbtour;
bool status2 = true;
@@ -918,11 +918,11 @@ void fmap(const char *FmapFile, long Nbx, long Nbz, 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 fz"
-// " dfx dfz D=log_10(sqrt(df_x^2+df_y^2))\n");
+// fprintf(outf,"# x[mm] z[mm] fx fy"
+// " dfx dfy D=log_10(sqrt(df_x^2+df_y^2))\n");
//
- fprintf(outf,"# xi[m] zi[m] fx fz"
- " dfx dfz\n");
+ fprintf(outf,"# xi[m] zi[m] fx fy"
+ " dfx dfy\n");
//file with lost particle information
@@ -940,17 +940,17 @@ if(printloss){
" xp[rad] z[m] zp[rad] delta ctau\n");
}
- if ((Nbx < 1) || (Nbz < 1))
- fprintf(stdout,"fmap: Error Nbx=%ld Nbz=%ld\n",Nbx,Nbz);
+ if ((Nbx < 1) || (Nby < 1))
+ fprintf(stdout,"fmap: Error Nbx=%ld Nby=%ld\n",Nbx,Nby);
// steps in both planes
xstep = xmax/sqrt((double)Nbx);
- zstep = zmax/sqrt((double)Nbz);
+ zstep = ymax/sqrt((double)Nby);
// double number of turn if diffusion to compute
if (diffusion) nturn = 2*Nbtour;
- // px and pz zeroed
+ // px and py zeroed
xp = xp0;
zp = zp0;
@@ -958,52 +958,52 @@ if(printloss){
for (i = 0; i <= Nbx; i++) {
x = x0 + sqrt((double)i)*xstep;
//fprintf(stdout,"\n");
- for (j = 0; j<= Nbz; j++) {
- z = z0 + sqrt((double)j)*zstep;
+ for (j = 0; j<= Nby; j++) {
+ y = z0 + sqrt((double)j)*zstep;
//print out the lost information
if(printloss) // tracking around closed orbit
- Trac_Simple4DCOD(x,xp,z,zp,energy,ctau,nturn,Tab,lastn, lastpos, x1, &status2);
+ Trac_Simple4DCOD(x,xp,y,zp,energy,ctau,nturn,Tab,lastn, lastpos, x1, &status2);
else
// tracking around closed orbit
- Trac_Simple4DCOD(x,xp,z,zp,energy,ctau,nturn,Tab,&status2);
+ Trac_Simple4DCOD(x,xp,y,zp,energy,ctau,nturn,Tab,&status2);
if (status2) { // if trajectory is stable
// gets frequency vectors
- Get_NAFF(NTERM2, Nbtour, Tab, fx, fz, nb_freq);
- Get_freq(fx,fz,&nux1,&nuz1); // gets nux and nuz
+ 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
Get_Tabshift(Tab,Tab0,Nbtour,Nbtour);
// gets frequency vectors
- Get_NAFF(NTERM2, Nbtour, Tab0, fx2, fz2, nb_freq);
- Get_freq(fx2,fz2,&nux2,&nuz2); // gets nux and nuz
+ Get_NAFF(NTERM2, Nbtour, Tab0, fx2, fy2, nb_freq);
+ Get_freq(fx2,fy2,&nux2,&nuy2); // gets nux and nuy
}
} // unstable trajectory
else { //zeroing output
- nux1 = 0.0; nuz1 = 0.0;
- nux2 = 0.0; nuz2 = 0.0;
+ nux1 = 0.0; nuy1 = 0.0;
+ nux2 = 0.0; nuy2 = 0.0;
}
// printout value
if (!diffusion){
fprintf(outf,"%+10.6e %+10.6e %+10.6e %+10.6e\n",
- x, z, nux1, nuz1);
+ x, y, nux1, nuy1);
//fprintf(stdout,"%+14.6e %+14.6e %+14.6e %+14.6e\n",
- // x, z, nux1, nuz1);
+ // x, y, nux1, nuy1);
}
else {
- dfx = nux1 - nux2; dfz = nuz1 - nuz2;
+ dfx = nux1 - nux2; dfy = nuy1 - nuy2;
fprintf(outf,"%+10.6e %+10.6e %+10.6e %+10.6e %+10.6e %+10.6e\n",
- x, z, nux1, nuz1, dfx, dfz);
+ x, y, nux1, nuy1, dfx, dfy);
//fprintf(stdout,"%+14.6e %+14.6e %+14.6e %+14.6e %+14.6e %+14.6e\n",
- // x, z, nux1, nuz1, dfx, dfz);
+ // 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",
- x, z, lastn,status.lossplane, Cell[lastpos].S, x1[0],x1[1],
+ x, y, lastn,status.lossplane, Cell[lastpos].S, x1[0],x1[1],
x1[2], x1[3], x1[4], x1[5]);
}
}
@@ -1016,16 +1016,16 @@ if(printloss){
#undef NTERM2
/****************************************************************************/
-/* void fmap_p(const char *FmapFile_p, long Nbx, long Nbz, long Nbtour, double xmax,
- double zmax, double energy, bool diffusion, int numprocs, int myid)
+/* void fmap_p(const char *FmapFile_p, long Nbx, long Nby, long Nbtour, double xmax,
+ double ymax, double energy, bool diffusion, int numprocs, int myid)
Purpose:
Parallelized version of fmap( ).
- Compute a frequency map of Nbx x Nbz points
+ Compute a frequency map of Nbx x Nby points
For each set of initial conditions the particle is tracked over
Nbtour for an energy offset dp
- Frequency map is based on fixed beam energy, trace x versus z,
+ Frequency map is based on fixed beam energy, trace x versus y,
or, tracking transverse dynamic aperture for fixed momentum
(usually, on-momentum) particle.
@@ -1038,7 +1038,7 @@ if(printloss){
Nbx horizontal step number
Nby vertical step number
xmax horizontal maximum amplitude
- zmax vertical maximum amplitude
+ ymax vertical maximum amplitude
Nbtour number of turn for tracking
energy particle energy offset
matlab set file print format for matlab post-process; specific for nsrl-ii
@@ -1063,19 +1063,19 @@ if(printloss){
Merged with the version written by Mao-Sen Qiu at Taiwan light source.
****************************************************************************/
#define NTERM2 10
-void fmap_p(const char *FmapFile_p, long Nbx, long Nbz, long Nbtour, double xmax,
- double zmax, double energy, bool diffusion, bool printloss, int numprocs, int myid)
+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)
{
FILE *outf; // output file
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], fz[NTERM2], fx2[NTERM2], fz2[NTERM2], dfx, dfz;
- double x = 0.0, xp = 0.0, z = 0.0, zp = 0.0;
+ double fx[NTERM2], fy[NTERM2], fx2[NTERM2], fy2[NTERM2], dfx, dfy;
+ double x = 0.0, xp = 0.0, y = 0.0, zp = 0.0;
double x0 = 1e-6, xp0 = 0.0, z0 = 1e-6, zp0 = 0.0;
const double ctau = 0.0;
double xstep = 0.0, zstep = 0.0;
- double nux1 = 0.0, nuz1 = 0.0, nux2 = 0.0, nuz2 = 0.0;
+ double nux1 = 0.0, nuy1 = 0.0, nux2 = 0.0, nuy2 = 0.0;
int nb_freq[2] = {0, 0};
long nturn = Nbtour;
bool status2 = true;
@@ -1084,7 +1084,7 @@ void fmap_p(const char *FmapFile_p, long Nbx, long Nbz, long Nbtour, double xmax
char FmapFile[max_str];
sprintf(FmapFile,"%d",myid);
strcat(FmapFile,FmapFile_p);
- printf("%s\n",FmapFile);
+ fprintf(stdout, "%s\n",FmapFile);
//variables of the returned the tracked particle
long lastn = 0;
@@ -1101,7 +1101,7 @@ void fmap_p(const char *FmapFile_p, long Nbx, long Nbz, long Nbtour, double xmax
time(&aclock); /* Get time in seconds */
newtime = localtime(&aclock); /* Convert time to struct */
- if (trace) printf("Entering fmap ... results in %s\n\n",FmapFile);
+ if (trace) fprintf(stdout, "Entering fmap ... results in %s\n\n",FmapFile);
/* Opening file */
if ((outf = fopen(FmapFile, "w")) == NULL){
@@ -1120,7 +1120,7 @@ void fmap_p(const char *FmapFile_p, long Nbx, long Nbz, long Nbtour, double xmax
if(myid==0){
fprintf(outf,"# TRACY III -- %s -- %s \n", FmapFile_p, asctime2(newtime));
fprintf(outf,"# nu = f(x) \n");
- fprintf(outf,"# x[m] z[m] fx fz dfx dfz\n");
+ fprintf(outf,"# x[m] z[m] fx fy dfx dfy\n");
if(printloss){
fprintf(outf_loss,"# TRACY III -- %s -- %s \n", FmapLossFile, asctime2(newtime));
@@ -1132,17 +1132,17 @@ void fmap_p(const char *FmapFile_p, long Nbx, long Nbz, long Nbtour, double xmax
}
}
- if ((Nbx < 1) || (Nbz < 1))
- fprintf(stdout,"fmap_p: Error Nbx=%ld Nbz=%ld\n",Nbx,Nbz);
+ if ((Nbx < 1) || (Nby < 1))
+ fprintf(stdout,"fmap_p: Error Nbx=%ld Nby=%ld\n",Nbx,Nby);
// steps in both planes
xstep = xmax/sqrt((double)Nbx);
- zstep = zmax/sqrt((double)Nbz);
+ zstep = ymax/sqrt((double)Nby);
// double number of turn if diffusion to compute
if (diffusion) nturn = 2*Nbtour;
- // px and pz zeroed
+ // px and py zeroed
xp = xp0;
zp = zp0;
@@ -1172,53 +1172,53 @@ void fmap_p(const char *FmapFile_p, long Nbx, long Nbz, long Nbtour, double xmax
x = x0 + sqrt((double)i)*xstep;
// fprintf(stdout,"\n");
- for (j = 0; j< Nbz; j++) {
- z = z0 + sqrt((double)j)*zstep;
+ for (j = 0; j< Nby; j++) {
+ y = z0 + sqrt((double)j)*zstep;
// tracking around closed orbit
if(printloss)
- Trac_Simple4DCOD(x,xp,z,zp,energy,ctau,nturn,Tab,lastn, lastpos, x1, &status2);
+ Trac_Simple4DCOD(x,xp,y,zp,energy,ctau,nturn,Tab,lastn, lastpos, x1, &status2);
else
- Trac_Simple4DCOD(x,xp,z,zp,energy,ctau,nturn,Tab,&status2);
+ Trac_Simple4DCOD(x,xp,y,zp,energy,ctau,nturn,Tab,&status2);
if (status2) // if trajectory is stable
{
// gets frequency vectors
- Get_NAFF(NTERM2, Nbtour, Tab, fx, fz, nb_freq);
- Get_freq(fx,fz,&nux1,&nuz1); // gets nux and nuz
+ 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
Get_Tabshift(Tab,Tab0,Nbtour,Nbtour);
// gets frequency vectors
- Get_NAFF(NTERM2, Nbtour, Tab0, fx2, fz2, nb_freq);
- Get_freq(fx2,fz2,&nux2,&nuz2); // gets nux and nuz
+ Get_NAFF(NTERM2, Nbtour, Tab0, fx2, fy2, nb_freq);
+ Get_freq(fx2,fy2,&nux2,&nuy2); // gets nux and nuy
}
} // unstable trajectory
else
{ //zeroing output
- nux1 = 0.0; nuz1 = 0.0;
- nux2 = 0.0; nuz2 = 0.0;
+ nux1 = 0.0; nuy1 = 0.0;
+ nux2 = 0.0; nuy2 = 0.0;
}
// printout value
if (!diffusion)
{
- fprintf(outf,"%+10.6e %+10.6e %+10.6e %+10.6e\n",x, z, nux1, nuz1);
- //fprintf(stdout,"%+14.6e %+14.6e %+14.6e %+14.6e\n",x, z, nux1, nuz1);
+ 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
{
dfx = nux1 - nux2;
- dfz = nuz1 - nuz2;
+ dfy = nuy1 - nuy2;
- fprintf(outf,"%+10.6e %+10.6e %+10.6e %+10.6e %+10.6e %+10.6e\n", x, z, nux1, nuz1, dfx, dfz);
- //fprintf(stdout,"%+14.6e %+14.6e %14.6e %14.6e %+14.6e %+14.6e\n", x, z, nux1, nuz1, dfx, dfz);
+ 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",
- x, z, lastn,status.lossplane, Cell[lastpos].S, x1[0],x1[1],
+ x, y, lastn,status.lossplane, Cell[lastpos].S, x1[0],x1[1],
x1[2], x1[3], x1[4], x1[5]);
}
}
@@ -1235,14 +1235,14 @@ void fmap_p(const char *FmapFile_p, long Nbx, long Nbz, long Nbtour, double xmax
/****************************************************************************/
/* void fmapdp(const char *FmapdpFile, long Nbx, long Nbe, long Nbtour, double xmax, double emax,
- double z, bool diffusion)
+ double y, bool diffusion)
Purpose:
- Compute a frequency map of Nbx x Nbz points
+ Compute a frequency map of Nbx x Nby points
For each set of initial conditions the particle is tracked over
Nbtour for an energy offset dp
- Frequency map is based on fixed vertical amplitude z, trace x versus energy,
+ Frequency map is based on fixed vertical amplitude y, trace x versus energy,
or, tracking x for off-momentum particle.
The stepsize follows a square root law
@@ -1256,7 +1256,7 @@ void fmap_p(const char *FmapFile_p, long Nbx, long Nbz, long Nbtour, double xmax
Nbtour number of turns for tracking
xmax horizontal maximum amplitude
emax maximum energy
- z vertical amplitude
+ y vertical amplitude
diffusion flag to calculate tune diffusion/ tune
difference between the first Nbtour and last Nbtour
matlab set file print format for matlab post-process; specific for nsrl-ii
@@ -1285,17 +1285,17 @@ void fmap_p(const char *FmapFile_p, long Nbx, long Nbz, long Nbtour, double xmax
****************************************************************************/
#define NTERM2 10
void fmapdp(const char *FmapdpFile, long Nbx, long Nbe, long Nbtour, double xmax, double emax,
- double z, bool diffusion, bool printloss)
+ double y, bool diffusion, bool printloss)
{
FILE * outf;
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], fz[NTERM2], fx2[NTERM2], fz2[NTERM2], dfx, dfz;
+ double fx[NTERM2], fy[NTERM2], fx2[NTERM2], fy2[NTERM2], dfx, dfy;
double x = 0.0, xp = 0.0, zp = 0.0, dp = 0.0, ctau = 0.0;
double x0 = 1e-6, xp0 = 0.0, zp0 = 0.0;
double xstep = 0.0, estep = 0.0;
- double nux1 = 0.0, nuz1 = 0.0, nux2 = 0.0, nuz2 = 0.0;
+ double nux1 = 0.0, nuy1 = 0.0, nux2 = 0.0, nuy2 = 0.0;
int nb_freq[2] = {0, 0};
long nturn = Nbtour;
@@ -1328,8 +1328,8 @@ void fmapdp(const char *FmapdpFile, long Nbx, long Nbe, long Nbtour, double xmax
fprintf(outf,"# TRACY III -- %s -- %s \n", FmapdpFile, asctime2(newtime));
fprintf(outf,"# nu = f(x) \n");
-// fprintf(outf,"# dp[%%] x[mm] fx fz dfx dfz\n");
-fprintf(outf,"# dp[m] x[m] fx fz dfx dfz\n");
+// fprintf(outf,"# dp[%%] x[mm] fx fy dfx dfy\n");
+fprintf(outf,"# dp[m] x[m] fx fy dfx dfy\n");
//file with lost particle information
if(printloss){
@@ -1371,35 +1371,35 @@ if(printloss){
}
if(printloss)
- Trac_Simple4DCOD(x,xp,z,zp,dp,ctau,nturn,Tab,lastn, lastpos, x1, &status2);
+ Trac_Simple4DCOD(x,xp,y,zp,dp,ctau,nturn,Tab,lastn, lastpos, x1, &status2);
else
- Trac_Simple4DCOD(x,xp,z,zp,dp,ctau,nturn,Tab,&status2);
+ Trac_Simple4DCOD(x,xp,y,zp,dp,ctau,nturn,Tab,&status2);
if (status2) {
- Get_NAFF(NTERM2, Nbtour, Tab, fx, fz, nb_freq);
- Get_freq(fx,fz,&nux1,&nuz1); // gets nux and nuz
+ Get_NAFF(NTERM2, Nbtour, Tab, fx, fy, nb_freq);
+ Get_freq(fx,fy,&nux1,&nuy1); // gets nux and nuy
if (diffusion) { // diffusion
Get_Tabshift(Tab,Tab0,Nbtour,Nbtour); // shift data for second round NAFF
- Get_NAFF(NTERM2, Nbtour, Tab0, fx2, fz2, nb_freq); // gets frequency vectors
- Get_freq(fx2,fz2,&nux2,&nuz2); // gets nux and nuz
+ Get_NAFF(NTERM2, Nbtour, Tab0, fx2, fy2, nb_freq); // gets frequency vectors
+ Get_freq(fx2,fy2,&nux2,&nuy2); // gets nux and nuy
}
} // unstable trajectory
else { //zeroing output
- nux1 = 0.0; nuz1 = 0.0;
- nux2 = 0.0; nuz2 = 0.0;
+ nux1 = 0.0; nuy1 = 0.0;
+ nux2 = 0.0; nuy2 = 0.0;
}
// printout value
if (!diffusion){
- fprintf(outf,"%+10.6e %+10.6e %+10.6e %+10.6e\n", dp, x, nux1, nuz1);
- //fprintf(stdout,"% 10.6e % 10.6e % 10.6e % 10.6e\n", dp, x, nux1, nuz1);
+ 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; dfz = nuz2 - nuz1;
+ dfx = nux2 - nux1; dfy = nuy2 - nuy1;
fprintf(outf,"%+10.6e %+10.6e %+10.6e %+10.6e %+10.6e %+10.6e\n",
- dp, x, nux1, nuz2, dfx, dfz);
+ dp, x, nux1, nuy2, dfx, dfy);
//fprintf(stdout,"% 10.6e % 10.6e % 10.6e % 10.6e % 10.6e % 10.6e\n",
- // dp, x, nux1, nuz2, dfx, dfz);
+ // dp, x, nux1, nuy2, dfx, dfy);
}
if(printloss)
@@ -1417,16 +1417,16 @@ if(printloss)
/****************************************************************************/
/* void fmapdp_p(const char *FmapdpFile_p, long Nbx, long Nbe, long Nbtour,double xmax,
- double emax, double z, bool diffusion, int numprocs, int myid)
+ double emax, double y, bool diffusion, int numprocs, int myid)
Purpose:
Parallel version of fmapdp( ).
- Compute a frequency map of Nbx x Nbz points
+ Compute a frequency map of Nbx x Nby points
For each set of initial conditions the particle is tracked over
Nbtour for an energy offset dp
- Frequency map is based on fixed vertical amplitude z, trace x versus energy,
+ Frequency map is based on fixed vertical amplitude y, trace x versus energy,
or, tracking x for off-momentum particle.
The stepsize follows a square root law
@@ -1440,7 +1440,7 @@ if(printloss)
Nbtour number of turns for tracking
xmax horizontal maximum amplitude
emax maximum energy
- z vertical amplitude
+ y vertical amplitude
diffusion flag to calculate tune diffusion
numprocs Number of processes used to do parallel computing
myid process used to do parallel computing
@@ -1467,17 +1467,17 @@ if(printloss)
****************************************************************************/
#define NTERM2 10
void fmapdp_p(const char *FmapdpFile_p, long Nbx, long Nbe, long Nbtour, double xmax,
- double emax, double z, bool diffusion, bool printloss, int numprocs, int myid)
+ 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
long i = 0L, j = 0L;
double Tab[DIM][NTURN], Tab0[DIM][NTURN];
- double fx[NTERM2], fz[NTERM2], fx2[NTERM2], fz2[NTERM2], dfx, dfz;
+ double fx[NTERM2], fy[NTERM2], fx2[NTERM2], fy2[NTERM2], dfx, dfy;
double x = 0.0, xp = 0.0, zp = 0.0, dp = 0.0, ctau = 0.0;
double x0 = 1e-6, xp0 = 0.0, zp0 = 0.0;
double xstep = 0.0, estep = 0.0;
- double nux1 = 0.0, nuz1 = 0.0, nux2 = 0.0, nuz2 = 0.0;
+ double nux1 = 0.0, nuy1 = 0.0, nux2 = 0.0, nuy2 = 0.0;
int nb_freq[2] = {0, 0};
long nturn = Nbtour;
@@ -1524,8 +1524,8 @@ if(printloss){
{
fprintf(outf,"# TRACY III -- %s -- %s \n", FmapdpFile_p, asctime2(newtime));
fprintf(outf,"# nu = f(x) \n");
- // fprintf(outf,"# dp[%%] x[mm] fx fz dfx dfz\n");
- fprintf(outf,"# dp[m] x[m] fx fz dfx dfz\n");
+ // fprintf(outf,"# dp[%%] x[mm] fx fy dfx dfy\n");
+ fprintf(outf,"# dp[m] x[m] fx fy dfx dfy\n");
if(printloss){
@@ -1593,35 +1593,35 @@ if ((Nbx <= 1) || (Nbe <= 1))
}
if(printloss)
- Trac_Simple4DCOD(x,xp,z,zp,dp,ctau,nturn,Tab,lastn, lastpos, x1, &status2);
+ Trac_Simple4DCOD(x,xp,y,zp,dp,ctau,nturn,Tab,lastn, lastpos, x1, &status2);
else
- Trac_Simple4DCOD(x,xp,z,zp,dp,ctau,nturn,Tab,&status2);
+ Trac_Simple4DCOD(x,xp,y,zp,dp,ctau,nturn,Tab,&status2);
if (status2) {
- Get_NAFF(NTERM2, Nbtour, Tab, fx, fz, nb_freq);
- Get_freq(fx,fz,&nux1,&nuz1); // gets nux and nuz
+ Get_NAFF(NTERM2, Nbtour, Tab, fx, fy, nb_freq);
+ Get_freq(fx,fy,&nux1,&nuy1); // gets nux and nuy
if (diffusion) { // diffusion
Get_Tabshift(Tab,Tab0,Nbtour,Nbtour); // shift data for second round NAFF
- Get_NAFF(NTERM2, Nbtour, Tab0, fx2, fz2, nb_freq); // gets frequency vectors
- Get_freq(fx2,fz2,&nux2,&nuz2); // gets nux and nuz
+ Get_NAFF(NTERM2, Nbtour, Tab0, fx2, fy2, nb_freq); // gets frequency vectors
+ Get_freq(fx2,fy2,&nux2,&nuy2); // gets nux and nuy
}
} // unstable trajectory
else { //zeroing output
- nux1 = 0.0; nuz1 = 0.0;
- nux2 = 0.0; nuz2 = 0.0;
+ nux1 = 0.0; nuy1 = 0.0;
+ nux2 = 0.0; nuy2 = 0.0;
}
// printout value
if (!diffusion){
- fprintf(outf,"% 10.6e % 10.6e % 10.6e % 10.6e\n", dp, x, nux1, nuz1);
- //fprintf(stdout,"% 10.6e % 10.6e % 10.6e % 10.6e\n", dp, x, nux1, nuz1);
+ 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; dfz = nuz2 - nuz1;
+ dfx = nux2 - nux1; dfy = nuy2 - nuy1;
fprintf(outf,"% 10.6e % 10.6e % 10.6e % 10.6e % 10.6e % 10.6e\n",
- dp, x, nux1, nuz2, dfx, dfz);
+ dp, x, nux1, nuy2, dfx, dfy);
//fprintf(stdout,"% 10.6e % 10.6e % 10.6e % 10.6e % 10.6e % 10.6e\n",
- // dp, x, nux1, nuz2, dfx, dfz);
+ // dp, x, nux1, nuy2, dfx, dfy);
}
if(printloss)
@@ -1676,10 +1676,10 @@ void TunesShiftWithEnergy(const char *NudpFile,long Nb, long Nbtour, double emax
long i = 0L;
// long lastpos = 0L;
double Tab[DIM][NTURN];
- double fx[NTERM], fz[NTERM];
- double x = 0.0, xp = 0.0, z = 0.0, zp = 0.0, ctau = 0.0, dp = 0.0;
+ double fx[NTERM], fy[NTERM];
+ double x = 0.0, xp = 0.0, y = 0.0, zp = 0.0, ctau = 0.0, dp = 0.0;
double x0 = 1e-6, xp0 = 0.0, z0 = 1e-6, zp0 = 0.0, ctau0 = 0.0, dp0 = 0.0;
- double nux1 = 0.0, nuz1 = 0.0;
+ double nux1 = 0.0, nuy1 = 0.0;
int nb_freq[2] = {0, 0};
bool status = true;
struct tm *newtime;
@@ -1687,7 +1687,7 @@ void TunesShiftWithEnergy(const char *NudpFile,long Nb, long Nbtour, double emax
/* Get time and date */
newtime = GetTime();
- if (!trace) printf("\n Entering TunesShiftWithEnergy ...\n\n");
+ if (!trace) fprintf(stdout, "\n Entering TunesShiftWithEnergy ...\n\n");
/* Opening file */
if ((outf = fopen(NudpFile, "w")) == NULL) {
@@ -1696,8 +1696,8 @@ void TunesShiftWithEnergy(const char *NudpFile,long Nb, long Nbtour, double emax
}
fprintf(outf,"# TRACY III -- %s -- %s \n", NudpFile, asctime2(newtime));
- fprintf(outf,"# dP/P fx fz xcod pxcod zcod pzcod\n");
- if (trace) fprintf(stdout,"# dP/P fx fz xcod pxcod zcod pzcod\n");
+ fprintf(outf,"# dP/P fx fy xcod pxcod zcod pycod\n");
+ if (trace) fprintf(stdout,"# dP/P fx fy xcod pxcod zcod pycod\n");
if (Nb <= 1L)
fprintf(stdout,"NuDp: Error Nb=%ld\n",Nb);
@@ -1709,17 +1709,17 @@ void TunesShiftWithEnergy(const char *NudpFile,long Nb, long Nbtour, double emax
dp = dp0 + i*emax/(Nb-1)*2;
x = x0 ;
xp = xp0 ;
- z = z0 ;
+ y = z0 ;
zp = zp0 ;
ctau = ctau0;
- Trac_Simple4DCOD(x,xp,z,zp,dp,ctau,Nbtour,Tab,&status); // tracking around closed orbit
+ Trac_Simple4DCOD(x,xp,y,zp,dp,ctau,Nbtour,Tab,&status); // tracking around closed orbit
if (status) {
- Get_NAFF(NTERM, Nbtour, Tab, fx, fz, nb_freq); // get frequency vectors
- Get_freq(fx,fz,&nux1,&nuz1); // gets nux and nuz
+ Get_NAFF(NTERM, Nbtour, Tab, fx, fy, nb_freq); // get frequency vectors
+ Get_freq(fx,fy,&nux1,&nuy1); // gets nux and nuy
}
else {
- nux1 = 0.0; nuz1 = 0.0;
+ nux1 = 0.0; nuy1 = 0.0;
status = true;
}
@@ -1728,11 +1728,11 @@ void TunesShiftWithEnergy(const char *NudpFile,long Nb, long Nbtour, double emax
fprintf(outf,"%14.6e %14.6e %14.6e %14.6e %14.6e %14.6e %14.6e\n",
- dp, nux1,nuz1, globval.CODvect[0], globval.CODvect[1],
+ dp, nux1,nuy1, globval.CODvect[0], globval.CODvect[1],
globval.CODvect[2], globval.CODvect[3]);
if (trace) fprintf(stdout,"%14.6e %14.6e %14.6e %14.6e %14.6e %14.6e %14.6e\n",
- dp, nux1,nuz1, globval.CODvect[0], globval.CODvect[1],
+ dp, nux1,nuy1, globval.CODvect[0], globval.CODvect[1],
globval.CODvect[2], globval.CODvect[3]);
}
@@ -1797,7 +1797,7 @@ void Phase(const char *phasefile, double x,double xp,double y, double yp,double
fprintf(outf,"# TRACY III -- %s -- %s \n", fic, asctime2(newtime));
fprintf(outf,"# Phase Space \n");
fprintf(outf,
- "# x xp z zp dp ctau\n");
+ "# x xp y zp dp ctau\n");
// initialization to zero (case where unstable
for (i = 0; i < Nbtour; i++) {
@@ -1818,7 +1818,7 @@ void Phase(const char *phasefile, double x,double xp,double y, double yp,double
}
/****************************************************************************/
-/* void PhasePoly(long pos, double x0,double px0, double z0, double pz0, double delta0,
+/* void PhasePoly(long pos, double x0,double px0, double z0, double py0, double delta0,
double ctau0, long Nbtour)
Purpose:
@@ -1846,16 +1846,16 @@ void Phase(const char *phasefile, double x,double xp,double y, double yp,double
none
****************************************************************************/
-void PhasePoly(long pos, double x0,double px0, double z0, double pz0, double delta0,
+void PhasePoly(long pos, double x0,double px0, double z0, double py0, double delta0,
double ctau0, long Nbtour)
{
FILE *outf;
const char *fic="phasepoly.out";
long lastpos = 0,lastn = 0;
int i,j;
- double x, z, px, pz, delta, ctau;
+ double x, y, px, py, delta, ctau;
double ex = 1368E-9, el = 1.78E-4;
- double betax = 9.0, /*betaz = 8.2, */betal = 45.5;
+ double betax = 9.0, /*betay = 8.2, */betal = 45.5;
Vector xsynch;
int nx = 1, ne = 400;
struct tm *newtime;
@@ -1864,7 +1864,7 @@ void PhasePoly(long pos, double x0,double px0, double z0, double pz0, double del
newtime = GetTime();
fprintf(stdout,"Closed orbit:\n");
- fprintf(stdout," x px z pz delta ctau\n");
+ fprintf(stdout," x px y py delta ctau\n");
fprintf(stdout,"% 12.8f % 12.8f % 12.8f % 12.8f % 12.8f % 12.8f\n",
globval.CODvect[0], globval.CODvect[1], globval.CODvect[2],
globval.CODvect[3], globval.CODvect[4], globval.CODvect[5]);
@@ -1885,19 +1885,19 @@ void PhasePoly(long pos, double x0,double px0, double z0, double pz0, double del
fprintf(outf,"# TRACY III -- %s -- %s \n", fic, asctime2(newtime));
fprintf(outf,"# 6D Phase Space \n");
fprintf(outf,
- "# num x xp z zp dp ctau\n");
+ "# num x xp y zp dp ctau\n");
for (j = 0; j < ne; j++){
for (i = 0; i < nx; i++){
x = x0 + xsynch[0] + sqrt(ex*betax)*cos(2.0*M_PI/nx*i)*0;
px = px0 + xsynch[1] + sqrt(ex/betax)*sin(2.0*M_PI/nx*i)*0;
- z = z0 + xsynch[2];
- pz = pz0 + xsynch[3];
+ y = z0 + xsynch[2];
+ py = py0 + xsynch[3];
delta = delta0 + xsynch[4] + sqrt(el/betal)*sin(2*M_PI/ne*j)*0 ;
ctau = ctau0 + xsynch[5] + sqrt(el*betal)*cos(2*M_PI/ne*j)*0 + j*0.002;
fprintf(outf, "%6ld %+10.5e %+10.5e %+10.5e %+10.5e %+10.5e %+10.5e",
- 0L, x, px, z, pz, delta, ctau);
- Trac(x,px,z,pz,delta,ctau, Nbtour,pos, lastn, lastpos, outf);
+ 0L, x, px, y, py, delta, ctau);
+ Trac(x,px,y,py,delta,ctau, Nbtour,pos, lastn, lastpos, outf);
fprintf(outf,"\n");
}
}
@@ -1905,7 +1905,7 @@ void PhasePoly(long pos, double x0,double px0, double z0, double pz0, double del
}
/****************************************************************************/
-/* void PhasePortrait(double x0,double px0,double z0, double pz0, double delta0,
+/* void PhasePortrait(double x0,double px0,double z0, double py0, double delta0,
double end, double Nb, long Nbtour, int num)
Purpose:
@@ -1913,7 +1913,7 @@ void PhasePoly(long pos, double x0,double px0, double z0, double pz0, double del
Results in phaseportrait.out
Input:
- x0, px0, z0, Pz0, delta0, starting position
+ x0, px0, z0, py0, delta0, starting position
num cooordinate to vary (0 is x and 4 is delta)
end is the last value for the varying coordinate
Nb is the number of orbits to draw
@@ -1935,7 +1935,7 @@ void PhasePoly(long pos, double x0,double px0, double z0, double pz0, double del
Change of tracking routine: do not use a tabular to store data
****************************************************************************/
-void PhasePortrait(double x0,double px0,double z0, double pz0, double delta0,
+void PhasePortrait(double x0,double px0,double z0, double py0, double delta0,
double ctau0, double end, long Nb, long Nbtour, int num)
{
double Tab[6][NTURN];
@@ -1943,7 +1943,7 @@ void PhasePortrait(double x0,double px0,double z0, double pz0, double delta0,
const char fic[] = "phaseportrait.out";
int i = 0, j = 0;
double start = 0.0, step = 0.0;
- double x = 0.0, px = 0.0, z = 0.0, pz = 0.0, delta = 0.0, ctau = 0.0;
+ double x = 0.0, px = 0.0, y = 0.0, py = 0.0, delta = 0.0, ctau = 0.0;
bool status = true;
struct tm *newtime;
@@ -1961,10 +1961,10 @@ void PhasePortrait(double x0,double px0,double z0, double pz0, double delta0,
}
fprintf(outf,"# TRACY III -- %s \n", asctime2(newtime));
- fprintf(outf,"# x xp z zp dp ctau\n#\n");
+ fprintf(outf,"# x xp y zp dp ctau\n#\n");
x = x0; px = px0;
- z = z0; pz = pz0;
+ y = z0; py = py0;
delta = delta0;
switch (num) {
@@ -1975,7 +1975,7 @@ void PhasePortrait(double x0,double px0,double z0, double pz0, double delta0,
case 2:
start = z0; break;
case 3:
- start = pz0; break;
+ start = py0; break;
case 4:
start = delta0; break;
case 5:
@@ -1992,9 +1992,9 @@ void PhasePortrait(double x0,double px0,double z0, double pz0, double delta0,
case 1:
px = start + j*step; break;
case 2:
- z = start + j*step; break;
+ y = start + j*step; break;
case 3:
- pz = start + j*step; break;
+ py = start + j*step; break;
case 4:
delta = start + j*step; break;
case 5:
@@ -2002,8 +2002,8 @@ void PhasePortrait(double x0,double px0,double z0, double pz0, double delta0,
}
fprintf(stdout,"% .5e % .5e % .5e % .5e % .5e % .5e\n",
- x,px,z,pz,delta,ctau);
- Trac_Simple4DCOD(x,px,z,pz,delta,ctau,Nbtour,Tab,&status);
+ x,px,y,py,delta,ctau);
+ Trac_Simple4DCOD(x,px,y,py,delta,ctau,Nbtour,Tab,&status);
for (i = 0; i < Nbtour; i++) {
fprintf(outf,"% .5e % .5e % .5e % .5e % .5e % .5e\n",
Tab[0][i],Tab[1][i],Tab[2][i],Tab[3][i],Tab[4][i],Tab[5][i]);
@@ -2058,7 +2058,7 @@ void Check_Trac(double x, double px, double y, double py, double dp)
x1[2] = y; x1[3] = py;
x1[4] = dp; x1[5] = 0e0;
- fprintf(outf,"# i x xp z zp delta cT \n");
+ fprintf(outf,"# i x xp y zp delta cT \n");
for (i = 1; i<= globval.Cell_nLoc; i++)
{
@@ -2107,7 +2107,7 @@ void Enveloppe(double x, double px, double y, double py, double dp, double nturn
/* Get cod the delta = energy*/
getcod(dp, lastpos);
- printf("Enveloppe: xcod=%.5e mm zcod=% .5e mm \n", globval.CODvect[0]*1e3, globval.CODvect[2]*1e3);
+ fprintf(stdout, "Enveloppe: xcod=%.5e mm zcod=% .5e mm \n", globval.CODvect[0]*1e3, globval.CODvect[2]*1e3);
if ((outf = fopen(fic, "w")) == NULL)
{
@@ -2119,7 +2119,7 @@ void Enveloppe(double x, double px, double y, double py, double dp, double nturn
x1[2] = y + globval.CODvect[2]; x1[3] = py + globval.CODvect[3];
x1[4] = dp; x1[5] = 0e0;
- fprintf(outf,"# i x xp z zp delta cT \n");
+ fprintf(outf,"# i x xp y zp delta cT \n");
for (j = 1; j <= nturn; j++)
{
@@ -2130,7 +2130,7 @@ void Enveloppe(double x, double px, double y, double py, double dp, double nturn
Cell_Pass(i,i+1, x1, lastpos);
if (lastpos != i+1)
{
- printf("Enveloppe: Unstable motion ...\n"); exit_(1);
+ fprintf(stdout, "Enveloppe: Unstable motion ...\n"); exit_(1);
}
fprintf(outf,"%6.2f %+.5e %+.5e %+.5e %+.5e %+.5e %+.5e\n",
@@ -2216,7 +2216,7 @@ void Multipole_thicksext(const char *fic_hcorr, const char *fic_vcorr, const cha
- printf("Enter multipole ... \n");
+ fprintf(stdout, "Enter multipole ... \n");
/* Make lists of dipoles, quadrupoles and sextupoles */
for (i = 0; i <= globval.Cell_nLoc; i++)
@@ -2229,37 +2229,37 @@ void Multipole_thicksext(const char *fic_hcorr, const char *fic_vcorr, const cha
{
dlist[ndip] = i;
ndip++;
- if (trace) printf("%s % f\n",Cell.Elem.PName, Cell.Elem.M->PB[0 + HOMmax]);
+ if (trace) fprintf(stdout, "%s % f\n",Cell.Elem.PName, Cell.Elem.M->PB[0 + HOMmax]);
}
else if (Cell.Elem.M->PBpar[2L + HOMmax] != 0.0)
{
qlist[nquad] = i;
nquad++;
- if (trace) printf("%s % f\n",Cell.Elem.PName, Cell.Elem.M->PBpar[2L + HOMmax]);
+ if (trace) fprintf(stdout, "%s % f\n",Cell.Elem.PName, Cell.Elem.M->PBpar[2L + HOMmax]);
}
else if (Cell.Elem.M->PBpar[3L + HOMmax] != 0.0)
{
slist[nsext] = i;
nsext++;
- if (trace) printf("%s % f\n",Cell.Elem.PName, Cell.Elem.M->PBpar[3L + HOMmax]);
+ if (trace) fprintf(stdout, "%s % f\n",Cell.Elem.PName, Cell.Elem.M->PBpar[3L + HOMmax]);
}
else if ( Cell.Elem.PName[0] == 'c' && Cell.Elem.PName[1] == 'h')
{
hcorrlist[nhcorr] = i;
nhcorr++;
- if (trace) printf("%s \n",Cell.Elem.PName);
+ if (trace) fprintf(stdout, "%s \n",Cell.Elem.PName);
}
else if ( Cell.Elem.PName[0] == 'c' && Cell.Elem.PName[1] == 'v')
{
vcorrlist[nvcorr] = i;
nvcorr++;
- if (trace) printf("%s \n",Cell.Elem.PName);
+ if (trace) fprintf(stdout, "%s \n",Cell.Elem.PName);
}
else if ( Cell.Elem.PName[0] == 'q' && Cell.Elem.PName[1] == 't')
{
qcorrlist[nqcorr] = i;
nqcorr++;
- if (trace) printf("%s \n",Cell.Elem.PName);
+ if (trace) fprintf(stdout, "%s \n",Cell.Elem.PName);
}
}
}
@@ -2369,7 +2369,7 @@ void Multipole_thicksext(const char *fic_hcorr, const char *fic_vcorr, const cha
}
- if (!trace) printf("Elements: ndip=%d nquad=%d nsext=%d nhcorr=%d nvcorr=%d nqcorr=%d\n",
+ if (!trace) fprintf(stdout, "Elements: ndip=%d nquad=%d nsext=%d nhcorr=%d nvcorr=%d nqcorr=%d\n",
ndip,nquad,nsext,nhcorr,nvcorr,nqcorr);
/***********************************************************************************/
@@ -2405,37 +2405,37 @@ void Multipole_thicksext(const char *fic_hcorr, const char *fic_vcorr, const cha
mKL += dBoB2*theta*x0i;
SetKLpar(Cell.Fnum, Cell.Knum, mOrder=2L, mKL);
- if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld theta=% e mKl=% e\n",i,
+ if (trace) fprintf(stdout, "num= %4d name = %s Fnum = %3ld, Knum=%3ld theta=% e mKl=% e\n",i,
Cell.Elem.PName,Cell.Fnum, Cell.Knum, theta, mKL);
/* sextupole error */
mKL = dBoB3*theta*x02i;
SetKLpar(Cell.Fnum, Cell.Knum, mOrder=3L, mKL);
- if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld theta=% e mKl=% e\n",i,
+ if (trace) fprintf(stdout, "num= %4d name = %s Fnum = %3ld, Knum=%3ld theta=% e mKl=% e\n",i,
Cell.Elem.PName,Cell.Fnum, Cell.Knum, theta, mKL);
/* octupole error */
mKL = dBoB4*theta*x03i;
SetKLpar(Cell.Fnum, Cell.Knum, mOrder=4L, mKL);
- if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld theta=% e mKl=% e\n",i,
+ if (trace) fprintf(stdout, "num= %4d name = %s Fnum = %3ld, Knum=%3ld theta=% e mKl=% e\n",i,
Cell.Elem.PName,Cell.Fnum, Cell.Knum, theta, mKL);
/* decapole error */
mKL = dBoB5*theta*x04i;
SetKLpar(Cell.Fnum, Cell.Knum, mOrder=5L, mKL);
- if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld theta=% e mKl=% e\n",i,
+ if (trace) fprintf(stdout, "num= %4d name = %s Fnum = %3ld, Knum=%3ld theta=% e mKl=% e\n",i,
Cell.Elem.PName,Cell.Fnum, Cell.Knum, theta, mKL);
/* 12-pole error */
mKL = dBoB6*theta*x05i;
SetKLpar(Cell.Fnum, Cell.Knum, mOrder=6L, mKL);
- if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld theta=% e mKl=% e\n",i,
+ if (trace) fprintf(stdout, "num= %4d name = %s Fnum = %3ld, Knum=%3ld theta=% e mKl=% e\n",i,
Cell.Elem.PName,Cell.Fnum, Cell.Knum, theta, mKL);
/* 14-pole error */
mKL = dBoB7*theta*x06i;
SetKLpar(Cell.Fnum, Cell.Knum, mOrder=7L, mKL);
- if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld theta=% e mKl=% e\n",i,
+ if (trace) fprintf(stdout, "num= %4d name = %s Fnum = %3ld, Knum=%3ld theta=% e mKl=% e\n",i,
Cell.Elem.PName,Cell.Fnum, Cell.Knum, theta, mKL);
}
@@ -2482,7 +2482,7 @@ void Multipole_thicksext(const char *fic_hcorr, const char *fic_vcorr, const cha
else
mKL= b2*dBoB6C*x04i;
SetKLpar(Cell.Fnum, Cell.Knum, mOrder=6L, mKL);
- if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld b2=% e mKl=% e\n",i,
+ if (trace) fprintf(stdout, "num= %4d name = %s Fnum = %3ld, Knum=%3ld b2=% e mKl=% e\n",i,
Cell.Elem.PName,Cell.Fnum, Cell.Knum, b2, mKL);
/* 20-pole multipole error */
@@ -2491,7 +2491,7 @@ void Multipole_thicksext(const char *fic_hcorr, const char *fic_vcorr, const cha
else
mKL= b2*dBoB10C*x08i;
SetKLpar(Cell.Fnum, Cell.Knum, mOrder=10L, mKL);
- if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld b2=% e mKl=% e\n",i,
+ if (trace) fprintf(stdout, "num= %4d name = %s Fnum = %3ld, Knum=%3ld b2=% e mKl=% e\n",i,
Cell.Elem.PName,Cell.Fnum, Cell.Knum, b2, mKL);
/* 28-pole multipole error */
@@ -2500,7 +2500,7 @@ void Multipole_thicksext(const char *fic_hcorr, const char *fic_vcorr, const cha
else
mKL= b2*dBoB14C*x012i;
SetKLpar(Cell.Fnum, Cell.Knum, mOrder=14L, mKL);
- if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld b2=% e mKl=% e\n",i,
+ if (trace) fprintf(stdout, "num= %4d name = %s Fnum = %3ld, Knum=%3ld b2=% e mKl=% e\n",i,
Cell.Elem.PName,Cell.Fnum, Cell.Knum, b2, mKL);
/* sextupole mesure quadrupoles longs*/
@@ -2509,7 +2509,7 @@ void Multipole_thicksext(const char *fic_hcorr, const char *fic_vcorr, const cha
else
mKL= b2*dBoB3C*x1i;
SetKLpar(Cell.Fnum, Cell.Knum, mOrder=3L, mKL);
- if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld b2=% e mKl=% e\n",i,
+ if (trace) fprintf(stdout, "num= %4d name = %s Fnum = %3ld, Knum=%3ld b2=% e mKl=% e\n",i,
Cell.Elem.PName,Cell.Fnum, Cell.Knum, b2, mKL);
/* octupole mesure quadrupoles longs*/
@@ -2518,7 +2518,7 @@ void Multipole_thicksext(const char *fic_hcorr, const char *fic_vcorr, const cha
else
mKL= b2*dBoB4C*x1i*x1i;
SetKLpar(Cell.Fnum, Cell.Knum, mOrder=4L, mKL);
- if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld b2=% e mKl=% e\n",i,
+ if (trace) fprintf(stdout, "num= %4d name = %s Fnum = %3ld, Knum=%3ld b2=% e mKl=% e\n",i,
Cell.Elem.PName,Cell.Fnum, Cell.Knum, b2, mKL);
}
@@ -2564,37 +2564,37 @@ void Multipole_thicksext(const char *fic_hcorr, const char *fic_vcorr, const cha
/* 10-pole multipole error */
mKL= b3*dBoB5*x02i;
SetKLpar(Cell.Fnum, Cell.Knum, mOrder=5L, mKL);
- if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld b3=% e mKl=% e\n",i,
+ if (trace) fprintf(stdout, "num= %4d name = %s Fnum = %3ld, Knum=%3ld b3=% e mKl=% e\n",i,
Cell.Elem.PName,Cell.Fnum, Cell.Knum, b3, mKL);
/* 14-pole multipole error */
mKL= b3*dBoB7*x04i;
SetKLpar(Cell.Fnum, Cell.Knum, mOrder=7L, mKL);
- if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld b3=% e mKl=% e\n",i,
+ if (trace) fprintf(stdout, "num= %4d name = %s Fnum = %3ld, Knum=%3ld b3=% e mKl=% e\n",i,
Cell.Elem.PName,Cell.Fnum, Cell.Knum, b3, mKL);
/* 18-pole multipole error */
mKL= b3*dBoB9*x06i;
SetKLpar(Cell.Fnum, Cell.Knum, mOrder=9L, mKL);
- if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld b3=% e mKl=% e\n",i,
+ if (trace) fprintf(stdout, "num= %4d name = %s Fnum = %3ld, Knum=%3ld b3=% e mKl=% e\n",i,
Cell.Elem.PName,Cell.Fnum, Cell.Knum, b3, mKL);
/* 30-pole multipole error */
mKL= b3*dBoB15*x012i;
SetKLpar(Cell.Fnum, Cell.Knum, mOrder=15L, mKL);
- if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld b3=% e mKl=% e\n",i,
+ if (trace) fprintf(stdout, "num= %4d name = %s Fnum = %3ld, Knum=%3ld b3=% e mKl=% e\n",i,
Cell.Elem.PName,Cell.Fnum, Cell.Knum, b3, mKL);
/* 42-pole multipole error */
mKL= b3*dBoB21*x018i;
SetKLpar(Cell.Fnum, Cell.Knum, mOrder=21L, mKL);
- if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld b3=% e mKl=% e\n",i,
+ if (trace) fprintf(stdout, "num= %4d name = %s Fnum = %3ld, Knum=%3ld b3=% e mKl=% e\n",i,
Cell.Elem.PName,Cell.Fnum, Cell.Knum, b3, mKL);
/* 54-pole multipole error */
mKL= b3*dBoB27*x024i;
SetKLpar(Cell.Fnum, Cell.Knum, mOrder=27L, mKL);
- if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld b3=% e mKl=% e\n",i,
+ if (trace) fprintf(stdout, "num= %4d name = %s Fnum = %3ld, Knum=%3ld b3=% e mKl=% e\n",i,
Cell.Elem.PName,Cell.Fnum, Cell.Knum, b3, mKL);
}
@@ -2642,14 +2642,14 @@ void Multipole_thicksext(const char *fic_hcorr, const char *fic_vcorr, const cha
mKL += dBoB5*corr_strength*x04i;
SetKLpar(Cell.Fnum, Cell.Knum, mOrder=5L, mKL);
- if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld BL/brho=% e mKl=% e\n",i,
+ if (trace) fprintf(stdout, "num= %4d name = %s Fnum = %3ld, Knum=%3ld BL/brho=% e mKl=% e\n",i,
Cell.Elem.PName,Cell.Fnum, Cell.Knum, corr_strength, mKL);
/* 14-pole error */
mKL = GetKLpar(Cell.Fnum, Cell.Knum, mOrder=7L);
mKL += dBoB7*corr_strength*x06i;
SetKLpar(Cell.Fnum, Cell.Knum, mOrder=7L, mKL);
- if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld BL/brho=% e mKl=% e\n",i,
+ if (trace) fprintf(stdout, "num= %4d name = %s Fnum = %3ld, Knum=%3ld BL/brho=% e mKl=% e\n",i,
Cell.Elem.PName,Cell.Fnum, Cell.Knum, corr_strength, mKL);
/* 22-pole error */
@@ -2657,7 +2657,7 @@ void Multipole_thicksext(const char *fic_hcorr, const char *fic_vcorr, const cha
mKL += dBoB11*corr_strength*x010i;
SetKLpar(Cell.Fnum, Cell.Knum, mOrder=11, mKL);
- if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld BL/brho=% e mKl=% e\n",i,
+ if (trace) fprintf(stdout, "num= %4d name = %s Fnum = %3ld, Knum=%3ld BL/brho=% e mKl=% e\n",i,
Cell.Elem.PName,Cell.Fnum, Cell.Knum, corr_strength, mKL);
}
@@ -2706,20 +2706,20 @@ void Multipole_thicksext(const char *fic_hcorr, const char *fic_vcorr, const cha
// mKL = dBoB5*corr_strength*x04i;
// SetKLpar(Cell.Fnum, Cell.Knum, mOrder=-5L, mKL);
//
-// if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld BL/brho=% e mKl=% e\n",i,
+// if (trace) fprintf(stdout, "num= %4d name = %s Fnum = %3ld, Knum=%3ld BL/brho=% e mKl=% e\n",i,
// Cell.Elem.PName,Cell.Fnum, Cell.Knum, corr_strength, mKL);
// /* skew 14-pole error */
// mKL = dBoB7*corr_strength*x06i;
// SetKLpar(Cell.Fnum, Cell.Knum, mOrder=-7L, mKL);
//
-// if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld BL/brho=% e mKl=% e\n",i,
+// if (trace) fprintf(stdout, "num= %4d name = %s Fnum = %3ld, Knum=%3ld BL/brho=% e mKl=% e\n",i,
// Cell.Elem.PName,Cell.Fnum, Cell.Knum, corr_strength, mKL);
//
// /* skew 22-pole error */
// mKL = dBoB11*corr_strength*x010i;
// SetKLpar(Cell.Fnum, Cell.Knum, mOrder=-11L, mKL);
//
-// if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld BL/brho=% e mKl=% e\n",i,
+// if (trace) fprintf(stdout, "num= %4d name = %s Fnum = %3ld, Knum=%3ld BL/brho=% e mKl=% e\n",i,
// Cell.Elem.PName,Cell.Fnum, Cell.Knum, corr_strength, mKL);
// }
@@ -2733,7 +2733,7 @@ void Multipole_thicksext(const char *fic_hcorr, const char *fic_vcorr, const cha
mKL += dBoB5*corr_strength*x04i;
SetKLpar(Cell.Fnum, Cell.Knum, mOrder=-5L, mKL);
- if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld BL/brho=% e mKl=% e\n",i,
+ if (trace) fprintf(stdout, "num= %4d name = %s Fnum = %3ld, Knum=%3ld BL/brho=% e mKl=% e\n",i,
Cell.Elem.PName,Cell.Fnum, Cell.Knum, corr_strength, mKL);
/* skew 14-pole error */
@@ -2741,7 +2741,7 @@ void Multipole_thicksext(const char *fic_hcorr, const char *fic_vcorr, const cha
mKL += dBoB7*corr_strength*x06i;
SetKLpar(Cell.Fnum, Cell.Knum, mOrder=-7L, mKL);
- if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld BL/brho=% e mKl=% e\n",i,
+ if (trace) fprintf(stdout, "num= %4d name = %s Fnum = %3ld, Knum=%3ld BL/brho=% e mKl=% e\n",i,
Cell.Elem.PName,Cell.Fnum, Cell.Knum, corr_strength, mKL);
/* skew 22-pole error */
@@ -2749,7 +2749,7 @@ void Multipole_thicksext(const char *fic_hcorr, const char *fic_vcorr, const cha
mKL += dBoB11*corr_strength*x010i;
SetKLpar(Cell.Fnum, Cell.Knum, mOrder=-11L, mKL);
- if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld BL/brho=% e mKl=% e\n",i,
+ if (trace) fprintf(stdout, "num= %4d name = %s Fnum = %3ld, Knum=%3ld BL/brho=% e mKl=% e\n",i,
Cell.Elem.PName,Cell.Fnum, Cell.Knum, corr_strength, mKL);
}
/***********************************************************************************
@@ -2801,7 +2801,7 @@ d file */
// mKL = dBoB4*qcorr[i]*conv/brho*x02i;
// SetKLpar(Cell.Fnum, Cell.Knum, mOrder=-4L, mKL);
//
-// if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld BL/brho=% e mKl=% e\n",i,
+// if (trace) fprintf(stdout, "num= %4d name = %s Fnum = %3ld, Knum=%3ld BL/brho=% e mKl=% e\n",i,
// Cell.Elem.PName,Cell.Fnum, Cell.Knum, corr_strength, mKL);
// }
for (i = 0; i < nqcorr; i++)
@@ -2813,7 +2813,7 @@ d file */
mKL += dBoB4*qcorr[i]*conv/brho*x02i;
SetKLpar(Cell.Fnum, Cell.Knum, mOrder=-4L, mKL);
- if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld BL/brho=% e mKl=% e\n",i,
+ if (trace) fprintf(stdout, "num= %4d name = %s Fnum = %3ld, Knum=%3ld BL/brho=% e mKl=% e\n",i,
Cell.Elem.PName,Cell.Fnum, Cell.Knum, corr_strength, mKL);
}
}
@@ -2888,7 +2888,7 @@ void Multipole_thinsext(const char *fic_hcorr, const char *fic_vcorr, const char
FILE *fi;
/*********************************************************/
- printf("Enter multipole ... \n");
+ fprintf(stdout, "Enter multipole ... \n");
/* Make lists of dipoles, quadrupoles and sextupoles */
for (i = 0; i <= globval.Cell_nLoc; i++)
@@ -2901,42 +2901,42 @@ void Multipole_thinsext(const char *fic_hcorr, const char *fic_vcorr, const char
{
dlist[ndip] = i;
ndip++;
- if (trace) printf("%s % f\n",Cell.Elem.PName, Cell.Elem.M->PB[0 + HOMmax]);
+ if (trace) fprintf(stdout, "%s % f\n",Cell.Elem.PName, Cell.Elem.M->PB[0 + HOMmax]);
}
else if (fabs(Cell.Elem.M->PBpar[2L + HOMmax]) > 0.0)
{
qlist[nquad] = i;
nquad++;
- if (trace) printf("%s % f\n",Cell.Elem.PName, Cell.Elem.M->PBpar[2L + HOMmax]);
+ if (trace) fprintf(stdout, "%s % f\n",Cell.Elem.PName, Cell.Elem.M->PBpar[2L + HOMmax]);
}
else if (fabs(Cell.Elem.M->PBpar[3L + HOMmax]) > 0.0)
{
slist[nsext] = i;
nsext++;
- if (trace) printf("%s % f\n",Cell.Elem.PName, Cell.Elem.M->PBpar[3L + HOMmax]);
+ if (trace) fprintf(stdout, "%s % f\n",Cell.Elem.PName, Cell.Elem.M->PBpar[3L + HOMmax]);
}
else if ( Cell.Elem.PName[0] == 'c' && Cell.Elem.PName[1] == 'h')
{
hcorrlist[nhcorr] = i;
nhcorr++;
- if (trace) printf("%s \n",Cell.Elem.PName);
+ if (trace) fprintf(stdout, "%s \n",Cell.Elem.PName);
}
else if ( Cell.Elem.PName[0] == 'c' && Cell.Elem.PName[1] == 'v')
{
vcorrlist[nvcorr] = i;
nvcorr++;
- if (trace) printf("%s \n",Cell.Elem.PName);
+ if (trace) fprintf(stdout, "%s \n",Cell.Elem.PName);
}
else if ( Cell.Elem.PName[0] == 'q' && Cell.Elem.PName[1] == 't')
{
qcorrlist[nqcorr] = i;
nqcorr++;
- if (trace) printf("%s \n",Cell.Elem.PName);
+ if (trace) fprintf(stdout, "%s \n",Cell.Elem.PName);
}
}
}
- if (!trace) printf("Elements: ndip=%d nquad=%d nsext=%d nhcorr=%d nvcorr=%d nqcorr=%d\n",
+ if (!trace) fprintf(stdout, "Elements: ndip=%d nquad=%d nsext=%d nhcorr=%d nvcorr=%d nqcorr=%d\n",
ndip,nquad,nsext,nhcorr,nvcorr,nqcorr);
/***********************************************************************************/
@@ -2971,7 +2971,7 @@ void Multipole_thinsext(const char *fic_hcorr, const char *fic_vcorr, const char
mKL = dBoB2*theta*x0i;
SetKLpar(Cell.Fnum, Cell.Knum, mOrder=2L, mKL);
- if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld theta=% e mKl=% e\n",i,
+ if (trace) fprintf(stdout, "num= %4d name = %s Fnum = %3ld, Knum=%3ld theta=% e mKl=% e\n",i,
Cell.Elem.PName,Cell.Fnum, Cell.Knum, theta, mKL);
/* sextupole error */
@@ -3036,7 +3036,7 @@ void Multipole_thinsext(const char *fic_hcorr, const char *fic_vcorr, const char
else
mKL= b2*dBoB6C*x04i;
SetKLpar(Cell.Fnum, Cell.Knum, mOrder=6L, mKL);
- if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld b2=% e mKl=% e\n",i,
+ if (trace) fprintf(stdout, "num= %4d name = %s Fnum = %3ld, Knum=%3ld b2=% e mKl=% e\n",i,
Cell.Elem.PName,Cell.Fnum, Cell.Knum, b2, mKL);
/* 20-pole multipole error */
@@ -3045,7 +3045,7 @@ void Multipole_thinsext(const char *fic_hcorr, const char *fic_vcorr, const char
else
mKL= b2*dBoB10C*x08i;
SetKLpar(Cell.Fnum, Cell.Knum, mOrder=10L, mKL);
- if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld b2=% e mKl=% e\n",i,
+ if (trace) fprintf(stdout, "num= %4d name = %s Fnum = %3ld, Knum=%3ld b2=% e mKl=% e\n",i,
Cell.Elem.PName,Cell.Fnum, Cell.Knum, b2, mKL);
/* 28-pole multipole error */
@@ -3054,7 +3054,7 @@ void Multipole_thinsext(const char *fic_hcorr, const char *fic_vcorr, const char
else
mKL= b2*dBoB14C*x012i;
SetKLpar(Cell.Fnum, Cell.Knum, mOrder=14L, mKL);
- if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld b2=% e mKl=% e\n",i,
+ if (trace) fprintf(stdout, "num= %4d name = %s Fnum = %3ld, Knum=%3ld b2=% e mKl=% e\n",i,
Cell.Elem.PName,Cell.Fnum, Cell.Knum, b2, mKL);
/* sextupole mesure quadrupoles longs*/
@@ -3063,7 +3063,7 @@ void Multipole_thinsext(const char *fic_hcorr, const char *fic_vcorr, const char
else
mKL= b2*dBoB3C*x1i;
SetKLpar(Cell.Fnum, Cell.Knum, mOrder=3L, mKL);
- if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld b2=% e mKl=% e\n",i,
+ if (trace) fprintf(stdout, "num= %4d name = %s Fnum = %3ld, Knum=%3ld b2=% e mKl=% e\n",i,
Cell.Elem.PName,Cell.Fnum, Cell.Knum, b2, mKL);
/* octupole mesure quadrupoles longs*/
@@ -3072,7 +3072,7 @@ void Multipole_thinsext(const char *fic_hcorr, const char *fic_vcorr, const char
else
mKL= b2*dBoB4C*x1i*x1i;
SetKLpar(Cell.Fnum, Cell.Knum, mOrder=4L, mKL);
- if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld b2=% e mKl=% e\n",i,
+ if (trace) fprintf(stdout, "num= %4d name = %s Fnum = %3ld, Knum=%3ld b2=% e mKl=% e\n",i,
Cell.Elem.PName,Cell.Fnum, Cell.Knum, b2, mKL);
}
@@ -3113,25 +3113,25 @@ void Multipole_thinsext(const char *fic_hcorr, const char *fic_vcorr, const char
/* 18-pole multipole error */
mKL= b3*dBoB9*x06i;
SetKLpar(Cell.Fnum, Cell.Knum, mOrder=9L, mKL);
- if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld b3=% e mKl=% e\n",i,
+ if (trace) fprintf(stdout, "num= %4d name = %s Fnum = %3ld, Knum=%3ld b3=% e mKl=% e\n",i,
Cell.Elem.PName,Cell.Fnum, Cell.Knum, b3, mKL);
/* 30-pole multipole error */
mKL= b3*dBoB15*x012i;
SetKLpar(Cell.Fnum, Cell.Knum, mOrder=15L, mKL);
- if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld b3=% e mKl=% e\n",i,
+ if (trace) fprintf(stdout, "num= %4d name = %s Fnum = %3ld, Knum=%3ld b3=% e mKl=% e\n",i,
Cell.Elem.PName,Cell.Fnum, Cell.Knum, b3, mKL);
/* 42-pole multipole error */
mKL= b3*dBoB21*x018i;
SetKLpar(Cell.Fnum, Cell.Knum, mOrder=21L, mKL);
- if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld b3=% e mKl=% e\n",i,
+ if (trace) fprintf(stdout, "num= %4d name = %s Fnum = %3ld, Knum=%3ld b3=% e mKl=% e\n",i,
Cell.Elem.PName,Cell.Fnum, Cell.Knum, b3, mKL);
/* 54-pole multipole error */
mKL= b3*dBoB27*x024i;
SetKLpar(Cell.Fnum, Cell.Knum, mOrder=27L, mKL);
- if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld b3=% e mKl=% e\n",i,
+ if (trace) fprintf(stdout, "num= %4d name = %s Fnum = %3ld, Knum=%3ld b3=% e mKl=% e\n",i,
Cell.Elem.PName,Cell.Fnum, Cell.Knum, b3, mKL);
}
@@ -3179,20 +3179,20 @@ void Multipole_thinsext(const char *fic_hcorr, const char *fic_vcorr, const char
mKL = dBoB5*corr_strength*x04i;
SetKLpar(Cell.Fnum, Cell.Knum, mOrder=5L, mKL);
- if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld BL/brho=% e mKl=% e\n",i,
+ if (trace) fprintf(stdout, "num= %4d name = %s Fnum = %3ld, Knum=%3ld BL/brho=% e mKl=% e\n",i,
Cell.Elem.PName,Cell.Fnum, Cell.Knum, corr_strength, mKL);
/* 14-pole error */
mKL = dBoB7*corr_strength*x06i;
SetKLpar(Cell.Fnum, Cell.Knum, mOrder=7L, mKL);
- if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld BL/brho=% e mKl=% e\n",i,
+ if (trace) fprintf(stdout, "num= %4d name = %s Fnum = %3ld, Knum=%3ld BL/brho=% e mKl=% e\n",i,
Cell.Elem.PName,Cell.Fnum, Cell.Knum, corr_strength, mKL);
/* 22-pole error */
mKL = dBoB11*corr_strength*x010i;
SetKLpar(Cell.Fnum, Cell.Knum, mOrder=11, mKL);
- if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld BL/brho=% e mKl=% e\n",i,
+ if (trace) fprintf(stdout, "num= %4d name = %s Fnum = %3ld, Knum=%3ld BL/brho=% e mKl=% e\n",i,
Cell.Elem.PName,Cell.Fnum, Cell.Knum, corr_strength, mKL);
}
@@ -3242,20 +3242,20 @@ void Multipole_thinsext(const char *fic_hcorr, const char *fic_vcorr, const char
mKL = dBoB5*corr_strength*x04i;
SetKLpar(Cell.Fnum, Cell.Knum, mOrder=-5L, mKL);
- if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld BL/brho=% e mKl=% e\n",i,
+ if (trace) fprintf(stdout, "num= %4d name = %s Fnum = %3ld, Knum=%3ld BL/brho=% e mKl=% e\n",i,
Cell.Elem.PName,Cell.Fnum, Cell.Knum, corr_strength, mKL);
/* skew 14-pole error */
mKL = dBoB7*corr_strength*x06i;
SetKLpar(Cell.Fnum, Cell.Knum, mOrder=-7L, mKL);
- if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld BL/brho=% e mKl=% e\n",i,
+ if (trace) fprintf(stdout, "num= %4d name = %s Fnum = %3ld, Knum=%3ld BL/brho=% e mKl=% e\n",i,
Cell.Elem.PName,Cell.Fnum, Cell.Knum, corr_strength, mKL);
/* skew 22-pole error */
mKL = dBoB11*corr_strength*x010i;
SetKLpar(Cell.Fnum, Cell.Knum, mOrder=-11, mKL);
- if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld BL/brho=% e mKl=% e\n",i,
+ if (trace) fprintf(stdout, "num= %4d name = %s Fnum = %3ld, Knum=%3ld BL/brho=% e mKl=% e\n",i,
Cell.Elem.PName,Cell.Fnum, Cell.Knum, corr_strength, mKL);
}
@@ -3308,7 +3308,7 @@ d file */
mKL = dBoB4*qcorr[i]*conv/brho*x02i;
SetKLpar(Cell.Fnum, Cell.Knum, mOrder=-4L, mKL);
- if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld BL/brho=% e mKl=% e\n",i,
+ if (trace) fprintf(stdout, "num= %4d name = %s Fnum = %3ld, Knum=%3ld BL/brho=% e mKl=% e\n",i,
Cell.Elem.PName,Cell.Fnum, Cell.Knum, corr_strength, mKL);
}
}
@@ -3363,7 +3363,7 @@ void SetSkewQuad(void)
{
qlist[nquad] = i;
nquad++;
- if (trace) printf("%s % f\n",Cell.Elem.PName,
+ if (trace) fprintf(stdout, "%s % f\n",Cell.Elem.PName,
Cell.Elem.M->PBpar[2L + HOMmax]);
}
}
@@ -3400,7 +3400,7 @@ void SetSkewQuad(void)
mKL = b2*cos(2*theta[i]);
SetKLpar(Cell.Fnum, Cell.Knum, mOrder=2L, mKL);
- if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld KL=% e, KtiltL=% e\n"
+ if (trace) fprintf(stdout, "num= %4d name = %s Fnum = %3ld, Knum=%3ld KL=% e, KtiltL=% e\n"
,i,
Cell.Elem.PName,Cell.Fnum, Cell.Knum,
Cell.Elem.M->PBpar[HOMmax+2],
@@ -3461,7 +3461,7 @@ void SetSkewQuad(void)
// getelem(globval.hcorr_list[i], &Cell);
// SetKLpar(Cell.Fnum, Cell.Knum, mOrder, mKL[i]);
-// if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld KL=% e, KtiltL=% e\n"
+// if (trace) fprintf(stdout, "num= %4d name = %s Fnum = %3ld, Knum=%3ld KL=% e, KtiltL=% e\n"
// ,i,
// Cell.Elem.PName,Cell.Fnum, Cell.Knum,
// Cell.Elem.M->PBpar[HOMmax+mOrder],
@@ -3524,31 +3524,34 @@ void SetSkewQuad(void)
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
- 27/06/11 fix the bug of the index in the tabz and tabpz when calling Trac( );
+ 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
****************************************************************************/
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 zmax)
+ double em_max, long nstepm, long nturn, double ymax)
{
double dP = 0.0, dp1 = 0.0, dp2 = 0.0;
long lastpos = 0L,lastn = 0L;
long i = 0L, j = 0L, pos = 0L;
CellType Cell, Clost;
- double x = 0.0, px = 0.0, z = 0.0, pz = 0.0, ctau0 = 0.0, delta = 0.0;
+ double x = 0.0, px = 0.0, y = 0.0, py = 0.0, ctau0 = 0.0, delta = 0.0;
Vector x0;
FILE *outf2, *outf1;
- double **tabz0, **tabpz0;
+ double **taby0, **tabpy0;
struct tm *newtime; // for time
Vector codvector[Cell_nLocMax];
+ ss_vect<double> x1;
bool cavityflag, radiationflag;
+ bool phasespaceflag = false;
x0.zero();
@@ -3560,20 +3563,24 @@ void MomentumAcceptance(char *MomAccFile, long deb, long fin,
/* File opening for writing */
- outf1 = fopen("phase.out", "w");
+ outf1 = fopen("phase.out", "w"); // empty file if trace = false
outf2 = fopen(MomAccFile, "w");
+ //outf2_loss = fopen(MomAccFile, "w");
fprintf(outf2,"# TRACY III -- %s \n", asctime2(newtime));
- fprintf(outf2,"# i s dp s_lost name_lost \n#\n");
-
- fprintf(outf1,"# TRACY III -- %s \n", asctime2(newtime));
- fprintf(outf1,"# i x xp z zp dp ctau\n#\n");
-
-
+ fprintf(outf2,"# i s0[m] name0 dp Lostplane s_lost name_lost"
+ " x[m] px y[m] py delta ctau\n"
+ "# Element index and loss particle informations\n");
+
+ // empty file if trace = false
+ if (phasespaceflag){ fprintf(outf1,"# TRACY III -- %s \n", asctime2(newtime));
+ fprintf(outf1,"# i x px y py dp"
+ "ctau\n#\n");
+ }
pos = deb; /* starting position or element index in the ring */
/***************************************************************/
- fprintf(stdout,"Computing initial conditions ... \n");
+ fprintf(stdout,"MomentumAcceptance: Computing initial conditions ... \n");
/***************************************************************/
// cod search has to be done in 4D since in 6D it is zero
@@ -3583,19 +3590,19 @@ void MomentumAcceptance(char *MomAccFile, long deb, long fin,
globval.radiation = false; /* radiation on/off */
// Allocation of an array of pointer array
- tabz0 = (double **)malloc((nstepp)*sizeof(double*));
- tabpz0 = (double **)malloc((nstepp)*sizeof(double*));
- if (tabz0 == NULL || tabpz0 == NULL){
- fprintf(stdout,"1 out of memory \n"); return;
+ taby0 = (double **)malloc((nstepp)*sizeof(double*));
+ tabpy0 = (double **)malloc((nstepp)*sizeof(double*));
+ if (taby0 == NULL || tabpy0 == NULL){
+ fprintf(stdout,"MomentumAcceptance: 1 out of memory \n"); return;
}
for (i = 1L; i <= nstepp; i++){ // loop over energy
// Dynamical allocation 0 to nstepp -1
- tabz0[i-1L] = (double *)malloc((fin+1L)*sizeof(double));
- tabpz0[i-1L] = (double *)malloc((fin+1L)*sizeof(double));
- if (tabz0[i-1L] == NULL || tabpz0[i-1L] == NULL)
+ 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");
+ fprintf(stdout,"MomentumAcceptance: 2 out of memory \n");
return;
}
@@ -3609,15 +3616,15 @@ void MomentumAcceptance(char *MomAccFile, long deb, long fin,
set_vectorcod(codvector, dP);
// coordinates around closed orbit specially useful for 6D
- x0[0] = codvector[0][0];
- x0[1] = codvector[0][1];
- x0[2] = codvector[0][2] + zmax;
- x0[3] = codvector[0][3];
- x0[4] = codvector[0][4];
- x0[5] = codvector[0][5];
+ x0[0] = codvector[0][x_];
+ x0[1] = codvector[0][px_];
+ x0[2] = codvector[0][y_] + ymax;
+ x0[3] = codvector[0][py_];
+ x0[4] = codvector[0][delta_];
+ x0[5] = codvector[0][ct_];
if (0) fprintf(stdout,"dP=% e : %e %e %e %e %e %e\n",
- dP,x0[0],x0[1],x0[2],x0[3],x0[4],x0[5]);
+ dP,x0[x_],x0[px_],x0[y_],x0[py_],x0[delta_],x0[ct_]);
// Store vertical initial conditions
// case where deb is not element 1
if (deb > 1L)
@@ -3627,20 +3634,20 @@ void MomentumAcceptance(char *MomAccFile, long deb, long fin,
if (lastpos != j)
{ // look if stable
- tabz0 [i- 1L][j] = 1.0;
- tabpz0[i- 1L][j] = 1.0;
+ taby0 [i- 1L][j] = 1.0;
+ tabpy0[i- 1L][j] = 1.0;
}
else
{ // stable case
- tabz0 [i - 1L][j] = x0[2] - codvector[deb-1L][2];
- tabpz0[i - 1L][j] = x0[3] - codvector[deb-1L][3];
+ taby0 [i - 1L][j] = x0[y_] - codvector[deb-1L][y_];
+ tabpy0[i - 1L][j] = x0[py_] - codvector[deb-1L][py_];
}
}
else
{ // case where deb is element 1
j = deb - 1L;
- tabz0 [i - 1L][j] = x0[2] - codvector[j][2];
- tabpz0[i - 1L][j] = x0[3] - codvector[j][3];
+ taby0 [i - 1L][j] = x0[y_] - codvector[j][y_];
+ tabpy0[i - 1L][j] = x0[py_] - codvector[j][py_];
}
for (j = deb; j < fin; j++)
@@ -3649,13 +3656,13 @@ void MomentumAcceptance(char *MomAccFile, long deb, long fin,
// Cell_Pass(j -1L, j, x0, lastpos);
if (lastpos != j){ // look if stable
- tabz0 [i - 1L][j] = 1.0;
- tabpz0[i - 1L][j] = 1.0;
+ taby0 [i - 1L][j] = 1.0;
+ tabpy0[i - 1L][j] = 1.0;
}
else{ // stable case
- tabz0 [i - 1L][j] = x0[2] - codvector[j][2];
- tabpz0[i - 1L][j] = x0[3] - codvector[j][3];
-// fprintf(stdout,"z0= % e pz0= % e\n", tabz0 [i - 1L][j], tabpz0 [i - 1L][j]);
+ taby0 [i - 1L][j] = x0[y_] - codvector[j][y_];
+ tabpy0[i - 1L][j] = x0[py_] - codvector[j][py_];
+// fprintf(stdout,"z0= % e py0= % e\n", taby0 [i - 1L][j], tabpy0 [i - 1L][j]);
}
}
}
@@ -3664,7 +3671,7 @@ void MomentumAcceptance(char *MomAccFile, long deb, long fin,
globval.radiation = radiationflag;
/***************************************************************/
- fprintf(stdout,"Computing positive momentum acceptance ... \n");
+ fprintf(stdout,"MomentumAcceptance: Computing positive momentum acceptance ... \n");
/***************************************************************/
do
@@ -3673,14 +3680,14 @@ void MomentumAcceptance(char *MomAccFile, long deb, long fin,
getelem(pos,&Cell);
// coordinates around closed orbit which is non zero for 6D tracking
- x = Cell.BeamPos[0];
- px = Cell.BeamPos[1];
- z = Cell.BeamPos[2];
- pz = 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",
- pos, x, px, z, pz, delta, ctau0);
+ 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;
@@ -3696,34 +3703,43 @@ void MomentumAcceptance(char *MomAccFile, long deb, long fin,
else
dp2 = ep_max;
- if (trace) printf("i=%4ld pos=%4ld dp=%6.4g\n",i,pos,dp2);
- if (0) fprintf(stdout,"pos=%4ld z0 =% 10.5f pz0 =% 10.5f \n", pos, tabz0[i-1L][pos-1L], tabpz0[i-1L][pos-1L]);
+ 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]);
- Trac(x, px, z+tabz0[i-1L][pos-1L], pz+tabpz0[i-1L][pos-1L], dp2+delta , ctau0, nturn, pos, lastn, lastpos, outf1);
+ 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 != nstepp));
if ((lastn) == nturn)
dp1 = dp2;
+ // lost particle information
getelem(lastpos,&Clost);
getelem(pos,&Cell);
- fprintf(stdout,"pos=%4ld z0 =% 10.5f pz0 =% 10.5f \n", pos, tabz0[i-1L][pos-1L], tabpz0[i-1L][pos-1L]);
- 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.5f %10.5f %10.5f %*s\n", pos,Cell.S,dp1,Clost.S,5,Clost.Elem.PName);
-
+ if (trace) fprintf(stdout,"pos=%4ld zi =%+10.5f pyi =%+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",
+ 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",
+ pos, Cell.S, Cell.Elem.PName, dp1, status.lossplane,
+ Clost.S, Clost.Elem.PName,
+ x1[x_], x1[px_], x1[y_], x1[py_], x1[delta_], x1[ct_]);
pos++;
}while(pos != fin);
// free memory
for (i = 1L; i <= nstepp; i++){
- free(tabz0 [i - 1L]);
- free(tabpz0[i - 1L]);
+ free(taby0 [i - 1L]);
+ free(tabpy0[i - 1L]);
}
- free(tabz0);
- free(tabpz0);
+ free(taby0);
+ free(tabpy0);
+ fflush(outf2);
/***************************************************************/
/***************************************************************/
@@ -3736,7 +3752,7 @@ void MomentumAcceptance(char *MomAccFile, long deb, long fin,
pos = deb; /* starting position in the ring */
/***************************************************************/
- fprintf(stdout,"Computing initial conditions ... \n");
+ fprintf(stdout,"MomentumAcceptance: Computing initial conditions ... \n");
/***************************************************************/
// cod search has to be done in 4D since in 6D it is zero
@@ -3746,17 +3762,17 @@ void MomentumAcceptance(char *MomAccFile, long deb, long fin,
globval.radiation = false; /* radiation on/off */
// Allocation of an array of pointer array
- tabz0 = (double **)malloc((nstepm)*sizeof(double*));
- tabpz0 = (double **)malloc((nstepm)*sizeof(double*));
- if (tabz0 == NULL || tabpz0 == NULL){
+ 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;
}
for (i = 1L; i <= nstepm; i++){ // loop over energy
// Dynamical allocation
- tabz0[i-1L] = (double *)malloc((fin+1L)*sizeof(double));
- tabpz0[i-1L] = (double *)malloc((fin+1L)*sizeof(double));
- if (tabz0[i-1L] == NULL || tabpz0[i-1L] == NULL){
+ 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;
}
@@ -3771,12 +3787,12 @@ void MomentumAcceptance(char *MomAccFile, long deb, long fin,
set_vectorcod(codvector, dP);
// coordinates around closed orbit specially usefull for 6D
- x0[0] = codvector[0][0];
- x0[1] = codvector[0][1];
- x0[2] = codvector[0][2] + zmax;
- x0[3] = codvector[0][3];
- x0[4] = codvector[0][4];
- x0[5] = codvector[0][5];
+ x0[x_] = codvector[0][x_];
+ x0[px_] = codvector[0][px_];
+ x0[y_] = codvector[0][y_] + ymax;
+ x0[py_] = codvector[0][py_];
+ x0[delta_] = codvector[0][delta_];
+ x0[ct_] = codvector[0][ct_];
// Store vertical initial conditions
// case where deb is not element 1
@@ -3784,32 +3800,33 @@ void MomentumAcceptance(char *MomAccFile, long deb, long fin,
Cell_Pass(1L, deb - 1L, x0, lastpos); // track from 1 to deb-1L element
j = deb -1L;
if (lastpos != j){ // look if stable
- tabz0 [i- 1L][j] = 1.0;
- tabpz0[i- 1L][j] = 1.0;
+ taby0 [i- 1L][j] = 1.0;
+ tabpy0[i- 1L][j] = 1.0;
}
else{ // stable case
- tabz0 [i - 1L][j] = x0[2] - codvector[deb-1L][2];
- tabpz0[i - 1L][j] = x0[3] - codvector[deb-1L][3];
+ taby0 [i - 1L][j] = x0[y_] - codvector[deb-1L][y_];
+ tabpy0[i - 1L][j] = x0[py_] - codvector[deb-1L][py_];
}
}
else { // case where deb is element 1
j = deb - 1L;
- tabz0 [i - 1L][j] = x0[2] - codvector[j][2];
- tabpz0[i - 1L][j] = x0[3] - codvector[j][3];
-// fprintf(stdout,"z0= % e pz0= % e\n", tabz0 [i - 1L][j], tabpz0 [i - 1L][j]);
+ taby0 [i - 1L][j] = x0[y_] - codvector[j][y_];
+ tabpy0[i - 1L][j] = x0[py_] - codvector[j][py_];
+// fprintf(stdout,"z0= % e py0= % e\n", taby0 [i - 1L][j], tabpy0 [i - 1L][j]);
}
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
- tabz0 [i - 1L][j] = 1.0;
- tabpz0[i - 1L][j] = 1.0;
+ taby0 [i - 1L][j] = 1.0;
+ tabpy0[i - 1L][j] = 1.0;
}
else{ // stable case
- tabz0 [i - 1L][j] = x0[2] - codvector[j][2];
- tabpz0[i - 1L][j] = x0[3] - codvector[j][3];
-// fprintf(stdout,"dP= % e pos= %ld z0= % e pz0= % e\n", dP, j, tabz0 [i - 1L][j], tabpz0 [i - 1L][j]);
+ taby0 [i - 1L][j] = x0[y_] - codvector[j][y_];
+ tabpy0[i - 1L][j] = x0[py_] - codvector[j][py_];
+// fprintf(stdout,"dP= % e pos= %ld z0= % e py0= % e\n", dP, j,
+// taby0 [i - 1L][j], tabpy0 [i - 1L][j]);
}
}
}
@@ -3818,7 +3835,7 @@ void MomentumAcceptance(char *MomAccFile, long deb, long fin,
globval.radiation = radiationflag;
/***************************************************************/
- fprintf(stdout,"Computing negative momentum acceptance ... \n");
+ fprintf(stdout,"MomentumAcceptance: Computing negative momentum acceptance ... \n");
/***************************************************************/
do {
@@ -3826,14 +3843,14 @@ void MomentumAcceptance(char *MomAccFile, long deb, long fin,
getelem(pos,&Cell);
// coordinates around closed orbit which is non zero for 6D tracking
- x = Cell.BeamPos[0];
- px = Cell.BeamPos[1];
- z = Cell.BeamPos[2];
- pz = 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",
- pos, x, px, z, pz, delta, ctau0);
+ 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;
@@ -3848,8 +3865,9 @@ void MomentumAcceptance(char *MomAccFile, long deb, long fin,
else {
dp2 = em_max;
}
- if (trace) printf("i=%4ld pos=%4ld dp=%6.4g\n",i,pos,dp2);
- Trac(x, px, z+tabz0[i-1L][pos-1L], pz+tabpz0[i-1L][pos-1L], dp2+delta , ctau0, nturn, pos, lastn, lastpos, outf1);
+ if (trace) fprintf(stdout, "i=%4ld pos=%4ld dp=%6.4g\n",i,pos,dp2);
+ 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));
@@ -3857,21 +3875,30 @@ void MomentumAcceptance(char *MomAccFile, long deb, long fin,
getelem(lastpos,&Clost);
getelem(pos,&Cell);
- if (!trace) printf("i=%4ld pos=%4ld dp=%6.4g\n",i,pos,dp2);
- fprintf(stdout,"pos=%4ld z0 =% 10.5f pz0 =% 10.5f \n", pos, tabz0[i-1L][pos-1L], tabpz0[i-1L][pos-1L]);
- 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.5f %10.5f %10.5f %*s\n", pos,Cell.S,dp1,Clost.S, 5, Clost.Elem.PName);
+ if (trace){
+ fprintf(stdout, "i=%4ld pos=%4ld dp=%6.4g\n",i,pos,dp2);
+ fprintf(stdout,"pos=%4ld y0 =%+10.5f py0 =%+10.5f \n", pos, taby0[i-1L][pos-1L],
+ tabpy0[i-1L][pos-1L]);
+ 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",
+ pos, Cell.S, Cell.Elem.PName, dp1, status.lossplane,
+ Clost.S, Clost.Elem.PName,
+ x1[x_], x1[px_], x1[y_], x1[py_], x1[delta_], x1[ct_]);
pos++;
}
while(pos != fin);
// free memory
for (i = 1L; i <= nstepp; i++){
- free(tabz0 [i - 1L]);
- free(tabpz0[i - 1L]);
+ free(taby0 [i - 1L]);
+ free(tabpy0[i - 1L]);
}
- free(tabz0);
- free(tabpz0);
+ free(taby0);
+ free(tabpy0);
+ fflush(outf2);
fflush(NULL); // force writing at the end (BUG??)
fclose(outf1);
@@ -3881,7 +3908,7 @@ void MomentumAcceptance(char *MomAccFile, long deb, long fin,
/****************************************************************************/
/*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 zmax, int numprocs,int myid)
+ long nstepm, long nturn, double ymax, int numprocs,int myid)
Purpose:
Parallel version of MomentumAcceptance( ).
@@ -3932,28 +3959,30 @@ Purpose:
****************************************************************************/
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 zmax, 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;
long i = 0L, j = 0L, pos = 0L;
CellType Cell, Clost;
- double x = 0.0, px = 0.0, z = 0.0, pz = 0.0, ctau0 = 0.0, delta = 0.0;
+ double x = 0.0, px = 0.0, y = 0.0, py = 0.0, ctau0 = 0.0, delta = 0.0;
Vector x0;
FILE *outf2, *outf1;
- double **tabz0, **tabpz0;
+ double **taby0, **tabpy0;
struct tm *newtime; // for time
Vector codvector[Cell_nLocMax];
bool cavityflag, radiationflag;
+ ss_vect<double> x1;
+
x0.zero();
// Get time and date
newtime = GetTime();
/************************/
- // Fin des declarations
+ // End of definitions
// File opening for writing
char PhaseFile[50];
@@ -3966,21 +3995,22 @@ void MomentumAcceptance_p(char *_MomAccFile, long deb, long fin, double ep_min,
if(myid==0)
{
fprintf(outf1,"# TRACY III -- %s \n", asctime2(newtime));
- fprintf(outf1,"# i x xp z zp dp ctau\n#\n");
+ fprintf(outf1,"# i x xp y zp dp ctau\n#\n");
}
char MomAccFile[50];
MomAccFile[0]='\0';
sprintf(MomAccFile,"%d",myid);
strcat(MomAccFile,_MomAccFile);
- printf("%s\n",MomAccFile);
+ fprintf(stdout, "%s\n",MomAccFile);
outf2 = fopen(MomAccFile, "w");
- if(myid==0)
- {
- fprintf(outf2,"# TRACY III -- %s \n", asctime2(newtime));
- fprintf(outf2,"# i s dp s_lost name_lost \n#\n");
+ if(myid==0){
+ fprintf(outf2,"# TRACY III -- %s \n", asctime2(newtime));
+ fprintf(outf2,"# i s0[m] name0 dp Lostplane s_lost name_lost"
+ " x[m] px y[m] py delta ctau\n"
+ "# Element index and loss particle informations\n");
}
// pos = deb; // starting position or element index in the ring
@@ -3996,18 +4026,18 @@ void MomentumAcceptance_p(char *_MomAccFile, long deb, long fin, double ep_min,
globval.radiation = false; /* radiation on/off */
// Memory allocation. Allocation of an array of pointer array
- tabz0 = (double **)malloc((nstepp)*sizeof(double*));
- tabpz0 = (double **)malloc((nstepp)*sizeof(double*));
- if (tabz0 == NULL || tabpz0 == NULL){
+ 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;
}
for (i = 1L; i <= nstepp; i++)
{ // loop over energy
// Dynamical allocation 0 to nstepp -1
- tabz0[i-1L] = (double *)malloc((fin+1L)*sizeof(double));
- tabpz0[i-1L] = (double *)malloc((fin+1L)*sizeof(double));
- if (tabz0[i-1L] == NULL || tabpz0[i-1L] == NULL)
+ 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;
@@ -4023,14 +4053,15 @@ void MomentumAcceptance_p(char *_MomAccFile, long deb, long fin, double ep_min,
set_vectorcod(codvector, dP);
// coordinates around closed orbit specially useful for 6D
- x0[0] = codvector[0][0];
- x0[1] = codvector[0][1];
- x0[2] = codvector[0][2] + zmax;
- x0[3] = codvector[0][3];
- x0[4] = codvector[0][4];
- x0[5] = codvector[0][5];
+ x0[x_] = codvector[0][x_];
+ x0[px_] = codvector[0][px_];
+ x0[y_] = codvector[0][y_] + ymax;
+ x0[py_] = codvector[0][py_];
+ x0[delta_] = codvector[0][delta_];
+ x0[ct_] = codvector[0][ct_];
- if (0) fprintf(stdout,"dP=% e : %e %e %e %e %e %e\n", dP,x0[0],x0[1],x0[2],x0[3],x0[4],x0[5]);
+ if (0) fprintf(stdout,"dP=% e : %e %e %e %e %e %e\n",
+ dP,x0[0],x0[1],x0[2],x0[3],x0[4],x0[5]);
// Store vertical initial conditions
// case where deb is not element 1
@@ -4041,20 +4072,20 @@ void MomentumAcceptance_p(char *_MomAccFile, long deb, long fin, double ep_min,
if (lastpos != j)
{ // look if stable
- tabz0 [i- 1L][j] = 1.0;
- tabpz0[i- 1L][j] = 1.0;
+ taby0 [i- 1L][j] = 1.0;
+ tabpy0[i- 1L][j] = 1.0;
}
else
{ // stable case
- tabz0 [i - 1L][j] = x0[2] - codvector[deb-1L][2];
- tabpz0[i - 1L][j] = x0[3] - codvector[deb-1L][3];
+ 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
j = deb - 1L;
- tabz0 [i - 1L][j] = x0[2] - codvector[j][2];
- tabpz0[i - 1L][j] = x0[3] - codvector[j][3];
+ taby0 [i - 1L][j] = x0[2] - codvector[j][2];
+ tabpy0[i - 1L][j] = x0[3] - codvector[j][3];
}
for (j = deb; j < fin; j++)
@@ -4063,13 +4094,13 @@ void MomentumAcceptance_p(char *_MomAccFile, long deb, long fin, double ep_min,
//Cell_Pass(j -1L, j, x0, lastpos);
if (lastpos != j){ // look if stable
- tabz0 [i - 1L][j] = 1.0;
- tabpz0[i - 1L][j] = 1.0;
+ taby0 [i - 1L][j] = 1.0;
+ tabpy0[i - 1L][j] = 1.0;
}
else{ // stable case
- tabz0 [i - 1L][j] = x0[2] - codvector[j][2];
- tabpz0[i - 1L][j] = x0[3] - codvector[j][3];
-// fprintf(stdout,"z0= % e pz0= % e\n", tabz0 [i - 1L][j], tabpz0 [i - 1L][j]);
+ 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]);
}
}
}
@@ -4088,14 +4119,16 @@ void MomentumAcceptance_p(char *_MomAccFile, long deb, long fin, double ep_min,
//the end element should not less than start element
if(fin < deb){
- printf("End element index %ld should be NOT smaller than the start element index %ld\n",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);
}
integer=(fin-deb+1)/numprocs;
residue=(fin-deb+1)-integer*numprocs;
- printf("myid:%d, integer:%d, resideu:%d, numprocs:%d\n",myid,integer,residue,numprocs);
+ 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
@@ -4118,17 +4151,19 @@ if(fin < deb){
getelem(pos,&Cell);
// coordinates around closed orbit which is non zero for 6D tracking
- x = Cell.BeamPos[0];
- px = Cell.BeamPos[1];
- z = Cell.BeamPos[2];
- pz = Cell.BeamPos[3];
- delta = Cell.BeamPos[4];
- ctau0 = Cell.BeamPos[5];
+ x = Cell.BeamPos[x_];
+ px = Cell.BeamPos[px_];
+ y = Cell.BeamPos[y_];
+ py = Cell.BeamPos[py_];
+ delta = Cell.BeamPos[delta_];
+ ctau0 = Cell.BeamPos[ct_];
- fprintf(stdout,"%3ld %6.4g %6.4g %6.4g %6.4g %6.4g %6.4g\n", pos, x, px, z, pz, delta, ctau0);
+ 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);
- 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]);
+ 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]);
dp1 = 0.0;
dp2 = 0.0;
@@ -4144,13 +4179,17 @@ if(fin < deb){
else
dp2 = ep_max;
- // if (trace) printf("i=%4ld pos=%4ld dp=%6.4g\n",i,pos,dp2);
+ // if (trace) fprintf(stdout, "i=%4ld pos=%4ld dp=%6.4g\n",i,pos,dp2);
if (trace)
- printf("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, z+tabz0[i-1L][pos-1L], pz+tabpz0[i-1L][pos-1L], dp2+delta, ctau0);
+ 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],
+ py+tabpy0[i-1L][pos-1L], dp2+delta, ctau0);
- if (0) fprintf(stdout,"pos=%4ld z0 =% 10.5f pz0 =% 10.5f \n", pos, tabz0[i-1L][pos-1L], tabpz0[i-1L][pos-1L]);
+ if (trace) fprintf(stdout,"pos=%4ld z0 =% 10.5f py0 =% 10.5f \n",
+ pos, taby0[i-1L][pos-1L], tabpy0[i-1L][pos-1L]);
- Trac(x, px, z+tabz0[i-1L][pos-1L], pz+tabpz0[i-1L][pos-1L], dp2+delta , ctau0, nturn, pos, lastn, lastpos, outf1);
+ Trac(x, px, y+taby0[i-1L][pos-1L], py+tabpy0[i-1L][pos-1L], dp2+delta ,
+ ctau0, nturn, pos, lastn, lastpos, outf1);
}while (((lastn) == nturn) && (i != nstepp));
@@ -4160,22 +4199,25 @@ if(fin < deb){
getelem(lastpos,&Clost);
getelem(pos,&Cell);
- fprintf(stdout,"pos=%4ld z0 =% 10.5f pz0 =% 10.5f \n", pos, tabz0[i-1L][pos-1L], tabpz0[i-1L][pos-1L]);
- 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.5f %10.5f %10.5f %*s\n", pos, Cell.S, dp1, Clost.S, 5, Clost.Elem.PName);
-
- // pos++;
-
+ 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",
+ 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",
+ pos, Cell.S, Cell.Elem.PName, dp1, status.lossplane,
+ Clost.S, Clost.Elem.PName,
+ x1[x_], x1[px_], x1[y_], x1[py_], x1[delta_], x1[ct_]);
}//while(pos != fin);
// free memory
for (i = 1L; i <= nstepp; i++)
{
- free(tabz0 [i - 1L]);
- free(tabpz0[i - 1L]);
+ free(taby0 [i - 1L]);
+ free(tabpy0[i - 1L]);
}
- free(tabz0);
- free(tabpz0);
+ free(taby0);
+ free(tabpy0);
/***************************************************************/
/***************************************************************/
@@ -4198,17 +4240,17 @@ if(fin < deb){
globval.radiation = false; // radiation on/off
// Allocation of an array of pointer array
- tabz0 = (double **)malloc((nstepm)*sizeof(double*));
- tabpz0 = (double **)malloc((nstepm)*sizeof(double*));
- if (tabz0 == NULL || tabpz0 == NULL){
+ 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;
}
for (i = 1L; i <= nstepm; i++){ // loop over energy
// Dynamical allocation
- tabz0[i-1L] = (double *)malloc((fin+1L)*sizeof(double));
- tabpz0[i-1L] = (double *)malloc((fin+1L)*sizeof(double));
- if (tabz0[i-1L] == NULL || tabpz0[i-1L] == NULL){
+ 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;
}
@@ -4223,12 +4265,12 @@ if(fin < deb){
set_vectorcod(codvector, dP);
// coordinates around closed orbit specially usefull for 6D
- x0[0] = codvector[0][0];
- x0[1] = codvector[0][1];
- x0[2] = codvector[0][2] + zmax;
- x0[3] = codvector[0][3];
- x0[4] = codvector[0][4];
- x0[5] = codvector[0][5];
+ x0[x_] = codvector[0][x_];
+ x0[px_] = codvector[0][px_];
+ x0[y_] = codvector[0][y_] + ymax;
+ x0[py_] = codvector[0][py_];
+ x0[delta_] = codvector[0][delta_];
+ x0[ct_] = codvector[0][ct_];
// Store vertical initial conditions
// case where deb is not element 1
@@ -4236,32 +4278,33 @@ if(fin < deb){
Cell_Pass(1L, deb - 1L, x0, lastpos); // track from 1 to deb-1L element
j = deb -1L;
if (lastpos != j){ // look if stable
- tabz0 [i- 1L][j] = 1.0;
- tabpz0[i- 1L][j] = 1.0;
+ taby0 [i- 1L][j] = 1.0;
+ tabpy0[i- 1L][j] = 1.0;
}
else{ // stable case
- tabz0 [i - 1L][j] = x0[2] - codvector[deb-1L][2];
- tabpz0[i - 1L][j] = x0[3] - codvector[deb-1L][3];
+ 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
j = deb - 1L;
- tabz0 [i - 1L][j] = x0[2] - codvector[j][2];
- tabpz0[i - 1L][j] = x0[3] - codvector[j][3];
-// fprintf(stdout,"z0= % e pz0= % e\n", tabz0 [i - 1L][j], tabpz0 [i - 1L][j]);
+ 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
Cell_Pass(j, j, x0, lastpos);
// Cell_Pass(j -1L, j, x0, lastpos);
if (lastpos != j){ // look if stable
- tabz0 [i - 1L][j] = 1.0;
- tabpz0[i - 1L][j] = 1.0;
+ taby0 [i - 1L][j] = 1.0;
+ tabpy0[i - 1L][j] = 1.0;
}
else{ // stable case
- tabz0 [i - 1L][j] = x0[2] - codvector[j][2];
- tabpz0[i - 1L][j] = x0[3] - codvector[j][3];
-// fprintf(stdout,"dP= % e pos= %ld z0= % e pz0= % e\n", dP, j, tabz0 [i - 1L][j], tabpz0 [i - 1L][j]);
+ 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]);
}
}
}
@@ -4282,11 +4325,12 @@ if(fin < deb){
// coordinates around closed orbit which is non zero for 6D tracking
x = Cell.BeamPos[0];
px = Cell.BeamPos[1];
- z = Cell.BeamPos[2];
- pz = Cell.BeamPos[3];
+ 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", pos, x, px, z, pz, delta, ctau0);
+ 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;
@@ -4302,8 +4346,10 @@ if(fin < deb){
dp2 = em_max;
}
if (trace)
- printf("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, z+tabz0[i-1L][pos-1L], pz+tabpz0[i-1L][pos-1L], dp2+delta, ctau0);
- Trac(x, px, z+tabz0[i-1L][pos-1L], pz+tabpz0[i-1L][pos-1L], dp2+delta , ctau0, nturn, pos, lastn, lastpos, outf1);
+ 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);
}
while (((lastn) == nturn) && (i != nstepm));
@@ -4311,10 +4357,16 @@ if(fin < deb){
getelem(lastpos,&Clost);
getelem(pos,&Cell);
- if (!trace) printf("i=%4ld pos=%4ld dp=%6.4g\n",i,pos,dp2);
- fprintf(stdout,"pos=%4ld z0 =% 10.5f pz0 =% 10.5f \n", pos, tabz0[i-1L][pos-1L], tabpz0[i-1L][pos-1L]);
- 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.5f %10.5f %10.5f %*s\n", pos, Cell.S, dp1, Clost.S, 5, Clost.Elem.PName);
+ 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",
+ 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",
+ pos, Cell.S, Cell.Elem.PName, dp1, status.lossplane,
+ Clost.S, Clost.Elem.PName,
+ x1[x_], x1[px_], x1[y_], x1[py_], x1[delta_], x1[ct_]);
// pos++;
}
//while(pos != fin);
@@ -4322,11 +4374,11 @@ if(fin < deb){
// free memory
for (i = 1L; i <= nstepp; i++)
{
- free(tabz0 [i - 1L]);
- free(tabpz0[i - 1L]);
+ free(taby0 [i - 1L]);
+ free(tabpy0[i - 1L]);
}
- free(tabz0);
- free(tabpz0);
+ free(taby0);
+ free(tabpy0);
fflush(NULL); // force writing at the end (BUG??)
fclose(outf1);
@@ -4386,11 +4438,11 @@ void set_vectorcod(Vector codvector[], double dP)
// LAURENT
/****************************************************************************/
-/* void spectrum(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
+/* void spectrum(long Nbx, long Nby, long Nbtour, double xmax, double ymax,
double energy, bool *status)
Purpose:
- Compute a frequency map of Nbx x Nbz points
+ Compute a frequency map of Nbx x Nby points
For each set of initial conditions the particle is tracked over
Nbtour for an energy offset dp
@@ -4402,7 +4454,7 @@ void set_vectorcod(Vector codvector[], double dP)
Nbx horizontal step number
Nby vertical step number
xmax horizontal maximum amplitude
- zmax vertical maximum amplitude
+ ymax vertical maximum amplitude
Nbtour number of turn for tracking
energy particle energy offset
@@ -4423,7 +4475,7 @@ void set_vectorcod(Vector codvector[], double dP)
15/10/03 run for the diffusion: nasty patch for retrieving the closed orbit
****************************************************************************/
-void spectrum(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
+void spectrum(long Nbx, long Nby, long Nbtour, double xmax, double ymax,
double energy, bool diffusion)
{
FILE *xoutf, *zoutf;
@@ -4431,8 +4483,8 @@ void spectrum(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
const char zfic[] = "zspectrum.out";
long i, j, k;
#define nterm2 20
- double Tab[6][NTURN], fx[nterm2], fz[nterm2], fx2[nterm2], fz2[nterm2];
- double x = 0.0, xp = 0.0, z = 0.0, zp = 0.0;
+ double Tab[6][NTURN], fx[nterm2], fy[nterm2], fx2[nterm2], fy2[nterm2];
+ double x = 0.0, xp = 0.0, y = 0.0, zp = 0.0;
double x0 = 1e-6, xp0 = 0.0, z0 = 1e-6, zp0 = 0.0;
double xstep = 0.0, zstep = 0.0;
int nb_freq[2] = {0, 0};
@@ -4447,7 +4499,7 @@ void spectrum(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
if (diffusion) nturn = 2*Nbtour;
-// if (trace) printf("Entering fmap ... results in %s\n\n",fic);
+// if (trace) fprintf(stdout, "Entering fmap ... results in %s\n\n",fic);
/* Opening file */
if ((xoutf = fopen(xfic, "w")) == NULL) {
@@ -4463,36 +4515,36 @@ void spectrum(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
fprintf(xoutf,"# TRACY II v. 2.6 -- %s -- %s \n", xfic, asctime2(newtime));
fprintf(zoutf,"# TRACY II v. 2.6 -- %s -- %s \n", zfic, asctime2(newtime));
// fprintf(outf,"# nu = f(x) \n");
-// fprintf(outf,"# x[m] z[m] fx fz dfx dfz\n");
+// fprintf(outf,"# x[m] z[m] fx fy dfx dfy\n");
- if ((Nbx <= 1) || (Nbz <= 1))
- fprintf(stdout,"fmap: Error Nbx=%ld Nbz=%ld\n",Nbx,Nbz);
+ if ((Nbx <= 1) || (Nby <= 1))
+ fprintf(stdout,"fmap: Error Nbx=%ld Nby=%ld\n",Nbx,Nby);
xp = xp0;
zp = zp0;
xstep = xmax/sqrt((double)Nbx);
- zstep = zmax/sqrt((double)Nbz);
+ zstep = ymax/sqrt((double)Nby);
for (i = 0; i <= Nbx; i++) {
x = x0 + sqrt((double)i)*xstep;
- for (j = 0; j<= Nbz; j++) {
- z = z0 + sqrt((double)j)*zstep;
- Trac_Simple4DCOD(x,xp,z,zp,energy,0.0,nturn,Tab,&status);
+ for (j = 0; j<= Nby; j++) {
+ y = z0 + sqrt((double)j)*zstep;
+ Trac_Simple4DCOD(x,xp,y,zp,energy,0.0,nturn,Tab,&status);
if (status) {
- Get_NAFF(nterm2, Nbtour, Tab, fx, fz, nb_freq);
+ Get_NAFF(nterm2, Nbtour, Tab, fx, fy, nb_freq);
}
else {
- fx[0] = 0.0; fz[0] = 0.0;
- fx2[0] = 0.0; fz2[0] = 0.0;
+ fx[0] = 0.0; fy[0] = 0.0;
+ fx2[0] = 0.0; fy2[0] = 0.0;
}
// printout value
if (!diffusion){
- fprintf(xoutf,"%14.6e %14.6e", x, z);
- fprintf(zoutf,"%14.6e %14.6e", x, z);
- fprintf(stdout,"%14.6e %14.6e", x, z);
+ fprintf(xoutf,"%14.6e %14.6e", x, y);
+ fprintf(zoutf,"%14.6e %14.6e", x, y);
+ fprintf(stdout,"%14.6e %14.6e", x, y);
for (k = 0; k < nb_freq[0]; k++){
fprintf(xoutf," %14.6e", fx[k]);
@@ -4500,8 +4552,8 @@ void spectrum(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
}
for (k = 0; k < nb_freq[1]; k++){
- fprintf(zoutf," %14.6e", fz[k]);
- fprintf(stdout," %14.6e", fz[k]);
+ fprintf(zoutf," %14.6e", fy[k]);
+ fprintf(stdout," %14.6e", fy[k]);
}
fprintf(stdout,"\n");
@@ -4510,9 +4562,9 @@ void spectrum(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
}
// else {
// fprintf(outf,"%14.6e %14.6e %14.6e %14.6e %14.6e %14.6e\n",
-// x, z, fx[0], fz[0], fx[0]-fx2[0], fz[0]-fz2[0]);
+// x, y, fx[0], fy[0], fx[0]-fx2[0], fy[0]-fy2[0]);
// fprintf(stdout,"%14.6e %14.6e %14.6e %14.6e %14.6e %14.6e\n",
-// x, z, fx[0], fz[0], fx[0]-fx2[0], fz[0]-fz2[0]);
+// x, y, fx[0], fy[0], fx[0]-fx2[0], fy[0]-fy2[0]);
// }
}
}
@@ -4571,14 +4623,13 @@ void TracCO(double x, double px, double y, double py, double dp, double ctau,
getcod(dp, lastpos);
getelem(pos-1,&Cell);
- if (!trace) printf("dp= % .5e %% xcod= % .5e mm zcod= % .5e mm \n",
- dp*1e2, Cell.BeamPos[0]*1e3, Cell.BeamPos[2]*1e3);
+ if (!trace) fprintf(stdout, "dp= % .5e %% xcod= % .5e mm zcod= % .5e mm \n",
+ dp*1e2, Cell.BeamPos[x_]*1e3, Cell.BeamPos[y_]*1e3);
/* Tracking coordinates around the closed orbit */
- x1[0] = x + Cell.BeamPos[0]; x1[1] = px + Cell.BeamPos[1];
- x1[2] = y + Cell.BeamPos[2]; x1[3] = py + Cell.BeamPos[3];
- x1[4] = dp; x1[5] = ctau; // line true in 4D tracking
-// x1[4] = dp + Cell.BeamPos[4]; x1[5] = ctau + Cell.BeamPos[5];
+ x1[x_] = x + Cell.BeamPos[x_]; x1[px_] = px + Cell.BeamPos[px_];
+ x1[y_] = y + Cell.BeamPos[y_]; x1[py_] = py + Cell.BeamPos[py_];
+ x1[delta_] = dp; x1[ct_] = ctau; // line true in 4D tracking
lastn = 0;
lostF = true;
@@ -4593,7 +4644,7 @@ void TracCO(double x, double px, double y, double py, double dp, double ctau,
if (!trace) { // print initial conditions
fprintf(outf1, "%6ld %+10.5e %+10.5e %+10.5e %+10.5e"
" %+10.5e %+10.5e \n",
- lastn, x1[0], x1[1], x1[2], x1[3], x1[4], x1[5]);
+ lastn, x1[x_], x1[px_], x1[y_], x1[py_], x1[delta_], x1[ct_]);
}
// Cell_Pass(pos-1L, globval.Cell_nLoc, x1, lastpos);
@@ -4604,10 +4655,10 @@ void TracCO(double x, double px, double y, double py, double dp, double ctau,
if (lastpos != pos-1L)
{
- printf("TracCO: Particle lost \n");
+ 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[0], x1[1], x1[2], x1[3], x1[4], x1[5]);
+ lastn, x1[x_], x1[px_], x1[y_], x1[py_], x1[delta_], x1[ct_]);
}
}
@@ -4655,7 +4706,7 @@ void getA4antidamping()
{
qlist[nquad] = i;
nquad++;
- if (!trace) printf("%s % f\n",Cell.Elem.PName, Cell.Elem.M->PBpar[2L + HOMmax]);
+ if (!trace) fprintf(stdout, "%s % f\n",Cell.Elem.PName, Cell.Elem.M->PBpar[2L + HOMmax]);
}
}
}
@@ -4675,11 +4726,11 @@ void getA4antidamping()
/****************************************************************************/
-/* void fmapfull(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
+/* void fmapfull(long Nbx, long Nby, long Nbtour, double xmax, double ymax,
double energy, bool *status)
Purpose:
- Compute a frequency map of Nbx x Nbz points
+ Compute a frequency map of Nbx x Nby points
For each set of initial conditions the particle is tracked over
Nbtour for an energy offset dp
@@ -4691,7 +4742,7 @@ void getA4antidamping()
Nbx horizontal step number
Nby vertical step number
xmax horizontal maximum amplitude
- zmax vertical maximum amplitude
+ ymax vertical maximum amplitude
Nbtour number of turn for tracking
energy particle energy offset
@@ -4713,19 +4764,19 @@ void getA4antidamping()
****************************************************************************/
#define NTERM 10
-void fmapfull(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
+void fmapfull(long Nbx, long Nby, long Nbtour, double xmax, double ymax,
double energy, bool diffusion)
{
FILE * outf;
const char fic[] = "fmapfull.out";
int i, j, k;
double Tab[DIM][NTURN], Tab0[DIM][NTURN];
- double fx[NTERM], fz[NTERM], fx2[NTERM], fz2[NTERM];
- double x = 0.0, xp = 0.0, z = 0.0, zp = 0.0;
+ double fx[NTERM], fy[NTERM], fx2[NTERM], fy2[NTERM];
+ double x = 0.0, xp = 0.0, y = 0.0, zp = 0.0;
double x0 = 1e-6, xp0 = 0.0, z0 = 1e-6, zp0 = 0.0;
double xstep = 0.0, zstep = 0.0;
int nb_freq[2] = {0, 0};
- double nux1[NTERM], nuz1[NTERM],nux2[NTERM], nuz2[NTERM];
+ double nux1[NTERM], nuy1[NTERM],nux2[NTERM], nuy2[NTERM];
long nturn = Nbtour;
bool status=true;
struct tm *newtime;
@@ -4738,7 +4789,7 @@ void fmapfull(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
if (diffusion) nturn = 2*Nbtour;
- if (trace) printf("Entering fmap ... results in %s\n\n",fic);
+ if (trace) fprintf(stdout, "Entering fmap ... results in %s\n\n",fic);
/* Opening file */
if ((outf = fopen(fic, "w")) == NULL) {
@@ -4755,7 +4806,7 @@ void fmapfull(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
fprintf(outf,"%s",name);
}
for (k = 0; k < NTERM; k++){
- sprintf(name,"f%dz ",k);
+ sprintf(name,"f%dy ",k);
fprintf(outf,"%s",name);
}
@@ -4768,98 +4819,98 @@ void fmapfull(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
fprintf(outf,"%s",name);
}
for (k = 0; k < NTERM; k++){
- sprintf(name,"df%dz ",k);
+ sprintf(name,"df%dy ",k);
fprintf(outf,"%s",name);
}
fprintf(outf,"\n");
}
- if ((Nbx <= 1) || (Nbz <= 1))
- fprintf(stdout,"fmap: Error Nbx=%ld Nbz=%ld\n",Nbx,Nbz);
+ if ((Nbx <= 1) || (Nby <= 1))
+ fprintf(stdout,"fmap: Error Nbx=%ld Nby=%ld\n",Nbx,Nby);
xp = xp0;
zp = zp0;
xstep = xmax/sqrt((double)Nbx);
- zstep = zmax/sqrt((double)Nbz);
+ zstep = ymax/sqrt((double)Nby);
for (i = 0; i <= Nbx; i++) {
x = x0 + sqrt((double)i)*xstep;
- for (j = 0; j<= Nbz; j++) {
- z = z0 + sqrt((double)j)*zstep;
- Trac_Simple4DCOD(x,xp,z,zp,energy,0.0,nturn,Tab,&status);
+ for (j = 0; j<= Nby; j++) {
+ y = z0 + sqrt((double)j)*zstep;
+ Trac_Simple4DCOD(x,xp,y,zp,energy,0.0,nturn,Tab,&status);
if (status) {
- Get_NAFF(NTERM, Nbtour, Tab, fx, fz, nb_freq);
+ Get_NAFF(NTERM, Nbtour, Tab, fx, fy, nb_freq);
for (k = 0; k < nb_freq[0]; k++){
nux1[k] = fx[k];
}
for (k = 0; k < nb_freq[1]; k++){
- nuz1[k] = fz[k];
+ nuy1[k] = fy[k];
}
for (k = nb_freq[0]; k < NTERM; k++){
nux1[k] = 0.0;
}
for (k = nb_freq[1]; k < NTERM; k++){
- nuz1[k] = 0.0;
+ nuy1[k] = 0.0;
}
if (diffusion){
Get_Tabshift(Tab,Tab0,Nbtour,Nbtour); // shift data for second round NAFF
- Get_NAFF(NTERM, Nbtour, Tab0, fx2, fz2, nb_freq); // gets frequency vectors
+ Get_NAFF(NTERM, Nbtour, Tab0, fx2, fy2, nb_freq); // gets frequency vectors
for (k = 0; k < nb_freq[0]; k++){
nux2[k] = fx2[k];
}
for (k = 0; k < nb_freq[1]; k++){
- nuz2[k] = fz2[k];
+ nuy2[k] = fy2[k];
}
for (k = nb_freq[0]; k < NTERM; k++){
nux2[k] = 0.0;
}
for (k = nb_freq[1]; k < NTERM; k++){
- nuz2[k] = 0.0;
+ nuy2[k] = 0.0;
}
}
}
else {
for (k = 0; k < NTERM; k++){
nux1[k] = 0.0;
- nuz1[k] = 0.0;
+ nuy1[k] = 0.0;
nux2[k] = 0.0;
- nuz2[k] = 0.0;
+ nuy2[k] = 0.0;
}
}
// printout value
if (!diffusion){
- fprintf(outf,"%14.6e %14.6e ", x, z);
+ fprintf(outf,"%14.6e %14.6e ", x, y);
for (k = 0; k < NTERM; k++){
fprintf(outf,"%14.6e ", nux1[k]);
}
for (k = 0; k < NTERM; k++){
- fprintf(outf,"%14.6e ", nuz1[k]);
+ fprintf(outf,"%14.6e ", nuy1[k]);
}
fprintf(outf,"\n");
-// fprintf(stdout,"%14.6e %14.6e %14.6e %14.6e\n", x, z, nux1, nuz1);
+// fprintf(stdout,"%14.6e %14.6e %14.6e %14.6e\n", x, y, nux1, nuy1);
}
else {
- fprintf(outf,"%14.6e %14.6e ", x, z);
+ fprintf(outf,"%14.6e %14.6e ", x, y);
for (k = 0; k < NTERM; k++){
fprintf(outf,"%14.6e ", nux1[k]);
}
for (k = 0; k < NTERM; k++){
- fprintf(outf,"%14.6e ", nuz1[k]);
+ fprintf(outf,"%14.6e ", nuy1[k]);
}
for (k = 0; k < NTERM; k++){
fprintf(outf,"%14.6e ", nux2[k]);
}
for (k = 0; k < NTERM; k++){
- fprintf(outf,"%14.6e ", nuz2[k]);
+ fprintf(outf,"%14.6e ", nuy2[k]);
}
fprintf(outf,"\n");
// fprintf(stdout,"%14.6e %14.6e %14.6e %14.6e %14.6e %14.6e\n",
-// x, z, nux1, nuz1, fx[0]-fx2[0], fz[0]-fz2[0]);
+// x, y, nux1, nuy1, fx[0]-fx2[0], fy[0]-fy2[0]);
}
}
}
@@ -4869,11 +4920,11 @@ void fmapfull(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
#undef NTERM
/****************************************************************************/
-/* void Dyna(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
+/* void Dyna(long Nbx, long Nby, long Nbtour, double xmax, double ymax,
double energy, bool *status)
Purpose:
- Compute a frequency map of Nbx x Nbz points
+ Compute a frequency map of Nbx x Nby points
For each set of initial conditions the particle is tracked over
Nbtour for an energy offset dp
@@ -4885,7 +4936,7 @@ void fmapfull(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
Nbx horizontal step number
Nby vertical step number
xmax horizontal maximum amplitude
- zmax vertical maximum maplitude
+ ymax vertical maximum maplitude
Nbtour number of turn for tracking
energy particle energy offset
@@ -4907,14 +4958,14 @@ void fmapfull(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
****************************************************************************/
#define NTERM2 2
-void Dyna(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
+void Dyna(long Nbx, long Nby, long Nbtour, double xmax, double ymax,
double energy, bool diffusion)
{
FILE * outf;
const char fic[] = "dyna.out";
long i, j;
- double Tab[6][NTURN], fx[NTERM2], fz[NTERM2];
- double x = 0.0, xp = 0.0, z = 0.0, zp = 0.0;
+ double Tab[6][NTURN], fx[NTERM2], fy[NTERM2];
+ double x = 0.0, xp = 0.0, y = 0.0, zp = 0.0;
double x0 = 1e-6, xp0 = 0.0, z0 = 1e-6, zp0 = 0.0;
double xstep = 0.0, zstep = 0.0;
int nb_freq[2] = {0, 0};
@@ -4927,7 +4978,7 @@ void Dyna(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
if (diffusion) nturn = 2*Nbtour;
- if (trace) printf("Entering fmap ... results in %s\n\n",fic);
+ if (trace) fprintf(stdout, "Entering fmap ... results in %s\n\n",fic);
/* Opening file moustache */
if ((outf = fopen(fic, "w")) == NULL)
@@ -4938,32 +4989,32 @@ void Dyna(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
fprintf(outf,"# TRACY II v. 2.6 -- %s -- %s \n", fic, asctime2(newtime));
fprintf(outf,"# nu = f(x) \n");
- fprintf(outf,"# x[m] z[m] fx fz \n");
+ fprintf(outf,"# x[m] z[m] fx fy \n");
- if ((Nbx <= 1) || (Nbz <= 1))
- fprintf(stdout,"fmap: Error Nbx=%ld Nbz=%ld\n",Nbx,Nbz);
+ if ((Nbx <= 1) || (Nby <= 1))
+ fprintf(stdout,"fmap: Error Nbx=%ld Nby=%ld\n",Nbx,Nby);
xp = xp0;
zp = zp0;
xstep = xmax/sqrt((double)Nbx);
- zstep = zmax/sqrt((double)Nbz);
+ zstep = ymax/sqrt((double)Nby);
for (i = 0; i <= Nbx; i++) {
x = x0 + sqrt((double)i)*xstep;
- for (j = 0; j<= Nbz; j++) {
- z = z0 + sqrt((double)j)*zstep;
- Trac_Simple4DCOD(x,xp,z,zp,energy,0.0,nturn,Tab,&status);
- if (status) Get_NAFF(NTERM2, Nbtour, Tab, fx, fz, nb_freq);
+ for (j = 0; j<= Nby; j++) {
+ y = z0 + sqrt((double)j)*zstep;
+ Trac_Simple4DCOD(x,xp,y,zp,energy,0.0,nturn,Tab,&status);
+ if (status) Get_NAFF(NTERM2, Nbtour, Tab, fx, fy, nb_freq);
else {
- fx[0] = 0.0; fz[0] = 0.0;
+ fx[0] = 0.0; fy[0] = 0.0;
}
- fprintf(outf,"%14.6e %14.6e %14.6e %14.6e %d\n", x, z, fx[0], fz[0], status);
- fprintf(stdout,"%14.6e %14.6e %14.6e %14.6e %d\n", x, z, fx[0], fz[0], status);
+ fprintf(outf,"%14.6e %14.6e %14.6e %14.6e %d\n", x, y, fx[0], fy[0], status);
+ fprintf(stdout,"%14.6e %14.6e %14.6e %14.6e %d\n", x, y, fx[0], fy[0], status);
if (diffusion) {
- if (status) Get_NAFF(NTERM2, Nbtour, Tab, fx, fz, nb_freq);
- fprintf(outf,"%14.6e %14.6e %14.6e %14.6e %d\n", x, z, fx[0], fz[0], status);
- fprintf(stdout,"%14.6e %14.6e %14.6e %14.6e %d\n", x, z, fx[0], fz[0], status);
+ if (status) Get_NAFF(NTERM2, Nbtour, Tab, fx, fy, nb_freq);
+ fprintf(outf,"%14.6e %14.6e %14.6e %14.6e %d\n", x, y, fx[0], fy[0], status);
+ fprintf(stdout,"%14.6e %14.6e %14.6e %14.6e %d\n", x, y, fx[0], fy[0], status);
}
}
}
@@ -4973,19 +5024,19 @@ void Dyna(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
for (i = 0; i <= Nbx; i++) {
x = x0 - sqrt((double)i)*xstep;
- for (j = 0; j<= Nbz; j++) {
- z = z0 + sqrt((double)j)*zstep;
- Trac_Simple4DCOD(x,xp,z,zp,energy,0.0,nturn,Tab,&status);
- if (status) Get_NAFF(NTERM2, Nbtour, Tab, fx, fz, nb_freq);
+ for (j = 0; j<= Nby; j++) {
+ y = z0 + sqrt((double)j)*zstep;
+ Trac_Simple4DCOD(x,xp,y,zp,energy,0.0,nturn,Tab,&status);
+ if (status) Get_NAFF(NTERM2, Nbtour, Tab, fx, fy, nb_freq);
else {
- fx[0] = 0.0; fz[0] =0.0;
+ fx[0] = 0.0; fy[0] =0.0;
}
- fprintf(outf,"%14.6e %14.6e %14.6e %14.6e %d\n", x, z, fx[0], fz[0], status);
- fprintf(stdout,"%14.6e %14.6e %14.6e %14.6e %d\n", x, z, fx[0], fz[0], status);
+ fprintf(outf,"%14.6e %14.6e %14.6e %14.6e %d\n", x, y, fx[0], fy[0], status);
+ fprintf(stdout,"%14.6e %14.6e %14.6e %14.6e %d\n", x, y, fx[0], fy[0], status);
if (diffusion) {
- if (status) Get_NAFF(NTERM2, Nbtour, Tab, fx, fz, nb_freq);
- fprintf(outf,"%14.6e %14.6e %14.6e %14.6e\n", x, z, fx[0], fz[0]);
- fprintf(stdout,"%14.6e %14.6e %14.6e %14.6e\n", x, z, fx[0], fz[0]);
+ if (status) Get_NAFF(NTERM2, Nbtour, Tab, fx, fy, nb_freq);
+ fprintf(outf,"%14.6e %14.6e %14.6e %14.6e\n", x, y, fx[0], fy[0]);
+ fprintf(stdout,"%14.6e %14.6e %14.6e %14.6e\n", x, y, fx[0], fy[0]);
}
}
}
@@ -5042,7 +5093,7 @@ void Phase2(long pos, double x,double px,double y, double py,double energy,
fprintf(outf,"# TRACY II v. 2.6 -- %s -- %s \n", fic, asctime2(newtime));
fprintf(outf,"# Phase Space \n");
fprintf(outf,
- "# num x xp z zp dp ctau\n");
+ "# num x xp y zp dp ctau\n");
Trac(x,px,y,py,energy,ctau, Nbtour,pos, lastn, lastpos, outf);
fclose(outf);
@@ -5070,7 +5121,7 @@ void Phase3(long pos, double x,double px,double y, double py,double energy,
fprintf(outf,"# TRACY II v. 2.6 -- %s -- %s \n", fic, asctime2(newtime));
fprintf(outf,"# Phase Space \n");
fprintf(outf,
- "# num x xp z zp dp ctau\n");
+ "# num x xp y zp dp ctau\n");
x1[0] = x; x1[1] = px; x1[2] = y;
x1[3] = py; x1[4] = energy; x1[5] = ctau;
@@ -5119,7 +5170,7 @@ void Phase3(long pos, double x,double px,double y, double py,double energy,
2 2 2 2 T
Given V = (<u >, <uu'>, <u' >, <v >, <uu'>, <v' >)
2 2 2 2 T
- and X = (<x >, <xx'>, <x' >, <z >, <zz'>, <z' >)
+ and X = (<x >, <xx'>, <x' >, <y >, <zz'>, <z' >)
Then X = U2T V where U2T if constucted using the uncoupling R matrix
@@ -5509,7 +5560,7 @@ void Coupling_Edwards_Teng(void)
if (!trace)
{
- fprintf(stdout,"Projected emittances: Ex = % 10.6e, Ez = % 10.6e, Ez/Ex = %10.6e\n",
+ fprintf(stdout,"Projected emittances: Ex = % 10.6e, Ey = % 10.6e, Ez/Ex = %10.6e\n",
W1, W2, W2/W1);
fprintf(stdout,"**************************************\n");
}
@@ -5601,7 +5652,7 @@ void PhaseLongitudinalHamiltonien(void)
exit_(1);
}
- printf("Last stable orbit %f\n", acos(1.0-T*E/eVRF*Hsynchrotron(0.0,-0.098)));
+ fprintf(stdout, "Last stable orbit %f\n", acos(1.0-T*E/eVRF*Hsynchrotron(0.0,-0.098)));
fprintf(outf,"# TRACY II v. 2.6 -- %s \n", asctime2(newtime));
fprintf(outf,"# i ctau dp DH/H H \n#\n");
@@ -5839,7 +5890,7 @@ void Enveloppe2(double x, double px, double y, double py, double dp, double ntur
// Envzm[i] = Cell.BeamPos[2]; Envzp[i] = Cell.BeamPos[2];
// }
- printf("xcod=%.5e mm zcod=% .5e mm \n",
+ fprintf(stdout, "xcod=%.5e mm zcod=% .5e mm \n",
globval.CODvect[0]*1e3, globval.CODvect[2]*1e3);
if ((outf = fopen(fic, "w")) == NULL) {
@@ -5860,7 +5911,7 @@ void Enveloppe2(double x, double px, double y, double py, double dp, double ntur
Cell_Pass(i,i+1, x1, lastpos);
if (lastpos != i+1)
{
- printf("Unstable motion ...\n"); exit_(1);
+ fprintf(stdout, "Unstable motion ...\n"); exit_(1);
}
Envxp[i] = x1[0]; Envxm[i] = x1[0]; Envzp[i] = x1[2]; Envzm[i] = x1[2];
@@ -5874,7 +5925,7 @@ void Enveloppe2(double x, double px, double y, double py, double dp, double ntur
Cell_Pass(i, i+1, x1, lastpos);
if (lastpos != i+1)
{
- printf("Unstable motion ...\n"); exit_(1);
+ fprintf(stdout, "Unstable motion ...\n"); exit_(1);
}
if (x1[0] >= Envxp[i]) Envxp[i] = x1[0];
if (x1[0] <= Envxm[i]) Envxm[i] = x1[0];
@@ -6018,7 +6069,7 @@ void ReadFieldErr(const char *FieldErrorFile)
inf = file_read(FieldErrorFile);
- printf("\n");
+ fprintf(stdout, "\n");
/* read lines*/
while ((line = fgets(in, max_str, inf))) {
@@ -6041,18 +6092,18 @@ void ReadFieldErr(const char *FieldErrorFile)
if (strcmp("seed", name) == 0) { // the line to set random seed
sscanf(line, "%*s %d", &seed_val);
- printf("ReadFieldErr: setting random seed to %d\n", seed_val);
+ fprintf(stdout, "ReadFieldErr: setting random seed to %d\n", seed_val);
set_rnd = true;
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) printf("\nsetting <%s> multipole error to: %-5s r0 = %7.1le\n",keywrd,name,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) {
- printf("ReadFieldErr: seed not defined\n");
+ fprintf(stdout, "ReadFieldErr: seed not defined\n");
exit(1);
}
@@ -6073,7 +6124,7 @@ void ReadFieldErr(const char *FieldErrorFile)
sscanf(prm, "%lf", &An);
if (prt)
- printf(" n = %2d, Bn = %9.1e, An = %9.1e\n", n, Bn, An);
+ fprintf(stdout, " n = %2d, Bn = %9.1e, An = %9.1e\n", n, Bn, An);
/* set multipole errors to horizontal correctors of soleil ring*/
@@ -6094,7 +6145,7 @@ void ReadFieldErr(const char *FieldErrorFile)
//end of the line
}else
continue;
- // printf("%s", line);
+ // fprintf(stdout, "%s", line);
}
fclose(inf);
@@ -6142,7 +6193,7 @@ void AddFieldErrors(const char *name,const char *keywrd, const double r0,
int Fnum = 0;
if (strcmp("all", name) == 0) {
- printf("all: not yet implemented\n");
+ fprintf(stdout, "all: not yet implemented\n");
} else if (strcmp("dip", name) == 0) {
AddFieldValues_type(Dip,keywrd, r0, n, Bn, An);
} else if (strcmp("quad", name) == 0) {
@@ -6154,7 +6205,7 @@ void AddFieldErrors(const char *name,const char *keywrd, const double r0,
if(Fnum > 0)
AddFieldValues_fam(Fnum,keywrd, r0, n, Bn, An);
else
- printf("SetFieldErrors: undefined element %s\n", name);
+ fprintf(stdout, "SetFieldErrors: undefined element %s\n", name);
}
}
@@ -6393,7 +6444,7 @@ void Add_bnL_sys_elem(const int Fnum, const int Knum, const char *keywrd,
Mpole_SetPB(Fnum, Knum, -n); //set for An component
if (prt)
- printf("add the %s error: n=%d component of %s, bnL = %e, %e, anL = %e, %e\n",
+ 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]);
}
@@ -6690,7 +6741,7 @@ void PrintTrack(const char *TrackFile, double x, double px, double y,double py,
xf[i] = x2[i];
if (prt) {
- printf("%4ld %4ld %s %6.4f %10.4f %10.4f %10.4f %10.4f %10.4f %10.4f %4ld \n",
+ 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);
}
@@ -6736,7 +6787,7 @@ void PrintTrack(const char *TrackFile, double x, double px, double y,double py,
conv = 93.88e-4; /*conversion des A en T*/
nqtcorr = GetnKid(ElemIndex("SQ"));
- printf("Number of virtual skew quadrupoles: %d\n",nqtcorr);
+ fprintf(stdout, "Number of virtual skew quadrupoles: %d\n",nqtcorr);
/* open virtual QT corrector file */
if ((fi = fopen(VirtualSkewQuadFile,"r")) == NULL)
@@ -6751,7 +6802,7 @@ void PrintTrack(const char *TrackFile, double x, double px, double y,double py,
//corr_strength = 20.*conv/brho;
SetKLpar(ElemIndex("SQ"),i, mOrder=2L, corr_strength);
if(trace)
- printf("virtual skew quadrupole: %d, qtcorr=%le, corr_strength=%le\n",i,qtcorr[i-1],corr_strength);
+ fprintf(stdout, "virtual skew quadrupole: %d, qtcorr=%le, corr_strength=%le\n",i,qtcorr[i-1],corr_strength);
}
fclose(fi);
}
@@ -6791,8 +6842,8 @@ void CODCorrect(const char *hcorr_file, const char *vcorr_file,int n_orbit,int n
fprintf(stdout, "V-plane %d singular values:\n\n",nwv);
// compute beam response matrix
- printf("\n");
- printf("Computing beam response matrix\n");
+ fprintf(stdout, "\n");
+ fprintf(stdout, "Computing beam response matrix\n");
//get the response matrix between bpm and correctors
gcmats(globval.bpm, globval.hcorr, 1, hcorrIdx);
gcmats(globval.bpm, globval.vcorr, 2, vcorrIdx);
@@ -6838,13 +6889,13 @@ void CODCorrect(const char *hcorr_file, const char *vcorr_file,int n_orbit,int n
cod = CorrectCODs(hOrbitFile, vOrbitFile, OrbScanFile, n_orbit, nwh, nwv, hcorrIdx, vcorrIdx);
/* cod = CorrectCOD_N(ae_file, n_orbit, n_scale, k);*/
- printf("\n");
+ fprintf(stdout, "\n");
if (cod){
- /* printf("done with orbit correction, now do coupling",
+ /* fprintf(stdout, "done with orbit correction, now do coupling",
" correction plus vert. disp\n");*/
- printf("Orbit correction succeeded\n");
+ fprintf(stdout, "Orbit correction succeeded\n");
}else{
fprintf(stdout, "!!! Orbit correction failed\n");
exit_(1);
--
GitLab