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

# Marshall-Palmer分布のパラメータ
N0 = 8000.0  # [m^-3 mm^-1]

def marshall_palmer_N(D, R_rain):
    """
    Marshall-Palmer粒径分布
    D: 雨滴直径 [mm]
    R_rain: 降水強度 [mm/h]
    """
    Lambda = 4.1 * R_rain**(-0.21)  # [mm^-1]
    return N0 * np.exp(-Lambda * D)

# 異なる降水強度での粒径分布
R_rains = [1, 5, 10, 25, 50]  # [mm/h]
D = np.linspace(0, 6, 500)     # 雨滴直径 [mm]

fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# 粒径分布
for R_r in R_rains:
    N = marshall_palmer_N(D, R_r)
    axes[0].semilogy(D, N, linewidth=2, label=f'R = {R_r} mm/h')

axes[0].set_xlabel('Drop Diameter D [mm]', fontsize=13)
axes[0].set_ylabel(r'$N(D)$ [$m^{-3} mm^{-1}$]', fontsize=13)
axes[0].set_title('Marshall-Palmer Drop Size Distribution', fontsize=13)
axes[0].set_xlim([0, 6])
axes[0].set_ylim([1, 1e4])
axes[0].legend(fontsize=11)
axes[0].grid(True, alpha=0.3)

# Z の計算（数値積分）
R_rain_range = np.logspace(-1, 2, 200)  # 0.1 ~ 100 mm/h
Z_values = np.zeros_like(R_rain_range)

for i, R_r in enumerate(R_rain_range):
    Lambda = 4.1 * R_r**(-0.21)
    # Z = ∫ N(D) D^6 dD = N0 * 720 / Lambda^7
    Z_values[i] = N0 * 720.0 / Lambda**7  # [mm^6/m^3]

# Z-R関係のフィット（Marshall-Palmer）
Z_mp = 200 * R_rain_range**1.6

# dBZ変換
Z_dBZ = 10 * np.log10(Z_values + 1e-30)
Z_mp_dBZ = 10 * np.log10(Z_mp + 1e-30)

axes[1].plot(R_rain_range, Z_dBZ, 'b-', linewidth=2,
             label='Numerical integration')
axes[1].plot(R_rain_range, Z_mp_dBZ, 'r--', linewidth=2,
             label=r'$Z = 200 R^{1.6}$')
axes[1].set_xlabel('Rainfall Rate R [mm/h]', fontsize=13)
axes[1].set_ylabel('Z [dBZ]', fontsize=13)
axes[1].set_title('Z-R Relationship', fontsize=13)
axes[1].set_xscale('log')
axes[1].legend(fontsize=12)
axes[1].grid(True, alpha=0.3, which='both')

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

def z_to_rain(Z_dBZ, a=200, b=1.6):
    """
    Z-R関係の逆変換: Z [dBZ] → R [mm/h]
    Z = a * R^b → R = (Z/a)^(1/b)
    """
    Z_linear = 10**(Z_dBZ / 10)  # dBZ → 線形値
    R = (Z_linear / a)**(1.0 / b)
    return R

# 仮想的なレーダー観測データ（PPIスキャン）
np.random.seed(42)
n_azimuth = 360   # 方位角の分割数
n_range = 200     # 距離ゲートの数
max_range = 100   # 最大距離 [km]

azimuth = np.linspace(0, 2 * np.pi, n_azimuth, endpoint=False)
range_gate = np.linspace(1, max_range, n_range)
R_grid, Az_grid = np.meshgrid(range_gate, azimuth)
X = R_grid * np.cos(Az_grid)
Y = R_grid * np.sin(Az_grid)

# 降水エコーの生成（2つのセルを重ね合わせ）
# セル1: ガウス状の降水域
x1, y1 = 30, 20
sigma1 = 15
cell1 = 45 * np.exp(-((X - x1)**2 + (Y - y1)**2) / (2 * sigma1**2))

# セル2: より強い対流セル
x2, y2 = -20, -30
sigma2 = 8
cell2 = 55 * np.exp(-((X - x2)**2 + (Y - y2)**2) / (2 * sigma2**2))

# 合成 Z [dBZ]
Z_dBZ = np.maximum(cell1, cell2)
Z_dBZ = np.maximum(Z_dBZ, 5 * np.random.randn(n_azimuth, n_range) + 10)
Z_dBZ = np.clip(Z_dBZ, 0, 65)

# 降水量に変換
R_rain = z_to_rain(Z_dBZ)

# PPIプロット
fig, axes = plt.subplots(1, 2, figsize=(16, 7))

