BeamStripping.cpp 12.5 KB
Newer Older
1
//
2 3 4 5
// Class BeamStripping
//   Defines the abstract interface environment for
//   beam stripping physics.
//
6
// Copyright (c) 2018-2021, Pedro Calvo, CIEMAT, Spain
7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
// All rights reserved
//
// Implemented as part of the PhD thesis
// "Optimizing the radioisotope production of the novel AMIT
// superconducting weak focusing cyclotron"
//
// This file is part of OPAL.
//
// OPAL is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// You should have received a copy of the GNU General Public License
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
22 23 24
//

#include "AbsBeamline/BeamStripping.h"
25 26

#include "AbsBeamline/BeamlineVisitor.h"
27 28 29 30 31 32 33 34 35 36 37
#include "AbsBeamline/Cyclotron.h"
#include "Algorithms/PartBunchBase.h"
#include "Fields/Fieldmap.h"
#include "Physics/Physics.h"
#include "Solvers/BeamStrippingPhysics.hh"
#include "Solvers/ParticleMatterInteractionHandler.hh"
#include "Utilities/LogicalError.h"
#include "Utilities/GeneralClassicException.h"

#include <fstream>
#include <iostream>
38
#include <cmath>
39
#include <cstdio>
40

41
#include "Utility/Inform.h"
42

