Cyclotron.cpp 50.2 KB
Newer Older
gsell's avatar
gsell committed
1 2
//
// Definitions for class: Cyclotron
gsell's avatar
gsell committed
3
// Defines the abstract interface for a cyclotron.
gsell's avatar
gsell committed
4 5
// Class category: AbsBeamline
//
gsell's avatar
gsell committed
6 7 8 9 10
// Copyright (c) 2008-2019
// Paul Scherrer Institut, Villigen PSI, Switzerland
// All rights reserved.
//
// OPAL is licensed under GNU GPL version 3.
gsell's avatar
gsell committed
11 12 13
//

#include "AbsBeamline/Cyclotron.h"
14

gsell's avatar
gsell committed
15
#include "AbsBeamline/BeamlineVisitor.h"
16 17
#include "Algorithms/PartBunchBase.h"
#include "Fields/Fieldmap.h"
gsell's avatar
gsell committed
18
#include "Physics/Physics.h"
19
#include "Structure/LossDataSink.h"
20
#include "TrimCoils/TrimCoil.h"
kraus's avatar
kraus committed
21 22
#include "Utilities/Options.h"
#include "Utilities/GeneralClassicException.h"
gsell's avatar
gsell committed
23 24
#include <fstream>

25
#define CHECK_CYC_FSCANF_EOF(arg) if(arg == EOF)\
kraus's avatar
kraus committed
26
throw GeneralClassicException("Cyclotron::getFieldFromFile",\
27
                              "fscanf returned EOF at " #arg);
28

gsell's avatar
gsell committed
29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46
extern Inform *gmsg;

using Physics::pi;

// Class Cyclotron
// ------------------------------------------------------------------------

Cyclotron::Cyclotron():
    Component() {
}


Cyclotron::Cyclotron(const Cyclotron &right):
    Component(right),
    fmapfn_m(right.fmapfn_m),
    rffrequ_m(right.rffrequ_m),
    rfphi_m(right.rfphi_m),
    escale_m(right.escale_m),
47
    superpose_m(right.superpose_m),
gsell's avatar
gsell committed
48 49 50 51
    symmetry_m(right.symmetry_m),
    rinit_m(right.rinit_m),
    prinit_m(right.prinit_m),
    phiinit_m(right.phiinit_m),
52 53
    zinit_m(right.zinit_m),
    pzinit_m(right.pzinit_m),
54
    spiral_flag_m(right.spiral_flag_m),
55
    trimCoilThreshold_m(right.trimCoilThreshold_m),
gsell's avatar
gsell committed
56 57 58
    type_m(right.type_m),
    harm_m(right.harm_m),
    bscale_m(right.bscale_m),
59
    trimcoils_m(right.trimcoils_m),
60 61 62 63
    minr_m(right.minr_m),
    maxr_m(right.maxr_m),
    minz_m(right.minz_m),
    maxz_m(right.maxz_m),
64 65
    fmLowE_m(right.fmLowE_m),
    fmHighE_m(right.fmHighE_m),
66
    RFfilename_m(right.RFfilename_m),
67
    RFFCoeff_fn_m(right.RFFCoeff_fn_m),
68
    RFVCoeff_fn_m(right.RFVCoeff_fn_m) {
gsell's avatar
gsell committed
69 70 71
}


72
Cyclotron::Cyclotron(const std::string &name):
gsell's avatar
gsell committed
73 74 75 76 77 78 79 80
    Component(name) {
}


Cyclotron::~Cyclotron() {
}


gsell's avatar
gsell committed
81 82 83
void Cyclotron::applyTrimCoil_m(const double r, const double z,
                                const double tet_rad,
                                double *br, double *bz) {
84
     for (auto trimcoil : trimcoils_m) {
85
         trimcoil->applyField(r,tet_rad,z,br,bz);
86
     }
87 88
}

gsell's avatar
gsell committed
89 90 91
void Cyclotron::applyTrimCoil(const double r, const double z,
                              const double tet_rad,
                              double& br, double& bz) {
frey_m's avatar
frey_m committed
92
    //Changed from > to >= to include case where bz == 0 and trimCoilThreshold_m == 0 -DW
93
    if (std::abs(bz) >= trimCoilThreshold_m)  {
94
        applyTrimCoil_m(r, z, tet_rad, &br, &bz);
95
    }
frey_m's avatar
frey_m committed
96 97 98
    else {
        // make sure to have a smooth transition
        double tmp_bz = 0.0;
99
        applyTrimCoil_m(r, z, tet_rad, &br, &tmp_bz);
frey_m's avatar
frey_m committed
100 101 102 103
        bz += tmp_bz * std::abs(bz) / trimCoilThreshold_m;
    }
}

gsell's avatar
gsell committed
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
void Cyclotron::accept(BeamlineVisitor &visitor) const {
    visitor.visitCyclotron(*this);
}

void Cyclotron::setRinit(double rinit) {
    rinit_m = rinit;
}

double Cyclotron::getRinit() const {
    return rinit_m;
}

void Cyclotron::setPRinit(double prinit) {
    prinit_m = prinit;
}

double Cyclotron::getPRinit() const {
    return prinit_m;
}

void Cyclotron::setPHIinit(double phiinit) {
    phiinit_m = phiinit;
}

double Cyclotron::getPHIinit() const {
    return phiinit_m;
}

132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147
void Cyclotron::setZinit(double zinit){
    zinit_m = zinit;
}

double Cyclotron::getZinit() const {
    return zinit_m;
}

void Cyclotron::setPZinit(double pzinit){
    pzinit_m = pzinit;
}

double Cyclotron::getPZinit() const {
    return pzinit_m;
}

148
void Cyclotron::setTrimCoilThreshold(double trimCoilThreshold) {
149
    trimCoilThreshold_m = 10.0 * trimCoilThreshold; // convert from Tesla to kGauss
150 151 152 153 154 155
}

double Cyclotron::getTrimCoilThreshold() const {
    return trimCoilThreshold_m;
}

156 157 158 159 160 161 162 163
void Cyclotron::setSpiralFlag(bool spiral_flag) {
    spiral_flag_m = spiral_flag;
}

bool Cyclotron::getSpiralFlag() const {
    return spiral_flag_m;
}

164
void Cyclotron::setFieldMapFN(std::string f) {
gsell's avatar
gsell committed
165 166 167
    fmapfn_m = f;
}

gsell's avatar
gsell committed
168
std::string Cyclotron::getFieldMapFN() const {
gsell's avatar
gsell committed
169 170 171
    return fmapfn_m;
}

gsell's avatar
gsell committed
172
void Cyclotron::setRfFieldMapFN(std::vector<std::string> f) {
gsell's avatar
gsell committed
173 174 175
    RFfilename_m = f;
}

gsell's avatar
gsell committed
176
void Cyclotron::setRFFCoeffFN(std::vector<std::string> f) {
177 178 179
    RFFCoeff_fn_m = f;
}

gsell's avatar
gsell committed
180
void Cyclotron::setRFVCoeffFN(std::vector<std::string> f) {
181 182 183
    RFVCoeff_fn_m = f;
}

184 185 186 187
int Cyclotron::getFieldFlag(const std::string& type) const {
    /**
     * To ease the initialise() function, set a integral parameter fieldflag internally.
     * Its value is  by the option "TYPE" of the element  "CYCLOTRON"
gsell's avatar
gsell committed
188 189 190 191 192 193
     * fieldflag = 1, read in PSI format measured field file (default)
     * fieldflag = 2, read in carbon cyclotron field file created by Jianjun Yang, TYPE=CARBONCYCL
     * fieldflag = 3, read in ANSYS format file for CYCIAE-100 created by Jianjun Yang, TYPE=CYCIAE
     * fieldflag = 4, read in AVFEQ format file for Riken cyclotrons
     * fieldflag = 5, read in FFA format file for MSU/FNAL FFA
     * fieldflag = 6, read in both median plane B field map and 3D E field map
gsell's avatar
gsell committed
194
     *                of RF cavity for compact cyclotron
195 196 197 198 199 200 201 202 203
     * fieldflag = 7, read in fields for Daniel's synchrocyclotron simulations
     */
    int  fieldflag;
    if(type == std::string("CARBONCYCL")) {
        fieldflag = 2;
    } else if(type == std::string("CYCIAE")) {
        fieldflag = 3;
    } else if(type == std::string("AVFEQ")) {
        fieldflag = 4;
ext-rogers_c's avatar
ext-rogers_c committed
204
    } else if(type == std::string("FFA")) {
205 206 207 208 209 210 211 212 213 214
        fieldflag = 5;
    } else if(type == std::string("BANDRF")) {
        fieldflag = 6;
    } else if(type == std::string("SYNCHROCYCLOTRON")) {
        fieldflag = 7;
    } else //(type == "RING")
        fieldflag = 1;
    return fieldflag;
}

gsell's avatar
gsell committed
215
void Cyclotron::setRfPhi(std::vector<double> f) {
gsell's avatar
gsell committed
216 217 218
    rfphi_m = f;
}

219 220 221 222 223 224 225 226 227
double Cyclotron::getRfPhi(unsigned int i) const {
    if (i < rfphi_m.size())
        return rfphi_m[i];
    else {
        throw GeneralClassicException("Cyclotron::getRfPhi",
                                      "RFPHI not defined for CYCLOTRON!");
    }
}

gsell's avatar
gsell committed
228
void Cyclotron::setRfFrequ(std::vector<double> f) {
gsell's avatar
gsell committed
229 230 231
    rffrequ_m = f;
}

232 233 234
double Cyclotron::getRfFrequ(unsigned int i) const {
    if (i < rffrequ_m.size())
        return rffrequ_m[i];
235 236 237 238
    else {
        throw GeneralClassicException("Cyclotron::getRfFrequ",
                                      "RFFREQ not defined for CYCLOTRON!");
    }
adelmann's avatar
adelmann committed
239 240
}

241
void Cyclotron::setSuperpose(std::vector<bool> flag) {
adelmann's avatar
adelmann committed
242
  superpose_m = flag;
243 244
}

245 246 247 248 249 250 251 252
bool Cyclotron::getSuperpose(unsigned int i) const {
    if (i < superpose_m.size())
        return superpose_m[i];
    else {
        throw GeneralClassicException("Cyclotron::getSuperpose",
                                      "SUPERPOSE not defined for CYCLOTRON!");
    }
}
253

gsell's avatar
gsell committed
254 255 256 257 258 259 260 261 262
void Cyclotron::setSymmetry(double s) {
    symmetry_m = s;
}

double Cyclotron::getSymmetry() const {
    return symmetry_m;
}


263
void Cyclotron::setType(std::string t) {
gsell's avatar
gsell committed
264 265 266
    type_m = t;
}

267
const std::string &Cyclotron::getCyclotronType() const {
gsell's avatar
gsell committed
268 269 270
    return type_m;
}

271 272 273 274
ElementBase::ElementType Cyclotron::getType() const {
    return CYCLOTRON;
}

gsell's avatar
gsell committed
275 276 277 278 279 280 281 282 283 284 285 286
void Cyclotron::setCyclHarm(double h) {
    harm_m = h;
}

void Cyclotron::setBScale(double s) {
    bscale_m = s;
}

double Cyclotron::getBScale() const {
    return bscale_m;
}

gsell's avatar
gsell committed
287
void Cyclotron::setEScale(std::vector<double> s) {
gsell's avatar
gsell committed
288 289 290
    escale_m = s;
}

291 292 293 294 295 296 297 298 299
double Cyclotron::getEScale(unsigned int i) const {
    if (i < escale_m.size())
        return escale_m[i];
    else {
        throw GeneralClassicException("Cyclotron::EScale",
                                      "EScale not defined for CYCLOTRON!");
    }
}

300
unsigned int Cyclotron::getNumberOfTrimcoils() const {
301
  return trimcoils_m.size();
302 303
}

gsell's avatar
gsell committed
304 305 306 307 308 309 310 311 312 313 314 315 316 317
double Cyclotron::getCyclHarm() const {
    return harm_m;
}

double Cyclotron::getRmin() const {
    return BP.rmin;
}


double Cyclotron::getRmax() const {
    return BP.rmin + (Bfield.nrad - 1) * BP.delr;
}


318
void Cyclotron::setMinR(double r) {
319 320 321
    // DW: This is to let the user keep using mm in the input file for now
    // while switching internally to m
    minr_m = 0.001 * r;
322 323 324
}

void Cyclotron::setMaxR(double r) {
325 326 327
    // DW: This is to let the user keep using mm in the input file for now
    // while switching internally to m
    maxr_m = 0.001 * r;
328
}
329 330 331 332
double Cyclotron::getMinR() const {
    return minr_m;
}

333 334 335 336 337
double Cyclotron::getMaxR() const {
    return maxr_m;
}

void  Cyclotron::setMinZ(double z) {
338 339 340
    // DW: This is to let the user keep using mm in the input file for now
    // while switching internally to m
    minz_m = 0.001 * z;
341 342 343 344 345
}
double Cyclotron::getMinZ() const {
    return minz_m;
}
void Cyclotron::setMaxZ(double z) {
346 347 348
    // DW: This is to let the user keep using mm in the input file for now
    // while switching internally to m
    maxz_m = 0.001 * z;
349 350 351 352
}
double Cyclotron::getMaxZ() const {
    return maxz_m;
}
gsell's avatar
gsell committed
353

frey_m's avatar
frey_m committed
354
void Cyclotron::setTrimCoils(const std::vector<TrimCoil*> &trimcoils) {
355
    trimcoils_m = trimcoils;
frey_m's avatar
frey_m committed
356 357
}

358 359 360 361 362 363 364
void Cyclotron::setFMLowE(double e) { fmLowE_m = e;}
double Cyclotron::getFMLowE() const { return fmLowE_m;}

void Cyclotron::setFMHighE(double e) { fmHighE_m = e;}
double Cyclotron::getFMHighE() const { return fmHighE_m;}


365
bool Cyclotron::apply(const size_t &id, const double &t, Vector_t &E, Vector_t &B) {
kraus's avatar
kraus committed
366

367 368 369 370 371 372 373 374 375 376 377 378
    bool flagNeedUpdate = false;

    const double rpos = std::hypot(RefPartBunch_m->R[id](0), RefPartBunch_m->R[id](1));
    const double zpos = RefPartBunch_m->R[id](2);

    Inform gmsgALL("OPAL", INFORM_ALL_NODES);
    if (zpos > maxz_m || zpos < minz_m || rpos > maxr_m || rpos < minr_m){
        flagNeedUpdate = true;
        gmsgALL << level4 << getName()
                << ": Particle " << id
                << " out of the global aperture of cyclotron!" << endl;
        gmsgALL << level4 << getName()
ext-calvo_p's avatar
ext-calvo_p committed
379
                << ": Coords: "<< RefPartBunch_m->R[id] << " m"  << endl;
380 381 382 383 384 385 386 387

    } else{
        flagNeedUpdate = apply(RefPartBunch_m->R[id], RefPartBunch_m->P[id], t, E, B);
        if(flagNeedUpdate){
            gmsgALL << level4 << getName()
                    << ": Particle "<< id
                    << " out of the field map boundary!" << endl;
            gmsgALL << level4 << getName()
ext-calvo_p's avatar
ext-calvo_p committed
388
                    << ": Coords: "<< RefPartBunch_m->R[id] << " m" << endl;
389 390 391 392 393 394 395 396 397 398
        }
    }

    if (flagNeedUpdate) {
        lossDs_m->addParticle(RefPartBunch_m->R[id], RefPartBunch_m->P[id],
                              id, t, 0, RefPartBunch_m->bunchNum[id]);
        RefPartBunch_m->Bin[id] = -1;
    }

    return flagNeedUpdate;
gsell's avatar
gsell committed
399 400
}

401
bool Cyclotron::apply(const Vector_t &R, const Vector_t &/*P*/,
gsell's avatar
gsell committed
402
                      const double &t, Vector_t &E, Vector_t &B) {
gsell's avatar
gsell committed
403

snuverink_j's avatar
snuverink_j committed
404
    const double rad   = std::hypot(R[0],R[1]);
gsell's avatar
gsell committed
405
    const double tempv = atan(R[1] / R[0]);
406
    double tet = tempv;
gsell's avatar
gsell committed
407 408 409

    /* if((R[0] > 0) && (R[1] >= 0)) tet = tempv;
       else*/
410
    if     ((R[0] < 0) && (R[1] >= 0)) tet = pi + tempv;
gsell's avatar
gsell committed
411 412 413 414 415 416 417
    else if((R[0] < 0) && (R[1] <= 0)) tet = pi + tempv;
    else if((R[0] > 0) && (R[1] <= 0)) tet = 2.0 * pi + tempv;
    else if((R[0] == 0) && (R[1] > 0)) tet = pi / 2.0;
    else if((R[0] == 0) && (R[1] < 0)) tet = 1.5 * pi;

    double tet_rad = tet;

418
    // the actual angle of particle in degree
gsell's avatar
gsell committed
419 420
    tet = tet / pi * 180.0;

421
    // Necessary for gap phase output -DW
422
    if (0 <= tet && tet <= 45) waiting_for_gap = 1;
423 424 425 426 427 428
    
    // dB_{z}/dr, dB_{z}/dtheta, B_{z}
    double brint = 0.0, btint = 0.0, bzint = 0.0;
    
    if ( this->interpolate(rad, tet_rad, brint, btint, bzint) ) {
        
gsell's avatar
gsell committed
429
        /* Br */
430 431
        double br = - brint * R[2];
        
gsell's avatar
gsell committed
432
        /* Btheta */
433 434 435
        double bt = - btint / rad * R[2];
        
        /* Bz */
436
        double bz = - bzint;
437

438
        this->applyTrimCoil(rad, R[2], tet_rad, br, bz);
439
        
440
        /* Br Btheta -> Bx By */
gsell's avatar
gsell committed
441 442 443 444
        B[0] = br * cos(tet_rad) - bt * sin(tet_rad);
        B[1] = br * sin(tet_rad) + bt * cos(tet_rad);
        B[2] = bz;
    } else {
445
        return true;
gsell's avatar
gsell committed
446
    }
447

448 449 450 451 452
    if(myBFieldType_m != SYNCHRO && myBFieldType_m != BANDRF) {
        return false;
    }

    //The RF field is supposed to be sampled on a cartesian grid
gsell's avatar
gsell committed
453 454 455 456 457 458 459 460
    std::vector<Fieldmap *>::const_iterator fi  = RFfields_m.begin();
    std::vector<double>::const_iterator rffi    = rffrequ_m.begin();
    std::vector<double>::const_iterator rfphii  = rfphi_m.begin();
    std::vector<double>::const_iterator escali  = escale_m.begin();
    std::vector<bool>::const_iterator superposei = superpose_m.begin();
    std::vector<std::vector<double>>::const_iterator rffci = rffc_m.begin();
    std::vector<std::vector<double>>::const_iterator rfvci = rfvc_m.begin();

461 462 463 464 465
    double xBegin(0), xEnd(0), yBegin(0), yEnd(0), zBegin(0), zEnd(0);
    int fcount = 0;

    for(; fi != RFfields_m.end(); ++fi, ++rffi, ++rfphii, ++escali, ++superposei) {
        (*fi)->getFieldDimensions(xBegin, xEnd, yBegin, yEnd, zBegin, zEnd);
gsell's avatar
gsell committed
466
        if (fcount > 0 && *superposei == false) continue;
467

468 469 470
        // Ok, this is a total patch job, but now that the internal cyclotron units are in m, we have to
        // change stuff here to match with the input units of mm in the fieldmaps. -DW
        const Vector_t temp_R = R * Vector_t(1000.0); //Keep this until we have transitioned fully to m -DW
471

gsell's avatar
gsell committed
472 473 474
        if (temp_R(0) < xBegin || temp_R(0) > xEnd ||
            temp_R(1) < yBegin || temp_R(1) > yEnd ||
            temp_R(2) < zBegin || temp_R(2) > zEnd) continue;
475

476
        Vector_t tmpE(0.0, 0.0, 0.0), tmpB(0.0, 0.0, 0.0);
gsell's avatar
gsell committed
477 478
        // out of bounds?
        if ((*fi)->getFieldstrength(temp_R, tmpE, tmpB)) continue;
479

480
        ++fcount;
481

482 483
        double frequency = (*rffi);   // frequency in MHz
        double ebscale = (*escali);   // E and B field scaling
484

485 486 487 488 489 490 491 492 493 494
        if(myBFieldType_m == SYNCHRO) {
            double powert = 1;
            for(const double fcoef : (*rffci)) {
                powert *= (t * 1e-9);
                frequency += fcoef * powert; // Add frequency ramp (in MHz/s^n)
            }
            powert = 1;
            for(const double vcoef : (*rfvci)) {
                powert *= (t * 1e-9);
                ebscale += vcoef * powert; // Add frequency ramp (in MHz/s^n)
495
            }
496
        }
497

498
        double phase = 2.0 * pi * 1.0E-3 * frequency * t + (*rfphii);  // f in [MHz], t in [ns]
499

500 501
        E += ebscale * cos(phase) * tmpE;
        B -= ebscale * sin(phase) * tmpB;
502

503 504
        if(myBFieldType_m != SYNCHRO)
            continue;
505

506 507 508 509
        // Some phase output -DW
        if (tet >= 90.0 && waiting_for_gap == 1) {
            double phase_print = 180.0 * phase / pi;
            phase_print = fmod(phase_print, 360) - 360.0;
510

511 512 513 514
            *gmsg << endl << "Gap 1 phase = " << phase_print << " Deg" << endl;
            *gmsg << "Gap 1 E-Field = (" << E[0] << "/" << E[1] << "/" << E[2] << ")" << endl;
            *gmsg << "Gap 1 B-Field = (" << B[0] << "/" << B[1] << "/" << B[2] << ")" << endl;
            *gmsg << "RF Frequency = " << frequency << " MHz" << endl;
515

516 517 518
            waiting_for_gap = 2;
        }
        else if (tet >= 270.0 && waiting_for_gap == 2) {
519

520 521
            double phase_print = 180.0 * phase / pi;
            phase_print = fmod(phase_print, 360) - 360.0;
522

523 524 525 526 527
            *gmsg << endl << "Gap 2 phase = " << phase_print << " Deg" << endl;
            *gmsg << "Gap 2 E-Field = (" << E[0] << "/" << E[1] << "/" << E[2] << ")" << endl;
            *gmsg << "Gap 2 B-Field = (" << B[0] << "/" << B[1] << "/" << B[2] << ")" << endl;
            *gmsg << "RF Frequency = " << frequency << " MHz" << endl;
            waiting_for_gap = 0;
528
        }
gsell's avatar
gsell committed
529 530 531
        if(myBFieldType_m == SYNCHRO) {
            ++rffci, ++rfvci;
        }
gsell's avatar
gsell committed
532 533 534 535
    }
    return false;
}

536 537 538 539
void Cyclotron::apply(const double& rad, const double& z,
                      const double& tet_rad, double& br,
                      double& bt, double& bz) {
    this->interpolate(rad, tet_rad, br, bt, bz);
540
    this->applyTrimCoil(rad, z, tet_rad, br, bz);
541 542
}

gsell's avatar
gsell committed
543
void Cyclotron::finalise() {
544

gsell's avatar
gsell committed
545
    online_m = false;
adelmann's avatar
adelmann committed
546
    lossDs_m->save();
adelmann's avatar
Cleanup  
adelmann committed
547
    *gmsg << "* Finalize cyclotron" << endl;
548

gsell's avatar
gsell committed
549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 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 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668
}

bool Cyclotron::bends() const {
    return true;
}

// calculate derivatives with 5-point lagrange's formula.
double Cyclotron::gutdf5d(double *f, double dx, const int kor, const int krl, const int lpr)

{
    double C[5][5][3], FAC[3];
    double result;
    int j;
    /* CALCULATE DERIVATIVES WITH 5-POINT LAGRANGE FORMULA
     * PARAMETERS:
     * F  STARTADDRESS FOR THE 5 SUPPORT POINTS
     * DX STEPWIDTH FOR ARGUMENT
     * KOR        ORDER OF DERIVATIVE (KOR=1,2,3).
     * KRL        NUMBER OF SUPPORT POINT, WHERE THE DERIVATIVE IS TO BE CALCULATED
     *  (USUALLY 3, USE FOR BOUNDARY 1 ,2, RESP. 4, 5)
     * LPR        DISTANCE OF THE 5 STORAGE POSITIONS (=1 IF THEY ARE NEIGHBORS OR LENGTH
     * OF COLUMNLENGTH OF A MATRIX, IF THE SUPPORT POINTS ARE ON A LINE).
     * ATTENTION! THE INDICES ARE NOW IN C-FORMAT AND NOT IN FORTRAN-FORMAT.*/

    /* COEFFICIENTS FOR THE 1ST DERIVATIVE: */
    C[0][0][0] = -50.0;
    C[1][0][0] = 96.0;
    C[2][0][0] = -72.0;
    C[3][0][0] = 32.0;
    C[4][0][0] = -6.0;
    C[0][1][0] = -6.0;
    C[1][1][0] = -20.0;
    C[2][1][0] = 36.0;
    C[3][1][0] = -12.0;
    C[4][1][0] =  2.0;
    C[0][2][0] =  2.0;
    C[1][2][0] = -16.0;
    C[2][2][0] =  0.0;
    C[3][2][0] = 16.0;
    C[4][2][0] = -2.0;
    C[0][3][0] = -2.0;
    C[1][3][0] = 12.0;
    C[2][3][0] = -36.0;
    C[3][3][0] = 20.0;
    C[4][3][0] =  6.0;
    C[0][4][0] =  6.0;
    C[1][4][0] = -32.0;
    C[2][4][0] = 72.0;
    C[3][4][0] = -96.0;
    C[4][4][0] = 50.0;

    /* COEFFICIENTS FOR THE 2ND DERIVATIVE: */
    C[0][0][1] = 35.0;
    C[1][0][1] = -104;
    C[2][0][1] = 114.0;
    C[3][0][1] = -56.0;
    C[4][0][1] = 11.0;
    C[0][1][1] = 11.0;
    C[1][1][1] = -20.0;
    C[2][1][1] =  6.0;
    C[3][1][1] =  4.0;
    C[4][1][1] = -1.0;
    C[0][2][1] = -1.0;
    C[1][2][1] = 16.0;
    C[2][2][1] = -30.0;
    C[3][2][1] = 16.0;
    C[4][2][1] = -1.0;
    C[0][3][1] = -1.0;
    C[1][3][1] =  4.0;
    C[2][3][1] =  6.0;
    C[3][3][1] = -20.0;
    C[4][3][1] = 11.0;
    C[0][4][1] = 11.0;
    C[1][4][1] = -56.0;
    C[2][4][1] = 114.0;
    C[3][4][1] = -104;
    C[4][4][1] = 35.0;


    /* COEFFICIENTS FOR THE 3RD DERIVATIVE: */
    C[0][0][2] = -10.0;
    C[1][0][2] = 36.0;
    C[2][0][2] = -48.0;
    C[3][0][2] = 28.0;
    C[4][0][2] = -6.0;
    C[0][1][2] = -6.0;
    C[1][1][2] = 20.0;
    C[2][1][2] = -24.0;
    C[3][1][2] = 12.0;
    C[4][1][2] = -2.0;
    C[0][2][2] = -2.0;
    C[1][2][2] =  4.0;
    C[2][2][2] =  0.0;
    C[3][2][2] = -4.0;
    C[4][2][2] =  2.0;
    C[0][3][2] =  2.0;
    C[1][3][2] = -12.0;
    C[2][3][2] = 24.0;
    C[3][3][2] = -20.0;
    C[4][3][2] =  6.0;
    C[0][4][2] =  6.0;
    C[1][4][2] = -28.0;
    C[2][4][2] = 48.0;
    C[3][4][2] = -36.0;
    C[4][4][2] = 10.0;

    /* FACTOR: */
    FAC[0] = 24.0;
    FAC[1] = 12.0;
    FAC[2] = 4.0;

    result = 0.0;
    for(j = 0; j < 5; j++) {
        result += C[j][krl][kor] * *(f + j * lpr);
    }

    return result / (FAC[kor] * pow(dx, (kor + 1)));
}


669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686
bool Cyclotron::interpolate(const double& rad,
                            const double& tet_rad,
                            double& brint,
                            double& btint,
                            double& bzint)
{
    const double xir = (rad - BP.rmin) / BP.delr;

    // ir : the number of path whose radius is less than the 4 points of cell which surround the particle.
    const int ir = (int)xir;

    // wr1 : the relative distance to the inner path radius
    const double wr1 = xir - (double)ir;
    // wr2 : the relative distance to the outer path radius
    const double wr2 = 1.0 - wr1;
    
    // the corresponding angle on the field map
    // Note: this does not work if the start point of field map does not equal zero.
687
    double tet_map = fmod(tet_rad / Physics::pi * 180.0, 360.0 / symmetry_m);
688 689 690 691 692 693 694 695
    
    double xit = tet_map / BP.dtet;

    int it = (int) xit;

    const double wt1 = xit - (double)it;
    const double wt2 = 1.0 - wt1;

gsell's avatar
gsell committed
696 697
    // it : the number of point on the inner path whose angle is less than
    // the particle' corresponding angle.
698
    // include zero degree point
gsell's avatar
gsell committed
699
    it++;
700 701 702 703 704 705 706

    int r1t1, r2t1, r1t2, r2t2;
    int ntetS = Bfield.ntet + 1;

    // r1t1 : the index of the "min angle, min radius" point in the 2D field array.
    // considering  the array start with index of zero, minus 1.

ext-rogers_c's avatar
ext-rogers_c committed
707
    if(myBFieldType_m != FFABF) {
708
        /*
ext-rogers_c's avatar
ext-rogers_c committed
709
          For FFA this does not work
710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729
        */
        r1t1 = it + ntetS * ir - 1;
        r1t2 = r1t1 + 1;
        r2t1 = r1t1 + ntetS;
        r2t2 = r2t1 + 1 ;

    } else {
        /*
          With this we have B-field AND this is far more
          intuitive for me ....
        */
        r1t1 = idx(ir, it);
        r2t1 = idx(ir + 1, it);
        r1t2 = idx(ir, it + 1);
        r2t2 = idx(ir + 1, it + 1);
    }

    if((it >= 0) && (ir >= 0) && (it < Bfield.ntetS) && (ir < Bfield.nrad)) {

        // B_{z}
730 731 732 733
        double bzf = Bfield.bfld[r1t1] * wr2 * wt2 +
                     Bfield.bfld[r2t1] * wr1 * wt2 +
                     Bfield.bfld[r1t2] * wr2 * wt1 +
                     Bfield.bfld[r2t2] * wr1 * wt1;
734

735
        bzint = /*- */bzf ;
736 737

        // dB_{z}/dr
738 739 740 741
        brint = Bfield.dbr[r1t1] * wr2 * wt2 +
                Bfield.dbr[r2t1] * wr1 * wt2 +
                Bfield.dbr[r1t2] * wr2 * wt1 +
                Bfield.dbr[r2t2] * wr1 * wt1;
742 743 744


        // dB_{z}/dtheta
745 746 747 748
        btint = Bfield.dbt[r1t1] * wr2 * wt2 +
                Bfield.dbt[r2t1] * wr1 * wt2 +
                Bfield.dbt[r1t2] * wr2 * wt1 +
                Bfield.dbt[r2t2] * wr1 * wt1;
749 750 751 752 753 754 755 756

        return true;
    }
    return false;
}


void Cyclotron::read(const int &fieldflag, const double &scaleFactor) {
ext-rogers_c's avatar
ext-rogers_c committed
757
    //    PSIBF, AVFEQBF, ANSYSBF, FFABF
758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780
    // for your own format field, you should add your own getFieldFromFile() function by yourself.

    if(fieldflag == 1) {
        //*gmsg<<"Read field data from PSI format field map file."<<endl;
        myBFieldType_m = PSIBF;
        getFieldFromFile(scaleFactor);

    } else if(fieldflag == 2) {
        // *gmsg<<"Read data from 450MeV Carbon cyclotron field file"<<endl
        myBFieldType_m = CARBONBF;
        getFieldFromFile_Carbon(scaleFactor);

    } else if(fieldflag == 3) {
        // *gmsg<<"Read data from 100MeV H- cyclotron CYCIAE-100 field file"<<endl;
        myBFieldType_m = ANSYSBF;
        getFieldFromFile_CYCIAE(scaleFactor);

    } else if(fieldflag == 4) {
        *gmsg << "* Read AVFEQ data (Riken) use bfield scale factor bs = " << getBScale() << endl;
        myBFieldType_m = AVFEQBF;
        getFieldFromFile_AVFEQ(scaleFactor);

    } else if(fieldflag == 5) {
ext-rogers_c's avatar
ext-rogers_c committed
781 782 783
        *gmsg << "* Read FFA data MSU/FNAL " << getBScale() << endl;
        myBFieldType_m = FFABF;
        getFieldFromFile_FFA(scaleFactor);
784 785

    } else if(fieldflag == 6) {
786
        *gmsg << "* Read both median plane B field map and 3D E field map of RF cavity for compact cyclotron" << endl;
787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804
        myBFieldType_m = BANDRF;
        getFieldFromFile_BandRF(scaleFactor);

    } else if(fieldflag == 7) {
        *gmsg << "* Read midplane B-field, 3D RF fieldmaps, and text files with RF frequency/Voltage coefficients for Synchrocyclotron. (Midplane scaling = " << getBScale() << ")" << endl;
        myBFieldType_m = SYNCHRO;
        getFieldFromFile_Synchrocyclotron(scaleFactor);

    } else
        ERRORMSG("* The field reading function of this TYPE of CYCLOTRON has not implemented yet!" << endl);

    // calculate the radii of initial grid.
    initR(BP.rmin, BP.delr, Bfield.nrad);

    // calculate the remaining derivatives
    getdiffs();
}

snuverink_j's avatar
snuverink_j committed
805
// evaluate other derivative of magnetic field.
gsell's avatar
gsell committed
806 807
void Cyclotron::getdiffs() {

808 809 810 811 812 813 814 815 816
    Bfield.dbr.resize(Bfield.ntot);
    Bfield.dbrr.resize(Bfield.ntot);
    Bfield.dbrrr.resize(Bfield.ntot);

    Bfield.dbrt.resize(Bfield.ntot);
    Bfield.dbrrt.resize(Bfield.ntot);
    Bfield.dbrtt.resize(Bfield.ntot);

    Bfield.f2.resize(Bfield.ntot);
817 818
    Bfield.f3.resize(Bfield.ntot);
    Bfield.g3.resize(Bfield.ntot);
gsell's avatar
gsell committed
819 820 821 822 823 824 825 826 827

    for(int i = 0; i < Bfield.nrad; i++) {

        for(int k = 0; k < Bfield.ntet; k++) {

            double dtheta = pi / 180.0 * BP.dtet;

            int kEdge;

gsell's avatar
gsell committed
828 829
            kEdge = std::max(k - 2, 0);
            kEdge = std::min(kEdge, Bfield.ntet - 5);
gsell's avatar
gsell committed
830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849

            int dkFromEdge = k - kEdge;
            int index = idx(i, k);
            int indexkEdge = idx(i, kEdge);


            Bfield.dbt[index]    = gutdf5d(&Bfield.bfld[indexkEdge], dtheta, 0, dkFromEdge, 1);
            Bfield.dbtt[index]   = gutdf5d(&Bfield.bfld[indexkEdge], dtheta, 1, dkFromEdge, 1);
            Bfield.dbttt[index]  = gutdf5d(&Bfield.bfld[indexkEdge], dtheta, 2, dkFromEdge, 1);
        }
    }



    for(int k = 0; k < Bfield.ntet; k++) {
        // inner loop varies R
        for(int i = 0; i < Bfield.nrad; i++) {
            double rac = BP.rarr[i];
            // define iredg, the reference index for radial interpolation
            // standard: i-2 minimal: 0 (not negative!)  maximal: nrad-4
gsell's avatar
gsell committed
850 851
            int iredg = std::max(i - 2, 0);
            iredg = std::min(iredg, Bfield.nrad - 5);
gsell's avatar
gsell committed
852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903
            int irtak = i - iredg;
            int index = idx(i, k);
            int indexredg = idx(iredg, k);


            Bfield.dbr[index]    = gutdf5d(&Bfield.bfld[indexredg], BP.delr, 0, irtak, Bfield.ntetS);
            Bfield.dbrr[index]   = gutdf5d(&Bfield.bfld[indexredg], BP.delr, 1, irtak, Bfield.ntetS);
            Bfield.dbrrr[index]  = gutdf5d(&Bfield.bfld[indexredg], BP.delr, 2, irtak, Bfield.ntetS);

            Bfield.dbrt[index]   = gutdf5d(&Bfield.dbt[indexredg], BP.delr, 0, irtak, Bfield.ntetS);
            Bfield.dbrrt[index]  = gutdf5d(&Bfield.dbt[indexredg], BP.delr, 1, irtak, Bfield.ntetS);
            Bfield.dbrtt[index]  = gutdf5d(&Bfield.dbtt[indexredg], BP.delr, 0, irtak, Bfield.ntetS);

            // fehlt noch!! f2,f3,g3,
            Bfield.f2[index] = (Bfield.dbrr[index]
                                + Bfield.dbr[index] / rac
                                + Bfield.dbtt[index] / rac / rac) / 2.0;

            Bfield.f3[index] = (Bfield.dbrrr[index]
                                + Bfield.dbrr[index] / rac
                                + (Bfield.dbrtt[index] - Bfield.dbr[index]) / rac / rac
                                - 2.0 * Bfield.dbtt[index] / rac / rac / rac) / 6.0;

            Bfield.g3[index] = (Bfield.dbrrt[index]
                                + Bfield.dbrt[index] / rac
                                + Bfield.dbttt[index] / rac / rac) / 6.0;
        } // Radius Loop
    } // Azimuth loop

    // copy 1st azimuth to last + 1 to always yield an interval
    for(int i = 0; i < Bfield.nrad; i++) {
        int iend = idx(i, Bfield.ntet);
        int istart = idx(i, 0);

        Bfield.bfld[iend]   = Bfield.bfld[istart];
        Bfield.dbt[iend]    = Bfield.dbt[istart];
        Bfield.dbtt[iend]   = Bfield.dbtt[istart];
        Bfield.dbttt[iend]  = Bfield.dbttt[istart];

        Bfield.dbr[iend]    = Bfield.dbr[istart];
        Bfield.dbrr[iend]   = Bfield.dbrr[istart];
        Bfield.dbrrr[iend]  = Bfield.dbrrr[istart];

        Bfield.dbrt[iend]   = Bfield.dbrt[istart];
        Bfield.dbrtt[iend]  = Bfield.dbrtt[istart];
        Bfield.dbrrt[iend]  = Bfield.dbrrt[istart];

        Bfield.f2[iend]     = Bfield.f2[istart];
        Bfield.f3[iend]     = Bfield.f3[istart];
        Bfield.g3[iend]     = Bfield.g3[istart];

    }
904

905
    /* debug
gsell's avatar
gsell committed
906 907 908 909

    for(int i = 0; i< Bfield.nrad; i++){
      for(int j = 0; j< Bfield.ntetS; j++){
    int index = idx(i,j);
910 911
    double x = (BP.rmin+i*BP.delr) * sin(j*BP.dtet*pi/180.0);
    double y = (BP.rmin+i*BP.delr) * cos(j*BP.dtet*pi/180.0);
gsell's avatar
gsell committed
912 913 914 915 916 917 918 919 920 921 922 923 924 925 926 927
    *gmsg<<"x= "<<x<<" y= "<<y<<" B= "<<Bfield.bfld[index]<<endl;
      }
    }
    */
}

// read field map from external file.
void Cyclotron::getFieldFromFile(const double &scaleFactor) {

    FILE *f = NULL;
    int lpar;
    char fout[100];
    double dtmp;

    *gmsg << "* ----------------------------------------------" << endl;
    *gmsg << "*             READ IN RING FIELD MAP            " << endl;
928
    *gmsg << "*      (The first data block is useless)        " << endl;
gsell's avatar
gsell committed
929 930 931 932 933
    *gmsg << "* ----------------------------------------------" << endl;

    BP.Bfact = scaleFactor;

    if((f = fopen(fmapfn_m.c_str(), "r")) == NULL) {
934 935
        throw GeneralClassicException("Cyclotron::getFieldFromField",
                                      "failed to open file '" + fmapfn_m + "', please check if it exists");
gsell's avatar
gsell committed
936 937
    }

938
    CHECK_CYC_FSCANF_EOF(fscanf(f, "%lf", &BP.rmin));
gsell's avatar
gsell committed
939
    *gmsg << "* Minimal radius of measured field map: " << BP.rmin << " [mm]" << endl;
940
    BP.rmin *= 0.001;  // mm --> m
gsell's avatar
gsell committed
941

942
    CHECK_CYC_FSCANF_EOF(fscanf(f, "%lf", &BP.delr));
943
    //if the value is negative, the actual value is its reciprocal.
944
    if(BP.delr < 0.0) BP.delr = 1.0 / (-BP.delr);
gsell's avatar
gsell committed
945
    *gmsg << "* Stepsize in radial direction: " << BP.delr << " [mm]" << endl;
946
    BP.delr *= 0.001;  // mm --> m
gsell's avatar
gsell committed
947

948
    CHECK_CYC_FSCANF_EOF(fscanf(f, "%lf", &BP.tetmin));
gsell's avatar
gsell committed
949 950
    *gmsg << "* Minimal angle of measured field map: " << BP.tetmin << " [deg.]" << endl;

951
    CHECK_CYC_FSCANF_EOF(fscanf(f, "%lf", &BP.dtet));
952
    //if the value is negative, the actual value is its reciprocal.
gsell's avatar
gsell committed
953 954 955
    if(BP.dtet < 0.0) BP.dtet = 1.0 / (-BP.dtet);
    *gmsg << "* Stepsize in azimuth direction: " << BP.dtet << " [deg.]" << endl;

956
    for(int i = 0; i < 13; i++)
957
        CHECK_CYC_FSCANF_EOF(fscanf(f, "%s", fout));
gsell's avatar
gsell committed
958

959
    CHECK_CYC_FSCANF_EOF(fscanf(f, "%d", &Bfield.nrad));
gsell's avatar
gsell committed
960 961
    *gmsg << "* Index in radial direction: " << Bfield.nrad << endl;

962
    CHECK_CYC_FSCANF_EOF(fscanf(f, "%d", &Bfield.ntet));
gsell's avatar
gsell committed
963 964 965 966 967 968
    *gmsg << "* Index in azimuthal direction: " << Bfield.ntet << endl;

    Bfield.ntetS = Bfield.ntet + 1;
    *gmsg << "* Accordingly, total grid point along azimuth:  " << Bfield.ntetS << endl;

    for(int i = 0; i < 5; i++) {
969
        CHECK_CYC_FSCANF_EOF(fscanf(f, "%s", fout));
gsell's avatar
gsell committed
970
    }
971
    CHECK_CYC_FSCANF_EOF(fscanf(f, "%d", &lpar));
gsell's avatar
gsell committed
972 973 974
    // msg<< "READ"<<lpar<<" DATA ENTRIES"<<endl;

    for(int i = 0; i < 4; i++) {
975
        CHECK_CYC_FSCANF_EOF(fscanf(f, "%s", fout));
gsell's avatar
gsell committed
976 977 978
    }

    for(int i = 0; i < lpar; i++) {
979
        CHECK_CYC_FSCANF_EOF(fscanf(f, "%16lE", &dtmp));
gsell's avatar
gsell committed
980 981
    }
    for(int i = 0; i < 6; i++) {
982
        CHECK_CYC_FSCANF_EOF(fscanf(f, "%s", fout));
gsell's avatar
gsell committed
983 984 985
    }
    //*gmsg << "* READ FILE DESCRIPTION..." <<endl;
    for(int i = 0; i < 10000; i++) {
986
        CHECK_CYC_FSCANF_EOF(fscanf(f, "%s", fout));
gsell's avatar
gsell committed
987 988 989 990
        if(strcmp(fout, "LREC=") == 0)break;
    }

    for(int i = 0; i < 5; i++) {
991
        CHECK_CYC_FSCANF_EOF(fscanf(f, "%s", fout));
gsell's avatar
gsell committed
992 993 994 995
    }
    Bfield.ntot = idx(Bfield.nrad - 1, Bfield.ntet) + 1;
    //jjyang
    *gmsg << "* Total stored grid point number ( ntetS * nrad ) : " << Bfield.ntot << endl;
996

997 998 999 1000
    Bfield.bfld.resize(Bfield.ntot);
    Bfield.dbt.resize(Bfield.ntot);
    Bfield.dbtt.resize(Bfield.ntot);
    Bfield.dbttt.resize(Bfield.ntot);
gsell's avatar
gsell committed
1001

1002
    *gmsg << "* Read-in loop one block per radius" << endl;
1003
    *gmsg << "* Rescaling of the magnetic fields with factor: " << BP.Bfact << endl;
gsell's avatar
gsell committed
1004 1005 1006 1007
    for(int i = 0; i < Bfield.nrad; i++) {

        if(i > 0) {
            for(int dummy = 0; dummy < 6; dummy++) {
1008
                CHECK_CYC_FSCANF_EOF(fscanf(f, "%s", fout)); // INFO-LINE
gsell's avatar
gsell committed
1009 1010 1011
            }
        }
        for(int k = 0; k < Bfield.ntet; k++) {
1012
            CHECK_CYC_FSCANF_EOF(fscanf(f, "%16lE", &(Bfield.bfld[idx(i, k)])));
gsell's avatar
gsell committed
1013 1014 1015
            Bfield.bfld[idx(i, k)] *= BP.Bfact;
        }
        for(int k = 0; k < Bfield.ntet; k++) {
1016
            CHECK_CYC_FSCANF_EOF(fscanf(f, "%16lE", &(Bfield.dbt[idx(i, k)])));
gsell's avatar
gsell committed
1017 1018 1019
            Bfield.dbt[idx(i, k)] *= BP.Bfact;
        }
        for(int k = 0; k < Bfield.ntet; k++) {
1020
            CHECK_CYC_FSCANF_EOF(fscanf(f, "%16lE", &(Bfield.dbtt[idx(i, k)])));
gsell's avatar
gsell committed
1021 1022 1023
            Bfield.dbtt[idx(i, k)] *= BP.Bfact;
        }
        for(int k = 0; k < Bfield.ntet; k++) {
1024
            CHECK_CYC_FSCANF_EOF(fscanf(f, "%16lE", &(Bfield.dbttt[idx(i, k)])));
gsell's avatar
gsell committed
1025 1026 1027 1028 1029 1030 1031 1032 1033
            Bfield.dbttt[idx(i, k)] *= BP.Bfact;
        }
    }
    fclose(f);


    *gmsg << "* Field Map read successfully!" << endl << endl;
}

1034 1035


1036
// Calculates Radii of initial grid.
1037
// dimensions in [m]!
gsell's avatar
gsell committed
1038
void Cyclotron::initR(double rmin, double dr, int nrad) {
1039
    BP.rarr.resize(nrad);
gsell's avatar
gsell committed
1040 1041 1042 1043 1044 1045
    for(int i = 0; i < nrad; i++) {
        BP.rarr[i] = rmin + i * dr;
    }
    BP.delr = dr;
}

1046
void Cyclotron::initialise(PartBunchBase<double, 3> *bunch, double &/*startField*/, double &/*endField*/) {
1047 1048 1049
    RefPartBunch_m = bunch;
    online_m = true;
}
gsell's avatar
gsell committed
1050

1051
void Cyclotron::initialise(PartBunchBase<double, 3> *bunch, const int &fieldflag, const double &scaleFactor) {
gsell's avatar
gsell committed
1052
    RefPartBunch_m = bunch;
adelmann's avatar
adelmann committed
1053
    lossDs_m = std::unique_ptr<LossDataSink>(new LossDataSink(getName(), !Options::asciidump));
kraus's avatar
kraus committed
1054

1055
    this->read(fieldflag, scaleFactor);
gsell's avatar
gsell committed
1056 1057 1058
}


1059
void Cyclotron::getFieldFromFile_FFA(const double &/*scaleFactor*/) {
gsell's avatar
gsell committed
1060 1061

    /*
1062
      Field is read in from ascii file (COSY output) in the order:
gsell's avatar
gsell committed
1063 1064
      R(m) theta(Deg) x(m) y(m) Bz(T).

1065
      Theta is the fast varying variable
gsell's avatar
gsell committed
1066 1067 1068 1069 1070 1071 1072 1073 1074 1075 1076 1077 1078 1079

      2.0000   0.0  2.0000  0.0000      0.0000000000000000
      2.0000   1.0  1.9997  0.0349      0.0000000000000000
      2.0000   2.0  1.9988  0.0698      0.0000000000000000
      2.0000   3.0  1.9973  0.1047      0.0000000000000000

      ......
      <blank line>

      2.1000   0.0  2.1000  0.0000      0.0000000000000000
      2.1000   1.0  2.0997  0.0367      0.0000000000000000
      2.1000   2.0  2.0987  0.0733      0.0000000000000000
    */

gsell's avatar
gsell committed
1080 1081 1082 1083 1084 1085
    std::vector<double> rv;
    std::vector<double> thv;
    std::vector<double> xv;
    std::vector<double> yv;
    std::vector<double> bzv;
    std::vector<double>::iterator vit;
gsell's avatar
gsell committed
1086 1087

    *gmsg << "* ----------------------------------------------" << endl;
ext-rogers_c's avatar
ext-rogers_c committed
1088
    *gmsg << "*    READ IN FFA FIELD MAP     " << endl;
gsell's avatar
gsell committed
1089 1090
    *gmsg << "* ----------------------------------------------" << endl;

ext-rogers_c's avatar
ext-rogers_c committed
1091
    BP.Bfact = -10.0; // T->kG and H- for the current FNAL FFA
gsell's avatar
gsell committed
1092

gsell's avatar
gsell committed
1093
    std::ifstream file_to_read(fmapfn_m.c_str());
gsell's avatar
gsell committed
1094 1095 1096 1097 1098 1099 1100
    const int max_num_of_char_in_a_line = 128;
    const int num_of_header_lines = 1;

    // STEP2: SKIP ALL THE HEADER LINES
    for(int i = 0; i < num_of_header_lines; ++i)
        file_to_read.ignore(max_num_of_char_in_a_line, '\n');