util_functions.py

import numpy as np

def createXYGrid(bbox, meshsize):
    x1, y1, x2, y2 = bbox
    xgrid, ygrid =np.meshgrid(np.linspace(x1, x2 ,meshsize), np.linspace(y1,y2,meshsize))

    zeroes = np.zeros(xgrid.shape)

    # create the individual position vectors
    vecs = np.stack((xgrid.flatten(), ygrid.flatten(), zeroes.flatten()), axis=-1).reshape(xgrid.shape + (3,))
    return xgrid, ygrid,vecs

def createXZGrid(bbox, meshsize):
    x1, z1, x2, z2 = bbox
    xgrid, zgrid =np.meshgrid(np.linspace(x1, x2 ,meshsize), np.linspace(z1,z2,meshsize))

    zeroes = np.zeros(xgrid.shape)

    # create the individual position vectors
    vecs = np.stack((xgrid.flatten(), zeroes.flatten(), zgrid.flatten()), axis=-1).reshape(xgrid.shape + (3,))
    return xgrid, zgrid, vecs

def createXYZGrid(bbox, meshsize):
    x1, y1, z1, x2, y2, z2 = bbox
    xgrid, ygrid, zgrid =np.meshgrid(np.linspace(x1, x2 ,meshsize), np.linspace(y1,y2,meshsize), np.linspace(z1,z2,meshsize))


    # create the individual position vectors
    vecs = np.stack((xgrid.flatten(), ygrid.flatten(), zgrid.flatten()), axis=-1).reshape(xgrid.shape + (3,))
    return xgrid, ygrid, zgrid, vecs


def wireIntegrate(meshgrid, l, closed=True):

    def process(meshgrid, vectorIntegral, l, dl):
        for lv, dlv in zip(l, dl):

            rp = meshgrid - lv



            # We want to make all the vectors unit vectors to have a magnitude of 1 but numpy doesn't seem to like how I am trying to do it
            # so I am duplicating the norm entries 3 times to make the matrices the same size.
            #
            # Because a nxn matrix is not broadcastable via a nxnx3 matrix make the smaller into a nxnx3 matrix by repeating
            # it's entries 3 times

            # Note we use the ellipsis notation here (below) to make this code generic enough to handle 2D or 3D grids.
            norms = np.repeat(np.linalg.norm(rp, axis=-1)[...,np.newaxis], 3, axis=-1)

            # Now compute the Biot Savart cross product term
            integrand = np.nan_to_num(np.cross(dlv, rp) /norms**3)

            # Add the dl magnetic contribution to the sum of all vector dl contributions
            vectorIntegral += integrand


    # create an array of the vector dl from l

    if type(l) == np.ndarray:
        dl = np.diff(l, n=1, axis=0)
        if closed:
            dl = np.append(dl, [l[0] - l[-1]], axis=0)
        # The result of the integral summation goes to integrand_vecum as a 2D array of 3-vectors.
        vectorIntegral = np.zeros_like(meshgrid, dtype=np.float32)
        # Now sum up the magnetic field contributions along all the wire lengths dl

        # The numpy way of doing this without looping but it tends to cause the kernel to die
        #lvdlv = np.stack((l, dl), axis=1)

        #rp = meshgrid[:,:,np.newaxis] - lvdlv[:,0,:]
        #norms = np.linalg.norm(rp, axis=-1)
        #cross = I*np.nan_to_num(np.cross(lvdlv[:,1,:], rp) / norms[... , np.newaxis] ** 3)
        #vectorIntegral = cross.sum(axis=2)

        process(meshgrid, vectorIntegral, l, dl)
    elif type(l) == tuple or type(l) == list:

        vectorIntegral = np.zeros_like(meshgrid, dtype=np.float32)
        for l1 in l:
            #print (l1.shape)
            dl = np.diff(l1, n=1, axis=0)
            if closed:
                dl = np.append(dl, [l1[0] - l1[-1]], axis=0)
            process(meshgrid, vectorIntegral, l1, dl)
    return vectorIntegral

def cutoffPercentile(magnitudes, cutoff):
    mags = np.copy(magnitudes)
    # Now this is the only tricky bit.  Because the magnetic values close or in the wire tend to blow up I am chopping those
    # values out and setting them to a maximum value.  This is because B = μ0 I / 2R causes a pole right on top of the wire.  If I don't do
    # this then the visualisation is swamped by the contrast by the extreme high values and the linear scale I am using to match values to grayscale colours.

    q = (np.percentile(mags, cutoff))
    mags[mags > q] = q
    #print("percentile ===", cutoff, " value = ", q)
    return mags

def toGreyScale(magnitudes):

    mags = magnitudes
    # Find the minimum and maximum values in our magnitudes.  Of course the maximum value will be our cutoff value.
    minv = np.amin(mags[~np.isnan(mags)])
    maxv = np.amax(mags[~np.isnan(mags)])

    print(minv, maxv)
    # Produce our array of greyscale values
    #print ("arrint")
    arrint = (((mags - minv)/(maxv - minv)*255).astype(int))
    #print(arrint)
    return arrint