diff --git a/tracy/tools/soltracy.cc b/tracy/tools/soltracy.cc
index a78f3a3b86e78653ae417e0582061d143a97e241..138d44f563c0f779aef3d54e61ed337a8dd86204 100644
--- a/tracy/tools/soltracy.cc
+++ b/tracy/tools/soltracy.cc
@@ -1,14 +1,13 @@
-/*
+/************************************
Current revision $Revision$
On branch $Name$
Latest change $Date$ by $Author$
-*/
-#define ORDER 1
+*************************************/
+#define ORDER 1
+//#define ORDER 4
int no_tps = ORDER; // arbitrary TPSA order is defined locally
-
-
extern bool freq_map;
#if HAVE_CONFIG_H
@@ -514,7 +513,8 @@ int main(int argc, char *argv[]) {
UserCommandFlag[i]._FmapFlag_xmax,
UserCommandFlag[i]._FmapFlag_ymax,
UserCommandFlag[i]._FmapFlag_delta,
- UserCommandFlag[i]._FmapFlag_diffusion);
+ UserCommandFlag[i]._FmapFlag_diffusion,
+ UserCommandFlag[i]._FmapFlag_printloss);
#endif
}
@@ -594,8 +594,9 @@ int main(int argc, char *argv[]) {
UserCommandFlag[i]._FmapdpFlag_xmax,
UserCommandFlag[i]._FmapdpFlag_emax,
UserCommandFlag[i]._FmapdpFlag_z,
- UserCommandFlag[i]._FmapdpFlag_diffusion);
- #endif
+ UserCommandFlag[i]._FmapdpFlag_diffusion,
+ UserCommandFlag[i]._FmapdpFlag_printloss);
+ #endif
}
// // if (CodeComparaisonFlag) {
@@ -831,7 +832,7 @@ int main(int argc, char *argv[]) {
UserCommandFlag[i]._Phase_delta,
UserCommandFlag[i]._Phase_ctau,
UserCommandFlag[i]._Phase_nturn);
- printf("the simulation time for phase space in tracy 3 is \n");
+ printf("6D phase space tracking: \n the simulation time for phase space in tracy 3 is \n");
stop = stampstop(start);
/* restore the initial values*/
@@ -954,7 +955,7 @@ Based on parts of functions get_param( ) & ID_corr(), etc in nsls-ii_lib.cc
}
//For debug
- if(!trace){
+ if(trace){
cout << "scl_dbetax = " << scl_dbetax << " scl_dbetay = " << scl_dbetay <<endl;
cout << "scl_dnux = " << scl_dnux << " scl_dnuy = " << scl_dnuy << endl;
cout << "scl_nux = " << scl_nux << " scl_nuy = " << scl_nuy << endl;
diff --git a/tracy/tracy/inc/naffutils.h b/tracy/tracy/inc/naffutils.h
index b7ade74289374fab46a56c9e85076ccc87ce9538..25f0e35a4d37415430833275bb70123a8583dbd0 100644
--- a/tracy/tracy/inc/naffutils.h
+++ b/tracy/tracy/inc/naffutils.h
@@ -24,7 +24,10 @@ void GetChromTrac(long Nb, long Nbtour, double emax, double *xix, double *xiz);
void GetTuneTrac(long Nbtour, double emax, double *nux, double *nuz);
void Trac_Simple4DCOD(double x, double px, double y, double py, double dp,
- double ctau, long nmax, double Tx[][NTURN], bool *status);
+ 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);
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/read_script.h b/tracy/tracy/inc/read_script.h
index ffdb8dcdfc8d733c66bc03db33f8f22c4485b01c..90c049667ed6b6e19285c0a824fd06015a1c7829 100644
--- a/tracy/tracy/inc/read_script.h
+++ b/tracy/tracy/inc/read_script.h
@@ -57,13 +57,16 @@
char _FmapFlag_fmap_file[max_str];
long _FmapFlag_nxpoint, _FmapFlag_nypoint, _FmapFlag_nturn;
double _FmapFlag_xmax, _FmapFlag_ymax, _FmapFlag_delta;
- bool _FmapFlag_diffusion;
+ bool _FmapFlag_diffusion;
+ bool _FmapFlag_printloss;
+
//extern bool FmapdpFlag;
char _FmapdpFlag_fmapdp_file[max_str];
long _FmapdpFlag_nxpoint, _FmapdpFlag_nepoint, _FmapdpFlag_nturn;
double _FmapdpFlag_xmax, _FmapdpFlag_emax, _FmapdpFlag_z;
- bool _FmapdpFlag_diffusion;
+ bool _FmapdpFlag_diffusion;
+ bool _FmapdpFlag_printloss;
//MomentumAccFlag;
@@ -109,14 +112,14 @@
strcpy(_FmapFlag_fmap_file,"fmap.out");
_FmapFlag_nxpoint=31L, _FmapFlag_nypoint=21L, _FmapFlag_nturn=516L;
_FmapFlag_xmax=0.025, _FmapFlag_ymax=0.005, _FmapFlag_delta=0.0;
- _FmapFlag_diffusion = true;
+ _FmapFlag_diffusion = true, _FmapFlag_printloss = false;
/*fmap for off momentum particle*/
strcpy(_FmapdpFlag_fmapdp_file,"fmapdp.out");
_FmapdpFlag_nxpoint=31L, _FmapdpFlag_nepoint=21L, _FmapdpFlag_nturn=516L;
_FmapdpFlag_xmax=0.025, _FmapdpFlag_emax=0.005, _FmapdpFlag_z=0.0;
- _FmapdpFlag_diffusion = true;
-
+ _FmapdpFlag_diffusion = true, _FmapdpFlag_printloss = false;
+
/* tune shift with amplitude*/
strcpy(_AmplitudeTuneShift_nudx_file,"nudx.out");
strcpy(_AmplitudeTuneShift_nudz_file,"nudz.out");
@@ -161,7 +164,9 @@ strcpy(_EnergyTuneShift_nudp_file,"nudp.out");
};
/***** file names************/
- //files with multipole errors; for soleil lattice
+extern char lat_file[max_str];
+
+//files with multipole errors; for soleil lattice
extern char fic_hcorr[max_str],fic_vcorr[max_str], fic_skew[max_str];
extern char multipole_file[max_str];
//file of source of coupling; for soleil lattice
diff --git a/tracy/tracy/inc/soleillib.h b/tracy/tracy/inc/soleillib.h
index c34ed541f7123548276edfd467e40976d3adf6d2..dab545916de16239f06322e24312359609df9212 100644
--- a/tracy/tracy/inc/soleillib.h
+++ b/tracy/tracy/inc/soleillib.h
@@ -57,13 +57,15 @@ void TunesShiftWithEnergy(const char *NudpFile,long Nb, long Nbtour, double emax
//void fmapdp(long Nbx, long Nbe, long Nbtour, double xmax, double emax,
// double z, bool diffusion, bool matlab);
void fmap(const char *FmapFile, long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
- double energy, bool diffusion);
+ double energy, bool diffusion);
+void fmap(const char *FmapFile, long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
+ double energy, bool diffusion, bool printloss);
void fmap_p(const char *FmapFile, long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
- double energy, bool diffusion, int numprocs, int myid);
+ double energy, bool diffusion, bool printloss, int numprocs, int myid);
void fmapdp(const char *FmapdpFile, long Nbx, long Nbe, long Nbtour, double xmax, double emax,
- double z, bool diffusion);
+ double z, bool diffusion, bool printloss);
void fmapdp_p(const char *FmapdpFile, long Nbx, long Nbe, long Nbtour, double xmax, double emax,
- double z, bool diffusion, int numprocs, int myid);
+ double z, bool diffusion, bool printloss, int numprocs, int myid);
void Nu_Naff(void);
diff --git a/tracy/tracy/src/naffutils.cc b/tracy/tracy/src/naffutils.cc
index d0d92979e0cf67cf2390e5619cfdf150d83a93ae..dfdac96d7eb6725a9de94eb0ad7c1dfccc6c3d5c 100644
--- a/tracy/tracy/src/naffutils.cc
+++ b/tracy/tracy/src/naffutils.cc
@@ -8,9 +8,9 @@
*/
/*
- Current revision $Revision: 1.3 $
+ Current revision $Revision: 1.4 $
On branch $Name: not supported by cvs2svn $
- Latest change $Date: 2010-12-01 15:29:26 $ by $Author: nadolski $
+ Latest change $Date: 2013-06-22 15:17:21 $ by $Author: nadolski $
*/
#include "complexeheader.h"
@@ -50,8 +50,8 @@ double pi = M_PI;
useful for connection with NAFF
19/01/03 tracking around the closed orbit
****************************************************************************/
-void Trac_Simple4DCOD(double x, double px, double y, double py, double dp,
- double ctau, long nmax, double Tx[][NTURN], bool *status2)
+ void Trac_Simple4DCOD(double x, double px, double y, double py, double dp,
+ double ctau, long nmax, double Tx[][NTURN], bool *status2)
{
bool lostF = false; /* Lost particle Flag */
ss_vect<double> x1; /* Tracking coordinates */
@@ -109,9 +109,125 @@ void Trac_Simple4DCOD(double x, double px, double y, double py, double dp,
if (lastpos != globval.Cell_nLoc)
{ /* Particle lost: Error message section */
*status2 = false;
- printf("Trac_Simple: Particle lost \n");
- fprintf(stdout, "turn=%5ld plane= %1d %+10.5g %+10.5g %+10.5g %+10.5g %+10.5g %+10.5g \n", lastn-1,
+ if (trace){
+ fprintf(stdout, "Trac_Simple: 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]);
+ }
+ }
+}
+
+/****************************************************************************/
+/* void Trac_Simple4DCOD(double x, double px, double y, double py, double dp,
+ double ctau, long nmax, double Tx[][NTURN], bool *status2)
+
+ Purpose:
+ Single particle tracking around the closed orbit for NTURN turns
+ The 6D phase trajectory is saved in a array, but the
+ tracked dp and ctau is not about the COD.
+
+ Input:
+ x, px, y, py 4 transverse coordinates
+ 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
+ Tx 6xNTURN matrix of phase trajectory
+
+ Return:
+ none
+
+ Global variables:
+ NTURN number of turn for tracking
+ globval
+
+ Specific functions:
+ Cell_Pass
+
+ Comments:
+ useful for connection with NAFF
+ 19/01/03 tracking around the closed orbit
+
+ 11/06/2013 Modified by Jianfeng Zhang @ LAL
+ The same feature as function
+ Trac_Simple4DCOD(double x, double px, double y, double py, double dp,
+ double ctau, long nmax, double Tx[][NTURN], bool *status2)
+ but add the feature to extract the position of the location of
+ the lost particle;
+ called by function
+ fmap(const char *FmapFile, long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
+ double energy, bool diffusion, bool loss)
+ in soleillib.cc.
+
+****************************************************************************/
+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)
+{
+ bool lostF = false; /* Lost particle Flag */
+ Vector2 aperture = {1.0, 1.0};
+
+ lastn = 0;
+ lastpos = globval.Cell_nLoc;
+ *status2 = true; /* stable */
+ x1[0]=0,x1[1]=0,x1[2]=0,x1[3]=0,x1[4]=0,x1[5]=0;
+
+ if (globval.MatMeth) Cell_Concat(dp);
+
+ /* Get closed orbit */
+
+ getcod(dp, lastpos);
+
+ if (status.codflag)
+ if (trace) fprintf(stdout,"Track_Simple4DCOD: dp= % .5e %% xcod= % .5e mm zcod= % .5e mm \n",
+ dp*1e2, globval.CODvect[0]*1e3, globval.CODvect[2]*1e3);
+
+ /* Tracking coordinates around the closed orbit */
+ x1[0] = x + globval.CODvect[0]; x1[1] = px + globval.CODvect[1];
+ x1[2] = y + globval.CODvect[2]; x1[3] = py + globval.CODvect[3];
+ x1[4] = dp; x1[5] = ctau;
+
+ Tx[0][lastn] = x1[0]; Tx[1][lastn] = x1[1];
+ Tx[2][lastn] = x1[2]; Tx[3][lastn] = x1[3];
+ Tx[4][lastn] = x1[4]; Tx[5][lastn] = x1[5];
+ lastn++;
+
+ do
+ { /* tracking through the ring */
+ if ((lastpos == globval.Cell_nLoc) &&
+ (fabs(x1[0]) < aperture[0]) && (fabs(x1[2]) < aperture[1]) && status.codflag)
+ {
+ if (globval.MatMeth)
+ Cell_fPass(x1, lastpos);
+ else
+ Cell_Pass(0, globval.Cell_nLoc, x1, lastpos);
+ Tx[0][lastn] = x1[0]; Tx[1][lastn] = x1[1];
+ Tx[2][lastn] = x1[2]; Tx[3][lastn] = x1[3];
+ Tx[4][lastn] = x1[4]; Tx[5][lastn] = x1[5];
+ }
+ else
+ {
+ if (trace) fprintf(stdout, "Trac_Simple4DCOD: Particle lost \n");
+ if (trace) fprintf(stdout, "%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]);
+ lostF = true;
+ *status2 = false;
+ }
+ lastn++;
+ } while ((lastn < nmax) && (lastpos == globval.Cell_nLoc) && (lostF == false));
+
+
+ if (lastpos != globval.Cell_nLoc)
+ { /* Particle lost: Error message section */
+ *status2 = false;
+ if (trace){
+ 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]);
+ }
}
}
@@ -311,7 +427,6 @@ void Get_NAFF(int nterm, long ndata, double Tab[DIM][NTURN],
naf_smoy(g_NAFVariable.ZTABS);
naf_prtabs(g_NAFVariable.KTABS,g_NAFVariable.ZTABS, 20);
-// trace=on;
naf_mftnaf(nterm,fabs(g_NAFVariable.FREFON)/g_NAFVariable.m_dneps);
/* fill up H-frequency vector */
diff --git a/tracy/tracy/src/physlib.cc b/tracy/tracy/src/physlib.cc
index 7a9db3bbc1b85e6eafa8b30a131ef58c51ac27cc..c2a1bf063cc87c99835ae710f15580e4885cb766 100644
--- a/tracy/tracy/src/physlib.cc
+++ b/tracy/tracy/src/physlib.cc
@@ -6,9 +6,9 @@
L. Nadolski SOLEIL 2002 Link to NAFF, Radia field maps
J. Bengtsson NSLS-II, BNL 2004 -
*/
-/* Current revision $Revision: 1.20 $
+/* Current revision $Revision: 1.21 $
On branch $Name: not supported by cvs2svn $
- Latest change by $Author: zhang $
+ Latest change by $Author: nadolski $
*/
/**************************/
@@ -81,7 +81,7 @@ uint32_t stampstart(void) {
// print detailed time in milliseconds
if (timedebug)
- printf("TIMESTAMP-START\t %d:%02d:%02d:%d (~%d ms)\n", tm->tm_hour,
+ fprintf(stdout, "TIMESTAMP-START\t %d:%02d:%02d:%d (~%d ms)\n", tm->tm_hour,
tm->tm_min, tm->tm_sec, tv.tv_usec, tm->tm_hour * 3600 * 1000
+ tm->tm_min * 60 * 1000 + tm->tm_sec * 1000
+ tv.tv_usec / 1000);
@@ -135,11 +135,11 @@ uint32_t stampstop(uint32_t start) {
// print detailed time in milliseconds
if (timedebug) {
- printf("TIMESTAMP-END\t %d:%02d:%02d:%d (~%d ms) \n", tm->tm_hour,
+ fprintf(stdout, "TIMESTAMP-END\t %d:%02d:%02d:%d (~%d ms) \n", tm->tm_hour,
tm->tm_min, tm->tm_sec, tv.tv_usec, tm->tm_hour * 3600 * 1000
+ tm->tm_min * 60 * 1000 + tm->tm_sec * 1000
+ tv.tv_usec / 1000);
- printf("ELAPSED\t %d ms\n", stop - start);
+ fprintf(stdout, "ELAPSED\t %d ms\n", stop - start);
}
uint32_t delta, hour, minute, second, millisecond;
@@ -150,7 +150,7 @@ uint32_t stampstop(uint32_t start) {
millisecond = delta - 3600000 * hour - minute * 60000 - second * 1000;
if (prt)
- printf("ELAPSED\t %d h %d min %d s %d ms\n", hour, minute, second,
+ fprintf(stdout, "ELAPSED\t %d h %d min %d s %d ms\n", hour, minute, second,
millisecond);
return (stop);
@@ -191,57 +191,59 @@ uint32_t stampstop(uint32_t start) {
****************************************************************************/
void printglob(void) {
- printf("\n***************************************************************"
+ fprintf(stdout, "\n***************************************************************"
"***************\n");
- printf("\n");
- printf(" dPcommon = %9.3e dPparticle = %9.3e"
+ fprintf(stdout, "\n");
+ fprintf(stdout, " dPcommon = %9.3e dPparticle = %9.3e"
" Energy [GeV] = %.3f\n", globval.dPcommon, globval.dPparticle,
globval.Energy);
- printf(" MaxAmplx [m] = %9.3e MaxAmply [m] = %9.3e"
+ fprintf(stdout, " MaxAmplx [m] = %9.3e MaxAmply [m] = %9.3e"
" RFAccept [%%] = +/- %4.2f\n", Cell[0].maxampl[X_][1],
Cell[0].maxampl[Y_][1], globval.delta_RF * 1e2);
- printf(" MatMeth = %s ", globval.MatMeth ? "TRUE " : "FALSE");
- printf(" Cavity_On = %s ", globval.Cavity_on ? "TRUE " : "FALSE");
- printf(" Radiation_On = %s \n", globval.radiation ? "TRUE " : "FALSE");
+ fprintf(stdout, " MatMeth = %s ", globval.MatMeth ? "TRUE " : "FALSE");
+ fprintf(stdout, " Cavity_On = %s ", globval.Cavity_on ? "TRUE " : "FALSE");
+ fprintf(stdout, " Radiation_On = %s \n", globval.radiation ? "TRUE " : "FALSE");
if(globval.bpm == 0)
- printf(" bpm = 0");
+ fprintf(stdout, " bpm = %3d", 0);
else
- printf(" bpm = %3d", GetnKid(globval.bpm));
+ fprintf(stdout, " bpm = %3d", GetnKid(globval.bpm));
if(globval.qt == 0)
- printf(" qt = 0 \n");
+ fprintf(stdout, " qt = 0 \n");
else
- printf(" qt = %3d \n", GetnKid(globval.qt));
+ fprintf(stdout, " qt = %3d \n",
+ GetnKid(globval.qt));
if(globval.hcorr == 0)
- printf(" hcorr = 0");
+ fprintf(stdout, " hcorr = %3d", 0);
else
- printf(" hcorr = %3d", GetnKid(globval.hcorr));
+ fprintf(stdout, " hcorr = %3d", GetnKid(globval.hcorr));
if(globval.vcorr == 0)
- printf(" vcorr = 0 \n");
+ fprintf(stdout, " vcorr = 0 \n");
else
- printf(" vcorr = %3d \n", GetnKid(globval.vcorr));
+ fprintf(stdout, " vcorr = %3d \n",
+ GetnKid(globval.vcorr));
- printf(" Chambre_On = %s \n", globval.Aperture_on ? "TRUE " : "FALSE");
+ fprintf(stdout, " Chambre_On = %s \n", globval.Aperture_on ? "TRUE " : "FALSE");
- printf(" alphac = %8.4e\n", globval.Alphac);
- printf(" nux = %8.6f nuy = %8.6f",
+ fprintf(stdout, " alphac = %8.4e\n", globval.Alphac);
+ fprintf(stdout, " nux = %+8.6f nuy = %+8.6f",
globval.TotalTune[X_], globval.TotalTune[Y_]);
if (globval.Cavity_on)
- printf(" omega = %13.9f\n", globval.Omega);
+ fprintf(stdout, " omega = %+13.9f\n", globval.Omega);
else {
- printf("\n");
- printf(" ksix = %8.6f ksiy = %8.6f\n",
+ fprintf(stdout, "\n");
+ fprintf(stdout, " ksix = %+8.6f ksiy = %+8.6f\n",
globval.Chrom[X_], globval.Chrom[Y_]);
}
- printf("\n");
- printf(" OneTurn matrix:\n");
- printf("\n");
+ fprintf(stdout, "\n");
+ fprintf(stdout, " OneTurn matrix:\n");
+ fprintf(stdout, "\n");
prtmat(ss_dim, globval.OneTurnMat);
- printf("\nTwiss parameters at entrance:\n");
+ fprintf(stdout, "\nTwiss parameters at entrance:\n");
printf(
- " Betax [m] = % 9.3e Alphax = % 9.3e Etax [m] = % 9.3e Etaxp = % 9.3e\n"
- " Betay [m] = % 9.3e Alphay = % 9.3e Etay [m] = % 9.3e Etayp = % 9.3e\n\n",
+ " Betax [m] = %+9.3e Alphax = %+9.3e Etax [m] = %+9.3e Etaxp = %+9.3e\n"
+ " Betay [m] = %+9.3e Alphay = %+9.3e Etay [m] = %+9.3e Etayp = %+9.3e\n\n",
Cell[1].Beta[X_], Cell[1].Alpha[X_], Cell[1].Eta[X_],
Cell[1].Etap[X_], Cell[1].Beta[Y_], Cell[1].Alpha[Y_],
Cell[1].Eta[Y_], Cell[1].Etap[Y_]);
@@ -743,9 +745,9 @@ void track(const char *file_name, double ic1, double ic2, double ic3,
lastn = 0;
if (prt) {
- printf("\n");
- printf("track:\n");
- printf("%4ld %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f\n", lastn, 1e3
+ fprintf(stdout, "\n");
+ fprintf(stdout, "track:\n");
+ fprintf(stdout, "%4ld %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f\n", lastn, 1e3
* x2[x_], 1e3 * x2[px_], 1e3 * x2[y_], 1e3 * x2[py_], 1e2
* x2[delta_], 1e3 * x2[ct_]);
}
@@ -939,7 +941,7 @@ void Trac(double x, double px, double y, double py, double dp, double ctau,
Cell_Pass(0L, pos, x1, lastpos);
if (lastpos != pos) {
- printf("Trac: Particle lost \n");
+ 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]);
@@ -994,9 +996,9 @@ bool chk_if_lost(double x0, double y0, double delta, long int nturn, bool floqs)
lastn = 0;
if (prt) {
- printf("\n");
- printf("chk_if_lost:\n");
- printf("%4ld %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f\n", lastn,
+ fprintf(stdout, "\n");
+ fprintf(stdout, "chk_if_lost:\n");
+ fprintf(stdout, "%4ld %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f\n", lastn,
1e3 * x[x_], 1e3 * x[px_], 1e3 * x[y_], 1e3 * x[py_], 1e2
* x[delta_], 1e3 * x[ct_]);
}
@@ -1007,7 +1009,7 @@ bool chk_if_lost(double x0, double y0, double delta, long int nturn, bool floqs)
else
Cell_Pass(0, globval.Cell_nLoc, x, lastpos);
if (prt)
- printf("%4ld %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f\n", lastn, 1e3
+ fprintf(stdout, "%4ld %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f\n", lastn, 1e3
* x[x_], 1e3 * x[px_], 1e3 * x[y_], 1e3 * x[py_], 1e2
* x[delta_], 1e3 * x[ct_]);
} while ((lastn != nturn) && (lastpos == globval.Cell_nLoc));
@@ -1051,7 +1053,7 @@ void getdynap(double &r, double phi, double delta, double eps, int nturn,
const double r_reset = 1e-3, r0 = 10e-3;
if (prt)
- printf("\n");
+ fprintf(stdout, "\n");
if (globval.MatMeth)
Cell_Concat(delta);
@@ -1063,14 +1065,14 @@ void getdynap(double &r, double phi, double delta, double eps, int nturn,
while (rmax - rmin >= eps) {
r = rmin + (rmax - rmin) / 2.0;
if (prt)
- printf("getdynap: %6.3f %6.3f %6.3f\n", 1e3 * rmin, 1e3 * rmax, 1e3
+ fprintf(stdout, "getdynap: %6.3f %6.3f %6.3f\n", 1e3 * rmin, 1e3 * rmax, 1e3
* r);
if (!chk_if_lost(r * cos(phi), r * sin(phi), delta, nturn, floqs))
rmin = r;
else
rmax = r;
if (rmin > rmax) {
- printf("getdynap: rmin > rmax\n");
+ fprintf(stdout, "getdynap: rmin > rmax\n");
exit_(0);
}
@@ -1128,7 +1130,7 @@ void getcsAscr(void) {
}
}
if (!InvMat(6, globval.Ascrinv))
- printf(" *** Ascr is singular\n");
+ fprintf(stdout, " *** Ascr is singular\n");
}
/****************************************************************************
@@ -1184,21 +1186,21 @@ void dynap(FILE *fp, double r, const double delta, const double eps,
if (floqs) {
Ring_GetTwiss(false, delta);
if (print) {
- printf("\n");
- printf("Dynamical Aperture (Floquet space):\n");
- printf(" x^ y^\n");
- printf("\n");
+ fprintf(stdout, "\n");
+ fprintf(stdout, "Dynamical Aperture (Floquet space):\n");
+ fprintf(stdout, " x^ y^\n");
+ fprintf(stdout, "\n");
}
fprintf(fp, "# Dynamical Aperture (Floquet space):\n");
fprintf(fp, "# x^ y^\n");
fprintf(fp, "#\n");
} else {
if (print) {
- printf("\n");
- printf("Dynamical Aperture:\n");
- printf(" x y\n");
- printf(" [mm] [mm]\n");
- printf("\n");
+ fprintf(stdout, "\n");
+ fprintf(stdout, "Dynamical Aperture:\n");
+ fprintf(stdout, " x y\n");
+ fprintf(stdout, " [mm] [mm]\n");
+ fprintf(stdout, "\n");
}
fprintf(fp, "# Dynamical Aperture:\n");
fprintf(fp, "# x y\n");
@@ -1225,19 +1227,19 @@ void dynap(FILE *fp, double r, const double delta, const double eps,
y_max = max(y[i], y_max);
if (!floqs) {
if (print)
- printf(" %6.2f %6.2f\n", 1e3 * x[i], 1e3 * y[i]);
+ fprintf(stdout, " %6.2f %6.2f\n", 1e3 * x[i], 1e3 * y[i]);
fprintf(fp, " %6.2f %6.2f\n", 1e3 * x[i], 1e3 * y[i]);
} else {
if (print)
- printf(" %10.3e %10.3e\n", x[i], y[i]);
+ fprintf(stdout, " %10.3e %10.3e\n", x[i], y[i]);
fprintf(fp, " %10.3e %10.3e\n", x[i], y[i]);
}
fflush(fp);
}
if (print) {
- printf("\n");
- printf(" x^: %6.2f - %5.2f y^: %6.2f - %5.2f mm\n", 1e3 * x_min, 1e3
+ fprintf(stdout, "\n");
+ fprintf(stdout, " x^: %6.2f - %5.2f y^: %6.2f - %5.2f mm\n", 1e3 * x_min, 1e3
* x_max, 1e3 * y_min, 1e3 * y_max);
}
}
@@ -1277,8 +1279,8 @@ double get_aper(int n, double x[], double y[]) {
A += x[n - 1] * y[0] - x[0] * y[n - 1];
// x2 from mid-plane symmetry
A = fabs(A);
- // printf("\n");
- // printf(" Dyn. Aper.: %5.1f mm^2\n", 1e6*A);
+ // fprintf(stdout, "\n");
+ // fprintf(stdout, " Dyn. Aper.: %5.1f mm^2\n", 1e6*A);
return (A);
}
@@ -1955,7 +1957,7 @@ void LoadTolerances(const char *TolFileName) {
if (Fnum > 0) {
SetTol(Fnum, dx, dy, dr);
} else {
- printf("LoadTolerances: undefined element %s\n", FamName);
+ fprintf(stdout, "LoadTolerances: undefined element %s\n", FamName);
exit_(1);
}
}
@@ -1984,7 +1986,7 @@ void ScaleTolerances(const char *TolFileName, const double scl) {
if (Fnum > 0) {
Scale_Tol(Fnum, scl * dx, scl * dy, scl * dr);
} else {
- printf("ScaleTolerances: undefined element %s\n", FamName);
+ fprintf(stdout, "ScaleTolerances: undefined element %s\n", FamName);
exit_(1);
}
}
@@ -2020,7 +2022,7 @@ void SetKLpar(int Fnum, int Knum, int Order, double kL) {
long int loc;
if (abs(Order) > HOMmax) {
- printf("SetKLPar: Error!....Multipole Order %d > HOMmax %d\n", Order,
+ fprintf(stdout, "SetKLPar: Error!....Multipole Order %d > HOMmax %d\n", Order,
HOMmax);
exit_(1);
}
@@ -2151,12 +2153,12 @@ void set_dbn_rel(const int type, const int n, const double dbn_rel) {
long int j;
double dbn;
- printf("\n");
- printf("Setting Db_%d/b_%d = %6.1e for:\n", n, type, dbn_rel);
- printf("\n");
+ fprintf(stdout, "\n");
+ fprintf(stdout, "Setting Db_%d/b_%d = %6.1e for:\n", n, type, dbn_rel);
+ fprintf(stdout, "\n");
for (j = 0; j <= globval.Cell_nLoc; j++)
if ((Cell[j].Elem.Pkind == Mpole) && (Cell[j].Elem.M->n_design == type)) {
- printf("%s\n", Cell[j].Elem.PName);
+ fprintf(stdout, "%s\n", Cell[j].Elem.PName);
dbn = dbn_rel * Cell[j].Elem.M->PBpar[type + HOMmax];
Cell[j].Elem.M->PBrms[n + HOMmax] = dbn;
Cell[j].Elem.M->PBrnd[n + HOMmax] = normranf();
@@ -2278,13 +2280,13 @@ void SetKL(int Fnum, int Order) {
void set_dx(const int type, const double sigma_x, const double sigma_y) {
long int j;
- printf("\n");
- printf("Setting sigma_x,y = (%6.1e, %6.1e) for b_%d:\n", sigma_x, sigma_y,
+ fprintf(stdout, "\n");
+ fprintf(stdout, "Setting sigma_x,y = (%6.1e, %6.1e) for b_%d:\n", sigma_x, sigma_y,
type);
- printf("\n");
+ fprintf(stdout, "\n");
for (j = 0; j <= globval.Cell_nLoc; j++)
if ((Cell[j].Elem.Pkind == Mpole) && (Cell[j].Elem.M->n_design == type)) {
- printf("%s\n", Cell[j].Elem.PName);
+ fprintf(stdout, "%s\n", Cell[j].Elem.PName);
Cell[j].Elem.M->PdSrms[X_] = sigma_x;
Cell[j].Elem.M->PdSrms[Y_] = sigma_y;
Cell[j].Elem.M->PdSrnd[X_] = normranf();
diff --git a/tracy/tracy/src/read_script.cc b/tracy/tracy/src/read_script.cc
index b7efda9b3c5e0d624e8da855f1b142169fd4c697..ef7d4619770ebebc50c15bb01b5f3b3caf7fc5f4 100644
--- a/tracy/tracy/src/read_script.cc
+++ b/tracy/tracy/src/read_script.cc
@@ -29,6 +29,7 @@
/* files */
+char lat_file[max_str]="voidlattice";
// girder file
char girder_file[max_str];
@@ -66,13 +67,13 @@ void read_script(const char *param_file_name, bool rd_lat, long& CommNo, UserCom
{
- char str[max_str]="voidstring", dummy[max_str]="voidstring", nextpara[max_str]="voidpara";
+ char str[max_str]="voidstring", dummy[max_str]="voidstring",dummy2[max_str]="voidstring", nextpara[max_str]="voidpara";
char in[max_str]; //temporary line with preceding white space
char *line; //line to store the command without preceding white space
char name[max_str]="voidname";
- char lat_file[max_str]="voidlattice";
+ // char lat_file[max_str]="voidlattice";
char EndName[]="void";
-
+
FILE *inf;
const bool prt = false; // for debugging printout each line of input file
long int LineNum=0L;
@@ -92,15 +93,12 @@ void read_script(const char *param_file_name, bool rd_lat, long& CommNo, UserCom
sprintf(full_param_file_name,"%s.prm",dummy);
if (prt) printf("\n reading in script file: %s\n",full_param_file_name);
-
//
inf = file_read(full_param_file_name);
-
// read parameter file, line by line
- // while (fgets(line, max_str, inf) != NULL) {
- while (line = fgets(in, max_str, inf)) {
-
+ // while (fgets(line, max_str, inf) != NULL) {
+ while ((line = fgets(in, max_str, inf))) {
/* kill preceding whitespace generated by "table" key
or "space" key, but leave \n
so we're guaranteed to have something*/
@@ -266,11 +264,12 @@ void read_script(const char *param_file_name, bool rd_lat, long& CommNo, UserCom
// FMA
else if (strcmp("FmapFlag", name) == 0){
strcpy(dummy, "");
-
+ strcpy(dummy2,"");
sscanf(line, "%*s %s",nextpara);
+ //if no definition of loss flag in the lattice file, then "dummy2" is empty string.
if(strcmp(nextpara,"voidpara")!=0)
- sscanf(line, "%*s %s %ld %ld %ld %lf %lf %lf %s",
+ sscanf(line, "%*s %s %ld %ld %ld %lf %lf %lf %s %s",
UserCommandFlag[CommNo]._FmapFlag_fmap_file,
&(UserCommandFlag[CommNo]._FmapFlag_nxpoint),
&(UserCommandFlag[CommNo]._FmapFlag_nypoint),
@@ -278,23 +277,35 @@ void read_script(const char *param_file_name, bool rd_lat, long& CommNo, UserCom
&(UserCommandFlag[CommNo]._FmapFlag_xmax),
&(UserCommandFlag[CommNo]._FmapFlag_ymax),
&(UserCommandFlag[CommNo]._FmapFlag_delta),
- dummy);
+ dummy,dummy2);
if(strcmp(dummy, "true") == 0)
UserCommandFlag[CommNo]._FmapFlag_diffusion = true;
else if(strcmp(dummy, "false") == 0)
UserCommandFlag[CommNo]._FmapFlag_diffusion = false;
+ if(strcmp(dummy2, "true") == 0)
+ UserCommandFlag[CommNo]._FmapFlag_printloss = true;
+ else if(strcmp(dummy2, "false") == 0)
+ UserCommandFlag[CommNo]._FmapFlag_printloss = false;
+ else
+ UserCommandFlag[CommNo]._FmapFlag_printloss = false;
+
+
+ //cout << "debug: " << line << endl;
+ //cout << "dummy2 = " << dummy2 << " lossflag = " << UserCommandFlag[CommNo]._FmapFlag_printloss << endl;
+
// FmapFlag = true;
strcpy(UserCommandFlag[CommNo].CommandStr,name);
}
// FMA dp
else if (strcmp("FmapdpFlag", name) == 0){
strcpy(dummy, "");
+ strcpy(dummy2,"");
sscanf(line, "%*s %s",nextpara);
-
+
if(strcmp(nextpara,"voidpara")!=0)
- sscanf(line, "%*s %s %ld %ld %ld %lf %lf %lf %s",
+ sscanf(line, "%*s %s %ld %ld %ld %lf %lf %lf %s %s",
UserCommandFlag[CommNo]._FmapdpFlag_fmapdp_file,
&(UserCommandFlag[CommNo]._FmapdpFlag_nxpoint),
&(UserCommandFlag[CommNo]._FmapdpFlag_nepoint),
@@ -302,13 +313,24 @@ void read_script(const char *param_file_name, bool rd_lat, long& CommNo, UserCom
&(UserCommandFlag[CommNo]._FmapdpFlag_xmax),
&(UserCommandFlag[CommNo]._FmapdpFlag_emax),
&(UserCommandFlag[CommNo]._FmapdpFlag_z),
- dummy);
+ dummy,dummy2);
if(strcmp(dummy, "true") == 0)
UserCommandFlag[CommNo]._FmapdpFlag_diffusion = true;
else if(strcmp(dummy, "false") == 0)
UserCommandFlag[CommNo]._FmapdpFlag_diffusion = false;
+ if(strcmp(dummy2, "true") == 0)
+ UserCommandFlag[CommNo]._FmapdpFlag_printloss = true;
+ else if(strcmp(dummy2, "false") == 0)
+ UserCommandFlag[CommNo]._FmapdpFlag_printloss = false;
+ else
+ UserCommandFlag[CommNo]._FmapdpFlag_printloss = false;
+
+
+ cout << "debug: " << line << endl;
+ cout << "dummy2 = " << dummy2 << " lossflag = " << UserCommandFlag[CommNo]._FmapdpFlag_printloss << endl;
+
// FmapdpFlag = true;
strcpy(UserCommandFlag[CommNo].CommandStr,name);
}
diff --git a/tracy/tracy/src/soleillib.cc b/tracy/tracy/src/soleillib.cc
index f0650c5361360485a1d2a90b9a22b19e5db50b1f..0f9fa2f7b4bd746b1ca00291cd9e4bcb860f06b9 100644
--- a/tracy/tracy/src/soleillib.cc
+++ b/tracy/tracy/src/soleillib.cc
@@ -476,7 +476,7 @@ void ReadCh(const char *AperFile)
printf("\n");
printf("Loading and setting vacuum apertures to lattice elements...\n");
- while (line=fgets(in, max_str, fp)) {
+ while ((line=fgets(in, max_str, fp))) {
/* kill preceding whitespace generated by "table" key
or "space" key, but leave \n
so we're guaranteed to have something*/
@@ -519,14 +519,14 @@ void ReadCh(const char *AperFile)
Kidnum2 = GetnKid(Fnum2);
if(Kidnum1 != Kidnum2){
printf("\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);
+ " 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"
- " %s %s dxmin = %e, dxmax = %e, dymin = %e, dymax = %e\n",
+ " %s %s dxmin = %+e, dxmax = %+e, dymin = %+e, dymax = %+e\n",
Name1, Name2, dxmin, dxmax, dymin, dymax);
@@ -814,11 +814,14 @@ 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,
- double energy, bool diffusion)
+ double energy, bool diffusion, bool loss)
Purpose:
+
+
Compute a frequency map of Nbx x Nbz points
For each set of initial conditions the particle is tracked over
Nbtour for an energy offset dp
@@ -833,16 +836,18 @@ double get_D(const double df_x, const double df_y)
Input:
FmapFile file to save calculated frequency map analysis
- Nbx horizontal step number
- Nby vertical step number
- xmax horizontal maximum amplitude
- zmax vertical maximum amplitude
- Nbtour number of turn for tracking
- energy particle energy offset
+ Nbx horizontal step number
+ Nby vertical step number
+ xmax horizontal maximum amplitude
+ zmax vertical maximum amplitude
+ Nbtour number of turn for tracking
+ energy particle energy offset
+ 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
Output:
- status true if stable
+ status2 true if stable
false otherwise
Return:
@@ -859,12 +864,18 @@ double get_D(const double df_x, const double df_y)
16/02/03 patch removed
19/07/11 add interface of file defined by user which is used to save calculated
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,
+ 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,
- double energy, bool diffusion)
+ double energy, bool diffusion, bool printloss)
{
- FILE * outf;
+ 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;
@@ -875,15 +886,28 @@ void fmap(const char *FmapFile, long Nbx, long Nbz, long Nbtour, double xmax, do
double nux1 = 0.0, nuz1 = 0.0, nux2 = 0.0, nuz2 = 0.0;
int nb_freq[2] = {0, 0};
long nturn = Nbtour;
- bool status = true;
+ bool status2 = true;
+
struct tm *newtime;
-
+
+//variables of the returned the tracked particle
+ long lastn = 0;
+ long lastpos = 1;
+ ss_vect<double> x1;
+
+
+ //if(loss){
+ char FmapLossFile[S_SIZE+5]=" ";
+ strcpy(FmapLossFile,FmapFile);
+ strcat(FmapLossFile,".loss");
+ //}
+
/* Get time and date */
time_t aclock;
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, "fmap: Entering fmap ... results in %s\n\n",FmapFile);
/* Opening file */
if ((outf = fopen(FmapFile, "w")) == NULL) {
@@ -891,15 +915,31 @@ void fmap(const char *FmapFile, long Nbx, long Nbz, long Nbtour, double xmax, do
exit_(1);
}
+
fprintf(outf,"# TRACY III -- %s -- %s \n", FmapFile, asctime2(newtime));
- fprintf(outf,"# nu = f(x) \n");
+ 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[m] z[m] fx fz"
+ fprintf(outf,"# xi[m] zi[m] fx fz"
" dfx dfz\n");
-
+
+//file with lost particle information
+if(printloss){
+ if ((outf_loss = fopen(FmapLossFile, "w")) == NULL) {
+ fprintf(stdout, "fmap: error while opening file %s\n", FmapLossFile);
+ exit_(1);
+ }
+
+ fprintf(outf_loss,"# TRACY III -- %s -- %s \n", FmapLossFile, asctime2(newtime));
+ fprintf(outf_loss,"# Information of the lost particle: "
+ " lostp = 0 (no particle lost /1(horizontal) /2(vertical) /3(longitudinal) plane \n"
+ " Starting values & Phase space at turn Nturn and position s\n");
+ fprintf(outf_loss,"# xi[m] zi[m] Nturn Plane s[m] x[m]"
+ " xp[rad] z[m] zp[rad] delta ctau\n");
+ }
+
if ((Nbx < 1) || (Nbz < 1))
fprintf(stdout,"fmap: Error Nbx=%ld Nbz=%ld\n",Nbx,Nbz);
@@ -916,25 +956,26 @@ void fmap(const char *FmapFile, long Nbx, long Nbz, long Nbtour, double xmax, do
// Tracking part + NAFF
for (i = 0; i <= Nbx; i++) {
- x = x0 + sqrt((double)i)*xstep;
-// for (i = -Nbx; i <= Nbx; i++) {
- // x = x0 + sgn(i)*sqrt((double)abs(i))*xstep;
-// if (!matlab) fprintf(outf,"\n");
- fprintf(stdout,"\n");
+ x = x0 + sqrt((double)i)*xstep;
+ //fprintf(stdout,"\n");
for (j = 0; j<= Nbz; j++) {
z = z0 + sqrt((double)j)*zstep;
- // for (j = -Nbz; j<= Nbz; j++) {
- // z = z0 + sgn(j)*sqrt((double)abs(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);
+ else
// tracking around closed orbit
- Trac_Simple4DCOD(x,xp,z,zp,energy,ctau,nturn,Tab,&status);
- if (status) { // if trajectory is stable
+ Trac_Simple4DCOD(x,xp,z,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
if (diffusion) { // diffusion
- // shift data for second round NAFF
+ // shift data for second round NAFF
Get_Tabshift(Tab,Tab0,Nbtour,Nbtour);
- // gets frequency vectors
+ // gets frequency vectors
Get_NAFF(NTERM2, Nbtour, Tab0, fx2, fz2, nb_freq);
Get_freq(fx2,fz2,&nux2,&nuz2); // gets nux and nuz
}
@@ -946,34 +987,34 @@ void fmap(const char *FmapFile, long Nbx, long Nbz, long Nbtour, double xmax, do
// printout value
if (!diffusion){
-// fprintf(outf,"%14.6e %14.6e %14.6e %14.6e\n",
-// 1e3*x, 1e3*z, nux1, nuz1);
-// fprintf(stdout,"%14.6e %14.6e %14.6e %14.6e\n",
-// 1e3*x, 1e3*z, nux1, nuz1);
- 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, z, nux1, nuz1);
+ //fprintf(stdout,"%+14.6e %+14.6e %+14.6e %+14.6e\n",
+ // x, z, nux1, nuz1);
}
else {
dfx = nux1 - nux2; dfz = nuz1 - nuz2;
-// fprintf(outf,"%14.6e %14.6e %14.6e %14.6e %14.6e %14.6e %14.6e\n",
-// 1e3*x, 1e3*z, nux1, nuz1, dfx, dfz, get_D(dfx, dfz));
-// fprintf(stdout,"%14.6e %14.6e %14.6e %14.6e %14.6e %14.6e %14.6e\n",
-// 1e3*x, 1e3*z, nux1, nuz1, dfx, dfz, get_D(dfx, dfz));
- 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",
+ 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);
}
+
+ //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],
+ x1[2], x1[3], x1[4], x1[5]);
}
}
+ // Closing output files
fclose(outf);
+ if(printloss)
+ fclose(outf_loss);
}
#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)
@@ -1023,9 +1064,10 @@ void fmap(const char *FmapFile, long Nbx, long Nbz, long Nbtour, double xmax, do
****************************************************************************/
#define NTERM2 10
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)
+ double zmax, double energy, bool diffusion, bool printloss, int numprocs, int myid)
{
- FILE * outf;
+ 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;
@@ -1036,7 +1078,7 @@ void fmap_p(const char *FmapFile_p, long Nbx, long Nbz, long Nbtour, double xmax
double nux1 = 0.0, nuz1 = 0.0, nux2 = 0.0, nuz2 = 0.0;
int nb_freq[2] = {0, 0};
long nturn = Nbtour;
- bool status = true;
+ bool status2 = true;
struct tm *newtime;
char FmapFile[max_str];
@@ -1044,6 +1086,16 @@ void fmap_p(const char *FmapFile_p, long Nbx, long Nbz, long Nbtour, double xmax
strcat(FmapFile,FmapFile_p);
printf("%s\n",FmapFile);
+//variables of the returned the tracked particle
+ long lastn = 0;
+ long lastpos = 1;
+ ss_vect<double> x1;
+
+ char FmapLossFile[S_SIZE+5]=" ";
+ strcpy(FmapLossFile,FmapFile);
+ strcat(FmapLossFile,".loss");
+
+
/* Get time and date */
time_t aclock;
time(&aclock); /* Get time in seconds */
@@ -1052,21 +1104,36 @@ void fmap_p(const char *FmapFile_p, long Nbx, long Nbz, long Nbtour, double xmax
if (trace) printf("Entering fmap ... results in %s\n\n",FmapFile);
/* Opening file */
- if ((outf = fopen(FmapFile, "w")) == NULL)
- {
- fprintf(stdout, "fmap: error while opening file %s\n", FmapFile);
+ if ((outf = fopen(FmapFile, "w")) == NULL){
+ fprintf(stdout, "fmap_p: error while opening file %s\n", FmapFile);
exit_(1);
}
- if(myid==0)
- {
+//file with lost particle information
+ if(printloss){
+ if ((outf_loss = fopen(FmapLossFile, "w")) == NULL) {
+ fprintf(stdout, "fmap_p: error while opening file %s\n", FmapLossFile);
+ exit_(1);
+ }
+ }
+
+ 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");
- }
+
+ if(printloss){
+ fprintf(outf_loss,"# TRACY III -- %s -- %s \n", FmapLossFile, asctime2(newtime));
+ fprintf(outf_loss,"# Information of the lost particle: "
+ " lostp = 0 (no particle lost /1(horizontal) /2(vertical) /3(longitudinal) plane \n"
+ " Starting values & Phase space at turn Nturn and position s\n");
+ fprintf(outf_loss,"# xi[m] zi[m] Nturn Plane s[m] x[m]"
+ " xp[rad] z[m] zp[rad] delta ctau\n");
+ }
+}
if ((Nbx < 1) || (Nbz < 1))
- fprintf(stdout,"fmap: Error Nbx=%ld Nbz=%ld\n",Nbx,Nbz);
+ fprintf(stdout,"fmap_p: Error Nbx=%ld Nbz=%ld\n",Nbx,Nbz);
// steps in both planes
xstep = xmax/sqrt((double)Nbx);
@@ -1082,38 +1149,38 @@ void fmap_p(const char *FmapFile_p, long Nbx, long Nbz, long Nbtour, double xmax
// Tracking part + NAFF
//for (i = 0; i< Nbx; i++)
//Each core or process calculate different region of fmap according to id number. MSChiu 2011/10/13
-int deb,fin;
+ int deb,fin;
int integer,residue;
integer=((int)Nbx)/numprocs;
residue=((int)Nbx)-integer*numprocs;
- printf("myid:%d, integer:%d, residue:%d, numprocs:%d, Nbx:%d\n\n",myid,integer,residue,numprocs,Nbx);
+ fprintf(stdout, "fmap_p: myid:%d, integer:%d, residue:%d, numprocs:%d, Nbx:%ld\n\n",myid,integer,residue,
+ numprocs,Nbx);
//split tracking region (x,z) for each process
- if(myid<residue)
- {
+ if(myid<residue){
deb=myid*(integer+1);
fin=(myid+1)*(integer+1);
}
- else
- {
+ else{
deb=residue*(integer+1)+(myid-residue)*integer;
fin=residue*(integer+1)+(myid+1-residue)*integer;
}
// tracking and FFT, and get the tunes for each particle starts from (x,z)
- for (i = deb; i < fin; i++)
- {
+ for (i = deb; i < fin; i++){
x = x0 + sqrt((double)i)*xstep;
- fprintf(stdout,"\n");
- for (j = 0; j< Nbz; j++)
- {
+ // fprintf(stdout,"\n");
+ for (j = 0; j< Nbz; j++) {
z = z0 + sqrt((double)j)*zstep;
-
// tracking around closed orbit
- Trac_Simple4DCOD(x,xp,z,zp,energy,ctau,nturn,Tab,&status);
- if (status) // if trajectory is stable
+ if(printloss)
+ Trac_Simple4DCOD(x,xp,z,zp,energy,ctau,nturn,Tab,lastn, lastpos, x1, &status2);
+ else
+ Trac_Simple4DCOD(x,xp,z,zp,energy,ctau,nturn,Tab,&status2);
+
+ if (status2) // if trajectory is stable
{
// gets frequency vectors
Get_NAFF(NTERM2, Nbtour, Tab, fx, fz, nb_freq);
@@ -1136,21 +1203,31 @@ int deb,fin;
// 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, z, nux1, nuz1);
+ //fprintf(stdout,"%+14.6e %+14.6e %+14.6e %+14.6e\n",x, z, nux1, nuz1);
}
else
{
dfx = nux1 - nux2;
dfz = nuz1 - nuz2;
- 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, 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);
}
+
+ //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],
+ x1[2], x1[3], x1[4], x1[5]);
}
}
fclose(outf);
+
+ if(printloss)
+ fclose(outf_loss);
+
}
#undef NTERM2
@@ -1180,10 +1257,11 @@ int deb,fin;
xmax horizontal maximum amplitude
emax maximum energy
z vertical amplitude
- diffusion flag to calculate tune diffusion
+ 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
Output:
- status true if stable
+ status2 true if stable
false otherwise
Return:
@@ -1200,12 +1278,17 @@ int deb,fin;
23/10/04 for 6D turn off diffusion automatically and horizontal amplitude
is negative for negative enrgy offset since this is true for the cod
19/07/11 add features to save calculated fmapdp in the user defined file
+
+ 18/06/2013 by Jianfeng Zhang @ LAL
+
+ Add bool flag to print out the last information of the tracked particle
****************************************************************************/
#define NTERM2 10
void fmapdp(const char *FmapdpFile, long Nbx, long Nbe, long Nbtour, double xmax, double emax,
- double z, bool diffusion)
+ double z, 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;
@@ -1216,9 +1299,18 @@ void fmapdp(const char *FmapdpFile, long Nbx, long Nbe, long Nbtour, double xmax
int nb_freq[2] = {0, 0};
long nturn = Nbtour;
- bool status=true;
+ bool status2=true;
struct tm *newtime;
+ //variables of the returned the tracked particle
+ long lastn = 0;
+ long lastpos = 1;
+ ss_vect<double> x1;
+
+ char FmapdpLossFile[S_SIZE+5]=" ";
+ strcpy(FmapdpLossFile,FmapdpFile);
+ strcat(FmapdpLossFile,".loss");
+
/* Get time and date */
time_t aclock;
time(&aclock); /* Get time in seconds */
@@ -1226,7 +1318,7 @@ void fmapdp(const char *FmapdpFile, long Nbx, long Nbe, long Nbtour, double xmax
if (diffusion && globval.Cavity_on == false) nturn = 2*Nbtour;
- if (trace) printf("Entering fmap ... results in %s\n\n",FmapdpFile);
+ if (trace) fprintf(stdout, "fmapdp: Entering fmapdp ... results in %s\n\n",FmapdpFile);
/* Opening file */
if ((outf = fopen(FmapdpFile, "w")) == NULL) {
@@ -1239,6 +1331,22 @@ void fmapdp(const char *FmapdpFile, long Nbx, long Nbe, long Nbtour, double xmax
// fprintf(outf,"# dp[%%] x[mm] fx fz dfx dfz\n");
fprintf(outf,"# dp[m] x[m] fx fz dfx dfz\n");
+//file with lost particle information
+if(printloss){
+
+ if ((outf_loss = fopen(FmapdpLossFile, "w")) == NULL) {
+ fprintf(stdout, "fmapdp: error while opening file %s\n", FmapdpLossFile);
+ exit_(1);
+ }
+
+ fprintf(outf_loss,"# TRACY III -- %s -- %s \n", FmapdpLossFile, asctime2(newtime));
+ fprintf(outf_loss,"# Information of the lost particle: "
+ " lostp = 0 (no particle lost /1(horizontal) /2(vertical) /3(longitudinal) plane \n"
+ " Starting values & Phase space at turn Nturn and position s\n");
+
+ fprintf(outf_loss,"# dpi xi[m] Nturn lostp s[m] x[m]"
+ " xp[rad] z[m] zp[rad] delta ctau\n");
+ }
if ((Nbx <= 1) || (Nbe <= 1))
fprintf(stdout,"fmapdp: Error Nbx=%ld Nbe=%ld\n",Nbx,Nbe);
@@ -1250,22 +1358,24 @@ fprintf(outf,"# dp[m] x[m] fx fz dfx
for (i = 0; i <= Nbe; i++) {
dp = -emax + i*estep;
-// if (!matlab) fprintf(outf,"\n");
- fprintf(stdout,"\n");
for (j = 0; j<= Nbx; j++) {
-// for (j = -Nbx; j<= Nbx; j++) {
- // IF 6D Tracking diffusion turn off and x negative for dp negative
+ // If 6D Tracking diffusion turn off and x negative for dp negative
if ((globval.Cavity_on == true) && (dp < 0.0)){
- // x = x0 - sgn(j)*sqrt((double)abs(j))*xstep;
x = x0 - sqrt((double)j)*xstep;
diffusion = false;
}
- else
+ else{
// x = x0 + sgn(j)*sqrt((double)abs(j))*xstep;
x = x0 + sqrt((double)j)*xstep;
- Trac_Simple4DCOD(x,xp,z,zp,dp,ctau,nturn,Tab,&status);
- if (status) {
+ }
+
+if(printloss)
+ Trac_Simple4DCOD(x,xp,z,zp,dp,ctau,nturn,Tab,lastn, lastpos, x1, &status2);
+ else
+ Trac_Simple4DCOD(x,xp,z,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
if (diffusion) { // diffusion
@@ -1281,28 +1391,27 @@ fprintf(outf,"# dp[m] x[m] fx fz dfx
// printout value
if (!diffusion){
-// fprintf(outf,"%14.6e %14.6e %14.6e %14.6e\n",
-// 1e2*dp, 1e3*x, nux1, nuz1);
-// fprintf(stdout,"%14.6e %14.6e %14.6e %14.6e\n",
-// 1e2*dp, 1e3*x, nux1, nuz1);
- 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, nuz1);
+ //fprintf(stdout,"% 10.6e % 10.6e % 10.6e % 10.6e\n", dp, x, nux1, nuz1);
}
else {
dfx = nux2 - nux1; dfz = nuz2 - nuz1;
-// fprintf(outf,"%14.6e %14.6e %14.6e %14.6e %14.6e %14.6e %14.6e\n",
-// 1e2*dp, 1e3*x, nux1, nuz2, dfx, dfz, get_D(dfx, dfz));
-// fprintf(stdout,"%14.6e %14.6e %14.6e %14.6e %14.6e %14.6e %14.6e\n",
-// 1e2*dp, 1e3*x, nux1, nuz2, dfx, dfz, get_D(dfx, dfz));
- fprintf(outf,"% 10.6e % 10.6e % 10.6e % 10.6e % 10.6e % 10.6e\n",
- dp, x, nux1, nuz2, dfx, dfz);
- fprintf(stdout,"% 10.6e % 10.6e % 10.6e % 10.6e % 10.6e % 10.6e\n",
+ fprintf(outf,"%+10.6e %+10.6e %+10.6e %+10.6e %+10.6e %+10.6e\n",
dp, x, nux1, nuz2, dfx, dfz);
+ //fprintf(stdout,"% 10.6e % 10.6e % 10.6e % 10.6e % 10.6e % 10.6e\n",
+ // dp, x, nux1, nuz2, dfx, dfz);
}
+
+ if(printloss)
+ fprintf(outf_loss," %+8.3e %+8.3e %8ld %2d %+12.3e %+10.6e %+10.6e %+10.6e %+10.6e %+10.6e %+10.6e \n",
+ dp, x, lastn, status.lossplane, Cell[lastpos].S, x1[0],x1[1], x1[2], x1[3], x1[4], x1[5]);
+
}
}
-
fclose(outf);
+
+if(printloss)
+ fclose(outf_loss);
}
#undef NTERM2
@@ -1337,7 +1446,7 @@ fprintf(outf,"# dp[m] x[m] fx fz dfx
myid process used to do parallel computing
Output:
- status true if stable
+ status2 true if stable
false otherwise
Return:
@@ -1352,12 +1461,16 @@ fprintf(outf,"# dp[m] x[m] fx fz dfx
Comments:
14/11/2011 add features to parallel calculate fmapdp.
Merged with the version written by Mao-sen Qiu at Taiwan light source.
+
+ 18/06/2013 by Jianfeng Zhang @ LAL
+ Add bool flag to print out the last information of the tracked particle
****************************************************************************/
#define NTERM2 10
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 z, 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;
@@ -1368,13 +1481,22 @@ void fmapdp_p(const char *FmapdpFile_p, long Nbx, long Nbe, long Nbtour, double
int nb_freq[2] = {0, 0};
long nturn = Nbtour;
- bool status=true;
+ bool status2=true;
struct tm *newtime;
char FmapdpFile[max_str];
sprintf(FmapdpFile,"%d",myid);
strcat(FmapdpFile,FmapdpFile_p);
- printf("%s\n",FmapdpFile);
+ fprintf(stdout, "fmap_dp_p: %s\n",FmapdpFile);
+
+//variables of the returned the tracked particle
+ long lastn = 0;
+ long lastpos = 1;
+ ss_vect<double> x1;
+
+ char FmapdpLossFile[S_SIZE+5]=" ";
+ strcpy(FmapdpLossFile,FmapdpFile);
+ strcat(FmapdpLossFile,".loss");
/* Get time and date */
time_t aclock;
@@ -1383,24 +1505,42 @@ char FmapdpFile[max_str];
if (diffusion && globval.Cavity_on == false) nturn = 2*Nbtour;
- if (trace) printf("Entering fmap ... results in %s\n\n",FmapdpFile);
+ if (trace) fprintf(stdout,"fmapdp_p: Entering fmap ... results in %s\n\n",FmapdpFile);
/* Opening file */
if ((outf = fopen(FmapdpFile, "w")) == NULL) {
- fprintf(stdout, "fmapdp: error while opening file %s\n", FmapdpFile);
+ fprintf(stdout, "fmapdp_p: error while opening file %s\n", FmapdpFile);
exit_(1);
}
+if(printloss){
+ if ((outf_loss = fopen(FmapdpLossFile, "w")) == NULL) {
+ fprintf(stdout, "fmapdp_p: error while opening file %s\n", FmapdpLossFile);
+ exit_(1);
+ }
+ }
+
if(myid==0)
{
fprintf(outf,"# TRACY III -- %s -- %s \n", FmapdpFile_p, asctime2(newtime));
fprintf(outf,"# nu = f(x) \n");
// fprintf(outf,"# dp[%%] x[mm] fx fz dfx dfz\n");
fprintf(outf,"# dp[m] x[m] fx fz dfx dfz\n");
- }
+
+
+ if(printloss){
+ fprintf(outf_loss,"# TRACY III -- %s -- %s \n", FmapdpLossFile, asctime2(newtime));
+ fprintf(outf_loss,"# Information of the lost particle: "
+ " lostp = 0 (no particle lost /1(horizontal) /2(vertical) /3(longitudinal) plane \n"
+ " Starting values & Phase space at turn Nturn and position s\n");
+
+ fprintf(outf_loss,"# dpi xi[m] Nturn lostp s[m] x[m]"
+ " xp[rad] z[m] zp[rad] delta ctau\n");
+ }
+ }
if ((Nbx <= 1) || (Nbe <= 1))
- fprintf(stdout,"fmapdp: Error Nbx=%ld Nbe=%ld\n",Nbx,Nbe);
+ fprintf(stdout,"fmapdp_p: Error Nbx=%ld Nbe=%ld\n",Nbx,Nbe);
xp = xp0;
zp = zp0;
@@ -1415,7 +1555,7 @@ if ((Nbx <= 1) || (Nbe <= 1))
integer=((int)Nbe)/numprocs;
residue=((int)Nbe)-integer*numprocs;
- printf("myid:%d, integer:%d, resideu:%d, numprocs:%d, Nbe:%d\n\n",myid,integer,residue,numprocs,Nbe);
+ fprintf(stdout, "fmapdp_p: myid:%d, integer:%d, resideu:%d, numprocs:%d, Nbe:%ld\n\n",myid,integer,residue,numprocs,Nbe);
//split tracking region for each process
if(myid<residue)
@@ -1437,7 +1577,7 @@ if ((Nbx <= 1) || (Nbe <= 1))
// for (i = 0; i <= Nbe; i++) {
// dp = -emax + i*estep;
// if (!matlab) fprintf(outf,"\n");
- fprintf(stdout,"\n");
+ //fprintf(stdout,"\n");
for (j = 0; j<= Nbx; j++) {
// for (j = -Nbx; j<= Nbx; j++) {
@@ -1447,11 +1587,17 @@ if ((Nbx <= 1) || (Nbe <= 1))
x = x0 - sqrt((double)j)*xstep;
diffusion = false;
}
- else
+ else{
// x = x0 + sgn(j)*sqrt((double)abs(j))*xstep;
x = x0 + sqrt((double)j)*xstep;
- Trac_Simple4DCOD(x,xp,z,zp,dp,ctau,nturn,Tab,&status);
- if (status) {
+ }
+
+ if(printloss)
+ Trac_Simple4DCOD(x,xp,z,zp,dp,ctau,nturn,Tab,lastn, lastpos, x1, &status2);
+ else
+ Trac_Simple4DCOD(x,xp,z,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
if (diffusion) { // diffusion
@@ -1467,28 +1613,29 @@ if ((Nbx <= 1) || (Nbe <= 1))
// printout value
if (!diffusion){
-// fprintf(outf,"%14.6e %14.6e %14.6e %14.6e\n",
-// 1e2*dp, 1e3*x, nux1, nuz1);
-// fprintf(stdout,"%14.6e %14.6e %14.6e %14.6e\n",
-// 1e2*dp, 1e3*x, nux1, nuz1);
- 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, nuz1);
+ //fprintf(stdout,"% 10.6e % 10.6e % 10.6e % 10.6e\n", dp, x, nux1, nuz1);
}
else {
dfx = nux2 - nux1; dfz = nuz2 - nuz1;
-// fprintf(outf,"%14.6e %14.6e %14.6e %14.6e %14.6e %14.6e %14.6e\n",
-// 1e2*dp, 1e3*x, nux1, nuz2, dfx, dfz, get_D(dfx, dfz));
-// fprintf(stdout,"%14.6e %14.6e %14.6e %14.6e %14.6e %14.6e %14.6e\n",
-// 1e2*dp, 1e3*x, nux1, nuz2, dfx, dfz, get_D(dfx, dfz));
- fprintf(outf,"% 10.6e % 10.6e % 10.6e % 10.6e % 10.6e % 10.6e\n",
- dp, x, nux1, nuz2, dfx, dfz);
- fprintf(stdout,"% 10.6e % 10.6e % 10.6e % 10.6e % 10.6e % 10.6e\n",
+ fprintf(outf,"% 10.6e % 10.6e % 10.6e % 10.6e % 10.6e % 10.6e\n",
dp, x, nux1, nuz2, dfx, dfz);
+ //fprintf(stdout,"% 10.6e % 10.6e % 10.6e % 10.6e % 10.6e % 10.6e\n",
+ // dp, x, nux1, nuz2, dfx, dfz);
}
+
+ if(printloss)
+ fprintf(outf_loss," %-+8.3f %-+8.3f %-8ld %2d %+12.3f %+10.6f %+10.6f %+10.6f %+10.6f %+10.6f %+10.6f \n",
+ dp, x, lastn, status.lossplane, Cell[lastpos].S, x1[0],x1[1], x1[2], x1[3], x1[4], x1[5]);
+
}
}
-
+
+ // Closing files
fclose(outf);
+
+ if(printloss)
+ fclose(outf_loss);
}
#undef NTERM2
@@ -1740,7 +1887,6 @@ void PhasePoly(long pos, double x0,double px0, double z0, double pz0, double del
fprintf(outf,
"# num x xp z zp dp ctau\n");
- trace = true;
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;
@@ -1961,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("xcod=%.5e mm zcod=% .5e mm \n", globval.CODvect[0]*1e3, globval.CODvect[2]*1e3);
+ printf("Enveloppe: xcod=%.5e mm zcod=% .5e mm \n", globval.CODvect[0]*1e3, globval.CODvect[2]*1e3);
if ((outf = fopen(fic, "w")) == NULL)
{
@@ -1984,10 +2130,10 @@ 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("Unstable motion ...\n"); exit_(1);
+ printf("Enveloppe: Unstable motion ...\n"); exit_(1);
}
- fprintf(outf,"%6.2f % .5e % .5e % .5e % .5e % .5e % .5e\n",
+ fprintf(outf,"%6.2f %+.5e %+.5e %+.5e %+.5e %+.5e %+.5e\n",
Cell.S, x1[0],x1[1],x1[2],x1[3],x1[4],x1[5]);
}
}
@@ -3403,7 +3549,6 @@ void MomentumAcceptance(char *MomAccFile, long deb, long fin,
struct tm *newtime; // for time
Vector codvector[Cell_nLocMax];
bool cavityflag, radiationflag;
- bool trace=true;
x0.zero();
@@ -3801,7 +3946,6 @@ void MomentumAcceptance_p(char *_MomAccFile, long deb, long fin, double ep_min,
struct tm *newtime; // for time
Vector codvector[Cell_nLocMax];
bool cavityflag, radiationflag;
- bool trace=true;
x0.zero();
@@ -3938,7 +4082,7 @@ void MomentumAcceptance_p(char *_MomAccFile, long deb, long fin, double ep_min,
/***************************************************************/
//split tracking element region for each process
-//Eace core or process calculate different region of fmap according to id number. MSChiu 2011/10/13
+//Each core or process calculate different region of fmap according to id number. MSChiu 2011/10/13
int debN,finN;
int integer,residue;
@@ -3951,7 +4095,7 @@ if(fin < deb){
integer=(fin-deb+1)/numprocs;
residue=(fin-deb+1)-integer*numprocs;
- printf("myid:%d, integer:%d, resideu:%d, numprocs:%d, Nbx:%d\n",myid,integer,residue,numprocs);
+ printf("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
@@ -4900,7 +5044,6 @@ void Phase2(long pos, double x,double px,double y, double py,double energy,
fprintf(outf,
"# num x xp z zp dp ctau\n");
- trace = true;
Trac(x,px,y,py,energy,ctau, Nbtour,pos, lastn, lastpos, outf);
fclose(outf);
}
@@ -4929,7 +5072,6 @@ void Phase3(long pos, double x,double px,double y, double py,double energy,
fprintf(outf,
"# num x xp z zp dp ctau\n");
- trace = true;
x1[0] = x; x1[1] = px; x1[2] = y;
x1[3] = py; x1[4] = energy; x1[5] = ctau;
Cell_Pass(0L, pos-1L, x1, lastpos);
@@ -5878,7 +6020,7 @@ void ReadFieldErr(const char *FieldErrorFile)
printf("\n");
/* read lines*/
- while (line=fgets(in, max_str, inf)) {
+ while ((line = fgets(in, max_str, inf))) {
/* kill preceding whitespace generated by "table" key
or "space" key, but leave \n
@@ -6533,7 +6675,7 @@ void PrintTrack(const char *TrackFile, double x, double px, double y,double py,
Cell_Pass(pos, pos, x2, lastpos);
//check whether particle is lost
if (!CheckAmpl(x2, pos)){
- fprintf(stderr,"Error!!! %d turn, Particle lost at element: %s!",
+ fprintf(stderr,"Error!!! %ld turn, Particle lost at element: %s!",
lastn, Cell[pos].Elem.PName);
exit_(1);
}
@@ -6552,7 +6694,7 @@ void PrintTrack(const char *TrackFile, double x, double px, double y,double py,
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);
}
- fprintf(outf,"%6d %s %8.4e %12.4e %12.4e %12.4e %12.4e %12.4e %12.4e %4ld \n",
+ fprintf(outf,"%6ld %s %8.4e %12.4e %12.4e %12.4e %12.4e %12.4e %12.4e %4ld \n",
pos,Cell[pos].Elem.PName,Cell[pos].S, 1e3*xf[x_], 1e3*xf[px_],
1e3*xf[y_], 1e3*xf[py_],1e2*xf[delta_], 1e3*xf[ct_],lastn);
@@ -6585,7 +6727,6 @@ void PrintTrack(const char *TrackFile, double x, double px, double y,double py,
void ReadVirtualSkewQuad(const char *VirtualSkewQuadFile)
{
int nqtcorr= 0; /* Number of virtual skew quadrupoles */
- int qtcorrlist[max_str]; /* virtual skew quad list */
int i=0;
long mOrder = 0L;
double qtcorr[max_str]; /* virtual skew quad value in Amps */