Notebook

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

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

ルンゲクッタによる軌道積分

ルンゲクッタを使って、軌道積分をしてみたでござる。コードは以下の通りに ござる。 scipy.integrate に含まれている solve_ivp () という関数を使っ たでござる。


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

#
# Time-stamp: <2022/05/01 16:14:03 (CST) daisuke>
#

# importing argparse module
import argparse

# importing sys module
import sys

# importing numpy module
import numpy

# importing scipy module
import scipy.integrate

# importing matplotlib module
import matplotlib.figure
import matplotlib.backends.backend_agg

# initialising a parser
desc   = 'Keplerian motion'
parser = argparse.ArgumentParser (description=desc)

# adding arguments
parser.add_argument ('-x0', '--x0', type=float, default=1.0,
                     help='initial position in X (default: 1 au)')
parser.add_argument ('-y0', '--y0', type=float, default=0.0,
                     help='initial position in Y (default: 0 au)')
parser.add_argument ('-vx0', '--vx0', type=float, default=0.0,
                     help='initial velocity in X (default: 0 sqrt(GM)')
parser.add_argument ('-vy0', '--vy0', type=float, default=1.0,
                     help='initial velocity in Y (default: 1 sqrt(GM))')
parser.add_argument ('-M', '--M', type=float, default=1.0,
                     help='mass of star (default: 1 solar mass)')
parser.add_argument ('-t', '--time', type=float, default=100.0,
                     help='time of simulation (default: 100 yr)')
parser.add_argument ('-i', '--interval', type=float, default=0.01,
                     help='time interval of image creation (default: 0.01 yr)')
parser.add_argument ('-r', '--resolution', type=int, default=225, \
                     help='resolution of output file (default: 225 dpi)')
parser.add_argument ('-o', '--output', default='', \
                     help='output PNG file prefix')

# parsing arguments
args = parser.parse_args ()

# parameters
qx0           = args.x0
qy0           = args.y0
vx0           = args.vx0
vy0           = args.vy0
Mstar         = args.M
time_end      = args.time
dt            = args.interval
resolution    = args.resolution
output_prefix = args.output

# check of output_prefix
if (output_prefix == ''):
    print ("ERROR: output file prefix must be given!")
    sys.exit ()

#
# constants
#

# unit of time: year
# unit of distance: au

# gravitational constant
GM = 4.0 * numpy.pi * numpy.pi

#
# equation of motion
#
def eqmo (t, y):
    dy = numpy.zeros_like (y)

    r_cubed = ( y[0]**2 + y[2]**2 )**1.5

    dy[0] = y[1]
    dy[1] = -GM * y[0] / r_cubed
    dy[2] = y[3]
    dy[3] = -GM * y[2] / r_cubed

    return dy

# time to write position and velocity
n_step = int (time_end / dt) + 1
t_eval = numpy.linspace (0.0, time_end, n_step)

# initial values
y_init = (qx0, vx0  * numpy.sqrt (GM), qy0, vy0 * numpy.sqrt (GM))

# orbital integration
sol = scipy.integrate.solve_ivp (eqmo, [0.0, time_end], y_init, \
                                 t_eval=t_eval, dense_output=True, \
                                 rtol=10**-6, atol=10**-9)

# results (positions and velocities)
qx = sol.y[0]
qy = sol.y[2]
vx = sol.y[1]
vy = sol.y[3]

# finding maximum and minimum values
qx_min = +10**10
qx_max = -10**10
qy_min = +10**10
qy_max = -10**10
for i in range ( len (qx) ):
    if (qx[i] < qx_min):
        qx_min = qx[i]
    if (qx[i] > qx_max):
        qx_max = qx[i]
    if (qy[i] < qy_min):
        qy_min = qy[i]
    if (qy[i] > qy_max):
        qy_max = qy[i]
list_maxmin = [abs (qx_min), abs (qx_max), abs (qy_min), abs (qy_max)]
sorted_maxmin = sorted (list_maxmin)
x_min = -1.0 * sorted_maxmin[-1] * 1.2
x_max = +1.0 * sorted_maxmin[-1] * 1.2
y_min = -1.0 * sorted_maxmin[-1] * 1.2
y_max = +1.0 * sorted_maxmin[-1] * 1.2

