Notebook

これは日々の作業を通して学んだことや毎日の生活で気づいたことをを記録しておく備忘録である。

HTML ファイル生成日時: 2024/12/05 07:26:42.662 (台灣標準時)

太陽系の立体構造

Matplotlib で三次元の図を作ることができると知ったので、太陽系内の惑星 と小惑星の軌道運動を可視化する動画を作ってみたでござる。 PNG 画像を作 る Python のプログラムは以下の通り。


#!/usr/pkg/bin/python3.9

#
# Time-stamp: <2022/10/02 23:05:00 (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_struct3'

# 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 = 5000

# number of asteroids to plot
n_asteroids = 5000

# 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'2022-07-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})...')
    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 ()
fig.subplots_adjust (left=0.0, right=1.0, bottom=0.0, top=1.0)

# making a canvas object
canvas = matplotlib.backends.backend_agg.FigureCanvasAgg (fig)

# making an axes object
ax = fig.add_subplot (111, 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 = 0.0

for i in range (n_steps):
    # clearing previous axes
    ax.cla ()
    
    # camera viewing angle
    if (i < 200):
        el = el0
        az = az0
    elif ( (i >= 200) and (i < 1400) ):
        el = el0 - (i - 200) * 0.1
        az = az0
    elif ( (i >= 1400) and (i < 1600) ):
        el = -30.0
        az = az0
    elif ( (i >= 1600) and (i < 1900) ):
        el = -30 + (i - 1600) * 0.1
        az = az0
    elif ( (i >= 1900) and (i < 2100) ):
        el = 0.0
        az = az0
    elif ( (i >= 2100) and (i < 2700) ):
        el = (i - 2100) * 0.1
        az = az0
    elif ( (i >= 2700) and (i < 3600) ):
        el = 60.0
        az = 360.0 - (i - 2700) * 0.1
    elif ( (i >= 3600) and (i < 3800) ):
        el = 60.0
        az = 270.0
    elif ( (i >= 3800) and (i < 4700) ):
        el = 60.0
        az = 270.0 - (i - 3800) * 0.1
    else:
        el = 60.0
        az = 180.0
    
    # 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_xlim (-6.0, +6.0)
    ax.set_ylim (-6.0, +6.0)
    ax.set_zlim (-2.0, +2.0)
    ax.set_box_aspect ( (6.0, 6.0, 2.0) )

    # viewing angles of camera
    ax.view_init (elev=el, azim=az)

    # using black background colour
    fig.set_facecolor ('black')
    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))

    # plotting the Sun
    sun = make_sphere (list_major[0][0][i], \
                       list_major[0][1][i], \
                       list_major[0][2][i], \
                       0.25, 'yellow')

    # 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):
        asteroid = ax.scatter (list_asteroids[j][0][i], \
                               list_asteroids[j][1][i], \
                               list_asteroids[j][2][i], \
                               s=0.1, \
                               color='saddlebrown')
    
    # title
    title = ax.text2D (0.5, 0.95, f'Inner Solar System', \
                       color='white', \
                       horizontalalignment='center', \
                       transform=ax.transAxes)

    # plotting the time
    time = ax.text2D (0.5, 0.05, f'Date/Time: {t} (UTC)', \
                      color='white', \
                      horizontalalignment='center', \
                      transform=ax.transAxes)

    # image file
    file_image = f'{file_prefix}_{i:06d}.{file_ext}'
    fig.savefig (file_image, dpi=255)

実行すると、多数の PNG ファイルができるので、 ffmpeg を使って動画にす ればよいでござる。


% ffmpeg5 -f image2 -start_number 0 -framerate 30 -i solsys_3d_struct3_%06d.png \
-an -vcodec libx264 -pix_fmt yuv420p -threads 8 solsys_3d_struct3.mp4

出来上がった動画は以下のものでござる。



Frequently accessed files

  1. Computer___Python/20220518_0.html
  2. Computer___Network/20230726_00.html
  3. Misc___Taiwan/20240207_00.html
  4. Computer___Network/20230516_00.html
  5. Computer___FreeBSD/20220621_0.html
  6. Computer___Python/20220715_0.html
  7. Computer___Network/20230508_00.html
  8. Food___Taiwan/20220429_0.html
  9. Computer___Python/20220410_0.html
  10. Computer___NetBSD/20220817_3.html
  11. Computer___Network/20240130_00.html
  12. Computer___Network/20240416_00.html
  13. Computer___Debian/20210223_1.html
  14. Computer___NetBSD/20230119_00.html
  15. Computer___Python/20221013_0.html
  16. Computer___Python/20210124_0.html
  17. Computer___NetBSD/20220428_0.html
  18. Computer___NetBSD/20220818_1.html
  19. Computer___NetBSD/20240101_02.html
  20. Science___Math/20220420_0.html
  21. Computer___NetBSD/20220808_0.html
  22. Computer___TeX/20230503_00.html
  23. Science___Astronomy/20220503_0.html
  24. Computer___NetBSD/20230515_00.html
  25. Computer___Python/20240101_00.html
  26. Computer___NetBSD/20210127_0.html
  27. Computer___Network/20220413_1.html
  28. Computer___Python/20220816_1.html
  29. Computer___NetBSD/20210204_0.html
  30. Travel___Taiwan/20220809_2.html


HTML file generated by Kinoshita Daisuke.