参考

# Program Name: small_angle_and_exponent_approximation.py
# Creation Date: 20250419
# Overview: A program to numerically verify small angle and exponent approximations
# Usage: To run the program, use the command `python small_angle_and_exponent_approximation.py` in the terminal

# --- 必要なライブラリのインストール / Install necessary libraries ---
# 以下のコマンドはColabやターミナルで実行してください
# !pip install numpy matplotlib

import numpy as np
import matplotlib.pyplot as plt

# --- パラメータ設定（この部分で数値をまとめて管理）/ Define parameters in one place ---
theta_vals = np.linspace(-0.1, 0.1, 200)   # θが0に近い範囲 / θ near 0
x_vals = np.linspace(-0.1, 0.1, 200)       # xが0に近い範囲 / x near 0
y_vals = np.linspace(-0.1, 0.1, 200)       # yが0に近い範囲 / y near 0
alpha = 2                                  # 任意の実数指数 / Arbitrary exponent α
a = 3                                      # 基数 a

# --- 三角関数の近似と誤差 / Small-angle approximations and errors ---
# sinθ ≃ θ
sin_true = np.sin(theta_vals)
sin_approx = theta_vals
sin_error = np.abs(sin_true - sin_approx)

# tanθ ≃ θ
tan_true = np.tan(theta_vals)
tan_approx = theta_vals
tan_error = np.abs(tan_true - tan_approx)

# cosθ ≃ 1
cos_true = np.cos(theta_vals)
cos_approx = np.ones_like(theta_vals)
cos_error = np.abs(cos_true - cos_approx)

# --- 指数近似と誤差 / Exponential approximations and errors ---
# (1 + x)^α ≃ 1 + αx
exp1_true = (1 + x_vals)**alpha
exp1_approx = 1 + alpha * x_vals
exp1_error = np.abs(exp1_true - exp1_approx)

# (a + x)^α ≃ a^α * (1 + (x/a) * α)
exp2_true = (a + x_vals)**alpha
exp2_approx = a**alpha * (1 + (x_vals / a) * alpha)
exp2_error = np.abs(exp2_true - exp2_approx)

# (1 + x)(1 + y) ≃ 1 + x + y
exp3_true = (1 + x_vals) * (1 + y_vals)
exp3_approx = 1 + x_vals + y_vals
exp3_error = np.abs(exp3_true - exp3_approx)

# --- グラフの描画 / Plotting ---
plt.figure(figsize=(15, 10))

# --- Small-angle approximations ---
plt.subplot(2, 2, 1)
plt.plot(theta_vals, sin_error, label="|sin(θ) - θ|")
plt.plot(theta_vals, tan_error, label="|tan(θ) - θ|")
plt.plot(theta_vals, cos_error, label="|cos(θ) - 1|")
plt.title("Error of Trigonometric Approximations")
plt.xlabel("θ (radians)")
plt.ylabel("Absolute Error")
plt.legend()
plt.grid(True)

# --- Exponential approximation (1 + x)^α ---
plt.subplot(2, 2, 2)
plt.plot(x_vals, exp1_error, label="|(1+x)^α - (1+αx)|", color='orange')
plt.title("Error of (1+x)^α Approximation")
plt.xlabel("x")
plt.ylabel("Absolute Error")
plt.legend()
plt.grid(True)

# --- Exponential approximation (a + x)^α ---
plt.subplot(2, 2, 3)
plt.plot(x_vals, exp2_error, label="|(a+x)^α - a^α(1+αx/a)|", color='green')
plt.title("Error of (a+x)^α Approximation")
plt.xlabel("x")
plt.ylabel("Absolute Error")
plt.legend()
plt.grid(True)

# --- Product approximation (1+x)(1+y) ---
plt.subplot(2, 2, 4)
plt.plot(x_vals, exp3_error, label="|(1+x)(1+y) - (1+x+y)|", color='red')
plt.title("Error of (1+x)(1+y) Approximation")
plt.xlabel("x, y")
plt.ylabel("Absolute Error")
plt.legend()
plt.grid(True)

