From 849e6bfce6caae9e080ec4fbbecec0451a9b542b Mon Sep 17 00:00:00 2001
From: Marija Boberg <m.boberg@uke.de>
Date: Tue, 2 May 2023 11:57:55 +0200
Subject: [PATCH 01/14] improved patch handling

---
 src/MPISphericalHarmonics.jl | 16 +++++++++++++---
 test/runtests.jl             |  9 ++++++++-
 2 files changed, 21 insertions(+), 4 deletions(-)

diff --git a/src/MPISphericalHarmonics.jl b/src/MPISphericalHarmonics.jl
index d5589ec..aaf8c4c 100644
--- a/src/MPISphericalHarmonics.jl
+++ b/src/MPISphericalHarmonics.jl
@@ -8,10 +8,10 @@ using MPIMagneticFields
 using SphericalHarmonicExpansions
 using MPIFiles
 
-import Base.write
+import Base.write, Base.length
 
 export MagneticFieldCoefficients
-export selectPatch
+export selectPatch, length
 
 # Spherical harmonic coefficients describing a magnetic field
 mutable struct MagneticFieldCoefficients
@@ -250,8 +250,18 @@ MPIMagneticFields.fieldType(::SphericalHarmonicsDefinedField) = OtherField()
 MPIMagneticFields.definitionType(::SphericalHarmonicsDefinedField) = SphericalHarmonicsDataBasedFieldDefinition()
 MPIMagneticFields.timeDependencyType(::SphericalHarmonicsDefinedField) = TimeConstant()
 
+# get field values
 MPIMagneticFields.value(field::SphericalHarmonicsDefinedField, r::PT) where {T <: Number, PT <: AbstractVector{T}} = [field.func[i, field.patch].(r...) for i=1:3]
 
-selectPatch(field::SphericalHarmonicsDefinedField, patchNum) = field.patch = patchNum
+# patches
+length(field::SphericalHarmonicsDefinedField) = size(field.func,2) # get number of patches
+function selectPatch(field::SphericalHarmonicsDefinedField, patchNum::Int) 
+  if 1 <= patchNum <= length(field) # test if patch exists
+    field.patch = patchNum
+  else
+    throw(DimensionMismatch("The field contains only $(length(field)) patches."))
+  end
+end
+
 
 end
diff --git a/test/runtests.jl b/test/runtests.jl
index 8bac765..025e63f 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -158,13 +158,17 @@ using Aqua
       ## Multi-patch setting: Second field with offset
       coeffsPatch = hcat(deepcopy(coeffs),deepcopy(coeffs)) # two patches
       for j=1:3 coeffsPatch[j,2][0,0] = 0.01 end # set offset
-      field = SphericalHarmonicsDefinedField(coeffsPatch) # test constructor on coefficients
+      coeffsPatch_MF = MagneticFieldCoefficients(coeffsPatch) # create MagneticFieldCoefficients
+      field = SphericalHarmonicsDefinedField(coeffsPatch_MF) # test constructor on coefficients
 
       # Test field types
       @test fieldType(field) isa OtherField
       @test definitionType(field) isa SphericalHarmonicsDataBasedFieldDefinition
       @test timeDependencyType(field) isa TimeConstant
 
+      # Test number of patches
+      @test length(field) == 2
+
       ## Test FFPs (for both patches)
       # First patch
       @test field.patch == 1
@@ -180,6 +184,9 @@ using Aqua
       # test FFP
       ffp = [0.01, 0.01, -0.005]
       @test isapprox(field[ffp...], zeros(3), atol=1e-10)
+
+      # Test wrong patch number
+      @test_throws DimensionMismatch selectPatch(field,3)
     end
   end
 

From e780760b04b23d1ac1cb1166bf2c6a8b21240287 Mon Sep 17 00:00:00 2001
From: Jonas Schumacher <github735@jonasschumacher.de>
Date: Tue, 16 May 2023 10:30:25 +0200
Subject: [PATCH 02/14] Update to MPIMagneticFields 0.0.4

---
 Project.toml                 |  4 ++--
 src/MPISphericalHarmonics.jl |  8 ++++----
 test/runtests.jl             | 16 +++++++++-------
 3 files changed, 15 insertions(+), 13 deletions(-)

diff --git a/Project.toml b/Project.toml
index 04f2e1c..303963a 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "MPISphericalHarmonics"
 uuid = "2527f7a4-0af4-4016-a944-036fbac19de9"
 authors = ["Marija Boberg <m.boberg@uke.de> and contributors"]
-version = "0.0.4"
+version = "0.0.5"
 
 [deps]
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
@@ -14,7 +14,7 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 [compat]
 DocStringExtensions = "0.8, 0.9"
 MPIFiles = "0.12"
-MPIMagneticFields = "0.0.2, 0.0.3"
+MPIMagneticFields = "0.0.4"
 SphericalHarmonicExpansions = "0.1"
 Unitful = "1.11, 1.12"
 julia = "1"
diff --git a/src/MPISphericalHarmonics.jl b/src/MPISphericalHarmonics.jl
index aaf8c4c..9adef14 100644
--- a/src/MPISphericalHarmonics.jl
+++ b/src/MPISphericalHarmonics.jl
@@ -246,12 +246,12 @@ end
 SphericalHarmonicsDefinedField(coeffs::Array{SphericalHarmonicCoefficients}) = SphericalHarmonicsDefinedField(func = fastfunc.(coeffs))
 SphericalHarmonicsDefinedField(coeffs_MF::MagneticFieldCoefficients) = SphericalHarmonicsDefinedField(coeffs_MF.coeffs)
 
-MPIMagneticFields.fieldType(::SphericalHarmonicsDefinedField) = OtherField()
-MPIMagneticFields.definitionType(::SphericalHarmonicsDefinedField) = SphericalHarmonicsDataBasedFieldDefinition()
-MPIMagneticFields.timeDependencyType(::SphericalHarmonicsDefinedField) = TimeConstant()
+MPIMagneticFields.FieldStyle(::SphericalHarmonicsDefinedField) = OtherField()
+MPIMagneticFields.FieldDefinitionStyle(::SphericalHarmonicsDefinedField) = SphericalHarmonicsDataBasedFieldDefinition()
+MPIMagneticFields.FieldTimeDependencyStyle(::SphericalHarmonicsDefinedField) = TimeConstant()
 
 # get field values
-MPIMagneticFields.value(field::SphericalHarmonicsDefinedField, r::PT) where {T <: Number, PT <: AbstractVector{T}} = [field.func[i, field.patch].(r...) for i=1:3]
+MPIMagneticFields.value_(field::SphericalHarmonicsDefinedField, r) = [field.func[i, field.patch].(r...) for i=1:3]
 
 # patches
 length(field::SphericalHarmonicsDefinedField) = size(field.func,2) # get number of patches
diff --git a/test/runtests.jl b/test/runtests.jl
index 025e63f..d20e4af 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -15,9 +15,9 @@ using Aqua
 
     ## Calculate some field values
     # load a spherical t-design
-    tDes = loadTDesign(6,26,42.0u"mm")
+    tDes = loadTDesign(6, 26, 42.0u"mm")
     # ideal gradient field with FFP
-    idealField = IdealFFP([-1,-1,2])
+    idealField = IdealFFP([-1, -1, 2])
     # get field values
     fieldValues = hcat([idealField[ustrip.(Unitful.m.(pos))...] for pos in tDes]...)
 
@@ -146,6 +146,8 @@ using Aqua
       @test isapprox(coeffsW.radius, 0.042, atol=ε) # radius
       @test isapprox(coeffsW.center, zeros(3), atol=ε) # center
       @test coeffsW.ffp == zeros(3,1) # FFP
+      coeffsW = nothing # This maybe masks an implementation error
+      GC.gc()
 
       # remove test files
       rm(filename2)
@@ -162,9 +164,9 @@ using Aqua
       field = SphericalHarmonicsDefinedField(coeffsPatch_MF) # test constructor on coefficients
 
       # Test field types
-      @test fieldType(field) isa OtherField
-      @test definitionType(field) isa SphericalHarmonicsDataBasedFieldDefinition
-      @test timeDependencyType(field) isa TimeConstant
+      @test FieldStyle(field) isa OtherField
+      @test FieldDefinitionStyle(field) isa SphericalHarmonicsDataBasedFieldDefinition
+      @test FieldTimeDependencyStyle(field) isa TimeConstant
 
       # Test number of patches
       @test length(field) == 2
@@ -178,9 +180,9 @@ using Aqua
       # Second patch
       selectPatch(field,2)   
       @test field.patch == 2
-      # test offset in (0,0,0)
+      # test offset in (0, 0, 0)
       offset = ones(3) .* 0.01
-      @test isapprox(field[0,0,0], offset, atol=1e-10)
+      @test isapprox(field[0, 0, 0], offset, atol=1e-10)
       # test FFP
       ffp = [0.01, 0.01, -0.005]
       @test isapprox(field[ffp...], zeros(3), atol=1e-10)

From cd21c08e06156606fd09ed7472a11e76f1733235 Mon Sep 17 00:00:00 2001
From: Marija Boberg <m.boberg@uke.de>
Date: Thu, 8 Jun 2023 11:54:04 +0200
Subject: [PATCH 03/14] apply fastfunc to coefficients

---
 src/MPISphericalHarmonics.jl | 15 ++++++---------
 test/runtests.jl             |  2 +-
 2 files changed, 7 insertions(+), 10 deletions(-)

diff --git a/src/MPISphericalHarmonics.jl b/src/MPISphericalHarmonics.jl
index 9adef14..5ca10f4 100644
--- a/src/MPISphericalHarmonics.jl
+++ b/src/MPISphericalHarmonics.jl
@@ -205,22 +205,20 @@ function loadTDesignCoefficients(filename::String)
   end
   @polyvar x y z
   tDes = MPIFiles.loadTDesign(Int(t),N,radius*u"m", center.*u"m")
-  coeffs, expansion, func = magneticField(tDes, field, x,y,z)
+  coeffs, = magneticField(tDes, field, x,y,z)
   for c=1:size(coeffs,2), j = 1:3
     coeffs[j, c] = SphericalHarmonicExpansions.translation(coeffs[j, c], correction[:, j])
-    expansion[j, c] = sphericalHarmonicsExpansion(coeffs[j, c], x, y, z);
-    func[j, c] = @fastfunc expansion[j, c]
   end
 
   coeffs_MF = MagneticFieldCoefficients(coeffs, radius, center)
 
-  return coeffs_MF, expansion, func
+  return coeffs_MF
 end
 
 ## SphericalHarmonicsDefinedField ##
 export SphericalHarmonicsDefinedField
 Base.@kwdef mutable struct SphericalHarmonicsDefinedField <: AbstractMagneticField
-  func::Array{Function, 2}
+  func::Array{SphericalHarmonicExpansions.StaticPolynomials.Polynomial, 2}
   patch::Integer = 1
 end
 
@@ -230,13 +228,12 @@ function SphericalHarmonicsDefinedField(filename::String)
     if haskey(file,"coeffs") 
       # load coefficients
       coeffs_MF = MagneticFieldCoefficients(filename)
-      func = fastfunc.(coeffs_MF.coeffs)
     else
       # load measured field
-      coeffs_MF, expansion, func = loadTDesignCoefficients(filename)
+      coeffs_MF = loadTDesignCoefficients(filename)
     end
 
-    return func
+    return fastfunc.(coeffs_MF.coeffs)
   end
 
   return SphericalHarmonicsDefinedField(func=func)