# making plots
for i in range ( len (qx) ):
    print ("%08d %8.4f %8.4f %8.4f %8.4f" % (i, qx[i], qy[i], vx[i], vy[i]) )

    # output file name
    file_output = "%s_%08d.png" % (output_prefix, i)
    
    #
    # plotting using Matplotlib
    #
    
    # making objects "fig" and "ax"
    fig = matplotlib.figure.Figure ()
    matplotlib.backends.backend_agg.FigureCanvasAgg (fig)
    ax = fig.add_subplot (111)

    # axes
    ax.set_title ('Planetary Motion')
    ax.set_xlabel ('X [au]')
    ax.set_ylabel ('Y [au]')
    ax.set_xlim (x_min, x_max)
    ax.set_ylim (y_min, y_max)
    ax.set_aspect ('equal')

    # plotting a figure
    ax.plot (0.0, 0.0, color='yellow', marker='o', markersize=10, \
             label='star')
    if (i < 100):
        for j in range (i):
            ax.plot (qx[j], qy[j], color='cyan', marker='o', markersize=3, \
                     alpha=j/i)
    else:
        for j in range (100):
            ax.plot (qx[i-j], qy[i-j], color='cyan', marker='o', markersize=3, \
                     alpha=1.0-j/100.0)
    ax.plot (qx[0:i], qy[0:i], color='black', marker=',', alpha=0.5)
    ax.plot (qx[i], qy[i], color='green', marker='o', markersize=5, \
             label='planet')
    text_time = "Time: %8.2f year" % (i * dt)
    text_initial = "Initial conditions"
    text_mass    = "mass of star = %5.2f solar mass" % (Mstar)
    text_iq      = "(qx0, qy0) = (%4.2f au, %4.2f au)" \
        % (qx0, qy0)
    text_iv      = "(vx0, vy0) = (%4.2f au/yr, %4.2f au/yr)" \
        % (vx0 * numpy.sqrt (GM), vy0 * numpy.sqrt (GM))
    ax.text (0.03, 0.95, text_time, transform=ax.transAxes)
    ax.text (0.03, 0.18, text_initial, transform=ax.transAxes)
    ax.text (0.05, 0.13, text_mass, transform=ax.transAxes)
    ax.text (0.05, 0.08, text_iq, transform=ax.transAxes)
    ax.text (0.05, 0.03, text_iv, transform=ax.transAxes)
    ax.legend ()

    # saving the figure to a file
    fig.savefig (file_output, dpi=resolution)

例えば、以下のように実行すればよいでござる。


% chmod a+x kepler_00.py
% ./kepler_00.py -h
usage: kepler_00.py [-h] [-x0 X0] [-y0 Y0] [-vx0 VX0] [-vy0 VY0] [-M M]
                    [-t TIME] [-i INTERVAL] [-r RESOLUTION] [-o OUTPUT]

Keplerian motion

optional arguments:
  -h, --help            show this help message and exit
  -x0 X0, --x0 X0       initial position in X (default: 1 au)
  -y0 Y0, --y0 Y0       initial position in Y (default: 0 au)
  -vx0 VX0, --vx0 VX0   initial velocity in X (default: 0 sqrt(GM)
  -vy0 VY0, --vy0 VY0   initial velocity in Y (default: 1 sqrt(GM))
  -M M, --M M           mass of star (default: 1 solar mass)
  -t TIME, --time TIME  time of simulation (default: 100 yr)
  -i INTERVAL, --interval INTERVAL
                        time interval of image creation (default: 0.01 yr)
  -r RESOLUTION, --resolution RESOLUTION
                        resolution of output file (default: 225 dpi)
  -o OUTPUT, --output OUTPUT
                        output PNG file prefix

% ./kepler_00.py -o orb -vy0 1.3

多数の PNG ファイルができるので、それらを使って動画を作れば完成にござる。


% ffmpeg5 -f image2 -start_number 0 -framerate 30 -i orb_%08d.png \
? -an -vcodec libx264 -pix_fmt yuv420p -threads 12 orb_no_audio.mp4
% ffmpeg5 -i orb_no_audio.mp4 -i foo.f251.webm -c:v copy orb_with_audio.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.