Fieldmap.cpp 29.5 KB
Newer Older
1 2
#include "Fields/Fieldmap.h"
#include "Fields/FM3DDynamic.h"
3
#include "Fields/FM3DH5BlockBase.h"
4 5 6
#include "Fields/FM3DH5Block.h"
#include "Fields/FM3DH5Block_nonscale.h"
#include "Fields/FM3DMagnetoStaticH5Block.h"
7 8
#include "Fields/FM3DMagnetoStatic.h"
#include "Fields/FM3DMagnetoStaticExtended.h"
9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27
#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
28
#include "Utilities/Options.h"
29
#include "Utilities/Util.h"
30
#include "AbstractObjects/OpalData.h"
31 32
#include "Physics/Physics.h"

gsell's avatar
gsell committed
33 34
#include "H5hut.h"

35 36
#include <boost/filesystem.hpp>

37
#include <cmath>
kraus's avatar
kraus committed
38 39 40
#include <iostream>
#include <fstream>
#include <ios>
41
#include <cassert>
42

43 44
namespace fs = boost::filesystem;

gsell's avatar
gsell committed
45 46 47
#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
48
Fieldmap *Fieldmap::getFieldmap(std::string Filename, bool fast) {
49
    std::map<std::string, FieldmapDescription>::iterator position = FieldmapDictionary.find(Filename);
50
    if (position != FieldmapDictionary.end()) {
gsell's avatar
gsell committed
51 52 53 54
        (*position).second.RefCounter++;
        return (*position).second.Map;
    } else {
        MapType type;
55
        std::pair<std::map<std::string, FieldmapDescription>::iterator, bool> position;
gsell's avatar
gsell committed
56 57
        type = readHeader(Filename);
        switch(type) {
58
        case T1DDynamic:
59
            if (fast) {
gsell's avatar
gsell committed
60 61 62 63 64
                position = FieldmapDictionary.insert(
                    std::make_pair(
                        Filename,
                        FieldmapDescription(
                            T1DDynamic, new FM1DDynamic_fast(Filename))));
65
            } else {
gsell's avatar
gsell committed
66 67 68 69 70
                position = FieldmapDictionary.insert(
                    std::make_pair(
                        Filename,
                        FieldmapDescription(
                            T1DDynamic, new FM1DDynamic(Filename))));
71 72 73 74 75
            }
            return (*position.first).second.Map;
            break;

        case TAstraDynamic:
76
            if (fast) {
gsell's avatar
gsell committed
77 78 79 80 81
                position = FieldmapDictionary.insert(
                    std::make_pair(
                        Filename,
                        FieldmapDescription(
                            TAstraDynamic, new Astra1DDynamic_fast(Filename))));
82
            } else {
gsell's avatar
gsell committed
83 84 85 86 87
                position = FieldmapDictionary.insert(
                    std::make_pair(
                        Filename,
                        FieldmapDescription(
                            TAstraDynamic, new Astra1DDynamic(Filename))));
88 89 90 91 92
            }
            return (*position.first).second.Map;
            break;

        case T1DElectroStatic:
93
            if (fast) {
gsell's avatar
gsell committed
94 95 96 97 98
                position = FieldmapDictionary.insert(
                    std::make_pair(
                        Filename,
                        FieldmapDescription(
                            T1DElectroStatic, new FM1DElectroStatic_fast(Filename))));
99
            } else {
gsell's avatar
gsell committed
100 101 102 103 104
                position = FieldmapDictionary.insert(
                    std::make_pair(
                        Filename,
                        FieldmapDescription(
                            T1DElectroStatic, new FM1DElectroStatic(Filename))));
105 106 107 108 109
            }
            return (*position.first).second.Map;
            break;

        case TAstraElectroStatic:
110
            if (fast) {
gsell's avatar
gsell committed
111 112 113 114 115
                position = FieldmapDictionary.insert(
                    std::make_pair(
                        Filename,
                        FieldmapDescription(
                            TAstraElectroStatic, new Astra1DElectroStatic_fast(Filename))));
116
            } else {
gsell's avatar
gsell committed
117 118 119 120 121
                position = FieldmapDictionary.insert(
                    std::make_pair(
                        Filename,
                        FieldmapDescription(
                            TAstraElectroStatic, new Astra1DElectroStatic(Filename))));
122 123 124 125 126
            }
            return (*position.first).second.Map;
            break;

        case T1DMagnetoStatic:
127
            if (fast) {
gsell's avatar
gsell committed
128 129 130 131 132
                position = FieldmapDictionary.insert(
                    std::make_pair(
                        Filename,
                        FieldmapDescription(
                            T1DMagnetoStatic, new FM1DMagnetoStatic_fast(Filename))));
133
            } else {
gsell's avatar
gsell committed
134 135 136 137 138
                position = FieldmapDictionary.insert(
                    std::make_pair(
                        Filename,
                        FieldmapDescription(
                            T1DMagnetoStatic, new FM1DMagnetoStatic(Filename))));
139 140 141 142 143
            }
            return (*position.first).second.Map;
            break;

        case TAstraMagnetoStatic:
144
            if (fast) {
gsell's avatar
gsell committed
145 146 147 148 149
                position = FieldmapDictionary.insert(
                    std::make_pair(
                        Filename,
                        FieldmapDescription(
                            TAstraMagnetoStatic, new Astra1DMagnetoStatic_fast(Filename))));
150
            } else {
gsell's avatar
gsell committed
151 152 153 154 155
                position = FieldmapDictionary.insert(
                    std::make_pair(
                        Filename,
                        FieldmapDescription(
                            TAstraMagnetoStatic, new Astra1DMagnetoStatic(Filename))));
156 157 158 159 160
            }
            return (*position.first).second.Map;
            break;

        case T1DProfile1:
gsell's avatar
gsell committed
161 162 163 164
            position = FieldmapDictionary.insert(
                std::make_pair(
                    Filename,
                    FieldmapDescription(T1DProfile1, new FM1DProfile1(Filename))));
165 166 167 168
            return (*position.first).second.Map;
            break;

        case T1DProfile2:
gsell's avatar
gsell committed
169 170 171 172 173
            position = FieldmapDictionary.insert(
                std::make_pair(
                    Filename,
                    FieldmapDescription(
                        T1DProfile2, new FM1DProfile2(Filename))));
174 175 176 177
            return (*position.first).second.Map;
            break;

        case T2DDynamic:
gsell's avatar
gsell committed
178 179 180 181 182
            position = FieldmapDictionary.insert(
                std::make_pair(
                    Filename,
                    FieldmapDescription(
                        T2DDynamic, new FM2DDynamic(Filename))));
183 184 185 186
            return (*position.first).second.Map;
            break;

        case T2DElectroStatic:
gsell's avatar
gsell committed
187 188 189 190 191
            position = FieldmapDictionary.insert(
                std::make_pair(
                    Filename,
                    FieldmapDescription(
                        T2DElectroStatic, new FM2DElectroStatic(Filename))));
192 193 194 195
            return (*position.first).second.Map;
            break;

        case T2DMagnetoStatic:
gsell's avatar
gsell committed
196 197 198 199 200
            position = FieldmapDictionary.insert(
                std::make_pair(
                    Filename,
                    FieldmapDescription(
                        T2DMagnetoStatic, new FM2DMagnetoStatic(Filename))));
201 202 203 204
            return (*position.first).second.Map;
            break;

        case T3DDynamic:
gsell's avatar
gsell committed
205 206 207 208 209
            position = FieldmapDictionary.insert(
                std::make_pair(
                    Filename,
                    FieldmapDescription(
                        T3DDynamic, new FM3DDynamic(Filename))));
210 211 212 213
            return (*position.first).second.Map;
            break;

        case T3DMagnetoStaticH5Block:
gsell's avatar
gsell committed
214 215 216 217
            position = FieldmapDictionary.insert(
                std::make_pair(
                    Filename, FieldmapDescription(
                        T3DMagnetoStaticH5Block, new FM3DMagnetoStaticH5Block(Filename))));
218 219 220
            return (*position.first).second.Map;
            break;

221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236
        case T3DMagnetoStatic:
            position = FieldmapDictionary.insert(
                std::make_pair(
                    Filename, FieldmapDescription(
                        T3DMagnetoStatic, new FM3DMagnetoStatic(Filename))));
            return (*position.first).second.Map;
            break;

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

237
        case T3DDynamicH5Block:
238
            if (fast) {
gsell's avatar
gsell committed
239 240 241 242 243
                position = FieldmapDictionary.insert(
                    std::make_pair(
                        Filename,
                        FieldmapDescription(
                            T3DDynamic, new FM3DH5Block_nonscale(Filename))));
244
            } else {
gsell's avatar
gsell committed
245 246 247 248 249
                position = FieldmapDictionary.insert(
                    std::make_pair(
                        Filename,
                        FieldmapDescription(
                            T3DDynamic, new FM3DH5Block(Filename))));
250 251 252 253 254
            }
            return (*position.first).second.Map;
            break;

        default:
255 256
            throw GeneralClassicException("Fieldmap::getFieldmap()",
                                          "Couldn't determine type of fieldmap in file \"" + Filename + "\"");
gsell's avatar
gsell committed
257 258 259 260
        }
    }
}