43
#define CHECK_BSTP_FSCANF_EOF(arg) if (arg == EOF)\
44 45 46
throw GeneralClassicException("BeamStripping::getPressureFromFile",\
                              "fscanf returned EOF at " #arg);

47
extern Inform* gmsg;
48 49 50 51 52 53


BeamStripping::BeamStripping():
    BeamStripping("")
{}

54
BeamStripping::BeamStripping(const BeamStripping& right):
55 56 57 58 59 60 61 62 63 64 65 66 67 68
    Component(right),
    gas_m(right.gas_m),
    pressure_m(right.pressure_m),
    pmapfn_m(right.pmapfn_m),
    pscale_m(right.pscale_m),
    temperature_m(right.temperature_m),
    stop_m(right.stop_m),
    minr_m(right.minr_m),
    maxr_m(right.maxr_m),
    minz_m(right.minz_m),
    maxz_m(right.maxz_m),
    parmatint_m(NULL) {
}

69
BeamStripping::BeamStripping(const std::string& name):
70
    Component(name),
71
    gas_m(ResidualGas::NOGAS),
72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88
    pressure_m(0.0),
    pmapfn_m(""),
    pscale_m(0.0),
    temperature_m(0.0),
    stop_m(true),
    minr_m(0.0),
    maxr_m(0.0),
    minz_m(0.0),
    maxz_m(0.0),
    parmatint_m(NULL) {
}

BeamStripping::~BeamStripping() {
    if (online_m)
        goOffline();
}

89 90

void BeamStripping::accept(BeamlineVisitor& visitor) const {
91 92 93 94 95 96 97 98
    visitor.visitBeamStripping(*this);
}

void BeamStripping::setPressure(double pressure) {
    pressure_m = pressure;
}

double BeamStripping::getPressure() const {
99
    if (pressure_m > 0.0)
100 101 102
        return pressure_m;
    else {
        throw LogicalError("BeamStripping::getPressure",
103
                           "Pressure must not be zero");
104 105 106 107 108 109 110 111
    }
}

void BeamStripping::setTemperature(double temperature) {
    temperature_m = temperature;
}

double BeamStripping::getTemperature() const {
112
    if (temperature_m > 0.0)
113 114 115 116 117 118 119 120 121 122 123
        return temperature_m;
    else {
        throw LogicalError("BeamStripping::getTemperature",
                           "Temperature must not be zero");
    }
}

void BeamStripping::setPressureMapFN(std::string p) {
    pmapfn_m = p;
}

124
std::string BeamStripping::getPressureMapFN() const {
125 126 127 128 129 130 131 132
    return pmapfn_m;
}

void BeamStripping::setPScale(double ps) {
    pscale_m = ps;
}

double BeamStripping::getPScale() const {
133
    if (pscale_m > 0.0)
134 135 136 137 138
        return pscale_m;
    else {
        throw LogicalError("BeamStripping::getPScale",
                           "PScale must be positive");
    }
139 140 141
}

void BeamStripping::setResidualGas(std::string gas) {
142 143 144 145 146
    if (gas == "AIR") {
        gas_m = ResidualGas::AIR;
    } else if (gas == "H2") {
        gas_m = ResidualGas::H2;
    }
147 148
}

149 150 151 152 153 154 155 156 157 158 159 160 161 162
ResidualGas BeamStripping::getResidualGas() const {
    return gas_m;
}

std::string BeamStripping::getResidualGasName() {
    switch (gas_m) {
    case ResidualGas::AIR:
        return "AIR";
    case ResidualGas::H2:
        return "H2";
    default:
       throw GeneralClassicException("BeamStripping::getResidualGasName",
                                      "Residual gas not found");
    }		
163 164 165 166 167 168 169 170 171 172 173
}

void BeamStripping::setStop(bool stopflag) {
    stop_m = stopflag;
}

bool BeamStripping::getStop() const {
    return stop_m;
}


174
bool BeamStripping::checkBeamStripping(PartBunchBase<double, 3>* bunch, Cyclotron* cycl) {
175 176 177 178 179 180 181 182 183

    bool flagNeedUpdate = false;

    Vector_t rmin, rmax;
    bunch->get_bounds(rmin, rmax);
    std::pair<Vector_t, double> boundingSphere;
    boundingSphere.first = 0.5 * (rmax + rmin);
    boundingSphere.second = euclidean_norm(rmax - boundingSphere.first);

184 185 186 187
    maxr_m = cycl->getMaxR();
    minr_m = cycl->getMinR();
    maxz_m = cycl->getMaxZ();
    minz_m = cycl->getMinZ();
188 189 190

    size_t tempnum = bunch->getLocalNum();
    for (unsigned int i = 0; i < tempnum; ++i) {
191
        int pflag = checkPoint(bunch->R[i](0), bunch->R[i](1), bunch->R[i](2));
192 193 194 195 196 197 198 199 200 201 202 203
        if ( (pflag != 0) && (bunch->Bin[i] != -1) )
            flagNeedUpdate = true;
    }
    reduce(&flagNeedUpdate, &flagNeedUpdate + 1, &flagNeedUpdate, OpBitwiseOrAssign());
    if (flagNeedUpdate && parmatint_m) {
        dynamic_cast<BeamStrippingPhysics*>(parmatint_m)->setCyclotron(cycl);
        parmatint_m->apply(bunch, boundingSphere);
    }
    return flagNeedUpdate;
}


204
void BeamStripping::initialise(PartBunchBase<double, 3>* bunch, double& startField, double& endField) {
205 206 207 208
    endField = startField + getElementLength();
    initialise(bunch, pscale_m);
}

209
void BeamStripping::initialise(PartBunchBase<double, 3>* bunch, const double& scaleFactor) {
210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225

    RefPartBunch_m = bunch;

    parmatint_m = getParticleMatterInteraction();

    if (pmapfn_m != "") {
        getPressureFromFile(scaleFactor);
        // calculate the radii of initial grid.
        initR(PP.rmin, PP.delr, PField.nrad);
    }

    goOnline(-1e6);

    // change the mass and charge to simulate real particles
    for (size_t i = 0; i < bunch->getLocalNum(); ++i) {
        bunch->M[i] = bunch->getM()*1E-9;
226
        bunch->Q[i] = bunch->getQ() * Physics::q_e;
227
        if (bunch->weHaveBins())
228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248
            bunch->Bin[bunch->getLocalNum()-1] = bunch->Bin[i];
    }
}

void BeamStripping::finalise() {
    if (online_m)
        goOffline();
    *gmsg << "* Finalize Beam Stripping" << endl;
}

void BeamStripping::goOnline(const double &) {
}

void BeamStripping::goOffline() {
    online_m = false;
}

bool BeamStripping::bends() const {
    return false;
}

249
void BeamStripping::getDimensions(double& /*zBegin*/, double& /*zEnd*/) const {}
250 251 252 253 254 255 256 257 258

ElementBase::ElementType BeamStripping::getType() const {
    return BEAMSTRIPPING;
}

std::string BeamStripping::getBeamStrippingShape() {
    return "BeamStripping";
}

259
int BeamStripping::checkPoint(const double& x, const double& y, const double& z) {
260
    int cn;
261
    double rpos = std::sqrt(x * x + y * y);
262
    if (z >= maxz_m || z <= minz_m || rpos >= maxr_m || rpos <= minr_m) {
263
        cn = 0;
264
    } else {
265
        cn = 1;
266
    }
267 268 269
    return (cn);  // 0 if out, and 1 if in
}

270
double BeamStripping::checkPressure(const double& x, const double& y) {
271 272 273

    double pressure = 0.0;

274
    if (pmapfn_m != "") {
275
        const double rad = std::sqrt(x * x + y * y);
276 277 278 279 280 281 282 283 284
        const double xir = (rad - PP.rmin) / PP.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;

285
        const double tempv = std::atan(y / x);
286
        double tet = tempv;
287 288 289 290 291
        if      ((x < 0) && (y >= 0)) tet = Physics::pi + tempv;
        else if ((x < 0) && (y <= 0)) tet = Physics::pi + tempv;
        else if ((x > 0) && (y <= 0)) tet = Physics::two_pi + tempv;
        else if ((x == 0) && (y > 0)) tet = Physics::pi / 2.0;
        else if ((x == 0) && (y < 0)) tet = 1.5 * Physics::pi;
292 293

        // the actual angle of particle
294
        tet = tet * Physics::rad2deg;
295 296 297

        // the corresponding angle on the field map
        // Note: this does not work if the start point of field map does not equal zero.
298
        double xit = tet / PP.dtet;
299 300 301 302 303
        int it = (int) xit;
        const double wt1 = xit - (double)it;
        const double wt2 = 1.0 - wt1;
        // it : the number of point on the inner path whose angle is less than the particle' corresponding angle.
        // include zero degree point
304 305
        it++;
        double epsilon = 0.06;
306
        if (tet > 360 - epsilon && tet < 360 + epsilon) it = 0;
307 308 309 310 311 312 313 314 315 316

        int r1t1, r2t1, r1t2, r2t2;
        // 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.

        r1t1 = idx(ir, it);
        r2t1 = idx(ir + 1, it);
        r1t2 = idx(ir, it + 1);
        r2t2 = idx(ir + 1, it + 1);

317
        if ((it >= 0) && (ir >= 0) && (it < PField.ntetS) && (ir < PField.nrad)) {
318
            pressure = (PField.pfld[r1t1] * wr2 * wt2 +
319 320 321
                        PField.pfld[r2t1] * wr1 * wt2 +
                        PField.pfld[r1t2] * wr2 * wt1 +
                        PField.pfld[r2t2] * wr1 * wt1);
322
            if (pressure <= 0.0) {
323 324 325 326
                *gmsg << level4 << getName() << ": Pressure data from file zero." << endl;
                *gmsg << level4 << getName() << ": Take constant value through BeamStripping::getPressure" << endl;
                pressure = getPressure();
            }
327
        } else if (ir >= PField.nrad) {
328 329 330
            *gmsg << level4 << getName() << ": Particle out of maximum radial position of pressure field map." << endl;
            *gmsg << level4 << getName() << ": Take constant value through BeamStripping::getPressure" << endl;
            pressure = getPressure();
331
        } else {
332 333 334
            throw GeneralClassicException("BeamStripping::checkPressure",
                                          "Pressure data not found");
        }
335
    } else {
336 337 338
        pressure = getPressure();
    }

339
    return pressure;
340 341
}

342
// Calculates radius of initial grid (dimensions in [m]!)
343 344
void BeamStripping::initR(double rmin, double dr, int nrad) {
    PP.rarr.resize(nrad);
345
    for (int i = 0; i < nrad; i++)
346 347 348 349 350 351
        PP.rarr[i] = rmin + i * dr;

    PP.delr = dr;
}

// Read pressure map from external file.
352
void BeamStripping::getPressureFromFile(const double& scaleFactor) {
353

354
    FILE* f = NULL;
355 356 357 358 359 360

    *gmsg << "* " << endl;
    *gmsg << "* Reading pressure field map " << pmapfn_m << endl;

    PP.Pfact = scaleFactor;

361
    if ((f = std::fopen(pmapfn_m.c_str(), "r")) == NULL) {
362
        throw GeneralClassicException("BeamStripping::getPressureFromFile",
363 364
                                      "failed to open file '" + pmapfn_m +
                                      "', please check if it exists");
365 366
    }

367
    CHECK_BSTP_FSCANF_EOF(std::fscanf(f, "%lf", &PP.rmin));
368 369
    *gmsg << "* --- Minimal radius of measured pressure map: " << PP.rmin << " [mm]" << endl;

370
    CHECK_BSTP_FSCANF_EOF(std::fscanf(f, "%lf", &PP.delr));
371
    //if the value is negative, the actual value is its reciprocal.
372
    if (PP.delr < 0.0) PP.delr = 1.0 / (-PP.delr);
373 374
    *gmsg << "* --- Stepsize in radial direction: " << PP.delr << " [mm]" << endl;

375
    CHECK_BSTP_FSCANF_EOF(std::fscanf(f, "%lf", &PP.tetmin));
376 377
    *gmsg << "* --- Minimal angle of measured pressure map: " << PP.tetmin << " [deg]" << endl;

378
    CHECK_BSTP_FSCANF_EOF(std::fscanf(f, "%lf", &PP.dtet));
379
    //if the value is negative, the actual value is its reciprocal.
380
    if (PP.dtet < 0.0) PP.dtet = 1.0 / (-PP.dtet);
381 382
    *gmsg << "* --- Stepsize in azimuthal direction: " << PP.dtet << " [deg]" << endl;

383
    CHECK_BSTP_FSCANF_EOF(std::fscanf(f, "%d", &PField.ntet));
384 385
    *gmsg << "* --- Grid points along azimuth (ntet): " << PField.ntet << endl;

386
    CHECK_BSTP_FSCANF_EOF(std::fscanf(f, "%d", &PField.nrad));
387
    *gmsg << "* --- Grid points along radius (nrad): " << PField.nrad << endl;
388 389 390
    *gmsg << "* --- Maximum radial position: " << PP.rmin + (PField.nrad-1)*PP.delr << " [mm]" << endl;
    PP.rmin *= 0.001;  // mm --> m
    PP.delr *= 0.001;  // mm --> m
391 392

    PField.ntetS = PField.ntet + 1;
393
    *gmsg << "* --- Adding a guard cell along azimuth" << endl;
394

395
    PField.ntot = PField.nrad * PField.ntetS;
396
    PField.pfld.resize(PField.ntot);
397
    *gmsg << "* --- Total stored grid point number ((ntet+1) * nrad) : " << PField.ntot << endl;
398 399
    *gmsg << "* --- Escale factor: " << PP.Pfact << endl;

400 401 402
    for (int i = 0; i < PField.nrad; i++) {
        for (int k = 0; k < PField.ntet; k++) {
            CHECK_BSTP_FSCANF_EOF(std::fscanf(f, "%16lE", &(PField.pfld[idx(i, k)])));
403 404 405 406 407
            PField.pfld[idx(i, k)] *= PP.Pfact;
        }
    }
    *gmsg << "*" << endl;

408
    std::fclose(f);
409 410
}

411
#undef CHECK_BSTP_FSCANF_EOF