Fieldmap.cpp 27.4 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
#include "Fields/Fieldmap.h"
#include "Fields/FM3DDynamic.h"
#include "Fields/FM3DH5Block.h"
#include "Fields/FM3DH5Block_nonscale.h"
#include "Fields/FM3DMagnetoStaticH5Block.h"
#include "Fields/FM2DDynamic.h"
#include "Fields/FM2DElectroStatic.h"
#include "Fields/FM2DMagnetoStatic.h"
#include "Fields/FM1DDynamic.h"
#include "Fields/FM1DDynamic_fast.h"
#include "Fields/Astra1DDynamic.h"
#include "Fields/Astra1DDynamic_fast.h"
#include "Fields/FM1DElectroStatic.h"
#include "Fields/FM1DElectroStatic_fast.h"
#include "Fields/Astra1DElectroStatic.h"
#include "Fields/Astra1DElectroStatic_fast.h"
#include "Fields/FM1DMagnetoStatic.h"
#include "Fields/FM1DMagnetoStatic_fast.h"
#include "Fields/Astra1DMagnetoStatic.h"
#include "Fields/Astra1DMagnetoStatic_fast.h"
#include "Fields/FM1DProfile1.h"
#include "Fields/FM1DProfile2.h"
#include "Fields/FMDummy.h"
#include "Utilities/GeneralClassicException.h"
kraus's avatar
kraus committed
25
#include "Utilities/Options.h"
26 27
#include "Physics/Physics.h"

gsell's avatar
gsell committed
28 29
#include "H5hut.h"

kraus's avatar
kraus committed
30 31 32 33
#include <iostream>
#include <fstream>
#include <ios>

gsell's avatar
gsell committed
34 35 36
#define REGISTER_PARSE_TYPE(X) template <> struct Fieldmap::TypeParseTraits<X> \
    { static const char* name; } ; const char* Fieldmap::TypeParseTraits<X>::name = #X