261 262
std::vector<std::string> Fieldmap::getListFieldmapNames() {
    std::vector<std::string> name_list;
263
    for (std::map<std::string, FieldmapDescription>::const_iterator it = FieldmapDictionary.begin(); it != FieldmapDictionary.end(); ++ it) {
gsell's avatar
gsell committed
264 265 266 267 268
        name_list.push_back((*it).first);
    }
    return name_list;
}

Steve Russell's avatar
Steve Russell committed
269
void Fieldmap::deleteFieldmap(std::string Filename) {
270
    freeMap(Filename);
gsell's avatar
gsell committed
271 272
}

273 274 275 276 277 278 279 280 281
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
282
MapType Fieldmap::readHeader(std::string Filename) {
gsell's avatar
gsell committed
283
    char magicnumber[5] = "    ";
Steve Russell's avatar
Steve Russell committed
284
    std::string buffer;
gsell's avatar
gsell committed
285 286 287
    int lines_read_m = 0;

    // Check for default map(s).
288
    if (Filename == "1DPROFILE1-DEFAULT")
gsell's avatar
gsell committed
289 290
        return T1DProfile1;

291 292 293 294
    if (!fs::exists(Filename))
        throw GeneralClassicException("Fieldmap::readHeader()",
                                      "File \"" + Filename + "\" doesn't exist");

295
    std::ifstream File(Filename.c_str());
296
    if (!File.good()) {
297
        std::cerr << "could not open file " << Filename << std::endl;
gsell's avatar
gsell committed
298 299 300 301
        return UNKNOWN;
    }

    getLine(File, lines_read_m, buffer);
302
    std::istringstream interpreter(buffer, std::istringstream::in);
gsell's avatar
gsell committed
303 304 305

    interpreter.read(magicnumber, 4);

306
    if (std::strcmp(magicnumber, "3DDy") == 0)
gsell's avatar
gsell committed
307 308
        return T3DDynamic;

309
    if (std::strcmp(magicnumber, "3DMa") == 0) {
310 311 312
        char tmpString[21] = "                    ";
        interpreter.read(tmpString, 20);

313
        if (std::strcmp(tmpString, "gnetoStatic_Extended") == 0)
314 315 316 317
            return T3DMagnetoStatic_Extended;
        else
            return T3DMagnetoStatic;
    }
gsell's avatar
gsell committed
318

319
    if (std::strcmp(magicnumber, "3DEl") == 0)
gsell's avatar
gsell committed
320 321
        return T3DElectroStatic;

322
    if (std::strcmp(magicnumber, "2DDy") == 0) {
323 324
        // char tmpString[14] = "             ";
        // interpreter.read(tmpString, 13);
325
        return T2DDynamic;
gsell's avatar
gsell committed
326 327
    }

328
    if (std::strcmp(magicnumber, "2DMa") == 0) {
329 330
        // char tmpString[20] = "                   ";
        // interpreter.read(tmpString, 19);
331
        return T2DMagnetoStatic;
gsell's avatar
gsell committed
332 333
    }

334
    if (std::strcmp(magicnumber, "2DEl") == 0) {
335 336
        // char tmpString[20] = "                   ";
        // interpreter.read(tmpString, 19);
337
        return T2DElectroStatic;
gsell's avatar
gsell committed
338 339
    }

340
    if (std::strcmp(magicnumber, "1DDy") == 0)
gsell's avatar
gsell committed
341 342
        return T1DDynamic;

343
    if (std::strcmp(magicnumber, "1DMa") == 0)
gsell's avatar
gsell committed
344 345
        return T1DMagnetoStatic;

346
    if (std::strcmp(magicnumber, "1DPr") == 0) {
347 348
        // char tmpString[7] = "      ";
        // interpreter.read(tmpString, 6);
349
        // if (strcmp(tmpString, "ofile1") == 0)
350
        return T1DProfile1;
351
        // if (strcmp(tmpString, "ofile2") == 0)
352
        //     return T1DProfile2;
gsell's avatar
gsell committed
353 354
    }

355
    if (std::strcmp(magicnumber, "1DEl") == 0)
gsell's avatar
gsell committed
356 357
        return T1DElectroStatic;

358
    if (std::strcmp(magicnumber, "\211HDF") == 0) {
gsell's avatar
gsell committed
359 360 361 362 363
        h5_err_t h5err = 0;
#if defined (NDEBUG)
        // mark variable as unused
        (void)h5err;
#endif
gsell's avatar
gsell committed
364
        char name[20];
365
        h5_size_t len_name = sizeof(name);
kraus's avatar
kraus committed
366
        h5_prop_t props = H5CreateFileProp ();
367 368 369 370
        MPI_Comm comm = Ippl::getComm();
        h5err = H5SetPropFileMPIOCollective (props, &comm);
        assert (h5err != H5_ERR);
        h5_file_t file = H5OpenFile (Filename.c_str(), H5_O_RDONLY, props);
kraus's avatar
kraus committed
371 372
        assert (file != (h5_file_t)H5_ERR);
        H5CloseProp (props);
373

kraus's avatar
kraus committed
374
        h5err = H5SetStep(file, 0);
375
        assert (h5err != H5_ERR);
gsell's avatar
gsell committed
376

kraus's avatar
kraus committed
377
        h5_int64_t num_fields = H5BlockGetNumFields(file);
378
        assert (num_fields != H5_ERR);
kraus's avatar
kraus committed
379
        MapType maptype = UNKNOWN;
380

381
        for (h5_ssize_t i = 0; i < num_fields; ++ i) {
gsell's avatar
gsell committed
382 383
            h5err = H5BlockGetFieldInfo(
                file, (h5_size_t)i, name, len_name, NULL, NULL, NULL, NULL);
384 385
            assert (h5err != H5_ERR);
            // using field name "Bfield" and "Hfield" to distinguish the type
386
            if (std::strcmp(name, "Bfield") == 0) {
387 388
                maptype = T3DMagnetoStaticH5Block;
                break;
389
            } else if (std::strcmp(name, "Hfield") == 0) {
390 391
                maptype = T3DDynamicH5Block;
                break;
gsell's avatar
gsell committed
392 393
            }
        }
394
        h5err = H5CloseFile(file);
gsell's avatar
gsell committed
395
        assert (h5err != H5_ERR);
396
        if (maptype != UNKNOWN)
397
            return maptype;
gsell's avatar
gsell committed
398
    }
399

400
    if (std::strcmp(magicnumber, "Astr") == 0) {
gsell's avatar
gsell committed
401 402
        char tmpString[3] = "  ";
        interpreter.read(tmpString, 2);
403
        if (std::strcmp(tmpString, "aE") == 0) {
gsell's avatar
gsell committed
404 405
            return TAstraElectroStatic;
        }
406
        if (std::strcmp(tmpString, "aM") == 0) {
gsell's avatar
gsell committed
407 408
            return TAstraMagnetoStatic;
        }
409
        if (std::strcmp(tmpString, "aD") == 0) {
gsell's avatar
gsell committed
410 411 412 413 414 415 416 417
            return TAstraDynamic;
        }
    }


    return UNKNOWN;
}

