diff --git a/gradientCalculationV2_1.py b/gradientCalculationV2_1.py new file mode 100644 index 0000000..e8bfd09 --- /dev/null +++ b/gradientCalculationV2_1.py @@ -0,0 +1,349 @@ +# -*- coding: utf-8 -*- +""" +Created on Mon Nov 4 17:02:30 2019 + +@author: to_reilly +""" + +import numpy as np +import matplotlib.pyplot as plt +from scipy import special +from scipy import spatial + +def calcP(mm, aa, bb, kk): + mu0 = 4 * np.pi * 10**(-7) + dBesselk_m = special.kvp(mm, abs(kk)*aa,1) + dBesseli_m = special.ivp(mm, abs(kk)*bb,1) + return aa * mu0 * kk * dBesseli_m * dBesselk_m + +def calcQ(mm, aa, bb, kk): + mu0 = 4 * np.pi * 10**(-7) + dBesselk_m = special.kvp(mm,abs(kk)*aa,1) + Besseli_m = special.iv(mm, abs(kk) * bb) + return mm * aa * mu0 / bb * abs(kk) / kk * Besseli_m * dBesselk_m + +def calculateContour(streamF, numWires, phi, z): + + phi2D, z2D= np.meshgrid(phi,z) + levels = np.linspace(np.min(streamF), np.max(streamF), numWires*2 + 4) + levels = levels[1:-1] + + # Wire should be laid along contours between the isolines, calculate midpoint between isolines + midPointLevels = [(levels[i]+levels[i+1])/2 for i in range(np.size(levels)-1)] + midPointLevels = np.array(midPointLevels)[np.abs(midPointLevels) >= 1e-6] #remove zeros, account for floating point error + + plt.ioff() + plt.figure() + contWires = plt.contour(phi2D,z2D,streamF,levels = midPointLevels) + + return contWires + +def halbachXgradient(linearLength = 140, coilRad = 135, coilLength = 350,numWires = 10, numHighOrders = 10, \ + linearityTerm = 16, apoTerm = .05, resolution = 1e-3): + + gradStrength = 1e-3 + a = coilRad*1e-3 #gradRad + b = 0.001*a + d = linearLength*1e-3 #Linear region + nShape = linearityTerm + res = resolution + Zmax = coilLength*1e-3 + N = numHighOrders + h = apoTerm + + Nsamp = np.int(2*Zmax/res) + z = np.linspace(-Zmax,Zmax,Nsamp) + phi = np.linspace(-1.5*np.pi,0.5*np.pi,Nsamp) + + + kfov = 1/res + k = np.linspace(0.00001,kfov,Nsamp).conj().T + + gradShape = np.divide(z,1+(z/d)**nShape) + gradShape_k = np.fft.fft(gradShape.conj().T); + + t_k = np.exp(-2*(h*k)**2) #apodisation term + + n_0 = 2*np.pi*gradShape_k*gradStrength/(calcQ(1,a,b,k)+calcP(1,a,b,k)) + streamF = 0 + + for n in range(N+1): + sign = (-1)**(n+1) + scale = 1 + for m in range(n+1): + scale *= np.divide((calcP(2*m-1,a,b,k)-calcQ(2*m-1,a,b,k)),(calcP(2*m+1,a,b,k)+calcQ(2*m+1,a,b,k))) + B_apo = sign*np.fft.ifft(np.divide(np.multiply(n_0,t_k*scale),k)) + streamF += np.outer(B_apo,np.cos((2*n+1)*phi)) + streamF = np.real(streamF) + + + return calculateContour(streamF, numWires, phi, z) + +def halbachYgradient(linearLength = 140, coilRad = 135, coilLength = 350, numWires = 10, numHighOrders = 10, \ + linearityTerm = 16, apoTerm = .05, resolution = 1e-3): + + gradStrength = 1e-3 + a = coilRad*1e-3 #gradRad + b = 0.001*a + d = linearLength*1e-3 #Linear region + Zmax = 2*coilLength*1e-3 + N = 0 + print("Number of higher order modes set to 0, not needed for Y and Z") + h = apoTerm + + Nsamp = np.int(2*Zmax/resolution) + z = np.linspace(-Zmax,Zmax,Nsamp) + phi = np.linspace(-np.pi,np.pi,Nsamp) + + k = np.linspace(0.00001,1/resolution,Nsamp).conj().T + + gradShape = 1/(1+(z/d)**linearityTerm) - 1/(1+((z+3.5*d)/(0.5*d))**linearityTerm)-1/(1+((z-3.5*d)/(0.5*d))**linearityTerm) + + gradShape_k = np.fft.fft(gradShape.conj().T); + + t_k = np.exp(-2*(h*k)**2) #apodisation term + + n_0 = b*2*np.pi*gradShape_k*gradStrength/(calcQ(2,a,b,k)+calcP(2,a,b,k)) + streamF = 0 + + for n in range(N+1): + amp = (-1)**(n+1) + scale = 1 + for m in range(1,n+1): + scale *= np.divide((calcP(2*m,a,b,k)-calcQ(2*m,a,b,k)),(calcP(2*m+2,a,b,k)+calcQ(2*m+2,a,b,k))) + B_apo = np.fft.ifft((2/np.pi)*amp*(np.divide(1,k))*n_0*scale*t_k) + streamF += np.outer(B_apo,np.sin((2*n+2)*phi)) + + #remove sidelobes + d_samp=d/resolution + streamF = streamF[int(Nsamp/2-1.5*d_samp):int(Nsamp/2+1.5*d_samp)] + z = z[int(Nsamp/2-1.5*d_samp):int(Nsamp/2+1.5*d_samp)] + phi2D, z2D= np.meshgrid(phi,z) + + return calculateContour(streamF.real, numWires, phi, z) + +def halbachZgradient(linearLength = 140, coilRad = 135, coilLength = 350, numWires = 10, numHighOrders = 10, \ + linearityTerm = 16, apoTerm = .05, resolution = 1e-3): + + gradStrength = 1e-3 + a = coilRad*1e-3 #gradRad + b = 0.001*a + d = linearLength*1e-3 #Linear region + Zmax = 2*coilLength*1e-3 + N = 0 + #print("Number of higher order modes set to 0, not needed for Y and Z") + h = apoTerm + + Nsamp = np.int(2*Zmax/resolution) + z = np.linspace(-Zmax,Zmax,Nsamp) + phi = np.linspace(-np.pi,np.pi,Nsamp)+np.pi/4 + + k = np.linspace(0.00001,1/resolution,Nsamp).conj().T + + gradShape = 1/(1+(z/d)**linearityTerm) - 1/(1+((z+3.5*d)/(0.5*d))**linearityTerm)-1/(1+((z-3.5*d)/(0.5*d))**linearityTerm) + gradShape_k = np.fft.fft(gradShape.conj().T); + + t_k = np.exp(-2*(h*k)**2) #apodisation term + + n_0 = b*2*np.pi*gradShape_k*gradStrength/(calcQ(2,a,b,k)+calcP(2,a,b,k)) + streamF = 0 + + for n in range(N+1): + amp = (-1)**(n+1) + scale = 1 + for m in range(1,n+1): + scale *= np.divide((calcP(2*m,a,b,k)-calcQ(2*m,a,b,k)),(calcP(2*m+2,a,b,k)+calcQ(2*m+2,a,b,k))) + B_apo = np.fft.ifft((2/np.pi)*amp*(np.divide(1,k))*n_0*scale*t_k) + streamF += np.outer(B_apo,np.cos((2*n+2)*phi)) + + #remove sideloabs + d_samp=d/resolution + streamF = streamF[int(Nsamp/2-1.5*d_samp):int(Nsamp/2+1.5*d_samp)] + streamF[0] = 0 + streamF[-1] = 0 + z = z[int(Nsamp/2-1.5*d_samp):int(Nsamp/2+1.5*d_samp)] + phi2D, z2D= np.meshgrid(phi,z) + np.save("streamF_z.npy", streamF) + np.save("phi2D.npy", streamF) + np.save("z2D.npy", z2D) + + return calculateContour(streamF.real, numWires, phi, z) + +def calculateBfield(contours, DSV, resolution, coilRad, direction): + import numexpr as ne + + radius = np.float32(DSV/2) + + phi = np.linspace(0, 2*np.pi, int(2*np.pi*radius/resolution), dtype = np.float32) + theta = np.linspace(0, np.pi, int(np.pi*radius/resolution), dtype = np.float32) + phiGrid, thetaGrid = np.meshgrid(phi,theta) + xSphere = radius*np.multiply(np.sin(thetaGrid), np.cos(phiGrid)) + ySphere = radius*np.multiply(np.sin(thetaGrid), np.sin(phiGrid)) + zSphere = radius*np.cos(thetaGrid) + + points = np.stack((np.ravel(xSphere),np.ravel(ySphere),np.ravel(zSphere)),axis=1) + + wireLevels = contours.allsegs + gradCurrent = np.float32(contours.levels[1] - contours.levels[0]) + + bField = np.zeros(np.shape(points)[0], dtype = np.float32) + + import time + startTime = time.time() + wireCounter = 1 + for wireLevel in wireLevels: + for wire in wireLevel: + print("Simulating wire %.0f of %0.f"%(wireCounter, np.size(wireLevels))) + wire = np.array(wire, dtype = np.float32) + wirePath3D = np.stack((np.cos(wire[:,0])*np.float32(coilRad),np.sin(wire[:,0])*np.float32(coilRad),wire[:,1]),axis=1) + idS = 1e-7*gradCurrent*(wirePath3D[1:,:] - wirePath3D[:-1,:]) + r = points[:,np.newaxis] - wirePath3D[:-1,:] + r3 = np.sum(np.square(r), axis = 2)[:,:,np.newaxis] + rNorm = r/(r3*np.sqrt(r3)) + bField += np.matmul(rNorm[:,:,2], idS[:,1]) - np.matmul(rNorm[:,:,1],idS[:,2]) + wireCounter += 1 + print("Execution time: %.2f seconds"%(time.time()-startTime)) + error = calculateError(points, bField, direction) + return [xSphere*1e3, ySphere*1e3, zSphere*1e3], np.reshape(bField, (np.size(theta), np.size(phi))), error + +def calculateError(coords, bField, direction): + if(direction == 0): + coordAxis = coords[:,2] + elif(direction == 1): + coordAxis = coords[:,1] + else: + coordAxis = coords[:,0] + + argMin = np.argmin(coordAxis) + argMax = np.argmax(coordAxis) + + posRange = np.max(coordAxis) - np.min(coordAxis) + + bRange = bField[argMax] - bField[argMin] + + efficiency = bRange/posRange + + coordAxis[np.abs(coordAxis) < 0.01*posRange] = 'NaN' + + return np.nanmax((bField - efficiency*coordAxis)/(efficiency*coordAxis)) + +def exportWires(contours, coilRad, direction, conjoined): + + wireNum = 0 + contourDict = {} + wireLevels = contours.allsegs + + if ((direction == 0) and conjoined): #for the X gradient the center of the smallest contour is needed for joining the wires + minLength = np.inf + for wireLevel in wireLevels: + for wire in wireLevel: + if(np.size(wire,0) < minLength): + centerHeight = np.abs(np.mean(wire[:,1])*1e3) + + for wireLevel in wireLevels: + for wire in wireLevel: + wirePath3D = np.stack((np.cos(wire[:,0])*coilRad,np.sin(wire[:,0])*coilRad,wire[:,1]*1e3),axis=1) + if(conjoined): + gapSize = 8 #gap in which the sections are joined + gapAngle = gapSize/coilRad + centerAngle = np.mean(wire[:,0]) + + if(direction == 0): + mask = (np.abs(wire[:,0] - centerAngle) > gapAngle) | (np.abs(wirePath3D[:,2]) < centerHeight) + else: + mask = (np.abs(wire[:,0] - centerAngle) > gapAngle) | (wirePath3D[:,2] < 0) + + while mask[0]: + mask = np.roll(mask,1) + wirePath3D = np.roll(wirePath3D, 1, axis = 0) + + contourDict[str(wireNum)] = np.stack((wirePath3D[mask, 0],wirePath3D[mask, 1],wirePath3D[mask, 2]),axis=1) + else: + contourDict[str(wireNum)] = wirePath3D + wireNum += 1 + + if(not conjoined): + return contourDict + else: + + ############################################# + # Join the wires with a gap in to one array # + ############################################# + + numCoilSegments = 4 #Number of quadrants + + joinedContour = {} + joinedContour[str(0)] = contourDict[str(0)] + joinedContour[str(1)] = contourDict[str(1)] + joinedContour[str(2)] = contourDict[str(int(2*wireNum/numCoilSegments))] + joinedContour[str(3)] = contourDict[str(int(2*wireNum/numCoilSegments)+1)] + + for idx in range(1,int(wireNum/numCoilSegments)): + joinedContour[str(0)] = np.append(joinedContour[str(0)], contourDict[str(2*idx)], axis = 0) + joinedContour[str(1)] = np.append(joinedContour[str(1)], contourDict[str(2*idx+1)], axis = 0) + joinedContour[str(2)] = np.append(joinedContour[str(2)], contourDict[str(int(2*wireNum/numCoilSegments) + idx*2 )], axis = 0) + joinedContour[str(3)] = np.append(joinedContour[str(3)], contourDict[str(int(2*wireNum/numCoilSegments) + idx*2 +1)], axis = 0) + + + ############################################ + # Check for consecutive identical elements # + ############################################ + tol = 1e-5 + for key in joinedContour: + delta = joinedContour[key][1:,:] - joinedContour[key][:-1,:] + delta = np.sum(np.square(delta), axis = 1) + zeroElements = delta < tol + joinedContour[key] = np.delete(joinedContour[key],np.nonzero(zeroElements), axis = 0) + + return joinedContour + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/gradientDesignTool.py b/gradientDesignTool.py index 8177e41..9564ac5 100644 --- a/gradientDesignTool.py +++ b/gradientDesignTool.py @@ -7,11 +7,9 @@ import numpy as np import tkinter as tk -from tkinter import ttk, font -from tkinter.ttk import Style +from tkinter import ttk import os -import gradientCalculation as gradCalc - +import gradientCalculationV2_1 as gradCalc from matplotlib.backends.backend_tkagg import (FigureCanvasTkAgg) from matplotlib.figure import Figure @@ -19,7 +17,6 @@ from mpl_toolkits.mplot3d import Axes3D from matplotlib import cm - class PopUpTextMessage(tk.Toplevel): def __init__(self, parent, textmsg): super(PopUpTextMessage, self).__init__(parent) @@ -34,358 +31,279 @@ def __init__(self, parent, textmsg): def cleanup(self): self.destroy() - class GDT_GUI(tk.Tk): - - def __init__(self, *args, **kwargs): + + def __init__ (self, *args, **kwargs): super(GDT_GUI, self).__init__(*args, **kwargs) # set name and icon # super().iconbitmap(self, default=None) - #super().configure(background='gray') super().wm_title('Gradient design tool') - - - '''INPUT FRAME''' - #s = Style() - #s.configure('My.TFrame', background='gray') - inputFrame = ttk.Frame(self, borderwidth=1)#,style='My.TFrame') - inputFrame.grid(row=0, column=0, columnspan=1, sticky=tk.W, - padx=10, pady=10) - - default_font = font.nametofont("TkDefaultFont") - default_font.configure(size=12) - - inputDescriptionLabel = ttk.Label(inputFrame, text='Design parameters', - font='Helvetica 18 bold') - inputDescriptionLabel.grid(row=0, column=0, columnspan=2, sticky=tk.W) + inputFrame = ttk.Frame(self,borderwidth=1) + inputFrame.grid(row=0,column=0,columnspan=1, sticky=tk.W, padx=10, pady=10) + + inputDescriptionLabel = ttk.Label(inputFrame, text='Design parameters', font='Helvetica 12 bold') + inputDescriptionLabel.grid(row=0,column=0, columnspan = 2, sticky=tk.W) + + #create variables for storing inputs + self.linearLength = tk.DoubleVar(inputFrame, value = 140) + self.coilRadius = tk.DoubleVar(inputFrame, value = 135) + self.coilLength = tk.DoubleVar(inputFrame, value = 400) + self.numWires = tk.IntVar(inputFrame, value = 14) + self.wireDiameter = tk.DoubleVar(inputFrame, value = 1.5) + self.apodisationTerm = tk.DoubleVar(inputFrame, value = 0.05) + self.linearityOrder = tk.IntVar(inputFrame, value = 16) + self.higherOrderTerms = tk.IntVar(inputFrame, value = 10) + self.simulationSize = tk.DoubleVar(inputFrame, value= 140) + self.simulationRes = tk.DoubleVar(inputFrame, value= 5) + self.DSV = tk.DoubleVar(inputFrame, value = 140) - - # create variables for storing inputs - self.linearLength = tk.DoubleVar(inputFrame, value=140) - self.coilRadius = tk.DoubleVar(inputFrame, value=135) - self.coilLength = tk.DoubleVar(inputFrame, value=400) - self.numWires = tk.IntVar(inputFrame, value=15) - self.wireDiameter = tk.DoubleVar(inputFrame, value=1.5) - self.apodisationTerm = tk.DoubleVar(inputFrame, value=0.05) - self.linearityOrder = tk.IntVar(inputFrame, value=16) - self.higherOrderTerms = tk.IntVar(inputFrame, value=0) - self.simulationSize = tk.DoubleVar(inputFrame, value=140) - self.simulationRes = tk.DoubleVar(inputFrame, value=5) - self.DSV = tk.DoubleVar(inputFrame, value=140) - - directionLabel = ttk.Label(inputFrame, text='Background field direction') - directionLabel.grid(row=1, column=0, pady=5, sticky=tk.W) - self.backgroundInput = ttk.Combobox(inputFrame, values=["Halbach", "Conventional"], - width=12, justify=tk.RIGHT, - state='readonly') - self.backgroundInput.current(0) - self.backgroundInput.grid(row=1, column=1, padx=20) - directionLabel = ttk.Label(inputFrame, text='Gradient direction') - directionLabel.grid(row=2, column=0, pady=5, sticky=tk.W) - self.directionInput = ttk.Combobox(inputFrame, values=["X", "Y", "Z"], - width=7, justify=tk.RIGHT, - state='readonly') + directionLabel.grid(row=1,column=0, pady=5,sticky=tk.W) + self.directionInput = ttk.Combobox(inputFrame, values=["X", "Y", "Z"], width = 7, justify=tk.RIGHT, state = 'readonly') self.directionInput.current(0) - self.directionInput.grid(row=2, column=1, padx=20) - + self.directionInput.grid(row=1,column=1, padx = 20) + targetRegionLabel = ttk.Label(inputFrame, text='Target linear region (mm)') - targetRegionLabel.grid(row=3, column=0, pady=5, sticky=tk.W) - targetRegionInput = ttk.Entry(inputFrame, textvariable=self.linearLength, width=10, justify=tk.RIGHT) - targetRegionInput.grid(row=3, column=1, padx=20) - + targetRegionLabel.grid(row=2,column=0, pady=5,sticky=tk.W) + targetRegionInput = ttk.Entry(inputFrame, textvariable=self.linearLength, width = 10, justify=tk.RIGHT) + targetRegionInput.grid(row=2,column=1,padx = 20) + coilRadiusLabel = ttk.Label(inputFrame, text='Coil radius (mm)') - coilRadiusLabel.grid(row=4, column=0, pady=5, sticky=tk.W) - coilRadiusInput = ttk.Entry(inputFrame, textvariable=self.coilRadius, width=10, justify=tk.RIGHT) - coilRadiusInput.grid(row=4, column=1, padx=20) - + coilRadiusLabel.grid(row=3,column=0, pady=5,sticky=tk.W) + coilRadiusInput = ttk.Entry(inputFrame, textvariable=self.coilRadius, width = 10, justify=tk.RIGHT) + coilRadiusInput.grid(row=3,column=1,padx = 20) + coilLengthLabel = ttk.Label(inputFrame, text='Coil length (mm)') - coilLengthLabel.grid(row=5, column=0, pady=5, sticky=tk.W) - coilLengthInput = ttk.Entry(inputFrame, textvariable=self.coilLength, width=10, justify=tk.RIGHT) - coilLengthInput.grid(row=5, column=1, padx=20) - + coilLengthLabel.grid(row=4,column=0, pady=5,sticky=tk.W) + coilLengthInput = ttk.Entry(inputFrame, textvariable=self.coilLength, width = 10, justify=tk.RIGHT) + coilLengthInput.grid(row=4,column=1,padx = 20) + numWiresLabel = ttk.Label(inputFrame, text='Wires per quadrant') - numWiresLabel.grid(row=6, column=0, pady=5, sticky=tk.W) - numWiresInput = ttk.Entry(inputFrame, textvariable=self.numWires, width=10, justify=tk.RIGHT) - numWiresInput.grid(row=6, column=1, padx=20) - - - # #Todo: remove higher order terms for longidutinal background fields - # numHigherTermsLabel = ttk.Label(inputFrame, text='Number of higher order terms') - # numHigherTermsLabel.grid(row=8, column=0, pady=5, sticky=tk.W) - # numHigherTermsInput = ttk.Entry(inputFrame, textvariable=self.higherOrderTerms, width=10, justify=tk.RIGHT) - # numHigherTermsInput.grid(row=8, column=1, padx=20) - - resolutionLabel = ttk.Label(inputFrame, text='Linearity order') - resolutionLabel.grid(row=7, column=0, pady=5, sticky=tk.W) - resolutionInput = ttk.Entry(inputFrame, textvariable=self.linearityOrder, width=10, justify=tk.RIGHT) - resolutionInput.grid(row=7, column=1, padx=20) - - apodisationTermLabel = ttk.Label(inputFrame, text='Apodisation term') - apodisationTermLabel.grid(row=8, column=0, pady=5, sticky=tk.W) - apodisationTermInput = ttk.Entry(inputFrame, textvariable=self.apodisationTerm, width=10, justify=tk.RIGHT) - apodisationTermInput.grid(row=8, column=1, padx=20) + numWiresLabel.grid(row=5,column=0, pady=5,sticky=tk.W) + numWiresInput = ttk.Entry(inputFrame, textvariable=self.numWires, width = 10, justify=tk.RIGHT) + numWiresInput.grid(row=5,column=1,padx = 20) - - '''SIMULATION FRAME''' - - inputDescriptionLabel = ttk.Label(inputFrame, text='Simulation parameters', - font='Helvetica 18 bold') - inputDescriptionLabel.grid(row=9, column=0, columnspan=2, sticky=tk.W) - - wireDiameterLabel = ttk.Label(inputFrame, text='Wire diameter (mm)') - wireDiameterLabel.grid(row=10, column=0, pady=5, sticky=tk.W) - wireDiameterInput = ttk.Entry(inputFrame, textvariable=self.wireDiameter, width=10, justify=tk.RIGHT) - wireDiameterInput.grid(row=10, column=1, padx=20) - - linearityOrderLabel = ttk.Label(inputFrame, text='Simulation volume (mm)') - linearityOrderLabel.grid(row=11, column=0, pady=5, sticky=tk.W) - linearityOrderInput = ttk.Entry(inputFrame, textvariable=self.simulationSize, width=10, justify=tk.RIGHT) - linearityOrderInput.grid(row=11, column=1, padx=20) + wireDiameterLabel.grid(row=6,column=0, pady=5,sticky=tk.W) + wireDiameterInput = ttk.Entry(inputFrame, textvariable=self.wireDiameter, width = 10, justify=tk.RIGHT) + wireDiameterInput.grid(row=6,column=1,padx = 20) + + numHigherTermsLabel = ttk.Label(inputFrame, text='Number of higher order terms') + numHigherTermsLabel.grid(row=7,column=0, pady=5,sticky=tk.W) + numHigherTermsInput = ttk.Entry(inputFrame, textvariable=self.higherOrderTerms, width = 10, justify=tk.RIGHT) + numHigherTermsInput.grid(row=7,column=1,padx = 20) + + linearityOrderLabel = ttk.Label(inputFrame, text='Linearity order') + linearityOrderLabel.grid(row=8,column=0, pady=5,sticky=tk.W) + linearityOrderInput = ttk.Entry(inputFrame, textvariable=self.linearityOrder, width = 10, justify=tk.RIGHT) + linearityOrderInput.grid(row=8,column=1,padx = 20) + + apodisationTermLabel = ttk.Label(inputFrame, text='Apodisation term') + apodisationTermLabel.grid(row=9,column=0, pady=5,sticky=tk.W) + apodisationTermInput = ttk.Entry(inputFrame, textvariable=self.apodisationTerm, width = 10, justify=tk.RIGHT) + apodisationTermInput.grid(row=9,column=1,padx=20) - simDimensionsLabel = ttk.Label(inputFrame, text='Simulation resolution (mm)') - simDimensionsLabel.grid(row=12, column=0, pady=5, sticky=tk.W) - simDimensionsInput = ttk.Entry(inputFrame, textvariable=self.simulationRes, width=10, justify=tk.RIGHT) - simDimensionsInput.grid(row=12, column=1, padx=20) + DSVLabel = ttk.Label(inputFrame, text='Simulation DSV (mm)') + DSVLabel.grid(row=10,column=0, pady=5,sticky=tk.W) + DSVInput = ttk.Entry(inputFrame, textvariable=self.DSV, width = 10, justify=tk.RIGHT) + DSVInput.grid(row=10,column=1,padx = 20) - DSVLabel = ttk.Label(inputFrame, text='Error DSV (mm)') - DSVLabel.grid(row=13, column=0, pady=5, sticky=tk.W) - DSVInput = ttk.Entry(inputFrame, textvariable=self.DSV, width=10, justify=tk.RIGHT) - DSVInput.grid(row=13, column=1, padx=20) + simResolutionLabel = ttk.Label(inputFrame, text='Simulation resolution (mm)') + simResolutionLabel.grid(row=11,column=0, pady=5,sticky=tk.W) + simResolutionInput = ttk.Entry(inputFrame, textvariable=self.simulationRes, width = 10, justify=tk.RIGHT) + simResolutionInput.grid(row=11,column=1,padx = 20) calculateWireButton = ttk.Button(inputFrame, text="Calculate wire pattern", command=self.calculateGradientWires) - calculateWireButton.grid(row=14, column=0, pady=5, sticky=tk.W) + calculateWireButton.grid(row=12, column=0, pady=5, sticky=tk.W) self.calculateFieldButton = ttk.Button(inputFrame, text="Simulate B-field", command=self.calculateGradientField) - self.calculateFieldButton.grid(row=14, column=1, padx=20) + self.calculateFieldButton.grid(row=12, column=1, padx=20) ''' OUTPUT FRAME''' outputFrame = ttk.Frame(self) - outputFrame.grid(row=1, column=0, columnspan=1, sticky=tk.NW, padx=10, pady=10) - - outputDescriptionLabel = ttk.Label(outputFrame, text='Simulation output', font='Helvetica 18 bold') - outputDescriptionLabel.grid(row=0, column=0, columnspan=1, sticky=tk.W) - - self.gradEfficiencyString = tk.StringVar(outputFrame, value="-") - self.gradErrorString = tk.StringVar(outputFrame, value="-") - self.resistanceString = tk.StringVar(outputFrame, value="-") - self.inductanceString = tk.StringVar(outputFrame, value="-") - self.wireLengthString = tk.StringVar(outputFrame, value="-") - self.DSVerrorString = tk.StringVar(outputFrame, value="Error over {:.0f} mm DSV:".format(self.DSV.get())) - self.exportConjoined = tk.BooleanVar(outputFrame, value=True) - + outputFrame.grid(row=1,column=0,columnspan=1, sticky=tk.NW, padx=10, pady=20) + + outputDescriptionLabel = ttk.Label(outputFrame, text='Simulation output', font='Helvetica 12 bold') + outputDescriptionLabel.grid(row=0,column=0, columnspan = 1, sticky=tk.W) + + self.gradEfficiencyString = tk.StringVar(outputFrame, value = "-") + self.gradErrorString = tk.StringVar(outputFrame, value = "-") + self.resistanceString = tk.StringVar(outputFrame, value = "-") + self.inductanceString = tk.StringVar(outputFrame, value = "-") + self.wireLengthString = tk.StringVar(outputFrame, value = "-") + self.exportConjoined = tk.BooleanVar(outputFrame, value = True) + self.zRangeString = tk.StringVar(outputFrame, value = "-") + efficiencyTextLabel = ttk.Label(outputFrame, text='Gradient efficiency:') - efficiencyTextLabel.grid(row=1, column=0, pady=5, sticky=tk.W) - efficiencyValueLabel = ttk.Label(outputFrame, textvariable=self.gradEfficiencyString, justify=tk.RIGHT) - efficiencyValueLabel.grid(row=1, column=1, padx=10, sticky=tk.E) - - # TODO: This doesn't update for a change in DSV value! - linearityTextLabel = ttk.Label(outputFrame, textvariable=self.DSVerrorString) - linearityTextLabel.grid(row=2, column=0, pady=5, sticky=tk.W) - linearityValueLabel = ttk.Label(outputFrame, textvariable=self.gradErrorString, justify=tk.RIGHT) - linearityValueLabel.grid(row=2, column=1, padx=10, sticky=tk.E) + efficiencyTextLabel.grid(row=1,column=0, pady=5,sticky=tk.W) + efficiencyValueLabel = ttk.Label(outputFrame, textvariable=self.gradEfficiencyString, justify = tk.RIGHT) + efficiencyValueLabel.grid(row=1,column=1, padx=10,sticky=tk.E) + + linearityTextLabel = ttk.Label(outputFrame, text='Error over %.0f mm DSV:'%(self.DSV.get())) + linearityTextLabel.grid(row=2,column=0, pady=5,sticky=tk.W) + linearityValueLabel = ttk.Label(outputFrame, textvariable=self.gradErrorString, justify = tk.RIGHT) + linearityValueLabel.grid(row=2,column=1, padx=10,sticky=tk.E) wireLengthTextLabel = ttk.Label(outputFrame, text='Wire length:') - wireLengthTextLabel.grid(row=3, column=0, pady=5, sticky=tk.W) - wireLengthValueLabel = ttk.Label(outputFrame, textvariable=self.wireLengthString, justify=tk.RIGHT) - wireLengthValueLabel.grid(row=3, column=1, padx=10, sticky=tk.E) + wireLengthTextLabel.grid(row=3,column=0, pady=5,sticky=tk.W) + wireLengthValueLabel = ttk.Label(outputFrame, textvariable=self.wireLengthString, justify = tk.RIGHT) + wireLengthValueLabel.grid(row=3,column=1, padx=10,sticky=tk.E) resistanceTextLabel = ttk.Label(outputFrame, text='Coil resistance:') - resistanceTextLabel.grid(row=4, column=0, pady=5, sticky=tk.W) - resistanceValueLabel = ttk.Label(outputFrame, textvariable=self.resistanceString, justify=tk.RIGHT) - resistanceValueLabel.grid(row=4, column=1, padx=10, sticky=tk.E) - - #inductanceTextLabel = ttk.Label(outputFrame, text='Coil inductance:') - #inductanceTextLabel.grid(row=5, column=0, pady=5, sticky=tk.W) - #inductanceValueLabel = ttk.Label(outputFrame, textvariable=self.inductanceString, justify=tk.RIGHT) - #inductanceValueLabel.grid(row=5, column=1, padx=10, sticky=tk.E) - - self.saveBfieldButton = ttk.Button(outputFrame, text="Export B-field", state=tk.DISABLED, command=self.exportBfield) + resistanceTextLabel.grid(row=4,column=0, pady=5,sticky=tk.W) + resistanceValueLabel = ttk.Label(outputFrame, textvariable=self.resistanceString, justify = tk.RIGHT) + resistanceValueLabel.grid(row=4,column=1, padx=10,sticky=tk.E) + + zRangeStringTextLabel = ttk.Label(outputFrame, text='Min/Max X: ') + zRangeStringTextLabel.grid(row=5,column=0, pady=5,sticky=tk.W) + zRangeStringValueLabel = ttk.Label(outputFrame, textvariable=self.zRangeString, justify = tk.RIGHT) + zRangeStringValueLabel.grid(row=5,column=1, padx=10,sticky=tk.E) + + self.saveBfieldButton = ttk.Button(outputFrame, text="Export B-field",state=tk.DISABLED, command=self.exportBfield) self.saveBfieldButton.grid(row=6, column=0, pady=5, sticky=tk.SW) - + saveWireButton = ttk.Button(outputFrame, text="Export wire to CSV", command=self.exportWireCSV) - saveWireButton.grid(row=6, column=1, pady=5, padx=20, sticky=tk.SE) - - conjoinedExport = ttk.Checkbutton(outputFrame, text="Join wires", variable=self.exportConjoined) - conjoinedExport.grid(row=7, column=1, pady=5, padx=20, sticky=tk.SE) - + saveWireButton.grid(row=6, column=1, pady = 5, padx=20,sticky=tk.SE) + + conjoinedExport = ttk.Checkbutton(outputFrame, text="Join wires", variable = self.exportConjoined) + conjoinedExport.grid(row=7, column=1, pady = 5, padx=20,sticky=tk.SE) + '''Position plots''' - - self.contourfig = Figure(figsize=(7, 7)) + + self.contourfig = Figure(figsize=(5,5)) self.contourfig.gca(projection='3d') self.contourfig.set_tight_layout(True) - - self.contourCanvas = FigureCanvasTkAgg(self.contourfig, master=self) + + self.contourCanvas = FigureCanvasTkAgg(self.contourfig, master = self) self.contourCanvas.draw() - self.contourCanvas.get_tk_widget().grid(row=0, column=1) - - self.fieldFig = Figure(figsize=(7, 7)) + self.contourCanvas.get_tk_widget().grid(row=0,column=1) + self.contourCanvas.mpl_connect('button_press_event', self.contourfig.gca()._button_press) + self.contourCanvas.mpl_connect('button_release_event', self.contourfig.gca()._button_release) + self.contourCanvas.mpl_connect('motion_notify_event', self.contourfig.gca()._on_move) + + self.fieldFig = Figure(figsize=(5,5)) self.fieldFig.gca(projection='3d') self.fieldFig.set_tight_layout(True) - - self.fieldCanvas = FigureCanvasTkAgg(self.fieldFig, master=self) + + self.fieldCanvas = FigureCanvasTkAgg(self.fieldFig, master = self) self.fieldCanvas.draw() - self.fieldCanvas.get_tk_widget().grid(row=0, column=2) - - + self.fieldCanvas.get_tk_widget().grid(row=1,column=1) + self.fieldCanvas.mpl_connect('button_press_event', self.fieldFig.gca()._button_press) + self.fieldCanvas.mpl_connect('button_release_event', self.fieldFig.gca()._button_release) + self.fieldCanvas.mpl_connect('motion_notify_event', self.fieldFig.gca()._on_move) self.calculateGradientWires() - + def calculateGradientWires(self): - - designParameters = {} - designParameters['length'] = self.linearLength.get()*1e-3 - designParameters['linearity'] = self.linearityOrder.get() - - designParameters['resolution'] = self.simulationRes.get()*1e-3 - # designParameters['gridSize'] = self.simulationSize.get()*1e-3 - designParameters['gridSize'] = self.coilLength.get()*1e-3 - - designParameters['B0'] = self.backgroundInput.get() - designParameters['coilRad'] = self.coilRadius.get()*1e-3 - # TODO: arbitrary scaling of b, why this value?? comment and win! - designParameters['fieldRad'] = 0.001*self.coilRadius.get()*1e-3 - designParameters['apo'] = self.apodisationTerm.get() - designParameters['nrModes'] = self.higherOrderTerms.get() - designParameters['nrWires'] = self.numWires.get() - designParameters['wireDiameter'] = self.wireDiameter.get()*1e-3 - - - # Halbach system: - if (self.backgroundInput.current() == 0): - - if (self.directionInput.current() == 0): - designParameters['type'] = "longitudinal" - designParameters['gradDir'] = "z" - - elif(self.directionInput.current() == 2): - designParameters['type'] = "transverse" - designParameters['gradDir'] = "x" - - elif(self.directionInput.current() == 1): - designParameters['type'] = "transverse" - designParameters['gradDir'] = "y" - - # Conventional MRI: - elif (self.backgroundInput.current() == 1): + + if (self.directionInput.current() == 0): + self.contourData = gradCalc.halbachXgradient(linearLength = self.linearLength.get(), coilRad = self.coilRadius.get(), coilLength = self.coilLength.get(),\ + numWires = self.numWires.get(), numHighOrders = self.higherOrderTerms.get(), linearityTerm = self.linearityOrder.get(),\ + apoTerm = self.apodisationTerm.get()) - if (self.directionInput.current() == 0): - designParameters['type'] = "transverse" - designParameters['gradDir'] = "x" - - elif(self.directionInput.current() == 1): - designParameters['type'] = "transverse" - designParameters['gradDir'] = "y" - - elif(self.directionInput.current() == 2): - designParameters['type'] = "longitudinal" - designParameters['gradDir'] = "z" + elif(self.directionInput.current() == 1): + self.contourData = gradCalc.halbachYgradient(linearLength = self.linearLength.get(), coilRad = self.coilRadius.get(), coilLength = self.coilLength.get(),\ + numWires = self.numWires.get(), numHighOrders = self.higherOrderTerms.get(), linearityTerm = self.linearityOrder.get(),\ + apoTerm = self.apodisationTerm.get()) - - - Coil = gradCalc.generateCoil(designParameters, gradStrength=1e-3, current=2) - - self.contourData = Coil['wires'] + elif(self.directionInput.current() == 2): + self.contourData = gradCalc.halbachZgradient(linearLength = self.linearLength.get(), coilRad = self.coilRadius.get(), coilLength = self.coilLength.get(),\ + numWires = self.numWires.get(), numHighOrders = self.higherOrderTerms.get(), linearityTerm = self.linearityOrder.get(),\ + apoTerm = self.apodisationTerm.get()) wireLevels = self.contourData.allsegs - self.inductanceString.set("{:.4g} \u03BCH".format(Coil['L']*1e6)) + self.gradEfficiencyString.set("%.4f mT/m/A"%(np.divide(1,self.contourData.levels[1] - self.contourData.levels[0]))) - self.contourfig.gca().clear() - - if (self.backgroundInput.current() == 0): - self.contourfig.gca().set_xlabel('Z/B$_{0}$ (mm)') - self.contourfig.gca().set_ylabel('X (mm)') - self.contourfig.gca().set_zlabel('Y (mm)') - else: - self.contourfig.gca().set_xlabel('X (mm)') - self.contourfig.gca().set_ylabel('Z/B$_{0}$ (mm)') - self.contourfig.gca().set_zlabel('Y (mm)') - - + self.contourfig.gca().set_xlabel('Z (mm)') + self.contourfig.gca().set_ylabel('Y (mm)') + self.contourfig.gca().set_zlabel('X (mm)') + np.warnings.filterwarnings('ignore', category=np.VisibleDeprecationWarning) + self.length = 0 for idx, wireLevel in enumerate(wireLevels): currentLevel = self.contourData.levels[idx] - for wire in range(np.size(wireLevel, 0)): - wirePath3D = np.stack((np.cos(wireLevel[wire][:, 0])*self.coilRadius.get(), - np.sin(wireLevel[wire][:, 0])*self.coilRadius.get(), - wireLevel[wire][:, 1]*1e3), axis=1) - if currentLevel > 0: - self.contourfig.gca().plot3D(wirePath3D[:, 0], - wirePath3D[:, 2], - wirePath3D[:, 1], 'red') + for wire in range (np.size(wireLevel,0)): + wirePath3D = np.stack((np.cos(wireLevel[wire][:,0])*self.coilRadius.get(), np.sin(wireLevel[wire][:,0])*self.coilRadius.get(),wireLevel[wire][:,1]*1e3),axis=1).astype(float) + if currentLevel>0: + self.contourfig.gca().plot3D(wirePath3D[:,0],wirePath3D[:,1],wirePath3D[:,2],'red') else: - self.contourfig.gca().plot3D(wirePath3D[:, 0], - wirePath3D[:, 2], - wirePath3D[:, 1], 'blue') - - self.length = Coil['length'] - self.wireLengthString.set("{:.2f} metres".format(self.length)) - self.contourfig.gca().mouse_init() - self.contourCanvas.draw() - self.contourCanvas.flush_events() - self.DSVerrorString.set('Error over {:.0f} mm DSV:'.format(self.DSV.get())) - self.resistanceString.set("{:.4f} \u03A9".format(Coil['R'])) - - def calculateGradientField(self): - self.calculateFieldButton['state'] = 'disabled' - super().update() + self.contourfig.gca().plot3D(wirePath3D[:,0],wirePath3D[:,1],wirePath3D[:,2],'blue') - if (self.backgroundInput.current() == 0): - self.fieldFig.gca().set_xlabel('Z/B$_{0}$ (mm)') - self.fieldFig.gca().set_ylabel('X (mm)') - self.fieldFig.gca().set_zlabel('Y (mm)') - else: - self.fieldFig.gca().set_xlabel('X (mm)') - self.fieldFig.gca().set_ylabel('Z/B$_{0}$ (mm)') - self.fieldFig.gca().set_zlabel('Y (mm)') + temp, indices = np.unique(wirePath3D.round(decimals=6), axis = 0, return_index = True) + wirePath3D = wirePath3D[sorted(indices)] + delta = wirePath3D[1:,:] - wirePath3D[:-1,:] + segLength = np.sqrt(np.sum(np.square(delta), axis = 1)) + self.length += np.sum(segLength)# + np.sum(np.square(wirePath3D[0,:] - wirePath3D[-1,:])) - self.coords, self.bField, self.efficiency, error= gradCalc.calculateBfield(self.contourData, self.DSV.get()*1e-3, self.simulationRes.get()*1e-3, self.coilRadius.get()*1e-3, self.directionInput.current(), self.backgroundInput.get()) + # self.contourfig.gca().mouse_init() + self.wireLengthString.set("%.2f meters"%(self.length*1e-3)) + self.contourCanvas.draw() + self.contourCanvas.flush_events() + self.wireDict = gradCalc.exportWires(self.contourData, self.coilRadius.get(), self.directionInput.current(), self.exportConjoined.get()) + minZ = 0 + maxZ = 0 + for segmentKey in self.wireDict: + try: + segment = self.wireDict[segmentKey] + minZSegment = np.min(segment[:,2]) + maxZSegment = np.max(segment[:,2]) + if minZSegment < minZ: + minZ = minZSegment + if maxZSegment > maxZ: + maxZ = maxZSegment + except: + print("Error") + + self.zRangeString.set("%.2f / %.2f mm"%(minZ, maxZ)) + self.calculateResistance() + + def calculateResistance(self): + copperResistivity = 1.68*1e-8 + resistance = copperResistivity*self.length*1e-3/(np.pi*(self.wireDiameter.get()*1e-3 /2)**2) + self.resistanceString.set("%.4f Ohms"%(resistance,)) + def calculateGradientField(self): + self.fieldFig.gca().clear() + self.fieldFig.gca().set_xlabel('Z (mm)') + self.fieldFig.gca().set_ylabel('Y (mm)') + self.fieldFig.gca().set_zlabel('X (mm)') + self.coords, self.bField, error = gradCalc.calculateBfield(self.contourData, self.DSV.get()*1e-3, self.simulationRes.get()*1e-3, self.coilRadius.get()*1e-3, self.directionInput.current()) bNorm = (self.bField - np.min(self.bField))/(np.max(self.bField) - np.min(self.bField)) - self.fieldFig.gca().plot_surface(self.coords[0], self.coords[2], self.coords[1], rstride=1, cstride=1, facecolors=cm.jet(bNorm)) - + self.fieldFig.gca().plot_surface(self.coords[0], self.coords[1], self.coords[2], rstride=1, cstride=1, facecolors=cm.jet(bNorm)) + self.fieldFig.gca().mouse_init() self.fieldCanvas.draw() self.fieldCanvas.flush_events() - self.gradEfficiencyString.set("{:.4f} mT/m/A".format(self.efficiency*1000)) - self.gradErrorString.set("{:.2f} %".format(error*100)) - self.fieldFig.gca().mouse_init() - - self.calculateFieldButton['state'] = 'normal' + self.gradErrorString.set("%.2f %%"%(error*100)) self.saveBfieldButton['state'] = 'normal' - super().update() - + def exportBfield(self): - bFieldOutputFile = tk.filedialog.asksaveasfilename(mode='w', defaultextension='.csv', filetypes=(("CSV file", "*.csv"), ("Text file", "*.txt"))) - if bFieldOutputFile is None: + bFieldOutputFile = tk.filedialog.asksaveasfilename(defaultextension = '.csv', filetypes=(("CSV file","*.csv"), ("Text file","*.txt"))) + if bFieldOutputFile == None: return header = 'X (mm),\tY (mm),\tZ (mm),\tBz (mT)' - delimiter = ',\t' - outputArray = np.zeros((np.size(self.bField), np.size(self.coords, 0)+1)) - for idx in range(np.size(self.coords, 0)): - outputArray[:, idx] = np.ravel(self.coords[idx]) - outputArray[:, -1] = np.ravel(self.bField) - np.savetxt(bFieldOutputFile.name, outputArray, delimiter=delimiter, header=header) + delimiter = ',\t' + outputArray = np.zeros((np.size(self.bField),np.size(self.coords,0)+1)) + for idx in range(np.size(self.coords,0)): + outputArray[:,idx] = np.ravel(self.coords[idx]) + outputArray[:,-1] = np.ravel(self.bField) + np.savetxt(bFieldOutputFile.name , outputArray,delimiter = delimiter ,header = header ) bFieldOutputFile.close() - + def exportWireCSV(self): - wireOutputFile = tk.filedialog.asksaveasfilename(defaultextension='.csv', filetypes=(("CSV file","*.csv"), ("Text file", "*.txt"))) + wireOutputFile = tk.filedialog.asksaveasfilename(defaultextension = '.csv', filetypes=(("CSV file","*.csv"), ("Text file","*.txt"))) if wireOutputFile == '': return folder, filename = os.path.split(wireOutputFile) file, extension = os.path.splitext(filename) - - self.wireDict = gradCalc.exportWires(self.contourData, - self.coilRadius.get(), - self.directionInput.current(), - self.exportConjoined.get()) - + + self.wireDict = gradCalc.exportWires(self.contourData, self.coilRadius.get(), self.directionInput.current(), self.exportConjoined.get()) + for key in self.wireDict: - filename = os.path.join(folder, file+key+extension) - np.savetxt(filename, self.wireDict[key], delimiter=",", fmt='%f') - - + filename = os.path.join(folder,file+key+extension) + np.savetxt(filename,self.wireDict[key],delimiter=",", fmt='%f') + + app = GDT_GUI() app.mainloop() -closePlots('all') +closePlots('all') \ No newline at end of file