Commit 972569ea authored by Otthorn's avatar Otthorn 💎

PEP8 and fixed OOB syntax and logic

parent c9014be6
......@@ -5,6 +5,9 @@
/dist/
/data/
# IDE file
/.idea
# Python egg metadata, regenerated from source files by setuptools.
/*.egg-info
......
......@@ -44,7 +44,7 @@ The observation grid is the collection of points where the B field will be evalu
Currently, there is a polygon class that can make regular polygons to approximate your coil shape. If you need a non-regular polygon shape, then create a class that inherits from the loop class and calculates the vertices of the non-regular shape. The new shape class then needs to call create_segments and pass in the vertices. Consider pushing your shape to the repo so others may use it too! Copy and paste the code block below into your shape class to do this.
```python
self.segments = self.create_segments(self.vertices)
self.segments = self.create_segments()
```
1. Look in the section called: *create coil*
......
#!/usr/bin/env python
from pylab import *
import matplotlib.pyplot as plt
import scipy.integrate as integrate
import staticbs as bs
import sys
import numpy as np
from plot import *
def save(data_dir='./data/', **kwargs):
......@@ -11,23 +11,28 @@ def save(data_dir='./data/', **kwargs):
print((fn, kwargs[fn]))
np.save(data_dir + fn, kwargs[fn])
# print dir(static_bs_module)
# observation grid setup
# grid units are in meters
x = linspace(-0.2, 0.2, 20)
y = linspace(-0.2, 0.2, 20)
x = linspace(-0.2, 0.2, 100)
y = linspace(-0.2, 0.2, 100)
enable_3D = False
#====2D case====
## Example for how to ignore the z axis
z = zeros((x.shape[0])) + 1e-2
zz = zeros(x.shape + y.shape) + 1e-2
xx, yy = meshgrid(x, y)
# ====2D case====
if not enable_3D:
# Example for how to ignore the z axis
z = zeros((x.shape[0])) + 1e-2
zz = zeros(x.shape + y.shape) + 1e-2
xx, yy = meshgrid(x, y)
#====3D case====
#z = linspace(-0.2, 0.2, 20)
#xx, yy, zz = meshgrid(x, y, z)
# ====3D case====
else:
z = linspace(-0.2, 0.2, 20)
xx, yy, zz = meshgrid(x, y, z)
observation_points = rollaxis(np.array((xx, yy, zz)), 0, xx.ndim+1)
......@@ -42,14 +47,17 @@ amp_turns = simulation_current * num_loops
velocity = np.array([0.001 ,0, 0]) # m/s
# create coil
coil_radius = 0.1 # meters
coil = bs.shapes.polygon(3, (0,0,0), coil_radius)
coil_radius = 0.05 # meters
coil_A = bs.shapes.Polygon(4, (-0.1,0,0), coil_radius)
coil_B = bs.shapes.Polygon(4, (0.1,0,0), coil_radius)
# make a list of all of the current segments
current_object = bs.iobjs.current_obj_list(coil.get_segments() )
# current_A = bs.iobjs.current_obj_list(coil_A.get_segments())
current_object = bs.iobjs.CurrentObjList(coil_A.get_segments())
#make e field object
e_field = bs.eobjs.e_field(amp_turns, current_object, velocity)
# make e field object
e_field = bs.eobjs.EField(amp_turns, current_object, velocity)
# calculate magnetic field at observation points
print("calculating B field")
......@@ -60,17 +68,23 @@ E_obs = e_field.E(observation_points)
# E field line integral
#result = integrate.quad(sim_obj.E_integrand, 1, 1, args=(0, 1))
#print(result)
# result = integrate.quad(sim_obj.E_integrand, 1, 1, args=(0, 1))
# print(result)
save(B_test=B_obs, E_test=E_obs, obs_points=observation_points)
from plot import *
if not enable_3D:
plot_B('./data/B_test.npy')
plot_E('./data/E_test.npy')
plot_B('./data/B_test.npy')
plot_E('./data/E_test.npy')
if enable_3D:
plot_3D('./data/B_test.npy')
import matplotlib.pyplot as plt
octagonal_coil = bs.shapes.Octagon(0.8, 0.2, 0.2, 0.2, 0.1)
# plot_coil(coil_A)
# plot_coil(coil_B)
plot_coil(octagonal_coil)
plt.show()
#!/usr/bin/env python
from pylab import *
from mpl_toolkits.mplot3d import axes3d
import numpy as np
from staticbs.common import mag
def plot_linearity(B_fn, obs_points_fn="./data/obs_points_linearity.npy"):
obs_points = np.load(obs_points_fn)
B = np.load(B_fn)
......@@ -64,6 +63,7 @@ def plot_3D(B_fn, obs_points_fn="./data/obs_points.npy"):
B[::4, ::4, 1]/float(np.amax(B[::4, ::4, 1])),
B[::4, ::4, 2]/float(np.amax(B[::4, ::4, 2])), color='blue', length=0.001)
def plot_B(B_fn, obs_points_fn="./data/obs_points.npy"):
"""Plot of 2 vector graph of the B field"""
obs_points = np.load(obs_points_fn)
......@@ -111,7 +111,7 @@ def plot_E(E_fn, obs_points_fn="./data/obs_points.npy"):
def plot_coil(coil):
"""Plot the verticies of a coil in a 3D graph"""
"""Plot the vertices of a coil in a 3D graph"""
vertices = coil.vertices
xx = vertices[:, 0]
......@@ -125,11 +125,10 @@ def plot_coil(coil):
if __name__ == '__main__':
#plot_B(sys.argv[1])
#plot_E(sys.argv[2])
#plot_3D(sys.argv[1])
# plot_B(sys.argv[1])
# plot_E(sys.argv[2])
# plot_3D(sys.argv[1])
scatter_3D(sys.argv[1])
#plot_linearity(sys.argv[1])
# plot_linearity(sys.argv[1])
plt.show()
plt.show()
\ No newline at end of file
import numpy as np
def vectorize_over_points(func):
def vectorized_func(a, *args, **kwargs):
r = np.array(tuple(func(p, *args, **kwargs) for p in a.reshape(-1, 3)))
return r.reshape(a.shape[:-1] + r.shape[1:])
return vectorized_func
def mag(vector):
"""Retun the list of the norms of all the vectors in a
"""Return the list of the norms of all the vectors in a
given array"""
return np.sqrt(np.sum(np.square(vector), axis=-1))
from .common import *
class e_field():
class EField:
def __init__(self, amp_turns, current_object, velocity):
"""calculates the e field given a current object"""
self.vectorized_mag = vectorize_over_points(mag)
......@@ -13,10 +14,12 @@ class e_field():
return np.cross(self.velocity, self.current_object.B(point, self.amp_turns))
def E_integrand(self, t, a, b):
return np.dot( self.E_field_at_point(self.r(t, a, b)), self.r_prime(t, a, b) )
return np.dot(self.E_field_at_point(self.r(t, a, b)), self.r_prime(a, b))
def r(self, t, a, b):
@staticmethod
def r(t, a, b):
return np.array([0, a + t*(b-a), 0])
def r_prime(self, t, a, b):
@staticmethod
def r_prime(a, b):
return np.array([0, (b-a), 0])
......@@ -2,36 +2,38 @@
from .common import *
import numpy as np
import scipy.constants as consts
import pprint
import copy
class current_obj_list(list):
class CurrentObjList(list):
def B(self, points, current):
return np.sum((i.B(points, current) for i in self), axis=0)
class loop:
class Loop:
"""creates straight wire segments from points. It closes the loop if it is left open.
Current flows in the order the points are given."""
def __init__(self, I, points):
def __init__(self, points):
self.vertices = points
self.segments = self.create_segments(self.vertices)
self.segments = self.create_segments()
def create_segments(self, points):
def create_segments(self):
points = self.vertices
wire_list = []
for i in range(0, len(points)):
if (i+1) == len(points):
wire_list.append(straight_wire(points[i], points[0]))
wire_list.append(StraightWire(points[i], points[0]))
else:
wire_list.append(straight_wire(points[i], points[i+1]))
wire_list.append(StraightWire(points[i], points[i+1]))
return copy.deepcopy(wire_list)
def get_segments(self):
return copy.deepcopy(self.segments)
class straight_wire:
class StraightWire:
def __init__(self, start, end):
self.start = start # xyz tuple
......@@ -39,11 +41,11 @@ class straight_wire:
self.unit_vector = (end-start)/mag(end-start)
self.B = vectorize_over_points(self.B_at_point)
def perpendicular_vector(self, observation_point, start, end):
def perpendicular_vector(self, observation_point):
# source point is origin
bold_r = observation_point - start
wire_vector = (start - end)
wv_unit_v = wire_vector / mag( wire_vector )
bold_r = observation_point - self.start
wire_vector = (self.start - self.end)
wv_unit_v = wire_vector / mag(wire_vector)
# in the direction of the unit vector
r_prime = np.dot(bold_r, wv_unit_v)
......@@ -69,13 +71,13 @@ class straight_wire:
start = self.start
end = self.end
r = self.perpendicular_vector(observation_point, start, end)
r = self.perpendicular_vector(observation_point)
r_mag = mag(r)
unit_vector = (end - start) / mag( end - start)
l1 = np.dot( (observation_point - start), unit_vector )
l2 = np.dot( (observation_point- end), -1*unit_vector )
unit_vector = (end - start) / mag(end - start)
l1 = np.dot(observation_point - start, unit_vector)
l2 = np.dot(observation_point - end, -1*unit_vector)
cos_theta_1 = l1/np.sqrt(r_mag**2 + np.square(l1))
cos_theta_2 = l2/np.sqrt(r_mag**2 + np.square(l2))
......@@ -83,9 +85,9 @@ class straight_wire:
B = (consts.mu_0/(4*consts.pi)) * \
(current/r_mag)*(cos_theta_2 + cos_theta_1)
return np.sqrt(B**2)
def get_wire_points(self, num_samples):
wire_points = make_3d_line(self.start, self.end, num_samples)
# where is the function make_3d_lines ?
return wire_points
import numpy as np
# useless import since cart2pol and pol2cart are not used
from .iobjs import *
def cart2pol(x, y):
"""Convert carthesian coordonates to polar coordonates"""
"""Convert cartesian coordinates to polar coordinates"""
rho = np.sqrt(x**2 + y**2)
phi = np.arctan2(y, x)
return(rho, phi)
return rho, phi
def pol2cart(rho, phi):
"""Convert polar coordonates to carthesian coordonates"""
"""Convert polar coordinates to cartesian coordinates"""
x = rho * np.cos(phi)
y = rho * np.sin(phi)
return(x, y)
return x, y
class polygon(loop):
class Polygon(Loop):
"""Class to define the geometry of an n-sized polygon"""
def __init__(self, n_sides, center, radius):
......@@ -28,29 +32,29 @@ class polygon(loop):
theta = self.phi
self.rot_matrix_arr = [np.matrix([[np.cos(theta), -np.sin(theta), 0],\
[np.sin(theta), np.cos(theta), 0],
[0, 0, 1]]) for theta in -np.arange(0, 2*np.pi, theta)]
self.rot_matrix_arr = [np.matrix([[np.cos(theta), -np.sin(theta), 0],
[np.sin(theta), np.cos(theta), 0],
[0, 0, 1]]) for theta in -np.arange(0, 2*np.pi, theta)]
self.vertices = np.array(self.rot_matrix_arr * self.prime_vector.T)
self.segments = self.create_segments(self.vertices)
self.segments = self.create_segments()
def area(self, vertices):
x = vertices[:, 0]
y = vertices[:, 1]
def area(self):
x = self.vertices[:, 0]
y = self.vertices[:, 1]
area = 0.0
num_points = len(x)
j = num_points - 1
for i in range(0, num_points):
area = area + (x[j] + x[i]) * (y[j] - y[i])
area = area + (x[j] + x[i]) * (y[j] - y[i])
j = i
return area / 2.0
class octagon(loop):
class Octagon(Loop):
def __init__(self, h_straight_length_percent, v_straight_length_percent, vertical_axis, horizontal_axis, height):
# height is the total height of the stack at the point where you start depositing the current layer
# thickness later on is the height of the layer being currently deposited
......@@ -64,34 +68,25 @@ class octagon(loop):
self.h_straight_length = self.horizontal_axis * self.h_straight_length_percent
self.v_straight_length = self.vertical_axis * self.v_straight_length_percent
# print "h_straight_length: ", self.h_straight_length
# print "v_straight_length: ", self.v_straight_length
self.vertices = self.calculate_vertices()
# print self.vertices
# print self.area(self.vertices)
#
# fig, ax = plt.subplots()
# plt.scatter(self.vertices[:,0], self.vertices[:,1])
# circle = plt.Circle((0,0), 0.1,fill=False)
# ax.add_artist(circle)
# plt.show()
self.segments = self.create_segments(self.vertices)
self.segments = self.create_segments()
def area(self, vertices):
x = vertices[:, 0]
y = vertices[:, 1]
@staticmethod
def area(self):
x = self.vertices[:, 0]
y = self.vertices[:, 1]
area = 0.0
num_points = len(x)
j = num_points - 1
for i in range(0, num_points):
area = area + (x[j] + x[i]) * (y[j] - y[i])
area = area + (x[j] + x[i]) * (y[j] - y[i])
j = i
return area / 2.0
def calculate_vertices(self):
center = (0, 0, self.height)
......@@ -109,5 +104,4 @@ class octagon(loop):
v8 = (x + self.h_straight_length / 2.0, y - self.vertical_axis / 2.0, z)
vertices = np.array([v1, v2, v3, v4, v5, v6, v7, v8])
# print "vertices: ", vertices
return vertices
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment