Code owners
Assign users and groups as approvers for specific file changes. Learn more.
soltracy.cc 38.48 KiB
/************************************
Current revision $Revision$
On branch $Name$
Latest change $Date$ by $Author$
*************************************/
#define ORDER 1
//#define ORDER 4
int no_tps = ORDER; // arbitrary TPSA order is defined locally
extern bool freq_map;
#if HAVE_CONFIG_H
#include <config.h>
#endif
#if MPI_EXEC
//library of parallel computation,must be before "stdio.h"
#include <mpi.h>
#endif
#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
//long j=0L; //initialize the for loop to read the files from parallel computing
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];
#if MPI_EXEC
//Initialize parallel computation
MPI_Init(&argc, &argv);
#endif
/************************************************************************
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.out");
// 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";
}
//deactive quadrupole fringe field
else if(strcmp(CommandStr,"QuadFringeOffFlag") == 0) { globval.quad_fringe = false;
cout << "\n";
cout << "globval.quad_fringe is " << globval.quad_fringe << "\n";
}
//set RF voltage
else if(strcmp(CommandStr,"RFvoltageFlag") == 0) {
fprintf(stdout, "\nSetting RF voltage:\n");
fprintf(stdout, " Old RF voltage is: %lf [MV]\n", get_RFVoltage(globval.cav)
/ 1e6);
set_RFVoltage(globval.cav, UserCommandFlag[i].RFvolt);
fprintf(stdout, " 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
// Based on the function error_and_correction() in nsls-ii_lib_templ.h
else if (strcmp(CommandStr,"ReadaefileFlag") == 0) {
// load misalignments
cout << "Read misalignment errors from file: " << UserCommandFlag[i].ae_file << endl;
LoadAlignTol(UserCommandFlag[i].ae_file, true, 1.0, true, 0);
Ring_GetTwiss(true, 0.0);
}
//do COD correction using SVD method.
else if(strcmp(CommandStr,"CODCorrectFlag") == 0) {
cout << "Begin Closed Orbit correction: " << endl;
if(n_orbit < 1){
cout << "interation number n_obit should NOT small than 1!!!" << endl;
exit_(1);
}
CODCorrect(hcorr_file,vcorr_file, n_orbit,nwh,nwv);
Ring_GetTwiss(true, 0.0);
prt_cod("summary_miserr_codcorr.out", globval.bpm, true);
}
// 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 multipole field errors from file: %s \n", UserCommandFlag[i].fe_file);
LoadFieldErrs(UserCommandFlag[i].fe_file, true, 1.0, true, 1);
Ring_GetTwiss(true, 0.0);
prtmfile("flat_file_fefile.out");
}
// read multipole errors from a file; specific for soleil lattice
else if(strcmp(CommandStr,"ReadMultipoleFlag") == 0) {
fprintf(stdout,"\n Read Multipoles field error for SOLEIL lattice with thick sextupoles, from file:\n");
fprintf(stdout," %s\n",multipole_file);
ReadFieldErr(multipole_file);
//first print the full lattice with error as a flat file
prtmfile("flat_file_errmultipole.out"); // writes flat file with all element errors /* very important file for debug*/
Ring_GetTwiss(chroma = true, 0.0); /* Compute and get Twiss parameters */
printglob();
}
// read the sources of coupling from a file; for soleil lattice
else if(strcmp(CommandStr,"ReadVirtualSkewquadFlag") == 0) {
fprintf(stdout,"\n Read virtual skew quadrupole setting of SOLEIL lattice from file:\n");
fprintf(stdout," %s \n",virtualskewquad_file);
ReadVirtualSkewQuad(virtualskewquad_file);
//first print the full lattice with sources of coupling as a flat file prtmfile("flat_file_skewquad.out"); /* very important file for debug*/
//get the coupling
Coupling_Edwards_Teng();
}
//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 closed orbit
else if(strcmp(CommandStr,"PrintCODFlag") == 0) {
cout << "\n";
cout << "print the closed 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*/
fprintf(stdout, "Before the fitting, the quadrupole field strengths are: \n");
fprintf(stdout, " %s = %f, %s = %f\n", UserCommandFlag[i].qf, qf_bn, UserCommandFlag[i].qd, qd_bn);
fprintf(stdout, "After the fitting, the quadrupole field strengths are: \n");
fprintf(stdout, " %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");
fprintf(stdout, " %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");
fprintf(stdout, " %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*/
fprintf(stdout, "Before the fitting, the sextupole field strengths are \n");
fprintf(stdout, " %s = %f, %s = %f\n", UserCommandFlag[i].sxm1, sxm1_bn, UserCommandFlag[i].sxm2, sxm2_bn);
fprintf(stdout, "After the fitting, the sextupole field strengths are: \n");
fprintf(stdout, " %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.out"); //print the elements without rotation errors
SetErr(UserCommandFlag[i].err_seed, UserCommandFlag[i].err_rms);
prtmfile("flat_file_errcoupling_full.out"); //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.out"); //print the elements without rotation errors
SetErr2(UserCommandFlag[i].err_seed, UserCommandFlag[i].err_rms);
prtmfile("flat_file_errcoupling_half.out"); //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) {
fprintf(stdout, "\n begin Fmap calculation for on momentum particles: \n");
#if MPI_EXEC
//initialization for parallel computing
int myid = 0, numprocs = 0;
int namelen = 0;
char processor_name[MPI_MAX_PROCESSOR_NAME]="";
MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
MPI_Comm_rank(MPI_COMM_WORLD, &myid);
MPI_Get_processor_name(processor_name, &namelen);
fprintf(stdout, "\n Beginning parallel computing: %d of %d is on %s\n", myid, numprocs, processor_name);
fprintf(stdout, "\n Begin Fmap calculation for on momentum particles: \n");
//Each core or process, characterized by myid, will track different region of fmap
fmap_p(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,
UserCommandFlag[i]._FmapFlag_printloss,
numprocs, myid);
//Synchronize all cores, till all cores finish fmap of different region,then all cores will proceed.
MPI_Barrier(MPI_COMM_WORLD);
//Collecting data from all files generated by cores, then write to fmap_p.out
if(myid == 0){
FILE *outf;
FILE *outf_loss;
char FmapLossFile[max_str+5]=" ";
// open frequency map file for final results
if ((outf = fopen(UserCommandFlag[i]._FmapFlag_fmap_file, "w")) == NULL){
fprintf(stdout, "psoltracy: error while opening file %s\n", UserCommandFlag[i]._FmapFlag_fmap_file);
exit_(1);
}
// open frequency map loss file for final results
if (UserCommandFlag[i]._FmapFlag_printloss){
strcpy(FmapLossFile, UserCommandFlag[i]._FmapFlag_fmap_file);
strcat(FmapLossFile, ".loss");
if ((outf_loss = fopen(FmapLossFile, "w")) == NULL){
fprintf(stdout, "psoltracy: error while opening file %s\n", UserCommandFlag[i]._FmapFlag_fmap_file);
exit_(1);
}
}
FILE *fp, *fp_loss;
char buffer[81];
char FmapFile[max_str];
for(int j = 0; j < numprocs; j++){
FmapFile[0]='\0';
sprintf(FmapFile,"%d",j);
strcat(FmapFile, UserCommandFlag[i]._FmapFlag_fmap_file);
if((fp = fopen(FmapFile, "r")) == NULL) {
fprintf(stdout, "%s: error while opening file.\n", FmapFile);
exit_(1);
}
// concatanation of the results, add current j-file to the master file
while(fgets(buffer, max_str, fp) != NULL){
fputs(buffer, outf);
buffer[0]='\0';
}
fclose(fp);
if (UserCommandFlag[i]._FmapFlag_printloss){
strcpy(FmapLossFile,FmapFile);
strcat(FmapLossFile,".loss");
if((fp_loss = fopen(FmapLossFile, "r")) == NULL) {
fprintf(stdout, "%s: error while opening file.\n", FmapLossFile);
exit_(1);
}
// concatanation of the results, add current j-file to the master file
while(fgets(buffer, max_str, fp_loss) != NULL){
fputs(buffer, outf_loss);
buffer[0]='\0';
}
fclose(fp_loss);
}
}
fclose(outf);
if (UserCommandFlag[i]._FmapFlag_printloss){
fclose(outf_loss);
}
}
MPI_Barrier(MPI_COMM_WORLD);
printf("Process %d finished for on momentum fmap.\n", myid);
#else // Single processor
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,
UserCommandFlag[i]._FmapFlag_printloss);
#endif // MPI
}
// Compute FMA dp
else if(strcmp(CommandStr,"FmapdpFlag") == 0) {
fprintf(stdout, "\n begin Fmap calculation for off momentum particles: \n");
#if MPI_EXEC
//initialization for parallel computing
int myid=0, numprocs=0;
int namelen=0;
char processor_name[MPI_MAX_PROCESSOR_NAME]="";
MPI_Comm_size(MPI_COMM_WORLD, &numprocs); MPI_Comm_rank(MPI_COMM_WORLD, &myid);
MPI_Get_processor_name(processor_name, &namelen);
printf("Beginning parallel computing: %d of %d is on %s\n", myid, numprocs, processor_name);
fprintf(stdout, "\n begin Fmap calculation for off momentum particles: \n");
fmapdp_p(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,
UserCommandFlag[i]._FmapdpFlag_printloss, numprocs, myid);
//Synchronize all cores, till all cores finish fmap of different region,then all cores will proceed.
MPI_Barrier(MPI_COMM_WORLD);
//Collecting data from all files generated by cores, then write to fmap_p.out
if(myid == 0){
FILE *outf, *outf_loss;
char FmapdpLossFile[max_str+5]=" ";
if ((outf = fopen(UserCommandFlag[i]._FmapdpFlag_fmapdp_file, "w")) == NULL){
fprintf(stdout, "psoltracy: error while opening file %s\n",
UserCommandFlag[i]._FmapdpFlag_fmapdp_file);
exit_(1);
}
// open frequency map loss file for final results
if (UserCommandFlag[i]._FmapdpFlag_printloss){
strcpy(FmapdpLossFile, UserCommandFlag[i]._FmapdpFlag_fmapdp_file);
strcat(FmapdpLossFile, ".loss");
if ((outf_loss = fopen(FmapdpLossFile, "w")) == NULL){
fprintf(stdout, "psoltracy: error while opening file %s\n",
UserCommandFlag[i]._FmapdpFlag_fmapdp_file);
exit_(1);
}
}
FILE *fp, *fp_loss;
char buffer[81];
char FmapdpFile[50];
for(int j = 0; j < numprocs; j++){
FmapdpFile[0]='\0';
sprintf(FmapdpFile,"%d",j);
strcat(FmapdpFile, UserCommandFlag[i]._FmapdpFlag_fmapdp_file);
if((fp = fopen(FmapdpFile, "r")) == NULL){
fprintf(stdout, "%s: error while opening file.\n", FmapdpFile);
exit_(1);
}
// concatanation of the results, add current j-file to the master file
while(fgets(buffer, max_str, fp) != NULL){
fputs(buffer, outf);
buffer[0]='\0';
}
fclose(fp);
if (UserCommandFlag[i]._FmapdpFlag_printloss){
strcpy(FmapdpLossFile, FmapdpFile);
strcat(FmapdpLossFile,".loss");
if((fp_loss = fopen(FmapdpLossFile, "r")) == NULL) {
fprintf(stdout, "%s: error while opening file.\n", FmapdpLossFile);
exit_(1);
}
// concatanation of the results, add current j-file to the master file
while(fgets(buffer, max_str, fp_loss) != NULL){ fputs(buffer, outf_loss);
buffer[0]='\0';
}
fclose(fp_loss);
} // if printloss
} // for loop
fclose(outf);
if (UserCommandFlag[i]._FmapFlag_printloss){
fclose(outf_loss);
}
} // if myid
MPI_Barrier(MPI_COMM_WORLD);
fprintf(stdout, "Process %d finished off momentum fmap.\n", myid);
#else //single processor
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,
UserCommandFlag[i]._FmapdpFlag_printloss);
#endif
}
// // 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*/
fprintf(stdout, "\n Calculate momentum acceptance: \n");
#if MPI_EXEC
/* calculate momentum acceptance*/
//initialization for parallel computing
int myid=0, numprocs=0;
int namelen=0;
char processor_name[MPI_MAX_PROCESSOR_NAME]="";
MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
MPI_Comm_rank(MPI_COMM_WORLD, &myid);
MPI_Get_processor_name(processor_name, &namelen);
printf("Beginning parallel computing: %d of %d is on %s\n", myid, numprocs, processor_name);
fprintf(stdout, "\n Calculate momentum acceptance: \n");
MomentumAcceptance_p(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,
numprocs,myid);
//merge the files
//Synchronize all cores, till finish cal of different region,
//then files from the cores will be processed.
MPI_Barrier(MPI_COMM_WORLD);
//Collect data from all files generated by different cores, then write to user defined file
if(myid==0){
//merge the phase.out from different cpus
FILE *outf1; //phase.out, for debugging
FILE *fp1;
char buffer1[81];
char PhaseFile[50];
if ((outf1 = fopen("phase_p.out", "w")) == NULL){
fprintf(stdout, "psoltracy: error while opening file %s\n", "phase_p.out");
exit_(1);
}
for(int j=0; j<numprocs; j++){
PhaseFile[0]='\0';
sprintf(PhaseFile,"%d",j);
strcat(PhaseFile,"phase.out");
if((fp1 = fopen(PhaseFile, "r")) == NULL){
fprintf(stdout, "%s: error while opening file.\n", PhaseFile);
exit_(1);
}
while(fgets(buffer1, 80, fp1) != NULL){
fputs(buffer1, outf1);
buffer1[0]='\0';
}
fclose(fp1);
} // for loop
//collect data for the momentum acceptance
FILE *outf2; //momentum acceptance
FILE *fp2;
char buffer2[81];
char MAFile[50];
if ((outf2 = fopen(UserCommandFlag[i]._MomentumAccFlag_momacc_file, "w")) == NULL){
fprintf(stdout, "psoltracy: error while opening file %s\n", UserCommandFlag[i]._MomentumAccFlag_momacc_file);
exit_(1);
}
//Collect data in Positive region
for(int j=0; j<numprocs; j++){
MAFile[0]='\0';
sprintf(MAFile,"%d",j);
strcat(MAFile,UserCommandFlag[i]._MomentumAccFlag_momacc_file);
if((fp2 = fopen(MAFile, "r")) == NULL){
fprintf(stdout, "%s: error while opening file.\n", MAFile);
exit_(1);
}
while(fgets(buffer2, 80, fp2) != NULL){
if(strstr(buffer2,"Negative")!=0) break;
fputs(buffer2, outf2);
buffer2[0]='\0'; //reset buffer to NULL
} buffer2[0]='\0';
fclose(fp2);
}
fprintf(outf2,"\n"); // A void line
//Collect data in Nagative region
for(int j=0; j<numprocs; j++){
MAFile[0]='\0';
sprintf(MAFile,"%d",j);
strcat(MAFile,UserCommandFlag[i]._MomentumAccFlag_momacc_file);
if((fp2 = fopen(MAFile, "r")) == NULL){
fprintf(stdout, "%s: error while opening file.\n", MAFile);
exit_(1);
}
bool found=false;
while(fgets(buffer2, 80, fp2) != NULL){
if((strstr(buffer2,"Negative")==0) && (!found)) { buffer2[0]='\0'; continue;}
if((strstr(buffer2,"Negative")!=0) ) { found=true; buffer2[0]='\0'; continue;}
fputs(buffer2, outf2);
buffer2[0]='\0';
} // while
buffer2[0]='\0';
fclose(fp2);
} // for
fclose(outf2);
}
MPI_Barrier(MPI_COMM_WORLD);
printf("process %d finished momentum acceptance.\n", myid);
#else // single processor coputation
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);
#endif
/* restore the initial values*/
globval.Cavity_on = cavityflag;
globval.radiation = radiationflag;
}
// induced amplitude
else if(strcmp(CommandStr,"InducedAmplitudeFlag") == 0) {
fprintf(stdout, "\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("6D phase space tracking: \n 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);
}
}
/*
To do ID correction.
Based on parts of functions get_param( ) & ID_corr(), etc in nsls-ii_lib.cc
*/ else if(strcmp(CommandStr,"IDCorrFlag") == 0) {
int k = 0;
cout << endl;
cout << "************ Begin ID matching: **********" <<endl;
// read the family index of quadrupoles used for ID correction
if (N_calls > 0) {
if (N_Fam > N_Fam_max) {
printf("get_param: N_Fam > N_Fam_max: %d (%d)\n",
N_Fam, N_Fam_max);
exit_(0);
}
for (k = 0; k < N_Fam; k++)
Q_Fam[k] = ElemIndex(IDCq_name[k]);
}
//For debug
if(trace){
cout << "scl_dbetax = " << scl_dbetax << " scl_dbetay = " << scl_dbetay <<endl;
cout << "scl_dnux = " << scl_dnux << " scl_dnuy = " << scl_dnuy << endl;
cout << "scl_nux = " << scl_nux << " scl_nuy = " << scl_nuy << endl;
cout << "ID_step = " << ID_step <<endl;
cout << "N_calls = " << N_calls << " N_steps = " << N_steps <<endl;
cout << "Number of quadrupole families used for ID correction: " << N_Fam << endl;
cout << "Quadrupoles used for ID correction: " << endl;
for (k = 0; k < N_Fam; k++)
fprintf(stdout, "%d\n",Q_Fam[k]);
}
//initialization
ini_ID_corr();
printlatt("linlat00.out");
ID_corr0();
printlatt("linlat01.out");
// exit(0);
ini_ID_corr();
printlatt("linlat0.out");
ID_corr(N_calls,N_steps);
printlatt("linlat1.out");
}
else{
printf("Command string %s is invalid!!!!!\n ",CommandStr);
exit_(1);
}//give an error message
}//end of looking for user defined flag
#if MPI_EXEC
MPI_Finalize(); //finish parallel computation
#endif
return 0;
}//end of main()