Steve Russell's avatar
Steve Russell committed
37
Fieldmap *Fieldmap::getFieldmap(std::string Filename, bool fast) {
38
    std::map<std::string, FieldmapDescription>::iterator position = FieldmapDictionary.find(Filename);
gsell's avatar
gsell committed
39 40 41 42 43
    if(position != FieldmapDictionary.end()) {
        (*position).second.RefCounter++;
        return (*position).second.Map;
    } else {
        MapType type;
44
        std::pair<std::map<std::string, FieldmapDescription>::iterator, bool> position;
gsell's avatar
gsell committed
45 46
        type = readHeader(Filename);
        switch(type) {
47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191
        case T1DDynamic:
            if(fast) {
                position = FieldmapDictionary.insert(std::make_pair(Filename,
                                                                    FieldmapDescription(T1DDynamic,
                                                                                        new FM1DDynamic_fast(Filename))));
            } else {
                position = FieldmapDictionary.insert(std::make_pair(Filename,
                                                                    FieldmapDescription(T1DDynamic,
                                                                                        new FM1DDynamic(Filename))));
            }
            return (*position.first).second.Map;
            break;

        case TAstraDynamic:
            if(fast) {
                position = FieldmapDictionary.insert(std::make_pair(Filename,
                                                                    FieldmapDescription(TAstraDynamic,
                                                                                        new Astra1DDynamic_fast(Filename))));
            } else {
                position = FieldmapDictionary.insert(std::make_pair(Filename,
                                                                    FieldmapDescription(TAstraDynamic,
                                                                                        new Astra1DDynamic(Filename))));
            }
            return (*position.first).second.Map;
            break;

        case T1DElectroStatic:
            if(fast) {
                position = FieldmapDictionary.insert(std::make_pair(Filename,
                                                                    FieldmapDescription(T1DElectroStatic,
                                                                                        new FM1DElectroStatic_fast(Filename))));
            } else {
                position = FieldmapDictionary.insert(std::make_pair(Filename,
                                                                    FieldmapDescription(T1DElectroStatic,
                                                                                        new FM1DElectroStatic(Filename))));
            }
            return (*position.first).second.Map;
            break;

        case TAstraElectroStatic:
            if(fast) {
                position = FieldmapDictionary.insert(std::make_pair(Filename,
                                                                    FieldmapDescription(TAstraElectroStatic,
                                                                                        new Astra1DElectroStatic_fast(Filename))));
            } else {
                position = FieldmapDictionary.insert(std::make_pair(Filename,
                                                                    FieldmapDescription(TAstraElectroStatic,
                                                                                        new Astra1DElectroStatic(Filename))));
            }
            return (*position.first).second.Map;
            break;

        case T1DMagnetoStatic:
            if(fast) {
                position = FieldmapDictionary.insert(std::make_pair(Filename,
                                                                    FieldmapDescription(T1DMagnetoStatic,
                                                                                        new FM1DMagnetoStatic_fast(Filename))));
            } else {
                position = FieldmapDictionary.insert(std::make_pair(Filename,
                                                                    FieldmapDescription(T1DMagnetoStatic,
                                                                                        new FM1DMagnetoStatic(Filename))));
            }
            return (*position.first).second.Map;
            break;

        case TAstraMagnetoStatic:
            if(fast) {
                position = FieldmapDictionary.insert(std::make_pair(Filename,
                                                                    FieldmapDescription(TAstraMagnetoStatic,
                                                                                        new Astra1DMagnetoStatic_fast(Filename))));
            } else {
                position = FieldmapDictionary.insert(std::make_pair(Filename,
                                                                    FieldmapDescription(TAstraMagnetoStatic,
                                                                                        new Astra1DMagnetoStatic(Filename))));
            }
            return (*position.first).second.Map;
            break;

        case T1DProfile1:
            position = FieldmapDictionary.insert(std::make_pair(Filename,
                                                                FieldmapDescription(T1DProfile1,
                                                                                    new FM1DProfile1(Filename))));
            return (*position.first).second.Map;
            break;

        case T1DProfile2:
            position = FieldmapDictionary.insert(std::make_pair(Filename,
                                                                FieldmapDescription(T1DProfile2,
                                                                                    new FM1DProfile2(Filename))));
            return (*position.first).second.Map;
            break;

        case T2DDynamic:
            position = FieldmapDictionary.insert(std::make_pair(Filename,
                                                                FieldmapDescription(T2DDynamic,
                                                                                    new FM2DDynamic(Filename))));
            return (*position.first).second.Map;
            break;

        case T2DElectroStatic:
            position = FieldmapDictionary.insert(std::make_pair(Filename,
                                                                FieldmapDescription(T2DElectroStatic,
                                                                                    new FM2DElectroStatic(Filename))));
            return (*position.first).second.Map;
            break;

        case T2DMagnetoStatic:
            position = FieldmapDictionary.insert(std::make_pair(Filename,
                                                                FieldmapDescription(T2DMagnetoStatic,
                                                                                    new FM2DMagnetoStatic(Filename))));
            return (*position.first).second.Map;
            break;

        case T3DDynamic:
            position = FieldmapDictionary.insert(std::make_pair(Filename,
                                                                FieldmapDescription(T3DDynamic,
                                                                                    new FM3DDynamic(Filename))));
            return (*position.first).second.Map;
            break;

        case T3DMagnetoStaticH5Block:
            position = FieldmapDictionary.insert(std::make_pair(Filename,
                                                                FieldmapDescription(T3DMagnetoStaticH5Block,
                                                                                    new FM3DMagnetoStaticH5Block(Filename))));
            return (*position.first).second.Map;
            break;

        case T3DDynamicH5Block:
            if(fast) {
                position = FieldmapDictionary.insert(std::make_pair(Filename,
                                                                    FieldmapDescription(T3DDynamic,
                                                                                        new FM3DH5Block_nonscale(Filename))));
            } else {
                position = FieldmapDictionary.insert(std::make_pair(Filename,
                                                                    FieldmapDescription(T3DDynamic,
                                                                                        new FM3DH5Block(Filename))));
            }
            return (*position.first).second.Map;
            break;

        default:
            position = FieldmapDictionary.insert(std::make_pair(Filename,
                                                                FieldmapDescription(UNKNOWN,
                                                                                    new FMDummy(Filename))));
            return (*position.first).second.Map;
gsell's avatar
gsell committed
192 193 194 195
        }
    }
}

