33 #include <boost/filesystem.hpp>
41 namespace fs = boost::filesystem;
43 #define REGISTER_PARSE_TYPE(X) template <> struct Fieldmap::TypeParseTraits<X> \
44 { static const char* name; } ; const char* Fieldmap::TypeParseTraits<X>::name = #X
49 (*position).second.RefCounter++;
50 return (*position).second.Map;
70 return (*position.first).second.Map;
87 return (*position.first).second.Map;
104 return (*position.first).second.Map;
121 return (*position.first).second.Map;
138 return (*position.first).second.Map;
155 return (*position.first).second.Map;
163 return (*position.first).second.Map;
172 return (*position.first).second.Map;
181 return (*position.first).second.Map;
190 return (*position.first).second.Map;
199 return (*position.first).second.Map;
208 return (*position.first).second.Map;
216 return (*position.first).second.Map;
224 return (*position.first).second.Map;
232 return (*position.first).second.Map;
249 return (*position.first).second.Map;
254 "Couldn't determine type of fieldmap in file \"" + Filename +
"\"");
260 std::vector<std::string> name_list;
262 name_list.push_back((*it).first);
274 delete it->second.Map;
275 it->second.Map = NULL;
281 char magicnumber[5] =
" ";
286 if(Filename ==
"1DPROFILE1-DEFAULT")
289 if (!fs::exists(Filename))
291 "File \"" + Filename +
"\" doesn't exist");
293 std::ifstream File(Filename.c_str());
295 std::cerr <<
"could not open file " << Filename <<
std::endl;
299 getLine(File, lines_read_m, buffer);
300 std::istringstream interpreter(buffer, std::istringstream::in);
302 interpreter.read(magicnumber, 4);
304 if(strcmp(magicnumber,
"3DDy") == 0)
307 if(strcmp(magicnumber,
"3DMa") == 0) {
308 char tmpString[21] =
" ";
309 interpreter.read(tmpString, 20);
311 if(strcmp(tmpString,
"gnetoStatic_Extended") == 0)
317 if(strcmp(magicnumber,
"3DEl") == 0)
320 if(strcmp(magicnumber,
"2DDy") == 0) {
326 if(strcmp(magicnumber,
"2DMa") == 0) {
332 if(strcmp(magicnumber,
"2DEl") == 0) {
338 if(strcmp(magicnumber,
"1DDy") == 0)
341 if(strcmp(magicnumber,
"1DMa") == 0)
344 if(strcmp(magicnumber,
"1DPr") == 0) {
353 if(strcmp(magicnumber,
"1DEl") == 0)
356 if(strcmp(magicnumber,
"\211HDF") == 0) {
363 h5_size_t len_name =
sizeof(
name);
364 h5_prop_t props = H5CreateFileProp ();
366 h5err = H5SetPropFileMPIOCollective (props, &comm);
367 assert (h5err != H5_ERR);
368 h5_file_t file = H5OpenFile (Filename.c_str(), H5_O_RDONLY, props);
369 assert (file != (h5_file_t)H5_ERR);
372 h5err = H5SetStep(file, 0);
373 assert (h5err != H5_ERR);
375 h5_int64_t num_fields = H5BlockGetNumFields(file);
376 assert (num_fields != H5_ERR);
379 for(h5_ssize_t i = 0; i < num_fields; ++ i) {
380 h5err = H5BlockGetFieldInfo(
381 file, (h5_size_t)i, name, len_name, NULL, NULL, NULL, NULL);
382 assert (h5err != H5_ERR);
384 if(strcmp(name,
"Bfield") == 0) {
387 }
else if(strcmp(name,
"Hfield") == 0) {
392 h5err = H5CloseFile(file);
393 assert (h5err != H5_ERR);
398 if(strcmp(magicnumber,
"Astr") == 0) {
399 char tmpString[3] =
" ";
400 interpreter.read(tmpString, 2);
401 if(strcmp(tmpString,
"aE") == 0) {
404 if(strcmp(tmpString,
"aM") == 0) {
407 if(strcmp(tmpString,
"aD") == 0) {
419 if(!(*position).second.read) {
420 (*position).second.Map->readMap();
421 (*position).second.read =
true;
431 if((*position).second.RefCounter > 0) {
432 (*position).second.RefCounter--;
435 if ((*position).second.RefCounter == 0) {
436 delete (*position).second.Map;
437 (*position).second.Map = NULL;
444 std::pair<double, double> fieldDimensions,
446 const std::vector<double> &fourierCoefficients,
447 gsl_spline *splineCoefficients,
448 gsl_interp_accel *splineAccelerator) {
449 double length = fieldDimensions.second - fieldDimensions.first;
450 unsigned int sizeSampling =
std::floor(length / deltaZ + 0.5);
451 std::vector<double> zSampling(sizeSampling);
452 zSampling[0] = fieldDimensions.first;
453 for (
unsigned int i = 1; i < sizeSampling; ++ i) {
454 zSampling[i] = zSampling[i-1] + deltaZ;
456 checkMap(accuracy, length, zSampling, fourierCoefficients, splineCoefficients, splineAccelerator);
461 const std::vector<double> &zSampling,
462 const std::vector<double> &fourierCoefficients,
463 gsl_spline *splineCoefficients,
464 gsl_interp_accel *splineAccelerator) {
466 double maxDiff = 0.0;
468 double ezSquare = 0.0;
469 size_t lastDot =
Filename_m.find_last_of(
".");
470 size_t lastSlash =
Filename_m.find_last_of(
"/");
471 lastSlash = (lastSlash == std::string::npos)? 0: lastSlash + 1;
476 out.open(
"data/" +
Filename_m.substr(lastSlash, lastDot) +
".check");
477 out <<
"# z original reproduced\n";
479 auto it = zSampling.begin();
480 auto end = zSampling.end();
481 for (; it != end; ++ it) {
483 double onAxisFieldCheck = fourierCoefficients[0];
485 for(
unsigned int l = 1; l < accuracy; ++l, n += 2) {
486 double coskzl =
cos(kz * l);
487 double sinkzl =
sin(kz * l);
489 onAxisFieldCheck += (fourierCoefficients[
n] * coskzl - fourierCoefficients[n + 1] * sinkzl);
491 double ez = gsl_spline_eval(splineCoefficients, *it, splineAccelerator);
492 double difference =
std::abs(ez - onAxisFieldCheck);
493 maxDiff = difference > maxDiff? difference: maxDiff;
499 out << std::setw(16) << std::setprecision(8) << *it
500 << std::setw(16) << std::setprecision(8) << ez
501 << std::setw(16) << std::setprecision(8) << onAxisFieldCheck
507 if (
sqrt(error / ezSquare) > 1
e-1 || maxDiff > 1
e-1 * ezMax) {
511 "Field map can't be reproduced properly with the given number of fourier components");
513 if (
sqrt(error / ezSquare) > 1
e-2 || maxDiff > 1
e-2 * ezMax) {
535 comment = buffer.find(
"#");
536 buffer = buffer.substr(0, comment);
540 }
while(!in.eof() && lastof == std::string::npos);
542 if(firstof != std::string::npos) {
543 buffer = buffer.substr(firstof, lastof - firstof + 1);
552 size_t comment = buffer.find_first_of(
"#");
553 buffer = buffer.substr(0, comment);
555 if(lasto != std::string::npos) {
564 const std::string &expecting,
565 const std::string &found) {
566 std::stringstream errormsg;
567 errormsg <<
"THERE SEEMS TO BE SOMETHING WRONG WITH YOUR FIELD MAP '" <<
Filename_m <<
"'.\n"
568 <<
"expecting: '" << expecting <<
"' on line " <<
lines_read_m <<
",\n"
569 <<
"found instead: '" << found <<
"'.";
576 const bool &read_all,
577 const std::string &expecting,
578 const std::string &found) {
579 std::string error_msg;
581 error_msg = std::string(
"Didn't find enough values!");
582 }
else if(state & std::ios_base::eofbit) {
583 error_msg = std::string(
"Found more values than expected!");
584 }
else if(state & std::ios_base::failbit) {
585 error_msg = std::string(
"Found wrong type of values!");
591 std::stringstream errormsg;
592 errormsg <<
"THERE SEEMS TO BE SOMETHING WRONG WITH YOUR FIELD MAP '" <<
Filename_m <<
"'.\n"
593 <<
"There are only " <<
lines_read_m - 1 <<
" lines in the file, expecting more.\n"
594 <<
"Please check the section about field maps in the user manual.";
601 std::stringstream errormsg;
602 errormsg <<
"THERE SEEMS TO BE SOMETHING WRONG WITH YOUR FIELD MAP '" <<
Filename_m <<
"'.\n"
603 <<
"There are too many lines in the file, expecting only " <<
lines_read_m <<
" lines.\n"
604 <<
"Please check the section about field maps in the user manual.";
611 std::stringstream errormsg;
612 errormsg <<
"DISABLING FIELD MAP '" +
Filename_m +
"' DUE TO PARSING ERRORS." ;
619 std::stringstream errormsg;
620 errormsg <<
"DISABLING FIELD MAP '" <<
Filename_m <<
"' SINCE FILE COULDN'T BE FOUND!";
627 std::stringstream errormsg;
628 errormsg <<
"IT SEEMS THAT YOU USE TOO FEW FOURIER COMPONENTS TO SUFFICIENTLY WELL\n"
629 <<
"RESOLVE THE FIELD MAP '" <<
Filename_m <<
"'.\n"
630 <<
"PLEASE INCREASE THE NUMBER OF FOURIER COMPONENTS!\n"
631 <<
"The ratio (e_i - E_i)^2 / E_i^2 is " << std::to_string(squareError) <<
" and\n"
632 <<
"the ratio (max_i(|e_i - E_i|) / max_i(|E_i|) is " << std::to_string(maxError) <<
".\n"
633 <<
"Here E_i is the field as in the field map and e_i is the reconstructed field.\n"
634 <<
"The lower limit for the two ratios is 1e-2\n"
635 <<
"Have a look into the directory data/ for a reconstruction of the field map.\n";
636 std::string errormsg_str =
typeset_msg(errormsg.str(),
"warning");
641 std::ofstream omsg(
"errormsg.txt", std::ios_base::app);
648 static std::string frame(
"* ******************************************************************************\n");
649 static unsigned int frame_width = frame.length() - 5;
650 static std::string closure(
" *\n");
652 std::string return_string(
"\n" + frame);
654 int remaining_length = msg.length();
655 unsigned int current_position = 0;
658 for(; ii < title.length(); ++ ii) {
661 return_string.replace(15 + 2 * ii, 1,
" ");
662 return_string.replace(16 + 2 * ii, 1, &c, 1);
664 return_string.replace(15 + 2 * ii, 1,
" ");
666 while(remaining_length > 0) {
667 size_t eol = msg.find(
"\n", current_position);
668 std::string next_to_process;
669 if(eol != std::string::npos) {
670 next_to_process = msg.substr(current_position, eol - current_position);
672 next_to_process = msg.substr(current_position);
676 if(eol - current_position < frame_width) {
677 return_string +=
"* " + next_to_process + closure.substr(eol - current_position + 2);
679 unsigned int last_space = next_to_process.rfind(
" ", frame_width);
681 if(last_space < frame_width) {
682 return_string +=
"* " + next_to_process.substr(0, last_space) + closure.substr(last_space + 2);
684 return_string +=
"* " + next_to_process.substr(0, last_space) +
" *\n";
686 if(next_to_process.length() - last_space + 1 < frame_width) {
687 return_string +=
"* " + next_to_process.substr(last_space + 1) + closure.substr(next_to_process.length() - last_space + 1);
689 return_string +=
"* " + next_to_process.substr(last_space + 1) +
" *\n";
692 return_string +=
"* " + next_to_process +
" *\n";
696 current_position = eol + 1;
697 remaining_length = msg.length() - current_position;
699 return_string += frame;
701 return return_string;
708 std::vector<double> &engeCoeffsExit) {
713 double &entranceParameter2,
714 double &entranceParameter3) {
719 double &exitParameter2,
720 double &exitParameter3) {
735 const std::pair<double, double> &xrange,
736 const std::pair<double, double> &yrange,
737 const std::pair<double, double> &zrange,
738 const std::vector<Vector_t> &ef,
739 const std::vector<Vector_t> &bf) {
740 const size_t numpoints = nx * ny * nz;
742 (ef.size() != numpoints && bf.size() != numpoints))
return;
744 size_t extensionStart =
Filename_m.find_last_of(
'.');
746 of.open (std::string (
"data/" +
Filename_m.substr(0,extensionStart) +
".vtk").c_str ());
747 assert (of.is_open ());
750 const double hx = (xrange.second - xrange.first) / (nx - 1);
751 const double hy = (yrange.second - yrange.first) / (ny - 1);
752 const double hz = (zrange.second - zrange.first) / (nz - 1);
754 of <<
"# vtk DataFile Version 2.0" <<
std::endl;
755 of <<
"generated by 3D fieldmaps" <<
std::endl;
757 of <<
"DATASET RECTILINEAR_GRID" <<
std::endl;
758 of <<
"DIMENSIONS " << nx <<
" " << ny <<
" " << nz <<
std::endl;
760 of <<
"X_COORDINATES " << nx <<
" float" <<
std::endl;
762 for (
unsigned int i = 1; i < nx - 1; ++ i) {
763 of <<
" " << xrange.first + i * hx;
767 of <<
"Y_COORDINATES " << ny <<
" float" <<
std::endl;
769 for (
unsigned int i = 1; i < ny - 1; ++ i) {
770 of <<
" " << yrange.first + i * hy;
774 of <<
"Z_COORDINATES " << nz <<
" float" <<
std::endl;
776 for (
unsigned int i = 1; i < nz - 1; ++ i) {
777 of <<
" " << zrange.first + i * hz;
781 of <<
"POINT_DATA " << numpoints <<
std::endl;
783 if (ef.size() == numpoints) {
784 of <<
"VECTORS EField float" <<
std::endl;
786 for (
size_t i = 0; i < numpoints; ++ i) {
787 of << ef[i](0) <<
" " << ef[i](1) <<
" " << ef[i](2) << std::endl;
792 if (bf.size() == numpoints) {
793 of <<
"VECTORS BField float" <<
std::endl;
795 for (
size_t i = 0; i < numpoints; ++ i) {
796 of << bf[i](0) <<
" " << bf[i](1) <<
" " << bf[i](2) << std::endl;
808 std::map<std::string, Fieldmap::FieldmapDescription>
Fieldmap::FieldmapDictionary = std::map<std::string, Fieldmap::FieldmapDescription>();
PETE_TUTree< FnAbs, typename T::PETE_Expr_t > abs(const PETE_Expr< T > &l)
constexpr double e
The value of .
static void clearDictionary()
static MapType readHeader(std::string Filename)
virtual void getOnaxisEz(std::vector< std::pair< double, double > > &onaxis)
Tps< T > sin(const Tps< T > &x)
Sine.
virtual double getFieldGap()
constexpr double two_pi
The value of .
void disableFieldmapWarning()
static std::string typeset_msg(const std::string &msg, const std::string &title)
static Fieldmap * getFieldmap(std::string Filename, bool fast=false)
static void deleteFieldmap(std::string Filename)
static OpalData * getInstance()
static char buffer_m[256]
void missingValuesWarning()
virtual void setFieldGap(double gap)
constexpr double c
The velocity of light in m/s.
static std::map< std::string, FieldmapDescription > FieldmapDictionary
bool interpreteEOF(std::ifstream &in)
virtual void setFieldLength(const double &)
Tps< T > pow(const Tps< T > &x, int y)
Integer power.
virtual void get1DProfile1EntranceParam(double &entranceParameter1, double &entranceParameter2, double &entranceParameter3)
static MPI_Comm getComm()
void interpreteWarning(const std::string &error_msg, const std::string &expecting, const std::string &found)
virtual void get1DProfile1EngeCoeffs(std::vector< double > &engeCoeffsEntry, std::vector< double > &engeCoeffsExit)
Tps< T > sqrt(const Tps< T > &x)
Square root.
static std::vector< std::string > getListFieldmapNames()
virtual void setEdgeConstants(const double &bendAngle, const double &entranceAngle, const double &exitAngle)
static std::string alpha_numeric
#define REGISTER_PARSE_TYPE(X)
Tps< T > cos(const Tps< T > &x)
Cosine.
void getLine(std::ifstream &in, std::string &buffer)
void exceedingValuesWarning()
std::string::iterator iterator
void checkMap(unsigned int accuracy, std::pair< double, double > fieldDimensions, double deltaZ, const std::vector< double > &fourierCoefficients, gsl_spline *splineCoefficients, gsl_interp_accel *splineAccelerator)
void lowResolutionWarning(double squareError, double maxError)
#define READ_BUFFER_LENGTH
void write3DField(unsigned int nx, unsigned int ny, unsigned int nz, const std::pair< double, double > &xrange, const std::pair< double, double > &yrange, const std::pair< double, double > &zrange, const std::vector< Vector_t > &ef, const std::vector< Vector_t > &bf)
PETE_TUTree< FnFloor, typename T::PETE_Expr_t > floor(const PETE_Expr< T > &l)
Inform & endl(Inform &inf)
virtual void get1DProfile1ExitParam(double &exitParameter1, double &exitParameter2, double &exitParameter3)