Steve Russell's avatar
Steve Russell committed
418
void Fieldmap::readMap(std::string Filename) {
419
    std::map<std::string, FieldmapDescription>::iterator position = FieldmapDictionary.find(Filename);
420 421
    if (position != FieldmapDictionary.end())
        if (!(*position).second.read) {
gsell's avatar
gsell committed
422
            (*position).second.Map->readMap();
423 424
            (*position).second.read = true;
        }
gsell's avatar
gsell committed
425 426
}

Steve Russell's avatar
Steve Russell committed
427
void Fieldmap::freeMap(std::string Filename) {
428
    std::map<std::string, FieldmapDescription>::iterator position = FieldmapDictionary.find(Filename);
429 430 431
    /*
      FIXME: find( ) make problem, crashes
    */
432 433
    if (position != FieldmapDictionary.end()) {
        if ((*position).second.RefCounter > 0) {
434 435 436 437 438 439
            (*position).second.RefCounter--;
        }

        if ((*position).second.RefCounter == 0) {
            delete (*position).second.Map;
            (*position).second.Map = NULL;
gsell's avatar
gsell committed
440 441 442 443 444
            FieldmapDictionary.erase(position);
        }
    }
}

445 446 447 448 449 450 451
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;
452
    unsigned int sizeSampling = std::round(length / deltaZ);