196 197 198
std::vector<std::string> Fieldmap::getListFieldmapNames() {
    std::vector<std::string> name_list;
    for(std::map<std::string, FieldmapDescription>::const_iterator it = FieldmapDictionary.begin(); it != FieldmapDictionary.end(); ++ it) {
gsell's avatar
gsell committed
199 200 201 202 203
        name_list.push_back((*it).first);
    }
    return name_list;
}

Steve Russell's avatar
Steve Russell committed
204
void Fieldmap::deleteFieldmap(std::string Filename) {
205
    freeMap(Filename);
gsell's avatar
gsell committed
206 207
}

208 209 210 211 212 213 214 215 216
void Fieldmap::clearDictionary() {
    std::map<std::string, FieldmapDescription>::iterator it = FieldmapDictionary.begin();
    for (;it != FieldmapDictionary.end(); ++ it) {
        delete it->second.Map;
        it->second.Map = NULL;
    }
    FieldmapDictionary.clear();
}

Steve Russell's avatar
Steve Russell committed
217
MapType Fieldmap::readHeader(std::string Filename) {
gsell's avatar
gsell committed
218
    char magicnumber[5] = "    ";
Steve Russell's avatar
Steve Russell committed
219
    std::string buffer;
gsell's avatar
gsell committed
220 221 222 223 224 225
    int lines_read_m = 0;

    // Check for default map(s).
    if(Filename == "1DPROFILE1-DEFAULT")
        return T1DProfile1;

226
    std::ifstream File(Filename.c_str());
gsell's avatar
gsell committed
227
    if(!File.good()) {
228
        std::cerr << "could not open file " << Filename << std::endl;
gsell's avatar
gsell committed
229 230 231 232
        return UNKNOWN;
    }

    getLine(File, lines_read_m, buffer);
233
    std::istringstream interpreter(buffer, std::istringstream::in);
gsell's avatar
gsell committed
234 235 236 237 238 239 240 241 242 243 244 245 246 247 248

    interpreter.read(magicnumber, 4);

    if(strcmp(magicnumber, "3DDy") == 0)
        return T3DDynamic;

    if(strcmp(magicnumber, "3DMa") == 0)
        return T3DMagnetoStatic;

    if(strcmp(magicnumber, "3DEl") == 0)
        return T3DElectroStatic;

    if(strcmp(magicnumber, "2DDy") == 0) {
        char tmpString[14] = "             ";
        interpreter.read(tmpString, 13);
249
        return T2DDynamic;
gsell's avatar
gsell committed
250 251 252 253 254
    }

    if(strcmp(magicnumber, "2DMa") == 0) {
        char tmpString[20] = "                   ";
        interpreter.read(tmpString, 19);
255
        return T2DMagnetoStatic;
gsell's avatar
gsell committed
256 257 258 259 260
    }

    if(strcmp(magicnumber, "2DEl") == 0) {
        char tmpString[20] = "                   ";
        interpreter.read(tmpString, 19);
261
        return T2DElectroStatic;
gsell's avatar
gsell committed
262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287
    }

    if(strcmp(magicnumber, "1DDy") == 0)
        return T1DDynamic;

    if(strcmp(magicnumber, "1DMa") == 0)
        return T1DMagnetoStatic;

    if(strcmp(magicnumber, "1DPr") == 0) {
        char tmpString[7] = "      ";
        interpreter.read(tmpString, 6);
        if(strcmp(tmpString, "ofile1") == 0)
            return T1DProfile1;
        if(strcmp(tmpString, "ofile2") == 0)
            return T1DProfile2;
    }

    if(strcmp(magicnumber, "1DEl") == 0)
        return T1DElectroStatic;

    if(strcmp(magicnumber, "\211HDF") == 0) {
        h5_err_t h5err;
        h5_size_t grid_rank;
        h5_size_t grid_dims[3];
        h5_size_t field_dims;
        char name[20];
288
        h5_size_t len_name = sizeof(name);
gsell's avatar
gsell committed
289 290
        h5_int64_t ftype;

291
        h5_file_t *file = H5OpenFile(Filename.c_str(), H5_O_RDONLY, Ippl::getComm());
292
        if(file != (void*)H5_ERR) {
gsell's avatar
gsell committed
293 294 295 296
            h5err = H5SetStep(file, 0);
            if(h5err != H5_SUCCESS)
                ERRORMSG("H5 rc= " << h5err << " in " << __FILE__ << " @ line " << __LINE__ << endl);
            h5_int64_t num_fields = H5BlockGetNumFields(file);
gsell's avatar
gsell committed
297
            MapType maptype = UNKNOWN;
gsell's avatar
gsell committed
298 299 300 301
            for(h5_ssize_t i = 0; i < num_fields; ++ i) {
                h5err = H5BlockGetFieldInfo(file, (h5_size_t)i, name, len_name, &grid_rank, grid_dims, &field_dims, &ftype);
                if(h5err != H5_SUCCESS)
                    ERRORMSG("H5 rc= " << h5err << " in " << __FILE__ << " @ line " << __LINE__ << endl);
gsell's avatar
gsell committed
302
                // using field name "Bfield" and "Hfield" to distinguish the type
303 304
                if(strcmp(name, "Bfield") == 0) {
                    maptype = T3DMagnetoStaticH5Block;
gsell's avatar
gsell committed
305
                    break;
gsell's avatar
gsell committed
306
                } else if(strcmp(name, "Hfield") == 0) {
307
                    maptype = T3DDynamicH5Block;
gsell's avatar
gsell committed
308
                    break;
gsell's avatar
gsell committed
309 310 311 312 313
                }
            }
            h5err = H5CloseFile(file);
            if(h5err != H5_SUCCESS)
                ERRORMSG("H5 rc= " << h5err << " in " << __FILE__ << " @ line " << __LINE__ << endl);
314 315
            if(maptype != UNKNOWN)
                return maptype;
gsell's avatar
gsell committed
316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335
        }
    }
    if(strcmp(magicnumber, "Astr") == 0) {
        char tmpString[3] = "  ";
        interpreter.read(tmpString, 2);
        if(strcmp(tmpString, "aE") == 0) {
            return TAstraElectroStatic;
        }
        if(strcmp(tmpString, "aM") == 0) {
            return TAstraMagnetoStatic;
        }
        if(strcmp(tmpString, "aD") == 0) {
            return TAstraDynamic;
        }
    }


    return UNKNOWN;
}

