diff --git a/libcon2020/modif/include/con2020.h b/libcon2020/modif/include/con2020.h new file mode 100644 index 0000000..29105d8 --- /dev/null +++ b/libcon2020/modif/include/con2020.h @@ -0,0 +1,670 @@ +#ifndef __LIBCON2020_H__ +#define __LIBCON2020_H__ +#include +#include +#include +#include +#include + +#define LIBCON2020_VERSION_MAJOR 0 +#define LIBCON2020_VERSION_MINOR 1 +#define LIBCON2020_VERSION_PATCH 0 + +#define M_PI 3.14159265358979323846 +#define USE_MATH_DEFINES +#define _USE_MATH_DEFINES +#define deg2rad M_PI/180.0 +#define rad2deg 180.0/M_PI + +extern "C" { + /* these wrappers can be used to get the magnetic field vectors */ + void Con2020FieldArray(int n, double *p0, double *p1, double *p2, + double *B0, double *B1, double *B2); + + void Con2020Field(double p0, double p1, double p2, + double *B0, double *B1, double *B2); + + + void GetCon2020Params(double *mui, double *irho, double *r0, double *r1, + double *d, double *xt, double *xp, char *eqtype, + bool *Edwards, bool *ErrChk, bool *CartIn, bool *CartOut, + bool *smooth, double *DeltaRho, double *DeltaZ, + double *g, char *azfunc, double *wO_open, double *wO_om, + double *thetamm, double *dthetamm, double *thetaoc, double *dthetaoc); + + + void SetCon2020Params(double mui, double irho, double r0, double r1, + double d, double xt, double xp, const char *eqtype, + bool Edwards, bool ErrChk, bool CartIn, bool CartOut, + bool smooth, double DeltaRho, double DeltaZ, + double g, const char *azfunc, double wO_open, double wO_om, + double thetamm, double dthetamm, double thetaoc, double dthetaoc); + + void Con2020AnalyticField( int n, double a, + double *rho, double *z, + double *Brho, double *Bz); + + void Con2020AnalyticFieldSmooth( int n, double a, + double *rho, double *z, + double *Brho, double *Bz); + + +/*************************************************************** +* +* NAME : ScalarPotentialSmallRho(rho,z,a,mui2,D) +* +* DESCRIPTION : Calcualte the small rho approximation +* of the scalar potential accoring to Edwards et al., +* 2001 (equation 8). +* +* INPUTS : +* double rho Cylindrical rho coordinate (in disc +* coordinate system, Rj) +* double z z-coordinate, Rj +* double a inner edge of semi-infinite current +* sheet, Rj +* double mui2 mu_0 I_0 /2 parameter (default 139.6 nT) +* double D Current sheet half-thickness, Rj +* +***************************************************************/ + double ScalarPotentialSmallRho( double rho, double z, double a, + double mui2, double D); + + +/*************************************************************** +* +* NAME : ScalarPotentialLargeRho(rho,z,a,mui2,D,deltaz) +* +* DESCRIPTION : Calcualte the large rho approximation +* of the scalar potential accoring to Edwards et al., +* 2001 (equation 12). +* +* INPUTS : +* double rho Cylindrical rho coordinate (in disc +* coordinate system, Rj) +* double z z-coordinate, Rj +* double a inner edge of semi-infinite current +* sheet, Rj +* double mui2 mu_0 I_0 /2 parameter (default 139.6 nT) +* double D Current sheet half-thickness, Rj +* double deltaz Scale length over which to smooth 4th +* term of the equation +* +***************************************************************/ + double ScalarPotentialLargeRho( double rho, double z, double a, + double mui2, double D, double deltaz); + + +/*************************************************************** +* +* NAME : ScalarPotential(rho,z,a,mui2,D,deltarho,deltaz) +* +* DESCRIPTION : Calculate the small/large rho approximation +* of the scalar potential accoring to Edwards et al., +* 2001 (equations 8 & 12). +* +* INPUTS : +* double rho Cylindrical rho coordinate (in disc +* coordinate system, Rj) +* double z z-coordinate, Rj +* double a inner edge of semi-infinite current +* sheet, Rj +* double mui2 mu_0 I_0 /2 parameter (default 139.6 nT) +* double D Current sheet half-thickness, Rj +* double deltarho Scale length to smoothly transition from +* small to large rho approx +* double deltaz Scale length over which to smooth 4th +* term of the equation +* +***************************************************************/ + double ScalarPotential( double rho, double z, double a, + double mui2, double D, + double deltarho, double deltaz); + + /************************************************************* + * + * NAME: f_theta(thetai) + * + * DESCRIPTION: Equation 5 of Cowley et al., 2008 + * + * INPUTS: + * double thetai colatitude of the ionospheric footprint + * in radians! + * + * RETURNS: + * double f_theta 1 + 0.25*tan^2 thetai + * + *************************************************************/ + double f_thetai(double thetai); + + /************************************************************* + * + * NAME: OmegaRatio(thetai,wO_open,wO_om,thetamm,dthetamm, + * thetaoc,dthetaoc) + * + * DESCRIPTION: Ratio of the angular velocity mapped to + * thetai to the planetary rotation. Equation 15 of + * Cowley et al., 2008. + * + * INPUTS: + * double thetai colatitude of the ionospheric footprint + * in radians! + * double wO_open angular velocity ratio of open flux to + * planetary spin + * double wO_om angular velocity ratio of outer magnetosphere + * to planetary spin + * double thetamm ionospheric footprint latitude of the + * middle magnetosphere (where plasma + * goes from rigid corotation to subcorotation) + * in radians. + * double dthetamm width of the middle magnetosphere in radians. + * double thetaoc ionospheric latitude of the open-closed field + * line boundary, in radians. + * double dthetaoc width of the open-closed field line boundary, + * in radians. + * + * RETURNS: + * double wO Ratio of plasma angular veloctiy to Jupiter + * spin. + * + *************************************************************/ + double OmegaRatio( double thetai, double wO_open, double wO_om, + double thetamm, double dthetamm, + double thetaoc, double dthetaoc); + + /************************************************************* + * + * NAME: PedersenCurrent(thetai,g,wO_open,wO_om,thetamm,dthetamm, + * thetaoc,dthetsoc) + * + * DESCRIPTION: Calculate the Pedersen current which maps to a + * given ionospheric latitude using equation 6 of Cowley et + * al., 2008. + * + * INPUTS: + * double thetai colatitude of the ionospheric footprint + * in radians! + * double g dipole coefficient, nT. + * double wO_open angular velocity ratio of open flux to + * planetary spin + * double wO_om angular velocity ratio of outer magnetosphere + * to planetary spin + * double thetamm ionospheric footprint latitude of the + * middle magnetosphere (where plasma + * goes from rigid corotation to subcorotation) + * in radians. + * double dthetamm width of the middle magnetosphere in radians. + * double thetaoc ionospheric latitude of the open-closed field + * line boundary, in radians. + * double dthetaoc width of the open-closed field line boundary, + * in radians. + * RETURNS: + * double Ihp Ionospheric Pedersen current. + * + *************************************************************/ + double PedersenCurrent( double thetai, double g, + double wO_open, double wO_om, + double thetamm, double dthetamm, + double thetaoc, double dthetaoc ); + + /************************************************************* + * + * NAME: ThetaIonosphere(r,theta,g,r0,r1,mui2,D,deltarho,deltaz) + * + * DESCRIPTION: Use the flux functions of the CAN model and a + * dipole field to map the current position to a position + * on the ionosphere. + * + * INPUTS: + * double r radial coordinate, Rj. + * double theta colatitude, radians. + * double g dipole coefficient, nT. + * double r0 Inner edge of the current sheet, Rj. + * double r1 Outer edge of the current sheet, Rj. + * double mui2 current parameter, nT. + * double D half-thickness of the current sheet, Rj. + * double deltarho scale distance of the smoothing between + * inner and outer approximations, Rj. + * double deltaz scale distance to smooth across the + * +/-D boundary, Rj. + * + * RETURNS: + * double thetai Ionospheric latitude in radians. + * + * + *************************************************************/ + double ThetaIonosphere( double r, double theta, double g, + double r0, double r1, + double mui2, double D, + double deltarho, double deltaz); + + /************************************************************* + * + * NAME: BphiLMIC(r,theta,g,r0,r1,mui2,D,deltarho,deltaz, + * wO_open,wO_om,thetamm,dthetamm, + * thetaom,dthetaom) + * + * DESCRIPTION: Calculate the azimuthal field using the LMIC + * model. + * + * INPUTS: + * double r radial coordinate, Rj. + * double theta colatitude, radians. + * double g dipole coefficient, nT. + * double r0 Inner edge of the current sheet, Rj. + * double r1 Outer edge of the current sheet, Rj. + * double mui2 current parameter, nT. + * double D half-thickness of the current sheet, Rj. + * double deltarho scale distance of the smoothing between + * inner and outer approximations, Rj. + * double deltaz scale distance to smooth across the + * +/-D boundary, Rj. + * double wO_open angular velocity ratio of open flux to + * planetary spin + * double wO_om angular velocity ratio of outer magnetosphere + * to planetary spin + * double thetamm ionospheric footprint latitude of the + * middle magnetosphere (where plasma + * goes from rigid corotation to subcorotation) + * in radians. + * double dthetamm width of the middle magnetosphere in radians. + * double thetaoc ionospheric latitude of the open-closed field + * line boundary, in radians. + * double dthetaoc width of the open-closed field line boundary, + * in radians. + * + * RETURNS: + * double Bphi Azimuthal field, nT. + * + *************************************************************/ + double BphiLMIC(double r, double theta, double g, + double r0, double r1, + double mui2, double D, + double deltarho, double deltaz, + double wO_open, double wO_om, + double thetamm, double dthetamm, + double thetaoc, double dthetaoc ); + + /************************************************************* + * + * NAME: BphiIonosphere(thetai,g,wO_open,wO_om,thetamm,dthetamm, + * thetaom,dthetaom) + * + * DESCRIPTION: Calculate the ionospheric azimuthal field using the LMIC + * model. + * + * INPUTS: + * double thetai ionospheric colatitude, radians. + * double g dipole coefficient, nT. + * double wO_open angular velocity ratio of open flux to + * planetary spin + * double wO_om angular velocity ratio of outer magnetosphere + * to planetary spin + * double thetamm ionospheric footprint latitude of the + * middle magnetosphere boundary (where plasma + * goes from rigid corotation to subcorotation) + * in radians. + * double dthetamm width of the middle magnetosphere boundary + * in radians. + * double thetaoc ionospheric latitude of the open-closed field + * line boundary, in radians. + * double dthetaoc width of the open-closed field line boundary, + * in radians. + * + * RETURNS: + * double Bphi Azimuthal field, nT. + * + *************************************************************/ + double BphiIonosphere( double thetai, double g, + double wO_open, double wO_om, + double thetamm, double dthetamm, + double thetaoc, double dthetaoc ); +} + + +/*********************************************************************** + * NAME : j0(x) + * + * DESCRIPTION : Fucntion to calculate an estimate of the Bessel + * function j0 using code based on the Cephes C library + * (Cephes Mathematical Functions Library, http://www.netlib.org/cephes/). + * + * INPUTS : + * double x position to calculate J0 at. + * + * RETURNS : + * double j j0 function evaluated at x. + * + * ********************************************************************/ +double j0(double x); + +/*********************************************************************** + * NAME : j1(x) + * + * DESCRIPTION : Fucntion to calculate an estimate of the Bessel + * function j1 using code based on the Cephes C library + * (Cephes Mathematical Functions Library, http://www.netlib.org/cephes/). + * + * INPUTS : + * double x position to calculate J1 at. + * + * RETURNS : + * double j j1 function evaluated at x. + * + * ********************************************************************/ +double j1(double x); + +/*********************************************************************** + * NAME : j0(n,x,j) + * + * DESCRIPTION : Fucntion to calculate an estimate of the Bessel + * function j0 using code based on the Cephes C library + * (Cephes Mathematical Functions Library, http://www.netlib.org/cephes/). + * + * INPUTS : + * int n Number of elements in x + * double *x position to calculate J0 at. + * + * OUTPUTS : + * double *j j0 function evaluated at x. + * + * ********************************************************************/ +void j0(int n, double *x, double *j); + +/*********************************************************************** + * NAME : j1(n,x,j) + * + * DESCRIPTION : Fucntion to calculate an estimate of the Bessel + * function j1 using code based on the Cephes C library + * (Cephes Mathematical Functions Library, http://www.netlib.org/cephes/). + * + * INPUTS : + * int n Number of elements in x + * double *x position to calculate J1 at. + * + * OUTPUTS : + * double *j j1 function evaluated at x. + * + * ********************************************************************/ +void j1(int n, double *x, double *j); + +/*********************************************************************** + * NAME : j0(n,x,multx,j) + * + * DESCRIPTION : Fucntion to calculate an estimate of the Bessel + * function j0 using code based on the Cephes C library + * (Cephes Mathematical Functions Library, http://www.netlib.org/cephes/). + * + * INPUTS : + * int n Number of elements in x + * double *x position to calculate J0(x*multx) at. + * double multx Constant to multiply x by + * + * OUTPUTS : + * double *j j0 function evaluated at x*multx. + * + * ********************************************************************/ +void j0(int n, double *x, double multx, double *j); + +/*********************************************************************** + * NAME : j1(n,x,multx,j) + * + * DESCRIPTION : Fucntion to calculate an estimate of the Bessel + * function j1 using code based on the Cephes C library + * (Cephes Mathematical Functions Library, http://www.netlib.org/cephes/). + * + * INPUTS : + * int n Number of elements in x + * double *x position to calculate J1(x*multx) at. + * double multx Constant to multiply x by + * + * OUTPUTS : + * double *j j1 function evaluated at x*multx. + * + * ********************************************************************/ +void j1(int n, double *x, double multx, double *j); + + +template T clip(T x, T mn, T mx) { + return std::min(mx,std::max(x,mn)); +} + +/* function pointer for input conversion */ +class Con2020; /*this is needed for the pointer below */ +typedef void (Con2020::*InputConvFunc)(int,double*,double*,double*, + double*,double*,double*,double*,double*, + double*,double*,double*,double*); +/* Output conversion */ +typedef void (Con2020::*OutputConvFunc)(int,double*,double*,double*, + double*,double*,double*,double*, + double*,double*,double*, + double*,double*,double*); + +/* Model function */ +typedef void (Con2020::*ModelFunc)(double,double,double,double*,double*,double*); + +/* analytical approximation equations */ +typedef void (Con2020::*Approx)(double,double,double,double,double,double*,double*); + +/* azimuthal function */ +typedef void (Con2020::*AzimFunc)(double,double,double,double*); + +class Con2020 { + public: + /* constructors */ + Con2020(); + Con2020(double,double,double,double,double,double,double,const char*,bool,bool,bool,bool); + + /* destructor */ + ~Con2020(); + + /* these functions will be used to set the equations used, if + * they need to be changed post-initialisation */ + + // BRE - 10/10/2023 - Add method to init integrals + void Init(); + //---------- + + void SetEdwardsEqs(bool); + void SetEqType(const char*); + void SetAzCurrentParameter(double); + void SetRadCurrentParameter(double); + void SetR0(double); + void SetR1(double); + void SetCSHalfThickness(double); + void SetCSTilt(double); + void SetCSTiltAzimuth(double); + void SetErrCheck(bool); + void SetCartIn(bool); + void SetCartOut(bool); + void SetSmooth(bool); + void SetDeltaRho(double); + void SetDeltaZ(double); + void SetOmegaOpen(double); + void SetOmegaOM(double); + void SetThetaMM(double); + void SetdThetaMM(double); + void SetThetaOC(double); + void SetdThetaOC(double); + void SetG(double); + void SetAzimuthalFunc(const char*); + + /* these mamber functions will be the "getter" version of the + * above setters */ + bool GetEdwardsEqs(); + void GetEqType(char*); + double GetAzCurrentParameter(); + double GetRadCurrentParameter(); + double GetR0(); + double GetR1(); + double GetCSHalfThickness(); + double GetCSTilt(); + double GetCSTiltAzimuth(); + bool GetErrCheck(); + bool GetCartIn(); + bool GetCartOut(); + bool GetSmooth(); + double GetDeltaRho(); + double GetDeltaZ(); + double GetOmegaOpen(); + double GetOmegaOM(); + double GetThetaMM(); + double GetdThetaMM(); + double GetThetaOC(); + double GetdThetaOC(); + double GetG(); + void GetAzimuthalFunc(char *); + + /* This function will be used to call the model, it is overloaded + * so that we have one for arrays, one for scalars */ + void Field(int,double*,double*,double*,double*,double*,double*); + void Field(double,double,double,double*,double*,double*); + + /* a function for testing purposes...*/ + void AnalyticField(double,double,double,double*,double*); + void AnalyticFieldSmooth(double,double,double,double*,double*); + + private: + /* model parameters */ + double mui_,irho_,r0_,r1_,d_,xt_,xp_,disctilt_,discshift_; + double r0sq_, r1sq_; + double cosxp_,sinxp_,cosxt_,sinxt_; + char eqtype_[9]; + bool Edwards_, ErrChk_; + bool CartIn_,CartOut_; + double deltaz_,deltarho_; + bool smooth_; + + /* LMIC parameters*/ + double wO_open_, wO_om_, thetamm_, dthetamm_, thetaoc_, dthetaoc_, g_; + char azfunc_[10]; + + /* Bessel function arrays - arrays prefixed with r and z are + * to be used for integrals which calcualte Brho and Bz, + * respectively */ + int *rnbes_; /* number of elements for each Bessel function (rho)*/ + int *znbes_; /* same as above for z */ + double **rlambda_;/* Lambda array to integrate over rho*/ + double **zlambda_;/* Lambda array to integrate over z*/ + double **rj0_lambda_r0_; /* j0(lambda*r0) */ + double **rj1_lambda_rho_;/* j1(lambda*rho) */ + double **zj0_lambda_r0_; /* j0(lambda*r0) */ + double **zj0_lambda_rho_;/* j0(lambda*rho) */ + + /* arrays to multiply be stuff to be integrated */ + /* these arrays will store the parts of equations 14, 15, 17 + * and 18 of Connerny 1981 which only need to be calculated once*/ + double **Eq14_; /* j0(lambda*r0)*sinh(lamba*d)/lambda */ + double **Eq15_; /* j0(lambda*r0)*sinh(lamba*d)/lambda */ + double **Eq17_; /* j0(lambda*r0)*exp(-lamba*d)/lambda */ + double **Eq18_; /* j0(lambda*r0)/lambda */ + double **ExpLambdaD_; + + /* integration step sizes */ + static constexpr double dlambda_ = 1e-4; + static constexpr double dlambda_brho_ = 1e-4; + static constexpr double dlambda_bz_ = 5e-5; + + /* Arrays containing maximum lambda values */ + double rlmx_array_[6]; + double zlmx_array_[6]; + + /* coordinate conversions for positions */ + InputConvFunc _ConvInput; + void _SysIII2Mag(int,double*,double*,double*, + double*,double*,double*,double*,double*, + double*,double*,double*,double*); + void _PolSysIII2Mag(int,double*,double*,double*, + double*,double*,double*,double*,double*, + double*,double*,double*,double*); + + /* coordinate conversion for magnetic field vector */ + OutputConvFunc _ConvOutput; + void _BMag2SysIII(int,double*,double*,double*, + double*,double*,double*,double*, + double*,double*,double*, + double*,double*,double*); + void _BMag2PolSysIII(int,double*,double*,double*, + double*,double*,double*,double*, + double*,double*,double*, + double*,double*,double*); + + /* Functions to update function pointers */ + void _SetIOFunctions(); + void _SetModelFunctions(); + ModelFunc _Model; + + /* Azimuthal field */ + AzimFunc _AzimuthalField; + void _BphiConnerney(int,double*,double*,double*,double*); + void _BphiConnerney(double,double,double,double*); + void _BphiLMIC(double,double,double,double*); + Approx _LargeRho; + Approx _SmallRho; + /* analytic equations */ + void _Analytic(double,double,double,double*,double*,double*); + void _AnalyticSmooth(double,double,double,double*,double*,double*); + void _SolveAnalytic(int,double*,double*,double,double*,double*); + void _AnalyticInner(double,double,double*,double*); + void _AnalyticOuter(double,double,double*,double*); + void _AnalyticInnerSmooth(double,double,double*,double*); + void _AnalyticOuterSmooth(double,double,double*,double*); + void _LargeRhoConnerney(double,double,double,double,double,double*,double*); + void _SmallRhoConnerney(double,double,double,double,double,double*,double*); + void _LargeRhoEdwards(double,double,double,double,double,double*,double*); + void _LargeRhoEdwardsSmooth(double,double,double,double,double,double*,double*); + void _SmallRhoEdwards(double,double,double,double,double,double*,double*); + + /* integral-related functions */ + void _Integral(double,double,double,double*,double*,double*); + void _IntegralInner(double, double, double, double*, double*); + void _InitIntegrals(); + void _RecalcIntegrals(); + void _DeleteIntegrals(); + void _IntegralChecks(int,double*,int*,int[]); + void _IntegralCheck(double,int*); + void _SolveIntegral(int,double*,double*,double*,double*,double*); + void _IntegrateEq14(int,double,double,double,double*); + void _IntegrateEq15(int,double,double,double*); + void _IntegrateEq17(int,double,double,double*); + void _IntegrateEq18(int,double,double,double*); + + /* hybrid */ + void _Hybrid(double,double,double,double*,double*,double*); +}; + +/* we want to initialize the model objects with its parameters */ +extern Con2020 con2020; + + +double polyeval(double x, double *c, int d); + +double pol1eval(double x, double *c, int d); + +/*********************************************************************** + * NAME : smoothd(z,dz,d) + * + * DESCRIPTION : Smooth fucntion for crossing the current sheet + * (replaces the last bit of equation 12 in Edwards et al 2000). + * + * INPUTS : + * double z z-coordinate in dipole coordinate system (Rj) + * double dz Scale of the transition to use (Rj) + * double d Half thickness of the current sheet. + * + * RETURNS : + * double out Smoothed function across the current sheet. + * + * ********************************************************************/ +double smoothd(double z, double dz, double d); + + +double trap(int n, double *x, double *y); +double trapc(int n, double dx, double *y); + +template T sgn(T x) { + return (x > 0) - (x < 0); +} + + +#endif diff --git a/libcon2020/modif/src/con2020.cc b/libcon2020/modif/src/con2020.cc new file mode 100644 index 0000000..0c1c12e --- /dev/null +++ b/libcon2020/modif/src/con2020.cc @@ -0,0 +1,1381 @@ +#include "con2020.h" + +Con2020::Con2020() { + /* set all model parameters to their default values */ + mui_ = 139.6; + irho_ = 16.7; + r0_ = 7.8; + r0sq_ = r0_*r0_; + r1_ = 51.4; + r1sq_ = r1_*r1_; + d_ = 3.6; + xt_ = 9.3; + xp_ = 155.8; + strcpy(eqtype_,"hybrid"); + Edwards_ = true; + ErrChk_ = true; + CartIn_ = true; + CartOut_ = true; + deltarho_ = 1.0; + deltaz_ = 0.1; + smooth_ = false; + wO_open_ = 0.1; + wO_om_ = 0.35; + thetamm_ = 16.1*deg2rad; + dthetamm_ = 0.5*deg2rad; + thetaoc_ = 10.716*deg2rad; + dthetaoc_ = 0.125*deg2rad; + g_ = 417659.3836476442; + strcpy(azfunc_,"connerney"); + + /* some other values which will only need calculating once */ + discshift_ = (xp_-180.0)*deg2rad; + disctilt_ = xt_*deg2rad; + cosxp_ = cos(discshift_); + sinxp_ = sin(discshift_); + cosxt_ = cos(disctilt_); + sinxt_ = sin(disctilt_); + + // BRE - 10/10/2023 - Do not init integrals in constructor to improve library load + /* initialize some things used for integration */ + init_ = false; + //_InitIntegrals(); + //------------------------- + + /* set appropriate coord conversion functions*/ + _SetIOFunctions(); + _SetModelFunctions(); + +} + +Con2020::Con2020(double mui, double irho, double r0, double r1, + double d, double xt, double xp, const char *eqtype, + bool Edwards, bool ErrChk, bool CartIn, bool CartOut) { + /* set non-boolean model parameters to their default values */ + mui_ = 139.6; + irho_ = 16.7; + r0_ = 7.8; + r0sq_ = r0_*r0_; + r1_ = 51.4; + r1sq_ = r1_*r1_; + d_ = 3.6; + xt_ = 9.3; + xp_ = 155.8; + strcpy(eqtype_,"hybrid"); + Edwards_ = Edwards; + ErrChk_ = ErrChk; + CartIn_ = CartIn; + CartOut_ = CartOut; + deltarho_ = 1.0; + deltaz_ = 0.01; + smooth_ = false; + wO_open_ = 0.1; + wO_om_ = 0.35; + thetamm_ = 16.1*deg2rad; + dthetamm_ = 0.5*deg2rad; + thetaoc_ = 10.716*deg2rad; + dthetaoc_ = 0.125*deg2rad; + g_ = 417659.3836476442; + strcpy(azfunc_,"connerney"); + + /* apply custom values if they are valid */ + SetAzCurrentParameter(mui); + SetRadCurrentParameter(irho); + SetR0(r0); + SetR1(r1); + SetCSHalfThickness(d); + SetCSTilt(xt); + SetCSTiltAzimuth(xp); + + + /* some other values which will only need calculating once */ + discshift_ = (xp_-180.0)*deg2rad; + disctilt_ = xt_*deg2rad; + + // BRE - 10/10/2023 - Do not init integrals in constructor to improve library load + /* initialize some things used for integration */ + init_ = false; + //_InitIntegrals(); + //------------------------- + + /* set appropriate coord conversion functions*/ + _SetIOFunctions(); + _SetModelFunctions(); + +} +Con2020::~Con2020() { + // BRE - 10/10/2023 - Delete integrals only if init + if (init_) + _DeleteIntegrals(); + //------------------------- +} + + +// BRE - 10/10/2023 - Add method to init integrals +void Con2020::Init() { + + if (init_) + return; + + /* initialize some things used for integration */ + _InitIntegrals(); + + init_ = true; +} +//------------------------- + +void Con2020::_SetIOFunctions() { + + if (CartIn_) { + _ConvInput = &Con2020::_SysIII2Mag; + } else { + _ConvInput = &Con2020::_PolSysIII2Mag; + } + if (CartOut_) { + _ConvOutput = &Con2020::_BMag2SysIII; + } else { + _ConvOutput = &Con2020::_BMag2PolSysIII; + } +} + +void Con2020::_SetModelFunctions() { + + /* firstly set whether we are using the Edwards or Connerney + * analytical functions */ + if (Edwards_) { + if (smooth_) { + _LargeRho = &Con2020::_LargeRhoEdwardsSmooth; + } else { + _LargeRho = &Con2020::_LargeRhoEdwards; + } + _SmallRho = &Con2020::_SmallRhoEdwards; + } else { + _LargeRho = &Con2020::_LargeRhoConnerney; + _SmallRho = &Con2020::_SmallRhoConnerney; + } + + /* now we need to set which model functions we will use + * (analytic, intergral or hybrid) */ + if (strcmp(eqtype_,"analytic") == 0) { + if (smooth_) { + _Model = &Con2020::_AnalyticSmooth; + } else { + _Model = &Con2020::_Analytic; + } + } else if (strcmp(eqtype_,"integral") == 0) { + _Model = &Con2020::_Integral; + } else if (strcmp(eqtype_,"hybrid") == 0) { + _Model = &Con2020::_Hybrid; + } else { + printf("What's going on here then?\n"); + } + + /* set the azimuthal function */ + if (strcmp(azfunc_,"connerney") == 0) { + _AzimuthalField = &Con2020::_BphiConnerney; + } else if (strcmp(azfunc_,"lmic") == 0) { + _AzimuthalField = &Con2020::_BphiLMIC; + } else { + printf("Azimuthal function %s not recognised\n",azfunc_); + } + +} + +void Con2020::_SysIII2Mag(int n, double *x0, double *y0, double *z0, + double *x1, double *y1, double *z1, double *rho, double *absz, + double *cost, double *sint, double *cosp, double *sinp) { + + + /* some temporary variables which get used more than once */ + double xt, theta, phi, r, rho0, rho0_sq; + + + int i; + for (i=0;i 0.0) { + Bphi[i] = -Bphi[i]; + } + + } +} + +void Con2020::_BphiConnerney(double rho, double absz, double z, double *Bphi) { + + Bphi[0] = 2.7975*irho_/rho; + + if (absz < d_) { + Bphi[0] = Bphi[0]*absz/d_; + } + + if (z > 0.0) { + Bphi[0] = -Bphi[0]; + } +} + +void Con2020::_BphiLMIC(double rho, double absz, double z, double *Bphi) { + + double r = sqrt(rho*rho + z*z); + double theta = asin(rho/r); + + Bphi[0] = BphiLMIC(r,theta,g_,r0_,r1_,mui_,d_,deltarho_,deltaz_, + wO_open_,wO_om_,thetamm_,dthetamm_,thetaoc_,dthetaoc_); + +} + +void Con2020::Field(double p0, double p1, double p2, + double *B0, double *B1, double *B2) { + + /* create a bunch of empty variables */ + double x, y, z, rho, absz; + double cost, sint, cosp, sinp; + double Brho, Bphi, Bz; + + /* convert the input coordinates */ + (this->*_ConvInput)(1,&p0,&p1,&p2,&x,&y,&z,&rho,&absz,&cost,&sint,&cosp,&sinp); + + /* calculate the mode field */ + (this->*_Model)(rho,absz,z,&Brho,&Bphi,&Bz); + + /*convert the output field */ + (this->*_ConvOutput)(1,&x,&y,&rho,&cost,&sint,&cosp,&sinp,&Brho,&Bphi,&Bz,B0,B1,B2); +} + + +void Con2020::Field(int n, double *p0, double *p1, double *p2, + double *B0, double *B1, double *B2) { + + int i; + + /* allocate arrays to store model coordinates and fields */ + double *x = new double[n]; + double *y = new double[n]; + double *z = new double[n]; + double *absz = new double[n]; + double *sint = new double[n]; + double *sinp = new double[n]; + double *cost = new double[n]; + double *cosp = new double[n]; + double *rho = new double[n]; + double *Brho = new double[n]; + double *Bphi = new double[n]; + double *Bz = new double[n]; + + /* convert input coordinates to the coordinate system used by the model */ + (this->*_ConvInput)(n,p0,p1,p2,x,y,z,rho,absz,cost,sint,cosp,sinp); + + /* call the model */ + for (i=0;i*_Model)(rho[i],absz[i],z[i],&Brho[i],&Bphi[i],&Bz[i]); + } + + /* convert B field back to appropriate coordinate system */ + (this->*_ConvOutput)(n,x,y,rho,cost,sint,cosp,sinp,Brho,Bphi,Bz,B0,B1,B2); + + /* delete temporary variables */ + delete[] x; + delete[] y; + delete[] z; + delete[] absz; + delete[] rho; + delete[] sint; + delete[] sinp; + delete[] cost; + delete[] cosp; + delete[] Brho; + delete[] Bphi; + delete[] Bz; + + +} + + + + + + +void Con2020::AnalyticFieldSmooth(double a, + double rho, double z, + double *Brho, double *Bz) { + + /* define a few required variables */ + double zpd = z + d_; + double zmd = z - d_; + + /* define some temporary variables */ + double Brho0, Brho1, Bz0, Bz1; + double tanhrho, C0, C1; + double asq = a*a; + + /* calculate small and large rho approximations + * NOTE: Add smoothed functions for small and large rho */ + (this->*_LargeRho)(rho,z,zmd,zpd,asq,&Brho1,&Bz1); + (this->*_SmallRho)(rho,z,zmd,zpd,asq,&Brho0,&Bz0); + + /* calculate the tanh smoothing parameters */ + tanhrho = tanh((rho-a)/deltarho_); + C0 = (1-tanhrho)/2.0; + C1 = (1+tanhrho)/2.0; + + /* splice together as suggested by Stan */ + *Brho = Brho0*C0 + Brho1*C1; + *Bz = Bz0*C0 + Bz1*C1; + + +} + +void Con2020::AnalyticField(double a, + double rho, double z, + double *Brho, double *Bz) { + + /* define a few required variables */ + double zpd = z + d_; + double zmd = z - d_; + + /* define some temporary variables */ + double asq = a*a; + + /* calculate small and large rho approximations */ + if (rho >= a) { + (this->*_LargeRho)(rho,z,zmd,zpd,asq,Brho,Bz); + } else { + (this->*_SmallRho)(rho,z,zmd,zpd,asq,Brho,Bz); + } + + + +} + +void Con2020::_Analytic(double rho, double absz, double z, + double *Brho, double *Bphi, double *Bz) { + + /* calculate the inner contribution to Brho and Bz */ + _AnalyticInner(rho,z,Brho,Bz); + + /* also the azimuthal field */ + (this->*_AzimuthalField)(rho,absz,z,Bphi); + + /* we need to calculate the outer edge contribution */ + double oBrho, oBz; + _AnalyticOuter(rho,z,&oBrho,&oBz); + + /*subtract it */ + Brho[0] -= oBrho; + Bz[0] -= oBz; + +} + +void Con2020::_AnalyticSmooth(double rho, double absz, double z, + double *Brho, double *Bphi, double *Bz) { + + /* calculate the inner contribution to Brho and Bz */ + _AnalyticInnerSmooth(rho,z,Brho,Bz); + + /* also the azimuthal field */ + (this->*_AzimuthalField)(rho,absz,z,Bphi); + + /* we need to calculate the outer edge contribution */ + double oBrho, oBz; + _AnalyticOuterSmooth(rho,z,&oBrho,&oBz); + + /*subtract it */ + Brho[0] -= oBrho; + Bz[0] -= oBz; + +} + +void Con2020::_AnalyticInner( double rho, double z, + double *Brho, double *Bz) { + + /* define a few required variables */ + double zpd = z + d_; + double zmd = z - d_; + + if (rho >= r0_) { + (this->*_LargeRho)(rho,z,zmd,zpd,r0sq_,Brho,Bz); + } else { + (this->*_SmallRho)(rho,z,zmd,zpd,r0sq_,Brho,Bz); + } + +} + +void Con2020::_AnalyticInnerSmooth( double rho, double z, + double *Brho, double *Bz) { + + /* define a few required variables */ + double zpd = z + d_; + double zmd = z - d_; + + /* define some temporary variables */ + double Brho0, Brho1, Bz0, Bz1; + double tanhrho, C0, C1; + + /* calculate small and large rho approximations + * NOTE: Add smoothed functions for small and large rho */ + (this->*_LargeRho)(rho,z,zmd,zpd,r0sq_,&Brho1,&Bz1); + (this->*_SmallRho)(rho,z,zmd,zpd,r0sq_,&Brho0,&Bz0); + + /* calculate the tanh smoothing parameters */ + tanhrho = tanh((rho-r0_)/deltarho_); + C0 = (1-tanhrho)/2.0; + C1 = (1+tanhrho)/2.0; + + /* splice together as suggested by Stan */ + *Brho = Brho0*C0 + Brho1*C1; + *Bz = Bz0*C0 + Bz1*C1; + + +} + +void Con2020::_AnalyticOuter( double rho, double z, + double *Brho, double *Bz) { + + /* define a few required variables */ + double zpd = z + d_; + double zmd = z - d_; + + if (rho >= r1_) { + (this->*_LargeRho)(rho,z,zmd,zpd,r1sq_,Brho,Bz); + } else { + (this->*_SmallRho)(rho,z,zmd,zpd,r1sq_,Brho,Bz); + } + +} + + +void Con2020::_AnalyticOuterSmooth( double rho, double z, + double *Brho, double *Bz) { + + /* define a few required variables */ + double zpd = z + d_; + double zmd = z - d_; + + /* define some temporary variables */ + double Brho0, Brho1, Bz0, Bz1; + double tanhrho, C0, C1; + + /* calculate small and large rho approximations + * NOTE: Add smoothed functions for small and large rho */ + (this->*_LargeRho)(rho,z,zmd,zpd,r1sq_,&Brho1,&Bz1); + (this->*_SmallRho)(rho,z,zmd,zpd,r1sq_,&Brho0,&Bz0); + + /* calculate the tanh smoothing parameters */ + tanhrho = tanh((rho-r1_)/deltarho_); + C0 = (1-tanhrho)/2.0; + C1 = (1+tanhrho)/2.0; + + /* splice together as suggested by Stan */ + *Brho = Brho0*C0 + Brho1*C1; + *Bz = Bz0*C0 + Bz1*C1; + + +} + + +void Con2020::_LargeRhoConnerney(double rho, double z, double zmd, + double zpd, double a2, double *Brho, double *Bz) { + + /* some common variables */ + double zmd2 = zmd*zmd; + double zpd2 = zpd*zpd; + double rho2 = rho*rho; + double f1 = sqrt(zmd2 + rho2); + double f2 = sqrt(zpd2 + rho2); + double f1cubed = f1*f1*f1; + double f2cubed = f2*f2*f2; + + /* Equation A7 */ + double termr0 = (1.0/rho)*(f1 -f2 + 2*clip(z,-d_,d_)); + double termr1 = (a2*rho/4.0)*(1/f1cubed - 1/f2cubed); + Brho[0] = mui_*(termr0 - termr1); + + /* Equation A8 */ + double termz0 = 2*d_/sqrt(z*z + rho*rho); + double termz1 = (a2/4.0)*(zmd/f1cubed - zpd/f2cubed); + Bz[0] = mui_*(termz0 - termz1); + +} + +void Con2020::_SmallRhoConnerney(double rho, double z, double zmd, double zpd, + double a2, double *Brho, double *Bz) { + + double zmd2 = zmd*zmd; + double zpd2 = zpd*zpd; + double f1 = sqrt(zmd2 + a2); + double f2 = sqrt(zpd2 + a2); + double f1cubed = f1*f1*f1; + double f2cubed = f2*f2*f2; + + /* Equation A1 */ + Brho[0] = mui_*(rho/2.0)*(1.0/f1 - 1.0/f2); + + /* Equation A2 */ + Bz[0] = mui_*(2*d_*(1/sqrt(z*z + a2)) - ((rho*rho)/4)*((zmd/f1cubed) - (zpd/f2cubed))); + +} + + +void Con2020::_LargeRhoEdwards(double rho, double z, double zmd, + double zpd, double a2,double *Brho, double *Bz) { + + + /* some common variables */ + double zmd2 = zmd*zmd; + double zpd2 = zpd*zpd; + double rho2 = rho*rho; + double f1 = sqrt(zmd2 + rho2); + double f2 = sqrt(zpd2 + rho2); + double f1cubed = f1*f1*f1; + double f2cubed = f2*f2*f2; + + /* equation 13a */ + double terma0 = (1.0/rho)*(f1 - f2); + double terma1 = (rho*a2/4)*(1.0/f2cubed - 1.0/f1cubed); + double terma2 = (2.0/rho)*clip(z,-d_,d_); + Brho[0] = mui_*(terma0 + terma1 + terma2); + + /* equation 13b */ + double termb0 = log((zpd + f2)/(zmd + f1)); + double termb1 = (a2/4.0)*(zpd/f2cubed - zmd/f1cubed); + Bz[0] = mui_*(termb0 + termb1); +} + +void Con2020::_LargeRhoEdwardsSmooth(double rho, double z, double zmd, + double zpd, double a2,double *Brho, double *Bz) { + + /* some common variables */ + double zmd2 = zmd*zmd; + double zpd2 = zpd*zpd; + double rho2 = rho*rho; + double f1 = sqrt(zmd2 + rho2); + double f2 = sqrt(zpd2 + rho2); + double f1cubed = f1*f1*f1; + double f2cubed = f2*f2*f2; + + /* equation 13a */ + double terma0 = (1.0/rho)*(f1 - f2); + double terma1 = (rho*a2/4)*(1.0/f2cubed - 1.0/f1cubed); + double terma2 = (2.0/rho)*smoothd(z,deltaz_,d_); + Brho[0] = mui_*(terma0 + terma1 + terma2); + + /* equation 13b */ + double termb0 = log((zpd + f2)/(zmd + f1)); + double termb1 = (a2/4.0)*(zpd/f2cubed - zmd/f1cubed); + Bz[0] = mui_*(termb0 + termb1); +} + + +void Con2020::_SmallRhoEdwards(double rho, double z, double zmd, double zpd, + double a2, double *Brho, double *Bz) { + + double zmd2 = zmd*zmd; + double zpd2 = zpd*zpd; + double f1 = sqrt(zmd2 + a2); + double f2 = sqrt(zpd2 + a2); + double f1cubed = f1*f1*f1; + double f2cubed = f2*f2*f2; + + /* calculate some of the common terms from equations 9a and 9b */ + double rhoov2 = rho/2.0; + double rho2ov4 = rhoov2*rhoov2; + double rho3ov16 = rho2ov4*rhoov2/2.0; + + /* equation 9a */ + double f3a = f1*f1; + double f4a = f2*f2; + double f3 = (a2 - 2*zmd2)/(f3a*f3a*f1); + double f4 = (a2 - 2*zpd2)/(f4a*f4a*f2); + + double terma0 = rhoov2*(1.0/f1 - 1.0/f2); + double terma1 = rho3ov16*(f3 - f4); + + Brho[0] = mui_*(terma0 + terma1); + + /* equation 9b */ + double termb0 = log((zpd + f2)/(zmd + f1)); + double termb1 = rho2ov4*(zpd/f2cubed - zmd/f1cubed); + Bz[0] = mui_*(termb0 + termb1); +} + + +void Con2020::_InitIntegrals() { + + /* the lambda max values for brho and bz */ + double lmx_brho, lmx_bz; + rlmx_array_[0] = 4.0; + rlmx_array_[1] = 4.0; + rlmx_array_[2] = 40.0; + rlmx_array_[3] = 40.0; + rlmx_array_[4] = 100.0; + rlmx_array_[5] = 100.0; + zlmx_array_[0] = 100.0; + zlmx_array_[1] = 20.0; + zlmx_array_[2] = 100.0; + zlmx_array_[3] = 20.0; + zlmx_array_[4] = 100.0; + zlmx_array_[5] = 20.0; + + + /* initialize the bessel functions which do not change */ + rnbes_ = new int[6]; + znbes_ = new int[6]; + rlambda_ = new double*[6]; + zlambda_ = new double*[6]; + rj0_lambda_r0_ = new double*[6]; + zj0_lambda_r0_ = new double*[6]; + + /*These functions do change and will need reassigning + * for each run of the integral function*/ + rj1_lambda_rho_ = new double*[6]; + zj0_lambda_rho_ = new double*[6]; + + /*these will store bits of the Connerney et al 1981 equations which + * don't change*/ + Eq14_ = new double*[6]; + Eq15_ = new double*[6]; + Eq17_ = new double*[6]; + Eq18_ = new double*[6]; + ExpLambdaD_ = new double*[6]; + + + /* initialize the second dimensions */ + int zcase, i; + double ld; + for (zcase=0;zcase<6;zcase++) { + /* calculate the length of the arrays for each case */ + rnbes_[zcase] = (int) (rlmx_array_[zcase]/dlambda_brho_) - 1; + znbes_[zcase] = (int) (zlmx_array_[zcase]/dlambda_bz_) - 1; + + /* allocate the second dimension */ + rlambda_[zcase] = new double[rnbes_[zcase]]; + zlambda_[zcase] = new double[znbes_[zcase]]; + rj0_lambda_r0_[zcase] = new double[rnbes_[zcase]]; + zj0_lambda_r0_[zcase] = new double[znbes_[zcase]]; + rj1_lambda_rho_[zcase] = new double[rnbes_[zcase]]; + zj0_lambda_rho_[zcase] = new double[znbes_[zcase]]; + Eq14_[zcase] = new double[rnbes_[zcase]]; + Eq15_[zcase] = new double[znbes_[zcase]]; + Eq17_[zcase] = new double[rnbes_[zcase]]; + Eq18_[zcase] = new double[znbes_[zcase]]; + ExpLambdaD_[zcase] = new double[znbes_[zcase]]; + } + + _RecalcIntegrals(); +} + +void Con2020::_RecalcIntegrals(){ + + int zcase, i; + double ld; + for (zcase=0;zcase<6;zcase++) { + /* initialize lambda for z and rho cases */ + for (i=0;i*_AzimuthalField)(rho,absz,z,Bphi); + + /* we need to calculate the outer edge contribution */ + double oBrho, oBz; + _AnalyticOuter(rho,z,&oBrho,&oBz); + + /*subtract it */ + Brho[0] -= oBrho; + Bz[0] -= oBz; + +} + + +void Con2020::_IntegralChecks(int n, double *absz, int *chind, int ncase[]) { + + int i; + double check1; + bool check2; + + for (i=0;i<6;i++) { + ncase[i] = 0; + } + + for (i=0;i= 0.7) { + chind[i] = 1; + } else if (check1 < 0.1) { + chind[i] = 5; + } else { + chind[i] = 3; + } + chind[i] -= ((int) check2); + ncase[chind[i]]++; + } + +} + +void Con2020::_IntegralCheck(double absz, int *chind) { + + int i; + double check1; + bool check2; + + + check1 = fabs(absz - d_); + check2 = (absz < d_*1.1); + + if (check1 >= 0.7) { + chind[0] = 1; + } else if (check1 < 0.1) { + chind[0] = 5; + } else { + chind[0] = 3; + } + chind[0] -= ((int) check2); + + +} + +void Con2020::_SolveIntegral(int n, double *rho, double *z, + double *absz, double *Brho, double *Bz) { + + int i, zcase, ncase[6]; + + + /* calculate the "checks" */ + int *chind = new int[n]; + _IntegralChecks(n,absz,chind,ncase); + + /* loop through each position */ + double br, bz; + for (i=0;i d_) { + /* in this case we want to use equations 14 and 15 outside + * of the finite current sheet */ + _IntegrateEq14(chind[i],rho[i],z[i],absz[i],&Brho[i]); + _IntegrateEq15(chind[i],rho[i],absz[i],&Bz[i]); + + } else { + /* here we are inside the current sheet and should integrate + * equations 17 and 18 */ + _IntegrateEq17(chind[i],rho[i],z[i],&Brho[i]); + _IntegrateEq18(chind[i],rho[i],z[i],&Bz[i]); + } + } + + delete[] chind; + +} + +void Con2020::_IntegralInner( double rho, double absz, double z, + double *Brho, double *Bz) { + //printf("\n_IntegralInner\n"); + int chind; + /* check which set of integral parameters we need to use*/ + _IntegralCheck(absz,&chind); + + if (absz > d_) { + /* in this case we want to use equations 14 and 15 outside + * of the finite current sheet */ + _IntegrateEq14(chind,rho,z,absz,Brho); + _IntegrateEq15(chind,rho,absz,Bz); + } else { + /* here we are inside the current sheet and should integrate + * equations 17 and 18 */ + _IntegrateEq17(chind,rho,z,Brho); + _IntegrateEq18(chind,rho,z,Bz); + } + +} + +void Con2020::_IntegrateEq14(int zcase, double rho, double z, double absz, double *Brho) { + + /* create an array to integrate over */ + int n = rnbes_[zcase]; + double *func = new double[n]; + double *j1lr = new double[n]; + + /* calculate the other bessel function */ + j1(n,&rlambda_[zcase][0],rho,j1lr); + + + /* calculate the function */ + int i; + for (i=0;i*_AzimuthalField)(rho,absz,z,Bphi); + + /* we need to calculate the outer edge contribution */ + double oBrho, oBz; + _AnalyticOuter(rho,z,&oBrho,&oBz); + + /*subtract it */ + Brho[0] -= oBrho; + Bz[0] -= oBz; + +} + + +void Con2020::SetEdwardsEqs(bool Edwards) { + Edwards_ = Edwards; + + /* reset function pointers */ + _SetModelFunctions(); +} + +bool Con2020::GetEdwardsEqs() { + return Edwards_; +} + +void Con2020::SetEqType(const char *eqtype) { + + if ( (strcmp(eqtype,"analytic") == 0) || + (strcmp(eqtype,"integral") == 0) || + (strcmp(eqtype,"hybrid") == 0)) { + /* this is a valid string - update it */ + strcpy(eqtype_,eqtype); + + _SetModelFunctions(); + } else { + printf("eqtype '%s' not recognised - ignoring\n",eqtype); + } + +} + +void Con2020::SetSmooth(bool smooth) { + smooth_ = smooth; + _SetModelFunctions(); +} + +bool Con2020::GetSmooth() { + + return smooth_; +} + + +void Con2020::GetEqType(char *eqtype) { + strcpy(eqtype,eqtype_); +} + +void Con2020::SetAzCurrentParameter(double mui) { + + if (std::isfinite(mui)) { + /* good value (hopefully) */ + mui_ = mui; + } else { + printf("Non-finite value - ignoring\n"); + } +} + +double Con2020::GetAzCurrentParameter() { + return mui_; +} + +void Con2020::SetRadCurrentParameter(double irho) { + if (std::isfinite(irho)) { + /* good value (hopefully) */ + irho_ = irho; + } else { + printf("Non-finite value - ignoring\n"); + } +} + +double Con2020::GetRadCurrentParameter() { + return irho_; +} + +void Con2020::SetR0(double r0) { + if (std::isfinite(r0) && (r0 >= 0.0)) { + /* good value (hopefully) */ + r0_ = r0; + r0sq_ = r0_*r0_; + } else if (!std::isfinite(r0)) { + printf("Non-finite value - ignoring\n"); + } else { + printf("r0 must have a positive value\n"); + } +} + +double Con2020::GetR0() { + return r0_; +} + +void Con2020::SetR1(double r1) { + if (std::isfinite(r1) && (r1 >= 0.0)) { + /* good value (hopefully) */ + r1_ = r1; + r1sq_ = r1_*r1_; + } else if (!std::isfinite(r1)) { + printf("Non-finite value - ignoring\n"); + } else { + printf("r1 must have a positive value\n"); + } +} + +double Con2020::GetR1() { + return r1_; +} + +void Con2020::SetCSHalfThickness(double d) { + + if (std::isfinite(d) && (d >= 0.0)) { + /* good value (hopefully) */ + d_ = d; + } else if (!std::isfinite(d)) { + printf("Non-finite value - ignoring\n"); + } else { + printf("d must have a positive value\n"); + } +} + +double Con2020::GetCSHalfThickness() { + return d_; +} + +void Con2020::SetCSTilt(double xt) { + + if (std::isfinite(xt)) { + /* good value (hopefully) */ + xt_ = xt; + disctilt_ = xt_*deg2rad; + cosxt_ = cos(disctilt_); + sinxt_ = sin(disctilt_); + } else { + printf("Non-finite value - ignoring\n"); + } +} + +double Con2020::GetCSTilt() { + return xt_; +} + +void Con2020::SetCSTiltAzimuth(double xp) { + if (std::isfinite(xp)) { + /* good value (hopefully) */ + xp_ = xp; + discshift_ = (xp_ - 180.0)*deg2rad; + cosxp_ = cos(discshift_); + sinxp_ = sin(discshift_); + } else { + printf("Non-finite value - ignoring\n"); + } +} + +double Con2020::GetCSTiltAzimuth() { + return xp_; +} + +void Con2020::SetErrCheck(bool ErrChk) { + ErrChk_ = ErrChk; +} + +bool Con2020::GetErrCheck() { + return ErrChk_; +} + +void Con2020::SetCartIn(bool CartIn) { + CartIn_ = CartIn; + _SetIOFunctions(); +} + +bool Con2020::GetCartIn() { + return CartIn_; +} + +void Con2020::SetCartOut(bool CartOut) { + CartOut_ = CartOut; + _SetIOFunctions(); +} + +bool Con2020::GetCartOut() { + return CartOut_; +} + +void Con2020::SetDeltaRho(double DeltaRho) { + deltarho_ = DeltaRho; +} + +double Con2020::GetDeltaRho() { + return deltarho_; +} + +void Con2020::SetDeltaZ(double DeltaZ) { + deltaz_ = DeltaZ; +} + +double Con2020::GetDeltaZ() { + return deltaz_; +} + +void Con2020::SetOmegaOpen(double OmegaOpen) { + wO_open_ = OmegaOpen; +} + +double Con2020::GetOmegaOpen() { + return wO_open_; +} + +void Con2020::SetOmegaOM(double OmegaOM) { + wO_om_ = OmegaOM; +} + +double Con2020::GetOmegaOM() { + return wO_om_; +} + +void Con2020::SetThetaMM(double ThetaMM) { + /* we need this to be in radians, but will accept it in degrees */ + thetamm_ = ThetaMM*deg2rad; +} + +double Con2020::GetThetaMM() { + /* convert back to degrees*/ + return thetamm_*rad2deg; +} + +void Con2020::SetdThetaMM(double dThetaMM) { + /* we need this to be in radians, but will accept it in degrees */ + dthetamm_ = dThetaMM*deg2rad; +} + +double Con2020::GetdThetaMM() { + /* convert back to degrees*/ + return dthetamm_*rad2deg; +} + +void Con2020::SetThetaOC(double ThetaOC) { + /* we need this to be in radians, but will accept it in degrees */ + thetaoc_ = ThetaOC*deg2rad; +} + +double Con2020::GetThetaOC() { + /* convert back to degrees*/ + return thetaoc_*rad2deg; +} + +void Con2020::SetdThetaOC(double dThetaOC) { + /* we need this to be in radians, but will accept it in degrees */ + dthetaoc_ = dThetaOC*deg2rad; +} + +double Con2020::GetdThetaOC() { + /* convert back to degrees*/ + return dthetaoc_*rad2deg; +} + +void Con2020::SetG(double g) { + g_ = g; +} + +double Con2020::GetG() { + return g_; +} + +void Con2020::SetAzimuthalFunc(const char* azfunc) { + if (strcmp(azfunc,"connerney") || strcmp(azfunc,"lmic")) { + strcpy(azfunc_,azfunc); + _SetModelFunctions(); + } else { + printf("Azimuthal function %s not recognised\n",azfunc); + } +} + +void Con2020::GetAzimuthalFunc(char* azfunc) { + strcpy(azfunc,azfunc_); +} diff --git a/libcon2020/modif/src/con2020.h b/libcon2020/modif/src/con2020.h new file mode 100644 index 0000000..92bde0e --- /dev/null +++ b/libcon2020/modif/src/con2020.h @@ -0,0 +1,227 @@ +#ifndef __CON2020_H__ +#define __CON2020_H__ +#include +#include +#include +#define _USE_MATH_DEFINES +#include +#include "bessel.h" +#include "sgn.h" +#include "clip.h" +#include "smoothd.h" +#include "trap.h" +#include "lmic.h" +#define deg2rad M_PI/180.0 +#define rad2deg 180.0/M_PI + + +/* function pointer for input conversion */ +class Con2020; /*this is needed for the pointer below */ +typedef void (Con2020::*InputConvFunc)(int,double*,double*,double*, + double*,double*,double*,double*,double*, + double*,double*,double*,double*); +/* Output conversion */ +typedef void (Con2020::*OutputConvFunc)(int,double*,double*,double*, + double*,double*,double*,double*, + double*,double*,double*, + double*,double*,double*); + +/* Model function */ +typedef void (Con2020::*ModelFunc)(double,double,double,double*,double*,double*); + +/* analytical approximation equations */ +typedef void (Con2020::*Approx)(double,double,double,double,double,double*,double*); + +/* azimuthal function */ +typedef void (Con2020::*AzimFunc)(double,double,double,double*); + +class Con2020 { + public: + /* constructors */ + Con2020(); + Con2020(double,double,double,double,double,double,double,const char*,bool,bool,bool,bool); + + /* destructor */ + ~Con2020(); + + // BRE - 10/10/2023 - Add method to init integrals + void Init(); + //---------- + + /* these functions will be used to set the equations used, if + * they need to be changed post-initialisation */ + void SetEdwardsEqs(bool); + void SetEqType(const char*); + void SetAzCurrentParameter(double); + void SetRadCurrentParameter(double); + void SetR0(double); + void SetR1(double); + void SetCSHalfThickness(double); + void SetCSTilt(double); + void SetCSTiltAzimuth(double); + void SetErrCheck(bool); + void SetCartIn(bool); + void SetCartOut(bool); + void SetSmooth(bool); + void SetDeltaRho(double); + void SetDeltaZ(double); + void SetOmegaOpen(double); + void SetOmegaOM(double); + void SetThetaMM(double); + void SetdThetaMM(double); + void SetThetaOC(double); + void SetdThetaOC(double); + void SetG(double); + void SetAzimuthalFunc(const char*); + + /* these mamber functions will be the "getter" version of the + * above setters */ + bool GetEdwardsEqs(); + void GetEqType(char*); + double GetAzCurrentParameter(); + double GetRadCurrentParameter(); + double GetR0(); + double GetR1(); + double GetCSHalfThickness(); + double GetCSTilt(); + double GetCSTiltAzimuth(); + bool GetErrCheck(); + bool GetCartIn(); + bool GetCartOut(); + bool GetSmooth(); + double GetDeltaRho(); + double GetDeltaZ(); + double GetOmegaOpen(); + double GetOmegaOM(); + double GetThetaMM(); + double GetdThetaMM(); + double GetThetaOC(); + double GetdThetaOC(); + double GetG(); + void GetAzimuthalFunc(char *); + + /* This function will be used to call the model, it is overloaded + * so that we have one for arrays, one for scalars */ + void Field(int,double*,double*,double*,double*,double*,double*); + void Field(double,double,double,double*,double*,double*); + + /* a function for testing purposes...*/ + void AnalyticField(double,double,double,double*,double*); + void AnalyticFieldSmooth(double,double,double,double*,double*); + + private: + + bool init_; + + /* model parameters */ + double mui_,irho_,r0_,r1_,d_,xt_,xp_,disctilt_,discshift_; + double r0sq_, r1sq_; + double cosxp_,sinxp_,cosxt_,sinxt_; + char eqtype_[9]; + bool Edwards_, ErrChk_; + bool CartIn_,CartOut_; + double deltaz_,deltarho_; + bool smooth_; + + /* LMIC parameters*/ + double wO_open_, wO_om_, thetamm_, dthetamm_, thetaoc_, dthetaoc_, g_; + char azfunc_[10]; + + /* Bessel function arrays - arrays prefixed with r and z are + * to be used for integrals which calcualte Brho and Bz, + * respectively */ + int *rnbes_; /* number of elements for each Bessel function (rho)*/ + int *znbes_; /* same as above for z */ + double **rlambda_;/* Lambda array to integrate over rho*/ + double **zlambda_;/* Lambda array to integrate over z*/ + double **rj0_lambda_r0_; /* j0(lambda*r0) */ + double **rj1_lambda_rho_;/* j1(lambda*rho) */ + double **zj0_lambda_r0_; /* j0(lambda*r0) */ + double **zj0_lambda_rho_;/* j0(lambda*rho) */ + + + /* arrays to multiply be stuff to be integrated */ + /* these arrays will store the parts of equations 14, 15, 17 + * and 18 of Connerny 1981 which only need to be calculated once*/ + double **Eq14_; /* j0(lambda*r0)*sinh(lamba*d)/lambda */ + double **Eq15_; /* j0(lambda*r0)*sinh(lamba*d)/lambda */ + double **Eq17_; /* j0(lambda*r0)*exp(-lamba*d)/lambda */ + double **Eq18_; /* j0(lambda*r0)/lambda */ + double **ExpLambdaD_; + + + /* integration step sizes */ + static constexpr double dlambda_ = 1e-4; + static constexpr double dlambda_brho_ = 1e-4; + static constexpr double dlambda_bz_ = 5e-5; + + /* Arrays containing maximum lambda values */ + double rlmx_array_[6]; + double zlmx_array_[6]; + + + /* coordinate conversions for positions */ + InputConvFunc _ConvInput; + void _SysIII2Mag(int,double*,double*,double*, + double*,double*,double*,double*,double*, + double*,double*,double*,double*); + void _PolSysIII2Mag(int,double*,double*,double*, + double*,double*,double*,double*,double*, + double*,double*,double*,double*); + + + /* coordinate conversion for magnetic field vector */ + OutputConvFunc _ConvOutput; + void _BMag2SysIII(int,double*,double*,double*, + double*,double*,double*,double*, + double*,double*,double*, + double*,double*,double*); + void _BMag2PolSysIII(int,double*,double*,double*, + double*,double*,double*,double*, + double*,double*,double*, + double*,double*,double*); + + /* Functions to update function pointers */ + void _SetIOFunctions(); + void _SetModelFunctions(); + ModelFunc _Model; + + /* Azimuthal field */ + AzimFunc _AzimuthalField; + void _BphiConnerney(int,double*,double*,double*,double*); + void _BphiConnerney(double,double,double,double*); + void _BphiLMIC(double,double,double,double*); + Approx _LargeRho; + Approx _SmallRho; + /* analytic equations */ + void _Analytic(double,double,double,double*,double*,double*); + void _AnalyticSmooth(double,double,double,double*,double*,double*); + void _SolveAnalytic(int,double*,double*,double,double*,double*); + void _AnalyticInner(double,double,double*,double*); + void _AnalyticOuter(double,double,double*,double*); + void _AnalyticInnerSmooth(double,double,double*,double*); + void _AnalyticOuterSmooth(double,double,double*,double*); + void _LargeRhoConnerney(double,double,double,double,double,double*,double*); + void _SmallRhoConnerney(double,double,double,double,double,double*,double*); + void _LargeRhoEdwards(double,double,double,double,double,double*,double*); + void _LargeRhoEdwardsSmooth(double,double,double,double,double,double*,double*); + void _SmallRhoEdwards(double,double,double,double,double,double*,double*); + + /* integral-related functions */ + void _Integral(double,double,double,double*,double*,double*); + void _IntegralInner(double, double, double, double*, double*); + void _InitIntegrals(); + void _RecalcIntegrals(); + void _DeleteIntegrals(); + void _IntegralChecks(int,double*,int*,int[]); + void _IntegralCheck(double,int*); + void _SolveIntegral(int,double*,double*,double*,double*,double*); + void _IntegrateEq14(int,double,double,double,double*); + void _IntegrateEq15(int,double,double,double*); + void _IntegrateEq17(int,double,double,double*); + void _IntegrateEq18(int,double,double,double*); + + /* hybrid */ + void _Hybrid(double,double,double,double*,double*,double*); +}; +#endif diff --git a/libcon2020/modif/src/libcon2020.cc b/libcon2020/modif/src/libcon2020.cc new file mode 100644 index 0000000..0ff9714 --- /dev/null +++ b/libcon2020/modif/src/libcon2020.cc @@ -0,0 +1,123 @@ +#include "libcon2020.h" + +// BRE - Due to this global variable, the Con2020 constructor is called on library load +// => the init integrals is a little bit long so the call of this method has been moved outside of the constructor +Con2020 con2020; + +// BRE - 10/10/2023 - Add method to init integrals +void Con2020Init() { + con2020.Init(); +} +//------------------------- + +void Con2020FieldArray(int n, double *p0, double *p1, double *p2, + double *B0, double *B1, double *B2) { + + /* could create a separate function for default model */ + con2020.Field(n,p0,p1,p2,B0,B1,B2); + +} + +void Con2020Field(double p0, double p1, double p2, + double *B0, double *B1, double *B2) { + + /* could create a separate function for default model */ + con2020.Field(p0,p1,p2,B0,B1,B2); + +} + +void GetCon2020Params(double *mui, double *irho, double *r0, double *r1, + double *d, double *xt, double *xp, char *eqtype, + bool *Edwards, bool *ErrChk, bool *CartIn, bool *CartOut, + bool *smooth, double *DeltaRho, double *DeltaZ, + double *g, char *azfunc, double *wO_open, double *wO_om, + double *thetamm, double *dthetamm, double *thetaoc, double *dthetaoc) { + + mui[0] = con2020.GetAzCurrentParameter(); + irho[0] = con2020.GetRadCurrentParameter(); + r0[0] = con2020.GetR0(); + r1[0] = con2020.GetR1(); + d[0] = con2020.GetCSHalfThickness(); + xt[0] = con2020.GetCSTilt(); + xp[0] = con2020.GetCSTiltAzimuth(); + Edwards[0] = con2020.GetEdwardsEqs(); + ErrChk[0] = con2020.GetErrCheck(); + CartIn[0] = con2020.GetCartIn(); + CartOut[0] = con2020.GetCartOut(); + con2020.GetEqType(eqtype); + smooth[0] = con2020.GetSmooth(); + DeltaRho[0] = con2020.GetDeltaRho(); + DeltaZ[0] = con2020.GetDeltaZ(); + + /* new LMIC parameters */ + g[0] = con2020.GetG(); + con2020.GetAzimuthalFunc(azfunc); + wO_open[0] = con2020.GetOmegaOpen(); + wO_om[0] = con2020.GetOmegaOM(); + thetamm[0] = con2020.GetThetaMM(); + dthetamm[0] = con2020.GetdThetaMM(); + thetaoc[0] = con2020.GetThetaOC(); + dthetaoc[0] = con2020.GetdThetaOC(); + + +} +void SetCon2020Params(double mui, double irho, double r0, double r1, + double d, double xt, double xp, const char *eqtype, + bool Edwards, bool ErrChk, bool CartIn, bool CartOut, + bool smooth, double DeltaRho, double DeltaZ, + double g, const char *azfunc, double wO_open, double wO_om, + double thetamm, double dthetamm, double thetaoc, double dthetaoc) { + + con2020.SetAzCurrentParameter(mui); + con2020.SetRadCurrentParameter(irho); + con2020.SetR0(r0); + con2020.SetR1(r1); + con2020.SetCSHalfThickness(d); + con2020.SetCSTilt(xt); + con2020.SetCSTiltAzimuth(xp); + con2020.SetEdwardsEqs(Edwards); + con2020.SetErrCheck(ErrChk); + con2020.SetCartIn(CartIn); + con2020.SetCartOut(CartOut); + con2020.SetEqType(eqtype); + con2020.SetSmooth(smooth); + con2020.SetDeltaRho(DeltaRho); + con2020.SetDeltaZ(DeltaZ); + + /*set LMIC parameters */ + con2020.SetG(g); + con2020.SetAzimuthalFunc(azfunc); + con2020.SetOmegaOpen(wO_open); + con2020.SetOmegaOM(wO_om); + con2020.SetThetaMM(thetamm); + con2020.SetdThetaMM(dthetamm); + con2020.SetThetaOC(thetaoc); + con2020.SetdThetaOC(dthetaoc); +} + +void Con2020AnalyticField( int n, double a, + double *rho, double *z, + double *Brho, double *Bz) { + + /* define a few required variables */ + int i; + + for (i=0;i +#include +#include "con2020.h" +#include + + + +/* we want to initialize the model objects with its parameters */ +extern Con2020 con2020; + +extern "C" { + // BRE - 10/10/2023 - Add method to init integrals + void Con2020Init(); + //---------- + + /* these wrappers can be used to get the magnetic field vectors */ + void Con2020FieldArray(int n, double *p0, double *p1, double *p2, + double *B0, double *B1, double *B2); + + void Con2020Field(double p0, double p1, double p2, + double *B0, double *B1, double *B2); + + + void GetCon2020Params(double *mui, double *irho, double *r0, double *r1, + double *d, double *xt, double *xp, char *eqtype, + bool *Edwards, bool *ErrChk, bool *CartIn, bool *CartOut, + bool *smooth, double *DeltaRho, double *DeltaZ, + double *g, char *azfunc, double *wO_open, double *wO_om, + double *thetamm, double *dthetamm, double *thetaoc, double *dthetaoc); + + + void SetCon2020Params(double mui, double irho, double r0, double r1, + double d, double xt, double xp, const char *eqtype, + bool Edwards, bool ErrChk, bool CartIn, bool CartOut, + bool smooth, double DeltaRho, double DeltaZ, + double g, const char *azfunc, double wO_open, double wO_om, + double thetamm, double dthetamm, double thetaoc, double dthetaoc); + + void Con2020AnalyticField( int n, double a, + double *rho, double *z, + double *Brho, double *Bz); + + void Con2020AnalyticFieldSmooth( int n, double a, + double *rho, double *z, + double *Brho, double *Bz); + +} +#endif -- libgit2 0.21.2