453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471
    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(".");
472 473 474
    size_t lastSlash = Filename_m.find_last_of("/");
    lastSlash = (lastSlash == std::string::npos)? 0: lastSlash + 1;

475
    auto opal = OpalData::getInstance();
476
    std::ofstream out;
477
    if (Ippl::myNode() == 0 && !opal->isOptimizerRun()) {
478
        std::string fname = Util::combineFilePath({
frey_m's avatar
frey_m committed
479
            opal->getAuxiliaryOutputDirectory(),
480 481 482
            Filename_m.substr(lastSlash, lastDot) + ".check"
        });
        out.open(fname);
483 484 485 486 487 488 489 490
        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;
491 492 493
        for (unsigned int l = 1; l < accuracy; ++l, n += 2) {
            double coskzl = std::cos(kz * l);
            double sinkzl = std::sin(kz * l);
494 495 496 497 498 499 500 501 502

            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);
503 504 505 506 507 508 509

        if (Ippl::myNode() == 0 && !opal->isOptimizerRun()) {
            out << std::setw(16) << std::setprecision(8) << *it
                << std::setw(16) << std::setprecision(8) << ez
                << std::setw(16) << std::setprecision(8) << onAxisFieldCheck
                << std::endl;
        }
510 511 512
    }
    out.close();