Steve Russell's avatar
Steve Russell committed
336
void Fieldmap::readMap(std::string Filename) {
337
    std::map<std::string, FieldmapDescription>::iterator position = FieldmapDictionary.find(Filename);
gsell's avatar
gsell committed
338 339 340 341 342
    if(position != FieldmapDictionary.end())
        if(!(*position).second.read)
            (*position).second.Map->readMap();
}

Steve Russell's avatar
Steve Russell committed
343
void Fieldmap::freeMap(std::string Filename) {
344
    std::map<std::string, FieldmapDescription>::iterator position = FieldmapDictionary.find(Filename);
345 346 347
    /*
      FIXME: find( ) make problem, crashes
    */
gsell's avatar
gsell committed
348
    if(position != FieldmapDictionary.end()) {
349 350 351 352 353 354 355
        if((*position).second.RefCounter > 0) {
            (*position).second.RefCounter--;
        }

        if ((*position).second.RefCounter == 0) {
            delete (*position).second.Map;
            (*position).second.Map = NULL;
gsell's avatar
gsell committed
356 357 358 359 360
            FieldmapDictionary.erase(position);
        }
    }
}

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
void Fieldmap::checkMap(unsigned int accuracy,
                        std::pair<double, double> fieldDimensions,
                        double deltaZ,
                        const std::vector<double> &fourierCoefficients,
                        gsl_spline *splineCoefficients,
                        gsl_interp_accel *splineAccelerator) {
    double length = fieldDimensions.second - fieldDimensions.first;
    unsigned int sizeSampling = std::floor(length / deltaZ + 0.5);
    std::vector<double> zSampling(sizeSampling);
    zSampling[0] = fieldDimensions.first;
    for (unsigned int i = 1; i < sizeSampling; ++ i) {
        zSampling[i] = zSampling[i-1] + deltaZ;
    }
    checkMap(accuracy, length, zSampling, fourierCoefficients, splineCoefficients, splineAccelerator);
}

