Editing Openai/6897769e-4ee4-800f-aba5-69cca34f701c (section)

=== Below is a full Python script. It is written to be robust: if you already have two CSV files, place them in the same folder and name them solar_spectrum.csv and water_n_lambda.csv. If you want me to fetch the online files, tell me and I’ll fetch and run. ===

Save the code as rainbow_weighted_angle.py and run with Python 3 (requires numpy, scipy, matplotlib, pandas).

<syntaxhighlight lang="python"># rainbow_weighted_angle.py
import numpy as np
import pandas as pd
from scipy.optimize import minimize_scalar
import matplotlib.pyplot as plt
from scipy.constants import c, h, k
from scipy.interpolate import interp1d

=== ---------- User inputs ---------- ===
=== If you already have files: ===
SOLAR_FILE = 'solar_spectrum.csv'      # columns: wavelength_nm, irradiance (any units)
WATER_N_FILE = 'water_n_lambda.csv'    # columns: wavelength_nm, n  (real refractive index)
=== If you don't have files, you can create a Planck spectrum below. ===

=== Emission lines to test (wavelength nm, amplitude relative to continuum, sigma nm) ===
TEST_LINES = [
    (656.28, 10.0, 0.5),  # H-alpha (strong test)
    (486.13, 5.0, 0.5),   # H-beta
    (589.3, 5.0, 0.4),    # Na D
]

=== Wavelength range for calculation ===
lam_min_nm = 350
lam_max_nm = 1100
dlam = 1.0   # nm resolution

=== --------- Helper functions ----------- ===
deg = 180.0/np.pi

def load_solar_spectrum(filename):
    df = pd.read_csv(filename)
    # expect columns named wavelength_nm and irradiance (if different adjust here)
    if 'wavelength_nm' not in df.columns:
        raise ValueError('solar_spectrum.csv must have wavelength_nm column')
    wl = np.array(df['wavelength_nm'], dtype=float)
    I = np.array(df.iloc[:,1], dtype=float)  # second column as irradiance
    f = interp1d(wl, I, kind='linear', bounds_error=False, fill_value=0.0)
    return f

def load_water_n(filename):
    df = pd.read_csv(filename)
    wl = np.array(df['wavelength_nm'], dtype=float)
    n = np.array(df['n'], dtype=float)
    fn = interp1d(wl, n, kind='cubic', bounds_error=False, fill_value='extrapolate')
    return fn

=== Planck spectral radiance (per nm) for convenience if user lacks measured solar file ===
def planck_lambda_nm(T, wl_nm):
    wl_m = wl_nm * 1e-9
    B = (2''h''c'''2) / (wl_m'''5) / (np.exp(h''c/(wl_m''k*T)) - 1.0)
    # convert to W sr^-1 m^-2 m^-1 -> convert m^-1 to nm^-1 by *1e-9
    return B * 1e-9

=== Deviation function D(theta_i, n) ===
def deviation(theta_i, n):
    # theta_i in radians
    s = np.sin(theta_i) / n
    # ensure domain [-1,1]
    s = np.clip(s, -0.9999999999, 0.9999999999)
    theta_t = np.arcsin(s)
    D = np.pi + 2.0''theta_i - 4.0''theta_t
    return D

=== get D_min for a wavelength (i.e. for a given n) ===
def D_min_for_n(n):
    # minimize D(theta) over theta in (0, pi/2)
    f = lambda th: deviation(th, n)
    # we want the minimum; use bracket from small angle to near 90deg
    res = minimize_scalar(f, bounds=(1e-6, np.pi/2 - 1e-6), method='bounded')
    if not res.success:
        raise RuntimeError('minimization failed')
    return res.fun  # this is D_min in radians

=== ---------- Main workflow ---------- ===
def compute_weighted_mean(wl_nm, Ddeg, spectrum):
    # wl_nm: array, Ddeg: array same shape (degrees), spectrum: irradiance function or array
    I = spectrum(wl_nm)
    # ensure non-negative
    I = np.clip(I, 0, None)
    num = np.trapz(Ddeg * I, wl_nm)
    den = np.trapz(I, wl_nm)
    return num/den if den>0 else np.nan

