polymachine/torus_raytrace
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

reduce to one axis to fix width scaling

Browse source
This commit is contained in:
Daniel Mevec 2023-02-11 22:28:02 +01:00
parent 908cd306c0
commit d25458f28a
1 changed files with 128 additions and 80 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

208
torus.py
View file

@ -83,13 +83,8 @@ class TorusWorld:
class NoInamge(): class NoInamge():
def __init__(self, rfrac_init=0.5, sun_init=np.pi/2): def __init__(self, rfrac_init=0.5, sun_init=np.pi/2):
fig, (ax_side, ax_top, ax_map) = plt.subplots(3, 1, gridspec_kw={'height_ratios': [2, 2, 1]}) self.fig, self.ax = plt.subplots()
self.fig = fig
self.ax = dict(
map=ax_map,
top=ax_top,
side=ax_side,
)
self.lines = dict() self.lines = dict()
self.sun_kwargs = dict(marker='o', color='r', markersize=10,) self.sun_kwargs = dict(marker='o', color='r', markersize=10,)
@ -98,85 +93,131 @@ class NoInamge():
self.torus = TorusWorld(rfrac_init) self.torus = TorusWorld(rfrac_init)
self.illumination = None self.illumination = None
self.res = (50, 100) self.res = 200
self.rectangular_map = True
self.update_illumination() self.update_illumination()
self.ax.set_aspect('equal')
self.ax.axis('off')
self.init_top_view() self.init_top_view()
self.init_side_view() self.init_side_view()
self.init_map_view() self.init_map_view()
def update_illumination(self): def update_illumination(self):
# TODO: refactor TorusWorld.illumination() do use vectoization!! # TODO: refactor TorusWorld.illumination() do use vectoization!!
phi = np.linspace(-np.pi, np.pi, self.res[0]) phi = np.linspace(-np.pi, np.pi, self.res)
theta = np.linspace(-np.pi, np.pi, self.res[1]) theta = np.linspace(-np.pi, np.pi, self.res)
self.illumination = np.array([[self.torus.illumination(ph, th) for th in theta] for ph in phi]) self.illumination = np.array([[self.torus.illumination(ph, th) for th in theta] for ph in phi])
def init_map_view(self): def init_map_view(self):
self.lines['map_border'], = self.ax['map'].plot(*self._mantle_map(), 'k') self.lines['map_border'], = self.ax.plot(*self._mantle_map(), 'k')
self.lines['pos_map'], = self.ax['map'].plot(*self._sunpos_map(), marker='o', color='r', markersize=10,) self.lines['pos_map'], = self.ax.plot(*self._sunpos_map(), marker='o', color='r', markersize=10,)
self.lines['dawnline_map'] = [ self.lines['dawnline_map'] = [
self.ax['map'].contourf(*self._contour_map(), self.levels, **self.contour_kwargs), self.ax.contourf(*self._contour_map(), self.levels, **self.contour_kwargs),
] ]
self.ax['map'].set_aspect('equal')
self.ax['map'].set_ylabel('Longitude')
self.ax['map'].set_xlabel('Lattitude')
self.ax['map'].axis('off')
def init_top_view(self): def init_top_view(self):
self.lines['circles_top'], = self.ax['top'].plot(*self._top_section(), 'k') self.lines['circles_top'], = self.ax.plot(*self._top_section(), 'k')
self.lines['path_top'], = self.ax['top'].plot(*self._sunpath_top(), 'k:') self.lines['path_top'], = self.ax.plot(*self._sunpath_top(), 'k:')
self.lines['pos_top'], = self.ax['top'].plot(*self._sunpos_top(), **self.sun_kwargs) self.lines['pos_top'], = self.ax.plot(*self._sunpos_top(), **self.sun_kwargs)
self.ax['top'].set_xlim(-2.05, 2.05) self.lines['dawnline_top'] = [
self.ax['top'].set_ylim(-2.05, None) self.ax.contourf(*self._contour_top(), self.levels, **self.contour_kwargs),
self.ax['top'].set_aspect('equal') ]
self.ax['top'].axis('off')
def init_side_view(self): def init_side_view(self):
self.lines['circles_side'], = self.ax['side'].plot(*self._crossection(), 'k') self.lines['circles_side'], = self.ax.plot(*self._crossection(), 'k')
self.lines['path_side'], = self.ax['side'].plot(*self._sunpath_side(), 'k:') self.lines['path_side'], = self.ax.plot(*self._sunpath_side(), 'k:')
self.lines['pos_side'], = self.ax['side'].plot(*self._sunpos_side(), **self.sun_kwargs) self.lines['pos_side'], = self.ax.plot(*self._sunpos_side(), **self.sun_kwargs)
self.ax['side'].set_xlim(-2.05, 2.05) self.lines['dawnline_side'] = [
self.ax['side'].set_aspect('equal') self.ax.contourf(*self._contour_side(), self.levels, **self.contour_kwargs),
# self.ax['side'].set_ylim(-self.torus.r_min, None) ]
self.ax['side'].axis('off')
def _offset_map(self):
return -self.torus.r_maj*1.1 - (1+np.pi*self._scale_map())*self.torus.r_min
def _offset_side(self):
return self.torus.r_maj*1.1
def _scale_map(self):
if self.rectangular_map:
return (self.torus.r_maj+self.torus.r_min)/(self.torus.r_maj*np.pi)
else:
return 1/np.pi
def _mantle_map(self, n=1000): def _mantle_map(self, n=1000):
phi = np.linspace(-np.pi, np.pi, n) if self.rectangular_map:
u = phi * self.torus.r_min a = self.torus.r_maj*np.pi
width = np.pi*(self.torus.r_maj - self.torus.r_min*np.cos(phi)) b = self.torus.r_min*np.pi
data = np.array([width, u]) data = np.transpose([[a, b], [-a, b], [-a, -b], [a, -b], [a, b]])
data2 = np.array([-width, u[::-1]]) else:
phi = np.linspace(-np.pi, np.pi, n)
u = phi * self.torus.r_min
width = np.pi*(self.torus.r_maj - self.torus.r_min*np.cos(phi))
data = np.array([width, u])
data2 = np.array([-width, u[::-1]])
data = np.append(data, data2, axis=1) data = np.append(data, data2, axis=1)
data = np.append(data, data[:, 0:1], axis=1) data = np.append(data, data[:, 0:1], axis=1)
data *= self._scale_map()
data[1, :] += self._offset_map()
return data return data
def _contour_map(self): def _contour_map(self):
phi = np.linspace(-np.pi, np.pi, self.res[0]) if self.rectangular_map:
y = np.array([[w]*self.res[1] for w in self.torus.r_min*phi]) a = self.torus.r_maj*np.pi
func = (np.pi)*(self.torus.r_maj - self.torus.r_min*np.cos(phi)) b = self.torus.r_min*np.pi
x = np.array([np.linspace(-f, f, self.res[1]) for f in func]) phi = np.linspace(-b, b, self.res)
# x, y = np.meshgrid(np.linspace(-func[0], func[0], n), rmin*theta) theta = np.linspace(-a, a, self.res)
x, y = np.meshgrid(theta, phi)
else:
phi = np.linspace(-np.pi, np.pi, self.res)
width = (np.pi)*(self.torus.r_maj - self.torus.r_min*np.cos(phi))
x = np.array([np.linspace(-w, w, self.res) for w in width])
y = np.array([[row]*self.res for row in self.torus.r_min*phi])
x, y = np.meshgrid(np.linspace(-width[0], width[0], self.res), self.torus.r_min*phi)
z = self.illumination z = self.illumination
x *= self._scale_map()
y *= self._scale_map()
return x, y + self._offset_map(), z
def _contour_top(self):
n = self.res
phi = np.linspace(0, np.pi, int(np.ceil(n/2)))
theta = np.linspace(0, -np.pi, int(np.ceil(n/2)))
phx, thy = np.meshgrid(phi, theta)
x, y, _ = self.torus.surface_point(np.transpose(phx), np.transpose(thy))
z = self.illumination[int(np.floor(n/2)):, int(np.floor(n/2)):]
return x, y, z return x, y, z
def _contour_side(self):
n = self.res
phi = np.linspace(-np.pi/2, np.pi/2, int(np.ceil(n/2)))
theta = np.linspace(np.pi, 0, int(np.ceil(n/2)))
phx, thy = np.meshgrid(phi, theta)
x, _, y = self.torus.surface_point(np.transpose(phx), np.transpose(thy))
z = self.illumination[int(np.floor(n/4)):int(np.ceil(n*3/4)), :int(np.floor(n/2))]
return x, y + self._offset_side(), z
def _sunpos_map(self): def _sunpos_map(self):
phi, theta = self.torus.sun_r if self.rectangular_map:
x = (self.torus.r_maj - self.torus.r_min*np.cos(phi)) * theta y, x = self.torus.sun_r * np.array([self.torus.r_min, self.torus.r_maj])
y = self.torus.r_min * phi else:
return x, y phi, theta = self.torus.sun_r
x = (self.torus.r_maj - self.torus.r_min*np.cos(phi)) * theta
y = self.torus.r_min * phi
x *= self._scale_map()
y *= self._scale_map()
return x, y + self._offset_map()
def _top_section(self, n=1000): def _top_section(self, n=1000):
phi = np.linspace(0, np.pi, n) theta = np.linspace(0, np.pi, n)
x1 = (self.torus.r_maj + self.torus.r_min) * np.cos(phi) x1 = (self.torus.r_maj + self.torus.r_min) * np.cos(theta)
y1 = -1 * (self.torus.r_maj + self.torus.r_min) * np.sin(phi) y1 = -1 * (self.torus.r_maj + self.torus.r_min) * np.sin(theta)
data1 = np.array([x1, y1]) data1 = np.array([x1, y1])
x2 = (self.torus.r_maj - self.torus.r_min) * np.cos(phi) x2 = (self.torus.r_maj - self.torus.r_min) * np.cos(theta)
y2 = -1 * (self.torus.r_maj - self.torus.r_min) * np.sin(phi) y2 = -1 * (self.torus.r_maj - self.torus.r_min) * np.sin(theta)
data2 = np.array([x2, y2]) data2 = np.array([x2, y2])
data = np.append(data1, [[np.nan,], [np.nan,]], axis=1) data = np.append(data1, [[np.nan,], [np.nan,]], axis=1)
@ -191,7 +232,7 @@ class NoInamge():
def _sunpos_top(self): def _sunpos_top(self):
phi, theta = self.torus.sun_r phi, theta = self.torus.sun_r
x = self.torus.r_maj + self.torus.r_maj*np.cos(phi) x = self.torus.r_maj - self.torus.r_maj*np.cos(phi)
y = 0 y = 0
return x, y return x, y
@ -203,56 +244,63 @@ class NoInamge():
data2 = np.array([-x, y]) data2 = np.array([-x, y])
data = np.append(data, [[np.nan,], [np.nan,]], axis=1) data = np.append(data, [[np.nan,], [np.nan,]], axis=1)
data = np.append(data, data2, axis=1) data = np.append(data, data2, axis=1)
data[1, :] += self._offset_side()
return data return data
def _sunpath_side(self, n=1000): def _sunpath_side(self, n=1000):
phi = np.linspace(0, 2*np.pi, n) phi = np.linspace(0, 2*np.pi, n)
x = self.torus.r_maj + self.torus.r_maj*np.cos(phi) x = self.torus.r_maj + self.torus.r_maj*np.cos(phi)
y = self.torus.r_maj*np.sin(phi) y = self.torus.r_maj*np.sin(phi)
return np.array([x, y]) return np.array([x, y+self._offset_side()])
def _sunpos_side(self): def _sunpos_side(self):
phi, theta = self.torus.sun_r phi, theta = self.torus.sun_r
x = self.torus.r_maj - self.torus.r_maj*np.cos(phi) x = self.torus.r_maj - self.torus.r_maj*np.cos(phi)
y = self.torus.r_maj*np.sin(phi) y = self.torus.r_maj*np.sin(phi)
return x, y return x, y+self._offset_side()
@staticmethod def redraw_plot(self, line, func):
def redraw_plot(ax, line, func):
x, y = func() x, y = func()
line.set_xdata(x) line.set_xdata(x)
line.set_ydata(y) line.set_ydata(y)
ax.relim() self.ax.relim()
ax.autoscale_view() self.ax.autoscale_view()
@staticmethod def redraw_contourf(self, container, func, levels=None, contour_kwargs=None):
def redraw_contourf(ax, container, func, levels=None, contour_kwargs=None):
levels = levels if levels else [-1, 0, 1] levels = levels if levels else [-1, 0, 1]
contour_kwargs = contour_kwargs if contour_kwargs else {} contour_kwargs = contour_kwargs if contour_kwargs else {}
for coll in container[0].collections: for coll in container[0].collections:
coll.remove() coll.remove()
container[0] = ax.contourf(*func(), levels, **contour_kwargs) container[0] = self.ax.contourf(*func(), levels, **contour_kwargs)
def update_torus(self, rfrac): def update_torus(self, rfrac):
self.torus.update(rfrac) self.torus.update(rfrac)
self.update_illumination() self.update_illumination()
self.redraw_plot(self.ax['map'], self.lines['map_border'], self._mantle_map) self.redraw_plot(self.lines['map_border'], self._mantle_map)
self.redraw_plot(self.ax['side'], self.lines['circles_side'], self._crossection) self.redraw_plot(self.lines['circles_side'], self._crossection)
self.redraw_plot(self.ax['side'], self.lines['path_side'], self._sunpath_side) self.redraw_plot(self.lines['path_side'], self._sunpath_side)
self.redraw_plot(self.ax['top'], self.lines['circles_top'], self._top_section) self.redraw_plot(self.lines['circles_top'], self._top_section)
self.redraw_plot(self.ax['top'], self.lines['path_top'], self._sunpath_top) self.redraw_plot(self.lines['path_top'], self._sunpath_top)
self.redraw_contourf(self.ax['map'], self.lines['dawnline_map'], self.redraw_contourf(self.lines['dawnline_map'],
self._contour_map, self.levels, self.contour_kwargs) self._contour_map, self.levels, self.contour_kwargs)
self.redraw_contourf(self.lines['dawnline_top'],
self._contour_top, self.levels, self.contour_kwargs)
self.redraw_contourf(self.lines['dawnline_side'],
self._contour_side, self.levels, self.contour_kwargs)
self.fig.canvas.draw_idle() self.fig.canvas.draw_idle()
def update_sun(self, phi, theta): def update_sun(self, phi, theta):
self.torus.put_sun(phi, theta) self.torus.put_sun(phi, theta)
self.update_illumination() self.update_illumination()
self.redraw_plot(self.ax['map'], self.lines['pos_map'], self._sunpos_map) self.redraw_plot(self.lines['pos_map'], self._sunpos_map)
self.redraw_plot(self.ax['side'], self.lines['pos_side'], self._sunpos_side) self.redraw_plot(self.lines['pos_side'], self._sunpos_side)
self.redraw_plot(self.ax['top'], self.lines['pos_top'], self._sunpos_top) self.redraw_plot(self.lines['pos_top'], self._sunpos_top)
self.redraw_contourf(self.ax['map'], self.lines['dawnline_map'], self.redraw_contourf(self.lines['dawnline_map'],
self._contour_map, self.levels, self.contour_kwargs) self._contour_map, self.levels, self.contour_kwargs)
self.redraw_contourf(self.lines['dawnline_top'],
self._contour_top, self.levels, self.contour_kwargs)
self.redraw_contourf(self.lines['dawnline_side'],
self._contour_side, self.levels, self.contour_kwargs)
self.fig.canvas.draw_idle() self.fig.canvas.draw_idle()
@ -272,19 +320,19 @@ class ImageInteractive(NoInamge):
def init_interactivity(self, rfrac_init, sun_init): def init_interactivity(self, rfrac_init, sun_init):
self.fig.subplots_adjust(left=0.25, bottom=0.25) self.fig.subplots_adjust(left=0.25, bottom=0.25)
self.ax['slider_sun'] = self.fig.add_axes([0.25, 0.1, 0.65, 0.03]) ax1 = self.fig.add_axes([0.25, 0.1, 0.65, 0.03])
self.ax['slider_rf'] = self.fig.add_axes([0.1, 0.25, 0.0225, 0.63]) ax2 = self.fig.add_axes([0.1, 0.25, 0.0225, 0.63])
self.ax['button_reset'] = self.fig.add_axes([0.8, 0.025, 0.1, 0.04]) ax3 = self.fig.add_axes([0.8, 0.025, 0.1, 0.04])
self.sliders = dict( self.sliders = dict(
sun_phi=Slider( sun_phi=Slider(
ax=self.ax['slider_sun'], ax=ax1,
label='Angle of Sun', label='Angle of Sun',
valmin=-np.pi, valmin=-np.pi,
valmax=np.pi, valmax=np.pi,
valinit=sun_init, valinit=sun_init,
), ),
rfrac=Slider( rfrac=Slider(
ax=self.ax['slider_rf'], ax=ax2,
label="Fraction of Radii (r/R)", label="Fraction of Radii (r/R)",
valmin=0, valmin=0,
valmax=1, valmax=1,
@ -295,7 +343,7 @@ class ImageInteractive(NoInamge):
) )
self.sliders['sun_phi'].on_changed(self._slider_update_sun) self.sliders['sun_phi'].on_changed(self._slider_update_sun)
self.sliders['rfrac'].on_changed(self._slider_update_torus) self.sliders['rfrac'].on_changed(self._slider_update_torus)
button = Button(self.ax['button_reset'], 'Reset', hovercolor='0.975') button = Button(ax3, 'Reset', hovercolor='0.975')
button.on_clicked(self._reset) button.on_clicked(self._reset)
def _slider_update_torus(self, val): def _slider_update_torus(self, val):