#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[])
{
const long seed = 1121;
const bool All = TRUE;
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;
// overview, on energy: 25-15
//const double x_max_FMA = 20e-3, y_max_FMA = 10e-3; //const x_max_FMA = 25e-3, y_max_FMA = 15e-3;
//const int n_x = 80, n_y = 80, n_tr = 2048;
// overview, off energy: 10-10
const double x_max_FMA = 10e-3, delta_FMA = 10e-2;
const int n_x = 80, n_dp = 80, n_tr = 2048;
//
// zoom, on energy: 8-2.5
//const double x_max_FMA = 8e-3, y_max_FMA = 2.5e-3;
//const int n_x = 64, n_y = 15, n_tr = 2048;
// zoom, off energy: 7-3
//const double x_max_FMA = 3e-3, delta_FMA = 7e-2;
//const int n_x = 28, n_dp = 56, n_tr = 2048;
double nux=0.0, nuz=0.0, ksix=0.0, ksiz=0.0;
bool chroma;
double dP = 0.0;
long lastpos = -1L;
char str1[S_SIZE];
/************************************************************************
start read in files and flags
*************************************************************************/
read_script(argv[1], true);
/************************************************************************
end read in files and flags
*************************************************************************/
// if (true)
// // Read_Lattice("/home/nadolski/codes/tracy/maille/soleil/solamor2_tracy3");
// Read_Lattice(argv[1]); //sets some globval params
// else
// rdmfile("flat_file.dat"); //instead of reading lattice file, get data from flat file
//no_sxt(); //turns off sextupoles
// Ring_GetTwiss(true, 0e-2); //gettwiss computes one-turn matrix arg=(w or w/o chromat, dp/p)
//get_matching_params_scl();
//get_alphac2();
//GetEmittance(ElemIndex("cav"), true);
//prt_lat("linlat.out", globval.bpm, true); /* print lattice file for nsrl-ii*/
prtmfile("flat_file.dat"); // writes flat file /* very important file for debug*/
//prt_chrom_lat(); //writes chromatic functions into chromlat.out
// printlatt(); /* print out lattice functions */
/* print lattice file */
// prt_lat("linlatBNL.out", globval.bpm, All); // BNL print for all elements
printlatt(); /* SOLEIL print out lattice functions */
printglob();
// Flag factory
// bool TuneTracFlag = true;
// bool ChromTracFlag = true;
//*************************************************************
//=============================================================
// 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("Apertures.dat"); /* read vacuum chamber from a file "Apertures.dat" , soleil version*/
//LoadApers("Apertures.dat", 1.0, 1.0); /* read vacuum chamber definition for bnl */
PrintCh();
// 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){
GetChromTrac(2L, 1026L, 1e-5, &ksix, &ksiz);
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 */
}
//SetKLpar(ElemIndex("QT"), 1, 2L, 0.001026770838382);
// coupling calculation
if (CouplingFlag == true){
Ring_GetTwiss(chroma=true, 0.0); /* Compute and get Twiss parameters */
printlatt(); /* dump linear lattice functions into "linlat.dat" */
// 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();
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(); /* 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(); /* for thick sextupoles */
Ring_GetTwiss(chroma=true, 0.0); /* Compute and get Twiss parameters */
printglob();
}
}
/* print parameter list */
// PX2 chicane
// if (PX2Flag ==true){
// setPX2chicane();
// //get closed orbit
// getcod (0.0, &lastpos);
// printcod();
// Ring_GetTwiss(chroma=true, 0.0); /* Compute and get Twiss parameters */
// printglob(); /* print parameter list */
// }
// Computes FMA
if (FmapFlag == true){
if (ChamberFlag == true ){
if (ExperimentFMAFlag == true)
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{
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){
// SOLEIL
fmap(100,50,1026,32e-3,7e-3,0.0,true);
//fmap(200,100,1026,-32e-3,7e-3,0.0,true);
}
if (MomentumAccFlag == true){
//MomentumAcceptance(10L, 28L, 0.01, 0.05, 4L, -0.01, -0.05, 4L);
MomentumAcceptance(1L, 28L, 0.01, 0.05, 40L, -0.01, -0.05, 40L);
// MomentumAcceptance(1L, 108L, 0.01, 0.05, 100L, -0.01, -0.05, 100L);
}
// computes Tuneshift with amplitudes
if (TuneShiftFlag == true){
if (ChamberFlag == true ){
TunesShiftWithAmplitude(31L,21L,516L,0.025,0.005,dP);
//NuDp(31L,516L,0.06);
//NuDp(31L,1026L,0.06);
}
else{
TunesShiftWithAmplitude(50L,30L,516L,0.035,0.02,dP);
TunesShiftWithEnergy(31L,1026L,0.06);
}
}
// if (SigmaFlag == true){printsigma();
// }
//
// 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);
}
}
//*********************************************************************************
//---------------------------------------------------------------------------------------------------------------------------------
// 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(20e-3, 10e-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);
}
//
// 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;
if (false) {
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;
}
//sigma_delta = 7.70e-04; //410:7.70e-4, 411:9.57e-4, 412:9.12e-4
//globval.eps[X_] = 0.326e-9; //410:3.26e-10, 411:2.63e-10, 412:2.01e-10
globval.eps[Y_] = 0.008e-9;
//sigma_s = 9.73e-3; //410:9.73e-3, 411:12.39e-3, 412:12.50e-3/10.33e-3
//globval.eps[Z_] = sigma_delta*sigma_s;
//globval.delta_RF given by cav voltage in lattice file
//globval.delta_RF = 6.20e-2; //410:6.196e-2, 411:5.285e-2, 412:4.046e-2/5.786e-2
n_turns = 490; //410:490(735), 411:503(755), 412:439(659)/529(794)
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
}
//globval.eps[X_] = 0.362e-9;
//sigma_delta = 1.04e-3;
//sigma_s = 14.8e-3;
Touschek(Qb, globval.delta_RF, globval.eps[X_], globval.eps[Y_],
sigma_delta, sigma_s);
// IBS
if (false) {
// 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
if (false) {
globval.Aperture_on = true;
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);
}
}
}