Simple Helmholtz Coils calculator (C++, demonstrates SWIG)

hhcoils.cxx

/* HelmholtzCoils - demonstrating separation of analysis and toolkit interface in NDIToolbox by creating a
simple magnetic field calculator

Chris Coughlin (TRI/Austin, Inc.)
*/
#include "hhcoils.h"

const double HelmholtzCoils::H(double position) const{
    double lh_geometry_factor = geometry_correction(lhcoil_position,position);
    double rh_geometry_factor = geometry_correction(rhcoil_position,position);
    double total_geometry_factor = lh_geometry_factor + rh_geometry_factor;
    return ((N*I)/(2*a))*total_geometry_factor;
}

const double HelmholtzCoils::geometry_correction(const double coil_position, const double position) const{
    return pow(1 + pow(coil_position - position,2) / pow(a,2), -1.5);
}

const double HelmholtzCoils::wirelength(void) const {
    double coil_circumference = M_PI*2*a;
    double wire_per_coil = N*coil_circumference;
    return 2*wire_per_coil;
}

const double HelmholtzCoils::awg_recommendation(void) const {
    double wire_gauges[] = {10,11,14,17,20,21,22,24,27,30,31,32,34,37,40};
    double currentrecs_awg[] = {15,10,5,2.5,1.5,1.0,0.75,0.5,0.25,0.125,0.100,0.075,0.050,0.025,0.0125};
    double delta = DBL_MAX;
    double epsilon = 0;
    double awg_rec = 0;
    if(I>=0.0125 && I<=15) {
        for(int iter=0;iter<16;iter++) {
            epsilon = currentrecs_awg[iter] - I;
            if(epsilon>=0 && epsilon<delta) {
                delta = epsilon;
                awg_rec = wire_gauges[iter];
            }
        }
    }
    return awg_rec;
}

hhcoils.h

/* HelmholtzCoils - demonstrating separation of analysis and toolkit interface in NDIToolbox by creating a
simple magnetic field calculator

Chris Coughlin (TRI/Austin, Inc.)
*/

#ifndef HELMHOLTZCOILS_H_
#define HELMHOLTZCOILS_H_

#define _USE_MATH_DEFINES // Required on some platforms to get mathematical constants
#include <cmath>
#include <cfloat>

static const double mu_0 = 4*M_PI*1e-7;

class HelmholtzCoils {
public:
    HelmholtzCoils(int turns_per_coil, double current_per_coil, double coil_radius):
        N(turns_per_coil), I(current_per_coil), a(coil_radius), lhcoil_position(-a/2), rhcoil_position(a/2) { }
    
    const double H(const double position) const; // Magnetic field at position (m) in A/m
    const double centerH(void) const { return H(0); } // Magnetic field at dead center of coils
    const double B(const double position) const { return mu_0*H(position)*1000; } // Flux density at position (m) in mT
    const double centerB(void) const { return B(0); } // Flux density at dead center of coils
    const double B_mG(const double position) const { return B(position) * 1e4; } // Flux density at position (m) in mG
    
    const double wirelength(void) const; // Length of wire (m) to make N turns of radius a
    const double awg_recommendation(void) const; // Lookup table to make AWG recommendations based on current I
        
private:
    int N; // Turns per coil
    double I; // Current (A) per coil
    double a; // Common radius (m) of coils
    double lhcoil_position; // Position of left coil (m), defined as -a/2
    double rhcoil_position; // Position of right coil (m), defined as a/2
    // Geometry correction for magnetic field calcs
    const double geometry_correction(const double coil_position, const double position) const; 
};

#endif // HELMHOLTZCOILS_H_

hhcoils.i

%module hhcoils
%{
/* Includes the header in the wrapper code */
#include "hhcoils.h"
%}
/* Include various STL interfaces - not really needed here but included as a demo*/
%include "std_string.i"
%include "std_vector.i"
/* For vectors it's necessary to specify what type of vector(s) you'll be using-here we're
using the Python "vector_float" as an alias for vector<float>*/
namespace std {
   %template(vector_float) vector<float>;
   %template(vector_vector_float) vector<std::vector<float> >;
};

/* Parse the header file to generate wrappers */
%include "hhcoils.h"