これは日々の作業を通して学んだことや毎日の生活で気づいたことをを記録しておく備忘録である。
HTML ファイル生成日時: 2024/11/21 17:40:55.112 (台灣標準時)
太陽、水星、金星、地球、火星、木星と小惑星二萬個の軌道運動をアニメーショ ンにしてみたでござる。
アニメーションを作るために使った Python スクリプトは以下の通りにござる。
#!/usr/pkg/bin/python3.9 # # Time-stamp: <2022/10/12 21:47:21 (CST) daisuke> # # importing sys module import sys # importing numpy module import numpy # importing astropy module import astropy import astropy.coordinates import astropy.time import astropy.units # importing astroquery module import astroquery.jplhorizons # importing matplotlib module import matplotlib.animation import matplotlib.backends.backend_agg import matplotlib.figure # output file name prefix file_prefix = 'solsys_3d_struct' # output file name extension file_ext = 'png' # units u_au = astropy.units.au u_hr = astropy.units.hour # number of steps to calculate n_steps = 10000 # number of asteroids to plot #n_asteroids = 30000 n_asteroids = 20000 #n_asteroids = 16000 #n_asteroids = 10 # step size in hr step_hr = 12 step_str = f'{step_hr}h' step = step_hr * u_hr # an empty list for storing asteroids positions list_asteroids = [] # date/time to start the simulation t_start_str = f'2020-01-01T00:00:00.000' # time to start the simulation in astropy.time object t_start = astropy.time.Time (t_start_str, format='isot', scale='utc') # time to stop the simulation in astropy.time object t_stop = t_start + step * n_steps # an empty list for storing major planets positions list_major = [] # major body names (Sun, Mercury, Venus, Earth, Mars, Jupiter) list_names = ['10', '199', '299', '399', '499', '599'] # getting positions of the Sun, Mercury, Venus, Earth, Mars, and Jupiter # from JPL/Horizons print (f'Now, getting positions of the Sun and planets...') for i in list_names: print (i) query = astroquery.jplhorizons.Horizons (id_type=None, id=f'{i}', \ location='@0', \ epochs={'start': t_start.iso, \ 'stop': t_stop.iso, \ 'step': step_str}) vec = query.vectors () print (vec) x = vec['x'] y = vec['y'] z = vec['z'] list_major.append ( [x, y, z] ) print (f'Finished getting positions of the Sun and planets!') # getting asteroids positions from JPL/Horizons print (f'Now, getting asteroids positions...') for i in range (1, n_asteroids + 1): #if (i % 10 == 0): #print (f' now, getting positions of asteroid ({i})...') print (f' now, getting positions of asteroid ({i})...') ast_query = astroquery.jplhorizons.Horizons (id_type='smallbody', \ id=f'{i}', \ location='@0', \ epochs={'start': t_start.iso, \ 'stop': t_stop.iso, \ 'step': step_str}) ast_vec = ast_query.vectors () x = ast_vec['x'] y = ast_vec['y'] z = ast_vec['z'] list_asteroids.append ( [x, y, z] ) print (f'Finished getting asteroids positions...') # making a fig object using object-oriented interface fig = matplotlib.figure.Figure (figsize=[15.36, 8.64]) fig.subplots_adjust (left=0.0, right=1.0, bottom=0.0, top=1.0, \ wspace=0.0, hspace=0.0) # making a canvas object canvas = matplotlib.backends.backend_agg.FigureCanvasAgg (fig) # making an axes object #ax = fig.add_subplot (111, projection='3d') ax = fig.add_axes ( (0, 0, 1, 1), projection='3d') # an empty list of frames for animation list_frame = [] # definition of a function for making a sphere def make_sphere (x_c, y_c, z_c, radius, colour): u = numpy.linspace (0, 2 * numpy.pi, 1000) v = numpy.linspace (0, numpy.pi, 1000) x = radius * numpy.outer (numpy.cos(u), numpy.sin(v)) + x_c y = radius * numpy.outer (numpy.sin(u), numpy.sin(v)) + y_c z = radius * numpy.outer (numpy.ones(numpy.size(u)), numpy.cos(v)) + z_c # plotting the surface sphere = ax.plot_surface (x, y, z, color=colour, antialiased=False, \ shade=True, rcount=100, ccount=100) return (sphere) # initial value of 'elev' angle el0 = 90.0 # initial value of 'azim' angle az0 = -90.0 # initial value of 'dist' dist0 = 10.0 for i in range (n_steps): # clearing previous axes ax.cla () # time t t = t_start + i * 12.0 * u_hr # printing positions of the Sun, planets, and asteroids if (i % 10 == 0): print (f'Now, making a plot for {t}...') # settings for plot ax.set_xlim3d (-7.5, +7.5) ax.set_ylim3d (-7.5, +7.5) ax.set_zlim3d (-2.0, +2.0) ax.set_box_aspect ( (7.5, 7.5, 2.0) ) # projection ax.set_proj_type ('persp') # using black background colour #fig.set_facecolor ('black') fig.set_facecolor ('white') ax.set_facecolor ('black') ax.grid (False) ax.w_xaxis.set_pane_color ((0.0, 0.0, 0.0, 0.0)) ax.w_yaxis.set_pane_color ((0.0, 0.0, 0.0, 0.0)) ax.w_zaxis.set_pane_color ((0.0, 0.0, 0.0, 0.0)) # camera viewing angle if (i < 300): el = el0 az = az0 dist = dist0 elif ( (i >= 300) and (i < 1200) ): el = el0 - (i - 300) * 0.1 az = az0 dist = dist0 elif ( (i >= 1200) and (i < 1500) ): el = 0.0 az = az0 dist = dist0 elif ( (i >= 1500) and (i < 2400) ): el = 0.0 az = az0 - (i - 1500) * 0.1 dist = dist0 elif ( (i >= 2400) and (i < 2700) ): el = 0.0 az = -180.0 dist = dist0 elif ( (i >= 2700) and (i < 3000) ): el = (i - 2700) * 0.1 az = -180.0 dist = dist0 elif ( (i >= 3000) and (i < 3300) ): el = 30.0 az = -180.0 dist = dist0 elif ( (i >= 3300) and (i < 3600) ): el = 30.0 az = -180.0 dist = dist0 - (i - 3300) * 0.01 elif ( (i >= 3600) and (i < 3900) ): el = 30.0 + (i - 3600) * 0.1 az = -180.0 dist = dist0 - 3.0 elif ( (i >= 3900) and (i < 4200) ): el = 60.0 az = -180.0 dist = dist0 - 3.0 elif ( (i >= 4200) and (i < 4500) ): el = 60.0 az = -180.0 dist = dist0 - 3.0 - (i - 4200) * 0.01 elif ( (i >= 4500) and (i < 5400) ): el = 60.0 az = -180.0 - (i - 4500) * 0.2 dist = dist0 - 6.0 elif ( (i >= 5400) and (i < 6000) ): el = 60.0 az = 0.0 dist = dist0 - 6.0 elif ( (i >= 6000) and (i < 6150) ): el = 60.0 - (i - 6000) * 0.1 az = 0.0 dist = dist0 - 6.0 elif ( (i >= 6150) and (i < 6350) ): el = 45.0 az = 0.0 dist = dist0 - 6.0 + (i - 6150) * 0.01 elif ( (i >= 6350) and (i < 6650) ): el = 45.0 az = 0.0 dist = dist0 - 4.0 elif ( (i >= 6650) and (i < 6950) ): el = 45.0 az = 0.0 - (i - 6650) * 0.3 dist = dist0 - 4.0 elif ( (i >= 6950) and (i < 7250) ): el = 45.0 az = -90.0 dist = dist0 - 4.0 elif ( (i >= 7250) and (i < 7550) ): el = 45.0 - (i - 7250) * 0.1 az = -90.0 dist = dist0 - 4.0 elif ( (i >= 7550) and (i < 7850) ): el = 15.0 az = -90.0 dist = dist0 - 4.0 elif ( (i >= 7850) and (i < 8150) ): el = 15.0 az = -90.0 dist = dist0 - 4.0 + (i - 7850) * 0.01 elif ( (i >= 8150) and (i < 8450) ): el = 15.0 az = -90.0 dist = dist0 - 1.0 elif ( (i >= 8450) and (i < 9050) ): el = 15.0 az = -90.0 dist = dist0 - 1.0 - (i - 8450) * 0.01 elif ( (i >= 9050) and (i < 9350) ): el = 15.0 az = -90.0 dist = dist0 - 7.0 elif ( (i >= 9350) and (i < 9650) ): el = 15.0 az = -90.0 dist = dist0 - 7.0 + (i - 9350) * 0.01 else: el = 15.0 az = -90.0 dist = dist0 - 4.0 # plotting the Sun sun = make_sphere (list_major[0][0][i], \ list_major[0][1][i], \ list_major[0][2][i], \ 0.25, 'orange') # plotting Mercury mercury = make_sphere (list_major[1][0][i], \ list_major[1][1][i], \ list_major[1][2][i], \ 0.05, 'cyan') # plotting Venus venus = make_sphere (list_major[2][0][i], \ list_major[2][1][i], \ list_major[2][2][i], \ 0.15, 'gold') # plotting Earth earth = make_sphere (list_major[3][0][i], \ list_major[3][1][i], \ list_major[3][2][i], \ 0.15, 'blue') # plotting Mars mars = make_sphere (list_major[4][0][i], \ list_major[4][1][i], \ list_major[4][2][i], \ 0.15, 'red') # plotting Jupiter jupiter = make_sphere (list_major[5][0][i], \ list_major[5][1][i], \ list_major[5][2][i], \ 0.15, 'bisque') # plotting asteroids for j in range (0, n_asteroids): ax.scatter (list_asteroids[j][0][i], \ list_asteroids[j][1][i], \ list_asteroids[j][2][i], \ s=1.0, \ color='saddlebrown', \ alpha=0.5) # title title = ax.text2D (0.5, 0.75, f'Inner Solar System', \ color='white', \ horizontalalignment='center', \ transform=ax.transAxes) # plotting the time time = ax.text2D (0.20, 0.23, f'Date/Time: {t} (UTC)', \ color='white', \ horizontalalignment='center', \ transform=ax.transAxes) # plotting the author name author = ax.text2D (0.85, 0.23, f'animation generated by Daisuke', \ color='white', \ horizontalalignment='center', \ transform=ax.transAxes) # viewing angles of camera ax.view_init (elev=el, azim=az) ax.dist = dist # image file file_image = f'{file_prefix}_{i:06d}.{file_ext}' fig.savefig (file_image, dpi=225)
一萬枚の画像から動画を作るためには、以下のコマンドを実行したでござる。
% ffmpeg5 -f image2 -start_number 0 -framerate 30 -i solsys_3d_struct_%06dc.png -an -vcodec libx264 -pix_fmt yuv420p -threads 12 solsys_3d_struct_no_audio.mp4
音楽をのせるには以下のようにしたでござる。
% ffmpeg5 -i solsys_3d_struct_no_audio.mp4 -i audio/fragments.f251.webm -c:v copy solsys_3d_struct_with_audio.mp4