void Fieldmap::checkMap(unsigned int accuracy,
                        double length,
                        const std::vector<double> &zSampling,
                        const std::vector<double> &fourierCoefficients,
                        gsl_spline *splineCoefficients,
                        gsl_interp_accel *splineAccelerator) {
    double error = 0.0;
    double maxDiff = 0.0;
    double ezMax = 0.0;
    double ezSquare = 0.0;
    size_t lastDot = Filename_m.find_last_of(".");
    std::ofstream out;
    if (Ippl::myNode() == 0) {
        out.open("data/" + Filename_m.substr(0, lastDot) + ".check");
        out << "# z  original reproduced\n";
    }
    auto it = zSampling.begin();
    auto end = zSampling.end();
    for (; it != end; ++ it) {
        const double kz = Physics::two_pi * (*it / length + 0.5);
        double onAxisFieldCheck = fourierCoefficients[0];
        unsigned int n = 1;
        for(unsigned int l = 1; l < accuracy; ++l, n += 2) {
            double coskzl = cos(kz * l);
            double sinkzl = sin(kz * l);

            onAxisFieldCheck += (fourierCoefficients[n] * coskzl - fourierCoefficients[n + 1] * sinkzl);
        }
        double ez = gsl_spline_eval(splineCoefficients, *it, splineAccelerator);
        double difference = std::abs(ez - onAxisFieldCheck);
        maxDiff = difference > maxDiff? difference: maxDiff;
        ezMax = std::abs(ez) > ezMax? std::abs(ez): ezMax;
        error += std::pow(difference, 2.0);
        ezSquare += std::pow(ez, 2.0);
        out << std::setw(16) << std::setprecision(8) << *it
            << std::setw(16) << std::setprecision(8) << ez
            << std::setw(16) << std::setprecision(8) << onAxisFieldCheck
            << std::endl;
    }
    out.close();

    if (sqrt(error / ezSquare) > 1e-1 || maxDiff > 1e-1 * ezMax) {
        lowResolutionWarning(sqrt(error / ezSquare), maxDiff / ezMax);

        throw GeneralClassicException("Astra2DDynamic_fast::readMap()",
                                      "Field map can't be reproduced properly with the given number of fourier components");
    }
    if (sqrt(error / ezSquare) > 1e-2 || maxDiff > 1e-2 * ezMax) {
        lowResolutionWarning(sqrt(error / ezSquare), maxDiff / ezMax);
    }
}

gsell's avatar
gsell committed
429 430 431 432 433 434
void Fieldmap::setEdgeConstants(const double &bendAngle, const double &entranceAngle, const double &exitAngle)
{};

void Fieldmap::setFieldLength(const double &)
{};

435
void Fieldmap::getLine(std::ifstream &in, int &lines_read, std::string &buffer) {
gsell's avatar
gsell committed
436 437 438 439 440 441 442 443
    size_t firstof = 0;
    size_t lastof;
    size_t comment;

    do {
        ++ lines_read;
        in.getline(buffer_m, READ_BUFFER_LENGTH);

Steve Russell's avatar
Steve Russell committed
444
        buffer = std::string(buffer_m);
gsell's avatar
gsell committed
445 446 447 448 449 450

        comment = buffer.find("#");
        buffer = buffer.substr(0, comment);

        lastof = buffer.find_last_of(alpha_numeric);
        firstof = buffer.find_first_of(alpha_numeric);
Steve Russell's avatar
Steve Russell committed
451
    } while(!in.eof() && lastof == std::string::npos);
gsell's avatar
gsell committed
452

Steve Russell's avatar
Steve Russell committed
453
    if(firstof != std::string::npos) {
gsell's avatar
gsell committed
454 455 456 457
        buffer = buffer.substr(firstof, lastof - firstof + 1);
    }
}

