From e7511dd16d4cc0ee365ef9659b28e16ecbeb7934 Mon Sep 17 00:00:00 2001 From: Marija Boberg Date: Wed, 7 Aug 2024 18:33:17 +0200 Subject: [PATCH] enable loading of older measurement data and findFFP(SHC) --- src/MPISphericalHarmonics.jl | 3 + src/MagneticFieldCoefficients.jl | 40 ++++++------ src/utils.jl | 108 +++++++++++++++++++++++++++++++ 3 files changed, 130 insertions(+), 21 deletions(-) create mode 100644 src/utils.jl diff --git a/src/MPISphericalHarmonics.jl b/src/MPISphericalHarmonics.jl index e0ae2ea..906fba2 100644 --- a/src/MPISphericalHarmonics.jl +++ b/src/MPISphericalHarmonics.jl @@ -72,4 +72,7 @@ function checkPatch(field::SphericalHarmonicsDefinedField, patch::Int) return nothing end + +# utilty functions +include("utils.jl") end diff --git a/src/MagneticFieldCoefficients.jl b/src/MagneticFieldCoefficients.jl index b8f9b48..82a23e4 100644 --- a/src/MagneticFieldCoefficients.jl +++ b/src/MagneticFieldCoefficients.jl @@ -92,10 +92,9 @@ function MagneticFieldCoefficients(path::String) # coefficients available return loadCoefficients(path) - elseif haskey(HDF5.root(file), "/fields") && - haskey(HDF5.root(file), "/positions/tDesign") + elseif haskey(HDF5.root(file), "/fields") # calculate coefficients - return loadTDesignCoefficients(path) + return calcTDesignCoefficients(path) else throw(ErrorException("Unknown file structure.")) @@ -127,20 +126,17 @@ function loadCoefficients(path::String) return coeffsMF end -# load coefficients from measurement data +# calculate 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 +function calcTDesignCoefficients(filename::String) + # load the measurement data + field, radius, N, t, center, correction = loadMagneticFieldMeasurementData(filename) + + # calculate the coefficients tDes = MPIFiles.loadTDesign(Int(t), N, radius * u"m", center .* u"m") coeffs = magneticField(tDes, field) + + # apply the correction of the sensors for c = 1:size(coeffs, 2), j = 1:3 coeffs[j, c] = SphericalHarmonicExpansions.translation(coeffs[j, c], correction[:, j]) end @@ -545,7 +541,7 @@ getVolumeIndex(vol::Union{String, Volume}) = getVolumeIndex(Symbol(vol)) *Description:* Newton method to find the FFPs of the magnetic fields\\ \\ *Input:* -- `coeffsMF` - MagneticFieldCoefficients +- `coeffsMF` - MagneticFieldCoefficients or Matrix{SphericalHarmonicCoefficients} **kwargs:** - `vol` - provide the volume where the FFP search should be done\\ options: $(instances(Volume))\\ @@ -558,14 +554,14 @@ getVolumeIndex(vol::Union{String, Volume}) = getVolumeIndex(Symbol(vol)) *Output:* - `ffp` - FFPs of the magnetic field """ -function findFFP(coeffsMF::MagneticFieldCoefficients; +function findFFP(coeffs::Matrix{SphericalHarmonicCoefficients}; vol::Union{Symbol, String, Volume} = xyz, start::Union{Vector{Vector{Float64}},Vector{Float64}} = [[0.0;0.0;0.0]], # start value returnasmatrix::Bool = true) # get spherical harmonic expansion @polyvar x y z - expansion = sphericalHarmonicsExpansion.(coeffsMF.coeffs, x, y, z) + expansion = sphericalHarmonicsExpansion.(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)) @@ -577,10 +573,10 @@ function findFFP(coeffsMF::MagneticFieldCoefficients; if length.(start) != 3 .* ones(length(start)) throw(DimensionMismatch("Length of start vector needs to be 3, not $(length.(start)).")) end - if length(start) == 1 && size(coeffsMF,2) !== 1 - start = repeat(start,size(coeffsMF,2)) - elseif length(start) != size(coeffsMF,2) - throw(DimensionMismatch("The numbers of start vectors ($(length(start))) has to be equal to the number of patches ($(size(coeffsMF,2))).")) + if length(start) == 1 && size(coeffs,2) !== 1 + start = repeat(start,size(coeffs,2)) + elseif length(start) != size(coeffs,2) + throw(DimensionMismatch("The numbers of start vectors ($(length(start))) has to be equal to the number of patches ($(size(coeffs,2))).")) end # get index set of the volume @@ -618,6 +614,8 @@ function findFFP(coeffsMF::MagneticFieldCoefficients; return ffp end +# MagneticFieldCoefficients as input +findFFP(coeffsMF::MagneticFieldCoefficients; kargs...) = findFFP(coeffsMF.coeffs; kargs...) """ ffp = findFFP!(coeffsMF::MagneticFieldCoefficients) diff --git a/src/utils.jl b/src/utils.jl new file mode 100644 index 0000000..048455f --- /dev/null +++ b/src/utils.jl @@ -0,0 +1,108 @@ +""" + loadMagneticFieldMeasurementData(filename::String) + +Load magnetic field measurement data by choosing the correct function for the version of the file and return all required data. +""" +function loadMagneticFieldMeasurementData(filename::String) + + measurementData = h5open(filename, "r") do file + + # test if the file contains the correct data + if haskey(file, "positions") && file["positions"] isa HDF5.Group # newest file structure + measurementData = loadMagneticFieldMeasurementDataV2(filename) + elseif haskey(file, "positions") && !(file["positions"] isa HDF5.Group)# oldest file structure + measurementData = loadMagneticFieldMeasurementDataV1(filename) + else # file structure not known + throw(ErrorException("Unknown file structure.")) + end + end + + # get required data + field = measurementData["fields"] + # positions + tDes = measurementData["positions"]["tDesign"] + radius = tDes["radius"] + N = tDes["N"] + t = tDes["t"] + center = tDes["center"] + + if haskey(measurementData, "sensor") + correction = measurementData["sensor"]["correctionTranslation"] + else + correction = Matrix{Float64}(I, 3, 3) + @warn "No correction of the sensor translations found. Using identity matrix." + end + + return field, radius, N, t, center, correction +end + +""" + loadMagneticFieldMeasurementDataV2(filename::String) + +Load magnetic field measurement data from hdf5-file (version v2). +""" +function loadMagneticFieldMeasurementDataV2(filename::String) + + measurementData = h5open(filename, "r") do file + measurementData = Dict{String, Any}() + + # magnetic field + measurementData["fields"] = read(file, "/fields") # measured field (size: 3 x #points x #patches) + + # positions + measurementData["positions"] = Dict{String, Any}() + measurementData["positions"]["tDesign"] = Dict{String, Any}() + tDes = measurementData["positions"]["tDesign"] + tDes["radius"] = read(file, "/positions/tDesign/radius") # radius of the measured ball + tDes["N"] = read(file, "/positions/tDesign/N") # number of points of the t-design + tDes["t"] = read(file, "/positions/tDesign/t") # t of the t-design + tDes["center"] = read(file, "/positions/tDesign/center") # center of the measured ball + + # sensor + measurementData["sensor"] = Dict{String, Any}() + measurementData["sensor"]["correctionTranslation"] = read(file, "/sensor/correctionTranslation") + + # optional data + if haskey(file, "currents") + measurementData["currents"] = read(file, "currents") # currents used for the measurement + end + + return measurementData + end + + return measurementData +end + +""" + loadMagneticFieldMeasurementDataV1(filename::String) + +Load magnetic field measurement data from hdf5-file (version v1). +""" +function loadMagneticFieldMeasurementDataV1(filename::String) + + measurementData = h5open(filename, "r") do file + measurementData = Dict{String, Any}() + + measurementData["fields"] = read(file, "fields") # measured field (size: 3 x #points x #patches) + measurementData["fieldsError"] = read(file,"fieldsError") # field error stemming from the measurement device + + # positions + measurementData["positions"] = Dict{String, Any}() + measurementData["positions"]["tDesign"] = Dict{String, Any}() + tDes = measurementData["positions"]["tDesign"] + tDes["positions"] = read(file,"positions") # measured positions (shifted and scaled t-design) + tDes["radius"] = read(file, "positionsTDesignRadius") # radius of the measured ball + tDes["N"] = read(file, "positionsTDesignN") # number of points of the t-design + tDes["t"] = read(file, "positionsTDesignT") # t of the t-design + tDes["center"] = read(file, "positionsCenter") # center of the measured ball + + # optional data + if haskey(file, "currents") + measurementData["currents"] = read(file, "currents") # currents used for the measurement + end + + return measurementData + end + + return measurementData +end \ No newline at end of file