513 514
    if (std::sqrt(error / ezSquare) > 1e-1 || maxDiff > 1e-1 * ezMax) {
        lowResolutionWarning(std::sqrt(error / ezSquare), maxDiff / ezMax);
515

516
        throw GeneralClassicException("Fieldmap::checkMap",
517 518
                                      "Field map can't be reproduced properly with the given number of fourier components");
    }
519 520
    if (std::sqrt(error / ezSquare) > 1e-2 || maxDiff > 1e-2 * ezMax) {
        lowResolutionWarning(std::sqrt(error / ezSquare), maxDiff / ezMax);
521 522 523
    }
}

524
void Fieldmap::setEdgeConstants(const double &/*bendAngle*/, const double &/*entranceAngle*/, const double &/*exitAngle*/)
gsell's avatar
gsell committed
525 526 527 528 529
{};

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

530
void Fieldmap::getLine(std::ifstream &in, int &lines_read, std::string &buffer) {
gsell's avatar
gsell committed
531 532 533 534 535 536 537
    size_t firstof = 0;
    size_t lastof;

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

Steve Russell's avatar
Steve Russell committed
538
        buffer = std::string(buffer_m);
gsell's avatar
gsell committed
539

540
        size_t comment = buffer.find("#");
gsell's avatar
gsell committed
541 542 543 544
        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
545
    } while(!in.eof() && lastof == std::string::npos);
gsell's avatar
gsell committed
546

547
    if (firstof != std::string::npos) {
gsell's avatar
gsell committed
548 549 550 551
        buffer = buffer.substr(firstof, lastof - firstof + 1);
    }
}

552
bool Fieldmap::interpreteEOF(std::ifstream &in) {
gsell's avatar
gsell committed
553 554 555
    while(!in.eof()) {
        ++lines_read_m;
        in.getline(buffer_m, READ_BUFFER_LENGTH);
Steve Russell's avatar
Steve Russell committed
556
        std::string buffer(buffer_m);
gsell's avatar
gsell committed
557 558 559
        size_t comment = buffer.find_first_of("#");
        buffer = buffer.substr(0, comment);
        size_t lasto = buffer.find_first_of(alpha_numeric);
560
        if (lasto != std::string::npos) {
gsell's avatar
gsell committed
561 562 563 564 565 566 567
            exceedingValuesWarning();
            return false;
        }
    }
    return true;
}

