diff --git a/gradientCalculationV2.py b/gradientCalculationV2.py new file mode 100644 index 0000000..3732b6a --- /dev/null +++ b/gradientCalculationV2.py @@ -0,0 +1,342 @@ +# -*- 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 = numHighOrders + h = apoTerm + + Nsamp = np.int(2*Zmax/resolution) + z = np.linspace(-Zmax,Zmax,Nsamp) + phi = np.linspace(-1*np.pi,1*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, 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 = numHighOrders + h = apoTerm + + Nsamp = np.int(2*Zmax/resolution) + z = np.linspace(-Zmax,Zmax,Nsamp) + phi = np.linspace(-0.75*np.pi,1.25*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.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)] + 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, 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, 2*np.pi*radius/resolution, dtype = np.float32) + theta = np.linspace(0, np.pi, 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 = ne.evaluate("sum(r**2, axis = 2)")[:,:,np.newaxis] + rNorm = ne.evaluate("r/r3**1.5") + 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 new file mode 100644 index 0000000..595b0f6 --- /dev/null +++ b/gradientDesignTool.py @@ -0,0 +1,295 @@ +# -*- coding: utf-8 -*- +""" +Created on Wed Oct 30 12:57:51 2019 + +@author: to_reilly +""" + +import numpy as np +import tkinter as tk +from tkinter import ttk +import os +import gradientCalculationV2 as gradCalc + +from matplotlib.backends.backend_tkagg import (FigureCanvasTkAgg) +from matplotlib.figure import Figure +from matplotlib.pyplot import close as closePlots +from mpl_toolkits.mplot3d import Axes3D +from matplotlib import cm + +class PopUpTextMessage(tk.Toplevel): + def __init__(self, parent, textmsg): + super(PopUpTextMessage, self).__init__(parent) + + self.wm_title('!') + label = ttk.Label(self, text=textmsg) + label.pack(side='top', fill='x', pady=10, padx=10) + + b1 = ttk.Button(self, text='Ok', command=self.cleanup) + b1.pack(pady=10, padx=10) + + def cleanup(self): + self.destroy() + +class GDT_GUI(tk.Tk): + + def __init__ (self, *args, **kwargs): + super(GDT_GUI, self).__init__(*args, **kwargs) + # set name and icon + # super().iconbitmap(self, default=None) + super().wm_title('Gradient design tool') + + 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) + + directionLabel = ttk.Label(inputFrame, text='Gradient direction') + 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=1,column=1, padx = 20) + + targetRegionLabel = ttk.Label(inputFrame, text='Target linear region (mm)') + 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=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=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=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) + + wireDiameterLabel = ttk.Label(inputFrame, text='Wire diameter (mm)') + 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) + + resolutionLabel = ttk.Label(inputFrame, text='Linearity order') + resolutionLabel.grid(row=8,column=0, pady=5,sticky=tk.W) + resolutionInput = ttk.Entry(inputFrame, textvariable=self.linearityOrder, width = 10, justify=tk.RIGHT) + resolutionInput.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) + + linearityOrderLabel = ttk.Label(inputFrame, text='Simulation volume (mm)') + linearityOrderLabel.grid(row=10,column=0, pady=5,sticky=tk.W) + linearityOrderInput = ttk.Entry(inputFrame, textvariable=self.simulationSize, width = 10, justify=tk.RIGHT) + linearityOrderInput.grid(row=10,column=1,padx = 20) + + simDimensionsLabel = ttk.Label(inputFrame, text='Simulation resolution (mm)') + simDimensionsLabel.grid(row=11,column=0, pady=5,sticky=tk.W) + simDimensionsInput = ttk.Entry(inputFrame, textvariable=self.simulationRes, width = 10, justify=tk.RIGHT) + simDimensionsInput.grid(row=11,column=1,padx = 20) + + DSVLabel = ttk.Label(inputFrame, text='Error DSV (mm)') + DSVLabel.grid(row=12,column=0, pady=5,sticky=tk.W) + DSVInput = ttk.Entry(inputFrame, textvariable=self.DSV, width = 10, justify=tk.RIGHT) + DSVInput.grid(row=12,column=1,padx = 20) + + calculateWireButton = ttk.Button(inputFrame, text="Calculate wire pattern", command=self.calculateGradientWires) + calculateWireButton.grid(row=13, column=0, pady=5, sticky=tk.W) + + self.calculateFieldButton = ttk.Button(inputFrame, text="Simulate B-field", command=self.calculateGradientField) + self.calculateFieldButton.grid(row=13, column=1, padx=20) + + ''' OUTPUT FRAME''' + outputFrame = ttk.Frame(self) + 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) + + 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) + + 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) + + 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 (beta):') + 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) + 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) + + '''Position plots''' + + 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.draw() + self.contourCanvas.get_tk_widget().grid(row=0,column=1) + + 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.draw() + self.fieldCanvas.get_tk_widget().grid(row=1,column=1) + + self.calculateGradientWires() + + + def calculateGradientWires(self): + + 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()) + + 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()) + + 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.gradEfficiencyString.set("%.4f mT/m/A"%(np.divide(1,self.contourData.levels[1] - self.contourData.levels[0]))) + + self.contourfig.gca().clear() + self.contourfig.gca().set_xlabel('Z (mm)') + self.contourfig.gca().set_ylabel('Y (mm)') + self.contourfig.gca().set_zlabel('X (mm)') + + 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[:,1],wirePath3D[:,2],'red') + else: + self.contourfig.gca().plot3D(wirePath3D[:,0],wirePath3D[:,1],wirePath3D[:,2],'blue') + 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.wireLengthString.set("%.2f meters"%(self.length*1e-3)) + self.contourCanvas.draw() + self.contourCanvas.flush_events() + + 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.calculateFieldButton['state'] = 'disabled' + super().update() + 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[1], self.coords[2], rstride=1, cstride=1, facecolors=cm.jet(bNorm)) + self.fieldCanvas.draw() + self.fieldCanvas.flush_events() + self.gradErrorString.set("%.2f %%"%(error*100)) + self.calculateFieldButton['state'] = 'normal' + 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 == 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 ) + bFieldOutputFile.close() + + def exportWireCSV(self): + 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()) + + for key in self.wireDict: + filename = os.path.join(folder,file+key+extension) + np.savetxt(filename,self.wireDict[key],delimiter=",", fmt='%f') + + +app = GDT_GUI() +app.mainloop() +closePlots('all') \ No newline at end of file