458
bool Fieldmap::interpreteEOF(std::ifstream &in) {
gsell's avatar
gsell committed
459 460 461
    while(!in.eof()) {
        ++lines_read_m;
        in.getline(buffer_m, READ_BUFFER_LENGTH);
Steve Russell's avatar
Steve Russell committed
462 463
        std::string buffer(buffer_m);
        std::string rest;
gsell's avatar
gsell committed
464 465 466
        size_t comment = buffer.find_first_of("#");
        buffer = buffer.substr(0, comment);
        size_t lasto = buffer.find_first_of(alpha_numeric);
Steve Russell's avatar
Steve Russell committed
467
        if(lasto != std::string::npos) {
gsell's avatar
gsell committed
468 469 470 471 472 473 474
            exceedingValuesWarning();
            return false;
        }
    }
    return true;
}

Steve Russell's avatar
Steve Russell committed
475 476 477 478 479
void Fieldmap::interpreteWarning(const std::string &error_msg,
                                 const std::string &expecting,
                                 const std::string &found) {
    std::stringstream errormsg;
    std::stringstream tmpmsg;
gsell's avatar
gsell committed
480 481 482
    errormsg << "THERE SEEMS TO BE SOMETHING WRONG WITH YOUR FIELD MAP '" << Filename_m << "'.\n"
             << "expecting: '" << expecting << "' on line " << lines_read_m << ",\n"
             << "found instead: '" << found << "'.";
Steve Russell's avatar
Steve Russell committed
483
    std::string errormsg_str = typeset_msg(errormsg.str(), "error");
kraus's avatar
kraus committed
484 485
    ERRORMSG(errormsg_str << "\n" << endl);

gsell's avatar
gsell committed
486
    if(Ippl::myNode() == 0) {
487 488
        std::ofstream omsg("errormsg.txt", std::ios_base::app);
        omsg << errormsg_str << std::endl;
gsell's avatar
gsell committed
489 490 491 492
        omsg.close();
    }
}

493
void Fieldmap::interpreteWarning(const std::ios_base::iostate &state,
gsell's avatar
gsell committed
494
                                 const bool &read_all,
Steve Russell's avatar
Steve Russell committed
495 496 497
                                 const std::string &expecting,
                                 const std::string &found) {
    std::string error_msg;
gsell's avatar
gsell committed
498
    if(!read_all) {
Steve Russell's avatar
Steve Russell committed
499
        error_msg = std::string("Didn't find enough values!");
500
    } else if(state & std::ios_base::eofbit) {
Steve Russell's avatar
Steve Russell committed
501
        error_msg = std::string("Found more values than expected!");
502
    } else if(state & std::ios_base::failbit) {
Steve Russell's avatar
Steve Russell committed
503
        error_msg = std::string("Found wrong type of values!");
gsell's avatar
gsell committed
504 505 506 507 508
    }
    interpreteWarning(error_msg, expecting, found);
}

void Fieldmap::missingValuesWarning() {
509
    std::stringstream errormsg;
gsell's avatar
gsell committed
510 511 512
    errormsg << "THERE SEEMS TO BE SOMETHING WRONG WITH YOUR FIELD MAP '" << Filename_m << "'.\n"
             << "There are only " << lines_read_m - 1 << " lines in the file, expecting more.\n"
             << "Please check the section about field maps in the user manual.";
Steve Russell's avatar
Steve Russell committed
513
    std::string errormsg_str = typeset_msg(errormsg.str(), "error");
gsell's avatar
gsell committed
514

kraus's avatar
kraus committed
515 516
    ERRORMSG(errormsg_str << "\n" << endl);

gsell's avatar
gsell committed
517
    if(Ippl::myNode() == 0) {
518 519
        std::ofstream omsg("errormsg.txt", std::ios_base::app);
        omsg << errormsg_str << std::endl;
gsell's avatar
gsell committed
520 521 522 523 524
        omsg.close();
    }
}

void Fieldmap::exceedingValuesWarning() {
525
    std::stringstream errormsg;
gsell's avatar
gsell committed
526 527 528
    errormsg << "THERE SEEMS TO BE SOMETHING WRONG WITH YOUR FIELD MAP '" << Filename_m << "'.\n"
             << "There are too many lines in the file, expecting only " << lines_read_m << " lines.\n"
             << "Please check the section about field maps in the user manual.";
Steve Russell's avatar
Steve Russell committed
529
    std::string errormsg_str = typeset_msg(errormsg.str(), "error");
kraus's avatar
kraus committed
530 531 532

    ERRORMSG(errormsg_str << "\n" << endl);

gsell's avatar
gsell committed
533
    if(Ippl::myNode() == 0) {
534 535
        std::ofstream omsg("errormsg.txt", std::ios_base::app);
        omsg << errormsg_str << std::endl;
gsell's avatar
gsell committed
536 537 538 539 540
        omsg.close();
    }
}

void Fieldmap::disableFieldmapWarning() {
541
    std::stringstream errormsg;
gsell's avatar
gsell committed
542
    errormsg << "DISABLING FIELD MAP '" + Filename_m + "' DUE TO PARSING ERRORS." ;
Steve Russell's avatar
Steve Russell committed
543
    std::string errormsg_str = typeset_msg(errormsg.str(), "error");
kraus's avatar
kraus committed
544 545 546

    ERRORMSG(errormsg_str << "\n" << endl);

gsell's avatar
gsell committed
547
    if(Ippl::myNode() == 0) {
548 549
        std::ofstream omsg("errormsg.txt", std::ios_base::app);
        omsg << errormsg_str << std::endl;
gsell's avatar
gsell committed
550 551 552 553 554
        omsg.close();
    }
}

void Fieldmap::noFieldmapWarning() {
555
    std::stringstream errormsg;
gsell's avatar
gsell committed
556
    errormsg << "DISABLING FIELD MAP '" << Filename_m << "' SINCE FILE COULDN'T BE FOUND!";
Steve Russell's avatar
Steve Russell committed
557
    std::string errormsg_str = typeset_msg(errormsg.str(), "error");
kraus's avatar
kraus committed
558 559 560

    ERRORMSG(errormsg_str << "\n" << endl);

gsell's avatar
gsell committed
561
    if(Ippl::myNode() == 0) {
562 563
        std::ofstream omsg("errormsg.txt", std::ios_base::app);
        omsg << errormsg.str() << std::endl;
gsell's avatar
gsell committed
564 565 566 567
        omsg.close();
    }
}

568 569 570 571 572 573 574 575 576 577 578
void Fieldmap::lowResolutionWarning(double squareError, double maxError) {
    std::stringstream errormsg;
    errormsg << "IT SEEMS THAT YOU USE TOO FEW FOURIER COMPONENTS TO SUFFICIENTLY WELL\n"
             << "RESOLVE THE FIELD MAP '" << Filename_m << "'.\n"
             << "PLEASE INCREASE THE NUMBER OF FOURIER COMPONENTS!\n"
             << "The ratio (e_i - E_i)^2 / E_i^2 is " << std::to_string(squareError) << " and\n"
             << "the ratio (max_i(|e_i - E_i|) / max_i(|E_i|) is " << std::to_string(maxError) << ".\n"
             << "Here E_i is the field as in the field map and e_i is the reconstructed field.\n"
             << "The lower limit for the two ratios is 1e-2\n"
             << "Have a look into the directory data/ for a reconstruction of the field map.\n";
    std::string errormsg_str = typeset_msg(errormsg.str(), "warning");
kraus's avatar
kraus committed
579 580 581

    ERRORMSG(errormsg_str << "\n" << endl);

582 583 584 585 586 587 588
    if(Ippl::myNode() == 0) {
        std::ofstream omsg("errormsg.txt", std::ios_base::app);
        omsg << errormsg.str() << std::endl;
        omsg.close();
    }
}