568
void Fieldmap::interpretWarning(const std::ios_base::iostate &state,
gsell's avatar
gsell committed
569
                                 const bool &read_all,
Steve Russell's avatar
Steve Russell committed
570 571
                                 const std::string &expecting,
                                 const std::string &found) {
572 573
    std::stringstream errormsg;
    errormsg << "THERE SEEMS TO BE SOMETHING WRONG WITH YOUR FIELD MAP '" << Filename_m << "'.\n";
574
    if (!read_all) {
575
        errormsg << "Didn't find enough values!" << std::endl;
576
    } else if (state & std::ios_base::eofbit) {
577
        errormsg << "Found more values than expected!" << std::endl;
578
    } else if (state & std::ios_base::failbit) {
579 580 581
        errormsg << "Found wrong type of values!" << "\n"
                 << "expecting: '" << expecting << "' on line " << lines_read_m << ",\n"
                 << "instead found: '" << found << "'." << std::endl;
gsell's avatar
gsell committed
582
    }
583 584
    throw GeneralClassicException("Fieldmap::interpretWarning()",
                                  errormsg.str());
gsell's avatar
gsell committed
585 586 587
}

void Fieldmap::missingValuesWarning() {
588
    std::stringstream errormsg;
gsell's avatar
gsell committed
589 590 591
    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.";
592

593 594
    throw GeneralClassicException("Fieldmap::missingValuesWarning()",
                                  errormsg.str());
gsell's avatar
gsell committed
595 596 597
}

void Fieldmap::exceedingValuesWarning() {
598
    std::stringstream errormsg;
gsell's avatar
gsell committed
599 600 601
    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.";
602

603 604
    throw GeneralClassicException("Fieldmap::exceedingValuesWarning()",
                                  errormsg.str());
gsell's avatar
gsell committed
605 606 607
}

void Fieldmap::disableFieldmapWarning() {
608
    std::stringstream errormsg;
gsell's avatar
gsell committed
609
    errormsg << "DISABLING FIELD MAP '" + Filename_m + "' DUE TO PARSING ERRORS." ;
610

611 612
    throw GeneralClassicException("Fieldmap::disableFieldmapsWarning()",
                                  errormsg.str());
gsell's avatar
gsell committed
613 614 615
}

void Fieldmap::noFieldmapWarning() {
616
    std::stringstream errormsg;
gsell's avatar
gsell committed
617
    errormsg << "DISABLING FIELD MAP '" << Filename_m << "' SINCE FILE COULDN'T BE FOUND!";
618

619 620
    throw GeneralClassicException("Fieldmap::noFieldmapsWarning()",
                                  errormsg.str());
gsell's avatar
gsell committed
621 622
}

623 624 625 626 627 628 629 630 631
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"
632 633 634
             << "Have a look into the directory "
             << OpalData::getInstance()->getAuxiliaryOutputDirectory()
             << " for a reconstruction of the field map.\n";
635
    std::string errormsg_str = typeset_msg(errormsg.str(), "warning");
kraus's avatar
kraus committed
636 637 638

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

639
    if (Ippl::myNode() == 0) {
640 641 642 643 644 645
        std::ofstream omsg("errormsg.txt", std::ios_base::app);
        omsg << errormsg.str() << std::endl;
        omsg.close();
    }
}

