diff --git a/tracy/tools/testtracy.cc b/tracy/tools/testtracy.cc deleted file mode 100644 index 84546f2079b159f3bc502580bd89d4b884b6a28a..0000000000000000000000000000000000000000 --- a/tracy/tools/testtracy.cc +++ /dev/null @@ -1,452 +0,0 @@ -/* 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); - } - -}