39 #include <boost/filesystem.hpp>
46 #define CHECK_CYC_FSCANF_EOF(arg) if (arg == EOF)\
47 throw GeneralClassicException("Cyclotron::getFieldFromFile",\
48 "fscanf returned EOF at " #arg);
59 fieldType_m(right.fieldType_m),
60 fmapfn_m(right.fmapfn_m),
61 rffrequ_m(right.rffrequ_m),
62 rfphi_m(right.rfphi_m),
63 escale_m(right.escale_m),
64 superpose_m(right.superpose_m),
65 symmetry_m(right.symmetry_m),
66 rinit_m(right.rinit_m),
67 prinit_m(right.prinit_m),
68 phiinit_m(right.phiinit_m),
69 zinit_m(right.zinit_m),
70 pzinit_m(right.pzinit_m),
71 spiralFlag_m(right.spiralFlag_m),
72 trimCoilThreshold_m(right.trimCoilThreshold_m),
73 typeName_m(right.typeName_m),
75 bscale_m(right.bscale_m),
76 trimcoils_m(right.trimcoils_m),
81 fmLowE_m(right.fmLowE_m),
82 fmHighE_m(right.fmHighE_m),
83 RFfilename_m(right.RFfilename_m),
84 RFFCoeff_fn_m(right.RFFCoeff_fn_m),
85 RFVCoeff_fn_m(right.RFVCoeff_fn_m) {
98 double* br,
double* bz) {
100 trimcoil->applyField(r,tet_rad,z,br,bz);
105 const double tet_rad,
106 double& br,
double& bz) {
183 if (boost::filesystem::exists(
fmapfn_m)) {
187 "Failed to open file '" +
fmapfn_m +
188 "', please check if it exists");
213 "RFPHI not defined for CYCLOTRON!");
226 "RFFREQ not defined for CYCLOTRON!");
239 "SUPERPOSE not defined for CYCLOTRON!");
284 "EScale not defined for CYCLOTRON!");
367 "TYPE is not defined for CYCLOTRON!");
368 }
else if (
typeName_m == std::string(
"RING")) {
370 }
else if (
typeName_m == std::string(
"CARBONCYCL")) {
372 }
else if (
typeName_m == std::string(
"CYCIAE")) {
374 }
else if (
typeName_m == std::string(
"AVFEQ")) {
376 }
else if (
typeName_m == std::string(
"FFA")) {
378 }
else if (
typeName_m == std::string(
"BANDRF")) {
380 }
else if (
typeName_m == std::string(
"SYNCHROCYCLOTRON")) {
387 bool flagNeedUpdate =
false;
394 flagNeedUpdate =
true;
396 <<
" out of the global aperture of cyclotron!" <<
endl;
402 if (flagNeedUpdate) {
404 <<
" out of the field map boundary!" <<
endl;
410 if (flagNeedUpdate) {
417 return flagNeedUpdate;
423 const double rad = std::hypot(
R[0],
R[1]);
430 else if ((
R[0] < 0) && (
R[1] <= 0)) tet =
Physics::pi + tempv;
435 double tet_rad = tet;
444 double brint = 0.0, btint = 0.0, bzint = 0.0;
446 if ( this->
interpolate(rad, tet_rad, brint, btint, bzint) ) {
449 double br = - brint *
R[2];
452 double bt = - btint /
rad *
R[2];
472 std::vector<Fieldmap *>::const_iterator fi =
RFfields_m.begin();
473 std::vector<double>::const_iterator rffi =
rffrequ_m.begin();
474 std::vector<double>::const_iterator rfphii =
rfphi_m.begin();
475 std::vector<double>::const_iterator escali =
escale_m.begin();
476 std::vector<bool>::const_iterator superposei =
superpose_m.begin();
477 std::vector<std::vector<double>>::const_iterator rffci =
rffc_m.begin();
478 std::vector<std::vector<double>>::const_iterator rfvci =
rfvc_m.begin();
480 double xBegin(0), xEnd(0), yBegin(0), yEnd(0), zBegin(0), zEnd(0);
483 for (; fi !=
RFfields_m.end(); ++fi, ++rffi, ++rfphii, ++escali, ++superposei) {
484 (*fi)->getFieldDimensions(xBegin, xEnd, yBegin, yEnd, zBegin, zEnd);
485 if (fcount > 0 && *superposei ==
false)
continue;
491 if (temp_R(0) < xBegin || temp_R(0) > xEnd ||
492 temp_R(1) < yBegin || temp_R(1) > yEnd ||
493 temp_R(2) < zBegin || temp_R(2) > zEnd)
continue;
495 Vector_t tmpE(0.0, 0.0, 0.0), tmpB(0.0, 0.0, 0.0);
497 if ((*fi)->getFieldstrength(temp_R, tmpE, tmpB))
continue;
501 double frequency = (*rffi);
502 double ebscale = (*escali);
506 for (
const double fcoef : (*rffci)) {
507 powert *= (t * 1
e-9);
508 frequency += fcoef * powert;
511 for (
const double vcoef : (*rfvci)) {
512 powert *= (t * 1
e-9);
513 ebscale += vcoef * powert;
519 E += ebscale *
std::cos(phase) * tmpE;
520 B -= ebscale *
std::sin(phase) * tmpB;
528 phase_print =
std::fmod(phase_print, 360) - 360.0;
530 *
gmsg <<
endl <<
"Gap 1 phase = " << phase_print <<
" Deg" <<
endl;
531 *
gmsg <<
"Gap 1 E-Field = (" << E[0] <<
"/" << E[1] <<
"/" << E[2] <<
")" <<
endl;
532 *
gmsg <<
"Gap 1 B-Field = (" << B[0] <<
"/" << B[1] <<
"/" << B[2] <<
")" <<
endl;
533 *
gmsg <<
"RF Frequency = " << frequency <<
" MHz" <<
endl;
537 phase_print =
std::fmod(phase_print, 360) - 360.0;
539 *
gmsg <<
endl <<
"Gap 2 phase = " << phase_print <<
" Deg" <<
endl;
540 *
gmsg <<
"Gap 2 E-Field = (" << E[0] <<
"/" << E[1] <<
"/" << E[2] <<
")" <<
endl;
541 *
gmsg <<
"Gap 2 B-Field = (" << B[0] <<
"/" << B[1] <<
"/" << B[2] <<
")" <<
endl;
542 *
gmsg <<
"RF Frequency = " << frequency <<
" MHz" <<
endl;
553 const double& tet_rad,
double& br,
554 double& bt,
double& bz) {
562 *
gmsg <<
"* Finalize cyclotron" <<
endl;
571 const int krl,
const int lpr) {
573 double C[5][5][3], FAC[3];
674 for (j = 0; j < 5; j++) {
675 result += C[j][krl][kor] * *(f + j * lpr);
678 return result / (FAC[kor] *
std::pow(dx, (kor + 1)));
683 const double& tet_rad,
691 const int ir = (
int)xir;
694 const double wr1 = xir - (double)ir;
696 const double wr2 = 1.0 - wr1;
705 const double wt1 = xit - (double)it;
706 const double wt2 = 1.0 - wt1;
713 int r1t1, r2t1, r1t2, r2t2;
723 r1t1 = it + ntetS * ir - 1;
734 r2t1 =
idx(ir + 1, it);
735 r1t2 =
idx(ir, it + 1);
736 r2t2 =
idx(ir + 1, it + 1);
771 *
gmsg <<
"* Read field data from PSI format field map file" <<
endl;
776 *
gmsg <<
"* Read data from 450MeV Carbon cyclotron field file" <<
endl;
781 *
gmsg <<
"* Read data from 100MeV H- cyclotron CYCIAE-100 field file" <<
endl;
786 *
gmsg <<
"* Read AVFEQ data (Riken)" <<
endl;
791 *
gmsg <<
"* Read FFA data MSU/FNAL" <<
endl;
796 *
gmsg <<
"* Read both median plane B field map and 3D E field map of RF cavity for compact cyclotron" <<
endl;
801 *
gmsg <<
"* Read midplane B-field, 3D RF fieldmaps, and text files with RF frequency/Voltage coefficients for Synchrocyclotron" <<
endl;
839 int dkFromEdge = k - kEdge;
840 int index =
idx(i, k);
841 int indexkEdge =
idx(i, kEdge);
857 int irtak = i - iredg;
858 int index =
idx(i, k);
859 int indexredg =
idx(iredg, k);
888 int istart =
idx(i, 0);
926 for (
int i = 0; i < nrad; i++) {
941 this->
read(scaleFactor);
953 *
gmsg <<
"* ----------------------------------------------" <<
endl;
954 *
gmsg <<
"* READ IN RING FIELD MAP " <<
endl;
955 *
gmsg <<
"* (The first data block is useless) " <<
endl;
956 *
gmsg <<
"* ----------------------------------------------" <<
endl;
960 f = std::fopen(
fmapfn_m.c_str(),
"r");
963 *
gmsg <<
"* Minimal radius of measured field map: " <<
BP_m.
rmin_m <<
" [mm]" <<
endl;
980 for (
int i = 0; i < 13; i++) {
993 for (
int i = 0; i < 5; i++) {
999 for (
int i = 0; i < 4; i++) {
1003 for (
int i = 0; i < lpar; i++) {
1006 for (
int i = 0; i < 6; i++) {
1010 for (
int i = 0; i < 10000; i++) {
1012 if (std::strcmp(fout,
"LREC=") == 0)
break;
1015 for (
int i = 0; i < 5; i++) {
1026 *
gmsg <<
"* Read-in loop one block per radius" <<
endl;
1030 for (
int dummy = 0; dummy < 6; dummy++) {
1053 *
gmsg <<
"* Field Map read successfully!" <<
endl <<
endl;
1078 std::vector<double> rv;
1079 std::vector<double> thv;
1080 std::vector<double> xv;
1081 std::vector<double> yv;
1082 std::vector<double> bzv;
1085 *
gmsg <<
"* ----------------------------------------------" <<
endl;
1086 *
gmsg <<
"* READ IN FFA FIELD MAP " <<
endl;
1087 *
gmsg <<
"* ----------------------------------------------" <<
endl;
1091 std::ifstream file_to_read(
fmapfn_m.c_str());
1092 const int max_num_of_char_in_a_line = 128;
1093 const int num_of_header_lines = 1;
1096 for (
int i = 0; i < num_of_header_lines; ++i)
1097 file_to_read.ignore(max_num_of_char_in_a_line,
'\n');
1100 while(!file_to_read.eof()) {
1101 double r, th, x, y, bz;
1102 file_to_read >> r >> th >> x >> y >> bz;
1103 if ((
int)th != 360) {
1112 double maxtheta = 360.0;
1115 double rmax = rv.back();
1118 for (vit = rv.begin(); *vit <=
BP_m.
rmin_m; ++vit) {}
1126 *
gmsg <<
"* Minimal radius of measured field map: " << 1000.0 *
BP_m.
rmin_m <<
" [mm]" <<
endl;
1127 *
gmsg <<
"* Maximal radius of measured field map: " << 1000.0 * rmax <<
" [mm]" <<
endl;
1128 *
gmsg <<
"* Stepsize in radial direction: " << 1000.0 *
BP_m.
delr_m <<
" [mm]" <<
endl;
1130 *
gmsg <<
"* Maximal angle of measured field map: " << maxtheta <<
" [deg]" <<
endl;
1159 *
gmsg <<
"* Number of elements read: " << count <<
endl;
1160 *
gmsg <<
"* Field Map read successfully!" <<
endl <<
endl;
1167 *
gmsg <<
"* ----------------------------------------------" <<
endl;
1168 *
gmsg <<
"* READ IN AVFEQ CYCLOTRON FIELD MAP " <<
endl;
1169 *
gmsg <<
"* ----------------------------------------------" <<
endl;
1187 f = std::fopen(
fmapfn_m.c_str(),
"r");
1190 *
gmsg <<
"* Minimal radius of measured field map: " <<
BP_m.
rmin_m <<
" [mm]" <<
endl;
1195 *
gmsg <<
"* Maximal radius of measured field map: " << rmax <<
" [mm]" <<
endl;
1207 *
gmsg <<
"* Maximal angle of measured field map: " << tetmax <<
" [deg]" <<
endl;
1223 *
gmsg <<
"* Total stored grid point number ( ntetS * nrad ): "
1249 *
gmsg <<
"* Number of elements read: " << count <<
endl;
1250 *
gmsg <<
"* Field Map read successfully!" <<
endl <<
endl;
1257 *
gmsg <<
"* ----------------------------------------------" <<
endl;
1258 *
gmsg <<
"* READ IN CARBON CYCLOTRON FIELD MAP " <<
endl;
1259 *
gmsg <<
"* ----------------------------------------------" <<
endl;
1263 f = std::fopen(
fmapfn_m.c_str(),
"r");
1266 *
gmsg <<
"* Minimal radius of measured field map: " <<
BP_m.
rmin_m <<
" [mm]" <<
endl;
1296 *
gmsg <<
"* Adding a guard cell along azimuth" <<
endl;
1317 *
gmsg <<
"* Field Map read successfully!" <<
endl <<
endl;
1327 *
gmsg <<
"* ----------------------------------------------" <<
endl;
1328 *
gmsg <<
"* READ IN CYCIAE-100 CYCLOTRON FIELD MAP " <<
endl;
1329 *
gmsg <<
"* ----------------------------------------------" <<
endl;
1333 f = std::fopen(
fmapfn_m.c_str(),
"r");
1336 *
gmsg <<
"* Minimal radius of measured field map: " <<
BP_m.
rmin_m <<
" [mm]" <<
endl;
1373 for (
int ii = 0; ii < 13; ii++) {
1376 for (
int k = 0; k < nHalfPoints; k++) {
1390 *
gmsg <<
"* Field Map read successfully!" <<
endl <<
endl;
1400 *
gmsg <<
"* Reading " << fm <<
endl;
1412 std::vector<std::string>::const_iterator fm =
RFfilename_m.begin();
1413 std::vector<std::string>::const_iterator rffcfni =
RFFCoeff_fn_m.begin();
1414 std::vector<std::string>::const_iterator rfvcfni =
RFVCoeff_fn_m.begin();
1421 *
gmsg <<
"* ------------------------------------------------------------" <<
endl;
1422 *
gmsg <<
"* READ IN 3D RF Fields and Frequency Coefficients " <<
endl;
1423 *
gmsg <<
"* ------------------------------------------------------------" <<
endl;
1425 for (; fm !=
RFfilename_m.end(); ++fm, ++rffcfni, ++rfvcfni, ++fcount) {
1432 *
gmsg <<
"RF Frequency Coefficient Filename: " << (*rffcfni) <<
endl;
1434 rffcf = std::fopen((*rffcfni).c_str(),
"r");
1436 if (rffcf == NULL) {
1438 "Cyclotron::getFieldFromFile_Synchrocyclotron",
1439 "failed to open file '" + *rffcfni +
"', please check if it exists");
1442 std::vector<double> fcoeff;
1447 *
gmsg <<
"* Number of coefficients in file: " << nc <<
endl;
1448 for (
int k = 0; k < nc; k++) {
1450 fcoeff.push_back(value);
1453 rffc_m.push_back(fcoeff);
1458 *
gmsg <<
"RF Voltage Coefficient Filename: " << (*rfvcfni) <<
endl;
1460 rfvcf = std::fopen((*rfvcfni).c_str(),
"r");
1461 if (rfvcf == NULL) {
1463 "Cyclotron::getFieldFromFile_Synchrocyclotron",
1464 "failed to open file '" + *rfvcfni +
"', please check if it exists");
1467 std::vector<double> vcoeff;
1470 *
gmsg <<
"* Number of coefficients in file: " << nc <<
endl;
1471 for (
int k = 0; k < nc; k++) {
1473 vcoeff.push_back(value);
1476 rfvc_m.push_back(vcoeff);
1495 fp1.open(fname, std::ios::out);
1513 fp2.open(fname, std::ios::out);
1517 Vector_t tmpE(0.0, 0.0, 0.0), tmpB(0.0, 0.0, 0.0);
1519 Vector_t E(0.0, 0.0, 0.0), B(0.0, 0.0, 0.0);
1520 if (!fi->getFieldstrength(tmpR, tmpE, tmpB)) {
1525 fp2 << tmpR <<
" \t E= " << tmpE <<
"\t B= " << tmpB <<
std::endl;
1528 *
gmsg <<
"\n* Writing 'gnu.out' and 'eb.out' files of cyclotron field maps\n" <<
endl;
1531 *
gmsg <<
"\n* Writing 'gnu.out' file of cyclotron field map\n" <<
endl;
1535 #undef CHECK_CYC_FSCANF_EOF
#define CHECK_CYC_FSCANF_EOF(arg)
Tps< T > cos(const Tps< T > &x)
Cosine.
Tps< T > pow(const Tps< T > &x, int y)
Integer power.
Tps< T > sin(const Tps< T > &x)
Sine.
T::PETE_Expr_t::PETE_Return_t max(const PETE_Expr< T > &expr, NDIndex< D > &loc)
T::PETE_Expr_t::PETE_Return_t min(const PETE_Expr< T > &expr, NDIndex< D > &loc)
PETE_TUTree< FnAbs, typename T::PETE_Expr_t > abs(const PETE_Expr< T > &l)
PETE_TBTree< FnFmod, PETE_Scalar< Vektor< T1, Dim > >, typename T2::PETE_Expr_t > fmod(const Vektor< T1, Dim > &l, const PETE_Expr< T2 > &r)
PETE_TUTree< FnArcTan, typename T::PETE_Expr_t > atan(const PETE_Expr< T > &l)
void fp2(BareField< T, 2U > &field, bool docomm=true)
void fp1(BareField< T, 1U > &field, bool docomm=true)
Inform & endl(Inform &inf)
Inform & level4(Inform &inf)
constexpr double two_pi
The value of.
constexpr double deg2rad
The conversion factor from degrees to radians.
constexpr double e
The value of.
constexpr double pi
The value of.
constexpr double rad2deg
The conversion factor from radians to degrees.
std::string::iterator iterator
T rad(T x)
Convert degrees to radians.
std::string combineFilePath(std::initializer_list< std::string > ilist)
ParticleAttrib< int > Bin
ParticleAttrib< double > M
ParticleAttrib< Vector_t > P
ParticleAttrib< double > Q
ParticleAttrib< short > bunchNum
static OpalData * getInstance()
std::string getAuxiliaryOutputDirectory() const
get the name of the the additional data directory
virtual void visitCyclotron(const Cyclotron &)=0
Apply the algorithm to a cyclotron.
Interface for a single beam element.
PartBunchBase< double, 3 > * RefPartBunch_m
std::vector< double > dbrrt_m
std::vector< double > g3_m
std::vector< double > dbtt_m
std::vector< double > dbrtt_m
std::vector< double > f3_m
std::vector< double > dbrr_m
std::vector< double > dbrrr_m
std::vector< double > dbttt_m
std::vector< double > dbrt_m
std::vector< double > dbr_m
std::vector< double > dbt_m
std::vector< double > bfld_m
std::vector< double > f2_m
std::vector< double > rarr_m
void setTrimCoilThreshold(double)
void setZinit(double zinit)
virtual double getCyclHarm() const
virtual void getDimensions(double &zBegin, double &zEnd) const
std::vector< std::string > RFfilename_m
virtual bool getSpiralFlag() const
virtual bool bends() const
void getFieldFromFile_CYCIAE(const double &scaleFactor)
virtual ElementBase::ElementType getType() const
Get element type std::string.
void setSpiralFlag(bool spiral_flag)
virtual double getPRinit() const
void getFieldFromFile_BandRF(const double &scaleFactor)
void setSuperpose(std::vector< bool > flag)
double getRfFrequ(unsigned int i) const
virtual double getTrimCoilThreshold() const
void setCyclotronType(std::string t)
void setRfPhi(std::vector< double > f)
void setRFVCoeffFN(std::vector< std::string > rfv_coeff_fn)
void setPZinit(double zinit)
virtual void accept(BeamlineVisitor &) const
Apply visitor to Cyclotron.
virtual double getPZinit() const
virtual std::string getFieldMapFN() const
virtual double getRmin() const
void read(const double &scaleFactor)
int idx(int irad, int ktet)
virtual double getMaxR() const
void setEScale(std::vector< double > bs)
void setSymmetry(double symmetry)
std::vector< bool > superpose_m
std::vector< double > escale_m
std::vector< std::string > RFVCoeff_fn_m
std::vector< double > rffrequ_m
void getFieldFromFile_AVFEQ(const double &scaleFactor)
virtual bool getSuperpose(unsigned int i) const
virtual double getRmax() const
virtual double getFMHighE() const
double gutdf5d(double *f, double dx, const int kor, const int krl, const int lpr)
double trimCoilThreshold_m
virtual double getMaxZ() const
virtual double getBScale() const
virtual double getZinit() const
void getFieldFromFile_Carbon(const double &scaleFactor)
virtual double getEScale(unsigned int i) const
void writeOutputFieldFiles()
virtual double getPHIinit() const
double getRfPhi(unsigned int i) const
std::vector< TrimCoil * > trimcoils_m
std::vector< std::string > RFFCoeff_fn_m
void setTrimCoils(const std::vector< TrimCoil * > &trimcoils)
virtual double getSymmetry() const
void setCyclHarm(double h)
virtual double getMinZ() const
void getFieldFromFile_Synchrocyclotron(const double &scaleFactor)
void applyTrimCoil_m(const double r, const double z, const double tet_rad, double *br, double *bz)
Apply trim coils (calculate field contributions)
virtual void initialise(PartBunchBase< double, 3 > *bunch, double &startField, double &endField)
void getFieldFromFile_FFA(const double &scaleFactor)
void initR(double rmin, double dr, int nrad)
void setPHIinit(double phiinit)
void setRfFrequ(std::vector< double > f)
void setPRinit(double prinit)
void applyTrimCoil(const double r, const double z, const double tet_rad, double &br, double &bz)
Apply trim coils (calculate field contributions) with smooth field transition.
virtual bool apply(const size_t &id, const double &t, Vector_t &E, Vector_t &B)
void setRfFieldMapFN(std::vector< std::string > rffmapfn)
bool interpolate(const double &rad, const double &tet_rad, double &br, double &bt, double &bz)
void setFMHighE(double e)
std::vector< Fieldmap * > RFfields_m
std::vector< double > rfphi_m
void setRinit(double rinit)
std::vector< std::vector< double > > rfvc_m
void setFieldMapFN(std::string fmapfn)
std::vector< std::vector< double > > rffc_m
virtual double getFMLowE() const
std::unique_ptr< LossDataSink > lossDs_m
unsigned int getNumberOfTrimcoils() const
virtual double getRinit() const
void getFieldFromFile_Ring(const double &scaleFactor)
virtual double getMinR() const
const std::string & getCyclotronType() const
void setRFFCoeffFN(std::vector< std::string > rff_coeff_fn)
void setBScale(double bs)
virtual const std::string & getName() const
Get element name.
std::string getOutputFN() const
Get output filename.
static Fieldmap * getFieldmap(std::string Filename, bool fast=false)
static void readMap(std::string Filename)
Vektor< double, 3 > Vector_t