Biot Savart in Julia

using LinearAlgebra
using Interpolations
using StatsBase
using NPZ

function createXZGrid(bbox, meshsize)
    x1, z1, x2, z2 = bbox
    xgrid, zgrid = meshgrid(range(x1, x2, length=meshsize), range(z1, z2, length=meshsize))

    zeroes = zeros(size(xgrid))

    # create the individual position vectors
    vecs = cat(reshape(xgrid, :, 1), reshape(zeroes, :, 1), reshape(zgrid, :, 1), dims=2)
    vecs = reshape(vecs, size(xgrid)..., 3)
    return xgrid, zgrid, vecs
end

function meshgrid(x, y)
    X = [i for j in y, i in x]
    Y = [j for j in y, i in x]
    return X, Y
end

function wireIntegrate(meshgrid::Array{Float64,3}, l::Matrix{Float64})


    dl = diff(l, dims=1)

    vectorIntegral = zeros(Float32, size(meshgrid))
    integrand = zeros(Float32, size(meshgrid))

    for (lv, dlv) in zip(eachrow(l[1:end-1,:]), eachrow(dl))
        rp = meshgrid .- reshape(lv, 1, 1, :)

        norms = norm.(eachslice(rp, dims=(1,2)))
        norms = reshape(norms, size(rp)[1:end-1])

        for i=1:size(meshgrid)[1], j=1:size(meshgrid)[2]
            integrand[i,j,:] = cross(dlv, rp[i,j,:]) / norms[i,j]^3
        end

        integrand = replace(integrand, NaN => 0.0)

        vectorIntegral += integrand
    end

    return vectorIntegral
end

function cutoffPercentile(magnitudes::Matrix{Float64}, cutoff::Real)
    mags = copy(magnitudes)

    q = quantile(vec(mags), cutoff / 100) # cutoff must be between 0 and 1.

    mags[mags .> q] .= q

    return mags
end

function createCanvasViaRotation(numCells::Int, newShape::Int, bField::Array{Float64, 3})
    # We create a 3D canvas of size newShape x newShape x newShape by rotating the 2D canvas x-z through 360 degrees and interpolate the values that don't fall on the grid.

    #  Create a 3D canvas of size newShape x newShape x newShape
    cellCanvas   = zeros(Float64, newShape, newShape, newShape, 3)

    for ii in 1:numCells
        for jj in 1:numCells
            for kk in 1:numCells
                d = sqrt((ii - 1)^2 + (jj - 1)^2) # Adjust for 1-based indexing
                if Int(floor(d)) >= numCells
                    break
                end
                if d == floor(d)
                    # The bField is a vertical slice through the 3D field.  The ii,jj
                    # indices are the x and y indices of the 2D canvas.  The kk index
                    # is the z index of the vertical slice of the 3D field.
                    cellCanvas[ii, jj, kk,:] = bField[kk, Int(floor(d)) + 1, :] # Adjust for 1-based indexing
                else
                    xp = [jj - 1, jj] # Adjust for 1-based indexing
                    if d < size(bField, 2) - 1
                        yps0 = [bField[kk, Int(floor(d)) + 1:Int(floor(d)) + 2, i] for i in 1:3] # Adjust for 1-based indexing

                        interpVector  = [LinearInterpolation([floor(d), floor(d) + 1], vs) for vs in yps0]
                        cellCanvas[ii, jj, kk,:] = map.(interpVector, repeat([d], 3))
                    end
                end
            end
        end
    end

    return cellCanvas
end

meshsize = 800  # The size of the mesh grid
I = 1 # 1 amp of current flowing through coil
mu0 = 4 *pi*10^(-7)  # magnetic permeability constant
radius = 1  # radius of coil in meters
# bounding box of visual area in meters
bbox = (0,0, 10, 10)  # bounding box of visual area in meters
xgrid, zgrid, vecs = createXZGrid(bbox, meshsize)
println(size(vecs))

rIdx = Int(meshsize/(bbox[3]-bbox[1]))

angnum = 360 # The number of evenly spaced angle values around a circumference of length 2*pi*R or part thereof
radius = 1.0 # Example radius

stop_value = pi/2
step_size = stop_value / angnum

angles = 0:step_size:(stop_value - step_size) # Calculate step size and create range

coss = cos.(angles)
sins = sin.(angles)
zz = zeros(size(angles))

l = hcat(coss, sins, zz) * radius

println(size(l))



@time vectorIntegral = wireIntegrate(vecs, l)

println("vectorIntegral = ", size(vectorIntegral))

doubleShape = Tuple(2 * s for s in size(vectorIntegral)[1:end-1])
doubleShape = (doubleShape..., size(vectorIntegral)[end])


cellsCanvas = zeros(Float64, doubleShape)
println(size(cellsCanvas))

cellsCanvas[meshsize+1:2*meshsize, meshsize+1:2*meshsize,:] += vectorIntegral
cellsCanvas[meshsize:-1:1, meshsize+1:2*meshsize,:] += vectorIntegral
cellsCanvas[meshsize:-1:1, meshsize:-1:1,:] += vectorIntegral
cellsCanvas[meshsize+1:2*meshsize, meshsize:-1:1,:] += vectorIntegral





# Take the norm of the vector field
magnitudes = dropdims(mapslices(norm, cellsCanvas, dims=3); dims=3)
cutoff = 95.0

result = cutoffPercentile(magnitudes, cutoff)



# Example usage:

bField = vectorIntegral # rand(numCells, numCells, 3)
newShape = size(bField,1)
numCells = newShape
@time result = createCanvasViaRotation(numCells, newShape, bField)

println(size(result))

npzwrite("cellCanvas-jl-$numCells.npy", result)