plt.tight_layout()
plt.show()

# Program Name: approximation_comparison.py
# Creation Date: 20250419
# Overview: A program to calculate and compare small-value approximations with actual values
# Usage: To run the program, use the command `python approximation_comparison.py` in the terminal

# --- 必要なライブラリのインストール / Install necessary libraries ---
# !pip install numpy matplotlib

import numpy as np
import matplotlib.pyplot as plt

# --- 数値設定（近似が有効な範囲）/ Define input range for comparison ---
x_vals = np.linspace(-0.1, 0.1, 200)  # |x| ≪ 1
θ_vals = np.linspace(-0.1, 0.1, 200)  # |θ| ≪ 1
a_vals = np.linspace(-0.1, 0.1, 200)  # |α| ≪ L
L = 1.0                               # 定数 L の値 / Value of L
n = 3                                 # 指数 n の値 / Exponent n

# --- 近似1: √(1 + x) ≃ 1 + x/2 ---
sqrt_true = np.sqrt(1 + x_vals)
sqrt_approx = 1 + x_vals / 2
sqrt_error = np.abs(sqrt_true - sqrt_approx)

# --- 近似2: (1 + x)^n ≃ 1 + nx ---
expn_true = (1 + x_vals)**n
expn_approx = 1 + n * x_vals
expn_error = np.abs(expn_true - expn_approx)

# --- 近似3: (L + α)^{1/2} ≃ L^{1/2}(1 + α/(2L)) ---
L_alpha_true = (L + a_vals)**0.5
L_alpha_approx = L**0.5 * (1 + a_vals / (2 * L))
L_alpha_error = np.abs(L_alpha_true - L_alpha_approx)

# --- 近似4: sinθ ≈ tanθ ≃ θ ---
sin_true = np.sin(θ_vals)
tan_true = np.tan(θ_vals)
sin_approx = θ_vals
tan_approx = θ_vals
sin_error = np.abs(sin_true - sin_approx)
tan_error = np.abs(tan_true - tan_approx)

# --- 近似5: cosθ ≃ 1 ---
cos_true = np.cos(θ_vals)
cos_approx = np.ones_like(θ_vals)
cos_error = np.abs(cos_true - cos_approx)

# --- グラフ表示 / Plotting ---
plt.figure(figsize=(15, 12))

# --- √(1+x) approximation ---
plt.subplot(3, 2, 1)
plt.plot(x_vals, sqrt_error, label="|√(1+x) - (1 + x/2)|", color="blue")
plt.title("Error of sqrt(1+x) Approximation")
plt.xlabel("x")
plt.ylabel("Absolute Error")
plt.grid(True)
plt.legend()

# --- (1+x)^n approximation ---
plt.subplot(3, 2, 2)
plt.plot(x_vals, expn_error, label="|(1+x)^n - (1 + nx)|", color="orange")
plt.title("Error of (1+x)^n Approximation")
plt.xlabel("x")
plt.ylabel("Absolute Error")
plt.grid(True)
plt.legend()

# --- (L+α)^{1/2} approximation ---
plt.subplot(3, 2, 3)
plt.plot(a_vals, L_alpha_error, label="|(L+α)^{1/2} - L^{1/2}(1 + α/(2L))|", color="green")
plt.title("Error of sqrt(L+α) Approximation")
plt.xlabel("α")
plt.ylabel("Absolute Error")
plt.grid(True)
plt.legend()

# --- sinθ ≃ θ and tanθ ≃ θ ---
plt.subplot(3, 2, 4)
plt.plot(θ_vals, sin_error, label="|sinθ - θ|", color="purple")
plt.plot(θ_vals, tan_error, label="|tanθ - θ|", color="red")
plt.title("Error of sinθ and tanθ Approximations")
plt.xlabel("θ (radians)")
plt.ylabel("Absolute Error")
plt.grid(True)
plt.legend()

# --- cosθ ≃ 1 ---
plt.subplot(3, 2, 5)
plt.plot(θ_vals, cos_error, label="|cosθ - 1|", color="brown")
plt.title("Error of cosθ ≃ 1 Approximation")
plt.xlabel("θ (radians)")
plt.ylabel("Absolute Error")
plt.grid(True)
plt.legend()

plt.tight_layout()
plt.show()

# Program Name: linear_approximation_difference.py
# Creation Date: 20250419
# Overview: A program to compare the exact and linear approximation of S1P - S2P using binomial expansion
# Usage: To run the program, use the command `python linear_approximation_difference.py` in the terminal

# --- 必要なライブラリのインストール / Install necessary libraries ---
# !pip install numpy matplotlib

import numpy as np
import matplotlib.pyplot as plt

# --- パラメータ定義 / Define parameters ---
L = 1.0                       # 底辺の長さ / Base length
d = 0.2                       # 中央をはさんだ距離 / Small displacement
x_vals = np.linspace(-0.1, 0.1, 200)  # xが0に近い範囲 / x close to 0

# --- 正確な式による評価 / Exact formula ---
term1 = np.sqrt(L**2 + (x_vals + d / 2)**2)
term2 = np.sqrt(L**2 + (x_vals - d / 2)**2)
exact_diff = term1 - term2

# --- 近似式 / Linear approximation using binomial expansion ---
# S1P - S2P ≒ (xd)/L
approx_diff = (x_vals * d) / L

# --- 誤差評価 / Error analysis ---
error = np.abs(exact_diff - approx_diff)

# --- グラフ描画 / Plotting ---
plt.figure(figsize=(12, 8))

# --- Exact vs Approximation ---
plt.subplot(2, 1, 1)
plt.plot(x_vals, exact_diff, label="Exact: $S_1P - S_2P$", color='blue')
plt.plot(x_vals, approx_diff, label="Approx: $\\frac{xd}{L}$", linestyle='--', color='orange')
plt.title("Exact vs Approximated $S_1P - S_2P$")
plt.xlabel("x")
plt.ylabel("Difference")
plt.legend()
plt.grid(True)

# --- Absolute error ---
plt.subplot(2, 1, 2)
plt.plot(x_vals, error, label="Absolute Error", color='red')
plt.title("Error of Approximation")
plt.xlabel("x")
plt.ylabel("Absolute Error")
plt.grid(True)
plt.legend()

plt.tight_layout()
plt.show()

# Program Name: pendulum_linear_approx_error.py
# Creation Date: 20250419
# Overview: Compare nonlinear and linear approximations of a pendulum and visualize the error
# Usage: Run with `python pendulum_linear_approx_error.py`

# --- 必要なライブラリのインストール / Install required libraries ---
# !pip install numpy matplotlib scipy

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# --- パラメータ定義 / Define physical parameters ---
g = 9.81       # 重力加速度 / gravitational acceleration [m/s^2]
L = 1.0        # 振り子の長さ / pendulum length [m]
theta0 = 0.2   # 初期角度（ラジアン）/ initial angle [rad]
omega0 = 0.0   # 初期角速度 / initial angular velocity [rad/s]
t_span = (0, 10)                          # シミュレーション時間 / simulation time
t_eval = np.linspace(*t_span, 1000)      # 時間分割 / time points

# --- 非線形モデル / Nonlinear pendulum model ---
def pendulum_nonlinear(t, y):
    theta, omega = y
    dtheta_dt = omega
    domega_dt = - (g / L) * np.sin(theta)
    return [dtheta_dt, domega_dt]

# --- 線形モデル / Linearized pendulum model ---
def pendulum_linear(t, y):
    theta, omega = y
    dtheta_dt = omega
    domega_dt = - (g / L) * theta
    return [dtheta_dt, domega_dt]

# --- 数値解法で解く / Solve both models ---
sol_nonlinear = solve_ivp(pendulum_nonlinear, t_span, [theta0, omega0], t_eval=t_eval)
sol_linear = solve_ivp(pendulum_linear, t_span, [theta0, omega0], t_eval=t_eval)

# --- 誤差の評価 / Evaluate error between the models ---
theta_error = np.abs(sol_nonlinear.y[0] - sol_linear.y[0])

# --- グラフ描画 / Plotting ---
plt.figure(figsize=(12, 8))

# θ(t)の比較
plt.subplot(2, 1, 1)
plt.plot(t_eval, sol_nonlinear.y[0], label="Nonlinear θ(t)", color='blue')
plt.plot(t_eval, sol_linear.y[0], label="Linearized θ(t)", linestyle='--', color='orange')
plt.title("Pendulum Angle θ(t): Nonlinear vs Linearized")
plt.xlabel("Time [s]")
plt.ylabel("Angle θ [rad]")
plt.legend()
plt.grid(True)

# 誤差のプロット
plt.subplot(2, 1, 2)
plt.plot(t_eval, theta_error, label="|θ_nonlinear - θ_linear|", color='red')
plt.title("Error Between Nonlinear and Linearized θ(t)")
plt.xlabel("Time [s]")
plt.ylabel("Absolute Error [rad]")
plt.legend()
plt.grid(True)

plt.tight_layout()
plt.show()

# Program Name: common_approximations_analysis.py
# Creation Date: 20250419
# Overview: Analyze common approximations (binomial, trigonometric, square root) for small x or θ
# Usage: Run this script using `python common_approximations_analysis.py`

# --- ライブラリインストール（必要に応じて）/ Install if needed ---
# !pip install numpy matplotlib

import numpy as np
import matplotlib.pyplot as plt

# --- パラメータ定義 / Define parameters ---
x_vals = np.linspace(-0.1, 0.1, 500)    # |x| ≪ 1
θ_vals = np.linspace(-0.1, 0.1, 500)    # |θ| ≪ 1
α_vals = np.linspace(-0.1, 0.1, 500)    # |α| ≪ L
L = 1.0                                 # 定数 L
n = 3                                   # 任意の自然数（べき乗の指数）

# --- 各近似式と誤差計算 / Compute approximations and errors ---

# ① 1 + x^2 ≈ 1
true1 = 1 + x_vals**2
approx1 = np.ones_like(x_vals)
error1 = np.abs(true1 - approx1)

# ② (1 + x)^n ≈ 1 + nx
true2 = (1 + x_vals)**n
approx2 = 1 + n * x_vals
error2 = np.abs(true2 - approx2)

# ③ sqrt(1 + x) ≈ 1 + x/2
true3 = np.sqrt(1 + x_vals)
approx3 = 1 + x_vals / 2
error3 = np.abs(true3 - approx3)

# ④ (L + α)^{1/2} ≈ L^{1/2}(1 + α / 2L)
true4 = (L + α_vals)**0.5
approx4 = L**0.5 * (1 + α_vals / (2 * L))
error4 = np.abs(true4 - approx4)

# ⑤ sinθ ≈ θ, tanθ ≈ θ
sin_error = np.abs(np.sin(θ_vals) - θ_vals)
tan_error = np.abs(np.tan(θ_vals) - θ_vals)

# ⑥ cosθ ≈ 1
cos_error = np.abs(np.cos(θ_vals) - 1)

# --- 可視化 / Visualization ---
plt.figure(figsize=(15, 12))

# 1 + x^2 ≈ 1
plt.subplot(3, 2, 1)
plt.plot(x_vals, error1, label="|1 + x² - 1|", color='blue')
plt.title("Error of 1 + x² ≈ 1")
plt.xlabel("x")
plt.ylabel("Absolute Error")
plt.grid(True)
plt.legend()

# (1 + x)^n ≈ 1 + nx
plt.subplot(3, 2, 2)
plt.plot(x_vals, error2, label="|(1+x)^n - (1+nx)|", color='orange')
plt.title("Error of (1 + x)^n ≈ 1 + nx")
plt.xlabel("x")
plt.ylabel("Absolute Error")
plt.grid(True)
plt.legend()

# sqrt(1 + x) ≈ 1 + x/2
plt.subplot(3, 2, 3)
plt.plot(x_vals, error3, label="|√(1+x) - (1 + x/2)|", color='green')
plt.title("Error of √(1 + x) ≈ 1 + x/2")
plt.xlabel("x")
plt.ylabel("Absolute Error")
plt.grid(True)
plt.legend()

# (L + α)^{1/2} ≈ L^{1/2}(1 + α / 2L)
plt.subplot(3, 2, 4)
plt.plot(α_vals, error4, label="|(L+α)^{1/2} - L^{1/2}(1 + α/(2L))|", color='purple')
plt.title("Error of (L+α)^{1/2} ≈ L^{1/2}(1 + α/(2L))")
plt.xlabel("α")
plt.ylabel("Absolute Error")
plt.grid(True)
plt.legend()

# sinθ ≈ θ, tanθ ≈ θ
plt.subplot(3, 2, 5)
plt.plot(θ_vals, sin_error, label="|sinθ - θ|", color='red')
plt.plot(θ_vals, tan_error, label="|tanθ - θ|", linestyle='--', color='orange')
plt.title("Error of sinθ ≈ θ and tanθ ≈ θ")
plt.xlabel("θ [rad]")
plt.ylabel("Absolute Error")
plt.grid(True)
plt.legend()

# cosθ ≈ 1
plt.subplot(3, 2, 6)
plt.plot(θ_vals, cos_error, label="|cosθ - 1|", color='brown')
plt.title("Error of cosθ ≈ 1")
plt.xlabel("θ [rad]")
plt.ylabel("Absolute Error")
plt.grid(True)
plt.legend()

plt.tight_layout()
plt.show()

# Program Name: infinitesimal_analysis_mechanics.py
# Creation Date: 20250419
# Overview: Analyze infinitesimal quantities and their approximations in spring and gravitational energy
# Usage: Run with `python infinitesimal_analysis_mechanics.py`

# --- ライブラリインストール（必要なら）/ Install libraries if needed ---
# !pip install numpy matplotlib

import numpy as np
import matplotlib.pyplot as plt

# --- 定数の設定 / Define physical constants ---
k = 2.0          # ばね定数 / spring constant [N/m]
G = 6.67e-11     # 万有引力定数 / gravitational constant
M = 5.97e24      # 質量M（地球）/ Mass of large object (e.g., Earth) [kg]
m = 1.0          # 質量m / small object [kg]

# --- ばねの仕事の微小量導出 / Spring energy via infinitesimal work ---
X_vals = np.linspace(0, 1, 1000)  # x from 0 to 1 meter
F_spring = k * X_vals             # F = kx
dW_spring = F_spring * (X_vals[1] - X_vals[0])  # 微小区間の仕事 dW
W_integral = 0.5 * k * X_vals**2  # 積分による理論解

# --- 万有引力による仕事 / Gravitational potential energy ---
r_vals = np.linspace(6.4e6, 1e7, 1000)  # 地表～高空 / radius from Earth center [m]
F_gravity = G * M * m / r_vals**2       # F = GMm/r²
dr = r_vals[1] - r_vals[0]
dW_gravity = -F_gravity * dr            # 微小仕事（負）
U_grav = -G * M * m / r_vals            # 理論位置エネルギー

# --- 微小量近似式と誤差比較 / Approximation formulas ---
x = np.linspace(-0.1, 0.1, 500)
alpha = x
theta = x
L = 1.0
n = 2

# 近似式一覧と誤差
sqrt_true = np.sqrt(1 + x)
sqrt_approx = 1 + x / 2
sqrt_error = np.abs(sqrt_true - sqrt_approx)

binomial_true = (1 + x)**n
binomial_approx = 1 + n * x
binomial_error = np.abs(binomial_true - binomial_approx)

trig_sin_error = np.abs(np.sin(theta) - theta)
trig_cos_error = np.abs(np.cos(theta) - 1)
trig_tan_error = np.abs(np.tan(theta) - theta)

# --- 描画 / Plotting ---
plt.figure(figsize=(15, 14))

# ばねのエネルギー
plt.subplot(3, 2, 1)
plt.plot(X_vals, W_integral, label="W = (1/2)kx²", color='blue')
plt.title("Spring Work Accumulated Over Distance")
plt.xlabel("Extension x [m]")
plt.ylabel("Elastic Potential Energy [J]")
plt.grid(True)
plt.legend()

# 万有引力のエネルギー
plt.subplot(3, 2, 2)
plt.plot(r_vals / 1e6, U_grav, label="U = -GMm/r", color='red')
plt.title("Gravitational Potential Energy")
plt.xlabel("Distance r [10⁶ m]")
plt.ylabel("Potential Energy [J]")
plt.grid(True)
plt.legend()

# √(1+x) 近似
plt.subplot(3, 2, 3)
plt.plot(x, sqrt_error, label="|√(1+x) - (1 + x/2)|", color='green')
plt.title("Error of sqrt(1 + x) ≈ 1 + x/2")
plt.xlabel("x")
plt.ylabel("Absolute Error")
plt.grid(True)
plt.legend()

# (1+x)^n 近似
plt.subplot(3, 2, 4)
plt.plot(x, binomial_error, label="|(1+x)^n - (1 + nx)|", color='orange')
plt.title("Error of (1 + x)^n ≈ 1 + nx")
plt.xlabel("x")
plt.ylabel("Absolute Error")
plt.grid(True)
plt.legend()

# 三角関数の小角近似
plt.subplot(3, 2, 5)
plt.plot(theta, trig_sin_error, label="|sinθ - θ|", color='purple')
plt.plot(theta, trig_tan_error, label="|tanθ - θ|", color='brown')
plt.title("Error of sinθ ≈ θ, tanθ ≈ θ")
plt.xlabel("θ [rad]")
plt.ylabel("Absolute Error")
plt.grid(True)
plt.legend()

# cosθ ≈ 1
plt.subplot(3, 2, 6)
plt.plot(theta, trig_cos_error, label="|cosθ - 1|", color='gray')
plt.title("Error of cosθ ≈ 1")
plt.xlabel("θ [rad]")
plt.ylabel("Absolute Error")
plt.grid(True)
plt.legend()

plt.tight_layout()
plt.show()

# Program Name: youngs_interference_pattern.py
# Creation Date: 20250419
# Overview: Simulates interference pattern from Young's double slit experiment and validates approximation formulas
# Usage: Run with `python youngs_interference_pattern.py`

# --- 必要なライブラリ / Required Libraries ---
# !pip install numpy matplotlib

import numpy as np
import matplotlib.pyplot as plt

# --- 定数設定 / Constants ---
λ = 600e-9            # 波長 [m] / wavelength (e.g., red light)
L = 1.0               # スリットからスクリーンまでの距離 [m]
d_values = [0.5e-3, 1e-3, 2e-3]  # スリット間隔 [m] 変化させて比較
m = np.arange(-10, 11)          # 明線の次数

# --- 干渉縞の間隔Δyを計算 / Compute fringe positions ---
plt.figure(figsize=(12, 6))
for d in d_values:
    y = m * λ * L / d  # Δy = mλL/d
    plt.plot(m, y * 1e3, 'o-', label=f"d = {d*1e3:.1f} mm")

plt.title("Interference Fringe Positions (Young's Experiment)")
plt.xlabel("Order m")
plt.ylabel("Fringe Position y [mm]")
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.show()

# --- 光路差 Δ = d sinθ の近似 / Path difference approximation ---
θ_vals = np.linspace(-0.01, 0.01, 500)  # [rad], small angle
d = 1e-3  # fix for now
delta_exact = d * np.sin(θ_vals)
delta_approx = d * θ_vals  # sinθ ≈ θ
error = np.abs(delta_exact - delta_approx)

# --- 誤差プロット / Error Plot ---
plt.figure(figsize=(12, 5))

plt.subplot(1, 2, 1)
plt.plot(θ_vals, delta_exact * 1e9, label="Exact Δ", color='blue')
plt.plot(θ_vals, delta_approx * 1e9, '--', label="Approx Δ ≈ dθ", color='orange')
plt.title("Path Difference Δ vs Angle")
plt.xlabel("θ [rad]")
plt.ylabel("Δ [nm]")
plt.grid(True)
plt.legend()