# Z [dBZ]
im1 = axes[0].pcolormesh(X, Y, Z_dBZ, cmap='jet', vmin=0, vmax=60, shading='auto')
axes[0].set_xlabel('X [km]', fontsize=12)
axes[0].set_ylabel('Y [km]', fontsize=12)
axes[0].set_title('Radar Reflectivity Z [dBZ]', fontsize=13)
axes[0].set_aspect('equal')
plt.colorbar(im1, ax=axes[0], label='Z [dBZ]')

# 降水量 [mm/h]
im2 = axes[1].pcolormesh(X, Y, R_rain, cmap='Blues', vmin=0, vmax=50, shading='auto')
axes[1].set_xlabel('X [km]', fontsize=12)
axes[1].set_ylabel('Y [km]', fontsize=12)
axes[1].set_title('Estimated Rainfall Rate [mm/h]', fontsize=13)
axes[1].set_aspect('equal')
plt.colorbar(im2, ax=axes[1], label='Rainfall Rate [mm/h]')

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# 真の風速プロファイル（高度に依存）
def true_wind(height):
    """
    高度に応じた真の風速 (u, v)
    u: 東西成分 [m/s], v: 南北成分 [m/s]
    """
    u = 5 + 15 * (height / 5000)   # 高度とともに増加
    v = -3 + 8 * (height / 5000)   # 高度とともに増加
    return u, v

# VADパラメータ
elevation = 5.0  # 仰角 [deg]
elev_rad = np.deg2rad(elevation)
n_heights = 20   # 推定する高度数
max_range = 50   # 最大距離 [km]

# 方位角（全方位走査）
n_azimuth = 360
azimuth = np.linspace(0, 360, n_azimuth, endpoint=False)
az_rad = np.deg2rad(azimuth)

# 各距離ゲート（= 高度）での風速推定
ranges = np.linspace(5, max_range, n_heights)
heights = ranges * 1000 * np.sin(elev_rad)  # 高度 [m]

u_est = np.zeros(n_heights)
v_est = np.zeros(n_heights)
u_true = np.zeros(n_heights)
v_true = np.zeros(n_heights)

for i, (r, h) in enumerate(zip(ranges, heights)):
    # 真の風速
    u_t, v_t = true_wind(h)
    u_true[i] = u_t
    v_true[i] = v_t

    # 落下速度（仮定）
    vt = -5.0  # [m/s]

    # 視線速度の生成（ノイズ付き）
    vr = (u_t * np.sin(az_rad) + v_t * np.cos(az_rad)) * np.cos(elev_rad) \
         + (0 + vt) * np.sin(elev_rad)
    noise = np.random.randn(n_azimuth) * 1.0  # 観測ノイズ
    vr_obs = vr + noise

    # フーリエ解析（1次のsin, cos成分を抽出）
    A1 = (2 / n_azimuth) * np.sum(vr_obs * np.sin(az_rad))
    B1 = (2 / n_azimuth) * np.sum(vr_obs * np.cos(az_rad))

    # 風速の推定
    u_est[i] = A1 / np.cos(elev_rad)
    v_est[i] = B1 / np.cos(elev_rad)

# プロット
fig, axes = plt.subplots(1, 3, figsize=(16, 6))

# 東西成分
axes[0].plot(u_true, heights / 1000, 'b-', linewidth=2, label='True')
axes[0].plot(u_est, heights / 1000, 'ro', markersize=6, label='VAD estimate')
axes[0].set_xlabel('u (E-W) [m/s]', fontsize=13)
axes[0].set_ylabel('Height [km]', fontsize=13)
axes[0].set_title('East-West Wind Component', fontsize=13)
axes[0].legend(fontsize=11)
axes[0].grid(True, alpha=0.3)

# 南北成分
axes[1].plot(v_true, heights / 1000, 'b-', linewidth=2, label='True')
axes[1].plot(v_est, heights / 1000, 'ro', markersize=6, label='VAD estimate')
axes[1].set_xlabel('v (N-S) [m/s]', fontsize=13)
axes[1].set_ylabel('Height [km]', fontsize=13)
axes[1].set_title('North-South Wind Component', fontsize=13)
axes[1].legend(fontsize=11)
axes[1].grid(True, alpha=0.3)

# 風速・風向
Vh_true = np.sqrt(u_true**2 + v_true**2)
Vh_est = np.sqrt(u_est**2 + v_est**2)

