gfield_multipoles.cc

Go to the documentation of this file.

// CLHEP units
#include "CLHEP/Units/PhysicalConstants.h"
 
// geant4 headers
#include "G4ThreeVector.hh"
 
// gfields
#include "gfield_multipoles.h"
#include "gfieldConventions.h"
 
// gemv
// #include "gutilities.h"
 
 
// tells the DLL how to create a GFieldFactory in each plugin .so/.dylib
extern "C" GField* GFieldFactory(const std::shared_ptr<GOptions>& g) { return static_cast<GField*>(new GField_MultipolesFactory(g)); }
 
 
// Implementation notes (non-Doxygen):
// - This follows common accelerator multipole conventions for transverse fields.
// - Strength is defined at a reference radius of 1 m (see r0 below).
// - References (background reading):
//   - https://cds.cern.ch/record/1333874/files/1.pdf
//   - https://uspas.fnal.gov/materials/12MSU/magnet_elements.pdf
//   - https://cas.web.cern.ch/sites/default/files/lectures/bruges-2009/wolski-1.pdf
void GField_MultipolesFactory::GetFieldValue(const double pos[3], G4double* bfield) const {
    // Evaluate an ideal multipole field at the query point in the lab frame.
    // The computation is performed in a magnet-local frame and then rotated back.
 
    // ======= Configuration / conventions =======
    // strength: Tesla at reference radius r0 for all multipole orders
    const double r0 = CLHEP::m; // reference radius; make this a member if you want
 
    // ======= Basic checks =======
    if (pole_number < 2 || (pole_number % 2) != 0) {
        log->error(
            gfields::ERR_WRONG_POLE_NUMBER,
            "Pole number must be an even integer >= 2 (2=dipole,4=quadrupole,...)");
    }
 
    // ======= Positions and local frame =======
    const G4ThreeVector x0(pos[0], pos[1], pos[2]);          // query point (lab)
    const G4ThreeVector x1(origin[0], origin[1], origin[2]); // magnet origin (lab)
 
    // shift to magnet-centered coordinates and "unroll" the magnet by -rotation_angle
    G4ThreeVector p = x0 - x1;
    if (rotaxis == 0) p.rotateX(-rotation_angle);
    else if (rotaxis == 1) p.rotateY(-rotation_angle);
    else if (rotaxis == 2) p.rotateZ(-rotation_angle);
 
    // Identify transverse plane (u,v) ⟂ to the axis, and the axial coordinate w (unused here)
    double u = 0.0, v = 0.0;
    switch (rotaxis) {
    case 0: u = p.y();
        v = p.z();
        break; // axis = X ⇒ transverse plane = (Y,Z)
    case 1: u = p.z();
        v = p.x();
        break; // axis = Y ⇒ transverse plane = (Z,X)
    case 2: u = p.x();
        v = p.y();
        break; // axis = Z ⇒ transverse plane = (X,Y)
    default: break;
    }
 
    const double r   = std::hypot(u, v); // transverse radius
    const double phi = std::atan2(v, u); // azimuth in transverse plane
 
    // ======= Axial (solenoid-like) mode if explicitly requested =======
    if (longitudinal) {
        // Uniform axial field aligned with rotaxis; not a multipole.
        G4ThreeVector B_local(0, 0, 0);
        switch (rotaxis) {
        case 0: B_local.setX(strength);
            break;
        case 1: B_local.setY(strength);
            break;
        case 2: B_local.setZ(strength);
            break;
        default: break;
        }
 
        // Roll back to lab
        G4ThreeVector B_lab = B_local;
        if (rotaxis == 0) B_lab.rotateX(+rotation_angle);
        else if (rotaxis == 1) B_lab.rotateY(+rotation_angle);
        else if (rotaxis == 2) B_lab.rotateZ(+rotation_angle);
 
        bfield[0] = B_lab.x();
        bfield[1] = B_lab.y();
        bfield[2] = B_lab.z();
 
        log->info(
            2, "Axial field mode (solenoid-like). Strength: %1 T, Field: (%2, %3, %4)", strength, bfield[0],
            bfield[1], bfield[2]);
        return;
    }
 
    // ======= Transverse multipole (standard accelerator definition) =======
    const int n = pole_number / 2; // 1=dipole, 2=quadrupole, 3=sextupole, ...
    const int a = n - 1;           // power of r
 
    // r^a scaling with reference radius r0; for dipole (a=0) this is 1
    double ra = 1.0;
    if (a > 0) {
        // If r=0 and a>0, |B|=0 (for ideal multipoles).
        if (r == 0.0) {
            bfield[0] = bfield[1] = bfield[2] = 0.0;
            return;
        }
        ra = std::pow(r / r0, static_cast<double>(a));
    }
 
    // "Normal" multipole angular dependence:
    // Bu = strength * r^(n-1) * cos((n-1)*phi)
    // Bv = strength * r^(n-1) * sin((n-1)*phi)
    // This yields a constant transverse dipole for n=1.
    const double Bu = strength * ra * std::cos(a * phi);
    const double Bv = strength * ra * std::sin(a * phi);
 
    // Place (Bu,Bv) into the correct transverse components of B_local; axial component = 0
    G4ThreeVector B_local(0, 0, 0);
    switch (rotaxis) {
    case 0: B_local.setY(Bu);
        B_local.setZ(Bv);
        break; // axis X → (Y,Z)
    case 1: B_local.setZ(Bu);
        B_local.setX(Bv);
        break; // axis Y → (Z,X)
    case 2: B_local.setX(Bu);
        B_local.setY(Bv);
        break; // axis Z → (X,Y)
    default: break;
    }
 
    // Rotate (roll) back to lab
    G4ThreeVector B_lab = B_local;
    if (rotaxis == 0) B_lab.rotateX(+rotation_angle);
    else if (rotaxis == 1) B_lab.rotateY(+rotation_angle);
    else if (rotaxis == 2) B_lab.rotateZ(+rotation_angle);
 
    // Output
    bfield[0] = B_lab.x();
    bfield[1] = B_lab.y();
    bfield[2] = B_lab.z();
 
    log->info(
        2,
        "Pole Number: %1, n: %2, Strength: %3, Requested at: (%4, %5, %6), Rotation angle: %7, Rotation "
        "axis: %8, longitudinal: %9, Field: (%10, %11, %12)", pole_number, n, strength, pos[0], pos[1], pos[2],
        rotation_angle, rotaxis, longitudinal, bfield[0], bfield[1], bfield[2]);
}
 
 
// Load and cache multipole parameters from the provided field definition (non-Doxygen summary).
void GField_MultipolesFactory::load_field_definitions(GFieldDefinition gfd) {
    gfield_definitions = gfd;
 
    pole_number          = get_field_parameter_int("pole_number");
    origin[0]            = get_field_parameter_double("vx");
    origin[1]            = get_field_parameter_double("vy");
    origin[2]            = get_field_parameter_double("vz");
    rotation_angle       = get_field_parameter_double("rotation_angle");
    auto rot_axis_option = gfield_definitions.field_parameters["rotaxis"];
    longitudinal         = false;
 
    // Longitudinal mode: return a uniform axial field aligned with rotaxis.
    if (gfield_definitions.field_parameters["longitudinal"] == "true") {
        longitudinal = true;
        log->info(1, "Longitudinal field");
    }
    else { log->info(1, "Transverse field"); }
 
    // Parse rotaxis as a single axis character.
    if (rot_axis_option == "X" || rot_axis_option == "x") { rotaxis = 0; }
    else if (rot_axis_option == "Y" || rot_axis_option == "y") { rotaxis = 1; }
    else if (rot_axis_option == "Z" || rot_axis_option == "z") { rotaxis = 2; }
    else {
        log->error(
            gfields::ERR_WRONG_FIELD_ROTATION,
            "GField_MultipolesFactory::load_field_definitions: Rotation axis %1 not supported. Exiting.",
            gfield_definitions.field_parameters["rotaxis"]);
    }
 
    log->info(1, "Rotation axis: %1", rotaxis);
    strength = get_field_parameter_double("strength");
}