/* Current revision $Revision$
On branch $Name$
Latest change $Date$ by $Author$
*/
#define ORDER 1
int no_tps = ORDER; // arbitrary TPSA order is defined locally
extern bool freq_map;
#include "tracy_lib.h"
void compare_cod(void)
{
const double dPmin = 1e-10, dPmax = 1e-6;
int k;
const int kmax = 20;
double dp, dpstep = dPmax/kmax;
FILE * outf;
const char fic[] = "compare_cod.out";
long lastpos = 0L;
Vector vector_findcod;
/* Opening file */
if ((outf = fopen(fic, "w")) == NULL) {
fprintf(stdout, "compare_cod: error while opening file %s\n", fic);
exit_(1);
}
fprintf(stdout,"\n");
for (k = 0; k <= kmax; k++){
dp = dPmin + k*dpstep;
getcod(dp, lastpos); // get cod for printout
CopyVec(6, globval.CODvect, vector_findcod);
findcod(dp);
// fprintf(stdout,"%14.6e %14.6e %14.6e %14.6e %14.6e %14.6e %14.6e %14.6e %14.6e\n",
// dp, globval.CODvect[0], vector_findcod[0], globval.CODvect[1],vector_findcod[1],
// globval.CODvect[2], vector_findcod[2], globval.CODvect[3], vector_findcod[3]);
fprintf(stdout,"%14.6e %14.6e %14.6e %14.6e %14.6e %14.6e %14.6e %14.6e %14.6e\n",
dp,
globval.CODvect[0], (globval.CODvect[0]-vector_findcod[0])/globval.CODvect[0],
globval.CODvect[1], (globval.CODvect[1]-vector_findcod[1])/globval.CODvect[1],
globval.CODvect[2], (globval.CODvect[2]-vector_findcod[2])/globval.CODvect[2],
globval.CODvect[3], (globval.CODvect[3]-vector_findcod[3])/globval.CODvect[3]);
}
fclose(outf);
}
#define LOG10 log(10.0)
void Getchrom2(double dP)
{
long lastpos = 0L;
double dPlocal = 0.0, expo = 0.0, TEMP = 0.0, TEMPX = 0.0, TEMPZ = 0.0;
Vector2 alpha={0.0,0.0}, beta={0.0,0.0}, gamma={0.0,0.0}, nu={0.0,0.0},
nu0={0.0,0.0}, Chrom={0.0,0.0};
trace = false;
if (dP != 0.0) {
fprintf(stderr,"Ring_Getchrom: Warning this is NOT the CHROMA, dP=%e\n",dP);
}
/* initialization */
globval.Chrom[0] = 1e38;
globval.Chrom[1] = 1e38;
expo = log(globval.dPcommon) / LOG10;
do
{
Chrom[0] = globval.Chrom[0];
Chrom[1] = globval.Chrom[1];
dPlocal = exp(expo*LOG10);
/* Get cod for energy dP - dPlocal*/
GetCOD(globval.CODimax, globval.CODeps, dP - dPlocal*0.5,
lastpos);
if (!status.codflag)
{ /* if no cod */
fprintf(stderr,"Ring_Getchrom: Lattice is unstable for dP - dPlocal=% .5e\n",
dP - dPlocal*0.5);
return;
}
/* get tunes for energy dP - dPlocal/2 from oneturn matrix */
Cell_GetABGN(globval.OneTurnMat, alpha, beta, gamma, nu0);
/* Get cod for energy dP + dPlocal*/
GetCOD(globval.CODimax, globval.CODeps, dP + dPlocal*0.5,
lastpos);
if (!status.codflag) { /* if no cod */
fprintf(stderr,"Ring_Getchrom Lattice is unstable for dP + dPlocal=% .5e \n",
dP + dPlocal*0.5);
return;
}
/* get tunes for energy dP + dPlocal/2 from oneturn matrix */
Cell_GetABGN(globval.OneTurnMat, alpha, beta, gamma, nu);
if (!globval.stable) {
printf("Ring_Getchrom: Lattice is unstable\n");
}
/* Get chromaticities by numerical differentiation*/
globval.Chrom[0] = (nu[0] - nu0[0]) / dPlocal;
globval.Chrom[1] = (nu[1] - nu0[1]) / dPlocal;
TEMP=sqrt(
fabs(globval.Chrom[0]*globval.Chrom[0]- Chrom[0]*Chrom[0])
+fabs(globval.Chrom[1]*globval.Chrom[1]- Chrom[1]*Chrom[1])
);
TEMPX = sqrt(fabs(globval.Chrom[0]*globval.Chrom[0]- Chrom[0]*Chrom[0]));
TEMPZ = sqrt(fabs(globval.Chrom[1]*globval.Chrom[1]- Chrom[1]*Chrom[1]));
// TEST CHROMA convergence
if (trace) {
fprintf(stdout,"\nexpo % e xix = % e xiz = % e TEMP = %e TEMPX %+e TEMPZ %+e\n",
expo,Chrom[0],Chrom[1],TEMP, TEMPX, TEMPZ);
fprintf(stdout,"expo % e nux = % e nuz = % e dPlocal= %+e\n",
expo,nu0[0],nu0[1],dP-0.5*dPlocal);
fprintf(stdout,"expo % e nux = % e nuz = % e dPlocal= %+e\n",
expo,nu[0],nu[1],dP+0.5*dPlocal);
}
fprintf(stdout, "%+e %+.12e %+.12e %+.12e %+.12e\n", dPlocal,
globval.Chrom[0], fabs(globval.Chrom[0]-Chrom[0])/Chrom[0],
globval.Chrom[1], fabs(globval.Chrom[1]-Chrom[1])/Chrom[1]);
expo += 0.1;
} while (expo<-2);
status.chromflag = true;
}
//***************************************************************************************
//
// MAIN CODE
//
//****************************************************************************************
int main(int argc, char *argv[])
{
const long seed = 1121;
iniranf(seed); setrancut(2.0);
// turn on globval.Cavity_on and globval.radiation to get proper synchr radiation damping
// IDs accounted too if: wiggler model and symplectic integrator (method = 1)
globval.H_exact = false;
globval.quad_fringe = false; // quadrupole fringe field
globval.radiation = false; // synchrotron radiation
globval.emittance = false; // emittance
globval.pathlength = false;
// globval.bpm = 0;
// const double x_max_FMA = 10e-3, delta_FMA = 10e-2;
// const int n_x = 801, n_dp = 80, n_tr = 2048;
double nux=0, nuz=0, ksix=0, ksiz=0;
bool chroma;
double dP = 0.0;
long lastpos = -1L;
char str1[S_SIZE];
// for time handling
uint32_t start, stop;
/************************************************************************
start read in files and flags
*************************************************************************/
if (argc > 1){
read_script(argv[1], true);}
else{
fprintf(stdout, "Not enough parameters\nSyntax is program parameterfile\n");
return 1;
}
/************************************************************************
end read in files and flags
*************************************************************************/
prtmfile("flat_file.dat"); // writes flat file /* very important file for debug*/
printlatt(); /* SOLEIL print out lattice functions */
printglob();
double V_RF;
V_RF = get_RFVoltage(ElemIndex("cav"));
set_RFVoltage(ElemIndex("cav"), 3e6); // 3 MV
V_RF = get_RFVoltage(ElemIndex("cav"));
// Flag factory
// Chamber factory
if (ChamberFlag == false)
ChamberOff(); // turn off vacuum chamber setting, use the default one
else if (ChamberNoU20Flag == true)
DefineChNoU20(); // using vacuum chamber setting but without undulator U20
else if (ReadChamberFlag == true)
ReadCh(chamber_file); /* read vacuum chamber from a file "Apertures.dat" , soleil version*/
PrintCh(); // print chamber into chamber.out
//get_matching_params_scl(); // get tunes and beta functions at entrance
get_alphac2(); // compute up to 3rd order mcf
//cout << endl << "computing tune shifts" << endl;
//dnu_dA(10e-3, 5e-3, 0.002);
//get_ksi2(0.0); // this gets the chromas and writes them into chrom2.out
// compute tunes by tracking (should be the same as by DA)
if (TuneTracFlag == true) {
GetTuneTrac(1026L, 0.0, &nux, &nuz);
fprintf(stdout,"From tracking: nux = % f nuz = % f \n",nux,nuz);
}
// compute chromaticities by tracking (should be the same as by DA)
if (ChromTracFlag == true){
start = stampstart();
GetChromTrac(2L, 1026L, 1e-5, &ksix, &ksiz);
stop = stampstop(start);
fprintf(stdout,"From tracking: ksix= % f ksiz= % f \n",ksix,ksiz);
}
if (FitTuneFlag == true){
fprintf(stdout, "\n Fitting tunes\n");
FitTune(ElemIndex("qp7"),ElemIndex("qp9"), targetnux, targetnuz);
Ring_GetTwiss(chroma=true, 0.0); /* Compute and get Twiss parameters */
printglob(); /* print parameter list */
}
if (FitChromFlag == true){
fprintf(stdout, "\n Fitting chromaticities\n");
FitChrom(ElemIndex("sx9"),ElemIndex("sx10"), targetksix, targetksiz);
Ring_GetTwiss(chroma=true, 0.0); /* Compute and get Twiss parameters */
printglob(); /* print parameter list */
}
// coupling calculation
if (CouplingFlag == true){
Ring_GetTwiss(chroma=true, 0.0); /* Compute and get Twiss parameters */
printlatt(); /* dump linear lattice functions into "linlat.dat" */
GetEmittance(ElemIndex("cav"), true);
Coupling_Edwards_Teng();
printglob(); /* print parameter list */
}
// add coupling by random rotating of the quadrupoles
if (ErrorCouplingFlag == true){
SetErr();
Ring_GetTwiss(chroma=true, 0.0); /* Compute and get Twiss parameters */
printlatt(); /* dump linear lattice functions into "linlat.dat" */
Coupling_Edwards_Teng();
GetEmittance(ElemIndex("cav"), true);
printglob(); /* print parameter list */
}
// WARNING Fit tunes and chromaticities before applying errors !!!!
//set multipoles in all magnets
if (MultipoleFlag == true ){
if (ThinsextFlag ==true){
fprintf(stdout, "\n Setting Multipoles for lattice with thin sextupoles \n");
Multipole_thinsext(fic_hcorr,fic_vcorr,fic_skew); /* for thin sextupoles */
Ring_GetTwiss(chroma=true, 0.0); /* Compute and get Twiss parameters */
printglob();
}
else{
fprintf(stdout, "\n Setting Multipoles for lattice with thick sextupoles \n");
Multipole_thicksext(fic_hcorr,fic_vcorr,fic_skew); /* for thick sextupoles */
Ring_GetTwiss(chroma=true, 0.0); /* Compute and get Twiss parameters */
printglob();
}
}
/******************************************************************************************/
// COMPUTATION PART after setting the model
/******************************************************************************************/
//first print the full lattice with error as a flat file
prtmfile("flat_file_error.dat"); // writes flat file /* very important file for debug*/
// computes TuneShift with amplitudes
if (AmplitudeTuneShiftFlag == true){
if (ChamberFlag == true ){
TunesShiftWithAmplitude(_AmplitudeTuneShift_nxpoint,
_AmplitudeTuneShift_nypoint, _AmplitudeTuneShift_nturn,
_AmplitudeTuneShift_xmax, _AmplitudeTuneShift_ymax,
_AmplitudeTuneShift_delta);
//NuDx(31L,21L,516L,0.025,0.005,dP);
}
else{ // Utility ?
TunesShiftWithAmplitude(50L,30L,516L,0.035,0.02,dP);
}
}
if (EnergyTuneShiftFlag == true){
if (ChamberFlag == true ){
TunesShiftWithEnergy(_EnergyTuneShift_npoint,
_EnergyTuneShift_nturn, _EnergyTuneShift_deltamax);
//NuDp(31L,1026L,0.06);
}
else{ // utility ?
TunesShiftWithEnergy(31L,1026L,0.06);
}
}
// Computes FMA
if (FmapFlag == true){
if (ChamberFlag == true ){
if (ExperimentFMAFlag == true)
fmap( _FmapFlag_nxpoint, _FmapFlag_nypoint,
_FmapFlag_nturn, _FmapFlag_xmax, _FmapFlag_ymax,
_FmapFlag_delta, _FmapFlag_diffusion);
//fmap(40,12,258,-20e-3,5e-3,0.0,true); // for experimental
if (DetailedFMAFlag == true)
fmap(100,50,1026,20e-3,5e-3,0.0,true);
}
else{ // Utility
if (ExperimentFMAFlag == true)
fmap(40,12,258,-32e-3,5e-3,0.0,true);
if (DetailedFMAFlag == true)
fmap(200,100,1026,32e-3,7e-3,0.0,true);
}
}
if (CodeComparaisonFlag){
fmap(200,100,1026,-32e-3,7e-3,0.0,true);
}
// MOMENTUM ACCEPTANCE
if (MomentumAccFlag == true){
MomentumAcceptance(
_MomentumAccFlag_istart, _MomentumAccFlag_istop,
_MomentumAccFlag_deltaminp, _MomentumAccFlag_deltamaxp, _MomentumAccFlag_nstepp,
_MomentumAccFlag_deltaminn, _MomentumAccFlag_deltamaxn, _MomentumAccFlag_nstepn);
// MomentumAcceptance(1L, 108L, 0.01, 0.05, 100L, -0.01, -0.05, 100L);
}
// induced amplitude
if (InducedAmplitudeFlag == true){
InducedAmplitude(193L);
}
if (EtaFlag == true){
// compute cod and twiss parameters for different energy offsets
for (int ii=0; ii<=40; ii++) {
dP = -0.02+ 0.001*ii;
Ring_GetTwiss(chroma=false, dP); /* Compute and get Twiss parameters */
printlatt(); /* dump linear lattice functions into "linlat.dat" */
getcod (dP, lastpos);
// printcod();
prt_cod("cod.out", globval.bpm, true);
//system("mv linlat.out linlat_ooo.out");
sprintf(str1, "mv cod.out cod_%02d.out", ii);
system(str1);
sprintf(str1, "mv linlat.out linlat_%02d.out", ii);
system(str1);
}
}
if (PhaseSpaceFlag == true){
start = stampstart();
Phase(0.001,0,0.001,0,0.001,0, 1);
printf("the simulation time for phase space in tracy 3 is \n");
stop = stampstop(start);
}
//compare_cod();
Getchrom2(0.0);
return 0;
/*********************************************************************************
Delicated for max4 lattice. To load alignment error files and do correction
----------------------------------------------------------------------------------
*********************************************************************************/
if (false) {
//prt_cod("cod.out", globval.bpm, true); //prints a specific closed orbit with corrector strengths
globval.bpm = ElemIndex("bpm_m"); //definition for max4 lattice, bpm
// globval.bpm = ElemIndex("bpm");
globval.hcorr = ElemIndex("corr_h"); globval.vcorr = ElemIndex("corr_v"); //definition for max4 lattice, correctors
// globval.hcorr = ElemIndex("ch"); globval.vcorr = ElemIndex("cv");
globval.gs = ElemIndex("GS"); globval.ge = ElemIndex("GE"); //definition for max4 lattice, girder maker
//compute response matrix (needed for OCO)
gcmat(globval.bpm, globval.hcorr, 1); gcmat(globval.bpm, globval.vcorr, 2);
//print response matrix (routine in lsoc.cc)
//prt_gcmat(globval.bpm, globval.hcorr, 1); prt_gcmat(globval.bpm, globval.vcorr, 2);
//gets response matrix, does svd, evaluates correction for N seed orbits
//get_cod_rms_scl(dx, dy, dr, n_seed)
//get_cod_rms_scl(100e-6, 100e-6, 1.0e-3, 100); //trim coils aren't reset when finished
//for alignments errors check LoadAlignTol (in nsls_ii_lib.cc) and AlignErr.dat
//CorrectCOD_N("/home/simon/projects/src/lattice/AlignErr.dat", 3, 1, 1);
//for field errors check LoadFieldErr(in nsls_ii_lib.cc) and FieldErr.dat
//LoadFieldErr("/home/simon/projects/src/lattice/FieldErr.dat", true, 1.0, true);
//for alignments errors check LoadAlignTol (in nsls_ii_lib.cc) and AlignErr.dat
//LoadAlignTol("/home/simon/projects/src/lattice/AlignErr.dat", true, 1.0, true, 1);
//LoadAlignTol("/home/simon/projects/out/20091126/AlignErr.dat", true, 1.0, true, 1);
//prt_cod("cod_err.out", globval.bpm, true); //prints a specific closed orbit with corrector strengths
// delicated for max4 lattice
//load alignment errors and field errors, correct orbit, repeat N times, and get statistics
get_cod_rms_scl_new(1); //trim coils aren't reset when finished
//for aperture limitations use LoadApers (in nsls_ii_lib.cc) and Apertures.dat
//globval.Aperture_on = true;
//LoadApers("/home/simon/projects/src/lattice/Apertures.dat", 1, 1);
}
/*******************************************************************************
Call nsls-ii_lib.cc
-----------------------------------------------------------------------------
*******************************************************************************/
//
// tune shift with amplitude
double delta = 0;
if (false) {
cout << endl;
cout << "computing tune shifts" << endl;
dnu_dA(25e-3, 5e-3, 0.0);
get_ksi2(delta); // this gets the chromas and writes them into chrom2.out
// get_ksi2(5.0e-2); // this gets the chromas and writes them into chrom2.out
}
if (false) {
//fmap(n_x, n_y, n_tr, x_max_FMA, y_max_FMA, 0.0, true, false);
// fmapdp(n_x, n_dp, n_tr, x_max_FMA, -delta_FMA, 1e-3, true, false); // always use -delta_FMA (+delta_FMA appears broken)
fmapdp(n_x, n_dp, n_tr, x_max_FMA, -delta_FMA, 1e-3, true); // always use -delta_FMA (+delta_FMA appears broken)
} else {
globval.Cavity_on = true; // this gives longitudinal motion
globval.radiation = false; // this adds ripple around long. ellipse (needs many turns to resolve damp.)
//globval.Aperture_on = true;
//LoadApers("/home/simon/projects/src/lattice/Apertures.dat", 1, 1);
//get_dynap_scl(delta, 512);
}
}