Skip to content

Commit

Permalink
Yet another major restructuring. Fields are now loaded into their own…
Browse files Browse the repository at this point in the history 
… module instead of some user module.
  • Loading branch information
@michakraus
michakraus committed Sep 21, 2020
1 parent c5499e6 commit 9c5b559
Show file tree
Hide file tree
Showing 16 changed files with 666 additions and 550 deletions.
7 changes: 3 additions & 4 deletions src/ElectromagneticFields.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,13 +9,13 @@ module ElectromagneticFields

export AnalyticField, AnalyticEquilibrium, AnalyticPerturbation, ZeroPerturbation

export CartesianField, CartesianEquilibrium, CartesianPerturbation

export load_equilibrium, periodicity

include("analytic/analytic_equilibrium.jl")
include("analytic/cartesian_equilibrium.jl")

using .AnalyticCartesianField

export ABC, AxisymmetricTokamakCartesian, AxisymmetricTokamakCylindrical,
AxisymmetricTokamakToroidal, EzCosZ,
Solovev, SolovevXpoint, SolovevQuadratic, SymmetricQuadratic, ThetaPinch
Expand All @@ -36,8 +36,7 @@ module ElectromagneticFields
include("analytic/axisymmetric_tokamak_cylindrical.jl")
include("analytic/axisymmetric_tokamak_toroidal.jl")
include("analytic/ezcosz.jl")
include("analytic/solovev_common.jl")
include("analytic/solovev_psi.jl")
include("analytic/solovev_abstract.jl")
include("analytic/solovev.jl")
include("analytic/solovev_quadratic.jl")
include("analytic/solovev_xpoint.jl")
Expand Down
59 changes: 38 additions & 21 deletions src/analytic/abc.jl
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
@doc raw"""
ABC equilibrium in (x,y,z) coordinates with covariant components of the vector
ABCEquilibrium equilibrium in (x,y,z) coordinates with covariant components of the vector
potential given by
```math
A (x,y,z) = \big( a \, \sin(z) + c \, \cos(y) , \, b \, \sin(x) + a \, \cos(z) , \, c \, \sin(y) + b \, \cos(x) \big)^T
Expand All @@ -8,32 +8,49 @@ resulting in the magnetic field ``B(x,y,z) = A(x,y,z)``.
Parameters: `a`, `b`, `c`
"""
struct ABC{T <: Number} <: CartesianEquilibrium
name::String
a::T
b::T
c::T
module ABC

ABC{T}(a::T, b::T, c::T) where T <: Number = new("ABCEquilibrium", a, b, c)
end
import ..ElectromagneticFields
import ..ElectromagneticFields: CartesianEquilibrium, ZeroPerturbation
import ..ElectromagneticFields: load_equilibrium, generate_equilibrium_code
import ..AnalyticCartesianField: X, Y, Z

ABC(a::T, b::T, c::T) where T <: Number = ABC{T}(a, b, c)
export ABCEquilibrium

struct ABCEquilibrium{T <: Number} <: CartesianEquilibrium
name::String
a::T
b::T
c::T

function Base.show(io::IO, equ::ABC)
print(io, "ABC Equilibrium with\n")
print(io, " A = ", equ.a, "\n")
print(io, " B = ", equ.b, "\n")
print(io, " C = ", equ.c)
end
ABCEquilibrium{T}(a::T, b::T, c::T) where T <: Number = new("ABCEquilibrium", a, b, c)
end

ABCEquilibrium(a::T, b::T, c::T) where T <: Number = ABCEquilibrium{T}(a, b, c)


function Base.show(io::IO, equ::ABCEquilibrium)
print(io, "ABC Equilibrium with\n")
print(io, " A = ", equ.a, "\n")
print(io, " B = ", equ.b, "\n")
print(io, " C = ", equ.c)
end


ElectromagneticFields.A₁(x::AbstractVector, equ::ABCEquilibrium) = equ.a * sin(x[3]) + equ.c * cos(x[2])
ElectromagneticFields.A₂(x::AbstractVector, equ::ABCEquilibrium) = equ.b * sin(x[1]) + equ.a * cos(x[3])
ElectromagneticFields.A₃(x::AbstractVector, equ::ABCEquilibrium) = equ.c * sin(x[2]) + equ.b * cos(x[1])