@@ -251,7 +248,7 @@ MPIMagneticFields.FieldDefinitionStyle(::SphericalHarmonicsDefinedField) = Spher
 MPIMagneticFields.FieldTimeDependencyStyle(::SphericalHarmonicsDefinedField) = TimeConstant()
 
 # get field values
-MPIMagneticFields.value_(field::SphericalHarmonicsDefinedField, r) = [field.func[i, field.patch].(r...) for i=1:3]
+MPIMagneticFields.value_(field::SphericalHarmonicsDefinedField, r) = [field.func[i, field.patch](r) for i=1:3]
 
 # patches
 length(field::SphericalHarmonicsDefinedField) = size(field.func,2) # get number of patches
diff --git a/test/runtests.jl b/test/runtests.jl
index d20e4af..3796bb2 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -99,7 +99,7 @@ using Aqua
       @test isapprox(field[0.01,0.01,0.01], [-0.01,-0.01,0.02], atol=ε)
 
       # get coefficients
-      coeffsMF, = MPISphericalHarmonics.loadTDesignCoefficients(filename)
+      coeffsMF = MPISphericalHarmonics.loadTDesignCoefficients(filename)
       @test isapprox(coeffsMF.radius, 0.042, atol=ε) # radius
       @test isapprox(coeffsMF.coeffs[1][1,1], -1.0, atol=1e-10) # gradient (x)
       @test isapprox(coeffsMF.coeffs[2][1,-1], -1.0, atol=1e-10) # gradient (y)

From 2085987fb44678edea7b99cf668e459d69f3f66b Mon Sep 17 00:00:00 2001
From: Jonas Schumacher <github735@jonasschumacher.de>
Date: Tue, 13 Jun 2023 11:42:20 +0200
Subject: [PATCH 04/14] Bump version compatability for `MPIMagneticFields`

---
 Project.toml | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/Project.toml b/Project.toml
index 303963a..f8fe85c 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "MPISphericalHarmonics"
 uuid = "2527f7a4-0af4-4016-a944-036fbac19de9"
 authors = ["Marija Boberg <m.boberg@uke.de> and contributors"]
-version = "0.0.5"
+version = "0.0.6"
 
 [deps]
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
@@ -14,7 +14,7 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 [compat]
 DocStringExtensions = "0.8, 0.9"
 MPIFiles = "0.12"
-MPIMagneticFields = "0.0.4"
+MPIMagneticFields = "0.0.4, 0.0.5"
 SphericalHarmonicExpansions = "0.1"
 Unitful = "1.11, 1.12"
 julia = "1"

From b6112fa563d35f9c0ba4218d730242a188414247 Mon Sep 17 00:00:00 2001
From: Marija Boberg <m.boberg@uke.de>
Date: Wed, 14 Jun 2023 16:34:54 +0200
Subject: [PATCH 05/14] added mathematical operators for coefficients

---
 src/MPISphericalHarmonics.jl     | 209 +-----------------------
 src/MagneticFieldCoefficients.jl | 269 +++++++++++++++++++++++++++++++
 test/runtests.jl                 |  44 +++++
 3 files changed, 318 insertions(+), 204 deletions(-)
 create mode 100644 src/MagneticFieldCoefficients.jl

diff --git a/src/MPISphericalHarmonics.jl b/src/MPISphericalHarmonics.jl
index 5ca10f4..ab38652 100644
--- a/src/MPISphericalHarmonics.jl
+++ b/src/MPISphericalHarmonics.jl
@@ -8,215 +8,16 @@ using MPIMagneticFields
 using SphericalHarmonicExpansions
 using MPIFiles
 
-import Base.write, Base.length
+import Base.length
 
+# load MagneticFieldCoefficients
+include("MagneticFieldCoefficients.jl")
 export MagneticFieldCoefficients
