192 lines
5.4 KiB
Python
Executable File
192 lines
5.4 KiB
Python
Executable File
#!/usr/bin/env python3
|
|
import numpy as np
|
|
import matplotlib.pyplot as plt
|
|
from alive_progress import alive_bar
|
|
|
|
img_res_x = 250
|
|
img_res_y = 250
|
|
total_pixels = img_res_x * img_res_y # so we don't gotta compute it every time
|
|
|
|
periods = 1
|
|
square_x = 0
|
|
square_y = 0
|
|
|
|
xmin = (-periods * np.pi) + (square_x * np.pi)
|
|
xmax = (periods * np.pi) + (square_x * np.pi)
|
|
ymin = (-periods * np.pi) + (square_y * np.pi)
|
|
ymax = (periods * np.pi) + (square_y * np.pi)
|
|
|
|
#xmin = -10
|
|
#xmax = 10
|
|
#ymin = -10
|
|
#ymax = 10
|
|
|
|
escape = 10000
|
|
iterations = 255*3
|
|
c_x = 2 * np.pi
|
|
c_y = 2 * np.pi
|
|
|
|
|
|
grid = np.meshgrid(np.linspace(ymin, ymax, 200), np.linspace(xmin, xmax, 200))
|
|
#image = np.meshgrid(np.linspace(ymin, ymax, img_res_y), np.linspace(xmin, xmax, img_res_x))
|
|
image = np.zeros([img_res_y, img_res_x], dtype=np.double)
|
|
|
|
|
|
class point_charge():
|
|
def __init__(self, x, y, c, mod):
|
|
self.x = x
|
|
self.y = y
|
|
self.c = c
|
|
self.mod = mod
|
|
def get_field(self, to_x, to_y):
|
|
if(self.mod):
|
|
to_x = (to_x % self.mod)
|
|
to_y = (to_y % self.mod)
|
|
return np.array([
|
|
((self.c * (self.x - to_x)) / ((self.x - to_x)**2 + (self.y - to_y)**2)**1.5),
|
|
((self.c * (self.y - to_y)) / ((self.x - to_x)**2 + (self.y - to_y)**2)**1.5)])
|
|
|
|
#will remove all the point charge code if it turns out to be good enough to be impliemnted into openCL
|
|
#point_charges = [point_charge(-5, -5, 100), point_charge(-5, 5, -100), point_charge(5, 0, 100)]
|
|
#point_charges = [point_charge(-1,-1, 100, 10), point_charge(1,1,-100, 0)]
|
|
point_charges = []
|
|
|
|
|
|
#plt.ion()
|
|
ax = plt.gca()
|
|
fig = plt.gcf()
|
|
#ax.set_autoscale_on(False)
|
|
#ax.set_xlim([xmin, xmax])
|
|
#ax.set_ylim([ymin, ymax])
|
|
|
|
vector_arrows = None
|
|
|
|
def show_field():
|
|
global vector_arrows
|
|
grid_f = np.zeros_like(grid)
|
|
for p in point_charges:
|
|
grid_f += p.get_field(grid[0], grid[1])
|
|
#plt.streamplot(grid[0], grid[1], grid_f[0], grid_f[1], density=5)
|
|
vector_arrows = plt.quiver(grid[0], grid[1], grid_f[0], grid_f[1])
|
|
plt.show(block=False)
|
|
plt.pause(.1)
|
|
|
|
|
|
#show_field()
|
|
|
|
timestep = .1
|
|
def test_sim():
|
|
particle_grid = np.meshgrid(np.linspace(ymin, ymax, 100), np.linspace(xmin, xmax, 100))
|
|
pos = particle_grid
|
|
acceleration = np.zeros_like(particle_grid)
|
|
velocity = np.zeros_like(particle_grid)
|
|
velocity = [np.ones_like(particle_grid[0]) * 1, np.ones_like(particle_grid[0]) * .5]
|
|
mass = 10
|
|
charge = 1
|
|
particle_plot = ax.plot(velocity[0], velocity[1], 'bo', animated=True)
|
|
#velocity += .1
|
|
|
|
background = fig.canvas.copy_from_bbox(ax.bbox)
|
|
ax.draw_artist(vector_arrows)
|
|
fig.canvas.blit(fig.bbox)
|
|
|
|
while True:
|
|
fig.canvas.restore_region(background)
|
|
field = np.zeros_like(particle_grid)
|
|
# TODO can make this quicker by skipping initilization
|
|
for p in point_charges:
|
|
field += p.get_field(pos[0], pos[1])
|
|
acceleration = ((charge * field) / mass) * timestep
|
|
#print(acceleration)
|
|
velocity += acceleration * timestep
|
|
pos += velocity * timestep
|
|
|
|
fig.canvas.restore_region(background)
|
|
particle_plot[0].set_data(pos[0],pos[1])
|
|
ax.draw_artist(particle_plot[0])
|
|
fig.canvas.blit(fig.bbox)
|
|
fig.canvas.flush_events()
|
|
plt.pause(1/60)
|
|
|
|
#fig.canvas.draw_idle()
|
|
#test_sim()
|
|
|
|
#exit(1)
|
|
|
|
|
|
max_timesteps = 10
|
|
|
|
def get_fractal_iter(img):
|
|
next_x = img[1][y][x] / np.sin(img[0][y][x])
|
|
img[0][y][x] = img[0][y][x] / np.sin(img[1][y][x])
|
|
img[1][y][x] = next_x
|
|
if (np.square(img[0][y][x]) + np.square(img[1][y][x])) >= escape:
|
|
z[y][x] = i
|
|
|
|
|
|
#meshgrid makes things slower as we can't test individual points for breaking to infinity;
|
|
#however, I will fix that later.
|
|
cmap = plt.cm.viridis
|
|
cmap.set_bad((0,0,0))
|
|
cmap.set_over((0,0,0))
|
|
cmap.set_under((0,0,0))
|
|
|
|
with alive_bar(img_res_y, bar = 'filling', spinner = 'waves') as bar:
|
|
for pix_y, y in enumerate(np.linspace(ymin, ymax, img_res_y)):
|
|
for pix_x, x in enumerate(np.linspace(xmin, xmax, img_res_x)):
|
|
on_x = x
|
|
on_y = y
|
|
for i in range(iterations):
|
|
completed_ratio = (((pix_x * pix_y * 1)) / total_pixels)
|
|
next_x = on_x/np.sin(on_y)
|
|
on_y = on_y/np.sin(on_x)
|
|
on_x = next_x
|
|
|
|
|
|
# do physics here - we could just use vectors but i keep rewriting things
|
|
timesteps = max_
|
|
acceleration = np.array([0, 0], dtype=np.double)
|
|
velocity = np.array([on_x, on_y], dtype=np.double) # maybe multiply by stuff
|
|
pos = np.array([0,0], dtype=np.double)
|
|
for t in range(timesteps):
|
|
for p in point_charges:
|
|
acceleration += p.get_field(on_x, on_y)
|
|
velocity += acceleration
|
|
pos += velocity * (timesteps)
|
|
on_x += pos[0]
|
|
on_y += pos[1]
|
|
|
|
if on_x**2 + on_y**2 > escape:
|
|
image[pix_x][pix_y] = i
|
|
break
|
|
bar()
|
|
|
|
ax.imshow(image, norm="log", aspect="auto", cmap=cmap)
|
|
plt.show()
|
|
|
|
|
|
# yeah, I shouldn't have switched to a meshgrid, oh well
|
|
#z = np.empty_like(image[0])
|
|
exit(0)
|
|
with alive_bar(img_res_x, bar = 'filling', spinner = 'waves') as bar:
|
|
for y in range(img_res_y):
|
|
for x in range(img_res_x):
|
|
for i in range(iterations):
|
|
if image[0][y][x] == np.NAN:
|
|
continue
|
|
next_x = image[1][y][x] / np.sin(image[0][y][x])
|
|
image[0][y][x] = image[0][y][x] / np.sin(image[1][y][x])
|
|
image[1][y][x] = next_x
|
|
if (np.square(image[0][y][x]) + np.square(image[1][y][x])) >= escape:
|
|
z[y][x] = i
|
|
image[0][y][x] = np.NAN
|
|
image[1][y][x] = np.NAN
|
|
break
|
|
# #do physics here
|
|
|
|
|
|
|
|
bar()
|
|
#image = np.clip(image, -escape, escape)
|
|
|