Steve Russell's avatar
Steve Russell committed
589 590
std::string Fieldmap::typeset_msg(const std::string &msg, const std::string &title) {
    static std::string frame("* ******************************************************************************\n");
gsell's avatar
gsell committed
591
    static unsigned int frame_width = frame.length() - 5;
Steve Russell's avatar
Steve Russell committed
592
    static std::string closure("                                                                               *\n");
gsell's avatar
gsell committed
593

Steve Russell's avatar
Steve Russell committed
594
    std::string return_string("\n" + frame);
gsell's avatar
gsell committed
595 596 597 598 599 600 601 602 603 604 605 606 607 608 609

    int remaining_length = msg.length();
    unsigned int current_position = 0;

    unsigned int ii = 0;
    for(; ii < title.length(); ++ ii) {
        char c = title[ii];
        c = toupper(c);
        return_string.replace(15 + 2 * ii, 1, " ");
        return_string.replace(16 + 2 * ii, 1, &c, 1);
    }
    return_string.replace(15 + 2 * ii, 1, " ");

    while(remaining_length > 0) {
        size_t eol = msg.find("\n", current_position);
Steve Russell's avatar
Steve Russell committed
610
        std::string next_to_process;
gsell's avatar
gsell committed
611 612
        if(eol + 1 == msg.length()) {
            break;
Steve Russell's avatar
Steve Russell committed
613
        } else if(eol != std::string::npos) {
gsell's avatar
gsell committed
614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648
            next_to_process = msg.substr(current_position, eol - current_position);
        } else {
            next_to_process = msg.substr(current_position);
            eol = msg.length();
        }

        if(eol - current_position < frame_width) {
            return_string += "* " + next_to_process + closure.substr(eol - current_position + 2);
        } else {
            unsigned int last_space = next_to_process.rfind(" ", frame_width);
            if(last_space > 0) {
                if(last_space < frame_width) {
                    return_string += "* " + next_to_process.substr(0, last_space) + closure.substr(last_space + 2);
                } else {
                    return_string += "* " + next_to_process.substr(0, last_space) + " *\n";
                }
                if(next_to_process.length() - last_space + 1 < frame_width) {
                    return_string += "* " + next_to_process.substr(last_space + 1) + closure.substr(next_to_process.length() - last_space + 1);
                } else {
                    return_string += "* " + next_to_process.substr(last_space + 1) + " *\n";
                }
            } else {
                return_string += "* " + next_to_process + " *\n";
            }
        }

        current_position = eol + 1;
        remaining_length = msg.length() - current_position;
    }

    return_string += frame;

    return return_string;
}

649
void Fieldmap::getOnaxisEz(std::vector<std::pair<double, double> > &onaxis)
gsell's avatar
gsell committed
650 651
{ }

652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675
void Fieldmap::Get1DProfile1EngeCoeffs(std::vector<double> &engeCoeffsEntry,
                                       std::vector<double> &engeCoeffsExit) {

}

void Fieldmap::Get1DProfile1EntranceParam(double &entranceParameter1,
        double &entranceParameter2,
        double &entranceParameter3) {

}

void Fieldmap::Get1DProfile1ExitParam(double &exitParameter1,
                                      double &exitParameter2,
                                      double &exitParameter3) {

}

double Fieldmap::GetFieldGap() {
    return 0.0;
}

void Fieldmap::SetFieldGap(double gap) {

}
gsell's avatar
gsell committed
676 677

REGISTER_PARSE_TYPE(int);
678
REGISTER_PARSE_TYPE(unsigned int);
gsell's avatar
gsell committed
679
REGISTER_PARSE_TYPE(double);
Steve Russell's avatar
Steve Russell committed
680
REGISTER_PARSE_TYPE(std::string);
gsell's avatar
gsell committed
681

Steve Russell's avatar
Steve Russell committed
682
std::string Fieldmap::alpha_numeric("abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ0123456789.-+\211");
683
std::map<std::string, Fieldmap::FieldmapDescription> Fieldmap::FieldmapDictionary = std::map<std::string, Fieldmap::FieldmapDescription>();
684
char Fieldmap::buffer_m[READ_BUFFER_LENGTH];