plt.subplot(1, 2, 2)
plt.plot(θ_vals, error * 1e9, color='red')
plt.title("Error of sinθ ≈ θ Approximation")
plt.xlabel("θ [rad]")
plt.ylabel("Absolute Error [nm]")
plt.grid(True)

plt.tight_layout()
plt.show()

# Program Name: linear_approximations_full.py
# Creation Date: 20250419
# Overview: Visualize true values, approximations, and errors for common linear approximations in high school physics
# Usage: Run this script with `python linear_approximations_full.py`

import numpy as np
import matplotlib.pyplot as plt

# --- 定義域とパラメータ設定 / Input ranges and constants ---
x = np.linspace(-0.1, 0.1, 500)
theta = np.linspace(-0.2, 0.2, 500)
alpha = np.linspace(-0.1, 0.1, 500)
L = 1.0
n = 2

# --- 各近似式とその比較 / Define and compute approximations ---

# ① sinθ ≈ θ
sin_true = np.sin(theta)
sin_approx = theta
sin_error = np.abs(sin_true - sin_approx)

# ② tanθ ≈ θ
tan_true = np.tan(theta)
tan_approx = theta
tan_error = np.abs(tan_true - tan_approx)

# ③ cosθ ≈ 1
cos_true = np.cos(theta)
cos_approx = np.ones_like(theta)
cos_error = np.abs(cos_true - cos_approx)

# ④ sqrt(1 + x) ≈ 1 + x/2
sqrt_true = np.sqrt(1 + x)
sqrt_approx = 1 + x / 2
sqrt_error = np.abs(sqrt_true - sqrt_approx)

# ⑤ (1 + x)^n ≈ 1 + nx
binom_true = (1 + x)**n
binom_approx = 1 + n * x
binom_error = np.abs(binom_true - binom_approx)

# ⑥ 1/(1 + x) ≈ 1 - x
inv_true = 1 / (1 + x)
inv_approx = 1 - x
inv_error = np.abs(inv_true - inv_approx)

# ⑦ sqrt(L + α) ≈ sqrt(L)(1 + α/(2L))
L_sqrt_true = np.sqrt(L + alpha)
L_sqrt_approx = np.sqrt(L) * (1 + alpha / (2 * L))
L_sqrt_error = np.abs(L_sqrt_true - L_sqrt_approx)

# --- 描画 / Plotting ---
plt.figure(figsize=(18, 20))

# ① sinθ
plt.subplot(4, 2, 1)
plt.plot(theta, sin_true, label="True: sinθ", color='blue')
plt.plot(theta, sin_approx, '--', label="Approx: θ", color='orange')
plt.plot(theta, sin_error, ':', label="Error |sinθ - θ|", color='red')
plt.title("① sinθ ≈ θ")
plt.xlabel("θ [rad]")
plt.ylabel("Value")
plt.grid(True)
plt.legend()

# ② tanθ
plt.subplot(4, 2, 2)
plt.plot(theta, tan_true, label="True: tanθ", color='blue')
plt.plot(theta, tan_approx, '--', label="Approx: θ", color='orange')
plt.plot(theta, tan_error, ':', label="Error |tanθ - θ|", color='red')
plt.title("② tanθ ≈ θ")
plt.xlabel("θ [rad]")
plt.ylabel("Value")
plt.grid(True)
plt.legend()

# ③ cosθ
plt.subplot(4, 2, 3)
plt.plot(theta, cos_true, label="True: cosθ", color='blue')
plt.plot(theta, cos_approx, '--', label="Approx: 1", color='orange')
plt.plot(theta, cos_error, ':', label="Error |cosθ - 1|", color='red')
plt.title("③ cosθ ≈ 1")
plt.xlabel("θ [rad]")
plt.ylabel("Value")
plt.grid(True)
plt.legend()

