From cb69732f68c0bd46c1574de16ce1aee6f38e439b Mon Sep 17 00:00:00 2001 From: Brett Weiland Date: Tue, 11 Jun 2024 14:50:14 -0500 Subject: restarting --- field_tests/field.py | 191 +++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 191 insertions(+) create mode 100755 field_tests/field.py (limited to 'field_tests/field.py') diff --git a/field_tests/field.py b/field_tests/field.py new file mode 100755 index 0000000..af2abbf --- /dev/null +++ b/field_tests/field.py @@ -0,0 +1,191 @@ +#!/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) + -- cgit v1.2.3