-export selectPatch, length
-
-# Spherical harmonic coefficients describing a magnetic field
-mutable struct MagneticFieldCoefficients
-  coeffs::Array{SphericalHarmonicCoefficients,2} # coefficients
-  radius::Float64 # radius of measured sphere
-  center::Vector{Float64} # center of measured sphere
-  ffp::Union{Array{Float64,2},Nothing} # field-free-point (if available)
-
-  function MagneticFieldCoefficients(coeffs::Array{SphericalHarmonicCoefficients,2}, radius::Float64,
-  				     center::Vector{Float64}, ffp::Union{Array{Float64,2},Nothing})
-    # test sizes of the arrays
-    if size(coeffs,1) != 3
-      throw(DimensionMismatch("The coefficient matrix needs 3 entries (x,y,z) in the first dimension, not $(size(coeffs,1))"))
-    elseif ffp !== nothing
-      if size(ffp,1) != 3
-        throw(DimensionMismatch("The FFP matrix needs 3 entries (x,y,z) in the first dimension, not $(size(coeffs,1))"))
-      elseif size(coeffs,2) != size(ffp,2)
-        throw(DimensionMismatch("The number of patches of the coefficients and FFPs does not match: $(size(coeffs,2)) != $(size(ffp,2))"))
-      end
-    end
-
-    return new(coeffs,radius,center,ffp)
-  end
-
-end
-
-# some other constructors
-MagneticFieldCoefficients(coeffs::Array{SphericalHarmonicCoefficients,2}) = MagneticFieldCoefficients(coeffs, 0.0)
-MagneticFieldCoefficients(coeffs::Array{SphericalHarmonicCoefficients,2}, radius::Float64) = MagneticFieldCoefficients(coeffs,radius,[0.0,0.0,0.0])
-MagneticFieldCoefficients(coeffs::Array{SphericalHarmonicCoefficients,2}, radius::Float64, center::Vector{Float64}) = MagneticFieldCoefficients(coeffs,radius,center,nothing)
-MagneticFieldCoefficients(coeffs::Array{SphericalHarmonicCoefficients,2}, radius::Float64, ffp::Array{Float64,2}) = MagneticFieldCoefficients(coeffs,radius,[0.0,0.0,0.0],ffp)
-
-function MagneticFieldCoefficients(L::Int)
-  if L<0
-    throw(DomainError(L,"Input vector needs to be of size (L+1)², where L ∈ ℕ₀."))
-  end
-  return MagneticFieldCoefficients(reshape([SphericalHarmonicCoefficients(L),SphericalHarmonicCoefficients(L),SphericalHarmonicCoefficients(L)],3,1))
-end
-
-# constructor using t-design
-MagneticFieldCoefficients(coeffs::Array{SphericalHarmonicCoefficients,2}, tDesign::SphericalTDesign, ffp=nothing) = MagneticFieldCoefficients(coeffs,ustrip(Unitful.m(tDesign.radius)), ustrip.(Unitful.m.(tDesign.center)), ffp)
-  
-  
-# read coefficients from an HDF5-file
-function MagneticFieldCoefficients(path::String)
-
-  # load spherical harmonic coefficients
-  shcoeffs = SphericalHarmonicCoefficients(path)
-  
-  coeffsMF = h5open(path,"r") do file
-    if haskey(HDF5.root(file), "/radius") 
-      # file contains all relevant information
-      radius = read(file, "/radius")
-      center = read(file, "/center")
-      if haskey(HDF5.root(file), "/ffp")
-        ffp = read(file, "/ffp")
-        return MagneticFieldCoefficients(shcoeffs, radius, center, ffp)
-      else
-        # field has not FFP -> ffp = nothing
-        return MagneticFieldCoefficients(shcoeffs, radius, center)
-      end
-    else
-      # convert file of SphericalHarmonicCoefficients into MagneticFieldCoefficients
-      # -> set all additional informations to 0 or nothing
-      # use radius = 0.01 as default value
-      return MagneticFieldCoefficients(shcoeffs, 0.01)
-    end
-  end
-
-  return coeffsMF
-end
-
-# write coefficients to an HDF5-file
-function write(path::String, coeffs::MagneticFieldCoefficients)
-
-  # save SphericalHarmonicCoefficients
-  write(path,coeffs.coeffs)
-  
-  # add field informations
-  radius = coeffs.radius
-  center = coeffs.center
-  ffp = coeffs.ffp
-
-  h5open(path,"cw") do file
-      write(file, "/radius", radius)
-      write(file, "/center", center)
-      if ffp !== nothing
-        write(file, "/ffp", ffp)
-      end
-  end
-end
-
-"""
-   magneticField(tDesign::SphericalTDesign, field::Union{AbstractArray{T,2},AbstractArray{T,3}}, 
-		       x::Variable, y::Variable, z::Variable;
-		       L::Int=Int(tDesign.T/2),
-		       calcSolid::Bool=true) where T <: Real
-*Description:*  Calculation of the spherical harmonic coefficients and expansion based on the measured t-design\\
- \\
-*Input:*
-- `tDesign`	- Measured t-design (type: SphericalTDesign)
-- `field`       - Measured field (size = (J,N,C)) with J <= 3
-- `x, y, z`     - Cartesian coordinates
-**kwargs:**
-- `L`           - Order up to which the coeffs be calculated (default: t/2)
-- `calcSolid`   - Boolean (default: true)\\
-    false -> spherical coefficients\\
-    true -> solid coefficients
-*Output:*
-- `coeffs`    - spherical/solid coefficients, type: Array{SphericalHarmonicCoefficients}(3,C)
-- `expansion` - related expansion (Cartesian polynomial), type: Array{AbstractPolynomialLike}(3,C)
-- `func`      - expansion converted to a function, type: Array{Function}(3,C)
-"""
-function magneticField(tDesign::SphericalTDesign, field::Union{AbstractArray{T,2},AbstractArray{T,3}}, 
-		       x::Variable, y::Variable, z::Variable;
-		       L::Int=floor(Int,tDesign.T/2),
-		       calcSolid::Bool=true) where T <: Real
-
-  # get tDesign positions [m] and removing the unit
-  # coordinates
-  coords = Float64.(ustrip.(Unitful.m.(hcat([p for p in tDesign]...))))
-
-  # radius
-  R = Float64(ustrip(Unitful.m(tDesign.radius)))
-
-  # center
-  center = Float64.(ustrip.(Unitful.m.(tDesign.center)))
-  
-  return magneticField(coords, field, 
-		       R, center, L,
-		       x, y, z,
-		       calcSolid)
-end
-
-function magneticField(coords::AbstractArray{T, 2}, field::Union{AbstractArray{T, 2}, AbstractArray{T, 3}}, 
-                       R::T, center::Vector{T}, L::Int,
-                       x::Variable, y::Variable, z::Variable, 
-                       calcSolid::Bool=true) where T <: Real
-
-  # transpose coords if its dimensions do not fit
-  if size(coords,1) != 3
-    coords = coords'
-  end
-
-  # test dimensions of field array
-  if size(field,1) > 3
-    throw(DimensionMismatch("The measured field has more than 3 entries in the first dimension: $(size(field,1))"))
-  elseif size(field,2) != size(coords,2)
-    throw(DimensionMismatch("The field vector does not match the size of the tdesign: $(size(field,2)) != $(size(coords,2))"))
-  end
-
-  func= Array{Function}(undef,size(field,1),size(field,3))
-  expansion = Array{Polynomial}(undef,size(field,1),size(field,3))
-  coeffs = Array{SphericalHarmonicCoefficients}(undef,size(field,1),size(field,3))
-
-  # rescale coordinates to t-design on unit sphere
-  coords = coords .- center
-  coords *= 1/R
-  for c in axes(field,3)
-    # calculation of the coefficients
-    for j in axes(field,1)
-      coeffs[j,c] = SphericalHarmonicExpansions.sphericalQuadrature(field[j,:,c],coords',L);
-      coeffs[j,c].R = R
-
-      normalize!(coeffs[j,c],R)
-
-      # convert spherical into solid coefficients
-      if calcSolid
-          solid!(coeffs[j,c])
-      end
-
-      # calculation of the expansion
-      expansion[j,c] = sphericalHarmonicsExpansion(coeffs[j,c],x,y,z) + 0*x;
-      func[j,c] = @fastfunc expansion[j,c]+0*x+0*y+0*z
-    end
-  end
-
-  return coeffs, expansion, func
-end
-
-#TODO: This should be merged with and moved to MPIFiles
-function loadTDesignCoefficients(filename::String)
-  field, radius, N, t, center, correction  = h5open(filename, "r") do file
-    field = read(file,"/fields") 		# measured field (size: 3 x #points x #patches)
-    radius = read(file,"/positions/tDesign/radius")	# radius of the measured ball
-    N = read(file,"/positions/tDesign/N")		# number of points of the t-design
-    t = read(file,"/positions/tDesign/t")		# t of the t-design
-    center = read(file,"/positions/tDesign/center")	# center of the measured ball
-    correction = read(file, "/sensor/correctionTranslation")
-    return field, radius, N, t, center, correction
-  end
-  @polyvar x y z
-  tDes = MPIFiles.loadTDesign(Int(t),N,radius*u"m", center.*u"m")
-  coeffs, = magneticField(tDes, field, x,y,z)
-  for c=1:size(coeffs,2), j = 1:3
-    coeffs[j, c] = SphericalHarmonicExpansions.translation(coeffs[j, c], correction[:, j])
-  end
-
-  coeffs_MF = MagneticFieldCoefficients(coeffs, radius, center)
-
-  return coeffs_MF
-end
 
 ## SphericalHarmonicsDefinedField ##
 export SphericalHarmonicsDefinedField
+export selectPatch, length
+
 Base.@kwdef mutable struct SphericalHarmonicsDefinedField <: AbstractMagneticField
   func::Array{SphericalHarmonicExpansions.StaticPolynomials.Polynomial, 2}
   patch::Integer = 1
diff --git a/src/MagneticFieldCoefficients.jl b/src/MagneticFieldCoefficients.jl
new file mode 100644
index 0000000..7e5bbfe
--- /dev/null
+++ b/src/MagneticFieldCoefficients.jl
@@ -0,0 +1,269 @@
+import Base.isapprox, Base.==,
+        Base.+, Base.-, Base.*, Base./, 
+        Base.write, Base.size
+
+# Spherical harmonic coefficients describing a magnetic field
+mutable struct MagneticFieldCoefficients
+    coeffs::Array{SphericalHarmonicCoefficients,2} # coefficients
+    radius::Float64 # radius of measured sphere
+    center::Vector{Float64} # center of measured sphere
+    ffp::Union{Array{Float64,2},Nothing} # field-free-point (if available)
+
+    function MagneticFieldCoefficients(coeffs::Array{SphericalHarmonicCoefficients,2}, radius::Float64,
+                                        center::Vector{Float64}, ffp::Union{Array{Float64,2},Nothing})
+        # test sizes of the arrays
+        if size(coeffs,1) != 3
+            throw(DimensionMismatch("The coefficient matrix needs 3 entries (x,y,z) in the first dimension, not $(size(coeffs,1))"))
+        elseif ffp !== nothing
+            if size(ffp,1) != 3
+                throw(DimensionMismatch("The FFP matrix needs 3 entries (x,y,z) in the first dimension, not $(size(coeffs,1))"))
+            elseif size(coeffs,2) != size(ffp,2)
+                throw(DimensionMismatch("The number of patches of the coefficients and FFPs does not match: $(size(coeffs,2)) != $(size(ffp,2))"))
+            end
+        end
+
+        return new(coeffs,radius,center,ffp)
+    end
+
+end
+
+# some other constructors
+MagneticFieldCoefficients(coeffs::Array{SphericalHarmonicCoefficients,2}) = MagneticFieldCoefficients(coeffs, 0.0)
+MagneticFieldCoefficients(coeffs::Array{SphericalHarmonicCoefficients,2}, radius::Float64) = MagneticFieldCoefficients(coeffs,radius,[0.0,0.0,0.0])
+MagneticFieldCoefficients(coeffs::Array{SphericalHarmonicCoefficients,2}, radius::Float64, center::Vector{Float64}) = MagneticFieldCoefficients(coeffs,radius,center,nothing)
+MagneticFieldCoefficients(coeffs::Array{SphericalHarmonicCoefficients,2}, radius::Float64, ffp::Array{Float64,2}) = MagneticFieldCoefficients(coeffs,radius,[0.0,0.0,0.0],ffp)
+
+function MagneticFieldCoefficients(L::Int)
+    if L<0
+        throw(DomainError(L,"Input vector needs to be of size (L+1)², where L ∈ ℕ₀."))
+    end
+    return MagneticFieldCoefficients(fill(SphericalHarmonicCoefficients(L), (3,1)))
+end
+
+# constructor using t-design
+MagneticFieldCoefficients(coeffs::Array{SphericalHarmonicCoefficients,2}, tDesign::SphericalTDesign, ffp=nothing) = MagneticFieldCoefficients(coeffs,ustrip(Unitful.m(tDesign.radius)), ustrip.(Unitful.m.(tDesign.center)), ffp)
+
+
+# read coefficients from an HDF5-file
+function MagneticFieldCoefficients(path::String)
+
+    # load spherical harmonic coefficients
+    shcoeffs = SphericalHarmonicCoefficients(path)
+
+    coeffsMF = h5open(path,"r") do file
+        if haskey(HDF5.root(file), "/radius") 
+            # file contains all relevant information
+            radius = read(file, "/radius")
+            center = read(file, "/center")
+            if haskey(HDF5.root(file), "/ffp")
+                ffp = read(file, "/ffp")
+                return MagneticFieldCoefficients(shcoeffs, radius, center, ffp)
+            else
+                # field has not FFP -> ffp = nothing
+                return MagneticFieldCoefficients(shcoeffs, radius, center)
+            end
+        else
+            # convert file of SphericalHarmonicCoefficients into MagneticFieldCoefficients
+            # -> set all additional informations to 0 or nothing
+            # use radius = 0.01 as default value
+            return MagneticFieldCoefficients(shcoeffs, 0.01)
+        end
+    end
+
+    return coeffsMF
+end
+
+# write coefficients to an HDF5-file
+function write(path::String, coeffs::MagneticFieldCoefficients)
+
+    # save SphericalHarmonicCoefficients
+    write(path,coeffs.coeffs)
+
+    # add field informations
+    radius = coeffs.radius
+    center = coeffs.center
+    ffp = coeffs.ffp
+
+    h5open(path,"cw") do file
+        write(file, "/radius", radius)
+        write(file, "/center", center)
+        if ffp !== nothing
+            write(file, "/ffp", ffp)
+        end
+    end
+end
+
+# Size
+size(mfc::MagneticFieldCoefficients) = size(mfc.coeffs)
+
+# Operations on MagneticFieldCoefficients
+function isapprox(mfc1::MagneticFieldCoefficients, mfc2::MagneticFieldCoefficients; kargs...)
+    val = all(isapprox.(mfc1.coeffs,mfc2.coeffs;kargs...)) && isapprox(mfc1.radius,mfc2.radius;kargs...) && isapprox(mfc1.center,mfc2.center;kargs...)
+    if mfc1.ffp === nothing && mfc2.ffp === nothing
+        return  val
+    elseif mfc1.ffp !== nothing || mfc2.ffp !== nothing
+        @info "Only one of the coefficients has FFPs. Applying isapprox to the other values yields $val."
+        return false
+    else
+        return val && isapprox(mfc1.ffp,mfc2.ffp;kargs...)
+    end
+end
+==(mfc1::MagneticFieldCoefficients, mfc2::MagneticFieldCoefficients) = 
+    mfc1.coeffs == mfc2.coeffs && mfc1.radius == mfc2.radius && mfc1.center == mfc2.center && mfc1.ffp == mfc2.ffp
+
+"""
+    +(mfc1::MagneticFieldCoefficients, mfc2::MagneticFieldCoefficients; force::Bool=false)
+
+`force = true` adds the coefficients even if the radius or center are not equal (set to values of the first coefficients).
+"""
+function +(mfc1::MagneticFieldCoefficients, mfc2::MagneticFieldCoefficients; force::Bool=false)
+    if force 
+        return MagneticFieldCoefficients(mfc1.coeffs .+ mfc2.coeffs, mfc1.radius, mfc1.center)
+    end
+    if mfc1.radius != mfc2.radius
+        throw(DomainError([mfc1.radius,mfc2.radius],"Coefficients do not have the same measurement radius."))
+    end
+    if mfc1.center != mfc2.center
+        throw(DomainError([mfc1.center,mfc2.center],"Coefficients do not have the same measurement center."))
+    end
+    return MagneticFieldCoefficients(mfc1.coeffs .+ mfc2.coeffs, mfc1.radius, mfc1.center)
+end
+
+"""
+    -(mfc1::MagneticFieldCoefficients, mfc2::MagneticFieldCoefficients; force::Bool=false)
+
+`force = true` subtracts the coefficients even if the radius or center are not equal (set to values of the first coefficients).
+"""
+function -(mfc1::MagneticFieldCoefficients, mfc2::MagneticFieldCoefficients; force::Bool=false)
+    if force 
+        return MagneticFieldCoefficients(mfc1.coeffs .- mfc2.coeffs, mfc1.radius, mfc1.center)
+    end
+    if mfc1.radius != mfc2.radius
+        throw(DomainError([mfc1.radius,mfc2.radius],"Coefficients do not have the same measurement radius."))
+    end
+    if mfc1.center != mfc2.center
+        throw(DomainError([mfc1.center,mfc2.center],"Coefficients do not have the same measurement center."))
+    end
+    return MagneticFieldCoefficients(mfc1.coeffs .- mfc2.coeffs, mfc1.radius, mfc1.center)
+end
+
++(mfc::MagneticFieldCoefficients, value::Number) = MagneticFieldCoefficients(mfc.coeffs .+ value, mfc.radius, mfc.center, mfc.ffp)
+-(mfc::MagneticFieldCoefficients, value::Number) = MagneticFieldCoefficients(mfc.coeffs .- value, mfc.radius, mfc.center, mfc.ffp)
+*(mfc::MagneticFieldCoefficients, value::Number) = MagneticFieldCoefficients(mfc.coeffs .* value, mfc.radius, mfc.center, mfc.ffp)
+/(mfc::MagneticFieldCoefficients, value::Number) = MagneticFieldCoefficients(mfc.coeffs ./ value, mfc.radius, mfc.center, mfc.ffp)
++(value::Number, mfc::MagneticFieldCoefficients) = +(mfc::MagneticFieldCoefficients, value)
+-(value::Number, mfc::MagneticFieldCoefficients) = MagneticFieldCoefficients(value .- mfc.coeffs, mfc.radius, mfc.center, mfc.ffp)
+*(value::Number, mfc::MagneticFieldCoefficients) = *(mfc::MagneticFieldCoefficients, value)
+
+
+## Load coefficients from t-design measurement ##
+"""
+    magneticField(tDesign::SphericalTDesign, field::Union{AbstractArray{T,2},AbstractArray{T,3}}, 
+		       x::Variable, y::Variable, z::Variable;
+		       L::Int=Int(tDesign.T/2),
+		       calcSolid::Bool=true) where T <: Real
+*Description:*  Calculation of the spherical harmonic coefficients and expansion based on the measured t-design\\
+ \\
+*Input:*
+- `tDesign`	- Measured t-design (type: SphericalTDesign)
+- `field`       - Measured field (size = (J,N,C)) with J <= 3
+- `x, y, z`     - Cartesian coordinates
+**kwargs:**
+- `L`           - Order up to which the coeffs be calculated (default: t/2)
+- `calcSolid`   - Boolean (default: true)\\
+    false -> spherical coefficients\\
+    true -> solid coefficients
+*Output:*
+- `coeffs`    - spherical/solid coefficients, type: Array{SphericalHarmonicCoefficients}(3,C)
+- `expansion` - related expansion (Cartesian polynomial), type: Array{AbstractPolynomialLike}(3,C)
+- `func`      - expansion converted to a function, type: Array{Function}(3,C)
+"""
+function magneticField(tDesign::SphericalTDesign, field::Union{AbstractArray{T,2},AbstractArray{T,3}}, 
+		       x::Variable, y::Variable, z::Variable;
+		       L::Int=floor(Int,tDesign.T/2),
+		       calcSolid::Bool=true) where T <: Real
+
+  # get tDesign positions [m] and removing the unit
+  # coordinates
+  coords = Float64.(ustrip.(Unitful.m.(hcat([p for p in tDesign]...))))
+
+  # radius
+  R = Float64(ustrip(Unitful.m(tDesign.radius)))
+
+  # center
+  center = Float64.(ustrip.(Unitful.m.(tDesign.center)))
+  
+  return magneticField(coords, field, 
+		       R, center, L,
+		       x, y, z,
+		       calcSolid)
+end
+
+function magneticField(coords::AbstractArray{T, 2}, field::Union{AbstractArray{T, 2}, AbstractArray{T, 3}}, 
+                       R::T, center::Vector{T}, L::Int,
+                       x::Variable, y::Variable, z::Variable, 
+                       calcSolid::Bool=true) where T <: Real
+
+  # transpose coords if its dimensions do not fit
+  if size(coords,1) != 3
+    coords = coords'
+  end
+
+  # test dimensions of field array
+  if size(field,1) > 3
+    throw(DimensionMismatch("The measured field has more than 3 entries in the first dimension: $(size(field,1))"))
+  elseif size(field,2) != size(coords,2)
+    throw(DimensionMismatch("The field vector does not match the size of the tdesign: $(size(field,2)) != $(size(coords,2))"))
+  end
+
+  func= Array{Function}(undef,size(field,1),size(field,3))
+  expansion = Array{Polynomial}(undef,size(field,1),size(field,3))
+  coeffs = Array{SphericalHarmonicCoefficients}(undef,size(field,1),size(field,3))
+
+  # rescale coordinates to t-design on unit sphere
+  coords = coords .- center
+  coords *= 1/R
+  for c in axes(field,3)
+    # calculation of the coefficients
+    for j in axes(field,1)
+      coeffs[j,c] = SphericalHarmonicExpansions.sphericalQuadrature(field[j,:,c],coords',L);
+      coeffs[j,c].R = R
+
+      normalize!(coeffs[j,c],R)
+
+      # convert spherical into solid coefficients
+      if calcSolid
+          solid!(coeffs[j,c])
+      end
+
+      # calculation of the expansion
+      expansion[j,c] = sphericalHarmonicsExpansion(coeffs[j,c],x,y,z) + 0*x;
+      func[j,c] = @fastfunc expansion[j,c]+0*x+0*y+0*z
+    end
+  end
+
+  return coeffs, expansion, func
+end
+
+#TODO: This should be merged with and moved to MPIFiles
+function loadTDesignCoefficients(filename::String)
+  field, radius, N, t, center, correction  = h5open(filename, "r") do file
+    field = read(file,"/fields") 		# measured field (size: 3 x #points x #patches)
+    radius = read(file,"/positions/tDesign/radius")	# radius of the measured ball
+    N = read(file,"/positions/tDesign/N")		# number of points of the t-design
+    t = read(file,"/positions/tDesign/t")		# t of the t-design
+    center = read(file,"/positions/tDesign/center")	# center of the measured ball
+    correction = read(file, "/sensor/correctionTranslation")
+    return field, radius, N, t, center, correction
+  end
+  @polyvar x y z
+  tDes = MPIFiles.loadTDesign(Int(t),N,radius*u"m", center.*u"m")
+  coeffs, = magneticField(tDes, field, x,y,z)
+  for c=1:size(coeffs,2), j = 1:3
+    coeffs[j, c] = SphericalHarmonicExpansions.translation(coeffs[j, c], correction[:, j])
+  end
+
+  coeffs_MF = MagneticFieldCoefficients(coeffs, radius, center)
+
+  return coeffs_MF
+end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index 3796bb2..b084d55 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -88,6 +88,50 @@ using Aqua
       # constructor with wrong FFP sizes
       @test_throws DimensionMismatch MagneticFieldCoefficients(coeffs, tDes, zeros(2,1))
       @test_throws DimensionMismatch MagneticFieldCoefficients(coeffs, tDes, zeros(3,2))
+
+      # test size
+      @test size(coeffsMF) == (3,1)
+    end
+
+    @testset "MagneticFieldCoefficient operations" begin
+      # coefficients for the tests
+      shc0 = fill(SphericalHarmonicCoefficients(zeros(4)), (3,1))
+      shc1 = fill(SphericalHarmonicCoefficients(ones(4)), (3,1))
+      shc2 = 2 .* shc1
+      c0 = MagneticFieldCoefficients(shc0,0.042,zeros(3)) 
+      c1 = MagneticFieldCoefficients(shc1,0.042,zeros(3)) 
+      c2 = MagneticFieldCoefficients(shc2,0.042,zeros(3)) 
+      c1R = MagneticFieldCoefficients(shc1,0.01,zeros(3)) # different radius
+      c1C = MagneticFieldCoefficients(shc1,0.042,ones(3)) # different center
+      c1F = MagneticFieldCoefficients(shc1,0.042,zeros(3),zeros(3,1)) # FFP 1
+      c1F2 = MagneticFieldCoefficients(shc1,0.042,zeros(3),ones(3,1)) # FFP 2
+
+      # isapprox 
+      @test !isapprox(c1, c2) # wrong coefficients
+      @test !isapprox(c1, c1R) # wrong radius
+      @test !isapprox(c1, c1C) # wrong center
+      @test !isapprox(c1, c1F) # one FFP
+      @test !isapprox(c1F, c1F2) # wrong FFP
+
+      # addition/subtraction # @test_throws DomainError mfc1+mfc2
+      @test_throws DomainError c1 + c1R # DomainError radius
+      @test_throws DomainError c2 - c1R # DomainError radius
+      @test_throws DomainError c1 + c1C # DomainError center
+      @test_throws DomainError c2 - c1C # DomainError center
+      @test isapprox(+(c1,c1R,force=true), c2) # force = true 
+      @test isapprox(-(c2,c1R,force=true), c1) # force = true 
+      @test isapprox(c1 + c1, c2) # correct
+      @test isapprox(c2 - c1, c1) # correct
+
+      # operations MFC and value
+      @test isapprox(1 + c1, c2)
+      @test isapprox(c0 + 2, c2)
+      @test isapprox(1 - c1, c0)
+      @test isapprox(c2 - 2, c0)
+      @test isapprox(0 * c1, c0)
+      @test isapprox(c1 * 2, c2)
+      @test isapprox(c2 / 2, c1)
+
     end
 
     @testset "Load/write data from/to file" begin

From eef1daf52c4a03af99794c7e21e5fdcfe0a07a00 Mon Sep 17 00:00:00 2001
From: Marija Boberg <m.boberg@uke.de>
Date: Fri, 7 Jul 2023 20:52:15 +0200
Subject: [PATCH 06/14] added shifting and FFP functions

---
 Project.toml                     |   1 +
 src/MPISphericalHarmonics.jl     |   2 +
 src/MagneticFieldCoefficients.jl | 155 ++++++++++++++++++++++++++++---
 test/runtests.jl                 |  10 +-
 4 files changed, 152 insertions(+), 16 deletions(-)

diff --git a/Project.toml b/Project.toml
index f8fe85c..f51020b 100644
--- a/Project.toml
+++ b/Project.toml
@@ -8,6 +8,7 @@ DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MPIFiles = "371237a9-e6c1-5201-9adb-3d8cfa78fa9f"
 MPIMagneticFields = "f6dda52a-86e9-4b50-afb7-f39836a99446"
+NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
 SphericalHarmonicExpansions = "504afa71-bae7-47b4-8ec9-3851161806ac"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
diff --git a/src/MPISphericalHarmonics.jl b/src/MPISphericalHarmonics.jl
index ab38652..732367e 100644
--- a/src/MPISphericalHarmonics.jl
+++ b/src/MPISphericalHarmonics.jl
@@ -7,12 +7,14 @@ using Unitful
 using MPIMagneticFields
 using SphericalHarmonicExpansions
 using MPIFiles
+using NLsolve # for Newton solver
 
 import Base.length
 
 # load MagneticFieldCoefficients
 include("MagneticFieldCoefficients.jl")
 export MagneticFieldCoefficients
+export shift, shift!, shiftFFP!
 
 ## SphericalHarmonicsDefinedField ##
 export SphericalHarmonicsDefinedField
diff --git a/src/MagneticFieldCoefficients.jl b/src/MagneticFieldCoefficients.jl
index 7e5bbfe..3ba39bd 100644
--- a/src/MagneticFieldCoefficients.jl
+++ b/src/MagneticFieldCoefficients.jl
@@ -6,15 +6,17 @@ import Base.isapprox, Base.==,
 mutable struct MagneticFieldCoefficients
     coeffs::Array{SphericalHarmonicCoefficients,2} # coefficients
     radius::Float64 # radius of measured sphere
-    center::Vector{Float64} # center of measured sphere
-    ffp::Union{Array{Float64,2},Nothing} # field-free-point (if available)
+    center::Array{Float64,2} # center of measured sphere (depends on expansion point)
+    ffp::Union{Array{Float64,2},Nothing} # field-free-point (if available, depends on expansion point)
 
     function MagneticFieldCoefficients(coeffs::Array{SphericalHarmonicCoefficients,2}, radius::Float64,
-                                        center::Vector{Float64}, ffp::Union{Array{Float64,2},Nothing})
+                                        center::Array{Float64,2}, ffp::Union{Array{Float64,2},Nothing})
         # test sizes of the arrays
         if size(coeffs,1) != 3
             throw(DimensionMismatch("The coefficient matrix needs 3 entries (x,y,z) in the first dimension, not $(size(coeffs,1))"))
-        elseif ffp !== nothing
+        elseif size(coeffs,2) != size(center,2)
+            throw(DimensionMismatch("The number of patches of the coefficients and center does not match: $(size(coeffs,2)) != $(size(center,2))"))
+        elseif !isnothing(ffp)
             if size(ffp,1) != 3
                 throw(DimensionMismatch("The FFP matrix needs 3 entries (x,y,z) in the first dimension, not $(size(coeffs,1))"))
             elseif size(coeffs,2) != size(ffp,2)
@@ -30,14 +32,14 @@ end
 # some other constructors
 MagneticFieldCoefficients(coeffs::Array{SphericalHarmonicCoefficients,2}) = MagneticFieldCoefficients(coeffs, 0.0)
 MagneticFieldCoefficients(coeffs::Array{SphericalHarmonicCoefficients,2}, radius::Float64) = MagneticFieldCoefficients(coeffs,radius,[0.0,0.0,0.0])
-MagneticFieldCoefficients(coeffs::Array{SphericalHarmonicCoefficients,2}, radius::Float64, center::Vector{Float64}) = MagneticFieldCoefficients(coeffs,radius,center,nothing)
-MagneticFieldCoefficients(coeffs::Array{SphericalHarmonicCoefficients,2}, radius::Float64, ffp::Array{Float64,2}) = MagneticFieldCoefficients(coeffs,radius,[0.0,0.0,0.0],ffp)
+MagneticFieldCoefficients(coeffs::Array{SphericalHarmonicCoefficients,2}, radius::Float64, center::VecOrMat{Float64}) = MagneticFieldCoefficients(coeffs,radius,center,nothing)
+MagneticFieldCoefficients(coeffs::Array{SphericalHarmonicCoefficients,2}, radius::Float64, center::Vector{Float64},ffp::Union{Matrix{Float64},Nothing}) = MagneticFieldCoefficients(coeffs,radius,hcat([center for p=1:size(coeffs,2)]...),ffp)
 
-function MagneticFieldCoefficients(L::Int)
+function MagneticFieldCoefficients(L::Int; numPatches::Int=1)
     if L<0
         throw(DomainError(L,"Input vector needs to be of size (L+1)², where L ∈ ℕ₀."))
     end
-    return MagneticFieldCoefficients(fill(SphericalHarmonicCoefficients(L), (3,1)))
+    return MagneticFieldCoefficients([SphericalHarmonicCoefficients(L) for j = 1:3, p=1:numPatches])
 end
 
 # constructor using t-design
@@ -94,14 +96,14 @@ function write(path::String, coeffs::MagneticFieldCoefficients)
 end
 
 # Size
-size(mfc::MagneticFieldCoefficients) = size(mfc.coeffs)
+size(mfc::MagneticFieldCoefficients, kargs...) = size(mfc.coeffs, kargs...)
 
 # Operations on MagneticFieldCoefficients
 function isapprox(mfc1::MagneticFieldCoefficients, mfc2::MagneticFieldCoefficients; kargs...)
     val = all(isapprox.(mfc1.coeffs,mfc2.coeffs;kargs...)) && isapprox(mfc1.radius,mfc2.radius;kargs...) && isapprox(mfc1.center,mfc2.center;kargs...)
-    if mfc1.ffp === nothing && mfc2.ffp === nothing
+    if isnothing(mfc1.ffp) && isnothing(mfc2.ffp)
         return  val
-    elseif mfc1.ffp !== nothing || mfc2.ffp !== nothing
+    elseif !isnothing(mfc1.ffp) || !isnothing(mfc2.ffp)
         @info "Only one of the coefficients has FFPs. Applying isapprox to the other values yields $val."
         return false
     else
@@ -266,4 +268,135 @@ function loadTDesignCoefficients(filename::String)
   coeffs_MF = MagneticFieldCoefficients(coeffs, radius, center)
 
   return coeffs_MF
+end
+
+# Shift coefficients into new expansion point
+function shift!(coeffsMF::MagneticFieldCoefficients, v::Union{Vector{T},Matrix{T}}) where T <: Real
+
+    # create matrix from vector
+    if typeof(v) <: Vector
+        v = hcat([v for i=1:size(coeffsMF,2)]...)
+    end
+
+    # test dimensions of shifting vector/array
+    if size(v,1) != 3
+        throw(DimensionMismatch("The shifting vector/matrix needs 3 entries but it has $(size(field,1))"))
+    elseif size(v,2) != size(coeffsMF,2)
+        throw(DimensionMismatch("The number of patches do not coincide: $(size(v,2)) != $(size(coeffsMF,2))"))
+    end
+
+    # shift coefficients by v
+    for p = 1:size(coeffsMF,2)
+        coeffsMF.coeffs[:,p] = SphericalHarmonicExpansions.translation.(coeffsMF.coeffs[:,p],[v[:,p]])
+    end
+
+    # change center and FFP position due to new expansion point
+    coeffsMF.center -= v
+    if !isnothing(coeffsMF.ffp)
+        coeffsMF.ffp -= v
+    end
+
+    return nothing
+end
+
+function shift(coeffsMF::MagneticFieldCoefficients,v::Union{Vector{T},Matrix{T}}) where T <: Real
+    coeffsMF_shifted = deepcopy(coeffsMF)
+    shift!(coeffsMF_shifted,v)
+    return coeffsMF_shifted
+end
+
+# shift coefficients into FFP
+"""
+    shiftFFP!(coeffsMF::MagneticFieldCoefficients)
+
+Shift magnetic-field coefficients into FFP (and calculate it if not available).
+"""
+function shiftFFP!(coeffsMF::MagneticFieldCoefficients)
+
+    # calculate FFP when not available
+    if isnothing(coeffsMF.ffp)
+        @info "Calculate FFP."
+        findFFP!(coeffsMF)
+    end
+
+    # shift coefficients into FFP
+    shift!(coeffsMF, coeffsMF.ffp)
+
+    return nothing
+end
+
+# Calculate the FFP of the magnetic field
+"""
+    findFFP(coeffsMF::MagneticFieldCoefficients; 
+            returnasmatrix::Bool=true)
+
+*Description:*  Newton method to find the FFPs of the magnetic fields\\
+ \\
+*Input:*
+- `coeffsMF`   - MagneticFieldCoefficients
+**kwargs:**
+- `returnasmatrix` - Boolean\\
+        true  -> return FFPs as Matrix with size (3,#Patches) (default)\\
+        false -> return FFPs as Array of NLsolve.SolverResults with size #Patches
+*Output:*
+- `ffp` - FFPs of the magnetic field
+"""
+function findFFP(coeffsMF::MagneticFieldCoefficients; 
+                 returnasmatrix::Bool=true)
+
+    # get spherical harmonic expansion
+    @polyvar x y z
+    expansion = sphericalHarmonicsExpansion.(coeffsMF.coeffs,x,y,z)
+
+    # return all FFPs in a matrix or as array with the solver results
+    ffp = returnasmatrix ? zeros(size(expansion)) : Array{NLsolve.SolverResults{Float64}}(undef,size(expansion,2))
+    for c in axes(expansion,2)
+
+        px = expansion[1,c]
+        py = expansion[2,c]
+        pz = expansion[3,c]
+        dpx = differentiate.(px,(x,y,z))
+        dpy = differentiate.(py,(x,y,z))
+        dpz = differentiate.(pz,(x,y,z))
+
+        function f!(fvec, xx)
+            fvec[1] = px((x,y,z)=>(xx[1],xx[2],xx[3]))
+            fvec[2] = py((x,y,z)=>(xx[1],xx[2],xx[3]))
+            fvec[3] = pz((x,y,z)=>(xx[1],xx[2],xx[3]))
+        end
+
+        function g!(fjac, xx)
+            fjac[1,1] = dpx[1]((x,y,z)=>(xx[1],xx[2],xx[3]))
+            fjac[1,2] = dpx[2]((x,y,z)=>(xx[1],xx[2],xx[3]))
+            fjac[1,3] = dpx[3]((x,y,z)=>(xx[1],xx[2],xx[3]))
+            fjac[2,1] = dpy[1]((x,y,z)=>(xx[1],xx[2],xx[3]))
+            fjac[2,2] = dpy[2]((x,y,z)=>(xx[1],xx[2],xx[3]))
+            fjac[2,3] = dpy[3]((x,y,z)=>(xx[1],xx[2],xx[3]))
+            fjac[3,1] = dpz[1]((x,y,z)=>(xx[1],xx[2],xx[3]))
+            fjac[3,2] = dpz[2]((x,y,z)=>(xx[1],xx[2],xx[3]))
+            fjac[3,3] = dpz[3]((x,y,z)=>(xx[1],xx[2],xx[3]))
+        end
+
+        if returnasmatrix
+            ffp[:,c] = nlsolve(f!,g!,[0.0;0.0;0.0],method=:newton,ftol=1e-16).zero
+        else
+            ffp[c] = nlsolve(f!,g!,[0.0;0.0;0.0],method=:newton,ftol=1e-16);
+        end
+    end
+
+    return ffp
+end
+
+"""
+    ffp = findFFP!(coeffsMF::MagneticFieldCoefficients)
+
+Find FFP and set it as coeffsMF.ffp.
+"""
+function findFFP!(coeffsMF::MagneticFieldCoefficients)
+    # find FFP
+    ffp = findFFP(coeffsMF, returnasmatrix=true)
+    # set FFP
+    coeffsMF.ffp = ffp
+
+    return ffp
 end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index b084d55..1a6bfa8 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -65,8 +65,8 @@ using Aqua
       coeffsMF = MagneticFieldCoefficients(coeffs)
       @test coeffsMF.coeffs == coeffs
       @test coeffsMF.radius == 0.0
-      @test coeffsMF.center == zeros(Float64,3)
-      @test coeffsMF.ffp === nothing
+      @test coeffsMF.center == zeros(Float64,3,1)
+      @test isnothing(coeffsMF.ffp)
 
       L = 1
       coeffsMF = MagneticFieldCoefficients(L)
@@ -76,13 +76,13 @@ using Aqua
       coeffsMF = MagneticFieldCoefficients(coeffs,0.042,zeros(3,1))
       @test coeffsMF.coeffs == coeffs
       @test coeffsMF.radius == 0.042
-      @test coeffsMF.center == zeros(Float64,3)
-      @test coeffsMF.ffp == zeros(3,1)
+      @test coeffsMF.center == zeros(Float64,3,1)
+      @test isnothing(coeffsMF.ffp)
 
       # constructor with t-design
       coeffsMF = MagneticFieldCoefficients(coeffs, tDes, zeros(Float64,3,1))
       @test coeffsMF.radius == 0.042
-      @test coeffsMF.center == zeros(Float64,3)
+      @test coeffsMF.center == zeros(Float64,3,1)
       @test coeffsMF.ffp == zeros(Float64,3,1) 
     
       # constructor with wrong FFP sizes

From ad2dc5f33f57f25d5c3eeb8b63ca8b439b00e0f3 Mon Sep 17 00:00:00 2001
From: Marija Boberg <m.boberg@uke.de>
Date: Thu, 13 Jul 2023 14:01:31 +0200
Subject: [PATCH 07/14] added getindex/setindex

---
 Project.toml                     |  1 +
 src/MagneticFieldCoefficients.jl | 53 ++++++++++++++++++++++++--------
 test/runtests.jl                 | 34 +++++++++++++++++---
 3 files changed, 71 insertions(+), 17 deletions(-)

diff --git a/Project.toml b/Project.toml
index f51020b..91679c2 100644
--- a/Project.toml
+++ b/Project.toml
@@ -16,6 +16,7 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 DocStringExtensions = "0.8, 0.9"
 MPIFiles = "0.12"
 MPIMagneticFields = "0.0.4, 0.0.5"
+NLsolve = "4.5"
 SphericalHarmonicExpansions = "0.1"
 Unitful = "1.11, 1.12"
 julia = "1"
diff --git a/src/MagneticFieldCoefficients.jl b/src/MagneticFieldCoefficients.jl
index 3ba39bd..7074370 100644
--- a/src/MagneticFieldCoefficients.jl
+++ b/src/MagneticFieldCoefficients.jl
@@ -1,4 +1,5 @@
-import Base.isapprox, Base.==,
+import Base.getindex, Base.setindex!,
+        Base.isapprox, Base.==,
         Base.+, Base.-, Base.*, Base./, 
         Base.write, Base.size
 
@@ -98,16 +99,36 @@ end
 # Size
 size(mfc::MagneticFieldCoefficients, kargs...) = size(mfc.coeffs, kargs...)
 
+# indexing
+getindex(mfc::MagneticFieldCoefficients, i) = MagneticFieldCoefficients(reshape(getindex(mfc.coeffs, :, i),3,length(i)), mfc.radius, reshape(getindex(mfc.center, :, i),3,length(i)), isnothing(mfc.ffp) ? nothing : reshape(getindex(mfc.ffp, :, i),3,length(i)))
+function setindex!(mfc1::MagneticFieldCoefficients, mfc2::MagneticFieldCoefficients, i) 
+    # setindex not possible in some cases 
+    if mfc1.radius != mfc2.radius
+        throw(DomainError([mfc1.radius,mfc2.radius],"Coefficients do not have the same measurement radius."))
+    elseif !isnothing(mfc1.ffp) && isnothing(mfc2.ffp)
+        throw(DomainError([mfc1.ffp,mfc2.ffp],"Coefficients do not provide an FFP."))
+    end
+    # set coefficients and center
+    setindex!(mfc1.coeffs, mfc2.coeffs, :, i) 
+    setindex!(mfc1.center, mfc2.center, :, i) 
+    # set FFP (ignore mfc2.ffp if mfc1.ffp === nothing)
+    if !isnothing(mfc1.ffp) && !isnothing(mfc2.ffp)
+        setindex!(mfc1.ffp, mfc2.ffp, :, i)
+    end
+    return nothing
+end
+
+
 # Operations on MagneticFieldCoefficients
 function isapprox(mfc1::MagneticFieldCoefficients, mfc2::MagneticFieldCoefficients; kargs...)
     val = all(isapprox.(mfc1.coeffs,mfc2.coeffs;kargs...)) && isapprox(mfc1.radius,mfc2.radius;kargs...) && isapprox(mfc1.center,mfc2.center;kargs...)
     if isnothing(mfc1.ffp) && isnothing(mfc2.ffp)
         return  val
-    elseif !isnothing(mfc1.ffp) || !isnothing(mfc2.ffp)
-        @info "Only one of the coefficients has FFPs. Applying isapprox to the other values yields $val."
-        return false
-    else
+    elseif !isnothing(mfc1.ffp) && !isnothing(mfc2.ffp)
         return val && isapprox(mfc1.ffp,mfc2.ffp;kargs...)
+    else
+        @info "Only one of the coefficients has FFPs. Applying isapprox to the other values yields $val."
+        return false   
     end
 end
 ==(mfc1::MagneticFieldCoefficients, mfc2::MagneticFieldCoefficients) = 
@@ -123,10 +144,10 @@ function +(mfc1::MagneticFieldCoefficients, mfc2::MagneticFieldCoefficients; for
         return MagneticFieldCoefficients(mfc1.coeffs .+ mfc2.coeffs, mfc1.radius, mfc1.center)
     end
     if mfc1.radius != mfc2.radius
-        throw(DomainError([mfc1.radius,mfc2.radius],"Coefficients do not have the same measurement radius."))
+        throw(DomainError([mfc1.radius,mfc2.radius],"Coefficients do not have the same measurement radius. (Use `force = true`.)"))
     end
     if mfc1.center != mfc2.center
-        throw(DomainError([mfc1.center,mfc2.center],"Coefficients do not have the same measurement center."))
+        throw(DomainError([mfc1.center,mfc2.center],"Coefficients do not have the same measurement center. (Use `force = true`.)"))
     end
     return MagneticFieldCoefficients(mfc1.coeffs .+ mfc2.coeffs, mfc1.radius, mfc1.center)
 end
@@ -271,12 +292,7 @@ function loadTDesignCoefficients(filename::String)
 end
 
 # Shift coefficients into new expansion point
-function shift!(coeffsMF::MagneticFieldCoefficients, v::Union{Vector{T},Matrix{T}}) where T <: Real
-
-    # create matrix from vector
-    if typeof(v) <: Vector
-        v = hcat([v for i=1:size(coeffsMF,2)]...)
-    end
+function shift!(coeffsMF::MagneticFieldCoefficients, v::Matrix{T}) where T <: Real
 
     # test dimensions of shifting vector/array
     if size(v,1) != 3
@@ -299,6 +315,17 @@ function shift!(coeffsMF::MagneticFieldCoefficients, v::Union{Vector{T},Matrix{T
     return nothing
 end
 
+function shift!(coeffsMF::MagneticFieldCoefficients, v::Vector{T}) where T <: Real
+
+    # create matrix from vector
+    v = hcat([v for i=1:size(coeffsMF,2)]...)
+
+    # shift coefficients
+    shift!(coeffsMF,v)
+
+    return nothing
+end
+
 function shift(coeffsMF::MagneticFieldCoefficients,v::Union{Vector{T},Matrix{T}}) where T <: Real
     coeffsMF_shifted = deepcopy(coeffsMF)
     shift!(coeffsMF_shifted,v)
diff --git a/test/runtests.jl b/test/runtests.jl
index 1a6bfa8..0f5ff7c 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -70,8 +70,12 @@ using Aqua
 
       L = 1
       coeffsMF = MagneticFieldCoefficients(L)
-      @test size(coeffsMF.coeffs) == (3,1)
+      @test size(coeffsMF) == (3,1)
       [@test coeffsMF.coeffs[j,1].L == L for j=1:3]
+      
+      # multiple patches
+      coeffsMF = MagneticFieldCoefficients(L, numPatches = 2)
+      @test size(coeffsMF) == (3,2)
 
       coeffsMF = MagneticFieldCoefficients(coeffs,0.042,zeros(3,1))
       @test coeffsMF.coeffs == coeffs
@@ -88,9 +92,6 @@ using Aqua
       # constructor with wrong FFP sizes
       @test_throws DimensionMismatch MagneticFieldCoefficients(coeffs, tDes, zeros(2,1))
       @test_throws DimensionMismatch MagneticFieldCoefficients(coeffs, tDes, zeros(3,2))
-
-      # test size
-      @test size(coeffsMF) == (3,1)
     end
 
     @testset "MagneticFieldCoefficient operations" begin
@@ -132,6 +133,31 @@ using Aqua
       @test isapprox(c1 * 2, c2)
       @test isapprox(c2 / 2, c1)
 
+      # indexing
+        # setup MagneticFieldCoefficients with multiple patches
+        csh = reshape(hcat(SphericalHarmonicCoefficients.([rand(16) for i=1:12])...), 3,4)
+        center = rand(3,4)
+        ffp = rand(3,4)
+        coeffsMF = MagneticFieldCoefficients(csh,0.01,center,ffp)
+        c1 = MagneticFieldCoefficients(csh[:,2:2],0.01,center[:,2:2],ffp[:,2:2])
+        c2 = MagneticFieldCoefficients(csh[:,3:4],0.01,center[:,3:4],ffp[:,3:4])
+
+        # test getindex
+        @test isapprox(coeffsMF[2], c1)
+        @test isapprox(coeffsMF[3:4], c2)
+
+        # test setindex
+        coeffsMF[1] = coeffsMF[2]
+        @test isapprox(coeffsMF[1], c1)
+        coeffsMF[1:2] = coeffsMF[3:4]
+        @test isapprox(coeffsMF[1:2], c2)
+
+        # test errors 
+        c3 = MagneticFieldCoefficients(csh[:,2:2],0.02,center[:,2:2],ffp[:,2:2]) # different radius
+        c4 = MagneticFieldCoefficients(csh[:,2:2],0.01,center[:,2:2]) # no FFP
+        
+        @test_throws DomainError coeffsMF[1] = c3 # different radius
+        @test_throws DomainError coeffsMF[1] = c4 # no FFP
     end
 
     @testset "Load/write data from/to file" begin

From da39c687eda0f63583c4f28fe3c42350fd2c0ca9 Mon Sep 17 00:00:00 2001
From: Marija Boberg <m.boberg@uke.de>
Date: Thu, 13 Jul 2023 17:59:35 +0200
Subject: [PATCH 08/14] improved file handling

---
 src/MagneticFieldCoefficients.jl | 95 ++++++++++++++++++--------------
 test/runtests.jl                 | 11 +++-
 2 files changed, 64 insertions(+), 42 deletions(-)

diff --git a/src/MagneticFieldCoefficients.jl b/src/MagneticFieldCoefficients.jl
index 7074370..5da11c7 100644
--- a/src/MagneticFieldCoefficients.jl
+++ b/src/MagneticFieldCoefficients.jl
@@ -47,35 +47,73 @@ end
 MagneticFieldCoefficients(coeffs::Array{SphericalHarmonicCoefficients,2}, tDesign::SphericalTDesign, ffp=nothing) = MagneticFieldCoefficients(coeffs,ustrip(Unitful.m(tDesign.radius)), ustrip.(Unitful.m.(tDesign.center)), ffp)
 
 
-# read coefficients from an HDF5-file
+# read coefficients/measurement data from an HDF5-file
 function MagneticFieldCoefficients(path::String)
 
+    coeffsMF = h5open(path,"r") do file
+        if haskey(HDF5.root(file),"/coeffs")
+            # coefficients available    
+            return loadCoefficients(path)
+
+        elseif haskey(HDF5.root(file),"/fields") && haskey(HDF5.root(file),"/positions/tDesign")
+            # calculate coefficients
+            return loadTDesignCoefficients(path)
+        
+        else
+            throw(ErrorException("Unknown file structure."))
+        end
+    end
+
+    return coeffsMF
+end
+
+# read coefficients from an HDF5-file
+function loadCoefficients(path::String)
+
     # load spherical harmonic coefficients
     shcoeffs = SphericalHarmonicCoefficients(path)
 
     coeffsMF = h5open(path,"r") do file
-        if haskey(HDF5.root(file), "/radius") 
-            # file contains all relevant information
-            radius = read(file, "/radius")
-            center = read(file, "/center")
-            if haskey(HDF5.root(file), "/ffp")
-                ffp = read(file, "/ffp")
-                return MagneticFieldCoefficients(shcoeffs, radius, center, ffp)
-            else
-                # field has not FFP -> ffp = nothing
-                return MagneticFieldCoefficients(shcoeffs, radius, center)
-            end
-        else
-            # convert file of SphericalHarmonicCoefficients into MagneticFieldCoefficients
-            # -> set all additional informations to 0 or nothing
-            # use radius = 0.01 as default value
-            return MagneticFieldCoefficients(shcoeffs, 0.01)
+        if !haskey(HDF5.root(file), "/radius") || !haskey(HDF5.root(file), "/center")
+            @warn "The file does not provide all necessary measurement informations."
         end
+
+        # set all informations not given in the file to 0 or nothing
+        radius = haskey(HDF5.root(file), "/radius") ? read(file, "/radius") : 0.0
+        center = haskey(HDF5.root(file), "/center") ? read(file, "/center") : zeros(size(shcoeffs))
+        ffp = haskey(HDF5.root(file), "/ffp") ? read(file, "/ffp") : nothing
+
+        return MagneticFieldCoefficients(shcoeffs, radius, center, ffp)
     end
 
     return coeffsMF
 end
 
+# load coefficients from measurement data
+# TODO: This should be merged with and moved to MPIFiles
+function loadTDesignCoefficients(filename::String)
+    field, radius, N, t, center, correction  = h5open(filename, "r") do file
+      field = read(file,"/fields") 		# measured field (size: 3 x #points x #patches)
+      radius = read(file,"/positions/tDesign/radius")	# radius of the measured ball
+      N = read(file,"/positions/tDesign/N")		# number of points of the t-design
+      t = read(file,"/positions/tDesign/t")		# t of the t-design
+      center = read(file,"/positions/tDesign/center")	# center of the measured ball
+      correction = read(file, "/sensor/correctionTranslation")
+      return field, radius, N, t, center, correction
+    end
+    @polyvar x y z
+    tDes = MPIFiles.loadTDesign(Int(t),N,radius*u"m", center.*u"m")
+    coeffs, = magneticField(tDes, field, x,y,z)
+    for c=1:size(coeffs,2), j = 1:3
+      coeffs[j, c] = SphericalHarmonicExpansions.translation(coeffs[j, c], correction[:, j])
+    end
+  
+    coeffs_MF = MagneticFieldCoefficients(coeffs, radius, center)
+  
+    return coeffs_MF
+end
+
+
 # write coefficients to an HDF5-file
 function write(path::String, coeffs::MagneticFieldCoefficients)
 
@@ -268,29 +306,6 @@ function magneticField(coords::AbstractArray{T, 2}, field::Union{AbstractArray{T
   return coeffs, expansion, func
 end
 
-#TODO: This should be merged with and moved to MPIFiles
-function loadTDesignCoefficients(filename::String)
-  field, radius, N, t, center, correction  = h5open(filename, "r") do file
-    field = read(file,"/fields") 		# measured field (size: 3 x #points x #patches)
-    radius = read(file,"/positions/tDesign/radius")	# radius of the measured ball
-    N = read(file,"/positions/tDesign/N")		# number of points of the t-design
-    t = read(file,"/positions/tDesign/t")		# t of the t-design
-    center = read(file,"/positions/tDesign/center")	# center of the measured ball
-    correction = read(file, "/sensor/correctionTranslation")
-    return field, radius, N, t, center, correction
-  end
-  @polyvar x y z
-  tDes = MPIFiles.loadTDesign(Int(t),N,radius*u"m", center.*u"m")
-  coeffs, = magneticField(tDes, field, x,y,z)
-  for c=1:size(coeffs,2), j = 1:3
-    coeffs[j, c] = SphericalHarmonicExpansions.translation(coeffs[j, c], correction[:, j])
-  end
-
-  coeffs_MF = MagneticFieldCoefficients(coeffs, radius, center)
-
-  return coeffs_MF
-end
-
 # Shift coefficients into new expansion point
 function shift!(coeffsMF::MagneticFieldCoefficients, v::Matrix{T}) where T <: Real
 
diff --git a/test/runtests.jl b/test/runtests.jl
index 0f5ff7c..8b51de4 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -169,7 +169,7 @@ using Aqua
       @test isapprox(field[0.01,0.01,0.01], [-0.01,-0.01,0.02], atol=ε)
 
       # get coefficients
-      coeffsMF = MPISphericalHarmonics.loadTDesignCoefficients(filename)
+      coeffsMF = MagneticFieldCoefficients(filename)
       @test isapprox(coeffsMF.radius, 0.042, atol=ε) # radius
       @test isapprox(coeffsMF.coeffs[1][1,1], -1.0, atol=1e-10) # gradient (x)
       @test isapprox(coeffsMF.coeffs[2][1,-1], -1.0, atol=1e-10) # gradient (y)
@@ -181,6 +181,7 @@ using Aqua
       filename3 = "Coeffs2.h5"
       filename4 = "Coeffs3.h5"
       filenameW = "CoeffsW.h5"
+      filenameE = "CoeffsE.h5"
       write(filename2, coeffsMF.coeffs)
       # add radius and center
       cp(filename2, filename3)
@@ -196,7 +197,7 @@ using Aqua
 
       # only coefficients (no further informations given)
       coeffsTest = MagneticFieldCoefficients(filename2)
-      @test isapprox(coeffsTest.radius, 0.01, atol=ε) # radius
+      @test isapprox(coeffsTest.radius, 0.0, atol=ε) # radius
  
       # with given radius & center
       coeffsTest = MagneticFieldCoefficients(filename3)
@@ -217,6 +218,11 @@ using Aqua
       @test isapprox(coeffsW.center, zeros(3), atol=ε) # center
       @test coeffsW.ffp == zeros(3,1) # FFP
       coeffsW = nothing # This maybe masks an implementation error
+
+      # test error
+      h5write(filenameE, "test", 0)
+      @test_throws ErrorException MagneticFieldCoefficients(filenameE)
+
       GC.gc()
 
       # remove test files
@@ -224,6 +230,7 @@ using Aqua
       rm(filename3)
       rm(filename4)
       rm(filenameW)
+      rm(filenameE)
     end
 
     @testset "SphericalHarmonicsDefinedField (multiple patches)" begin

From ba30fa5a03074c255a8500a52448d27a5059ee98 Mon Sep 17 00:00:00 2001
From: Marija Boberg <m.boberg@uke.de>
Date: Wed, 19 Jul 2023 18:54:13 +0200
Subject: [PATCH 09/14] return only coefficients with magneticField()

---
 src/MPISphericalHarmonics.jl     | 15 ++++-----------
 src/MagneticFieldCoefficients.jl | 22 ++++------------------
 test/runtests.jl                 |  9 ++++-----
 3 files changed, 12 insertions(+), 34 deletions(-)

diff --git a/src/MPISphericalHarmonics.jl b/src/MPISphericalHarmonics.jl
index 732367e..ceb983e 100644
--- a/src/MPISphericalHarmonics.jl
+++ b/src/MPISphericalHarmonics.jl
@@ -27,17 +27,10 @@ end
 
 function SphericalHarmonicsDefinedField(filename::String)
 
-  func = h5open(filename,"r") do file
-    if haskey(file,"coeffs") 
-      # load coefficients
-      coeffs_MF = MagneticFieldCoefficients(filename)
-    else
-      # load measured field
-      coeffs_MF = loadTDesignCoefficients(filename)
-    end
-
-    return fastfunc.(coeffs_MF.coeffs)
-  end
+  # load coefficients
+  mfc = MagneticFieldCoefficients(filename)
+  # get field function
+  func = fastfunc.(mfc.coeffs)
 
   return SphericalHarmonicsDefinedField(func=func)
 end
diff --git a/src/MagneticFieldCoefficients.jl b/src/MagneticFieldCoefficients.jl
index 5da11c7..c138f79 100644
--- a/src/MagneticFieldCoefficients.jl
+++ b/src/MagneticFieldCoefficients.jl
@@ -101,9 +101,8 @@ function loadTDesignCoefficients(filename::String)
       correction = read(file, "/sensor/correctionTranslation")
       return field, radius, N, t, center, correction
     end
-    @polyvar x y z
     tDes = MPIFiles.loadTDesign(Int(t),N,radius*u"m", center.*u"m")
-    coeffs, = magneticField(tDes, field, x,y,z)
+    coeffs = magneticField(tDes, field)
     for c=1:size(coeffs,2), j = 1:3
       coeffs[j, c] = SphericalHarmonicExpansions.translation(coeffs[j, c], correction[:, j])
     end
@@ -219,16 +218,14 @@ end
 
 ## Load coefficients from t-design measurement ##
 """
-    magneticField(tDesign::SphericalTDesign, field::Union{AbstractArray{T,2},AbstractArray{T,3}}, 
-		       x::Variable, y::Variable, z::Variable;
+    magneticField(tDesign::SphericalTDesign, field::Union{AbstractArray{T,2},AbstractArray{T,3}};
 		       L::Int=Int(tDesign.T/2),
 		       calcSolid::Bool=true) where T <: Real
-*Description:*  Calculation of the spherical harmonic coefficients and expansion based on the measured t-design\\
+*Description:*  Calculation of the spherical harmonic coefficients based on the measured t-design\\
  \\
 *Input:*
 - `tDesign`	- Measured t-design (type: SphericalTDesign)
 - `field`       - Measured field (size = (J,N,C)) with J <= 3
-- `x, y, z`     - Cartesian coordinates
 **kwargs:**
 - `L`           - Order up to which the coeffs be calculated (default: t/2)
 - `calcSolid`   - Boolean (default: true)\\
@@ -236,11 +233,8 @@ end
     true -> solid coefficients
 *Output:*
 - `coeffs`    - spherical/solid coefficients, type: Array{SphericalHarmonicCoefficients}(3,C)
-- `expansion` - related expansion (Cartesian polynomial), type: Array{AbstractPolynomialLike}(3,C)
-- `func`      - expansion converted to a function, type: Array{Function}(3,C)
 """
 function magneticField(tDesign::SphericalTDesign, field::Union{AbstractArray{T,2},AbstractArray{T,3}}, 
-		       x::Variable, y::Variable, z::Variable;
 		       L::Int=floor(Int,tDesign.T/2),
 		       calcSolid::Bool=true) where T <: Real
 
@@ -256,13 +250,11 @@ function magneticField(tDesign::SphericalTDesign, field::Union{AbstractArray{T,2
   
   return magneticField(coords, field, 
 		       R, center, L,
-		       x, y, z,
 		       calcSolid)
 end
 
 function magneticField(coords::AbstractArray{T, 2}, field::Union{AbstractArray{T, 2}, AbstractArray{T, 3}}, 
                        R::T, center::Vector{T}, L::Int,
-                       x::Variable, y::Variable, z::Variable, 
                        calcSolid::Bool=true) where T <: Real
 
   # transpose coords if its dimensions do not fit
@@ -277,8 +269,6 @@ function magneticField(coords::AbstractArray{T, 2}, field::Union{AbstractArray{T
     throw(DimensionMismatch("The field vector does not match the size of the tdesign: $(size(field,2)) != $(size(coords,2))"))
   end
 
-  func= Array{Function}(undef,size(field,1),size(field,3))
-  expansion = Array{Polynomial}(undef,size(field,1),size(field,3))
   coeffs = Array{SphericalHarmonicCoefficients}(undef,size(field,1),size(field,3))
 
   # rescale coordinates to t-design on unit sphere
@@ -296,14 +286,10 @@ function magneticField(coords::AbstractArray{T, 2}, field::Union{AbstractArray{T
       if calcSolid
           solid!(coeffs[j,c])
       end
-
-      # calculation of the expansion
-      expansion[j,c] = sphericalHarmonicsExpansion(coeffs[j,c],x,y,z) + 0*x;
-      func[j,c] = @fastfunc expansion[j,c]+0*x+0*y+0*z
     end
   end
 
-  return coeffs, expansion, func
+  return coeffs
 end
 
 # Shift coefficients into new expansion point
diff --git a/test/runtests.jl b/test/runtests.jl
index 8b51de4..f5e23fd 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -22,8 +22,7 @@ using Aqua
     fieldValues = hcat([idealField[ustrip.(Unitful.m.(pos))...] for pos in tDes]...)
 
     ## Calculate coefficients
-    MPISphericalHarmonics.@polyvar x y z
-    coeffs, = MPISphericalHarmonics.magneticField(tDes, fieldValues, x,y,z)
+    coeffs = MPISphericalHarmonics.magneticField(tDes, fieldValues)
 
     ## Test coefficients
     numCoeffs = length(coeffs[1,1].c)
@@ -43,14 +42,14 @@ using Aqua
       L = floor(Int,tDes.T/2)
   
       # transposed positions
-      coeffsTest, = MPISphericalHarmonics.magneticField(coords', fieldValues, R, center, L, x, y, z)
+      coeffsTest = MPISphericalHarmonics.magneticField(coords', fieldValues, R, center, L)
       for j=1:3
         @test isapprox(coeffs[j,1].c, coeffsTest[j,1].c, atol = 1e-10)
       end
 
       # Errors
-      @test_throws DimensionMismatch MPISphericalHarmonics.magneticField(coords, zeros(4,length(tDes)), R, center, L, x, y, z) # >3 field values in the first dimension 
-      @test_throws DimensionMismatch MPISphericalHarmonics.magneticField(coords, fieldValues[:,1:end-1], R, center, L, x, y, z) # number of field values != number of measured positions
+      @test_throws DimensionMismatch MPISphericalHarmonics.magneticField(coords, zeros(4,length(tDes)), R, center, L) # >3 field values in the first dimension 
+      @test_throws DimensionMismatch MPISphericalHarmonics.magneticField(coords, fieldValues[:,1:end-1], R, center, L) # number of field values != number of measured positions
     end
 
 

From a993ed7ca5c31f5cea269259032a2b2859066f59 Mon Sep 17 00:00:00 2001
From: "github-actions[bot]"
 <41898282+github-actions[bot]@users.noreply.github.com>
Date: Wed, 26 Jul 2023 09:09:33 +0200
Subject: [PATCH 10/14] CompatHelper: bump compat for MPIMagneticFields to
 0.0.6, (keep existing compat) (#12)

Co-authored-by: CompatHelper Julia <compathelper_noreply@julialang.org>
Co-authored-by: Jonas Schumacher <github735@jonasschumacher.de>
---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 91679c2..d26acdf 100644
--- a/Project.toml
+++ b/Project.toml
@@ -15,7 +15,7 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 [compat]
 DocStringExtensions = "0.8, 0.9"
 MPIFiles = "0.12"
-MPIMagneticFields = "0.0.4, 0.0.5"
+MPIMagneticFields = "0.0.4, 0.0.5, 0.0.6"
 NLsolve = "4.5"
 SphericalHarmonicExpansions = "0.1"
 Unitful = "1.11, 1.12"

From fae627cc9fb3264a33c6851e2f7dfb1d0dbc9068 Mon Sep 17 00:00:00 2001
From: "github-actions[bot]"
 <41898282+github-actions[bot]@users.noreply.github.com>
Date: Wed, 26 Jul 2023 09:10:52 +0200
Subject: [PATCH 11/14] CompatHelper: add new compat entry for NLsolve at
 version 4, (keep existing compat) (#13)

Co-authored-by: CompatHelper Julia <compathelper_noreply@julialang.org>
Co-authored-by: Jonas Schumacher <github735@jonasschumacher.de>
---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index d26acdf..48499ac 100644
--- a/Project.toml
+++ b/Project.toml
@@ -16,7 +16,7 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 DocStringExtensions = "0.8, 0.9"
 MPIFiles = "0.12"
 MPIMagneticFields = "0.0.4, 0.0.5, 0.0.6"
-NLsolve = "4.5"
+NLsolve = "4"
 SphericalHarmonicExpansions = "0.1"
 Unitful = "1.11, 1.12"
 julia = "1"

From 2a53edc71cf29899cdae13fbbe27a78db50cd401 Mon Sep 17 00:00:00 2001
From: "github-actions[bot]"
 <41898282+github-actions[bot]@users.noreply.github.com>
Date: Tue, 26 Sep 2023 10:46:19 +0200
Subject: [PATCH 12/14] CompatHelper: bump compat for MPIFiles to 0.13, (keep
 existing compat) (#14)

Co-authored-by: CompatHelper Julia <compathelper_noreply@julialang.org>
---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 48499ac..80e20ef 100644
--- a/Project.toml
+++ b/Project.toml
@@ -14,7 +14,7 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [compat]
 DocStringExtensions = "0.8, 0.9"
-MPIFiles = "0.12"
+MPIFiles = "0.12, 0.13"
 MPIMagneticFields = "0.0.4, 0.0.5, 0.0.6"
 NLsolve = "4"
 SphericalHarmonicExpansions = "0.1"

From a2f46061f9c0b618f10506499213895f4557f159 Mon Sep 17 00:00:00 2001
From: Jonas Schumacher <github735@jonasschumacher.de>
Date: Tue, 26 Sep 2023 12:53:04 +0200
Subject: [PATCH 13/14] Remove nightly and older Julia versions from CI

---
 .github/workflows/CI.yml | 6 ++----
 1 file changed, 2 insertions(+), 4 deletions(-)

diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 510c110..fffa84f 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -18,9 +18,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.7'
-          - '1.8'
-          - 'nightly'
+          - '1'
         os:
           - ubuntu-latest
           - macOS-latest
@@ -62,4 +60,4 @@ jobs:
       - run: julia --project=docs docs/make.jl
         env:
           GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
-          DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }}
\ No newline at end of file
+          DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }}

From 97dfbb27da90958f38c9363a6906174dbb93b53b Mon Sep 17 00:00:00 2001
From: Marija Boberg <m.boberg@uke.de>
Date: Wed, 27 Sep 2023 11:38:25 +0200
Subject: [PATCH 14/14] added functions to get field-specific values

---
 src/MPISphericalHarmonics.jl     |  3 +-
 src/MagneticFieldCoefficients.jl | 52 ++++++++++++++++++++++++++++++--
 test/runtests.jl                 |  5 +++
 3 files changed, 57 insertions(+), 3 deletions(-)

diff --git a/src/MPISphericalHarmonics.jl b/src/MPISphericalHarmonics.jl
index ceb983e..8d4b0a1 100644
--- a/src/MPISphericalHarmonics.jl
+++ b/src/MPISphericalHarmonics.jl
@@ -14,6 +14,7 @@ import Base.length
 # load MagneticFieldCoefficients
 include("MagneticFieldCoefficients.jl")
 export MagneticFieldCoefficients
+export getOffset, getGradient, getJacobian
 export shift, shift!, shiftFFP!
 
 ## SphericalHarmonicsDefinedField ##
@@ -21,7 +22,7 @@ export SphericalHarmonicsDefinedField
 export selectPatch, length
 
 Base.@kwdef mutable struct SphericalHarmonicsDefinedField <: AbstractMagneticField
-  func::Array{SphericalHarmonicExpansions.StaticPolynomials.Polynomial, 2}
+  func::Array{Union{Function,SphericalHarmonicExpansions.StaticPolynomials.Polynomial}, 2}
   patch::Integer = 1
 end
 
diff --git a/src/MagneticFieldCoefficients.jl b/src/MagneticFieldCoefficients.jl
index c138f79..682005e 100644
--- a/src/MagneticFieldCoefficients.jl
+++ b/src/MagneticFieldCoefficients.jl
@@ -36,11 +36,13 @@ MagneticFieldCoefficients(coeffs::Array{SphericalHarmonicCoefficients,2}, radius
 MagneticFieldCoefficients(coeffs::Array{SphericalHarmonicCoefficients,2}, radius::Float64, center::VecOrMat{Float64}) = MagneticFieldCoefficients(coeffs,radius,center,nothing)
 MagneticFieldCoefficients(coeffs::Array{SphericalHarmonicCoefficients,2}, radius::Float64, center::Vector{Float64},ffp::Union{Matrix{Float64},Nothing}) = MagneticFieldCoefficients(coeffs,radius,hcat([center for p=1:size(coeffs,2)]...),ffp)
 
-function MagneticFieldCoefficients(L::Int; numPatches::Int=1)
+function MagneticFieldCoefficients(L::Int, R::Float64=1.0, solid::Bool=true;
+                                   radius::Float64=0.0, center::VecOrMat{Float64}=[0.0,0.0,0.0],
+                                   numPatches::Int=1)
     if L<0
         throw(DomainError(L,"Input vector needs to be of size (L+1)², where L ∈ ℕ₀."))
     end
-    return MagneticFieldCoefficients([SphericalHarmonicCoefficients(L) for j = 1:3, p=1:numPatches])
+    return MagneticFieldCoefficients([SphericalHarmonicCoefficients(L,R,solid) for j = 1:3, p=1:numPatches],radius,center)
 end
 
 # constructor using t-design
@@ -215,6 +217,51 @@ end
 -(value::Number, mfc::MagneticFieldCoefficients) = MagneticFieldCoefficients(value .- mfc.coeffs, mfc.radius, mfc.center, mfc.ffp)
 *(value::Number, mfc::MagneticFieldCoefficients) = *(mfc::MagneticFieldCoefficients, value)
 
+## get some field-specific values
+"""
+    getOffset(mfc::MagneticFieldCoefficients, idx::AbstractUnitRange{Int64}=axes(mfc.coeffs,2))
+
+Get the offset of the field described by mfc[idx].
+"""
+getOffset(mfc::MagneticFieldCoefficients, idx::AbstractUnitRange{Int64}=axes(mfc.coeffs,2)) = [c[0,0] for c in mfc.coeffs[idx]]
+"""
+    getGradient(mfc::MagneticFieldCoefficients, idx::AbstractUnitRange{Int64}=axes(mfc.coeffs,2))
+
+Get the gradient of the field described by mfc[idx].
+"""
+function getGradient(mfc::MagneticFieldCoefficients, idx::AbstractUnitRange{Int64}=axes(mfc.coeffs,2))
+    # ideal gradient would be [[mfc.coeffs[1,p][1,1], mfc.coeffs[2,p][1,-1], mfc.coeffs[3,p][1,0]] for p in idx]
+
+    G = Vector{Float64}[]
+    for p in idx
+        # calculate jacobian matrix
+        @polyvar x y z
+        expansion = sphericalHarmonicsExpansion.(mfc.coeffs[:,p],[x],[y],[z])
+        jexp = differentiate(expansion,[x,y,z]);
+        push!(G,[jexp[i,i]((x,y,z) => [0.0,0.0,0.0]) for i=1:3]) # gradient = diag(jacobian matrix)
+    end
+
+    return G
+end
+"""
+    getJacobian(mfc::MagneticFieldCoefficients, idx::AbstractUnitRange{Int64}=axes(mfc.coeffs,2))
+
+Get the Jacobian matrix of the field described by mfc[idx].
+"""
+function getJacobian(mfc::MagneticFieldCoefficients, idx::AbstractUnitRange{Int64}=axes(mfc.coeffs,2))
+    # ideal Jacobian matrix would be [vcat([[c[1,1], c[1,-1], c[1,0]]' for c in mfc.coeffs[:,p]]...) for p in idx]
+
+    J = Matrix{Float64}[]
+    for p in idx
+        # calculate jacobian matrix
+        @polyvar x y z
+        expansion = sphericalHarmonicsExpansion.(mfc.coeffs[:,p],[x],[y],[z])
+        jexp = differentiate(expansion,[x,y,z]);
+        push!(J,[jexp[i,j]((x,y,z) => [0.0,0.0,0.0]) for i=1:3, j=1:3]) # jacobian matrix
+    end
+
+    return J
+end
 
 ## Load coefficients from t-design measurement ##
 """
@@ -379,6 +426,7 @@ function findFFP(coeffsMF::MagneticFieldCoefficients;
     # return all FFPs in a matrix or as array with the solver results
     ffp = returnasmatrix ? zeros(size(expansion)) : Array{NLsolve.SolverResults{Float64}}(undef,size(expansion,2))
     for c in axes(expansion,2)
+        print("$c ")
 
         px = expansion[1,c]
         py = expansion[2,c]
diff --git a/test/runtests.jl b/test/runtests.jl
index f5e23fd..92aef06 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -91,6 +91,11 @@ using Aqua
       # constructor with wrong FFP sizes
       @test_throws DimensionMismatch MagneticFieldCoefficients(coeffs, tDes, zeros(2,1))
       @test_throws DimensionMismatch MagneticFieldCoefficients(coeffs, tDes, zeros(3,2))
+
+      # test offset/gradient/Jacobian
+      @test isapprox(getOffset(coeffsMF),[0.0], atol=1e-10)
+      @test isapprox(getGradient(coeffsMF), [[-1.0,-1.0,2.0]], atol=1e-10)
+      @test isapprox(getJacobian(coeffsMF), [[-1.0 0.0 0.0; 0.0 -1.0 0.0; 0.0 0.0 2.0]], atol=1e-10)
     end
 
     @testset "MagneticFieldCoefficient operations" begin