axes[2].plot(Vh_true, heights / 1000, 'b-', linewidth=2, label='True')
axes[2].plot(Vh_est, heights / 1000, 'ro', markersize=6, label='VAD estimate')
axes[2].set_xlabel('Wind Speed [m/s]', fontsize=13)
axes[2].set_ylabel('Height [km]', fontsize=13)
axes[2].set_title('Horizontal Wind Speed Profile', fontsize=13)
axes[2].legend(fontsize=11)
axes[2].grid(True, alpha=0.3)

plt.suptitle(f'VAD Wind Profile Retrieval (Elevation = {elevation}°)', fontsize=14)
plt.tight_layout()
plt.show()

# 推定精度の評価
rmse_u = np.sqrt(np.mean((u_est - u_true)**2))
rmse_v = np.sqrt(np.mean((v_est - v_true)**2))
print(f"RMSE (u): {rmse_u:.2f} m/s")
print(f"RMSE (v): {rmse_v:.2f} m/s")

import numpy as np
import matplotlib.pyplot as plt

# 各降水種別の偏波パラメータの統計量（平均, 標準偏差）
hydrometeor_types = {
    'Drizzle': {'Zh_mean': 15, 'Zh_std': 5,
                'Zdr_mean': 0.3, 'Zdr_std': 0.2,
                'rhohv_mean': 0.995, 'rhohv_std': 0.005,
                'color': 'cyan', 'n': 200},
    'Light Rain': {'Zh_mean': 28, 'Zh_std': 5,
                   'Zdr_mean': 0.8, 'Zdr_std': 0.3,
                   'rhohv_mean': 0.99, 'rhohv_std': 0.008,
                   'color': 'blue', 'n': 300},
    'Heavy Rain': {'Zh_mean': 48, 'Zh_std': 5,
                   'Zdr_mean': 2.5, 'Zdr_std': 0.8,
                   'rhohv_mean': 0.985, 'rhohv_std': 0.01,
                   'color': 'green', 'n': 200},
    'Dry Snow': {'Zh_mean': 25, 'Zh_std': 8,
                 'Zdr_mean': 0.3, 'Zdr_std': 0.3,
                 'rhohv_mean': 0.985, 'rhohv_std': 0.015,
                 'color': 'gray', 'n': 200},
    'Wet Snow': {'Zh_mean': 35, 'Zh_std': 7,
                 'Zdr_mean': 1.5, 'Zdr_std': 0.8,
                 'rhohv_mean': 0.9, 'rhohv_std': 0.04,
                 'color': 'orange', 'n': 200},
    'Hail': {'Zh_mean': 55, 'Zh_std': 8,
             'Zdr_mean': 0.2, 'Zdr_std': 1.0,
             'rhohv_mean': 0.88, 'rhohv_std': 0.05,
             'color': 'red', 'n': 150},
}

fig, axes = plt.subplots(1, 3, figsize=(18, 5))

for name, params in hydrometeor_types.items():
    n = params['n']
    Zh = np.random.normal(params['Zh_mean'], params['Zh_std'], n)
    Zdr = np.random.normal(params['Zdr_mean'], params['Zdr_std'], n)
    rhohv = np.clip(
        np.random.normal(params['rhohv_mean'], params['rhohv_std'], n),
        0.6, 1.0
    )

    axes[0].scatter(Zh, Zdr, c=params['color'], s=10, alpha=0.4, label=name)
    axes[1].scatter(Zh, rhohv, c=params['color'], s=10, alpha=0.4, label=name)
    axes[2].scatter(Zdr, rhohv, c=params['color'], s=10, alpha=0.4, label=name)

axes[0].set_xlabel('$Z_H$ [dBZ]', fontsize=13)
axes[0].set_ylabel('$Z_{DR}$ [dB]', fontsize=13)
axes[0].set_title('$Z_H$ vs $Z_{DR}$', fontsize=13)
axes[0].legend(fontsize=9, markerscale=3)
axes[0].grid(True, alpha=0.3)

axes[1].set_xlabel('$Z_H$ [dBZ]', fontsize=13)
axes[1].set_ylabel(r'$\rho_{HV}$', fontsize=13)
axes[1].set_title(r'$Z_H$ vs $\rho_{HV}$', fontsize=13)
axes[1].grid(True, alpha=0.3)

axes[2].set_xlabel('$Z_{DR}$ [dB]', fontsize=13)
axes[2].set_ylabel(r'$\rho_{HV}$', fontsize=13)
axes[2].set_title(r'$Z_{DR}$ vs $\rho_{HV}$', fontsize=13)
axes[2].grid(True, alpha=0.3)

plt.suptitle('Dual-Polarization Parameters by Hydrometeor Type', fontsize=14)
plt.tight_layout()
plt.show()