Steve Russell's avatar
Steve Russell committed
646 647
std::string Fieldmap::typeset_msg(const std::string &msg, const std::string &title) {
    static std::string frame("* ******************************************************************************\n");
gsell's avatar
gsell committed
648
    static unsigned int frame_width = frame.length() - 5;
Steve Russell's avatar
Steve Russell committed
649
    static std::string closure("                                                                               *\n");
gsell's avatar
gsell committed
650

Steve Russell's avatar
Steve Russell committed
651
    std::string return_string("\n" + frame);
gsell's avatar
gsell committed
652 653 654 655 656

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

    unsigned int ii = 0;
657
    for (; ii < title.length(); ++ ii) {
gsell's avatar
gsell committed
658
        char c = title[ii];
659
        c = std::toupper(c);
gsell's avatar
gsell committed
660 661 662 663 664 665 666
        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
667
        std::string next_to_process;
668
        if (eol != std::string::npos) {
gsell's avatar
gsell committed
669 670 671 672 673 674
            next_to_process = msg.substr(current_position, eol - current_position);
        } else {
            next_to_process = msg.substr(current_position);
            eol = msg.length();
        }

675
        if (eol - current_position < frame_width) {
gsell's avatar
gsell committed
676 677 678
            return_string += "* " + next_to_process + closure.substr(eol - current_position + 2);
        } else {
            unsigned int last_space = next_to_process.rfind(" ", frame_width);
679 680
            if (last_space > 0) {
                if (last_space < frame_width) {
gsell's avatar
gsell committed
681 682 683 684
                    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";
                }
685
                if (next_to_process.length() - last_space + 1 < frame_width) {
gsell's avatar
gsell committed
686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702
                    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;
}

703
void Fieldmap::getOnaxisEz(std::vector<std::pair<double, double> > &/*onaxis*/)
gsell's avatar
gsell committed
704 705
{ }

706 707
void Fieldmap::get1DProfile1EngeCoeffs(std::vector<double> &/*engeCoeffsEntry*/,
                                       std::vector<double> &/*engeCoeffsExit*/) {
708 709 710

}

711 712 713
void Fieldmap::get1DProfile1EntranceParam(double &/*entranceParameter1*/,
                                          double &/*entranceParameter2*/,
                                          double &/*entranceParameter3*/) {
714 715 716

}

717 718 719
void Fieldmap::get1DProfile1ExitParam(double &/*exitParameter1*/,
                                      double &/*exitParameter2*/,
                                      double &/*exitParameter3*/) {
720 721 722

}

723
double Fieldmap::getFieldGap() {
724 725 726
    return 0.0;
}

727
void Fieldmap::setFieldGap(double /*gap*/) {
728 729

}
gsell's avatar
gsell committed
730

731 732 733 734 735 736 737 738 739 740 741 742
void Fieldmap::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) {
    const size_t numpoints = nx * ny * nz;
    if (Ippl::myNode() != 0 ||
        (ef.size() != numpoints && bf.size() != numpoints)) return;

743
    boost::filesystem::path p(Filename_m);
744 745
    std::string fname = Util::combineFilePath({
        OpalData::getInstance()->getAuxiliaryOutputDirectory(),
746
        p.stem().string() + ".vtk"
747
    });
748
    std::ofstream of;
749
    of.open (fname);
750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804
    assert (of.is_open ());
    of.precision (6);

    const double hx = (xrange.second - xrange.first) / (nx - 1);
    const double hy = (yrange.second - yrange.first) / (ny - 1);
    const double hz = (zrange.second - zrange.first) / (nz - 1);

    of << "# vtk DataFile Version 2.0" << std::endl;
    of << "generated by 3D fieldmaps" << std::endl;
    of << "ASCII" << std::endl << std::endl;
    of << "DATASET RECTILINEAR_GRID" << std::endl;
    of << "DIMENSIONS " << nx << " " << ny << " " << nz << std::endl;

    of << "X_COORDINATES " << nx << " float" << std::endl;
    of << xrange.first;
    for (unsigned int i = 1; i < nx - 1; ++ i) {
        of << " " << xrange.first + i * hx;
    }
    of << " " << xrange.second << std::endl;

    of << "Y_COORDINATES " << ny << " float" << std::endl;
    of << yrange.first;
    for (unsigned int i = 1; i < ny - 1; ++ i) {
        of << " " << yrange.first + i * hy;
    }
    of << " " << yrange.second << std::endl;

    of << "Z_COORDINATES " << nz << " float" << std::endl;
    of << zrange.first;
    for (unsigned int i = 1; i < nz - 1; ++ i) {
        of << " " << zrange.first + i * hz;
    }
    of << " " << zrange.second << std::endl;

    of << "POINT_DATA " << numpoints << std::endl;

    if (ef.size() == numpoints) {
        of << "VECTORS EField float" << std::endl;
        // of << "LOOKUP_TABLE default" << std::endl;
        for (size_t i = 0; i < numpoints; ++ i) {
            of << ef[i](0) << " " << ef[i](1) << " " << ef[i](2) << std::endl;
        }
        // of << std::endl;
    }

    if (bf.size() == numpoints) {
        of << "VECTORS BField float" << std::endl;
        // of << "LOOKUP_TABLE default" << std::endl;
        for (size_t i = 0; i < numpoints; ++ i) {
            of << bf[i](0) << " " << bf[i](1) << " " << bf[i](2) << std::endl;
        }
        // of << std::endl;
    }
}

gsell's avatar
gsell committed
805
REGISTER_PARSE_TYPE(int);
806
REGISTER_PARSE_TYPE(unsigned int);
gsell's avatar
gsell committed
807
REGISTER_PARSE_TYPE(double);
Steve Russell's avatar
Steve Russell committed
808
REGISTER_PARSE_TYPE(std::string);
gsell's avatar
gsell committed
809

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