diff --git a/docs/apps/FourierTransforms/app.py b/docs/apps/FourierTransforms/app.py index 2e24842..05e88b1 100644 --- a/docs/apps/FourierTransforms/app.py +++ b/docs/apps/FourierTransforms/app.py @@ -71,8 +71,6 @@ def server(input, output, session): click_data = reactive.value(None) - blue = 'navy' - red = 'firebrick' x_axis = np.linspace(-20, 20, 1024) freq_axis = np.fft.fftfreq(len(x_axis), d=(x_axis[1] - x_axis[0])) diff --git a/docs/apps/LagrangianPoints/app.md b/docs/apps/LagrangianPoints/app.md new file mode 100644 index 0000000..c7f2569 --- /dev/null +++ b/docs/apps/LagrangianPoints/app.md @@ -0,0 +1,16 @@ +--- +authors: + - ptuemmler +categories: + - Physics +tags: + - Draft +date: 2025-11-23 +hide: + - toc +--- + +# Lagrangian Points + +{{embed_app("100%", "830px")}} +Underlying potential taken from: [https://doi.org/10.1086/381315](https://doi.org/10.1086/381315) diff --git a/docs/apps/LagrangianPoints/app.py b/docs/apps/LagrangianPoints/app.py new file mode 100644 index 0000000..dcdd527 --- /dev/null +++ b/docs/apps/LagrangianPoints/app.py @@ -0,0 +1,68 @@ +import numpy as np +import plotly.graph_objects as go +from shiny import App, Inputs, Outputs, Session, render, ui + +app_ui = ui.page_sidebar( + ui.sidebar( + ui.input_slider("q", "Mass Ratio (q)", min=0.0, max=1.0, value=0.05, step=0.01), + # ui.input_slider("x", "x", min=0.0, max=2.0, value=1.0, step=0.01), + ui.input_dark_mode(id='dark_mode'), + ), + ui.output_ui("plot") +) + +def server(input: Inputs, output: Outputs, session: Session): + @render.ui + def plot(): + # Set template based on dark mode + if input.dark_mode() == "dark": + template = "plotly_dark" + else: + template = "plotly_white" + + fig = go.Figure() + # Add isosurface representing the potential field + x, y = np.meshgrid(np.linspace(-2, 2, 256), + np.linspace(-2, 2, 256)) + z = np.zeros_like(x) + + q = input.q() + # x_pos = input.x() + values = (x - q/(1+q))**2 + y**2 + 2/((1+q)*np.sqrt(x**2 + y**2 + z**2)) + 2*q/((1+q)*np.sqrt((x-1)**2 + y**2 + z**2)) + + + values -= np.min(values) + min_value = 1.5 + values[values > min_value] = np.nan # Mask values above a certain threshold + + fig.add_trace(go.Surface( + contours = { + "z": {"show": True, "start": -min_value, "end": 0.01, "size": 0.05} + }, + x=x, + y=y, + z=-values, # Offset for better visibility + colorscale='turbo', + colorbar_title_text='f(x,y)', + showscale=True, + hoverinfo='skip' + )) + + # Set layout + fig.update_layout( + template=template, + scene=dict( + xaxis_title='X-axis', + yaxis_title='Y-axis', + zaxis_title='Z-axis', + xaxis_showspikes=False, + yaxis_showspikes=False, + zaxis_showspikes=False, + ), + height=700, + margin=dict(l=0, r=0, t=0, b=0), + ) + + return ui.HTML(fig.to_html()) + +app = App(app_ui, server) diff --git a/docs/apps/LagrangianPoints/requirements.txt b/docs/apps/LagrangianPoints/requirements.txt new file mode 100644 index 0000000..51e83f6 --- /dev/null +++ b/docs/apps/LagrangianPoints/requirements.txt @@ -0,0 +1,2 @@ +numpy +plotly \ No newline at end of file diff --git a/docs/apps/MinkowskiSpaceTime/app.md b/docs/apps/MinkowskiSpaceTime/app.md index 5041331..63fe42e 100644 --- a/docs/apps/MinkowskiSpaceTime/app.md +++ b/docs/apps/MinkowskiSpaceTime/app.md @@ -16,4 +16,6 @@ draft: true --- # Minkowski Space-Time -{{embed_app("100%", "800px")}} +{{embed_app("100%", "800px", "stationary")}} +{{embed_app("100%", "800px", "doppler")}} +{{embed_app("100%", "800px", "motion")}} diff --git a/docs/apps/MinkowskiSpaceTime/doppler/app.py b/docs/apps/MinkowskiSpaceTime/doppler/app.py new file mode 100644 index 0000000..c30c753 --- /dev/null +++ b/docs/apps/MinkowskiSpaceTime/doppler/app.py @@ -0,0 +1,151 @@ +import matplotlib.pyplot as plt +from matplotlib.colors import LinearSegmentedColormap, Normalize +import numpy as np +import pandas as pd +from shiny import App, render, ui, reactive + +berlin_cmap = [[0.621082, 0.690182, 0.999507],[0.612157, 0.689228, 0.995374],[0.603202, 0.688250, 0.991239],[0.594200, 0.687257, 0.987092],[0.585165, 0.686248, 0.982922],[0.576088, 0.685222, 0.978733],[0.566961, 0.684166, 0.974524],[0.557791, 0.683098, 0.970288],[0.548590, 0.681992, 0.966016],[0.539327, 0.680859, 0.961704],[0.530034, 0.679691, 0.957350],[0.520687, 0.678484, 0.952942],[0.511295, 0.677230, 0.948466],[0.501863, 0.675908, 0.943923],[0.492368, 0.674526, 0.939297],[0.482832, 0.673075, 0.934574],[0.473239, 0.671530, 0.929751],[0.463610, 0.669898, 0.924806],[0.453931, 0.668152, 0.919735],[0.444213, 0.666275, 0.914518],[0.434440, 0.664271, 0.909136],[0.424645, 0.662120, 0.903586],[0.414818, 0.659791, 0.897845],[0.404975, 0.657289, 0.891905],[0.395137, 0.654579, 0.885750],[0.385296, 0.651674, 0.879368],[0.375493, 0.648536, 0.872757],[0.365742, 0.645164, 0.865903],[0.356059, 0.641552, 0.858801],[0.346453, 0.637692, 0.851451],[0.336982, 0.633574, 0.843855],[0.327642, 0.629189, 0.836017],[0.318487, 0.624551, 0.827937],[0.309539, 0.619657, 0.819628],[0.300784, 0.614497, 0.811108],[0.292309, 0.609115, 0.802379],[0.284098, 0.603485, 0.793470],[0.276205, 0.597634, 0.784386],[0.268595, 0.591580, 0.775143],[0.261308, 0.585335, 0.765780],[0.254368, 0.578908, 0.756296],[0.247753, 0.572328, 0.746719],[0.241464, 0.565596, 0.737066],[0.235515, 0.558748, 0.727351],[0.229842, 0.551802, 0.717600],[0.224503, 0.544750, 0.707805],[0.219485, 0.537628, 0.697998],[0.214694, 0.530433, 0.688190],[0.210172, 0.523193, 0.678377],[0.205889, 0.515897, 0.668578],[0.201771, 0.508598, 0.658787],[0.197878, 0.501258, 0.649030],[0.194172, 0.493903, 0.639287],[0.190556, 0.486541, 0.629572],[0.187112, 0.479181, 0.619898],[0.183752, 0.471826, 0.610241],[0.180500, 0.464474, 0.600622],[0.177365, 0.457117, 0.591037],[0.174264, 0.449788, 0.581483],[0.171224, 0.442474, 0.571966],[0.168242, 0.435172, 0.562486],[0.165292, 0.427884, 0.553021],[0.162439, 0.420608, 0.543603],[0.159545, 0.413370, 0.534210],[0.156739, 0.406147, 0.524856],[0.153905, 0.398932, 0.515524],[0.151122, 0.391757, 0.506230],[0.148346, 0.384591, 0.496972],[0.145641, 0.377462, 0.487751],[0.142879, 0.370343, 0.478544],[0.140138, 0.363257, 0.469389],[0.137466, 0.356204, 0.460239],[0.134777, 0.349162, 0.451147],[0.132079, 0.342150, 0.442085],[0.129401, 0.335173, 0.433042],[0.126735, 0.328195, 0.424036],[0.124090, 0.321259, 0.415071],[0.121456, 0.314347, 0.406144],[0.118899, 0.307460, 0.397234],[0.116316, 0.300608, 0.388376],[0.113731, 0.293781, 0.379546],[0.111187, 0.286980, 0.370748],[0.108613, 0.280217, 0.362004],[0.106159, 0.273497, 0.353280],[0.103670, 0.266776, 0.344594],[0.101183, 0.260108, 0.335952],[0.098776, 0.253467, 0.327342],[0.096347, 0.246850, 0.318783],[0.094059, 0.240264, 0.310267],[0.091788, 0.233727, 0.301758],[0.089506, 0.227245, 0.293318],[0.087341, 0.220800, 0.284914],[0.085142, 0.214360, 0.276576],[0.083069, 0.207981, 0.268249],[0.081098, 0.201631, 0.259992],[0.079130, 0.195361, 0.251781],[0.077286, 0.189136, 0.243589],[0.075571, 0.182943, 0.235502],[0.073993, 0.176835, 0.227434],[0.072410, 0.170785, 0.219433],[0.071045, 0.164795, 0.211500],[0.069767, 0.158901, 0.203628],[0.068618, 0.153040, 0.195818],[0.067560, 0.147319, 0.188124],[0.066665, 0.141671, 0.180452],[0.065923, 0.136076, 0.172917],[0.065339, 0.130695, 0.165458],[0.064911, 0.125349, 0.158169],[0.064636, 0.120132, 0.150946],[0.064517, 0.115070, 0.143889],[0.064554, 0.110222, 0.136957],[0.064749, 0.105427, 0.130230],[0.065100, 0.100849, 0.123569],[0.065383, 0.096469, 0.117170],[0.065574, 0.092338, 0.111008],[0.065892, 0.088201, 0.104982],[0.066388, 0.084134, 0.099288],[0.067108, 0.080051, 0.093829],[0.068193, 0.076099, 0.088470],[0.069720, 0.072283, 0.083025],[0.071639, 0.068654, 0.077544],[0.073978, 0.065058, 0.072110],[0.076596, 0.061657, 0.066651],[0.079637, 0.058550, 0.061133],[0.082963, 0.055666, 0.055745],[0.086537, 0.052997, 0.050336],[0.090315, 0.050699, 0.045040],[0.094260, 0.048753, 0.039773],[0.098319, 0.047041, 0.034683],[0.102458, 0.045624, 0.030074],[0.106732, 0.044705, 0.026012],[0.110986, 0.043972, 0.022379],[0.115245, 0.043596, 0.019150],[0.119547, 0.043567, 0.016299],[0.123812, 0.043861, 0.013797],[0.128105, 0.044459, 0.011588],[0.132315, 0.045229, 0.009531],[0.136451, 0.046164, 0.007895],[0.140635, 0.047374, 0.006502],[0.144884, 0.048634, 0.005327],[0.149230, 0.049836, 0.004346],[0.153685, 0.050997, 0.003537],[0.158309, 0.052130, 0.002882],[0.163014, 0.053218, 0.002363],[0.167811, 0.054240, 0.001963],[0.172736, 0.055172, 0.001669],[0.177801, 0.056018, 0.001469],[0.182863, 0.056820, 0.001340],[0.188058, 0.057574, 0.001262],[0.193233, 0.058514, 0.001226],[0.198463, 0.059550, 0.001227],[0.203778, 0.060501, 0.001260],[0.209092, 0.061486, 0.001322],[0.214470, 0.062710, 0.001412],[0.219897, 0.063823, 0.001529],[0.225345, 0.065027, 0.001675],[0.230856, 0.066297, 0.001853],[0.236422, 0.067645, 0.002068],[0.242016, 0.069092, 0.002325],[0.247681, 0.070458, 0.002632],[0.253390, 0.071986, 0.002998],[0.259176, 0.073640, 0.003435],[0.264997, 0.075237, 0.003955],[0.270934, 0.076965, 0.004571],[0.276928, 0.078822, 0.005301],[0.283017, 0.080819, 0.006161],[0.289196, 0.082879, 0.007171],[0.295466, 0.085075, 0.008349],[0.301858, 0.087460, 0.009726],[0.308387, 0.089912, 0.011455],[0.315024, 0.092530, 0.013324],[0.321806, 0.095392, 0.015413],[0.328738, 0.098396, 0.017780],[0.335805, 0.101580, 0.020449],[0.343036, 0.104977, 0.023440],[0.350413, 0.108640, 0.026771],[0.357947, 0.112564, 0.030456],[0.365629, 0.116658, 0.034571],[0.373470, 0.120971, 0.039115],[0.381463, 0.125606, 0.043693],[0.389583, 0.130457, 0.048471],[0.397845, 0.135474, 0.053136],[0.406220, 0.140795, 0.057848],[0.414690, 0.146274, 0.062715],[0.423229, 0.151979, 0.067685],[0.431837, 0.157906, 0.073044],[0.440444, 0.164028, 0.078620],[0.449085, 0.170269, 0.084644],[0.457704, 0.176666, 0.090869],[0.466314, 0.183213, 0.097335],[0.474900, 0.189888, 0.104064],[0.483420, 0.196677, 0.111039],[0.491910, 0.203516, 0.118190],[0.500322, 0.210433, 0.125501],[0.508690, 0.217425, 0.132983],[0.516977, 0.224432, 0.140623],[0.525197, 0.231543, 0.148349],[0.533349, 0.238624, 0.156261],[0.541440, 0.245755, 0.164233],[0.549481, 0.252923, 0.172265],[0.557462, 0.260091, 0.180403],[0.565378, 0.267255, 0.188640],[0.573272, 0.274461, 0.196924],[0.581112, 0.281673, 0.205237],[0.588920, 0.288894, 0.213625],[0.596716, 0.296114, 0.222054],[0.604484, 0.303345, 0.230529],[0.612228, 0.310617, 0.239052],[0.619976, 0.317867, 0.247618],[0.627708, 0.325132, 0.256189],[0.635438, 0.332443, 0.264815],[0.643173, 0.339745, 0.273490],[0.650917, 0.347064, 0.282179],[0.658661, 0.354395, 0.290887],[0.666419, 0.361751, 0.299640],[0.674194, 0.369121, 0.308415],[0.681975, 0.376518, 0.317219],[0.689783, 0.383920, 0.326043],[0.697596, 0.391354, 0.334929],[0.705434, 0.398794, 0.343796],[0.713288, 0.406271, 0.352720],[0.721158, 0.413757, 0.361662],[0.729054, 0.421259, 0.370618],[0.736968, 0.428796, 0.379616],[0.744900, 0.436349, 0.388639],[0.752851, 0.443923, 0.397680],[0.760831, 0.451512, 0.406747],[0.768821, 0.459124, 0.415838],[0.776844, 0.466756, 0.424962],[0.784879, 0.474407, 0.434092],[0.792935, 0.482080, 0.443269],[0.801009, 0.489763, 0.452465],[0.809110, 0.497486, 0.461672],[0.817222, 0.505207, 0.470910],[0.825358, 0.512962, 0.480170],[0.833517, 0.520732, 0.489445],[0.841692, 0.528527, 0.498763],[0.849885, 0.536335, 0.508096],[0.858092, 0.544161, 0.517448],[0.866324, 0.552013, 0.526825],[0.874568, 0.559879, 0.536218],[0.882829, 0.567761, 0.545643],[0.891110, 0.575670, 0.555082],[0.899407, 0.583585, 0.564550],[0.907716, 0.591530, 0.574038],[0.916031, 0.599492, 0.583552],[0.924368, 0.607473, 0.593095],[0.932714, 0.615460, 0.602649],[0.941076, 0.623483, 0.612229],[0.949447, 0.631512, 0.621832],[0.957832, 0.639563, 0.631467],[0.966219, 0.647628, 0.641113],[0.974619, 0.655718, 0.650792],[0.983030, 0.663823, 0.660487],[0.991448, 0.671939, 0.670216],[0.999873, 0.680072, 0.679950]] +#Taken from F. Crameri's scientific-colour-maps version 8.0.1. https://doi.org/10.5281/zenodo.1243862 (included by default in newer matplotlib versions) +berlin = LinearSegmentedColormap.from_list('berlin', berlin_cmap, N=256) +max_iterations = 1e6 + + +app_ui = ui.page_sidebar( + ui.sidebar( + ui.input_slider( + "velocity", + "Velocity (units of c)", + min=-0.9, + max=0.9, + value=0.5, + step=0.01, + animate=True + ), + ui.input_slider("period", "Period of signal (self time)", min=1, max=3, value=2, step=0.1, animate=True), + ui.input_radio_buttons("signal_type", "Signal type", choices={"ping": "ping", "periodic": "periodic"}, selected="periodic", inline=True), + ui.input_dark_mode(id='dark_mode'), + open='always' + ), + ui.output_plot( + "plot", + click=True, + width="100%", height="700px" + ), +) + + +def lorentz_transform(x, t, v): + c = 1 # velocity of light v will be in units of c + gamma = 1 / (1 - (v**2 / c**2))**0.5 + x_prime = gamma * (x + v * t) + t_prime = gamma * (t + (v * x) / c**2) + return x_prime, t_prime + + +def mirrored(maxval, inc=1): + x = np.arange(inc, maxval, inc) + return np.r_[-x[::-1], 0, x] + + +def server(input, output, session): + tlimits = (-10, 10) + points_per_period = 12 + color_norm = Normalize(vmin=-1, vmax=1) + + @render.plot() + def plot(): + if input.dark_mode() == "dark": + style_label = 'dark_background' + blue = 'lightsteelblue' + red = 'lightcoral' + cmap = berlin + else: + style_label = 'seaborn-v0_8' + blue = 'navy' + red = 'firebrick' + cmap = plt.get_cmap('RdBu_r') + + with plt.style.context(style_label): + fig, ax = plt.subplot_mosaic([['minkowski', 'sender'], ['minkowski', 'receiver']]) + + v = input.velocity() + period = input.period() + + signal_times = mirrored(tlimits[1], inc=period) + if input.signal_type() == 'ping': + signal_amplitude = np.ones_like(signal_times) + else: # input.signal_type() == 'periodic': + signal_times = np.linspace(signal_times[0], signal_times[-1], points_per_period * (len(signal_times) - 1) + 1) + signal_amplitude = np.cos(2 * np.pi * signal_times / period) + + for indx, signal in enumerate(signal_times): + + if input.signal_type() == 'ping': + color = 'grey' + else: # input.signal_type() == 'periodic': + # Color line according to amplitude by mapping to colormap between -1 and 1 + color = cmap(color_norm(signal_amplitude[indx])) + + ax['sender'].plot([signal, signal], [0, signal_amplitude[indx]], color=color) + + ct_axis_transformed, x_axis_transformed = lorentz_transform(signal, 0, v) + + sign = np.sign(signal) * np.sign(v) + + intersection = ct_axis_transformed + sign * x_axis_transformed + + # Extend signal to receiver worldline with slope 1 + ax['minkowski'].plot([x_axis_transformed , 0], [ct_axis_transformed, intersection], color=color) + ax['receiver'].plot([intersection, intersection], [0, signal_amplitude[indx]], color=color) + + + ax['minkowski'].axvline(0, color=blue, label='Reciever worldline') + ax['minkowski'].axline((0, 0), slope=np.tan(np.pi/2 - np.atan(v)), color=red, label='Sender worldline') + + ax['sender'].axhline(0, color=red) + ax['sender'].axvline(0, color=red) + ax['receiver'].axhline(0, color=blue) + ax['receiver'].axvline(0, color=blue) + + ax['sender'].set_title("Sender Frame (in motion)") + ax['receiver'].set_title("Receiver Frame (at rest)") + ax['minkowski'].set_title("Minkowski Spacetime Diagram") + + ax['sender'].set_xlabel("sender self time in seconds") + ax['sender'].set_ylabel("signal amplitude") + ax['receiver'].set_xlabel("receiver self time in seconds") + ax['receiver'].set_ylabel("signal amplitude") + # + ax['minkowski'].set_xlabel("x in light-seconds") + ax['minkowski'].set_ylabel("ct in light-seconds") + ax['minkowski'].set_aspect('equal') + + # # Add self time ticks according to x lims of sender and receiver plots + # sender_xticks = ax['sender'].get_xticks() + # for time in sender_xticks: + # sign = np.sign(time) * np.sign(v) + # ct_axis_transformed, x_axis_transformed = lorentz_transform(time, 0, v) + # if sign < 0: + # ax['minkowski'].text(x_axis_transformed, ct_axis_transformed, f'τ={time:.0f}s ', color=red, verticalalignment='center', + # horizontalalignment='right', ) + # else: + # ax['minkowski'].text(x_axis_transformed, ct_axis_transformed, f' τ={time:.0f}s', color=red, verticalalignment='center', + # horizontalalignment='left') + # + # receiver_xticks = ax['receiver'].get_xticks() + # for time in receiver_xticks: + # sign = np.sign(time) * np.sign(v) + # if sign > 0: + # ax['minkowski'].text(0, time, f'τ={time:.0f}s ', color=blue, verticalalignment='center', + # horizontalalignment='right', ) + # else: + # ax['minkowski'].text(0, time, f' τ={time:.0f}s', color=blue, verticalalignment='center', + # horizontalalignment='left') + + xlims = ax['minkowski'].get_xlim() + new_lim = max(5, xlims[1]) + ax['minkowski'].set_xlim(-new_lim, new_lim) + return fig + + +app = App(app_ui, server, debug=True) diff --git a/docs/apps/MinkowskiSpaceTime/requirements.txt b/docs/apps/MinkowskiSpaceTime/doppler/requirements.txt similarity index 100% rename from docs/apps/MinkowskiSpaceTime/requirements.txt rename to docs/apps/MinkowskiSpaceTime/doppler/requirements.txt diff --git a/docs/apps/MinkowskiSpaceTime/motion/app.py b/docs/apps/MinkowskiSpaceTime/motion/app.py new file mode 100644 index 0000000..f3b4173 --- /dev/null +++ b/docs/apps/MinkowskiSpaceTime/motion/app.py @@ -0,0 +1,291 @@ +import matplotlib.pyplot as plt +from matplotlib.colors import LinearSegmentedColormap +import numpy as np +import pandas as pd +from shiny import App, render, ui, reactive + +berlin_cmap = [[0.621082, 0.690182, 0.999507],[0.612157, 0.689228, 0.995374],[0.603202, 0.688250, 0.991239],[0.594200, 0.687257, 0.987092],[0.585165, 0.686248, 0.982922],[0.576088, 0.685222, 0.978733],[0.566961, 0.684166, 0.974524],[0.557791, 0.683098, 0.970288],[0.548590, 0.681992, 0.966016],[0.539327, 0.680859, 0.961704],[0.530034, 0.679691, 0.957350],[0.520687, 0.678484, 0.952942],[0.511295, 0.677230, 0.948466],[0.501863, 0.675908, 0.943923],[0.492368, 0.674526, 0.939297],[0.482832, 0.673075, 0.934574],[0.473239, 0.671530, 0.929751],[0.463610, 0.669898, 0.924806],[0.453931, 0.668152, 0.919735],[0.444213, 0.666275, 0.914518],[0.434440, 0.664271, 0.909136],[0.424645, 0.662120, 0.903586],[0.414818, 0.659791, 0.897845],[0.404975, 0.657289, 0.891905],[0.395137, 0.654579, 0.885750],[0.385296, 0.651674, 0.879368],[0.375493, 0.648536, 0.872757],[0.365742, 0.645164, 0.865903],[0.356059, 0.641552, 0.858801],[0.346453, 0.637692, 0.851451],[0.336982, 0.633574, 0.843855],[0.327642, 0.629189, 0.836017],[0.318487, 0.624551, 0.827937],[0.309539, 0.619657, 0.819628],[0.300784, 0.614497, 0.811108],[0.292309, 0.609115, 0.802379],[0.284098, 0.603485, 0.793470],[0.276205, 0.597634, 0.784386],[0.268595, 0.591580, 0.775143],[0.261308, 0.585335, 0.765780],[0.254368, 0.578908, 0.756296],[0.247753, 0.572328, 0.746719],[0.241464, 0.565596, 0.737066],[0.235515, 0.558748, 0.727351],[0.229842, 0.551802, 0.717600],[0.224503, 0.544750, 0.707805],[0.219485, 0.537628, 0.697998],[0.214694, 0.530433, 0.688190],[0.210172, 0.523193, 0.678377],[0.205889, 0.515897, 0.668578],[0.201771, 0.508598, 0.658787],[0.197878, 0.501258, 0.649030],[0.194172, 0.493903, 0.639287],[0.190556, 0.486541, 0.629572],[0.187112, 0.479181, 0.619898],[0.183752, 0.471826, 0.610241],[0.180500, 0.464474, 0.600622],[0.177365, 0.457117, 0.591037],[0.174264, 0.449788, 0.581483],[0.171224, 0.442474, 0.571966],[0.168242, 0.435172, 0.562486],[0.165292, 0.427884, 0.553021],[0.162439, 0.420608, 0.543603],[0.159545, 0.413370, 0.534210],[0.156739, 0.406147, 0.524856],[0.153905, 0.398932, 0.515524],[0.151122, 0.391757, 0.506230],[0.148346, 0.384591, 0.496972],[0.145641, 0.377462, 0.487751],[0.142879, 0.370343, 0.478544],[0.140138, 0.363257, 0.469389],[0.137466, 0.356204, 0.460239],[0.134777, 0.349162, 0.451147],[0.132079, 0.342150, 0.442085],[0.129401, 0.335173, 0.433042],[0.126735, 0.328195, 0.424036],[0.124090, 0.321259, 0.415071],[0.121456, 0.314347, 0.406144],[0.118899, 0.307460, 0.397234],[0.116316, 0.300608, 0.388376],[0.113731, 0.293781, 0.379546],[0.111187, 0.286980, 0.370748],[0.108613, 0.280217, 0.362004],[0.106159, 0.273497, 0.353280],[0.103670, 0.266776, 0.344594],[0.101183, 0.260108, 0.335952],[0.098776, 0.253467, 0.327342],[0.096347, 0.246850, 0.318783],[0.094059, 0.240264, 0.310267],[0.091788, 0.233727, 0.301758],[0.089506, 0.227245, 0.293318],[0.087341, 0.220800, 0.284914],[0.085142, 0.214360, 0.276576],[0.083069, 0.207981, 0.268249],[0.081098, 0.201631, 0.259992],[0.079130, 0.195361, 0.251781],[0.077286, 0.189136, 0.243589],[0.075571, 0.182943, 0.235502],[0.073993, 0.176835, 0.227434],[0.072410, 0.170785, 0.219433],[0.071045, 0.164795, 0.211500],[0.069767, 0.158901, 0.203628],[0.068618, 0.153040, 0.195818],[0.067560, 0.147319, 0.188124],[0.066665, 0.141671, 0.180452],[0.065923, 0.136076, 0.172917],[0.065339, 0.130695, 0.165458],[0.064911, 0.125349, 0.158169],[0.064636, 0.120132, 0.150946],[0.064517, 0.115070, 0.143889],[0.064554, 0.110222, 0.136957],[0.064749, 0.105427, 0.130230],[0.065100, 0.100849, 0.123569],[0.065383, 0.096469, 0.117170],[0.065574, 0.092338, 0.111008],[0.065892, 0.088201, 0.104982],[0.066388, 0.084134, 0.099288],[0.067108, 0.080051, 0.093829],[0.068193, 0.076099, 0.088470],[0.069720, 0.072283, 0.083025],[0.071639, 0.068654, 0.077544],[0.073978, 0.065058, 0.072110],[0.076596, 0.061657, 0.066651],[0.079637, 0.058550, 0.061133],[0.082963, 0.055666, 0.055745],[0.086537, 0.052997, 0.050336],[0.090315, 0.050699, 0.045040],[0.094260, 0.048753, 0.039773],[0.098319, 0.047041, 0.034683],[0.102458, 0.045624, 0.030074],[0.106732, 0.044705, 0.026012],[0.110986, 0.043972, 0.022379],[0.115245, 0.043596, 0.019150],[0.119547, 0.043567, 0.016299],[0.123812, 0.043861, 0.013797],[0.128105, 0.044459, 0.011588],[0.132315, 0.045229, 0.009531],[0.136451, 0.046164, 0.007895],[0.140635, 0.047374, 0.006502],[0.144884, 0.048634, 0.005327],[0.149230, 0.049836, 0.004346],[0.153685, 0.050997, 0.003537],[0.158309, 0.052130, 0.002882],[0.163014, 0.053218, 0.002363],[0.167811, 0.054240, 0.001963],[0.172736, 0.055172, 0.001669],[0.177801, 0.056018, 0.001469],[0.182863, 0.056820, 0.001340],[0.188058, 0.057574, 0.001262],[0.193233, 0.058514, 0.001226],[0.198463, 0.059550, 0.001227],[0.203778, 0.060501, 0.001260],[0.209092, 0.061486, 0.001322],[0.214470, 0.062710, 0.001412],[0.219897, 0.063823, 0.001529],[0.225345, 0.065027, 0.001675],[0.230856, 0.066297, 0.001853],[0.236422, 0.067645, 0.002068],[0.242016, 0.069092, 0.002325],[0.247681, 0.070458, 0.002632],[0.253390, 0.071986, 0.002998],[0.259176, 0.073640, 0.003435],[0.264997, 0.075237, 0.003955],[0.270934, 0.076965, 0.004571],[0.276928, 0.078822, 0.005301],[0.283017, 0.080819, 0.006161],[0.289196, 0.082879, 0.007171],[0.295466, 0.085075, 0.008349],[0.301858, 0.087460, 0.009726],[0.308387, 0.089912, 0.011455],[0.315024, 0.092530, 0.013324],[0.321806, 0.095392, 0.015413],[0.328738, 0.098396, 0.017780],[0.335805, 0.101580, 0.020449],[0.343036, 0.104977, 0.023440],[0.350413, 0.108640, 0.026771],[0.357947, 0.112564, 0.030456],[0.365629, 0.116658, 0.034571],[0.373470, 0.120971, 0.039115],[0.381463, 0.125606, 0.043693],[0.389583, 0.130457, 0.048471],[0.397845, 0.135474, 0.053136],[0.406220, 0.140795, 0.057848],[0.414690, 0.146274, 0.062715],[0.423229, 0.151979, 0.067685],[0.431837, 0.157906, 0.073044],[0.440444, 0.164028, 0.078620],[0.449085, 0.170269, 0.084644],[0.457704, 0.176666, 0.090869],[0.466314, 0.183213, 0.097335],[0.474900, 0.189888, 0.104064],[0.483420, 0.196677, 0.111039],[0.491910, 0.203516, 0.118190],[0.500322, 0.210433, 0.125501],[0.508690, 0.217425, 0.132983],[0.516977, 0.224432, 0.140623],[0.525197, 0.231543, 0.148349],[0.533349, 0.238624, 0.156261],[0.541440, 0.245755, 0.164233],[0.549481, 0.252923, 0.172265],[0.557462, 0.260091, 0.180403],[0.565378, 0.267255, 0.188640],[0.573272, 0.274461, 0.196924],[0.581112, 0.281673, 0.205237],[0.588920, 0.288894, 0.213625],[0.596716, 0.296114, 0.222054],[0.604484, 0.303345, 0.230529],[0.612228, 0.310617, 0.239052],[0.619976, 0.317867, 0.247618],[0.627708, 0.325132, 0.256189],[0.635438, 0.332443, 0.264815],[0.643173, 0.339745, 0.273490],[0.650917, 0.347064, 0.282179],[0.658661, 0.354395, 0.290887],[0.666419, 0.361751, 0.299640],[0.674194, 0.369121, 0.308415],[0.681975, 0.376518, 0.317219],[0.689783, 0.383920, 0.326043],[0.697596, 0.391354, 0.334929],[0.705434, 0.398794, 0.343796],[0.713288, 0.406271, 0.352720],[0.721158, 0.413757, 0.361662],[0.729054, 0.421259, 0.370618],[0.736968, 0.428796, 0.379616],[0.744900, 0.436349, 0.388639],[0.752851, 0.443923, 0.397680],[0.760831, 0.451512, 0.406747],[0.768821, 0.459124, 0.415838],[0.776844, 0.466756, 0.424962],[0.784879, 0.474407, 0.434092],[0.792935, 0.482080, 0.443269],[0.801009, 0.489763, 0.452465],[0.809110, 0.497486, 0.461672],[0.817222, 0.505207, 0.470910],[0.825358, 0.512962, 0.480170],[0.833517, 0.520732, 0.489445],[0.841692, 0.528527, 0.498763],[0.849885, 0.536335, 0.508096],[0.858092, 0.544161, 0.517448],[0.866324, 0.552013, 0.526825],[0.874568, 0.559879, 0.536218],[0.882829, 0.567761, 0.545643],[0.891110, 0.575670, 0.555082],[0.899407, 0.583585, 0.564550],[0.907716, 0.591530, 0.574038],[0.916031, 0.599492, 0.583552],[0.924368, 0.607473, 0.593095],[0.932714, 0.615460, 0.602649],[0.941076, 0.623483, 0.612229],[0.949447, 0.631512, 0.621832],[0.957832, 0.639563, 0.631467],[0.966219, 0.647628, 0.641113],[0.974619, 0.655718, 0.650792],[0.983030, 0.663823, 0.660487],[0.991448, 0.671939, 0.670216],[0.999873, 0.680072, 0.679950]] +#Taken from F. Crameri's scientific-colour-maps version 8.0.1. https://doi.org/10.5281/zenodo.1243862 (included by default in newer matplotlib versions) +berlin = LinearSegmentedColormap.from_list('berlin', berlin_cmap, N=256) +max_iterations = 1e6 + +# Generate random events for background defined in stationary frame +random_generator = np.random.default_rng(42) + +# num_random_events = int(1e3) # Reasonable for rest frame +# random_t = random_generator.uniform(-30, 200, num_random_events) +# random_x = random_generator.uniform(-30, 30, num_random_events) + + +num_random_events = int(1e6) +random_t = random_generator.uniform(-1000, 1000, num_random_events) +random_x = random_generator.uniform(-1000, 1000, num_random_events) +random_points = np.array([random_x, random_t]) + +app_ui = ui.page_sidebar( + ui.sidebar( + ui.input_radio_buttons("frame_of_reference", "Frame of reference", choices={"rest": "rest", "motion": "motion"}, selected="rest", inline=True), + ui.input_slider("frame", "Self time", min=0, max=1, value=0, step=0.5, animate=True), + ui.input_slider("turning_point", "Turning point (in light-seconds)", min=1, max=20, value=10, step=0.1), + ui.input_slider("acceleration", "Proper acceleration (in c/s)", min=0.01, max=0.25, value=0.15, step=0.001), + ui.input_checkbox('show_light_cones', 'Show light cones', value=True), + ui.input_dark_mode(id='dark_mode'), + ui.accordion( + ui.accordion_panel('Advanced settings', + ui.input_slider("dtau", "Proper time step (log10)", min=-5, max=-1, value=-2, step=0.1), + ), + open=False, + ), + open='always' + ), + ui.output_plot( + "plot", + click=True, + width="100%", height="700px" + ), +) + + +def hypTStep(dt, v0, x0, tau0, g): + ## Hyperbolic step. + ## If an object has proper acceleration g and starts at position x0 with speed v0 and proper time tau0 + ## as seen from an inertial frame, then return the new v, x, tau after time dt has elapsed. + if g == 0: + return v0, x0 + v0 * dt, tau0 + dt * (1. - v0 ** 2) ** 0.5 + + tinit = v0 / (g * (1 - v0 ** 2) ** 0.5) + B = (1 + (g ** 2 * (dt + tinit) ** 2)) ** 0.5 + + v1 = g * (dt + tinit) / B + x1 = x0 + (1.0 / g) * (B - 1. / (1. - v0 ** 2) ** 0.5) + tau1 = tau0 + (np.arcsinh(g * (dt + tinit)) - np.arcsinh(g * tinit)) / g + return v1, x1, tau1 + +def tauStep(dtau, v0, x0, t0, g): + ## linear step in proper time of clock. + ## If an object has proper acceleration g and starts at position x0 with speed v0 at time t0 + ## as seen from an inertial frame, then return the new v, x, t after proper time dtau has elapsed. + + ## Compute how much t will change given a proper-time step of dtau + gamma = (1. - v0 ** 2) ** -0.5 + if g == 0: + dt = dtau * gamma + else: + v0g = v0 * gamma + dt = (np.sinh(dtau * g + np.arcsinh(v0g)) - v0g) / g + + # return v0 + dtau * g, x0 + v0*dt, t0 + dt + v1, x1, t1 = hypTStep(dt, v0, x0, t0, g) + return v1, x1, t0 + dt + +def lorentz_transform(x, t, v): + c = 1 # velocity of light v will be in units of c + gamma = 1 / (1 - (v**2 / c**2))**0.5 + x_prime = gamma * (x + v * t) + t_prime = gamma * (t + (v * x) / c**2) + return x_prime, t_prime + +def ceil_to_next_multiple(value, multiple): + return np.ceil(value / multiple) * multiple + +def server(input, output, session): + simulation_data = reactive.Value(pd.DataFrame()) + random_points = reactive.Value(np.array([0, 0])) + + @reactive.effect + def _(): + dtau = 10 ** input.dtau() + g = input.acceleration() # proper acceleration of the moving observer + + with ui.Progress(min=0, max=1) as p: + p.set(message="Calculation in progress") + # Initialize arrays to hold data + times = [0.0] + positions = [0.0] + velocities = [0.0] + proper_times = [0.0] + + reached_first_turning_point = False + reached_second_turning_point = False + # Simulate motion of the accelerating observer + count = 0 + while g != 0 and count < max_iterations: + v_new, x_new, t_new = tauStep(dtau, velocities[-1], positions[-1], times[-1], g) + if x_new >= input.turning_point()/2 and not reached_first_turning_point and not reached_second_turning_point: + reached_first_turning_point = True + g = -g # Reverse acceleration at turning point + elif x_new <= input.turning_point()/2 and reached_first_turning_point and not reached_second_turning_point: + reached_second_turning_point = True + g = -g # Reverse acceleration to head back toward origin + if x_new <= 0 and reached_second_turning_point: + g = 0 # Stop acceleration at origin + v_new, x_new, t_new = [0.0, 0.0, times[-1] + (dtau * (1. - velocities[-1] ** 2) ** 0.5)] + max_time = t_new + ceiled_max_time = ceil_to_next_multiple(max_time, 5) + + velocities.append(v_new) + positions.append(x_new) + times.append(t_new) + proper_times.append(proper_times[-1] + dtau) + count += 1 + + if count >= max_iterations: + ui.notification_show("Maximum number of iterations reached. Consider increasing the proper time step.", type="error") + return + + # Keep appending zeros until max_time is reached + while proper_times[-1] < ceiled_max_time + dtau: + velocities.append(0.0) + positions.append(0.0) + times.append(times[-1] + dtau) + proper_times.append(proper_times[-1] + dtau) + + final_count = len(times) + + # Fill trivial case of no acceleration + stationary_times = np.arange(0, final_count, 1) * dtau + stationary_positions = np.zeros(final_count) + stationary_velocities = np.zeros(final_count) + stationary_proper_times = np.arange(0, final_count, 1) * dtau + + data = pd.DataFrame({ + 'm_time': times, + 'm_position': positions, + 'm_velocity': velocities, + 'm_proper_time': proper_times, + 's_time': stationary_times, + 's_position': stationary_positions, + 's_velocity': stationary_velocities, + 's_proper_time': stationary_proper_times, + }) + ui.update_slider("frame", max=ceiled_max_time, value=0) + simulation_data.set(data) + + @render.plot() + def plot(): + if input.dark_mode() == "dark": + style_label = 'dark_background' + blue = 'lightsteelblue' + red = 'lightcoral' + cmap = berlin + else: + style_label = 'default' + blue = 'navy' + red = 'firebrick' + cmap = 'RdBu_r' + + with plt.style.context(style_label): + fig, ax = plt.subplots() + data = simulation_data() + + # Draw axes of moving observer of current frame + frametime = input.frame() + frame_data = lambda: data.iloc[(data['m_proper_time'] - frametime).abs().argsort()[:1]].squeeze() + v_frame = frame_data()['m_velocity'] + times = np.arange(0, ceil_to_next_multiple(data['m_proper_time'].max(), 5) + 1, 5) + + if input.frame_of_reference() == 'rest': + shift_pos, shift_time = frame_data()['s_position'], frame_data()['s_time'] + + ax.scatter(data['m_position'] - shift_pos, data['m_time'] - shift_time, c=data['m_velocity'], vmin=-1, vmax=1, cmap=cmap) + ax.plot(data['m_position'] - shift_pos, data['m_time'] - shift_time, color='grey') + ax.plot(data['s_position'] - shift_pos, data['s_time'] - shift_time, color=blue) + + cbar = fig.colorbar(plt.cm.ScalarMappable(cmap=cmap, norm=plt.Normalize(vmin=-1, vmax=1)), ax=ax) + + ax.scatter([frame_data()['m_position'] - shift_pos], [frame_data()['m_time'] - shift_time], color='grey', marker='o', zorder=5) + ax.scatter([frametime - shift_time], [0 - shift_pos], color='grey', marker='o', zorder=5) + + if input.show_light_cones(): + ax.axline([frame_data()['m_position'] - shift_pos, frame_data()['m_time'] - shift_time], color=red, slope=1, linestyle='--', alpha=0.5) + ax.axline([frame_data()['m_position'] - shift_pos, frame_data()['m_time'] - shift_time], color=red, slope=-1, linestyle='--', alpha=0.5) + + ax.axline([frametime - shift_time, 0 - shift_pos], color=blue, slope=1, linestyle='--', alpha=0.5) + ax.axline([frametime - shift_time, 0 - shift_pos], color=blue, slope=-1, linestyle='--', alpha=0.5) + + # # Draw axes + # ax.axline((frame_data()['m_position'] - shift_pos, frame_data()['m_time'] - shift_time), slope=frame_data()['m_velocity'], color=red) + # ax.axline((frame_data()['m_position'] - shift_pos, frame_data()['m_time'] - shift_time), slope=np.tan(np.pi/2 - np.atan(frame_data()['m_velocity'])), color=red) + # ax.axhline(frametime - shift_time, color=blue) + # ax.axvline(0 - shift_pos, color=blue) + + # ax.axvline(input.turning_point() - shift_pos, color='grey', linestyle='--', zorder=-1) + + for prefix, color in [('s', blue), ('m', red)]: + # Interpolate positions of moving observer at these times to plot self times + times_at_this_proper_time = np.interp(times, data[prefix + '_proper_time'], data[prefix + '_time'], left=np.nan, right=np.nan) + positions_at_times = np.interp(times_at_this_proper_time, data[prefix + '_time'], data[prefix + '_position'], left=np.nan, right=np.nan) + for t, x, tau in zip(times_at_this_proper_time, positions_at_times, times): + if prefix is 's': + ax.text(x - shift_pos, t - shift_time, f'τ={tau:.0f}s ', color=color, verticalalignment='center', horizontalalignment='right', clip_on=True) + else: + ax.text(x - shift_pos, t - shift_time, f' τ={tau:.0f}s', color=color, verticalalignment='center', horizontalalignment='left', clip_on=True) + ax.scatter(x - shift_pos, t - shift_time, color=color, marker='x') + + # Plot random events in background + ax.scatter(random_x - shift_pos, random_t - shift_time, color='grey', s=1, alpha=0.5, zorder=-1) + + else: # frame of reference is motion + # For this we have to Lorentz transform all points into the moving frame at the selected self time + shift_pos, shift_time = frame_data()['m_position'], frame_data()['m_time'] + m_x_prime, m_t_prime = lorentz_transform(data['m_position'] - shift_pos, data['m_time'] - shift_time, -v_frame) + s_x_prime, s_t_prime = lorentz_transform(data['s_position'] - shift_pos, data['s_time'] - shift_time, -v_frame) + ax.scatter(m_x_prime, m_t_prime, c=data['m_velocity'], vmin=-1, vmax=1, cmap=cmap) + + cbar = fig.colorbar(plt.cm.ScalarMappable(cmap=cmap, norm=plt.Normalize(vmin=-1, vmax=1)), ax=ax) + ax.plot(m_x_prime, m_t_prime, color='grey') + ax.plot(s_x_prime, s_t_prime, color=blue) + + ax.scatter([0], [0], color='grey', marker='o', zorder=5) + stationary_time = np.interp(frametime, data['m_proper_time'], s_t_prime, left=np.nan, right=np.nan) + stationary_pos = np.interp(stationary_time, s_t_prime, s_x_prime, left=np.nan, right=np.nan) + + ax.scatter([stationary_pos], [stationary_time], color='grey', marker='o', zorder=5) + + if input.show_light_cones(): + ax.axline([0, 0], color=red, slope=1, linestyle='--', alpha=0.5) + ax.axline([0, 0], color=red, slope=-1, linestyle='--', alpha=0.5) + + ax.axline([stationary_pos, stationary_time], color=blue, slope=1, linestyle='--', alpha=0.5) + ax.axline([stationary_pos, stationary_time], color=blue, slope=-1, linestyle='--', alpha=0.5) + + # # Draw axes + # ax.axhline(0, color=red) + # ax.axvline(0, color=red) + + for (pos, time, proper_time), color in [((s_x_prime, s_t_prime, data['s_proper_time']), blue), ((m_x_prime, m_t_prime, data['m_proper_time']), red)]: + # Interpolate positions of moving observer at these times to plot self times + times_at_this_proper_time = np.interp(times, proper_time, time, left=np.nan, right=np.nan) + positions_at_times = np.interp(times_at_this_proper_time, time, pos, left=np.nan, right=np.nan) + for t, x, tau in zip(times_at_this_proper_time, positions_at_times, times): + if color is blue: # Shift not needed here, as the transformed data is already centered + ax.text(x, t, f'τ={tau:.0f}s ', color=color, verticalalignment='center', horizontalalignment='right', clip_on=True) + else: + ax.text(x, t, f' τ={tau:.0f}s', color=color, verticalalignment='center', horizontalalignment='left', clip_on=True) + ax.scatter(x, t, color=color, marker='x') + + # Also transform random points + r_x_prime, r_t_prime = lorentz_transform(random_x - shift_pos, random_t - shift_time, -v_frame) + ax.scatter(r_x_prime, r_t_prime, color='grey', s=1, alpha=0.5, zorder=-1) + + top_speed = np.max(np.abs(data['m_velocity'])) + cbar.ax.axhline(top_speed, color='grey', linestyle='--') + cbar.ax.text(0, top_speed, 'max ', color='grey', verticalalignment='center', horizontalalignment='right') + cbar.ax.axhline(-top_speed, color='grey', linestyle='--') + cbar.ax.text(0, -top_speed, 'max ', color='grey', verticalalignment='center', horizontalalignment='right') + cbar.ax.axhline(v_frame, color='grey') + + cbar.set_label('Velocity (in units of c)') + + if frametime == 0: + xlim = np.max([times[-1] * 0.5, input.turning_point() * 1.5]) + ax.set_xlim(-xlim, xlim) + ax.set_ylim(-5, times[-1] + 5) + else: + ax.set_xlim(-input.turning_point() * 1.3, input.turning_point() * 1.3) + ax.set_ylim(-input.turning_point() * 1.3, input.turning_point() * 1.3) + ax.axis('off') + ax.set_aspect('equal') + return fig + +app = App(app_ui, server, debug=True) diff --git a/docs/apps/MinkowskiSpaceTime/motion/requirements.txt b/docs/apps/MinkowskiSpaceTime/motion/requirements.txt new file mode 100644 index 0000000..fb0f9d3 --- /dev/null +++ b/docs/apps/MinkowskiSpaceTime/motion/requirements.txt @@ -0,0 +1,3 @@ +numpy +pandas +matplotlib \ No newline at end of file diff --git a/docs/apps/MinkowskiSpaceTime/app.py b/docs/apps/MinkowskiSpaceTime/stationary/app.py similarity index 100% rename from docs/apps/MinkowskiSpaceTime/app.py rename to docs/apps/MinkowskiSpaceTime/stationary/app.py diff --git a/docs/apps/MinkowskiSpaceTime/stationary/requirements.txt b/docs/apps/MinkowskiSpaceTime/stationary/requirements.txt new file mode 100644 index 0000000..fb0f9d3 --- /dev/null +++ b/docs/apps/MinkowskiSpaceTime/stationary/requirements.txt @@ -0,0 +1,3 @@ +numpy +pandas +matplotlib \ No newline at end of file diff --git a/docs/apps/NablaShowcase/app.md b/docs/apps/NablaShowcase/app.md index 1170076..d9d4180 100644 --- a/docs/apps/NablaShowcase/app.md +++ b/docs/apps/NablaShowcase/app.md @@ -18,4 +18,4 @@ hide: {{embed_app("100%", "830px", "gradient")}} ## Curl -{{embed_app("100%", "830px", "Curl")}} \ No newline at end of file +{{embed_app("100%", "830px", "curl")}} \ No newline at end of file diff --git a/docs/apps/OrbitIntegrator/app.md b/docs/apps/OrbitIntegrator/app.md new file mode 100644 index 0000000..8f08461 --- /dev/null +++ b/docs/apps/OrbitIntegrator/app.md @@ -0,0 +1,15 @@ +--- +authors: + - ptuemmler +categories: + - Physics +tags: + - Draft +date: 2025-11-23 +hide: + - toc +--- + +# Oribit Integrator + +{{embed_app("100%", "830px")}} diff --git a/docs/apps/OrbitIntegrator/app.py b/docs/apps/OrbitIntegrator/app.py new file mode 100644 index 0000000..4c87712 --- /dev/null +++ b/docs/apps/OrbitIntegrator/app.py @@ -0,0 +1,328 @@ +import numpy as np +import plotly.graph_objects as go +from shiny import App, Inputs, Outputs, Session, render, ui, reactive +from copy import deepcopy + + +app_ui = ui.page_sidebar( + ui.sidebar( + ui.input_select('ref_frame', 'Reference Frame', choices=['Heliocentric', 'Geocentric'], selected='Geocentric'), + ui.input_selectize('bodies', 'Bodies to Simulate', {'Sun': 'Sun', 'Moon': 'Moon'}, + selected=['Sun', 'Earth', 'Moon'], multiple=True), + ui.input_slider('time_step', 'Time Step (hours)', min=1, max=24, value=5, step=0.01), + ui.input_slider('total_time', 'Total Simulation Time (days)', min=0, max=365, value=5, step=0.005), + ui.input_slider('max_frames', 'Maximum Animation Frames', min=100, max=10000, value=1000, step=100), + ui.input_dark_mode(id='dark_mode'), + ), + ui.output_ui("plot", height='100%') +) + +class Body: + def __init__(self, mass: float, position, velocity, name: str=None): + self.name = name + self.mass = mass + self.pos = position + self.velo = velocity + + def return_vec(self): + return np.concatenate((self.pos, self.velo)) + + +class Simulation: + def __init__(self, bodies: list[Body]): + self.calc_diff_eqs: callable = None + self.diff_eq_kwargs: dict = {} + + self.times = None + self.history = None + self.bodies = bodies + self.N_bodies = len(self.bodies) + + self.quant_vec = np.concatenate(np.array([i.return_vec() for i in self.bodies])) + self.mass_vec = np.array([i.mass for i in self.bodies]) + self.name_vec = [i.name for i in self.bodies] + + def set_diff_eq(self, calc_diff_eqs: callable, **kwargs): + """ + Method which assigns an external solver function as the diff-eq solver for RK4. + For N-body or gravitational setups, this is the function which calculates accelerations. + --------------------------------- + Params: + calc_diff_eqs: A function which returns a [y] vector for RK4 + **kwargs: Any additional inputs/hyperparameters the external function requires + """ + self.diff_eq_kwargs = kwargs + self.calc_diff_eqs = calc_diff_eqs + + def rk4(self, t:float, dt:float): + """ + RK4 integrator. Calculates the K values and returns a new y vector + -------------------------------- + Params: + t: a time. + dt: timestep. + """ + k1 = dt * self.calc_diff_eqs(t,self.quant_vec,self.mass_vec,**self.diff_eq_kwargs) + k2 = dt * self.calc_diff_eqs(t + 0.5*dt,self.quant_vec+0.5*k1,self.mass_vec,**self.diff_eq_kwargs) + k3 = dt * self.calc_diff_eqs(t + 0.5*dt,self.quant_vec+0.5*k2,self.mass_vec,**self.diff_eq_kwargs) + k4 = dt * self.calc_diff_eqs(t + dt,self.quant_vec + k2,self.mass_vec,**self.diff_eq_kwargs) + + y_new = self.quant_vec + ((k1 + 2*k2 + 2*k3 + k4) / 6.0) + + return y_new + + + def run(self, T, dt, t0=0): + """ + Method which runs the simulation on a given set of bodies. + --------------------- + Params: + T: total time (in simulation units) to run the simulation. + dt: timestep (in simulation units) to advance the simulation. + t0 (optional): set a non-zero start time to the simulation. + + Returns: + None, but leaves an attribute history accessed via + 'simulation.history' which contains all y vectors for the simulation. + These are of shape (Nstep,Nbodies * 6), so the x and y positions of particle 1 are + simulation.history[:,0], simulation.history[:,1], while the same for particle 2 are + simulation.history[:,6], simulation.history[:,7]. Velocities are also extractable. + """ + + if self.calc_diff_eqs is None: + raise AttributeError('You must set a differential equation to solve first.') + + self.history = [self.quant_vec] + nsteps = int((T - t0) / dt) + for step in range(nsteps): + y_new = self.rk4(step, dt) + self.history.append(y_new) + self.quant_vec = y_new + self.history = np.array(self.history) + self.times = np.arange(nsteps) * dt + + def translate_history_to_dict(self): + """ + Method which translates the history array into a dictionary of arrays for easier access. + ---------------------- + Params: + None + Returns: + A dictionary where each key is the name of a body, and each value is another dictionary + with keys 'pos' and 'velo' containing the position and velocity arrays over time. + """ + + if self.history is None: + raise AttributeError('You must run the simulation first.') + history_dict = {} + for i, name in enumerate(self.name_vec): + history_dict[name] = {'pos': self.history[:, i*6:i*6+3], + 'velo': self.history[:, i*6+3:i*6+6]} + return history_dict + + +grav_constant = 6.67430e-11 # m^3 kg^-1 s^-2 +def nbody_gravity(t, y, masses, progress_callback=None): + """ + y: 1D array with length 6*N in the layout [x,y,z,vx,vy,vz, x,y,z,vx,vy,vz, ...] + masses: array-like of length N + """ + + y = np.asarray(y) + masses = np.asarray(masses) + N = len(masses) + assert y.size == 6 * N, "y must have length 6*N" + + # reshape into (N,6): columns 0:3 => pos, 3:6 => vel + state = y.reshape(N, 6) + pos = state[:, 0:3] # (N,3) + vel = state[:, 3:6] # (N,3) + + # pairwise displacement: r_ij = pos[i] - pos[j] + rij = pos[:, None, :] - pos[None, :, :] # (N, N, 3) + + # distances + dist2 = np.sum(rij * rij, axis=2) # (N, N) + + # compute |r|^3 with safety for self-terms + dist3 = dist2 * np.sqrt(dist2) # (N, N) + + # avoid division by zero on diagonal by setting diagonal to inf so contribution is 0 + np.fill_diagonal(dist3, np.inf) + + # acceleration: -G * sum_j m_j * r_ij / |r_ij|^3 + # einsum sums over j: result shape (N,3) + acc = -grav_constant * np.einsum('ijk,ij->ik', rij, 1.0 / dist3 * masses[None, :]) + + # assemble derivative in same interleaved format + dydt = np.empty_like(state) + dydt[:, 0:3] = vel + dydt[:, 3:6] = acc + + if progress_callback is not None: + progress_callback(t) + + return dydt.reshape(-1) + +def server(input: Inputs, output: Outputs, session: Session): + Sun = Body(mass=1.989e30, # kg + position=np.array([0, 0, 0]), # m + velocity=np.array([0, 0, 0]), # m/s + name='Sun') + + Earth = Body(mass=5.972e24, # kg + position=np.array([1.496e11, 0., 0.]), # m + velocity=np.array([0, 29780, 0]), # m/s + name='Earth') + + Moon = Body(mass=7.347e22, # kg + position=np.array([1.496e11 + 3.844e8, 0., 0.]), # m + velocity=np.array([0, 29780 + 1022, 0]), # m/s + name='Moon') + + bodies = [Sun, Earth, Moon] + + @reactive.calc + def sim(): + with ui.Progress(min=0, max=input.total_time() * 86400) as p: + p.set(message='Solving...') + selected_bodies = input.bodies() + simulation = Simulation([Earth, *[body for body in [Sun, Earth, Moon] if body.name in selected_bodies]]) + simulation.set_diff_eq(nbody_gravity) + simulation.run(T=input.total_time() * 86400, dt=input.time_step() * 3600) + ui.notification_show("Simulation complete!", type="message") + return simulation.translate_history_to_dict() + + @render.ui + def plot(): + # Set template based on dark mode + if input.dark_mode() == "dark": + template = "plotly_dark" + else: + template = "plotly_white" + + + result = deepcopy(sim()) + if input.ref_frame() == 'Geocentric': + earth_pos = result['Earth']['pos'] + for body_name, data in result.items(): + data['pos'] = data['pos'] - earth_pos + + fig = go.Figure() + + # Initialize empty scatter for animation + for body_name, data in result.items(): + fig.add_trace(go.Scatter3d( + x=[], y=[], z=[], + mode="markers", marker=dict(size=10), name=body_name + )) # current position + + + # Create frames for animation + frames = [] + N_steps = result[body_name]['pos'].shape[0] + N_bodies = len(result.keys()) + + if N_steps > input.max_frames(): + indices = np.linspace(0, N_steps, input.max_frames(), endpoint=False).astype(int) + else: + indices = list(range(N_steps)) + + with ui.Progress(min=0, max=len(indices)) as p: + p.set(message='Creating frames...') + for k in indices: + frame_data = [] + for body_name, data in result.items(): + frame_data.append(go.Scatter3d( + x=[data['pos'][k,0]], + y=[data['pos'][k,1]], + z=[data['pos'][k,2]] + )) + frames.append(go.Frame(data=frame_data, traces=list(range(N_bodies)), name=f'frame{k}')) + p.set(k) + + fig.update(frames=frames) + + for body_name, data in result.items(): + fig.add_trace(go.Scatter3d( + x=data['pos'][:,0], + y=data['pos'][:,1], + z=data['pos'][:,2], + mode="lines", line=dict(color="grey", width=2), showlegend=False, + )) # full trajectory + + def frame_args(duration): + return { + "frame": {"duration": duration}, + "mode": "immediate", + "fromcurrent": True, + "transition": {"duration": duration, "easing": "linear"}, + } + + + sliders = [ + {"pad": {"b": 10, "t": 60}, + "len": 0.9, + "x": 0.1, + "y": 0, + + "steps": [ + {"args": [[f.name], frame_args(0)], + "label": str(k), + "method": "animate", + } for k, f in enumerate(fig.frames) + ] + } + ] + + fig.update_layout( + updatemenus = [{"buttons":[ + { + "args": [None, frame_args(50)], + "label": "Play", + "method": "animate", + }, + { + "args": [[None], frame_args(0)], + "label": "Pause", + "method": "animate", + }], + + "direction": "left", + "pad": {"r": 10, "t": 70}, + "type": "buttons", + "x": 0.1, + "y": 0, + } + ], + sliders=sliders + ) + + # Collect max x and y values for setting axis limits + max_range = 0 + for body_name, data in result.items(): + pos = data['pos'] + max_range = max(max_range, + np.max(np.abs(pos[:,0])), + np.max(np.abs(pos[:,1])), + np.max(np.abs(pos[:,2]))) + axis_limit = max_range * 1.2 + + fig.update_layout( + scene=dict( + xaxis_title='X Position (m)', + yaxis_title='Y Position (m)', + zaxis_title='Z Position (m)', + xaxis_showspikes=False, + yaxis_showspikes=False, + zaxis_showspikes=False, + xaxis=dict(range=[-axis_limit, axis_limit]), + yaxis=dict(range=[-axis_limit, axis_limit]), + ), + height=700, + template=template, + sliders = sliders + ) + return ui.HTML(fig.to_html()) + +app = App(app_ui, server) diff --git a/docs/apps/OrbitIntegrator/requirements.txt b/docs/apps/OrbitIntegrator/requirements.txt new file mode 100644 index 0000000..5c84b3c --- /dev/null +++ b/docs/apps/OrbitIntegrator/requirements.txt @@ -0,0 +1,2 @@ +numpy +plotly diff --git a/docs/apps/TaylorExpansion/app.py b/docs/apps/TaylorExpansion/app.py index fabd317..55ca186 100644 --- a/docs/apps/TaylorExpansion/app.py +++ b/docs/apps/TaylorExpansion/app.py @@ -54,7 +54,6 @@ def server(input: Inputs, output: Outputs, session: Session): - @reactive.calc @reactive.calc def clean_function(): """Clean and validate the input function""" diff --git a/docs/contributing.md b/docs/contributing.md index a2c5610..a9a344b 100644 --- a/docs/contributing.md +++ b/docs/contributing.md @@ -109,7 +109,7 @@ If you have multiple apps in a single tool, as is the case in the example above, To test your tool locally, you should build the full documentation using the following command in the root directory of your repository: ```bash mkdocs build --clean -python -m http.server --directory ./site --bind localhost 8008 +python -m http.server --directory ./site --bind localhost 8008 ``` This will create a `site` directory containing the static HTML files for the documentation, including your new tool. You can then open your browser and navigate to [localhost http://[::1]:8008/](http://[::1]:8008/) to see the documentation and test if everything works.