{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Orbit plots in 3D\n", "\n", "Illustrate how 3D orbits can be animated. Note that the calculation of the orbits is purely for illustration purposes; they are not right! \n", "\n", "Some more information on animation is available here: https://alexgude.com/blog/matplotlib-blitting-supernova/" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Date and time 2023-03-22 15:24:53.788287\n", " \n", "Number of planets 7, time simulated 3.39\n", " \n", "Orbit 0\n", "Radius 0.60, angular freq. 13.52\n", "Euler angles -0.09, 0.05 and 0.00\n", "Start time 0.39\n", "Rotation matrix \n", " [[ 0.999 -0.001 0.053]\n", " [-0.003 0.996 0.086]\n", " [-0.053 -0.086 0.995]]\n", " \n", "Orbit 1\n", "Radius 0.80, angular freq. 8.78\n", "Euler angles -0.19, 0.02 and -0.00\n", "Start time 0.56\n", "Rotation matrix \n", " [[ 1. 0.001 0.017]\n", " [-0.004 0.983 0.184]\n", " [-0.016 -0.184 0.983]]\n", " \n", "Orbit 2\n", "Radius 1.10, angular freq. 5.45\n", "Euler angles -0.11, -0.08 and -0.00\n", "Start time 0.66\n", "Rotation matrix \n", " [[ 0.997 0.001 -0.081]\n", " [ 0.008 0.994 0.106]\n", " [ 0.081 -0.107 0.991]]\n", " \n", "Orbit 3\n", "Radius 1.20, angular freq. 4.78\n", "Euler angles 0.17, -0.09 and 0.00\n", "Start time 0.87\n", "Rotation matrix \n", " [[ 0.996 -0. -0.087]\n", " [-0.015 0.985 -0.17 ]\n", " [ 0.086 0.171 0.982]]\n", " \n", "Orbit 4\n", "Radius 1.50, angular freq. 3.42\n", "Euler angles -0.06, -0.17 and -0.00\n", "Start time 0.87\n", "Rotation matrix \n", " [[ 0.985 0. -0.172]\n", " [ 0.011 0.998 0.062]\n", " [ 0.172 -0.063 0.983]]\n", " \n", "Orbit 5\n", "Radius 1.90, angular freq. 2.40\n", "Euler angles -0.07, -0.02 and 0.00\n", "Start time 0.61\n", "Rotation matrix \n", " [[ 1. -0. -0.021]\n", " [ 0.002 0.997 0.072]\n", " [ 0.02 -0.072 0.997]]\n", " \n", "Orbit 6\n", "Radius 2.20, angular freq. 1.93\n", "Euler angles 0.19, 0.15 and -0.00\n", "Start time 0.38\n", "Rotation matrix \n", " [[ 0.988 0.001 0.151]\n", " [ 0.028 0.982 -0.189]\n", " [-0.149 0.191 0.97 ]]\n", "Date and time 2023-03-22 15:25:56.026732\n", "Time since last check is 0:01:02.238445\n" ] } ], "source": [ "import datetime\n", "now = datetime.datetime.now()\n", "print(\"Date and time\",str(now))\n", "#\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import matplotlib.animation as animation\n", "#\n", "# Can't use inline matplotlib backend for animations, so try TkAgg\n", "import matplotlib\n", "matplotlib.use('TkAgg')\n", "#\n", "# Set seed if want to have same sequence of random numbers every run\n", "#np.random.seed(137)\n", "#\n", "# Rotation matrix elements\n", "def Rx(theta):\n", " return np.matrix([[ 1, 0 , 0 ],\n", " [ 0, np.cos(theta),-np.sin(theta)],\n", " [ 0, np.sin(theta), np.cos(theta)]])\n", "#\n", "def Ry(theta):\n", " return np.matrix([[ np.cos(theta), 0, np.sin(theta)],\n", " [ 0 , 1, 0 ],\n", " [-np.sin(theta), 0, np.cos(theta)]])\n", "# \n", "def Rz(theta):\n", " return np.matrix([[ np.cos(theta), -np.sin(theta), 0 ],\n", " [ np.sin(theta), np.cos(theta) , 0 ],\n", " [ 0 , 0 , 1 ]])\n", "#\n", "# Calculate n_planets orbits\n", "n_frames = 50\n", "radii = np.array([0.6, 0.8, 1.1, 1.2, 1.5, 1.9, 2.2])\n", "#radii = np.array([0.6, 1.5])\n", "#\n", "# Make sure run for enough time to complete largest orbit (don't leave a gap!)\n", "time = (n_frames + 2)*np.amax(np.sqrt(radii**3))/n_frames\n", "n_planets = len(radii)\n", "print(\" \")\n", "print(f\"Number of planets {n_planets}, time simulated {time:.2f}\")\n", "#\n", "# Define figure\n", "fig = plt.figure(figsize = (8, 8))\n", "ax = fig.add_subplot(1, 1, 1, projection = '3d')\n", "#\n", "# Set the axis properties\n", "ax.set_xlim(-np.amax(radii), np.amax(radii))\n", "ax.set_ylim(-np.amax(radii), np.amax(radii))\n", "ax.set_zlim(-1.0, 1.0)\n", "ax.set_xlabel('x')\n", "ax.set_xlabel('y')\n", "ax.set_zlabel('z')\n", "#\n", "# Add sun to plot\n", "ax.plot(0, 0, 0, marker = '.', markersize = 20, color = 'yellow')\n", "#\n", "def update_lines_zip(num, orbits, lines):\n", " '''\n", " Update the lines in the plot, using zip to ensure that the orbits and lines are coreclty associated.\n", " '''\n", " for line, orbit in zip(lines, orbits):\n", " line.set_data(orbit[:num, :2].T)\n", " line.set_3d_properties(orbit[:num, 2])\n", " return lines\n", "#\n", "def update_lines(num, orbits, lines):\n", " '''\n", " Update the lines in the plot, using a for loop to ensure that the orbits and lines are coreclty associated.\n", " '''\n", " for n in range(0, len(lines)):\n", " lines[n].set_data(orbits[n][:num, :2].T)\n", " lines[n].set_3d_properties(orbits[n][:num, 2])\n", " return lines\n", "#\n", "def orbit(n_orbit, n_frames = 10, time = 2, radius = 1.0):\n", " '''\n", " Return a 3D \"orbit\" in an array with dimensions (n_frames, 3)\n", " '''\n", " t_step = time/n_frames\n", " period = np.sqrt(radius**3)\n", " omega = 2*np.pi/period\n", " #\n", " start = np.random.random(1)[0]\n", " alpha = 2*np.pi*(np.random.random(1)[0] - 0.5)/15.0\n", " beta = 2*np.pi*(np.random.random(1)[0] - 0.5)/15.0\n", " gamma = 2*np.pi*(np.random.random(1)[0] - 0.5)/2000.0\n", " #\n", " R = Rx(alpha)*Ry(beta)*Rz(gamma)\n", " #\n", " print(\" \")\n", " print(f\"Orbit {n_orbit}\")\n", " print(f\"Radius {radius:.2f}, angular freq. {omega:.2f}\")\n", " print(f\"Euler angles {alpha:.2f}, {beta:.2f} and {gamma:.2f}\")\n", " print(f\"Start time {start:.2f}\")\n", " print(\"Rotation matrix \\n\",np.around(R, 3))\n", " #\n", " path = np.zeros((n_frames, 3))\n", " for n in range(0, n_frames):\n", " t = n*t_step\n", " flat = np.array([radius*np.cos(omega*(t - start)), radius*np.sin(omega*(t - start)), 0.0])\n", " tilt = np.matmul(R, flat)\n", " #\n", " # matmul returns a 1D matrix \n", " path[n, :] = tilt[0, :]\n", " #\n", " return path\n", "#\n", "#orbit_test = orbit(0)\n", "#print(\" \")\n", "#print(\"orbit_test \\n\",np.around(orbit_test, 2))\n", "#\n", "# Create empty lines and make orbits \n", "lines = []\n", "orbits = []\n", "for n in range(0, n_planets):\n", " lines.append(ax.plot([], [], [])[0])\n", " orbits.append(orbit(n, n_frames, time, radii[n]))\n", "#\n", "# Make animation and save as gif. Arguments here are: \n", "# figure name, function called for each frame, number of frames, \n", "# arguments passed to function, delay between frames (in ms)\n", "ani = animation.FuncAnimation(fig, update_lines, fargs = (orbits, lines), \n", " frames = range(0, n_frames), interval = 100, blit = True)\n", "#\n", "# Save the animation (need to have imagemagick installed!)\n", "ani.save('orbit.gif', writer = 'imagemagick')\n", "#\n", "plt.show()\n", "#\n", "then = now\n", "now = datetime.datetime.now()\n", "print(\"Date and time\",str(now))\n", "print(\"Time since last check is\",str(now - then))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.16" } }, "nbformat": 4, "nbformat_minor": 4 }