import numpy as np
import matplotlib.pyplot as plt
from math import pi, cos, sin, sqrt

class BallisticSimulator:
    def __init__(self, dt=0.01):
        self.dt = dt
        self.g = 9.80665

    def get_accelerations(self, vx, vy, k, am):
        """
        現在の速度と設定から、x方向とy方向の加速度を返す
        k: 空気抵抗係数, am: 磁気浮上加速度
        """
        v = sqrt(vx**2 + vy**2)
        # x方向: 空気抵抗のみ
        ax = -k * vx * v
        # y方向: 重力(下) + 磁気浮上(上) + 空気抵抗(上下)
        ay = -self.g + am - (k * vy * v)
        return ax, ay

    def simulate(self, mode="normal", v0_kmh=1000.0, deg=0.0, ht=1.5, mas_g=9.0, area_cm2=0.5):
        # モードごとの環境設定
        config = {
            "normal":      {"rho": 1.225, "am": 0.0},      # 1. 通常（空気あり、浮力なし）
            "vacuum":      {"rho": 0.0,   "am": 0.0},      # 2. 真空（空気なし、浮力なし）
            "maglev":      {"rho": 0.0,   "am": 9.80665}   # 3. 真空 + 磁気浮上
        }
        env = config.get(mode, config["normal"])
        v0 = v0_kmh / 3.6
        rad = deg * pi / 180
        # 物理定数の計算
        mas = mas_g / 1000.0
        s_area = area_cm2 / 10000.0
        cd = 0.3
        k = 0.5 * s_area * cd * env["rho"] / mas
        # 初期値
        x, y = 0.0, ht
        vx, vy = v0 * cos(rad), v0 * sin(rad)
        res_x, res_y = [x], [y]
        for _ in range(1000000):
            k1vx, k1vy = self.get_accelerations(vx, vy, k, env["am"])
            k2vx, k2vy = self.get_accelerations(vx + k1vx*self.dt/2, vy + k1vy*self.dt/2, k, env["am"])
            k3vx, k3vy = self.get_accelerations(vx + k2vx*self.dt/2, vy + k2vy*self.dt/2, k, env["am"])
            k4vx, k4vy = self.get_accelerations(vx + k3vx*self.dt, vy + k3vy*self.dt, k, env["am"])
            vx += (k1vx + 2*k2vx + 2*k3vx + k4vx) / 6 * self.dt
            vy += (k1vy + 2*k2vy + 2*k3vy + k4vy) / 6 * self.dt
            x += vx * self.dt
            y += vy * self.dt
            if y < 0 or x >= 300000: break
            res_x.append(x)
            res_y.append(y)
        return res_x, res_y

# --- 実行と可視化 ---
sim = BallisticSimulator()
modes = ["normal", "vacuum", "maglev"]
mode = modes[1]

rx, ry = sim.simulate(mode=mode, v0_kmh=1000.0, deg=0)
print(f"飛行時間  : {(len(rx) - 1) * sim.dt:.2f} 秒")
print(f"移動距離  : {rx[-1]:.2f} メートル, {rx[-1]/1000:.5f} キロメートル")
plt.figure(figsize=(16, 3))
plt.plot(rx, ry, label=f"Mode: {mode}")
plt.axhline(0, color='black', linewidth=1)
plt.title("Comparison of Bullet Trajectories (v=1000km/h)")
plt.xlabel("Distance (m)")
plt.ylabel("Height (m)")
plt.legend()
plt.grid(True)
plt.show()

タグ ：

#プログラミング

#Python

#名探偵コナン

#緋色の弾丸

#弾道シミュレーション

brew install ffmpeg py3cairo pkg-config

brew install --cask mactex-no-gui

pip install manim

manim --version

manim -pqh ファイル名 シーン名

# hello.py
from manim import *

class Hello(Scene):
    def construct(self):
        text = Text("どうもこんばんわ、乙pyです")
        self.play(Write(text))
        self.wait()

manim -pqh hello.py Hello

from manim import *

class MyAnimation(Scene):
    def construct(self):

変数名 = オブジェクト名(引数など)

self.play(): 
アニメーションを実行

.animate: 
オブジェクトのプロパティ（位置や色）を滑らかに変化

ReplacementTransform: 
前の物体を消して、次の物体へモーフィング

# basicflow.py
# manim -pqh basicflow.py BasicFlow
from manim import *

class BasicFlow(Scene):
    def construct(self):
        # 1. オブジェクトの生成
        title = Text("Manimの基本操作", font_size=40)
        square = Square(color=BLUE, fill_opacity=0.5)
        circle = Circle(color=RED, fill_opacity=0.5)

        # 2. 配置の設定
        title.to_edge(UP) # 画面上端へ

        # 3. アニメーションの実行
        self.play(Write(title))        # 文字を書く
        self.play(DrawBorderThenFill(square)) # 枠線から塗りつぶし
        self.wait(1)

        # 4. 移動と変形 (Transform)
        self.play(
            square.animate.shift(LEFT * 2), # 左に動かしながら
            run_time=1
        )
        # 四角を円に変形させる
        self.play(ReplacementTransform(square, circle))
        self.play(circle.animate.shift(RIGHT * 4))
        
        self.wait(2)

タグ ：

#Python

#プログラミング

#Manim

#数式

#アニメーション

地震の発生確率の計算方法

prob = (f_future - f_t) / (1 - f_t)

import numpy as np
from scipy.stats import norm

def bpt_cdf_30(t, mu, alpha, delta_t=30):
    """
    BPT分布に基づき、今後 delta_t 年以内の地震発生確率を計算する
    
    Parameters:
    t: 前回からの経過年数
    mu: 平均再来間隔
    alpha: ばらつき係数
    delta_t: 評価期間（デフォルトは30年）
    """
    if t <= 0:
        return 0.0
    
    # 現在までの間に、既に地震が起きているはずだった累積確率
    inv_alpha_sq = 2 / (alpha**2)
    term1_t = norm.cdf((t / mu - 1) / (alpha * np.sqrt(t / mu)))
    term2_t = np.exp(inv_alpha_sq) * norm.cdf(-(t / mu + 1) / (alpha * np.sqrt(t / mu)))
    f_t = term1_t + term2_t
    
    # delta_t 年後までの間に、地震が起きる累積確率
    t_future = t + delta_t
    term1_f = norm.cdf((t_future / mu - 1) / (alpha * np.sqrt(t_future / mu)))
    term2_f = np.exp(inv_alpha_sq) * norm.cdf(-(t_future / mu + 1) / (alpha * np.sqrt(t_future / mu)))
    f_future = term1_f + term2_f
    
    if f_t >= 1.0:
        return 1.0
        
    prob = (f_future - f_t) / (1 - f_t)
    return max(0.0, min(prob, 1.0))

mu = 97.5       # 平均再来間隔 (南海トラフの直近データ 1361 ~ 1946)
alpha = 0.24    # ばらつき係数 (地震調査委員会採用値)
last_event = 1946  # 前回の発生年 (昭和南海地震)
target_year = 2026
t = target_year - last_event

# 30年分を予測
for delta_t in range(1, 31):
  prob = bpt_cdf_30(t, mu, alpha, delta_t)
  print(f"今後 {delta_t} 年以内の発生確率: {prob:.2%}")

タグ ：

#Python

#プログラミング

#南海トラフ地震

#シミュレーション

#BPT分布