ElectromagneticFields.get_functions(::ABCEquilibrium) = (X=X, Y=Y, Z=Z)

A₁(x::AbstractVector, equ::ABC) = equ.a * sin(x[3]) + equ.c * cos(x[2])
A₂(x::AbstractVector, equ::ABC) = equ.b * sin(x[1]) + equ.a * cos(x[3])
A₃(x::AbstractVector, equ::ABC) = equ.c * sin(x[2]) + equ.b * cos(x[1])
macro abc_equilibrium(a, b, c)
generate_equilibrium_code(ABCEquilibrium(a, b, c); output=false)
end

get_functions(::ABC) = (X=X, Y=Y, Z=Z)
function init(a, b, c; perturbation=ZeroPerturbation())
equilibrium = ABCEquilibrium(a, b, c)
load_equilibrium(equilibrium, perturbation; target_module=ABC)
return equilibrium
end

macro abc_equilibrium(a, b, c)
generate_equilibrium_code(ABC(a, b, c); output=false)
end
26 changes: 13 additions & 13 deletions src/analytic/analytic_equilibrium.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,15 +8,15 @@ abstract type AnalyticPerturbation <: AnalyticField end

get_functions(::AnalyticField) = (;)

x¹(::AbstractVector{T}, ::AnalyticField) where {T} = error("x¹() not implemented for ", ET)
x²(::AbstractVector{T}, ::AnalyticField) where {T} = error("x²() not implemented for ", ET)
x³(::AbstractVector{T}, ::AnalyticField) where {T} = error("x³() not implemented for ", ET)
x¹(::AbstractVector, ::ET) where {ET <: AnalyticField} = error("x¹() not implemented for ", ET)
x²(::AbstractVector, ::ET) where {ET <: AnalyticField} = error("x²() not implemented for ", ET)
x³(::AbstractVector, ::ET) where {ET <: AnalyticField} = error("x³() not implemented for ", ET)

ξ¹(::AbstractVector{T}, ::AnalyticField) where {T} = error("ξ¹() not implemented for ", ET)
ξ²(::AbstractVector{T}, ::AnalyticField) where {T} = error("ξ²() not implemented for ", ET)
ξ³(::AbstractVector{T}, ::AnalyticField) where {T} = error("ξ³() not implemented for ", ET)
ξ¹(::AbstractVector, ::ET) where {ET <: AnalyticField} = error("ξ¹() not implemented for ", ET)
ξ²(::AbstractVector, ::ET) where {ET <: AnalyticField} = error("ξ²() not implemented for ", ET)
ξ³(::AbstractVector, ::ET) where {ET <: AnalyticField} = error("ξ³() not implemented for ", ET)

J(::AbstractVector, ::AnalyticField) = error("J() not implemented for ", ET)
J(::AbstractVector, ::ET) where {ET <: AnalyticField} = error("J() not implemented for ", ET)

g₁₁(::AbstractVector{T}, ::AnalyticField) where {T} = one(T)
g₁₂(::AbstractVector{T}, ::AnalyticField) where {T} = zero(T)
Expand All @@ -33,9 +33,9 @@ periodicity(x::AbstractVector, ::AnalyticField) = zero(x)
from_cartesian(x::AbstractVector, equ::AnalyticField) = [ξ¹(x,equ), ξ²(x,equ), ξ³(x,equ)]
to_cartesian(ξ::AbstractVector, equ::AnalyticField) = [x¹(ξ,equ), x²(ξ,equ), x³(ξ,equ)]

A₁(::AbstractVector{T}, ::AnalyticField) where {T} = zero(T)
A₂(::AbstractVector{T}, ::AnalyticField) where {T} = zero(T)
A₃(::AbstractVector{T}, ::AnalyticField) where {T} = zero(T)
A₁(::AbstractVector, ::ET) where {ET <: AnalyticField} = error("A₁() not implemented for ", ET)
A₂(::AbstractVector, ::ET) where {ET <: AnalyticField} = error("A₂() not implemented for ", ET)
A₃(::AbstractVector, ::ET) where {ET <: AnalyticField} = error("A₃() not implemented for ", ET)

φ(::AbstractVector{T}, ::AnalyticField) where {T} = zero(T)

Expand All @@ -55,9 +55,9 @@ struct ZeroPerturbation <: AnalyticPerturbation
ZeroPerturbation() = new("ZeroPerturbation")
end

