From d3f396ce1c0e3646b57ab7e2f6a8738c6fafb10f Mon Sep 17 00:00:00 2001 From: Darren Garnier Date: Fri, 30 Jun 2023 01:17:21 -0400 Subject: [PATCH] finished docstrings --- .flake8 | 5 +- pyproject.toml | 2 + src/pysol/filaments.py | 173 +++++++++++ src/pysol/inductance.py | 642 ++++++++++++++++++++++++---------------- tests/test_elliptics.py | 5 + tests/test_ldx.py | 111 +++++-- 6 files changed, 657 insertions(+), 281 deletions(-) create mode 100644 src/pysol/filaments.py diff --git a/.flake8 b/.flake8 index 1e747b4..0721e16 100644 --- a/.flake8 +++ b/.flake8 @@ -1,8 +1,9 @@ [flake8] select = B,B9,C,D,DAR,E,F,N,RST,W -ignore = E203,E501,RST201,RST203,RST301,W503 -max-line-length = 80 +ignore = E203,E501,RST201,RST203,RST301,W503,N802,N803,N806 +max-line-length = 88 max-complexity = 10 docstring-convention = google rst-roles = class,const,func,meth,mod,ref rst-directives = deprecated +exclude = src/pysol/elliptics.py diff --git a/pyproject.toml b/pyproject.toml index d2a338d..3b13090 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -35,3 +35,5 @@ build-backend = "poetry_dynamic_versioning.backend" extend-exclude = ''' ( src/elliptics.py ) ''' +[tool.isort] +profile = "black" diff --git a/src/pysol/filaments.py b/src/pysol/filaments.py new file mode 100644 index 0000000..c9739e6 --- /dev/null +++ b/src/pysol/filaments.py @@ -0,0 +1,173 @@ +"""Filament coil calculations. + +This is a dictionary based API, making it easier to define coils. +probably it should be made into a proper object class +but really I only use it for benchmarking against LDX values +and I wanted to copy from old Mathematica routines and other tests. + +Filaments are defined as an numpy 3 element vector + - r, z, and n. + +A coil is defined as an dictionary with +r1, r2, z1, and z2 and nt turns +""" + +from .inductance import ( + L_long_solenoid_butterworth, + L_lyle4, + L_lyle6, + L_lyle6_appendix, + L_maxwell, + L_thin_wall_babic_akyel, + L_thin_wall_lorentz, + dLdR_lyle6, + filament_coil, + mutual_inductance_of_filaments, + radial_force_of_filaments, + self_inductance_by_filaments, + vertical_force_of_filaments, +) + + +def FilamentCoil(C: dict, nr: int, nz: int): + """Break a coil into filaments, as nr by nz array. Adds the filaments to the coil. + + Args: + C (dict): Coil dictionary (class?) + nr (int): number of columns to break the coil into + nz (int): number of rows to break the coil into + + Returns: + dict: Coil dictionary with filaments added + """ + C["fil"] = filament_coil( + float((C["r1"] + C["r2"]) / 2), + float((C["z2"] + C["z1"]) / 2), + float(C["r2"] - C["r1"]), + float(C["z2"] - C["z1"]), + C["nt"], + nr, + nz, + ) + C["dr"] = (C["r2"] - C["r1"]) / nr + C["dz"] = (C["z2"] - C["z1"]) / nz + return C + + +def TotalM0(C1, C2): + """Mutual inductance of two coils by filamentation.""" + return mutual_inductance_of_filaments(C1["fil"], C2["fil"]) + + +def TotalFz(C1, C2): + """Vertical force of two coils by filamentation.""" + F_a2 = vertical_force_of_filaments(C1["fil"], C2["fil"]) + return C1["at"] / C1["nt"] * C2["at"] / C2["nt"] * F_a2 + + +def TotalFrOn1(C1, C2): + """Radial force of on coil 1 by filamentation from second coil and self force.""" + Fr_11 = ( + 0.5 + * C1["at"] + / C1["nt"] + * dLdR_lyle6( + (C1["r2"] + C1["r1"]) / 2.0, + C1["r2"] - C1["r1"], + C1["z2"] - C1["z1"], + C1["nt"], + ) + ) + + Fr_12 = ( + C1["at"] + / C1["nt"] + * C2["at"] + / C2["nt"] + * radial_force_of_filaments(C1["fil"], C2["fil"]) + ) + return Fr_11 + Fr_12 + + +def LMaxwell(C): + """Inductance by Maxwell's formula.""" + return L_maxwell( + float((C["r1"] + C["r2"]) / 2), + float(C["r2"] - C["r1"]), + float(C["z2"] - C["z1"]), + C["nt"], + ) + + +def LLyle4(C): + """Inductance by Lyle's formula, 4th order.""" + return L_lyle4( + float((C["r1"] + C["r2"]) / 2), + float(C["r2"] - C["r1"]), + float(C["z2"] - C["z1"]), + C["nt"], + ) + + +def LLyle6(C): + """Inductance by Lyle's formula, 6th order.""" + return L_lyle6( + float((C["r1"] + C["r2"]) / 2), + float(C["r2"] - C["r1"]), + float(C["z2"] - C["z1"]), + C["nt"], + ) + + +def LLyle6A(C): + """Inductance by Lyle's formula, 6th order, appendix.""" + return L_lyle6_appendix( + float((C["r1"] + C["r2"]) / 2), + float(C["r2"] - C["r1"]), + float(C["z2"] - C["z1"]), + C["nt"], + ) + + +def Lfil(Coil): + """Inductance by filamentation. + + Args: + Coil (dict): coil dictionary, must have run FilamentCoil first. + + Returns: + float: inductance of coil in Henries + """ + return self_inductance_by_filaments( + Coil["fil"], conductor="rect", dr=Coil["dr"], dz=Coil["dz"] + ) + + +def LLS(C): + """Inductance by Butterworth's formula for a long solenoid.""" + return L_long_solenoid_butterworth( + float((C["r1"] + C["r2"]) / 2), + float(C["r2"] - C["r1"]), + float(C["z2"] - C["z1"]), + C["nt"], + ) + + +def LBA(C): + """Inductance by Babic and Akyel's formula for a thin wall solenoid.""" + return L_thin_wall_babic_akyel( + float((C["r1"] + C["r2"]) / 2), + float(C["r2"] - C["r1"]), + float(C["z2"] - C["z1"]), + C["nt"], + ) + + +def LL(C): + """Inductance by Lorentz's formula for a thin wall solenoid.""" + return L_thin_wall_lorentz( + float((C["r1"] + C["r2"]) / 2), + float(C["r2"] - C["r1"]), + float(C["z2"] - C["z1"]), + C["nt"], + ) diff --git a/src/pysol/inductance.py b/src/pysol/inductance.py index 45b6843..387202b 100644 --- a/src/pysol/inductance.py +++ b/src/pysol/inductance.py @@ -1,25 +1,23 @@ -""" - inductance.py - - author: Darren Garnier +"""Inductance calculations for coils. - basic equations come from analytic approximations of old. - tested in real life with the LDX Fcoil / Ccoil / Lcoil system. +author: Darren Garnier - One from Maxwell himself, and better ones from: +basic equations come from analytic approximations of old. +tested in real life with the LDX Fcoil / Ccoil / Lcoil system. - Lyle, T. R. "On the Self-inductance of Circular Coils of - Rectangular Section". Roy. Soc. London. A. V213 (1914) pp 421-435. - https://doi.org/10.1098/rsta.1914.0009 +One from Maxwell himself, and better ones from: - Unfortunately, Lyle doesn't work that well with large dz/R coils. Other approximations - are also included. +Lyle, T. R. "On the Self-inductance of Circular Coils of + Rectangular Section". Roy. Soc. London. A. V213 (1914) pp 421-435. + https://doi.org/10.1098/rsta.1914.0009 - This code now uses numba to do just-in-time compiliation and parallel execution - to greatly increase speed. Also requires numba-scipy for elliptical functions. - numba-scipy can be fragile and sometimes needs to be "updated" before installation - to the newest version of numba. +Unfortunately, Lyle doesn't work that well with large dz/R coils. Other +approximations are also included. +This code now uses numba to do just-in-time compiliation and parallel execution +to greatly increase speed. Also requires numba-scipy for elliptical functions. +numba-scipy can be fragile and sometimes needs to be "updated" before installation +to the newest version of numba. """ import numpy as np @@ -37,7 +35,7 @@ ) warn_explicit(warning, RuntimeWarning, "inductance.py", 0) - def jit(*args, **kwargs): + def _jit(*args, **kwargs): if len(args) == 1 and len(kwargs) == 0 and callable(args[0]): # called as @decorator return args[0] @@ -45,10 +43,7 @@ def jit(*args, **kwargs): # called as @decorator(*args, **kwargs) return lambda f: f - njit = jit - prange = range - - def guvectorize(*args, **kwds): + def _guvectorize(*args, **kwds): def fake_decorator(f): warning = f"{f.__name__} requires Numba to be installed." warn_explicit(warning, RuntimeWarning, "inductance.py", 0) @@ -56,19 +51,24 @@ def fake_decorator(f): return fake_decorator + guvectorize = _guvectorize + njit = _jit + prange = range + jit = _jit + from .elliptics import ellipke try: - from mpmath import fp, mp + from mpmath import mp USE_MPMATH = True -except (ImportError, ModuleNotFoundError): +except ImportError: USE_MPMATH = False @njit -def lyle_terms(b, c): +def _lyle_terms(b, c): # helper functions for most self inductance equations. # b : cylindrical height of coil # c : radial width of coil @@ -85,7 +85,7 @@ def lyle_terms(b, c): if USE_MPMATH: - def lyle_terms_mp(b, c): + def _lyle_terms_mp(b, c): # b : cylindrical height of coil # c : radial width of coil d = mp.sqrt(b**2 + c**2) # diagnonal length @@ -100,38 +100,74 @@ def lyle_terms_mp(b, c): @njit def L_maxwell(r, dr, dz, n): - # calculation of parameters used in approximations - # inputs are: - # r : Major Radius - # dz : cylindrical height of coil - # dr : radial width of coil - # n : number of turns + """Self inductance of a rectangular coil with constant current density by Maxwell. + + Args: + r (float): coil centerline radius + dr (float): coil radial width + dz (float): coil height + n (int): number of turns + + Returns: + float: coil self inductance in Henrys + """ a = float(r) b = float(dz) c = float(dr) - d, u, v, w, wp, phi, GMD = lyle_terms(b, c) + d, u, v, w, wp, phi, GMD = _lyle_terms(b, c) L = 4e-7 * np.pi * (n**2) * a * (np.log(8 * a / GMD) - 2) return L @njit -def L_circle(R, r, n): - L = n**2 * 4e-7 * np.pi * (n**2) * R * (np.log(8 * R / r) - 1.75) +def L_round(r, a, n): + """Self inductance of a round conductor coil with constant current density. + + Args: + r (float): coil centerline radius + a (float): coil conductor radius + n (int): number of turns + + Returns: + float: coil self inductance in Henrys + """ + L = 4e-7 * np.pi * (n**2) * r * (np.log(8 * r / a) - 1.75) + return L + + +@njit +def L_hollow_round(r, a, n): + """Self inductance of a round conductor coil with skin current. + + Args: + r (float): coil centerline radius + a (float): coil conductor radius + n (int): number of turns + + Returns: + float: coil self inductance in Henrys + """ + L = 4e-7 * np.pi * (n**2) * r * (np.log(8 * r / a) - 2) return L @njit def L_lyle4_eq3(r, dr, dz, n): - # calculation of parameters used in approximations - # inputs are: - # r : Major Radius - # dz : cylindrical height of coil - # dr : radial width of coil - # n : number of turns + """Self inductance of a rectangular coil via Lyle to 4th order, Eq3. + + Args: + r (float): coil centerline radius + dr (float): coil radial width + dz (float): coil height + n (int): number of turns + + Returns: + float: coil self inductance in Henrys + """ a = float(r) b = float(dz) c = float(dr) - d, u, v, w, wp, phi, GMD = lyle_terms(b, c) + d, u, v, w, wp, phi, GMD = _lyle_terms(b, c) p2 = 1 / (2**5 * 3 * d**2) * (3 * b**2 + c**2) q2 = ( 1 @@ -184,16 +220,21 @@ def L_lyle4_eq3(r, dr, dz, n): # equation 3 above matches what I did for the 6th order. @njit def L_lyle4_eq4(r, dr, dz, n): - # calculation of parameters used in approximations - # inputs are: - # r : Major Radius - # dz : cylindrical height of coil - # dr : radial width of coil - # n : number of turns + """Self inductance of a rectangular coil via Lyle to 4th order, Eq4. + + Args: + r (float): coil centerline radius + dr (float): coil radial width + dz (float): coil height + n (int): number of turns + + Returns: + float: coil self inductance in Henrys + """ a = float(r) b = float(dz) c = float(dr) - d, u, v, w, wp, phi, GMD = lyle_terms(b, c) + d, u, v, w, wp, phi, GMD = _lyle_terms(b, c) p2 = 1 / (2**5 * 3 * d**2) * (3 * b**2 + c**2) q2 = ( 1 @@ -246,16 +287,21 @@ def L_lyle4_eq4(r, dr, dz, n): @njit def L_lyle6(r, dr, dz, n): - # calculation of parameters used in approximations - # inputs are: - # r : Major Radius - # dz : cylindrical height of coil - # dr : radial width of coil - # n : number of turns + """Self inductance of a rectangular coil via Lyle to 6th order. + + Args: + r (float): coil centerline radius + dr (float): coil radial width + dz (float): coil height + n (int): number of turns + + Returns: + float: coil self inductance in Henrys + """ a = float(r) b = float(dz) c = float(dr) - d, u, v, w, ww, phi, GMD = lyle_terms(b, c) + d, u, v, w, ww, phi, GMD = _lyle_terms(b, c) bd2 = (b / d) ** 2 cd2 = (c / d) ** 2 da2 = (d / a) ** 2 @@ -303,17 +349,21 @@ def L_lyle6(r, dr, dz, n): @njit def dLdR_lyle6(r, dr, dz, n): - # calculation of parameters used in approximations - # inputs are: - # r : Major Radius - # dz : cylindrical height of coil - # dr : radial width of coil - # n : number of turns + """Radial derivative of self inductance of a rectangular coil via Lyle to 6th order. + Args: + r (float): coil centerline radius + dr (float): coil radial width + dz (float): coil height + n (int): number of turns + + Returns: + float: radial derivative of inductance in Henrys/meter + """ a = float(r) b = float(dz) c = float(dr) - d, u, v, w, ww, phi, GMD = lyle_terms(b, c) + d, u, v, w, ww, phi, GMD = _lyle_terms(b, c) bd2 = (b / d) ** 2 cd2 = (c / d) ** 2 da2 = (d / a) ** 2 @@ -361,16 +411,21 @@ def dLdR_lyle6(r, dr, dz, n): @njit def L_lyle6_appendix(r, dr, dz, n): - # calculation of parameters used in approximations - # inputs are: - # r : Major Radius - # dz : cylindrical height of coil - # dr : radial width of coil - # n : number of turns + """Self inductance of a rectangular coil via Lyle to 6th order, appendix. + + Args: + r (float): coil centerline radius + dr (float): coil radial width + dz (float): coil height + n (int): number of turns + + Returns: + float: coil self inductance in Henrys + """ a = float(r) b = float(dz) c = float(dr) - d, u, v, w, wp, phi, GMD = lyle_terms(b, c) + d, u, v, w, wp, phi, GMD = _lyle_terms(b, c) p6 = ( 1 @@ -405,18 +460,25 @@ def L_lyle6_appendix(r, dr, dz, n): if USE_MPMATH: def L_lyle6_mp(r, dr, dz, n): - # use arbitrary precision library - # calculation of parameters used in approximations - # inputs are: - # r : Major Radius - # dz : cylindrical height of coil - # dr : radial width of coil - # n : number of turns + """Self inductance of a rectangular coil via Lyle to 6th order. + + This function uses the mpmath arbitrary precision library to + calculate to arbitrary precision. The default precision is 30. + + Args: + r (float): coil centerline radius + dr (float): coil radial width + dz (float): coil height + n (int): number of turns + + Returns: + float: coil self inductance in Henrys + """ mp.dps = 30 a = mp.mpf(r) b = mp.mpf(dz) c = mp.mpf(dr) - d, u, v, w, wp, phi, GMD = lyle_terms_mp(b, c) + d, u, v, w, wp, phi, GMD = _lyle_terms_mp(b, c) ML = mp.log(8 * a / d) f = ( @@ -475,6 +537,17 @@ def L_lyle6_mp(r, dr, dz, n): @njit def L_long_solenoid(r, dr, dz, n): + """Self inductance of a long solenoid. + + Args: + r (float): coil centerline radius + dr (float): coil radial width (ignored) + dz (float): coil height + n (int): number of turns + + Returns: + float: coil self inductance in Henrys + """ a = float(r) b = float(dz) # c = float(dr) @@ -484,6 +557,17 @@ def L_long_solenoid(r, dr, dz, n): @njit def L_long_solenoid_butterworth(r, dr, dz, n): + """Self inductance of a long solenoid by Butterworth's formula. + + Args: + r (float): coil centerline radius + dr (float): coil radial width + dz (float): coil height + n (int): number of turns + + Returns: + float: coil self inductance in Henrys + """ a = float(r) b = float(dz) c = float(dr) @@ -531,6 +615,7 @@ def L_long_solenoid_butterworth(r, dr, dz, n): - 1.0 / 96 * (c / a) ** 3 * (a / b) ** 2 * (1 - 39.0 / 10 * (a / b) ** 2) ) ) + D # hmm.. not using D, wonder why I bothered to calculate it.. lets look back at this function L = 4e-7 * np.pi * (n**2) * (a / b) * (K) return L @@ -538,10 +623,20 @@ def L_long_solenoid_butterworth(r, dr, dz, n): @njit def L_thin_wall_babic_akyel(r, dr, dz, n): - a = float(dz) / 2 - R1 = float(r) - float(dr) / 2 - R2 = float(r) + float(dr) / 2 + """Self inductance thin wall solenoid by Babic and Akyel's formula. + Args: + r (float): coil centerline radius + dr (float): coil radial width + dz (float): coil height + n (int): number of turns + + Returns: + float: coil self inductance in Henrys + """ + a = float(dz) / 2 + # R1 = float(r) - float(dr) / 2 + # R2 = float(r) + float(dr) / 2 # alpha = R2 / R1 # lets double check why I'm not using alpha... beta = a / r @@ -560,7 +655,18 @@ def L_thin_wall_babic_akyel(r, dr, dz, n): @njit def L_thin_wall_lorentz(r, dr, dz, n): - a = float(dz) / 2 + """Self inductance of a thin wall solenoid by Lorentz's formula. + + Args: + r (float): coil centerline radius + dr (float): coil radial width + dz (float): coil height + n (int): number of turns + + Returns: + float: coil self inductance in Henrys + """ + # a = float(dz) / 2 # beta = a / r k2 = 4 * r**2 / (4 * r**2 + dz**2) @@ -604,8 +710,15 @@ def L_thin_wall_lorentz(r, dr, dz, n): @njit def mutual_inductance_fil(rzn1, rzn2): - # calculatate the mutual inductance of two filaments defined as r, z, and n. + """Mutual inductance of two filaments. + + Args: + rzn1 (array): (r, z, n) of first filament + rzn2 (array): (r, z, n) of second filament + Returns: + float: mutual inductance in Henrys + """ r1 = rzn1[0] r2 = rzn2[0] z1 = rzn1[1] @@ -622,9 +735,15 @@ def mutual_inductance_fil(rzn1, rzn2): @njit def vertical_force_fil(rzn1, rzn2): - # calculatate the vertical force per conductor amp. - # of two filaments defined as r, z, and n. + """Vertical force between two filaments per conductor amp. + + Args: + rzn1 (array): (r, z, n) of first filament + rzn2 (array): (r, z, n) of second filament + Returns: + float: force in Newtons/Amp^2 + """ r1 = rzn1[0] r2 = rzn2[0] z1 = rzn1[1] @@ -639,9 +758,15 @@ def vertical_force_fil(rzn1, rzn2): @njit def radial_force_fil(rzn1, rzn2): - # calculatate the radial force per conductor amp. - # of two filaments defined as r, z, and n. + """Radial force between two filaments per conductor amp. + Args: + rzn1 (array): (r, z, n) of first filament + rzn2 (array): (r, z, n) of second filament + + Returns: + float: force in Newtons/Amp^2 + """ r1 = rzn1[0] r2 = rzn2[0] z1 = rzn1[1] @@ -661,8 +786,17 @@ def radial_force_fil(rzn1, rzn2): @njit def AGreen(r, z, a, b): - # PSI green's function of a filament with radius a, and z postion b, - # calculated at r,z + """Psi at position r, z, due to a filament with radius a, and z postion b. + + Args: + r (float): radial position of a point + z (float): vertical position of a point + a (float): radius of a filament + b (float): vertical position of a filament + + Returns: + float: Psi at r, z in Weber / Amp + """ k2 = 4 * a * r / ((r + a) ** 2 + (z - b) ** 2) elk, ele = ellipke(k2) amp = 4e-7 * a / np.sqrt((r + a) ** 2 + (z - b) ** 2) @@ -671,8 +805,17 @@ def AGreen(r, z, a, b): @njit def BrGreen(r, z, a, b): - # Br green's function of a filament with radius a, and z postion b, - # calculated at r,z + """Br at position r, z, due to a filament with radius a, and z postion b. + + Args: + r (float): radial position of a point + z (float): vertical position of a point + a (float): radius of a filament + b (float): vertical position of a filament + + Returns: + float: Br at r, z in Tesla / Amp + """ k2 = 4 * a * r / ((r + a) ** 2 + (z - b) ** 2) elk, ele = ellipke(k2) amp = -2e-7 / np.sqrt((r + a) ** 2 + (z - b) ** 2) @@ -682,8 +825,17 @@ def BrGreen(r, z, a, b): @njit def BzGreen(r, z, a, b): - # Bz green's function of a filament with radius a, and z postion b, - # calculated at r,z + """Bz at position r, z, due to a filament with radius a, and z postion b. + + Args: + r (float): radial position of a point + z (float): vertical position of a point + a (float): radius of a filament + b (float): vertical position of a filament + + Returns: + float: Bz at r, z in Tesla / Amp + """ k2 = 4 * a * r / ((r + a) ** 2 + (z - b) ** 2) elk, ele = ellipke(k2) amp = -2e-7 / np.sqrt((r + a) ** 2 + (z - b) ** 2) @@ -693,8 +845,17 @@ def BzGreen(r, z, a, b): @njit(parallel=True) def green_sum_over_filaments(gfunc, fil, r, z): - # calculate the greens function over grid of r, z - # can take different shaped r & z inputs + """Calculate a greens function at position r, z, to an array of filaments. + + Args: + gfunc (function): Green's function to use + fil (array): filament array N x (r, z, n) + r (float): radial position of a point + z (float): vertical position of a point + + Returns: + float: sum of the greens function at r, z, due to all filaments + """ tmp = np.zeros_like(r) tshp = tmp.shape tmp = tmp.reshape((tmp.size,)) @@ -708,6 +869,17 @@ def green_sum_over_filaments(gfunc, fil, r, z): def segment_path(pts, ds=0, close=False): + """Segment a path into equal length segments. + + Args: + pts (array): N x 3 array of points + ds (float, optional): length between points. Defaults to minimun in pts. + close (bool, optional): Should the path close the points. Defaults to False. + + Returns: + array: M x 3 array of points along the path + array: M array of length along the path + """ # this will make random points have equal spaces along path # from scipy.interpolate import interp1d # assume path is x,y,z and the 0 axis is the segment @@ -723,17 +895,13 @@ def segment_path(pts, ds=0, close=False): slen = s[-1] # length snew = np.linspace(0, slen, int(slen / ds)) - # oldway with interp1d - # segf = interp1d(s, pts, axis=0) - # segs = segf(snew) - # use np.interp instead of scipy.interpolation.interp1d segs = np.column_stack([np.interp(snew, s, pts[:, i]) for i in range(3)]) return segs, snew @njit -def loop_segmented_mutual(r, z, pts): +def _loop_segmented_mutual(r, z, pts): # segments is array of n x 3 (x,y,z) # segments should contain first point at start AND end. # r, z is r & z of loop @@ -751,33 +919,46 @@ def loop_segmented_mutual(r, z, pts): @njit(parallel=True) def mutual_filaments_segmented(fils, pts): + """Mutual inductance between a set of axisymmetric filaments and a segmented path.""" M = float(0) for i in prange(fils.shape[0]): - M += fils[i, 2] * loop_segmented_mutual(fils[i, 0], fils[i, 1], pts) + M += fils[i, 2] * _loop_segmented_mutual(fils[i, 0], fils[i, 1], pts) return M def M_filsol_path(fil, pts, n, ds=0): - # + """Mutual inductance between a set of axisymmetric filaments and a path from pts.""" segs, s = segment_path(pts, ds) return n * mutual_filaments_segmented(fil, segs) @njit(parallel=True) def segmented_self_inductance(pts, s, a): - # Neumann's formula.. double curve integral of - # mu_0/(4*pi) * int * int ( dx1 . dx2 / norm(x1-x2))) - # do all points except where x1-x2 blows up. - # instead follow Dengler https://doi.org/10.7716/aem.v5i1.331 - # which makes a approximation for that is good O(mu_0*a) - # lets just assume the thing is broken into small pieces and neglect end points - # - # this code doesn't work very well.. comparison test sorta fails - # phi = np.linspace(0,2*np.pi,100) - # test_xyz = np.array([[np.cos(p), np.sin(p), 0] for p in phi]) - # L_maxwell(1, .01, .01, 1), L_lyle6(1, .01, .01, 1), L_approx_path_rect(test_xyz, .01, .01, 1, .1) - # (6.898558527897293e-06, 6.8985981243869525e-06, 6.907313505254537e-06) # with Y=1/4 hmmm... - # (6.898558527897293e-06, 6.8985981243869525e-06, 7.064366776069971e-06) # with Y=1/2 + """Self inductance of a segmented path by double integral. + + Args: + pts (N x 3 array): points along the path + s (array): length along the path + a (float): radius of wire + + Returns: + float: self inductance of the path + + Neumann's formula.. double curve integral of + mu_0/(4*pi) * int * int ( dx1 . dx2 / norm(x1-x2))) + do all points except where x1-x2 blows up. + instead follow Dengler https://doi.org/10.7716/aem.v5i1.331 + which makes a approximation for that is good O(mu_0*a) + lets just assume the thing is broken into small pieces and neglect end points + + this code doesn't work very well.. comparison test sorta fails + phi = np.linspace(0,2*np.pi,100) + test_xyz = np.array([[np.cos(p), np.sin(p), 0] for p in phi]) + L_seg = L_approx_path_rect(test_xyz, .01, .01, 1, .1) + L_maxwell(1, .01, .01, 1), L_lyle6(1, .01, .01, 1), L_seg + (6.898558527897293e-06, 6.8985981243869525e-06, 6.907313505254537e-06) # Y=1/4 hmmm... + (6.898558527897293e-06, 6.8985981243869525e-06, 7.064366776069971e-06) # Y=1/2 + """ ds = s[1] - s[0] # the real ds and thus the real b dx = pts[1:] - pts[:-1] x = (pts[1:] + pts[:-1]) / 2 @@ -795,9 +976,12 @@ def segmented_self_inductance(pts, s, a): def L_approx_path_rect(pts, b, c, n, ds=1): - # take a path of points n x 3, with a cross section b x c - # and approximate self inductance using Dengler - _, _, _, _, _, _, a = lyle_terms( + """Approximate self inductance of a path of points with a rectangular cross section. + + take a path of points n x 3, with a cross section b x c + and approximate self inductance using Dengler + """ + _, _, _, _, _, _, a = _lyle_terms( b, c ) # get Maxwell mean radius to approximate "wire" radius ds *= a @@ -808,14 +992,20 @@ def L_approx_path_rect(pts, b, c, n, ds=1): # for some reason, these are slightly slower than the above. # but they work for any shape r & z as long as they are the same - - @guvectorize( ["void(float64[:,:], float64[:], float64[:], float64[:])"], "(p, q),(n),(n)->(n)", target="parallel", ) def BrGreenFil(fil, r, z, gr): + """Radial magnetic field greens functions from a set of filaments over a set of points. + + Args: + fil (_type_): filament array (r,z,n) + r (_type_): radial positions of points + z (_type_): vertical positions of points + gr (_type_): greens function values at points + """ for j in range(r.shape[0]): tmp = 0.0 for i in range(len(fil)): @@ -829,6 +1019,14 @@ def BrGreenFil(fil, r, z, gr): target="parallel", ) def BzGreenFil(fil, r, z, gr): + """Vertical magnetic field greens functions from a set of filaments over a set of points. + + Args: + fil (_type_): filament array (r,z,n) + r (_type_): radial positions of points + z (_type_): vertical positions of points + gr (_type_): greens function values at points + """ for j in range(r.shape[0]): tmp = 0.0 for i in range(len(fil)): @@ -842,6 +1040,14 @@ def BzGreenFil(fil, r, z, gr): target="parallel", ) def AGreenFil(fil, r, z, gr): + """Psi greens functions from a set of filaments over a set of points. + + Args: + fil (_type_): filament array (r,z,n) + r (_type_): radial positions of points + z (_type_): vertical positions of points + gr (_type_): greens function values at points + """ for j in range(r.shape[0]): tmp = 0.0 for i in range(len(fil)): @@ -851,8 +1057,7 @@ def AGreenFil(fil, r, z, gr): # cant jit this... meshgrid not allowed def filament_coil(r, z, dr, dz, nt, nr, nz, rot=0): - """Create an array of filaments, each with its own - radius, height, and amperage. + """Create an array of filaments, each with its own radius, height, and amperage. r : Major radius of coil center. z : Vertical center of coil. @@ -879,15 +1084,37 @@ def filament_coil(r, z, dr, dz, nt, nr, nz, rot=0): @njit(parallel=True) def sum_over_filaments(func, f1, f2): - S = float(0) + """Apply a function and sum over all combinations of two sets of filaments. + + Args: + func (function): function to apply to each pair of filaments + f1 (array): first filament array + f2 (array): second filament array + + Returns: + float: result of summing over all combinations of func(f1[i], f2[j]) + """ + result = float(0) for i in prange(f1.shape[0]): for j in range(f2.shape[0]): - S += func(f1[i], f2[j]) - return S + result += func(f1[i], f2[j]) + return result @njit(parallel=True) def mutual_inductance_of_filaments(f1, f2): + """Mutual inductance of two sets of filaments. + + These should not be the same filament array, + not setup to handle self inductance of filament. + + Args: + f1 (array): first filament array + f2 (array): second filament array + + Returns: + float: Mutual inductance of filament sets in Henries + """ M = float(0) for i in prange(f1.shape[0]): for j in range(f2.shape[0]): @@ -896,18 +1123,47 @@ def mutual_inductance_of_filaments(f1, f2): @njit(parallel=True) -def self_inductance_by_filaments(f, conductor="round", rc=0, hc=0, wc=0): +def self_inductance_by_filaments(f, conductor="round", a=0.01, dr=0.01, dz=0.01): + """Self inductance of filament set. + + Args: + f (array): first filament array + conductor (str, optional): conductor shape. Defaults to "round". + a (float, optional): conductor radius. Defaults to 0.01. + dr (float, optional): conductor radial width. Defaults to 0.01 + dz (float, optional): conductor vertical height. Defaults to 0.01 + + Returns: + float: self inductance of filament set in Henries + """ L = float(0) for i in prange(f.shape[0]): for j in range(f.shape[0]): if i != j: L += mutual_inductance_fil(f[i, :], f[j, :]) - L += self_inductance_fil(f[i, :], conductor) + if conductor == "round": + L += L_round(f[i, 0], a, f[i, 2]) + elif conductor == "hollow_round": + L += L_hollow_round(f[i, 0], a, f[i, 2]) + elif conductor == "rect": + L += L_lyle6(f[i, 0], dr, dz, f[i, 2]) return L @njit(parallel=True) def vertical_force_of_filaments(f1, f2): + """Vertical force between two sets of filaments. + + These should not be the same filament array, + not setup to handle self inductance of filament. + + Args: + f1 (array): first filament array + f2 (array): second filament array + + Returns: + float: Vertical force in Newtons/Amp^2 + """ F = float(0) for i in prange(f1.shape[0]): for j in range(f2.shape[0]): @@ -917,130 +1173,20 @@ def vertical_force_of_filaments(f1, f2): @njit(parallel=True) def radial_force_of_filaments(f1, f2): + """Radial force between two sets of filaments. + + These should not be the same filament array, + not setup to handle self inductance of filament. + + Args: + f1 (array): first filament array + f2 (array): second filament array + + Returns: + float: Radial force in Newtons/Amp^2 + """ F = float(0) for i in prange(f1.shape[0]): for j in range(f2.shape[0]): F += radial_force_fil(f1[i, :], f2[j, :]) return F - - -# This is a dictionary based API, making it easier to define coils. -# probably it should be made into a proper object class -# but really I only use it for benchmarking against LDX values -# and I wanted to copy from old Mathematica routines -# and other testing -# assumes coil is defined with r1, r2, z1, z2 with nt turns. -# one can then filament the coil nr x nz times - - -def FilamentCoil(C, nr, nz): - return filament_coil( - float((C["r1"] + C["r2"]) / 2), - float((C["z2"] + C["z1"]) / 2), - float(C["r2"] - C["r1"]), - float(C["z2"] - C["z1"]), - C["nt"], - nr, - nz, - ) - - -def TotalM0(C1, C2): - return mutual_inductance_of_filaments(C1["fil"], C2["fil"]) - - -def TotalFz(C1, C2): - F_a2 = vertical_force_of_filaments(C1["fil"], C2["fil"]) - return C1["at"] / C1["nt"] * C2["at"] / C2["nt"] * F_a2 - - -def TotalFrOn1(C1, C2): - Fr_11 = ( - 0.5 - * C1["at"] - / C1["nt"] - * dLdR_lyle6( - (C1["r2"] + C1["r1"]) / 2.0, - C1["r2"] - C1["r1"], - C1["z2"] - C1["z1"], - C1["nt"], - ) - ) - - Fr_12 = ( - C1["at"] - / C1["nt"] - * C2["at"] - / C2["nt"] - * radial_force_of_filaments(C1["fil"], C2["fil"]) - ) - return Fr_11 + Fr_12 - - -def LMaxwell(C): - return L_maxwell( - float((C["r1"] + C["r2"]) / 2), - float(C["r2"] - C["r1"]), - float(C["z2"] - C["z1"]), - C["nt"], - ) - - -def LLyle4(C): - return L_lyle4( - float((C["r1"] + C["r2"]) / 2), - float(C["r2"] - C["r1"]), - float(C["z2"] - C["z1"]), - C["nt"], - ) - - -def LLyle6(C): - return L_lyle6( - float((C["r1"] + C["r2"]) / 2), - float(C["r2"] - C["r1"]), - float(C["z2"] - C["z1"]), - C["nt"], - ) - - -def LLyle6A(C): - return L_lyle6_appendix( - float((C["r1"] + C["r2"]) / 2), - float(C["r2"] - C["r1"]), - float(C["z2"] - C["z1"]), - C["nt"], - ) - - -def LLyle6m(C): - return L_lyle6_mp( - ((C["r1"] + C["r2"]) / 2), (C["r2"] - C["r1"]), (C["z2"] - C["z1"]), C["nt"] - ) - - -def LLS(C): - return L_long_solenoid_butterworth( - float((C["r1"] + C["r2"]) / 2), - float(C["r2"] - C["r1"]), - float(C["z2"] - C["z1"]), - C["nt"], - ) - - -def LBA(C): - return L_thin_wall_babic_akyel( - float((C["r1"] + C["r2"]) / 2), - float(C["r2"] - C["r1"]), - float(C["z2"] - C["z1"]), - C["nt"], - ) - - -def LL(C): - return L_thin_wall_lorentz( - float((C["r1"] + C["r2"]) / 2), - float(C["r2"] - C["r1"]), - float(C["z2"] - C["z1"]), - C["nt"], - ) diff --git a/tests/test_elliptics.py b/tests/test_elliptics.py index 2be843b..cc16da4 100644 --- a/tests/test_elliptics.py +++ b/tests/test_elliptics.py @@ -1,10 +1,15 @@ +"""Test the elliptics module.""" + import unittest from pysol.elliptics import ellipke class TestElliptics(unittest.TestCase): + """Test the elliptics module.""" + def test_ellipke(self): + """Test the ellipke function.""" scipy_data = [ (0.00, 1.5707963267948966, 1.5707963267948966), (0.05, 1.591003453790792, 1.5509733517804725), diff --git a/tests/test_ldx.py b/tests/test_ldx.py index bbed21b..19d28df 100644 --- a/tests/test_ldx.py +++ b/tests/test_ldx.py @@ -1,10 +1,80 @@ +"""Test filamentation and inducatance of LDX coils.""" + +import unittest + import numpy as np -from pysol.inductance import (FilamentCoil, LLyle4, LLyle6, LLyle6A, LMaxwell, - TotalFz, TotalM0) +from pysol.filaments import ( + FilamentCoil, + Lfil, + LLyle4, + LLyle6, + LLyle6A, + LMaxwell, + TotalFz, + TotalM0, +) + + +class TestLDXInductance(unittest.TestCase): + """Test the inductance of the LDX coils. + + We know the coil weighed about 590kg, and had about 0.4H inductance. + """ + + def setUp(self) -> None: + """Set up the LDX coils.""" + ( + self.HTS_LC, + self.LC, + self.CC, + self.F1, + self.F2, + self.F3, + self.FC, + ) = define_ldx_coils() + return super().setUp() + + def test_self_inductances(self): + """Test the self-inductances routines.""" + self.assertAlmostEqual(LMaxwell(self.LC), 0.005397519485316731) + self.assertAlmostEqual(LLyle4(self.LC), 0.005623887363536865) + self.assertAlmostEqual(LLyle6(self.LC), 0.005626208118367906) + self.assertAlmostEqual(LLyle6A(self.LC), 0.005626208118367906) + self.assertAlmostEqual(Lfil(self.LC), 0.005625117051066614) + self.assertAlmostEqual(LMaxwell(self.CC), 85.91637858501646) + self.assertAlmostEqual(LLyle4(self.CC), 90.90254315310752) + self.assertAlmostEqual(LLyle6(self.CC), 90.90927053789774) + self.assertAlmostEqual(LLyle6A(self.CC), 90.90927053789773) + self.assertAlmostEqual(Lfil(self.CC), 90.89118772211005) -def run_ldx_tests(): + def test_mutual_inductances(self): + """Test the mutual inductances routines.""" + MFL = TotalM0(self.FC, self.LC) + MFC = TotalM0(self.FC, self.CC) + self.assertAlmostEqual(MFL, 0.611729e-3) + self.assertAlmostEqual(MFC, 1.686291576361124) + + def test_fcoil_selfinductance(self): + """Test the self-inductance of the F coils.""" + LF1 = LLyle6(self.F1) + LF2 = LLyle6(self.F2) + LF3 = LLyle6(self.F3) + MF12 = TotalM0(self.F1, self.F2) + MF13 = TotalM0(self.F1, self.F3) + MF23 = TotalM0(self.F2, self.F3) + LFC = LF1 + LF2 + LF3 + 2 * MF12 + 2 * MF13 + 2 * MF23 + self.assertAlmostEqual(LFC, 0.38692937382836284) + + def test_levitation_force(self): + """Test the levitation force.""" + FZ = TotalFz(self.FC, self.LC) / 9.81 + self.assertAlmostEqual(FZ, 589.263842397435) + + +def define_ldx_coils(): + """Define the LDX coils.""" HTS_LC = {} HTS_LC["r1"] = 0.41 / 2 HTS_LC["r2"] = 1.32 / 2 @@ -21,7 +91,7 @@ def run_ldx_tests(): LC["z2"] = LC["z1"] + 0.1 LC["nt"] = 80 LC["at"] = 3500 * 80 - LC["fil"] = FilamentCoil(LC, 20, 4) + LC = FilamentCoil(LC, 20, 4) CC = {} CC["r1"] = 0.645 @@ -30,7 +100,7 @@ def run_ldx_tests(): CC["z2"] = -0.002 + 0.750 / 2 CC["nt"] = 8388 CC["at"] = 8388 * 420 - CC["fil"] = FilamentCoil(CC, 3, 7) + CC = FilamentCoil(CC, 3, 7) F1 = {} F1["r1"] = 0.2717 - 0.01152 / 2 @@ -39,7 +109,7 @@ def run_ldx_tests(): F1["z2"] = +0.0694 / 2 F1["nt"] = 26.6 F1["at"] = 26.6 * 1629 - F1["fil"] = FilamentCoil(F1, 2, 4) + F1 = FilamentCoil(F1, 2, 4) F2 = {} F2["r1"] = 0.28504 - 0.01508 / 2 @@ -48,7 +118,7 @@ def run_ldx_tests(): F2["z2"] = +0.125 / 2 F2["nt"] = 81.7 F2["at"] = 81.7 * 1629 - F2["fil"] = FilamentCoil(F2, 3, 7) + F2 = FilamentCoil(F2, 3, 7) F3 = {} F3["r1"] = 0.33734 - 0.08936 / 2 @@ -57,36 +127,15 @@ def run_ldx_tests(): F3["z2"] = +0.1615 / 2 F3["nt"] = 607.7 F3["at"] = 607.7 * 1629 - F3["fil"] = FilamentCoil(F3, 10, 15) + F3 = FilamentCoil(F3, 10, 15) FC = {} FC["fil"] = np.concatenate((F1["fil"], F2["fil"], F3["fil"])) FC["nt"] = F1["nt"] + F2["nt"] + F3["nt"] FC["at"] = F1["at"] + F2["at"] + F3["at"] - print("Inductances of LC by approximations") - print(LMaxwell(LC), LLyle4(LC), LLyle6(LC), LLyle6A(LC)) - print("Inductances of CC by approximations") - print(LMaxwell(CC), LLyle4(CC), LLyle6(CC), LLyle6A(CC)) - - print("Self-Inductance of FC by approximations and mutuals") - LF1 = LLyle6(F1) - LF2 = LLyle6(F2) - LF3 = LLyle6(F3) - MF12 = TotalM0(F1, F2) - MF13 = TotalM0(F1, F3) - MF23 = TotalM0(F2, F3) - LFC = LF1 + LF2 + LF3 + 2 * MF12 + 2 * MF13 + 2 * MF23 - print("LFC %f H" % (LFC)) - - MFL = TotalM0(FC, LC) - MFC = TotalM0(FC, CC) - FFL = TotalFz(FC, LC) - - print("Mutual of F & L: %f mH" % (MFL * 1000)) - print("Mutual of F & C: %f H" % (MFC)) - print("Force between F & L in %6.2f kg" % (FFL / 9.81)) + return HTS_LC, LC, CC, F1, F2, F3, FC if __name__ == "__main__": - run_ldx_tests() + unittest.main()