def add_lines_to_spectrum(base_I, wl, lines):
    I = base_I.copy()
    for (lam0, amp, sigma) in lines:
        I += amp '' np.exp(-0.5''((wl - lam0)/sigma)**2)
    return I

def main():
    wl = np.arange(lam_min_nm, lam_max_nm + dlam, dlam)
    # --- load or create spectrum ---
    try:
        solar_f = load_solar_spectrum(SOLAR_FILE)
        base_I = solar_f(wl)
        print("Loaded solar_spectrum.csv")
    except Exception as e:
        print("Could not load solar_spectrum.csv; using 5800K Planck as fallback.")
        base_I = planck_lambda_nm(5800.0, wl)
        # normalize for convenience
        base_I /= np.max(base_I)

    # --- load refractive index for water ---
    try:
        n_f = load_water_n(WATER_N_FILE)
        n_vals = n_f(wl)
        print("Loaded water_n_lambda.csv")
    except Exception as e:
        print("Could not load water_n_lambda.csv; using simple Cauchy approx.")
        # simple Cauchy (not high precision): n = A + B/(wl_um)^2
        wl_um = wl*1e-3
        A = 1.322
        B = 0.003  # small value
        n_vals = A + B/(wl_um**2)

    # compute D_min(lambda) numerically for each lambda
    Dmins_rad = np.zeros_like(wl, dtype=float)
    for i, n in enumerate(n_vals):
        if n <= 1.0:
            Dmins_rad[i] = np.nan
            continue
        Dmins_rad[i] = D_min_for_n(n)
    Dmins_deg = Dmins_rad * deg

    # compute weighted mean for base spectrum
    base_interp = lambda x: np.interp(x, wl, base_I)
    mean_base = compute_weighted_mean(wl, Dmins_deg, base_interp)
    print("Weighted mean D (base):", mean_base, "deg")

    # add lines and recompute
    I_with_lines = add_lines_to_spectrum(base_I, wl, TEST_LINES)
    mean_lines = compute_weighted_mean(wl, Dmins_deg, lambda x: np.interp(x, wl, I_with_lines))
    print("Weighted mean D (with lines):", mean_lines, "deg")

    # display plot
    plt.figure(figsize=(10,6))
    plt.subplot(2,1,1)
    plt.plot(wl, Dmins_deg, label='D_min(λ) (deg)')
    plt.axhline(137.5078, color='k', linestyle='--', label='Golden-angle 137.5078°')
    plt.xlabel('Wavelength (nm)')
    plt.ylabel('D_min (deg)')
    plt.legend()
    plt.subplot(2,1,2)
    plt.plot(wl, base_I/np.max(base_I), label='Base spectrum (norm)')
    plt.plot(wl, I_with_lines/np.max(I_with_lines), label='Spectrum + lines (norm)')
    plt.xlabel('Wavelength (nm)')
    plt.ylabel('Relative irradiance')
    plt.legend()
    plt.tight_layout()
    plt.show()

    # print also anti-solar angle (rainbow apparent angle = 180 - D)
    print("Apparent rainbow angle (base) = ", 180.0 - mean_base, "deg")
    print("Apparent rainbow angle (with lines) = ", 180.0 - mean_lines, "deg")

if __name__ == '__main__':
    main()

</syntaxhighlight>

Notes for using the script:
* solar_spectrum.csv should be a two-column CSV (wavelength_nm, irradiance). NREL’s ASTM G173 files can be easily converted or downloaded as CSV. If not present the script falls back to a 5800-K Planck curve.
* water_n_lambda.csv should be two columns (wavelength_nm, n). Refractiveindex.info exports or Hale & Querry tables are ideal for this.
* TEST_LINES contains emission lines you requested; you can edit amplitudes and widths freely to see how strong / narrow lines pull the mean.