A₁(::AbstractVector{T}, ::ZeroPerturbation) where {T} = zero(T)
A₂(::AbstractVector{T}, ::ZeroPerturbation) where {T} = zero(T)
A₃(::AbstractVector{T}, ::ZeroPerturbation) where {T} = zero(T)
A₁(::AbstractVector{T}, ::AnalyticPerturbation) where {T} = zero(T)
A₂(::AbstractVector{T}, ::AnalyticPerturbation) where {T} = zero(T)
A₃(::AbstractVector{T}, ::AnalyticPerturbation) where {T} = zero(T)


"Returns the i-th component of the vector corresponding to the one-form α"
Expand Down
82 changes: 51 additions & 31 deletions src/analytic/axisymmetric_tokamak_cartesian.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,45 +11,65 @@ B (x,y,z) = \frac{B_0}{q_0} \, \bigg( - \frac{q_0 R_0 y + x z}{R^2} , \, \frac{q
where $R = \sqrt{ x^2 + y^2 }$ and $r = \sqrt{ (R - R_0)^2 + z^2 }$.
Parameters:
* `R₀`: position of magnetic axis
* `B₀`: B-field at magnetic axis
* `q₀`: safety factor at magnetic axis
* `R₀`: position of magnetic axis
* `B₀`: B-field at magnetic axis
* `q₀`: safety factor at magnetic axis
"""
struct AxisymmetricTokamakCartesian{T <: Number} <: CartesianEquilibrium
name::String
R₀::T
B₀::T
q₀::T

function AxisymmetricTokamakCartesian{T}(R₀::T, B₀::T, q₀::T) where T <: Number
new("AxisymmetricTokamakCartesianEquilibrium", R₀, B₀, q₀)
module AxisymmetricTokamakCartesian

import ..ElectromagneticFields
import ..ElectromagneticFields: CartesianEquilibrium, ZeroPerturbation
import ..ElectromagneticFields: load_equilibrium, generate_equilibrium_code
import ..AnalyticCartesianField: X, Y, Z

export AxisymmetricTokamakCartesianEquilibrium

const DEFAULT_R₀ = 1.0
const DEFAULT_B₀ = 1.0
const DEFAULT_q₀ = 2.0

struct AxisymmetricTokamakCartesianEquilibrium{T <: Number} <: CartesianEquilibrium
name::String
R₀::T
B₀::T
q₀::T

function AxisymmetricTokamakCartesianEquilibrium{T}(R₀::T, B₀::T, q₀::T) where T <: Number
new("AxisymmetricTokamakCartesianEquilibrium", R₀, B₀, q₀)
end
end
end

AxisymmetricTokamakCartesian(R₀::T=1.0, B₀::T=1.0, q₀::T=2.0) where T <: Number = AxisymmetricTokamakCartesian{T}(R₀, B₀, q₀)
AxisymmetricTokamakCartesianEquilibrium(R₀::T=DEFAULT_R₀, B₀::T=DEFAULT_B₀, q₀::T=DEFAULT_q₀) where T <: Number = AxisymmetricTokamakCartesianEquilibrium{T}(R₀, B₀, q₀)

function Base.show(io::IO, equ::AxisymmetricTokamakCartesianEquilibrium)
print(io, "Axisymmetric Tokamak Equilibrium in (x,y,z) Coordinates with\n")
print(io, " R₀ = ", equ.R₀, "\n")
print(io, " B₀ = ", equ.B₀, "\n")
print(io, " q₀ = ", equ.q₀)
end

function Base.show(io::IO, equ::AxisymmetricTokamakCartesian)
print(io, "Axisymmetric Tokamak Equilibrium in (x,y,z) Coordinates with\n")
print(io, " R₀ = ", equ.R₀, "\n")
print(io, " B₀ = ", equ.B₀, "\n")
print(io, " q₀ = ", equ.q₀)
end

R²(x::AbstractVector, equ::AxisymmetricTokamakCartesianEquilibrium) = X(x,equ)^2 + Y(x,equ)^2
r²(x::AbstractVector, equ::AxisymmetricTokamakCartesianEquilibrium) = (R(x,equ) - equ.R₀)^2 + Z(x,equ)^2
R(x::AbstractVector, equ::AxisymmetricTokamakCartesianEquilibrium) = sqrt(R²(x,equ))
r(x::AbstractVector, equ::AxisymmetricTokamakCartesianEquilibrium) = sqrt(r²(x,equ))
θ(x::AbstractVector, equ::AxisymmetricTokamakCartesianEquilibrium) = atan(Z(x,equ), R(x,equ) - equ.R₀)
ϕ(x::AbstractVector, equ::AxisymmetricTokamakCartesianEquilibrium) = atan(Y(x,equ), X(x,equ))

R²(x::AbstractVector, equ::AxisymmetricTokamakCartesian) = X(x,equ)^2 + Y(x,equ)^2
r²(x::AbstractVector, equ::AxisymmetricTokamakCartesian) = (R(x,equ) - equ.R₀)^2 + Z(x,equ)^2
R(x::AbstractVector, equ::AxisymmetricTokamakCartesian) = sqrt(R²(x,equ))
r(x::AbstractVector, equ::AxisymmetricTokamakCartesian) = sqrt(r²(x,equ))
θ(x::AbstractVector, equ::AxisymmetricTokamakCartesian) = atan(Z(x,equ), R(x,equ) - equ.R₀)
ϕ(x::AbstractVector, equ::AxisymmetricTokamakCartesian) = atan(Y(x,equ), X(x,equ))
ElectromagneticFields.A₁(x::AbstractVector, equ::AxisymmetricTokamakCartesianEquilibrium) = + equ.B₀ * (equ.R₀ * X(x,equ) * Z(x,equ) - r²(x,equ) * Y(x,equ) / equ.q₀ ) / R²(x,equ) / 2
ElectromagneticFields.A₂(x::AbstractVector, equ::AxisymmetricTokamakCartesianEquilibrium) = + equ.B₀ * (equ.R₀ * Y(x,equ) * Z(x,equ) + r²(x,equ) * X(x,equ) / equ.q₀ ) / R²(x,equ) / 2
ElectromagneticFields.A₃(x::AbstractVector, equ::AxisymmetricTokamakCartesianEquilibrium) = - equ.B₀ * equ.R₀ * log(R(x,equ) / equ.R₀) / 2

A₁(x::AbstractVector, equ::AxisymmetricTokamakCartesian) = + equ.B₀ * (equ.R₀ * X(x,equ) * Z(x,equ) - r²(x,equ) * Y(x,equ) / equ.q₀ ) / R²(x,equ) / 2
A₂(x::AbstractVector, equ::AxisymmetricTokamakCartesian) = + equ.B₀ * (equ.R₀ * Y(x,equ) * Z(x,equ) + r²(x,equ) * X(x,equ) / equ.q₀ ) / R²(x,equ) / 2
A₃(x::AbstractVector, equ::AxisymmetricTokamakCartesian) = - equ.B₀ * equ.R₀ * log(R(x,equ) / equ.R₀) / 2
ElectromagneticFields.get_functions(::AxisymmetricTokamakCartesianEquilibrium) = (X=X, Y=Y, Z=Z, R=R, r=r, θ=θ, ϕ=ϕ, R²=R², r²=r²)

get_functions(::AxisymmetricTokamakCartesian) = (X=X, Y=Y, Z=Z, R=R, r=r, θ=θ, ϕ=ϕ, R²=R², r²=r²)
macro axisymmetric_tokamak_equilibrium_cartesian(R₀, B₀, q₀)
generate_equilibrium_code(AxisymmetricTokamakCartesianEquilibrium(R₀, B₀, q₀); output=false)
end

function init(R₀=DEFAULT_R₀, B₀=DEFAULT_B₀, q₀=DEFAULT_q₀; perturbation=ZeroPerturbation())
equilibrium = AxisymmetricTokamakCartesianEquilibrium(R₀, B₀, q₀)
load_equilibrium(equilibrium, perturbation; target_module=AxisymmetricTokamakCartesian)
return equilibrium
end

macro axisymmetric_tokamak_equilibrium_cartesian(R₀, B₀, q₀)
generate_equilibrium_code(AxisymmetricTokamakCartesian(R₀, B₀, q₀); output=false)
end
111 changes: 63 additions & 48 deletions src/analytic/axisymmetric_tokamak_cylindrical.jl
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,3 @@
using SymEngine

@doc raw"""
Axisymmetric tokamak equilibrium in (R,Z,ϕ) coordinates with covariant
components of the vector potential given by
Expand All @@ -17,65 +15,82 @@ Parameters:
* `B₀`: B-field at magnetic axis
* `q₀`: safety factor at magnetic axis
"""
struct AxisymmetricTokamakCylindrical{T <: Number} <: AnalyticEquilibrium
name::String
R₀::T
B₀::T
q₀::T

function AxisymmetricTokamakCylindrical{T}(R₀::T, B₀::T, q₀::T) where T <: Number
new("AxisymmetricTokamakCylindricalEquilibrium", R₀, B₀, q₀)
end
end
module AxisymmetricTokamakCylindrical

AxisymmetricTokamakCylindrical(R₀::T=1.0, B₀::T=1.0, q₀::T=2.0) where T <: Number = AxisymmetricTokamakCylindrical{T}(R₀, B₀, q₀)
import ..ElectromagneticFields
import ..ElectromagneticFields: AnalyticEquilibrium, ZeroPerturbation
import ..ElectromagneticFields: load_equilibrium, generate_equilibrium_code

export AxisymmetricTokamakCylindricalEquilibrium

function Base.show(io::IO, equ::AxisymmetricTokamakCylindrical)
print(io, "Axisymmetric Tokamak Equilibrium in (R,Z,ϕ) Coordinates with\n")
print(io, " R₀ = ", equ.R₀, "\n")
print(io, " B₀ = ", equ.B₀, "\n")
print(io, " q₀ = ", equ.q₀)
end
const DEFAULT_R₀ = 1.0
const DEFAULT_B₀ = 1.0
const DEFAULT_q₀ = 2.0

struct AxisymmetricTokamakCylindricalEquilibrium{T <: Number} <: AnalyticEquilibrium
name::String
R₀::T
B₀::T
q₀::T

R(x::AbstractVector, equ::AxisymmetricTokamakCylindrical) = x[1]
Z(x::AbstractVector, equ::AxisymmetricTokamakCylindrical) = x[2]
ϕ(x::AbstractVector, equ::AxisymmetricTokamakCylindrical) = x[3]
r²(x::AbstractVector, equ::AxisymmetricTokamakCylindrical) = (R(x,equ) - equ.R₀)^2 + Z(x,equ)^2
r(x::AbstractVector, equ::AxisymmetricTokamakCylindrical) = sqrt(r²(x, equ))
X(x::AbstractVector, equ::AxisymmetricTokamakCylindrical) = R(x,equ) * cos(ϕ(x,equ))
Y(x::AbstractVector, equ::AxisymmetricTokamakCylindrical) = R(x,equ) * sin(ϕ(x,equ))
θ(x::AbstractVector, equ::AxisymmetricTokamakCylindrical) = atan(Z(x,equ), R(x,equ) - equ.R₀)
function AxisymmetricTokamakCylindricalEquilibrium{T}(R₀::T, B₀::T, q₀::T) where T <: Number
new("AxisymmetricTokamakCylindricalEquilibrium", R₀, B₀, q₀)
end
end

J(x::AbstractVector, equ::AxisymmetricTokamakCylindrical) = R(x,equ)
AxisymmetricTokamakCylindricalEquilibrium(R₀::T=DEFAULT_R₀, B₀::T=DEFAULT_B₀, q₀::T=DEFAULT_q₀) where T <: Number = AxisymmetricTokamakCylindricalEquilibrium{T}(R₀, B₀, q₀)

A₁(x::AbstractVector, equ::AxisymmetricTokamakCylindrical) = + equ.B₀ * equ.R₀ * Z(x,equ) / R(x,equ) / 2
A₂(x::AbstractVector, equ::AxisymmetricTokamakCylindrical) = - equ.B₀ * equ.R₀ * log(R(x,equ) / equ.R₀) / 2
A₃(x::AbstractVector, equ::AxisymmetricTokamakCylindrical) = - equ.B₀ * r²(x,equ) / equ.q₀ / 2
function Base.show(io::IO, equ::AxisymmetricTokamakCylindricalEquilibrium)
print(io, "Axisymmetric Tokamak Equilibrium in (R,Z,ϕ) Coordinates with\n")
print(io, " R₀ = ", equ.R₀, "\n")
print(io, " B₀ = ", equ.B₀, "\n")
print(io, " q₀ = ", equ.q₀)
end

x¹(ξ::AbstractVector, equ::AxisymmetricTokamakCylindrical) = X(ξ,equ)
x²(ξ::AbstractVector, equ::AxisymmetricTokamakCylindrical) = Y(ξ,equ)
x³(ξ::AbstractVector, equ::AxisymmetricTokamakCylindrical) = Z(ξ,equ)

ξ¹(x::AbstractVector, equ::AxisymmetricTokamakCylindrical) = sqrt(x[1]^2 + x[2]^2)
ξ²(x::AbstractVector, equ::AxisymmetricTokamakCylindrical) = x[3]
ξ³(x::AbstractVector, equ::AxisymmetricTokamakCylindrical) = atan(x[2], x[1])
R(x::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = x[1]
Z(x::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = x[2]
ϕ(x::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = x[3]
r²(x::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = (R(x,equ) - equ.R₀)^2 + Z(x,equ)^2
r(x::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = sqrt(r²(x, equ))
X(x::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = R(x,equ) * cos(ϕ(x,equ))
Y(x::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = R(x,equ) * sin(ϕ(x,equ))
θ(x::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = atan(Z(x,equ), R(x,equ) - equ.R₀)

g₁₁(x::AbstractVector, equ::AxisymmetricTokamakCylindrical) = one(eltype(x))
g₂₂(x::AbstractVector, equ::AxisymmetricTokamakCylindrical) = one(eltype(x))
g₃₃(x::AbstractVector, equ::AxisymmetricTokamakCylindrical) = R(x, equ)^2
ElectromagneticFields.J(x::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = R(x,equ)

ElectromagneticFields.A₁(x::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = + equ.B₀ * equ.R₀ * Z(x,equ) / R(x,equ) / 2
ElectromagneticFields.A₂(x::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = - equ.B₀ * equ.R₀ * log(R(x,equ) / equ.R₀) / 2
ElectromagneticFields.A₃(x::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = - equ.B₀ * r²(x,equ) / equ.q₀ / 2

get_functions(::AxisymmetricTokamakCylindrical) = (X=X, Y=Y, Z=Z, R=R, r=r, θ=θ, ϕ=ϕ, r²=r²)
ElectromagneticFields.x¹(ξ::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = X(ξ,equ)
ElectromagneticFields.x²(ξ::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = Y(ξ,equ)
ElectromagneticFields.x³(ξ::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = Z(ξ,equ)

function periodicity(x::AbstractVector, ::AxisymmetricTokamakCylindrical)
p = zero(x)
p[3] = 2π
return p
end
ElectromagneticFields.ξ¹(x::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = sqrt(x[1]^2 + x[2]^2)
ElectromagneticFields.ξ²(x::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = x[3]
ElectromagneticFields.ξ³(x::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = atan(x[2], x[1])

ElectromagneticFields.g₁₁(x::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = one(eltype(x))
ElectromagneticFields.g₂₂(x::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = one(eltype(x))
ElectromagneticFields.g₃₃(x::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = R(x, equ)^2

ElectromagneticFields.get_functions(::AxisymmetricTokamakCylindricalEquilibrium) = (X=X, Y=Y, Z=Z, R=R, r=r, θ=θ, ϕ=ϕ, r²=r²)

function ElectromagneticFields.periodicity(x::AbstractVector, ::AxisymmetricTokamakCylindricalEquilibrium)
p = zero(x)
p[3] = 2π
return p
end

macro axisymmetric_tokamak_equilibrium_cylindrical(R₀, B₀, q₀)
generate_equilibrium_code(AxisymmetricTokamakCylindricalEquilibrium(R₀, B₀, q₀); output=false)
end

function init(R₀=DEFAULT_R₀, B₀=DEFAULT_B₀, q₀=DEFAULT_q₀; perturbation=ZeroPerturbation())
equilibrium = AxisymmetricTokamakCylindricalEquilibrium(R₀, B₀, q₀)
load_equilibrium(equilibrium, perturbation; target_module=AxisymmetricTokamakCylindrical)
return equilibrium
end

macro axisymmetric_tokamak_equilibrium_cylindrical(R₀, B₀, q₀)
generate_equilibrium_code(AxisymmetricTokamakCylindrical(R₀, B₀, q₀); output=false)
end
Loading

0 comments on commit 9c5b559

Please sign in to comment.