Select Git revision
soltracy.cc
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
soltracy.cc 28.35 KiB
/*
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"
//***************************************************************************************
//
// MAIN CODE
//
//****************************************************************************************
int main(int argc, char *argv[]) {
/* for time handling */
uint32_t start, stop;
/*set the random value for random generator*/
const long seed = 1121; //the default random seed number
iniranf(seed); //initialize the seed
setrancut(2.0); //default value of the normal cut for the normal distribution
// turn on globval.Cavity_on and globval.radiation to get proper synchro radiation damping
// IDs accounted too if: wiggler model and symplectic integrator (method = 1)
globval.H_exact = false;
/* parameters to read the user input script .prm */
long i=0L; //initialize the for loop to read command string
char CommandStr[max_str];
double nux = 0.0, nuy = 0.0, ksix = 0.0, ksiy = 0.0;
bool chroma=true;
double dP = 0.0;
long lastpos = -1L;
char str1[S_SIZE];
// paramters to read the command line in user input script
long CommandNo = -1L; //the number of commands, since this value is also
// the index of array UserCommandFlag[], so this
// value is always less than 1 in the really case.
const int NCOMMAND = 500;
UserCommand UserCommandFlag[NCOMMAND];
/************************************************************************
read in files and flags
*************************************************************************/
if (argc > 1) {
read_script(argv[1], true,CommandNo, UserCommandFlag);
} else {
fprintf(stdout, "Not enough parameters\nSyntax is program parameter file\n");
exit_(1);
}
/* get the family index of the special elements, prepare for printglob()*/
if(strcmp(bpm_name,"")==0)
globval.bpm = 0;
else
globval.bpm = ElemIndex(bpm_name);
if(strcmp(skew_quad_name,"")==0)
globval.qt = 0;
else
globval.qt = ElemIndex(skew_quad_name); if(strcmp(hcorr_name,"")==0)
globval.hcorr = 0;
else
globval.hcorr = ElemIndex(hcorr_name);
if(strcmp(vcorr_name,"")==0)
globval.vcorr = 0;
else
globval.vcorr = ElemIndex(vcorr_name);
if(strcmp(gs_name,"")==0)
globval.gs = 0;
else
globval.gs = ElemIndex(gs_name);
if(strcmp(ge_name,"")==0)
globval.ge = 0;
else
globval.ge = ElemIndex(ge_name);
// globval.g = ElemIndex("g"); /* get family index of girder*/
/* print the summary of the element in lattice */
printglob();
/************************************************************************
print files, very important file for debug
*************************************************************************/
//print flat file with all the design values of the lattice,
prtmfile("flat_file.dat");
// print location, twiss parameters and close orbit at all elements position to a file
getcod(dP, lastpos);
prt_cod("cod.out", globval.bpm, true);
//get_matching_params_scl(); // get tunes and beta functions at entrance
// compute up to 3rd order momentum compact factor
get_alphac2();
//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
/*********************************************************************
Flag factory
*********************************************************************/
for(i=0; i<=CommandNo; i++){//read user defined command by sequence
//assign user defined command
strcpy(CommandStr,UserCommandFlag[i].CommandStr);
//turn on flag for quadrupole fringe field
if(strcmp(CommandStr,"QuadFringeOnFlag") == 0) {
globval.quad_fringe = true;
cout << "\n";
cout << "globval.quad_fringe is " << globval.quad_fringe << "\n";
}
//set RF voltage
else if(strcmp(CommandStr,"RFvoltageFlag") == 0) {
printf("\nSetting RF voltage:\n");
printf(" Old RF voltage is: %lf [MV]\n", get_RFVoltage(globval.cav)
/ 1e6);
set_RFVoltage(globval.cav, UserCommandFlag[i].RFvolt);
printf(" New RF voltage is: %lf [MV]\n", get_RFVoltage(globval.cav)
/ 1e6);
}
// Chamber factory
else if(strcmp(CommandStr,"ReadChamberFlag") == 0) { ReadCh(UserCommandFlag[i].chamber_file); /* read vacuum chamber from a file "Apertures.dat" , soleil version*/
PrintCh(); // print chamber into chamber.out
}
// read the misalignment errors to the elements, then do COD correction
// using SVD method.
// Based on the function error_and_correction() in nsls-ii_lib_templ.h
// else if(strcmp(ReadaefileFlag, "") != 0){
else if (strcmp(CommandStr,"ReadaefileFlag") == 0) {
// ***** Apply corrections and output flatfile for n_stat sets of random #'s
bool cod = true;
int k, icod=0;
FILE *hOrbitFile, *vOrbitFile ;
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
if (n_orbit!=0 || n_scale!=0) {
// 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(OrbScanFileName);
//write files with orbits at all element locations
hOrbitFile = file_write(hOrbitFileName);
vOrbitFile = file_write(vOrbitFileName);
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");
for (k = 1; k <= n_stat; k++) {// loop for n_stat sets of errors
if (n_orbit!=0 || n_scale!=0) {
cod =
CorrectCOD_Ns(hOrbitFile, vOrbitFile, UserCommandFlag[i].ae_file,
n_orbit,n_scale,k,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");*/
if(n_orbit == 0)
printf("iter # %3d n_obit = 0, no orbit correction \n",k);
else
printf("iter # %3d Orbit correction succeeded\n", k);
}
else
if(!cod){
icod = icod + 1;
fprintf(stdout, "!!! iter # %3d error_and_correction\n",k);
}
//chk_cod(cod, "iter # %3d error_and_correction");
}
// Ring_GetTwiss(chroma, dp);
Ring_GetTwiss(true, 0.0);
// for debugging
//print flat lattice
//sprintf(mfile_name, "flat_file.%03d.dat",k);
//prtmfile(mfile_name);
}
fprintf(stdout, "Number of unstable orbits %d/%d", icod, n_stat);
prt_cod("corr_after.out", globval.bpm, true);
// close file giving orbit at BPM location
fclose(hOrbitFile);
fclose(vOrbitFile);
fclose(OrbScanFile);
}
// set the field error into the lattice
// The corresponding field error is replaced by the new value.
// This feature is generic, works for all lattices
else if (strcmp(CommandStr,"ReadfefileFlag") == 0) {
fprintf(stdout,"\n Read field error from fe_file: \n");
LoadFieldErrs(UserCommandFlag[i].fe_file, true, 1.0, true, 1);
Ring_GetTwiss(true, 0.0);
prtmfile("flat_file_fefile.dat");
}
// read multipole errors from a file; specific for soleil lattice
else if(strcmp(CommandStr,"ReadMultipoleFlag") == 0) {
fprintf(stdout,"\n Read Multipoles file for lattice with thick sextupoles, specific for SOLEIL lattice: \n");
ReadFieldErr(multipole_file);
//first print the full lattice with error as a flat file
prtmfile("flat_file_errmultipole.dat"); // writes flat file with all element errors /* very important file for debug*/
Ring_GetTwiss(chroma = true, 0.0); /* Compute and get Twiss parameters */
printglob();
}
//print the twiss paramaters in a file defined by the name
else if(strcmp(CommandStr,"PrintTwissFlag") == 0) {
cout << "\n";
cout << "print the twiss parameters to file: "
<< UserCommandFlag[i]._PrintTwiss_twiss_file << "\n";
printlatt(UserCommandFlag[i]._PrintTwiss_twiss_file);
}
//print the close orbit else if(strcmp(CommandStr,"PrintCODFlag") == 0) {
cout << "\n";
cout << "print the close orbit to file: "
<< UserCommandFlag[i]._PrintCOD_cod_file << "\n";
getcod(dP, lastpos);
prt_cod(UserCommandFlag[i]._PrintCOD_cod_file, globval.bpm, true);
}
//print coordinates at lattice elements obtained by tracking
else if(strcmp(CommandStr,"PrintTrackFlag") == 0) {
cout << "\n";
cout << "print the tracked coordinates to file: "
<< UserCommandFlag[i]._PrintTrack_track_file << "\n";
PrintTrack(UserCommandFlag[i]._PrintTrack_track_file,
UserCommandFlag[i]._PrintTrack_x, UserCommandFlag[i]._PrintTrack_px,
UserCommandFlag[i]._PrintTrack_y, UserCommandFlag[i]._PrintTrack_py,
UserCommandFlag[i]._PrintTrack_delta,UserCommandFlag[i]._PrintTrack_ctau,
UserCommandFlag[i]._PrintTrack_nmax);
}
//print the girder
// else if(strcmp(CommandStr,"PrintGirderFlag") == 0) {
// cout << "\n";
// cout << "print the information of girder to file: "<< girder_file << "\n";
// getcod(dP, lastpos);
// prt_cod(cod_file, globval.bpm, true);
// }
// compute tunes by tracking (should be the same as by tps)
else if (strcmp(CommandStr,"TuneTracFlag") == 0) {
GetTuneTrac(1026L, 0.0, &nux, &nuy);
fprintf(stdout, "From tracking: nux = % f nuz = % f \n", nux, nuy);
}
// compute chromaticities by tracking (should be the same as by DA)
else if(strcmp(CommandStr,"ChromTracFlag") == 0) {
start = stampstart();
GetChromTrac(2L, 1026L, 1e-5, &ksix, &ksiy);
stop = stampstop(start);
fprintf(stdout, "From tracking: ksix= % f ksiz= % f \n", ksix, ksiy);
}
//generic function, to fit tunes using 1 family of quadrupoles
// if (FitTuneFlag == true) {
else if(strcmp(CommandStr,"FitTuneFlag") == 0) {
double qf_bn = 0.0, qd_bn = 0.0;
double qf_bn_fit = 0.0, qd_bn_fit = 0.0;
/* quadrupole field strength before fitting*/
qf_bn = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].qf), 1)].Elem.M->PB[HOMmax + 2];
qd_bn = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].qd), 1)].Elem.M->PB[HOMmax + 2];
/* fitting tunes*/
fprintf(stdout,
"\n Fitting tunes: %s %s, targetnux = %f, targetnuz = %f \n", UserCommandFlag[i].qf,
UserCommandFlag[i].qd,UserCommandFlag[i].targetnux, UserCommandFlag[i].targetnuz);
FitTune(ElemIndex(UserCommandFlag[i].qf), ElemIndex(UserCommandFlag[i].qd),
UserCommandFlag[i].targetnux, UserCommandFlag[i].targetnuz);
/* integrated field strength after fitting*/
qf_bn_fit = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].qf), 1)].Elem.M->PB[HOMmax + 2];
qd_bn_fit = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].qd), 1)].Elem.M->PB[HOMmax + 2];
/* print out the quadrupole strengths before and after the fitting*/
printf("Before the fitting, the quadrupole field strengths are: \n");
printf(" %s = %f, %s = %f\n", UserCommandFlag[i].qf, qf_bn, UserCommandFlag[i].qd, qd_bn); printf("After the fitting, the quadrupole field strengths are: \n");
printf(" %s = %f, %s = %f\n", UserCommandFlag[i].qf, qf_bn_fit, UserCommandFlag[i].qd, qd_bn_fit);
/* Compute and get Twiss parameters */
Ring_GetTwiss(chroma = true, 0.0);
printglob(); /* print parameter list */
}
/* specific for soleil ring in which the quadrupole is cut into 2 parts*/
//if (FitTune4Flag == true) {
else if(strcmp(CommandStr,"FitTune4Flag") == 0) {
double qf1_bn = 0.0, qf2_bn = 0.0, qd1_bn = 0.0, qd2_bn = 0.0;
double qf1_bn_fit = 0.0, qf2_bn_fit = 0.0, qd1_bn_fit = 0.0, qd2_bn_fit =
0.0;
/* quadrupole field strength before fitting*/
qf1_bn = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].qf1), 1)].Elem.M->PB[HOMmax + 2];
qf2_bn = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].qf2), 1)].Elem.M->PB[HOMmax + 2];
qd1_bn = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].qd1), 1)].Elem.M->PB[HOMmax + 2];
qd2_bn = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].qd2), 1)].Elem.M->PB[HOMmax + 2];
/* fitting tunes*/
fprintf(
stdout,
"\n Fitting tunes for Soleil ring: %s %s %s %s, targetnux = %f, targetnuz = %f \n",
UserCommandFlag[i].qf1, UserCommandFlag[i].qf2, UserCommandFlag[i].qd1, UserCommandFlag[i].qd2,
UserCommandFlag[i].targetnux, UserCommandFlag[i].targetnuz);
FitTune4(ElemIndex(UserCommandFlag[i].qf1), ElemIndex(UserCommandFlag[i].qf2),
ElemIndex(UserCommandFlag[i].qd1), ElemIndex(UserCommandFlag[i].qd2),
UserCommandFlag[i].targetnux, UserCommandFlag[i].targetnuz);
/* integrated field strength after fitting*/
qf1_bn_fit = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].qf1), 1)].Elem.M->PB[HOMmax + 2];
qf2_bn_fit = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].qf2), 1)].Elem.M->PB[HOMmax + 2];
qd1_bn_fit = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].qd1), 1)].Elem.M->PB[HOMmax + 2];
qd2_bn_fit = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].qd2), 1)].Elem.M->PB[HOMmax + 2];
/* print out the quadrupole strengths before and after the fitting*/
printf("Before the fitting, the quadrupole field strengths are: \n");
printf(" %s = %f, %s = %f, %s = %f, %s = %f\n", UserCommandFlag[i].qf1, qf1_bn,
UserCommandFlag[i].qf2, qf2_bn, UserCommandFlag[i].qd1, qd1_bn, UserCommandFlag[i].qd2, qd2_bn);
printf("After the fitting, the quadrupole field strengths are: \n");
printf(" %s = %f, %s = %f, %s = %f, %s = %f\n", UserCommandFlag[i].qf1,
qf1_bn_fit, UserCommandFlag[i].qf2, qf2_bn_fit, UserCommandFlag[i].qd1, qd1_bn_fit,
UserCommandFlag[i].qd2, qd2_bn_fit);
/* Compute and get Twiss parameters */
Ring_GetTwiss(chroma = true, 0.0);
printglob(); /* print parameter list */
}
/* fit the chromaticities*/
else if(strcmp(CommandStr,"FitChromFlag") == 0) {
double sxm1_bn = 0.0, sxm2_bn = 0.0;
double sxm1_bn_fit = 0.0, sxm2_bn_fit = 0.0;
sxm1_bn = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].sxm1), 1)].Elem.M->PB[HOMmax + 3];
sxm2_bn = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].sxm2), 1)].Elem.M->PB[HOMmax + 3];
fprintf(stdout,"\n Fitting chromaticities: %s %s, targetksix = %f, targetksiz = %f\n",
UserCommandFlag[i].sxm1, UserCommandFlag[i].sxm2, UserCommandFlag[i].targetksix,
UserCommandFlag[i].targetksiz);
FitChrom(ElemIndex(UserCommandFlag[i].sxm1), ElemIndex(UserCommandFlag[i].sxm2),
UserCommandFlag[i].targetksix, UserCommandFlag[i].targetksiz);
sxm1_bn_fit = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].sxm1), 1)].Elem.M->PB[HOMmax + 3];
sxm2_bn_fit = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].sxm2), 1)].Elem.M->PB[HOMmax + 3];
/* print out the sextupole strengths before and after the fitting*/
printf("Before the fitting, the sextupole field strengths are \n");
printf(" %s = %f, %s = %f\n", UserCommandFlag[i].sxm1, sxm1_bn, UserCommandFlag[i].sxm2, sxm2_bn); printf("After the fitting, the sextupole field strengths are: \n");
printf(" %s = %f, %s = %f\n", UserCommandFlag[i].sxm1, sxm1_bn_fit, UserCommandFlag[i].sxm2, sxm2_bn_fit);
Ring_GetTwiss(chroma = true, 0.0); /* Compute and get Twiss parameters */
printglob(); /* print parameter list */
}
// coupling calculation
else if(strcmp(CommandStr,"CouplingFlag") == 0) {
Ring_GetTwiss(chroma = true, 0.0); /* Compute and get Twiss parameters */
printlatt("linlat_coupling.out"); /* dump linear lattice functions into "linlat.dat" */
// GetEmittance(ElemIndex("cav"), true); //only for test
Coupling_Edwards_Teng();
Ring_GetTwiss(chroma = true, 0.0); /* Compute and get Twiss parameters */
printglob(); /* print parameter list */
}
// add coupling by random rotating of the full quadrupole magnets
//if (ErrorCouplingFlag == true) {
else if(strcmp(CommandStr,"ErrorCouplingFlag") == 0) {
prtmfile("flat_file.dat"); //print the elements without rotation errors
SetErr(UserCommandFlag[i].err_seed, UserCommandFlag[i].err_rms);
prtmfile("flat_file_errcoupling_full.dat"); //print the elements with rotation errors
Ring_GetTwiss(chroma = true, 0.0); /* Compute and get Twiss parameters */
printlatt("linlat_errcoupling.out"); /* dump linear lattice functions into "linlat.dat" */
Coupling_Edwards_Teng();
// GetEmittance(ElemIndex("cav"), true); //only for test
Ring_GetTwiss(chroma = true, 0.0); /* Compute and get Twiss parameters */
// printlatt();
printglob(); /* print parameter list */
}
// add coupling by random rotating of the half quadrupole magnets, delicated for soleil
else if(strcmp(CommandStr,"ErrorCoupling2Flag") == 0) {
prtmfile("flat_file.dat"); //print the elements without rotation errors
SetErr2(UserCommandFlag[i].err_seed, UserCommandFlag[i].err_rms);
prtmfile("flat_file_errcoupling_half.dat"); //print the elements with rotation errors
Ring_GetTwiss(chroma = true, 0.0); /* Compute and get Twiss parameters */
printlatt("linlat_errcoupling2.out"); /* dump linear lattice functions into "linlat.dat" */
Coupling_Edwards_Teng();
// GetEmittance(ElemIndex("cav"), true);
Ring_GetTwiss(chroma = true, 0.0); /* Compute and get Twiss parameters */
printglob(); /* print parameter list */
}
/******************************************************************************************/
// COMPUTATION PART after setting the model
/******************************************************************************************/
// computes TuneShift with amplitudes
else if(strcmp(CommandStr,"AmplitudeTuneShiftFlag") == 0) {
TunesShiftWithAmplitude(UserCommandFlag[i]._AmplitudeTuneShift_nudx_file,
UserCommandFlag[i]._AmplitudeTuneShift_nudz_file,
UserCommandFlag[i]._AmplitudeTuneShift_nxpoint,
UserCommandFlag[i]._AmplitudeTuneShift_nypoint,
UserCommandFlag[i]._AmplitudeTuneShift_nturn,
UserCommandFlag[i]._AmplitudeTuneShift_xmax,
UserCommandFlag[i]._AmplitudeTuneShift_ymax,
UserCommandFlag[i]._AmplitudeTuneShift_delta);
}
// compute tuneshift with energy
else if(strcmp(CommandStr,"EnergyTuneShiftFlag") == 0) {
TunesShiftWithEnergy(UserCommandFlag[i]._EnergyTuneShift_nudp_file,
UserCommandFlag[i]._EnergyTuneShift_npoint,
UserCommandFlag[i]._EnergyTuneShift_nturn,
UserCommandFlag[i]._EnergyTuneShift_deltamax);
}
// Computes FMA
else if(strcmp(CommandStr,"FmapFlag") == 0) {
printf("\n begin Fmap calculation for on momentum particles: \n");
fmap(UserCommandFlag[i]._FmapFlag_fmap_file,
UserCommandFlag[i]._FmapFlag_nxpoint,
UserCommandFlag[i]._FmapFlag_nypoint,
UserCommandFlag[i]._FmapFlag_nturn,
UserCommandFlag[i]._FmapFlag_xmax,
UserCommandFlag[i]._FmapFlag_ymax,
UserCommandFlag[i]._FmapFlag_delta,
UserCommandFlag[i]._FmapFlag_diffusion);
}
// Compute FMA dp
else if(strcmp(CommandStr,"FmapdpFlag") == 0) {
printf("\n begin Fmap calculation for off momentum particles: \n");
fmapdp(UserCommandFlag[i]._FmapdpFlag_fmapdp_file,
UserCommandFlag[i]._FmapdpFlag_nxpoint,
UserCommandFlag[i]._FmapdpFlag_nepoint,
UserCommandFlag[i]._FmapdpFlag_nturn,
UserCommandFlag[i]._FmapdpFlag_xmax,
UserCommandFlag[i]._FmapdpFlag_emax,
UserCommandFlag[i]._FmapdpFlag_z,
UserCommandFlag[i]._FmapdpFlag_diffusion);
}
// // if (CodeComparaisonFlag) {
// else if(strcmp(CommandStr,"CodeComparaisonFlag") == 0) {
// fmap(200, 100, 1026, -32e-3, 7e-3, 0.0, true);
// }
// MOMENTUM ACCEPTANCE
else if(strcmp(CommandStr,"MomentumAccFlag") == 0) {
bool cavityflag, radiationflag;
/* record the initial values*/
cavityflag = globval.Cavity_on;
radiationflag = globval.radiation;
if (strcmp(UserCommandFlag[i]._MomentumAccFlag_TrackDim,"6D") == 0) {
globval.Cavity_on = true;
globval.radiation = true;
}else if (strcmp(UserCommandFlag[i]._MomentumAccFlag_TrackDim,"4D") == 0) {
globval.Cavity_on = false;
globval.radiation = false;
} else {
printf("MomentumAccFlag: Error!!! DimTrack must be '4D' or '6D'\n");
exit_(1);
};
/* calculate momentum acceptance*/
printf("\n Calculate momentum acceptance: \n");
MomentumAcceptance(UserCommandFlag[i]._MomentumAccFlag_momacc_file,
UserCommandFlag[i]._MomentumAccFlag_istart,
UserCommandFlag[i]._MomentumAccFlag_istop,
UserCommandFlag[i]._MomentumAccFlag_deltaminp,
UserCommandFlag[i]._MomentumAccFlag_deltamaxp,
UserCommandFlag[i]._MomentumAccFlag_nstepp,
UserCommandFlag[i]._MomentumAccFlag_deltaminn,
UserCommandFlag[i]._MomentumAccFlag_deltamaxn,
UserCommandFlag[i]._MomentumAccFlag_nstepn,
UserCommandFlag[i]._MomentumAccFlag_nturn,
UserCommandFlag[i]._MomentumAccFlag_zmax);
/* restore the initial values*/
globval.Cavity_on = cavityflag;
globval.radiation = radiationflag;
}
// induced amplitude
else if(strcmp(CommandStr,"InducedAmplitudeFlag") == 0) {
printf("\n Calculate induced amplitude: \n");
InducedAmplitude(193L);
}
else if(strcmp(CommandStr,"EtaFlag") == 0) {
// 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("linlat_eta.out"); /* 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);
}
}
// track to get phase space
else if(strcmp(CommandStr,"PhaseSpaceFlag") == 0) {
bool cavityflag, radiationflag;
/* record the initial values*/
cavityflag = globval.Cavity_on;
radiationflag = globval.radiation;
/* set the dimension for the momentum tracking*/
if (strcmp("6D", UserCommandFlag[i]._Phase_Dim) == 0) {
globval.Cavity_on = true;
} else if (strcmp("4D", UserCommandFlag[i]._Phase_Dim) == 0) {
globval.Cavity_on = false;
} else {
printf("MomentumAccFlag: Error!!! _Phase_Dim must be '4D' or '6D'\n");
exit_(1);
};
/* setting damping */
if ( UserCommandFlag[i]._Phase_Damping == true) {
globval.radiation = true;
} else {
globval.radiation = false;
}
start = stampstart();
Phase(UserCommandFlag[i]._Phase_phase_file,
UserCommandFlag[i]._Phase_X,
UserCommandFlag[i]._Phase_Px,
UserCommandFlag[i]._Phase_Y,
UserCommandFlag[i]._Phase_Py,
UserCommandFlag[i]._Phase_delta,
UserCommandFlag[i]._Phase_ctau,
UserCommandFlag[i]._Phase_nturn);
printf("the simulation time for phase space in tracy 3 is \n");
stop = stampstop(start);
/* restore the initial values*/
globval.Cavity_on = cavityflag;
globval.radiation = radiationflag;
}
else if (strcmp(CommandStr,"TouschekFlag") == 0 ||strcmp(CommandStr,"IBSFlag") == 0 ||
strcmp(CommandStr,"TousTrackFlag") == 0 ){
// ????????????? NSRL version, Check with Soleil version "MomentumAcceptance" // IBS & TOUSCHEK
int k, n_turns;
double sigma_s, sigma_delta, tau, alpha_z, beta_z, gamma_z, eps[3];
FILE *outf;
const double Qb = 5e-9; // e charge in one bunch
// else if(strcmp(CommandStr,"TouschekFlag") == 0) {
double sum_delta[globval.Cell_nLoc + 1][2];
double sum2_delta[globval.Cell_nLoc + 1][2];
GetEmittance(globval.cav, true);
// initialize momentum aperture arrays
for (k = 0; k <= globval.Cell_nLoc; k++) {
sum_delta[k][0] = 0.0;
sum_delta[k][1] = 0.0;
sum2_delta[k][0] = 0.0;
sum2_delta[k][1] = 0.0;
}
globval.eps[Y_] = 0.008e-9;
n_turns = 40;
// get the twiss parameters
alpha_z = -globval.Ascr[ct_][ct_] * globval.Ascr[delta_][ct_]
- globval.Ascr[ct_][delta_] * globval.Ascr[delta_][delta_];
beta_z = sqr(globval.Ascr[ct_][ct_]) + sqr(globval.Ascr[ct_][delta_]);
gamma_z = (1 + sqr(alpha_z)) / beta_z;
sigma_delta = sqrt(gamma_z * globval.eps[Z_]);
sigma_s = sqrt(beta_z * globval.eps[Z_]);//50e-3
beta_z = sqr(sigma_s) / globval.eps[Z_];
alpha_z = sqrt(beta_z * gamma_z - 1);
// INCLUDE LC (LC changes sigma_s and eps_z, but has no influence on sigma_delta)
if (false) {
double newLength, bunchLengthening;
newLength = 50e-3;
bunchLengthening = newLength / sigma_s;
sigma_s = newLength;
globval.eps[Z_] = globval.eps[Z_] * bunchLengthening;
beta_z = beta_z * bunchLengthening;
gamma_z = gamma_z / bunchLengthening;
alpha_z = sqrt(beta_z * gamma_z - 1); // this doesn't change
}
Touschek(Qb, globval.delta_RF, globval.eps[X_], globval.eps[Y_],
sigma_delta, sigma_s);
// Intra Beam Scattering(IBS)
if (strcmp(CommandStr,"IBSFlag") == 0) {
// initialize eps_IBS with eps_SR
for (k = 0; k < 3; k++)
eps[k] = globval.eps[k];
for (k = 0; k < 20; k++) //prototype (looping because IBS routine doesn't check convergence)
IBS(Qb, globval.eps, eps, alpha_z, beta_z);
}
// TOUSCHEK TRACKING
// Calculate Touschek lifetime
// with the momentum acceptance which is determined by
// the RF acceptance delta_RF and the momentum aperture
// at each element location which is tracked over n turns,
//the vacuum chamber is read from the file "chamber_file"
// finally, write the momentum acceptance to file "mom_aper.out".
if (strcmp(CommandStr,"TousTrackFlag") == 0) {
globval.Aperture_on = true;
ReadCh(UserCommandFlag[i].chamber_file);
// LoadApers("/home/simon/projects/src/lattice/Apertures.dat", 1, 1); tau = Touschek(Qb, globval.delta_RF, false, globval.eps[X_],
globval.eps[Y_], sigma_delta, sigma_s, n_turns, true, sum_delta,
sum2_delta); //the TRUE flag requires apertures loaded
printf("Touschek lifetime = %10.3e hrs\n", tau / 3600.0);
outf = file_write("mom_aper.out");
for (k = 0; k <= globval.Cell_nLoc; k++)
fprintf(outf, "%4d %7.2f %5.3f %6.3f\n", k, Cell[k].S, 1e2
* sum_delta[k][0], 1e2 * sum_delta[k][1]);
fclose(outf);
}
}
else
printf("Wrong!!!!!");
}//end of looking for user defined flag
return 0;
}//end of main()