From 1c95b76cb8b17666354f690cf686cd67e57c3225 Mon Sep 17 00:00:00 2001
From: zhang <zhang@9a6e40ed-f3a0-4838-9b4a-bf418f78e88d>
Date: Thu, 22 Dec 2011 13:06:13 +0000
Subject: [PATCH] 1) Add feature to read virtual source of coupling for soleil
lattice. 2) Add feature to call COD correction. 3) Add fmap_p( ), fmapdp_p(
), momacceptance_p( ), to do parallel computation of fmap, fmapdp,
momacceptance.
---
tracy/tracy/inc/soleillib.h | 15 +-
tracy/tracy/src/soleillib.cc | 1333 +++++++++++++++++++++++++++++-----
2 files changed, 1171 insertions(+), 177 deletions(-)
diff --git a/tracy/tracy/inc/soleillib.h b/tracy/tracy/inc/soleillib.h
index 3ad771b..c34ed54 100644
--- a/tracy/tracy/inc/soleillib.h
+++ b/tracy/tracy/inc/soleillib.h
@@ -36,6 +36,9 @@ void Multipole_thinsext(const char *fic_hcorr, const char *fic_vcorr, const cha
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);
+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);
void Trac_Tab(double x, double px, double y, double py, double dp,
long nmax, long pos, long *lastn, long *lastpos, FILE *outf1, double Tx[][NTURN]);
@@ -55,8 +58,12 @@ void TunesShiftWithEnergy(const char *NudpFile,long Nb, long Nbtour, 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);
+void fmap_p(const char *FmapFile, long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
+ double energy, bool diffusion, int numprocs, int myid);
void fmapdp(const char *FmapdpFile, long Nbx, long Nbe, long Nbtour, double xmax, double emax,
double z, bool diffusion);
+void fmapdp_p(const char *FmapdpFile, long Nbx, long Nbe, long Nbtour, double xmax, double emax,
+ double z, bool diffusion, int numprocs, int myid);
void Nu_Naff(void);
@@ -94,6 +101,9 @@ double get_RFVoltage(const int Fnum);
void set_RFVoltage(const int Fnum, const double V_RF);
+/* close orbit correction */
+void CODCorrect(const char *hcorr_file,const char *vcorr_file,int n_orbit,int nwh,int nwv);
+
/* Read multipole errors from a file for soleil*/
void ReadFieldErr(const char *FieldErrorFile);
@@ -111,7 +121,10 @@ void Add_bnL_sys_elem(const int Fnum, const int Knum, const char *keywrd,
void AddCorrQtErr_fam(char const *fic, const int Fnum, const double conv, const char *keywrd, const double r0,
const int n, const double Bn, const double An);
-
+
+/* Read the setting of skew quadrupoles from a file; for soleil lattice */
+void ReadVirtualSkewQuad(const char *VirtualSkewQuadFile);
+
/* fit tunes for soleil lattice, in which each quadrupole is cut into two parts*/
void FitTune4(long qf1,long qf2, long qd1, long qd2, double nux, double nuy);
//print the coordinates at lattice elements
diff --git a/tracy/tracy/src/soleillib.cc b/tracy/tracy/src/soleillib.cc
index b774ce4..f0650c5 100644
--- a/tracy/tracy/src/soleillib.cc
+++ b/tracy/tracy/src/soleillib.cc
@@ -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) printf("%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);
@@ -973,6 +973,189 @@ void fmap(const char *FmapFile, long Nbx, long Nbz, long Nbtour, double xmax, do
}
#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)
+
+ Purpose:
+ Parallelized version of fmap( ).
+ 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
+
+ Frequency map is based on fixed beam energy, trace x versus z,
+ or, tracking transverse dynamic aperture for fixed momentum
+ (usually, on-momentum) particle.
+
+ The stepsize follows a square root law
+
+ Results in fmap.out
+
+ Input:
+ FmapFile_p 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
+ matlab set file print format for matlab post-process; specific for nsrl-ii
+ numprocs number of processes used to do parallel computing
+ myid process used to do parallel computing
+
+ Output:
+ status true if stable
+ false otherwise
+
+ Return:
+ none
+
+ Global variables:
+ none
+
+ Specific functions:
+ Trac_Simple, Get_NAFF
+
+ Comments:
+ 14/11/2011 add feature to do parallel computing of frequency map analysis.
+ 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, int numprocs, int myid)
+{
+ FILE * outf;
+ 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 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;
+ int nb_freq[2] = {0, 0};
+ long nturn = Nbtour;
+ bool status = true;
+ struct tm *newtime;
+
+ char FmapFile[max_str];
+ sprintf(FmapFile,"%d",myid);
+ strcat(FmapFile,FmapFile_p);
+ printf("%s\n",FmapFile);
+
+ /* 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);
+
+ /* Opening file */
+ if ((outf = fopen(FmapFile, "w")) == NULL)
+ {
+ fprintf(stdout, "fmap: error while opening file %s\n", FmapFile);
+ 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 ((Nbx < 1) || (Nbz < 1))
+ fprintf(stdout,"fmap: Error Nbx=%ld Nbz=%ld\n",Nbx,Nbz);
+
+ // steps in both planes
+ xstep = xmax/sqrt((double)Nbx);
+ zstep = zmax/sqrt((double)Nbz);
+
+ // double number of turn if diffusion to compute
+ if (diffusion) nturn = 2*Nbtour;
+
+ // px and pz zeroed
+ xp = xp0;
+ zp = zp0;
+
+// 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 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);
+
+ //split tracking region (x,z) for each process
+ if(myid<residue)
+ {
+ deb=myid*(integer+1);
+ fin=(myid+1)*(integer+1);
+ }
+ 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++)
+ {
+ x = x0 + sqrt((double)i)*xstep;
+
+ 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
+ {
+ // 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
+ 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
+ }
+ } // unstable trajectory
+ else
+ { //zeroing output
+ nux1 = 0.0; nuz1 = 0.0;
+ nux2 = 0.0; nuz2 = 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);
+ }
+ 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);
+ }
+ }
+ }
+
+ fclose(outf);
+}
+#undef NTERM2
+
+
+
/****************************************************************************/
/* void fmapdp(const char *FmapdpFile, long Nbx, long Nbe, long Nbtour, double xmax, double emax,
double z, bool diffusion)
@@ -1124,104 +1307,291 @@ fprintf(outf,"# dp[m] x[m] fx fz dfx
#undef NTERM2
/****************************************************************************/
-/* void TunesShiftWithEnergy(long Nb, long Nbtour, double emax)
+/* 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)
+
Purpose:
- Computes tunes versus energy offset by tracking
- by linear energy step between -emax and emax
+ Parallel version of fmapdp( ).
+ 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
+
+ Frequency map is based on fixed vertical amplitude z, trace x versus energy,
+ or, tracking x for off-momentum particle.
+
+ The stepsize follows a square root law
- Input:
- Nb+1 numbers of points
- NbTour number of turns for tracking
- emax maximum energy
+ Results in fmapdp.out
+
+ Input:
+ FmapdpFile_p file to save calculated frequency map analysis
+ Nbx horizontal step number
+ Nbe energy step number
+ Nbtour number of turns for tracking
+ xmax horizontal maximum amplitude
+ emax maximum energy
+ z vertical amplitude
+ diffusion flag to calculate tune diffusion
+ numprocs Number of processes used to do parallel computing
+ myid process used to do parallel computing
Output:
- none
+ status true if stable
+ false otherwise
Return:
none
Global variables:
- trace
+ none
Specific functions:
Trac_Simple, Get_NAFF
Comments:
- none
-
+ 14/11/2011 add features to parallel calculate fmapdp.
+ Merged with the version written by Mao-sen Qiu at Taiwan light source.
****************************************************************************/
-#define NTERM 4
-void TunesShiftWithEnergy(const char *NudpFile,long Nb, long Nbtour, double emax)
+#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)
{
- FILE * outf;
+ FILE * outf;
+ 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, 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;
- 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 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;
- int nb_freq[2] = {0, 0};
- bool status = true;
- struct tm *newtime;
+ int nb_freq[2] = {0, 0};
+ long nturn = Nbtour;
+ bool status=true;
+ struct tm *newtime;
- /* Get time and date */
- newtime = GetTime();
+char FmapdpFile[max_str];
+ sprintf(FmapdpFile,"%d",myid);
+ strcat(FmapdpFile,FmapdpFile_p);
+ printf("%s\n",FmapdpFile);
- if (!trace) printf("\n Entering TunesShiftWithEnergy ...\n\n");
+ /* Get time and date */
+ time_t aclock;
+ time(&aclock); /* Get time in seconds */
+ newtime = localtime(&aclock); /* Convert time to struct */
- /* Opening file */
- if ((outf = fopen(NudpFile, "w")) == NULL) {
- fprintf(stdout, "NuDp: error while opening file %s\n", NudpFile);
- exit_(1);
- }
+ if (diffusion && globval.Cavity_on == false) nturn = 2*Nbtour;
- 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");
+ if (trace) printf("Entering fmap ... results in %s\n\n",FmapdpFile);
- if (Nb <= 1L)
- fprintf(stdout,"NuDp: Error Nb=%ld\n",Nb);
+ /* Opening file */
+ if ((outf = fopen(FmapdpFile, "w")) == NULL) {
+ fprintf(stdout, "fmapdp: error while opening file %s\n", FmapdpFile);
+ exit_(1);
+ }
- // start loop over energy
- dp0 = -emax;
+ 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 ((Nbx <= 1) || (Nbe <= 1))
+ fprintf(stdout,"fmapdp: Error Nbx=%ld Nbe=%ld\n",Nbx,Nbe);
- for (i = 0L; i < Nb; i++) {
- dp = dp0 + i*emax/(Nb-1)*2;
- x = x0 ;
- xp = xp0 ;
- z = z0 ;
- zp = zp0 ;
- ctau = ctau0;
-
- Trac_Simple4DCOD(x,xp,z,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
- }
- else {
- nux1 = 0.0; nuz1 = 0.0;
- status = true;
- }
-
- long lastpos=0L;
- getcod(dp, lastpos); // get cod for printout
+ xp = xp0;
+ zp = zp0;
+ xstep = xmax/sqrt((double)Nbx);
+ estep = 2.0*emax/Nbe;
- 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],
- globval.CODvect[2], globval.CODvect[3]);
+// for (i = 0; i < Nbe; i++)
+// Eace core or process calculate different region of fmapdp according to id number. MSChiu 2011/10/13
+ int deb,fin;
+ int integer,residue;
+ integer=((int)Nbe)/numprocs;
+ residue=((int)Nbe)-integer*numprocs;
- 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],
- globval.CODvect[2], globval.CODvect[3]);
- }
+ printf("myid:%d, integer:%d, resideu:%d, numprocs:%d, Nbe:%d\n\n",myid,integer,residue,numprocs,Nbe);
- fclose(outf);
-}
-#undef NTERM
+ //split tracking region for each process
+ if(myid<residue)
+ {
+ deb=myid*(integer+1);
+ fin=(myid+1)*(integer+1);
+ }
+ else
+ {
+ deb=residue*(integer+1)+(myid-residue)*integer;
+ fin=residue*(integer+1)+(myid+1-residue)*integer;
+ }
+
+ //begin tracking and FFT, and get tunes for the particle starts from (x,p)
+ for (i = deb; i < fin; i++)
+ {
+ dp = -emax + i*estep;
+
+ // for (i = 0; i <= Nbe; i++) {
+ // dp = -emax + i*estep;
+// if (!matlab) fprintf(outf,"\n");
+ fprintf(stdout,"\n");
+ for (j = 0; j<= Nbx; j++) {
+// for (j = -Nbx; j<= Nbx; j++) {
+
+ // IF 6D Tracking diffusion turn off and x negative for dp negative
+ if ((globval.Cavity_on == true) && (dp < 0.0)){
+ // x = x0 - sgn(j)*sqrt((double)abs(j))*xstep;
+ x = x0 - sqrt((double)j)*xstep;
+ diffusion = false;
+ }
+ 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) {
+ Get_NAFF(NTERM2, Nbtour, Tab, fx, fz, nb_freq);
+ Get_freq(fx,fz,&nux1,&nuz1); // gets nux and nuz
+ 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
+ }
+ } // unstable trajectory
+ else { //zeroing output
+ nux1 = 0.0; nuz1 = 0.0;
+ nux2 = 0.0; nuz2 = 0.0;
+ }
+
+ // 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);
+ }
+ 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",
+ dp, x, nux1, nuz2, dfx, dfz);
+ }
+ }
+ }
+
+ fclose(outf);
+}
+#undef NTERM2
+
+
+/****************************************************************************/
+/* void TunesShiftWithEnergy(long Nb, long Nbtour, double emax)
+
+ Purpose:
+ Computes tunes versus energy offset by tracking
+ by linear energy step between -emax and emax
+
+ Input:
+ Nb+1 numbers of points
+ NbTour number of turns for tracking
+ emax maximum energy
+
+ Output:
+ none
+
+ Return:
+ none
+
+ Global variables:
+ trace
+
+ Specific functions:
+ Trac_Simple, Get_NAFF
+
+ Comments:
+ none
+
+****************************************************************************/
+#define NTERM 4
+void TunesShiftWithEnergy(const char *NudpFile,long Nb, long Nbtour, double emax)
+{
+ FILE * outf;
+
+ 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 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;
+ int nb_freq[2] = {0, 0};
+ bool status = true;
+ struct tm *newtime;
+
+ /* Get time and date */
+ newtime = GetTime();
+
+ if (!trace) printf("\n Entering TunesShiftWithEnergy ...\n\n");
+
+ /* Opening file */
+ if ((outf = fopen(NudpFile, "w")) == NULL) {
+ fprintf(stdout, "NuDp: error while opening file %s\n", NudpFile);
+ exit_(1);
+ }
+
+ 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");
+
+ if (Nb <= 1L)
+ fprintf(stdout,"NuDp: Error Nb=%ld\n",Nb);
+
+ // start loop over energy
+ dp0 = -emax;
+
+ for (i = 0L; i < Nb; i++) {
+ dp = dp0 + i*emax/(Nb-1)*2;
+ x = x0 ;
+ xp = xp0 ;
+ z = z0 ;
+ zp = zp0 ;
+ ctau = ctau0;
+
+ Trac_Simple4DCOD(x,xp,z,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
+ }
+ else {
+ nux1 = 0.0; nuz1 = 0.0;
+ status = true;
+ }
+
+ long lastpos=0L;
+ getcod(dp, lastpos); // get cod for printout
+
+
+ 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],
+ 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],
+ globval.CODvect[2], globval.CODvect[3]);
+ }
+
+ fclose(outf);
+}
+#undef NTERM
@@ -2913,51 +3283,463 @@ void SetSkewQuad(void)
Specific functions:
GetElem, SetKLpar, GetKpar
- Comments:
- none
+ Comments:
+ none
+
+****************************************************************************/
+// void SetDecapole(void)
+// {
+// FILE *fi;
+// const char fic_deca[] ="/home/nadolski/soltracy/deca.dat";
+// int i;
+// double mKL[56]; /* array for skew quad tilt*/
+// CellType Cell;
+// long mOrder=5L;
+
+
+// /* open skew quad file */
+// if ((fi = fopen(fic_deca,"r")) == NULL){
+// fprintf(stderr, "Error while opening file %s \n",fic_deca);
+// exit(1);
+// }
+
+// /* read decapole strength */
+// for (i = 0; i < globval.hcorr; i++){
+// fscanf(fi,"%le \n", &mKL[i]);
+// }
+// fclose(fi);
+
+// for (i = 0; i < globval.hcorr; i++){
+// if (trace) fprintf(stdout,"%le \n", mKL[i]);
+
+// 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"
+// ,i,
+// Cell.Elem.PName,Cell.Fnum, Cell.Knum,
+// Cell.Elem.M->PBpar[HOMmax+mOrder],
+// Cell.Elem.M->PBpar[HOMmax-mOrder]);
+// }
+// }
+
+/****************************************************************************/
+/* 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)
+ Purpose:
+ Compute momemtum acceptance along the ring, track the particle with
+ different energy, momentum acceptance is the energy when the particle
+ is lost or the last energy if the particle is not lost.
+
+ Based on the version in tracy 2.
+
+ Input:
+ MomAccFile file to save calculated momentum compact factor
+ deb first element for momentum acceptance,"debut" is beginning in French
+ fin last element for momentum acceptance,"fin" is end in French
+
+ ep_min minimum energy deviation for positive momentum acceptance
+ ep_max maximum energy deviation for positive momentum acceptance
+ nstepp number of energy steps for positive momentum acceptance
+
+ em_min minimum energy deviation for negative momentum acceptance
+ em_max maximum energy deviation for negative momentum acceptance
+ nstepm number of energy steps for negative momentum acceptance
+
+
+ * 1 grande section droite
+ * 13 entree premier bend
+ * 22 sortie SX4
+ * 41 section droite moyenne
+ * 173 fin superperiode
+
+ Output:
+ output file soleil.out : file of results
+ output file phase.out : file of tracking results during the process
+
+ Return:
+ none
+
+ Global variables:
+ none
+
+ specific functions:
+ set_vectorcod
+
+ Comments:
+ 30/06/03 add fflush(NULL) to force writing at the end to correct
+ unexpected bug: rarely the output file is not finished
+ 31/07/03 add closed orbit a element: useful for 6D tracking
+ delta_closed_orbite = dp(cavity)/2
+ 21/10/03 add array for vertical initial conditions using tracking
+ removed choice of tracking: now this should be done outside
+
+ 23/07/10 modify the call variable to the Cell_Pass( ): j-1L --> j (L3435, L3590)
+ since the Cell_Pass( ) is tracking from element i0 to i1(tracy 3), and
+ the Cell_Pass( ) is tracking from element i0+1L to i1(tracy 2).
+ 17/04/11 add number of turn
+ 27/06/11 fix the bug of the index in the tabz and tabpz 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.
+
+****************************************************************************/
+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 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;
+ Vector x0;
+ FILE *outf2, *outf1;
+
+ double **tabz0, **tabpz0;
+ struct tm *newtime; // for time
+ Vector codvector[Cell_nLocMax];
+ bool cavityflag, radiationflag;
+ bool trace=true;
+
+ x0.zero();
+
+ /* Get time and date */
+ newtime = GetTime();
+
+ /************************/
+ /* Fin des declarations */
+
+ /* File opening for writing */
+
+ outf1 = fopen("phase.out", "w");
+ outf2 = 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");
+
+
+ pos = deb; /* starting position or element index in the ring */
+
+ /***************************************************************/
+ fprintf(stdout,"Computing initial conditions ... \n");
+ /***************************************************************/
+
+ // cod search has to be done in 4D since in 6D it is zero
+ cavityflag = globval.Cavity_on;
+ radiationflag = globval.radiation;
+ globval.Cavity_on = false; /* Cavity on/off */
+ 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;
+ }
+
+ 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)
+ {
+ fprintf(stdout,"2 out of memory \n");
+ return;
+ }
+
+ // compute dP
+ if (nstepp != 1L)
+ dP = ep_max - (nstepp - i)*(ep_max - ep_min)/(nstepp - 1L);
+ else
+ dP = ep_max;
+
+ // find and store closed orbit for dP energy offset
+ 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];
+
+ 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
+ if (deb > 1L)
+ {
+ Cell_Pass(1L, deb - 1L, x0, lastpos); // track from 1 to deb-1L element
+ j = deb -1L;
+
+ if (lastpos != j)
+ { // look if stable
+ tabz0 [i- 1L][j] = 1.0;
+ tabpz0[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];
+ }
+ }
+ 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];
+ }
+
+ 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;
+ }
+ 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]);
+ }
+ }
+ }
+
+ globval.Cavity_on = cavityflag;
+ globval.radiation = radiationflag;
+
+ /***************************************************************/
+ fprintf(stdout,"Computing positive momentum acceptance ... \n");
+ /***************************************************************/
+
+ do
+ {
+ getcod(dP=0.0, lastpos); /* determine closed orbit */
+
+ getelem(pos,&Cell);
+ // coordinates around closed orbit which is non zero for 6D tracking
+ x = Cell.BeamPos[0];
+ px = Cell.BeamPos[1];
+ 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);
+
+ dp1 = 0.0;
+ dp2 = 0.0;
+ i = 0L;
+
+ do /* Tracking over nturn */
+ {
+ i++;
+ dp1 = dp2;
+
+ if (nstepp != 1L)
+ dp2= ep_max - (nstepp - i)*(ep_max - ep_min)/(nstepp - 1L);
+ 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]);
+
+ Trac(x, px, z+tabz0[i-1L][pos-1L], pz+tabpz0[i-1L][pos-1L], dp2+delta , ctau0, nturn, pos, lastn, lastpos, outf1);
+
+ }while (((lastn) == nturn) && (i != nstepp));
+
+ if ((lastn) == nturn)
+ dp1 = dp2;
+
+ 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++;
+
+ }while(pos != fin);
+
+ // free memory
+ for (i = 1L; i <= nstepp; i++){
+ free(tabz0 [i - 1L]);
+ free(tabpz0[i - 1L]);
+ }
+ free(tabz0);
+ free(tabpz0);
+
+ /***************************************************************/
+ /***************************************************************/
+ // NEGATIVE MOMENTUM ACCEPTANCE
+ /***************************************************************/
+ /***************************************************************/
+
+ fprintf(outf2,"\n"); /* A void line */
+
+ pos = deb; /* starting position in the ring */
+
+ /***************************************************************/
+ fprintf(stdout,"Computing initial conditions ... \n");
+ /***************************************************************/
+
+ // cod search has to be done in 4D since in 6D it is zero
+ cavityflag = globval.Cavity_on;
+ radiationflag = globval.radiation;
+ globval.Cavity_on = false; /* Cavity on/off */
+ 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){
+ 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){
+ fprintf(stdout,"2 out of memory \n"); return;
+ }
+
+ // compute dP
+ if (nstepm != 1L) {
+ dP = em_max - (nstepm - i)*(em_max - em_min)/(nstepm - 1L);
+ }
+ else {
+ dP = em_max;
+ }
+ // store closed orbit
+ 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];
+
+ // Store vertical initial conditions
+ // case where deb is not element 1
+ if (deb > 1L){
+ Cell_Pass(1L, deb - 1L, x0, lastpos); // track from 1 to deb-1L element
+ j = deb -1L;
+ if (lastpos != j){ // look if stable
+ tabz0 [i- 1L][j] = 1.0;
+ tabpz0[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];
+ }
+ }
+ 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]);
+ }
+
+ 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;
+ }
+ 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]);
+ }
+ }
+ }
-****************************************************************************/
-// void SetDecapole(void)
-// {
-// FILE *fi;
-// const char fic_deca[] ="/home/nadolski/soltracy/deca.dat";
-// int i;
-// double mKL[56]; /* array for skew quad tilt*/
-// CellType Cell;
-// long mOrder=5L;
+ globval.Cavity_on = cavityflag;
+ globval.radiation = radiationflag;
+ /***************************************************************/
+ fprintf(stdout,"Computing negative momentum acceptance ... \n");
+ /***************************************************************/
+
+ do {
+ getcod(dP=0.0, lastpos); /* determine closed orbit */
-// /* open skew quad file */
-// if ((fi = fopen(fic_deca,"r")) == NULL){
-// fprintf(stderr, "Error while opening file %s \n",fic_deca);
-// exit(1);
-// }
+ 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);
-// /* read decapole strength */
-// for (i = 0; i < globval.hcorr; i++){
-// fscanf(fi,"%le \n", &mKL[i]);
-// }
-// fclose(fi);
+ dp1 = 0.0;
+ dp2 = 0.0;
+ i = 0L;
+ do /* Tracking over nturn */
+ {
+ i++;
+ dp1 = dp2;
+ if (nstepm != 1L) {
+ dp2= em_max - (nstepm - i)*(em_max - em_min)/(nstepm - 1L);
+ }
+ 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);
+ }
+ while (((lastn) == nturn) && (i != nstepm));
-// for (i = 0; i < globval.hcorr; i++){
-// if (trace) fprintf(stdout,"%le \n", mKL[i]);
+ if ((lastn) == nturn) dp1 = dp2;
-// getelem(globval.hcorr_list[i], &Cell);
-// SetKLpar(Cell.Fnum, Cell.Knum, mOrder, mKL[i]);
+ 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);
+ pos++;
+ }
+ while(pos != fin);
-// if (trace) printf("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],
-// Cell.Elem.M->PBpar[HOMmax-mOrder]);
-// }
-// }
+ // free memory
+ for (i = 1L; i <= nstepp; i++){
+ free(tabz0 [i - 1L]);
+ free(tabpz0[i - 1L]);
+ }
+ free(tabz0);
+ free(tabpz0);
+
+ fflush(NULL); // force writing at the end (BUG??)
+ fclose(outf1);
+ fclose(outf2);
+}
/****************************************************************************/
-/* 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)
- 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)
+
+Purpose:
+ Parallel version of MomentumAcceptance( ).
Compute momemtum acceptance along the ring, track the particle with
different energy, momentum acceptance is the energy when the particle
is lost or the last energy if the particle is not lost.
@@ -2976,8 +3758,10 @@ void SetSkewQuad(void)
em_min minimum energy deviation for negative momentum acceptance
em_max maximum energy deviation for negative momentum acceptance
nstepm number of energy steps for negative momentum acceptance
-
-
+
+ numprocs number of processes used to do parallel computation.
+ myid process used to do parallel computation.
+
* 1 grande section droite
* 13 entree premier bend
* 22 sortie SX4
@@ -2998,28 +3782,12 @@ void SetSkewQuad(void)
set_vectorcod
Comments:
- 30/06/03 add fflush(NULL) to force writing at the end to correct
- unexpected bug: rarely the output file is not finished
- 31/07/03 add closed orbit a element: useful for 6D tracking
- delta_closed_orbite = dp(cavity)/2
- 21/10/03 add array for vertical initial conditions using tracking
- removed choice of tracking: now this should be done outside
-
- 23/07/10 modify the call variable to the Cell_Pass( ): j-1L --> j (L3435, L3590)
- since the Cell_Pass( ) is tracking from element i0 to i1(tracy 3), and
- the Cell_Pass( ) is tracking from element i0+1L to i1(tracy 2).
- 17/04/11 add number of turn
- 27/06/11 fix the bug of the index in the tabz and tabpz 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.
-
-****************************************************************************/
-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)
+ 19/07/11 Add feature to do parallel calculation of momentum compact factor.
+ Merged with the version written by Mao-Sen Qiu at Taiwan light source.
+****************************************************************************/
+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)
{
double dP = 0.0, dp1 = 0.0, dp2 = 0.0;
long lastpos = 0L,lastn = 0L;
@@ -3037,25 +3805,41 @@ void MomentumAcceptance(char *MomAccFile, long deb, long fin,
x0.zero();
- /* Get time and date */
+ // Get time and date
newtime = GetTime();
/************************/
- /* Fin des declarations */
+ // Fin des declarations
- /* File opening for writing */
+ // File opening for writing
+ char PhaseFile[50];
+ PhaseFile[0]='\0';
+ sprintf(PhaseFile,"%d",myid);
+ strcat(PhaseFile,"phase.out");
+ //printf("%s\n",PhaseFile);
+ outf1 = fopen(PhaseFile, "w");
- outf1 = fopen("phase.out", "w");
- outf2 = fopen(MomAccFile, "w");
+ if(myid==0)
+ {
+ fprintf(outf1,"# TRACY III -- %s \n", asctime2(newtime));
+ fprintf(outf1,"# i x xp z zp dp ctau\n#\n");
+ }
- fprintf(outf2,"# TRACY III -- %s \n", asctime2(newtime));
- fprintf(outf2,"# i s dp s_lost name_lost \n#\n");
+ char MomAccFile[50];
+ MomAccFile[0]='\0';
+ sprintf(MomAccFile,"%d",myid);
+ strcat(MomAccFile,_MomAccFile);
+ printf("%s\n",MomAccFile);
- fprintf(outf1,"# TRACY III -- %s \n", asctime2(newtime));
- fprintf(outf1,"# i x xp z zp dp ctau\n#\n");
-
+ outf2 = fopen(MomAccFile, "w");
- pos = deb; /* starting position or element index in the ring */
+ if(myid==0)
+ {
+ fprintf(outf2,"# TRACY III -- %s \n", asctime2(newtime));
+ fprintf(outf2,"# i s dp s_lost name_lost \n#\n");
+ }
+
+ // pos = deb; // starting position or element index in the ring
/***************************************************************/
fprintf(stdout,"Computing initial conditions ... \n");
@@ -3067,14 +3851,15 @@ void MomentumAcceptance(char *MomAccFile, long deb, long fin,
globval.Cavity_on = false; /* Cavity on/off */
globval.radiation = false; /* radiation on/off */
- // Allocation of an array of pointer array
+ // 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){
fprintf(stdout,"1 out of memory \n"); return;
}
- for (i = 1L; i <= nstepp; i++){ // loop over energy
+ 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));
@@ -3101,8 +3886,8 @@ void MomentumAcceptance(char *MomAccFile, long deb, long fin,
x0[4] = codvector[0][4];
x0[5] = codvector[0][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]);
+ 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
if (deb > 1L)
@@ -3131,7 +3916,7 @@ void MomentumAcceptance(char *MomAccFile, long deb, long fin,
for (j = deb; j < fin; j++)
{ // loop over elements
Cell_Pass(j, j, x0, lastpos);
- // Cell_Pass(j -1L, j, x0, lastpos);
+ //Cell_Pass(j -1L, j, x0, lastpos);
if (lastpos != j){ // look if stable
tabz0 [i - 1L][j] = 1.0;
@@ -3140,7 +3925,7 @@ void MomentumAcceptance(char *MomAccFile, long deb, long fin,
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]);
+// fprintf(stdout,"z0= % e pz0= % e\n", tabz0 [i - 1L][j], tabpz0 [i - 1L][j]);
}
}
}
@@ -3152,9 +3937,40 @@ void MomentumAcceptance(char *MomAccFile, long deb, long fin,
fprintf(stdout,"Computing positive momentum acceptance ... \n");
/***************************************************************/
- do
+//split tracking element region for each process
+//Eace core or process calculate different region of fmap according to id number. MSChiu 2011/10/13
+ int debN,finN;
+ int integer,residue;
+
+ //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);
+ exit_(1);
+}
+
+ 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);
+
+ //split tracking element region for each process
+ //the start element is from deb
+ if(myid<residue)
+ {
+ debN=myid*(integer+1) +deb;
+ finN=(myid+1)*(integer+1) +deb;
+ }
+ else
+ {
+ debN=residue*(integer+1)+(myid-residue)*integer +deb;
+ finN=residue*(integer+1)+(myid+1-residue)*integer +deb;
+ }
+
+
+ //do
+ for(pos=debN;pos<finN;pos++)
{
- getcod(dP=0.0, lastpos); /* determine closed orbit */
+ getcod(dP=0.0, lastpos); // determine closed orbit
getelem(pos,&Cell);
// coordinates around closed orbit which is non zero for 6D tracking
@@ -3164,14 +3980,17 @@ void MomentumAcceptance(char *MomAccFile, long deb, long fin,
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);
+
+ 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, globval.CODvect[0], globval.CODvect[1],globval.CODvect[2],
+ globval.CODvect[3], globval.CODvect[4],globval.CODvect[5]);
+
dp1 = 0.0;
dp2 = 0.0;
i = 0L;
- do /* Tracking over nturn */
+ do // Tracking over nturn
{
i++;
dp1 = dp2;
@@ -3181,7 +4000,10 @@ 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 (trace) printf("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);
+
if (0) fprintf(stdout,"pos=%4ld z0 =% 10.5f pz0 =% 10.5f \n", pos, tabz0[i-1L][pos-1L], tabpz0[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);
@@ -3195,15 +4017,16 @@ void MomentumAcceptance(char *MomAccFile, long deb, long fin,
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);
+ 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++;
+ // pos++;
- }while(pos != fin);
+ }//while(pos != fin);
// free memory
- for (i = 1L; i <= nstepp; i++){
+ for (i = 1L; i <= nstepp; i++)
+ {
free(tabz0 [i - 1L]);
free(tabpz0[i - 1L]);
}
@@ -3216,9 +4039,9 @@ void MomentumAcceptance(char *MomAccFile, long deb, long fin,
/***************************************************************/
/***************************************************************/
- fprintf(outf2,"\n"); /* A void line */
+ fprintf(outf2,"Negative\n"); // A void line
- pos = deb; /* starting position in the ring */
+ // pos = deb; // starting position in the ring
/***************************************************************/
fprintf(stdout,"Computing initial conditions ... \n");
@@ -3227,8 +4050,8 @@ void MomentumAcceptance(char *MomAccFile, long deb, long fin,
// cod search has to be done in 4D since in 6D it is zero
cavityflag = globval.Cavity_on;
radiationflag = globval.radiation;
- globval.Cavity_on = false; /* Cavity on/off */
- globval.radiation = false; /* radiation on/off */
+ globval.Cavity_on = false; // Cavity on/off
+ globval.radiation = false; // radiation on/off
// Allocation of an array of pointer array
tabz0 = (double **)malloc((nstepm)*sizeof(double*));
@@ -3281,7 +4104,7 @@ void MomentumAcceptance(char *MomAccFile, long deb, long fin,
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]);
+// fprintf(stdout,"z0= % e pz0= % e\n", tabz0 [i - 1L][j], tabpz0 [i - 1L][j]);
}
for (j = deb; j < fin; j++){ // loop over elements
@@ -3306,10 +4129,12 @@ void MomentumAcceptance(char *MomAccFile, long deb, long fin,
fprintf(stdout,"Computing negative momentum acceptance ... \n");
/***************************************************************/
- do {
- getcod(dP=0.0, lastpos); /* determine closed orbit */
+ // do
+ for(pos=debN;pos<finN;pos++)
+{
+ getcod(dP=0.0, lastpos); // determine closed orbit
- getelem(pos,&Cell);
+ getelem(pos,&Cell);
// coordinates around closed orbit which is non zero for 6D tracking
x = Cell.BeamPos[0];
px = Cell.BeamPos[1];
@@ -3317,13 +4142,12 @@ void MomentumAcceptance(char *MomAccFile, long deb, long fin,
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);
+ fprintf(stdout,"%3ld %6.4g %6.4g %6.4g %6.4g %6.4g %6.4g\n", pos, x, px, z, pz, delta, ctau0);
dp1 = 0.0;
dp2 = 0.0;
i = 0L;
- do /* Tracking over nturn */
+ do // Tracking over nturn
{
i++;
dp1 = dp2;
@@ -3333,7 +4157,8 @@ 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);
+ 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);
}
while (((lastn) == nturn) && (i != nstepm));
@@ -3344,14 +4169,15 @@ void MomentumAcceptance(char *MomAccFile, long deb, long fin,
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);
- pos++;
+ 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++;
}
- while(pos != fin);
+ //while(pos != fin);
// free memory
- for (i = 1L; i <= nstepp; i++){
+ for (i = 1L; i <= nstepp; i++)
+ {
free(tabz0 [i - 1L]);
free(tabpz0[i - 1L]);
}
@@ -3362,7 +4188,8 @@ void MomentumAcceptance(char *MomAccFile, long deb, long fin,
fclose(outf1);
fclose(outf2);
}
-
+
+
/****************************************************************************/
/* set_vectorcod(double codvector[Cell_nLocMax][6], double dP)
@@ -5739,4 +6566,158 @@ void PrintTrack(const char *TrackFile, double x, double px, double y,double py,
fclose(outf);
}
+/*******************************************************************************************************
+ void ReadVirtualSkewQuad(const char *SkewQuadFile)
+
+ Purpose:
+ Set sources of coupling on SOLEIL storage ring, to simulate localized skew gradient.
+ Then calculate the coupling.
+
+ Comments:
+ Based on the tracy-2.7 file VirtualSkewQuad(void).
+
+ 21/12/2011 Jianfeng Zhang@soleil
+ Fix the bug to call Elem_Getpos() in
+ SetKLpar(ElemIndex("SQ"),i+1, mOrder=2L, corr_strength),
+ the start kid index of the family starts from 1 !!!!!
+
+*****************************************************************************************************/
+ 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 */
+ double conv = 0.0, brho = 0.0, corr_strength =0.0;
+
+ brho = globval.Energy/0.299792458; /* magnetic rigidity */
+ conv = 93.88e-4; /*conversion des A en T*/
+
+ nqtcorr = GetnKid(ElemIndex("SQ"));
+ printf("Number of virtual skew quadrupoles: %d\n",nqtcorr);
+
+ /* open virtual QT corrector file */
+ if ((fi = fopen(VirtualSkewQuadFile,"r")) == NULL)
+ {
+ fprintf(stderr, "Error while opening file %s \n",VirtualSkewQuadFile);
+ exit(1);
+ }
+
+ for (i = 1; i <= nqtcorr; i++){
+ fscanf(fi,"%le \n", &qtcorr[i-1]);
+ corr_strength = qtcorr[i-1]*conv/brho;
+ //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);
+ }
+ fclose(fi);
+}
+/*********************************************************************************
+void CODCorrect(FILE *hcorr_file,FILE *vcorr_file,int n_orbit,int nwh,int nwv)
+
+ Purpose:
+ Get the response matrix between bpms and h/v correctors; then call CorrectCODs( )
+ to do orbit correction; finally saved the summary of the orbit correction.
+
+
+
+ Comments:
+ Written by Jianfeng Zhang@soleil 20/12/2011
+
+*********************************************************************************/
+void CODCorrect(const char *hcorr_file, const char *vcorr_file,int n_orbit,int nwh,int nwv)
+{
+ bool cod = false;
+ int k=0;
+ FILE *hOrbitFile, *vOrbitFile, *OrbScanFile;
+ int hcorrIdx[nCOR], vcorrIdx[nCOR]; //list of corr for orbit correction
+
+ //initialize the corrector list
+ for ( k = 0; k < nCOR; k++){
+ hcorrIdx[k] = -1;
+ vcorrIdx[k] = -1;
+ }
+
+
+ //Get response matrix between bpm and correctors, and then print the SVD setting to the files
+ // select correctors to be used
+ readCorrectorList(hcorr_file, vcorr_file, hcorrIdx, vcorrIdx);
+
+ fprintf(stdout, "\n\nSVD correction setting:\n");
+ fprintf(stdout, "H-plane %d singular values:\n", nwh);
+ fprintf(stdout, "V-plane %d singular values:\n\n",nwv);
+
+ // compute beam response matrix
+ printf("\n");
+ printf("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);
+ /* gcmat(globval.bpm, globval.hcorr, 1);
+ gcmat(globval.bpm, globval.vcorr, 2);*/
+
+ // print response matrices to files 'svdh.out' and file 'svdv.out'
+ prt_gcmat(globval.bpm, globval.hcorr, 1);
+ prt_gcmat(globval.bpm, globval.vcorr, 2);
+
+ //print the statistics of orbit in file 'OrbScanFile.out'
+ OrbScanFile = file_write("OrbScanFile.out");
+
+ //write files with orbits at all element locations
+ hOrbitFile = file_write("horbit.out");
+ vOrbitFile = file_write("vorbit.out");
+
+ fprintf(hOrbitFile, "# First line: s-location (m) \n");
+ fprintf(hOrbitFile, "# After orbit correction: Horizontal closed orbit at all element locations (with %3d BPMs) at different loop\n", GetnKid(globval.bpm));
+ fprintf(vOrbitFile, "# First line s-location (m) \n");
+ fprintf(vOrbitFile, "# After orbit correction: Vertical closed orbit at all element locations (with %3d BPMs) at different loop\n", GetnKid(globval.bpm));
+
+ for (k = 0; k < globval.Cell_nLoc; k++){
+ fprintf(hOrbitFile, "% 9.3e ", Cell[k].S);
+ fprintf(vOrbitFile, "% 9.3e ", Cell[k].S);
+ } // end for
+
+ fprintf(hOrbitFile, "\n");
+ fprintf(vOrbitFile, "\n");
+
+
+ //prepare for the orbit correction
+ // Clear trim setpoints
+ set_bnL_design_fam(globval.hcorr, Dip, 0.0, 0.0);
+ set_bnL_design_fam(globval.vcorr, Dip, 0.0, 0.0);
+
+ // Beam based alignment
+ if (bba) {
+ Align_BPM2quad(Quad);
+ }
+
+ //orbit correction
+ cod = CorrectCODs(hOrbitFile, vOrbitFile, OrbScanFile, n_orbit, nwh, nwv, hcorrIdx, vcorrIdx);
+
+ /* cod = CorrectCOD_N(ae_file, n_orbit, n_scale, k);*/
+ printf("\n");
+
+ if (cod){
+ /* printf("done with orbit correction, now do coupling",
+ " correction plus vert. disp\n");*/
+
+ printf("Orbit correction succeeded\n");
+ }else{
+ fprintf(stdout, "!!! Orbit correction failed\n");
+ exit_(1);
+ }
+ //chk_cod(cod, "iter # %3d error_and_correction");
+
+ // for debugging
+ //print flat lattice
+ //sprintf(mfile_name, "flat_file.%03d.dat",k);
+ //prtmfile(mfile_name);
+
+ // close file giving orbit at BPM location
+ fclose(hOrbitFile);
+ fclose(vOrbitFile);
+ fclose(OrbScanFile);
+}
--
GitLab