Newer
Older
#define ORDER 1
int no_tps = ORDER; // arbitrary TPSA order is defined locally
//***************************************************************************************
//
// 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;
// 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;
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
// 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 */
// prt_lat("linlatBNL.out", globval.bpm, All); // BNL print for all elements
printlatt(); /* SOLEIL print out lattice functions */
//*************************************************************
//=============================================================
// 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){
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" */
// 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();
}
// 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 (CodeComparaisonFlag){
// SOLEIL
fmap(100,50,1026,32e-3,7e-3,0.0,true);
//fmap(200,100,1026,-32e-3,7e-3,0.0,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 ){
//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
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
//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
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);
}
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
//
// 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(ElemIndex("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);
}
}
}