# ④ sqrt(1+x)
plt.subplot(4, 2, 4)
plt.plot(x, sqrt_true, label="True: √(1 + x)", color='blue')
plt.plot(x, sqrt_approx, '--', label="Approx: 1 + x/2", color='orange')
plt.plot(x, sqrt_error, ':', label="Error", color='red')
plt.title("④ √(1 + x) ≈ 1 + x/2")
plt.xlabel("x")
plt.ylabel("Value")
plt.grid(True)
plt.legend()

# ⑤ (1 + x)^n
plt.subplot(4, 2, 5)
plt.plot(x, binom_true, label="True: (1 + x)^n", color='blue')
plt.plot(x, binom_approx, '--', label="Approx: 1 + nx", color='orange')
plt.plot(x, binom_error, ':', label="Error", color='red')
plt.title("⑤ (1 + x)^n ≈ 1 + nx")
plt.xlabel("x")
plt.ylabel("Value")
plt.grid(True)
plt.legend()

# ⑥ 1 / (1 + x)
plt.subplot(4, 2, 6)
plt.plot(x, inv_true, label="True: 1 / (1 + x)", color='blue')
plt.plot(x, inv_approx, '--', label="Approx: 1 - x", color='orange')
plt.plot(x, inv_error, ':', label="Error", color='red')
plt.title("⑥ 1 / (1 + x) ≈ 1 - x")
plt.xlabel("x")
plt.ylabel("Value")
plt.grid(True)
plt.legend()

# ⑦ sqrt(L + α)
plt.subplot(4, 2, 7)
plt.plot(alpha, L_sqrt_true, label="True: √(L + α)", color='blue')
plt.plot(alpha, L_sqrt_approx, '--', label="Approx: √L(1 + α/(2L))", color='orange')
plt.plot(alpha, L_sqrt_error, ':', label="Error", color='red')
plt.title("⑦ √(L + α) ≈ √L(1 + α/(2L))")
plt.xlabel("α")
plt.ylabel("Value")
plt.grid(True)
plt.legend()

plt.tight_layout()
plt.show()

# Program Name: newtons_rings_path_difference.py
# Creation Date: 20250419
# Overview: Calculate and visualize the air gap thickness in Newton's rings using exact and approximate formulas
# Usage: Run with `python newtons_rings_path_difference.py`

import numpy as np
import matplotlib.pyplot as plt

# --- パラメータ定義 / Define parameters ---
R = 1.0         # レンズの曲率半径 [m]
r_vals = np.linspace(0, 0.01, 500)  # 中心からの距離 r [m]（最大1cm）

# --- (1) 正確な式 / Exact thickness: R² = (R - t)² + r² → t = R - √(R² - r²)
t_exact = R - np.sqrt(R**2 - r_vals**2)

# --- (2) 近似式 / Approximate thickness: t ≈ r² / (2R)
t_approx = r_vals**2 / (2 * R)

# --- 誤差評価 / Absolute error between exact and approximate
error = np.abs(t_exact - t_approx)

# --- プロット / Plotting results ---
plt.figure(figsize=(14, 6))

# 厚みの比較 / Thickness Comparison
plt.subplot(1, 2, 1)
plt.plot(r_vals * 1000, t_exact * 1e6, label="Exact: $t = R - \sqrt{R^2 - r^2}$", color='blue')
plt.plot(r_vals * 1000, t_approx * 1e6, '--', label="Approx: $t ≈ \\frac{r^2}{2R}$", color='orange')
plt.title("Newton's Rings: Air Gap Thickness")
plt.xlabel("Distance from Center r [mm]")
plt.ylabel("Air Gap Thickness t [μm]")
plt.legend()
plt.grid(True)

# 誤差プロット / Error
plt.subplot(1, 2, 2)
plt.plot(r_vals * 1000, error * 1e9, label="Absolute Error [nm]", color='red')
plt.title("Error of Approximation")
plt.xlabel("Distance from Center r [mm]")
plt.ylabel("Error [nm]")
plt.grid(True)
plt.legend()

plt.tight_layout()
plt.show()