增加模块；增加主调用命令

prosail_method/modules/Jfunc.py

@@ -0,0 +1,54 @@
+import numpy as np
+
+def Jfunc1(k, l, t):
+    """
+    Computes J1(k, l, t) = (exp(-lt) - exp(-kt)) / (k - l)
+    Uses Taylor expansion when k ≈ l.
+    
+    k: scalar or array
+    l: array
+    t: scalar or array
+    """
+    # 广播所有变量到相同 shape（自动支持 k 是标量、l 是向量）
+    k, l, t_array = np.broadcast_arrays(k, l, t)
+
+    delta = (k - l) * t_array
+    Jout = np.zeros_like(delta)
+
+    mask_far = np.abs(delta) > 1e-3
+    mask_near = ~mask_far
+
+    # 正常计算
+    Jout[mask_far] = (
+        np.exp(-l[mask_far] * t_array[mask_far]) -
+        np.exp(-k[mask_far] * t_array[mask_far])
+    ) / (k[mask_far] - l[mask_far])
+
+    # k ≈ l 的情况，使用二阶展开
+    d = delta[mask_near]
+    Jout[mask_near] = 0.5 * t_array[mask_near] * (
+        np.exp(-k[mask_near] * t_array[mask_near]) +
+        np.exp(-l[mask_near] * t_array[mask_near])
+    ) * (1 - d * d / 12)
+
+    return Jout
+
+
+
+
+def Jfunc2(k, l, t):
+    k, l, t_array = np.broadcast_arrays(k, l, t)
+    delta = (k + l) * t_array
+    Jout = np.zeros_like(k)
+
+    mask_far = np.abs(delta) > 1e-3
+    mask_near = ~mask_far
+
+    Jout[mask_far] = (1 - np.exp(-delta[mask_far])) / (k[mask_far] + l[mask_far])
+    d = delta[mask_near]
+    Jout[mask_near] = t_array[mask_near] * (1 - 0.5 * d + d ** 2 / 6)
+
+    return Jout
+
+def Jfunc3(k, l, t):
+    return Jfunc2(k, l, t)

prosail_method/modules/calctav.py

@@ -0,0 +1,46 @@
+import numpy as np
+
+def calctav(alfa, nr):
+    """
+    Calculate average transmittance of isotropic light across a dielectric surface.
+    Translated from MATLAB version of calctav.m.
+
+    Parameters:
+        alfa : incidence angle in degrees (scalar or array)
+        nr   : refractive index of second medium (scalar or array of same shape as alfa)
+
+    Returns:
+        tav : average transmittance (same shape as alfa)
+    """
+    alfa = np.asarray(alfa)
+    nr = np.asarray(nr)
+    rd = np.pi / 180.0
+
+    n2 = nr ** 2
+    np_ = n2 + 1
+    nm = n2 - 1
+    a = ((nr + 1) ** 2) / 2
+    k = -((n2 - 1) ** 2) / 4
+    sa = np.sin(alfa * rd)
+
+    b1 = np.where(alfa != 90,
+                  np.sqrt((sa ** 2 - np_ / 2) ** 2 + k),
+                  0.0)
+    b2 = sa ** 2 - np_ / 2
+    b = b1 - b2
+    b3 = b ** 3
+    a3 = a ** 3
+
+    ts = (k ** 2 / (6 * b3) + k / b - b / 2) - (k ** 2 / (6 * a3) + k / a - a / 2)
+
+    tp1 = -2 * n2 * (b - a) / (np_ ** 2)
+    tp2 = -2 * n2 * np_ * np.log(b / a) / (nm ** 2)
+    tp3 = n2 * (1 / b - 1 / a) / 2
+    tp4 = 16 * n2 ** 2 * (n2 ** 2 + 1) * np.log((2 * np_ * b - nm ** 2) / (2 * np_ * a - nm ** 2)) / (np_ ** 3 * nm ** 2)
+    tp5 = 16 * n2 ** 3 * (1 / (2 * np_ * b - nm ** 2) - 1 / (2 * np_ * a - nm ** 2)) / (np_ ** 3)
+
+    tp = tp1 + tp2 + tp3 + tp4 + tp5
+    tav = (ts + tp) / (2 * sa ** 2)
+
+    return tav
+

prosail_method/modules/campbell.py

@@ -0,0 +1,53 @@
+import numpy as np
+
+def campbell(ala):
+    """
+    Campbell's leaf angle distribution function (ellipsoidal).
+    
+    Parameters:
+        ala : average leaf angle in degrees
+
+    Returns:
+        freq : relative frequency for 13 inclination bins
+        litab: bin centers (degrees)
+    """
+    tx1 = np.array([10., 20., 30., 40., 50., 60., 70., 80., 82., 84., 86., 88., 90.])
+    tx2 = np.array([0., 10., 20., 30., 40., 50., 60., 70., 80., 82., 84., 86., 88.])
+    litab = (tx1 + tx2) / 2
+    n = len(litab)
+
+    tl1 = np.radians(tx1)
+    tl2 = np.radians(tx2)
+
+    excent = np.exp(-1.6184e-5 * ala**3 + 2.1145e-3 * ala**2 - 0.12390 * ala + 3.2491)
+
+    freq = np.zeros(n)
+    for i in range(n):
+        x1 = excent / np.sqrt(1 + (excent**2) * np.tan(tl1[i])**2)
+        x2 = excent / np.sqrt(1 + (excent**2) * np.tan(tl2[i])**2)
+
+        if excent == 1:
+            freq[i] = abs(np.cos(tl1[i]) - np.cos(tl2[i]))
+        else:
+            alpha = excent / np.sqrt(abs(1 - excent**2))
+            alpha2 = alpha**2
+            x1_sq = x1**2
+            x2_sq = x2**2
+
+            if excent > 1:
+                alpx1 = np.sqrt(alpha2 + x1_sq)
+                alpx2 = np.sqrt(alpha2 + x2_sq)
+                dum1 = x1 * alpx1 + alpha2 * np.log(x1 + alpx1)
+                dum2 = x2 * alpx2 + alpha2 * np.log(x2 + alpx2)
+                freq[i] = abs(dum1 - dum2)
+            else:
+                almx1 = np.sqrt(alpha2 - x1_sq)
+                almx2 = np.sqrt(alpha2 - x2_sq)
+                dum1 = x1 * almx1 + alpha2 * np.arcsin(x1 / alpha)
+                dum2 = x2 * almx2 + alpha2 * np.arcsin(x2 / alpha)
+                freq[i] = abs(dum1 - dum2)
+
+    freq_sum = np.sum(freq)
+    freq_normalized = freq / freq_sum if freq_sum != 0 else freq
+
+    return freq_normalized, litab

prosail_method/modules/dcum.py

@@ -0,0 +1,29 @@
+import numpy as np
+
+def dcum(a, b, t):
+    """
+    计算 PROSAIL 中使用的累积分布函数 dcum
+    
+    Parameters:
+    - a: float
+    - b: float
+    - t: float, in degrees
+
+    Returns:
+    - f: float
+    """
+    rd = np.pi / 180  # degree to radian
+    if a >= 1:
+        f = 1 - np.cos(rd * t)
+    else:
+        eps = 1e-8
+        x = 2 * rd * t
+        p = x
+        delx = 1.0
+        while delx >= eps:
+            y = a * np.sin(x) + 0.5 * b * np.sin(2 * x)
+            dx = 0.5 * (y - x + p)
+            x += dx
+            delx = abs(dx)
+        f = (2 * y + p) / np.pi
+    return f

prosail_method/modules/dladgen.py

@@ -0,0 +1,35 @@
+import numpy as np
+from modules.dcum import dcum  # 确保你有这个函数
+
+def dladgen(a, b):
+    """
+    Generate Leaf Inclination Distribution Function (LIDF)
+    using a two-parameter beta distribution.
+
+    Parameters:
+        a, b : shape parameters of beta distribution
+
+    Returns:
+        freq  : frequency (probability) per inclination class (length 13)
+        litab : central inclination angles (degrees)
+    """
+    litab = np.array([5., 15., 25., 35., 45., 55., 65., 75., 81., 83., 85., 87., 89.])
+    freq = np.zeros_like(litab)
+
+    # First 8 bins: 10° steps
+    for i in range(8):
+        t = (i + 1) * 10
+        freq[i] = dcum(a, b, t)
+
+    # Next 4 bins: finer steps from 80° to 89°
+    for i in range(8, 12):
+        t = 80.0 + (i - 8 + 1) * 2.0
+        freq[i] = dcum(a, b, t)
+
+    freq[12] = 1.0  # total cumulative = 1
+
+    # Convert cumulative to class probabilities
+    for i in range(12, 0, -1):
+        freq[i] -= freq[i - 1]
+
+    return freq, litab

prosail_method/modules/load_coefficients.py

@@ -0,0 +1,32 @@
+import numpy as np
+
+def load_coefficients(filepath):
+    """
+    Load spectral coefficients from a full-format PROSAIL data file (e.g., dataSpec_PDB.txt).
+
+    Returns:
+        wl     : ndarray, wavelength (nm)
+        nr     : ndarray, refractive index of leaf material
+        Kab    : ndarray, absorption coefficient of chlorophyll (a+b)
+        Kcar   : ndarray, absorption coefficient of carotenoids
+        Kant   : ndarray, absorption coefficient of anthocyanins
+        Kbrown : ndarray, absorption coefficient of brown pigments
+        Kw     : ndarray, absorption coefficient of water
+        Km     : ndarray, absorption coefficient of dry matter
+        Rsoil1 : ndarray, dry soil reflectance
+        Rsoil2 : ndarray, wet soil reflectance
+    """
+    data = np.loadtxt(filepath, comments='%', delimiter=None)
+
+    wl     = data[:, 0]
+    nr     = data[:, 1]
+    Kab    = data[:, 2]
+    Kcar   = data[:, 3]
+    Kant   = data[:, 4]
+    Kbrown = data[:, 5]
+    Kw     = data[:, 6]
+    Km     = data[:, 7]
+    Rsoil1 = data[:, 10]
+    Rsoil2 = data[:, 11]
+
+    return wl, nr, Kab, Kcar, Kant, Kbrown, Kw, Km, Rsoil1, Rsoil2

prosail_method/modules/prosail.py

@@ -0,0 +1,23 @@
+from modules.prospect_d import prospect_d
+from modules.sail import sail
+
+def prosail(N, Cab, Car, Ant, Cbrown, Cw, Cm,
+            LIDFa, LIDFb, TypeLidf, lai, hspot,
+            tts, tto, psi, rsoil,
+            wl, nr, Kab, Kcar, Kant, Kbrown, Kw, Km):
+    """
+    PROSAIL simulator using PROSPECT-D + 4SAIL.
+
+    Returns:
+        refl: canopy directional reflectance
+        LRT : leaf reflectance/transmittance spectra
+    """
+    # Step 1: simulate leaf optical properties using PROSPECT-D
+    LRT = prospect_d(N, Cab, Car, Ant, Cbrown, Cw, Cm, wl, nr, Kab, Kcar, Kant, Kbrown, Kw, Km)
+    rho = LRT[:, 1]
+    tau = LRT[:, 2]
+
+    # Step 2: simulate canopy reflectance using 4SAIL
+    refl = sail(rho, tau, LIDFa, LIDFb, TypeLidf, lai, hspot, tts, tto, psi, rsoil, wl)
+
+    return wl, refl, LRT

prosail_method/modules/prospect_d.py

@@ -0,0 +1,80 @@
+import numpy as np
+from scipy.special import exp1
+from modules.calctav import calctav  # 你自己已有的函数
+
+def prospect_d(N, Cab, Car, Ant, Brown, Cw, Cm, wl, nr, Kab, Kcar, Kant, Kbrown, Kw, Km):
+    """
+    Accurate PROSPECT-D implementation matching the original MATLAB model.
+
+    Parameters:
+        N       : Structure parameter
+        Cab     : Chlorophyll content (µg/cm²)
+        Car     : Carotenoid content (µg/cm²)
+        Ant     : Anthocyanins content (µg/cm²)
+        Brown   : Brown pigment content (a.u.)
+        Cw      : Equivalent water thickness (cm)
+        Cm      : Dry matter content (g/cm²)
+        wl      : Wavelength array (nm)
+        nr      : Real refractive index
+        Kab,... : Absorption coefficients for components (same length as wl)
+
+    Returns:
+        LRT : ndarray of shape (n, 3) → [wavelength, reflectance, transmittance]
+    """
+
+    # 1. Compute total absorption coefficient
+    Kall = (Cab * Kab + Car * Kcar + Ant * Kant + Brown * Kbrown + Cw * Kw + Cm * Km) / N
+
+    # 2. Compute tau (internal transmittance) with expint for K > 0
+    t1 = (1 - Kall) * np.exp(-Kall)
+    t2 = (Kall ** 2) * exp1(Kall)
+    tau = np.ones_like(Kall)
+    j = Kall > 0
+    tau[j] = t1[j] + t2[j]
+
+    # 3. Fresnel interface properties
+    talf = calctav(40, nr)
+    ralf = 1 - talf
+    t12 = calctav(90, nr)
+    r12 = 1 - t12
+    t21 = t12 / (nr ** 2)
+    r21 = 1 - t21
+
+    # 4. Top surface
+    denom = 1 - (r21 ** 2) * (tau ** 2)
+    Ta = talf * tau * t21 / denom
+    Ra = ralf + r21 * tau * Ta
+
+    # 5. Bottom surface
+    t = t12 * tau * t21 / denom
+    r = r12 + r21 * tau * t
+
+    # 6. Stokes multilayer reflectance and transmittance
+    D = np.sqrt((1 + r + t) * (1 + r - t) * (1 - r + t) * (1 - r - t))
+    rq = r ** 2
+    tq = t ** 2
+    a = (1 + rq - tq + D) / (2 * r)
+    b = (1 - rq + tq + D) / (2 * t)
+
+    bNm1 = b ** (N - 1)
+    bN2 = bNm1 ** 2
+    a2 = a ** 2
+    denom2 = a2 * bN2 - 1
+    Rsub = a * (bN2 - 1) / denom2
+    Tsub = bNm1 * (a2 - 1) / denom2
+
+    # 7. Zero absorption edge-case correction
+    j0 = (r + t >= 1)
+    Tsub[j0] = t[j0] / (t[j0] + (1 - t[j0]) * (N - 1))
+    Rsub[j0] = 1 - Tsub[j0]
+
+    # 8. Final reflectance and transmittance
+    denom3 = 1 - Rsub * r
+    tran = Ta * Tsub / denom3
+    refl = Ra + Ta * Rsub * t / denom3
+
+    # 9. Clip values to [0, 1]
+    refl = np.clip(refl, 0, 1)
+    tran = np.clip(tran, 0, 1)
+
+    return np.column_stack((wl, refl, tran))

prosail_method/modules/sail.py

@@ -0,0 +1,141 @@
+import numpy as np
+from modules.volscatt import volscatt
+from modules.dladgen import dladgen
+from modules.campbell import campbell
+from modules.Jfunc import Jfunc1, Jfunc2, Jfunc3
+
+def sail(rho, tau, LIDFa, LIDFb, TypeLidf, lai, q, tts, tto, psi, rsoil,wl):
+    rd = np.pi / 180
+    cts = np.cos(rd * tts)
+    cto = np.cos(rd * tto)
+    ctscto = cts * cto
+    tants = np.tan(rd * tts)
+    tanto = np.tan(rd * tto)
+    cospsi = np.cos(rd * psi)
+    dso = np.sqrt(tants**2 + tanto**2 - 2 * tants * tanto * cospsi)
+
+    # Generate LIDF
+    if TypeLidf == 1:
+        lidf, litab = dladgen(LIDFa, LIDFb)
+    elif TypeLidf == 2:
+        lidf, litab = campbell(LIDFa)
+
+    # Initialize geometric coefficients
+    ks = ko = bf = sob = sof = 0
+    na = len(litab)
+
+    for i in range(na):
+        ttl = litab[i]
+        ctl = np.cos(rd * ttl)
+        chi_s, chi_o, frho, ftau = volscatt(tts, tto, psi, ttl)
+
+        ksli = chi_s / cts
+        koli = chi_o / cto
+        sobli = frho * np.pi / ctscto
+        sofli = ftau * np.pi / ctscto
+        bfli = ctl ** 2
+
+        ks += ksli * lidf[i]
+        ko += koli * lidf[i]
+        bf += bfli * lidf[i]
+        sob += sobli * lidf[i]
+        sof += sofli * lidf[i]
+
+    # Structure factors
+    sdb = 0.5 * (ks + bf)
+    sdf = 0.5 * (ks - bf)
+    dob = 0.5 * (ko + bf)
+    dof = 0.5 * (ko - bf)
+    ddb = 0.5 * (1 + bf)
+    ddf = 0.5 * (1 - bf)
+
+    sigb = ddb * rho + ddf * tau
+    sigf = ddf * rho + ddb * tau
+    att = 1 - sigf
+    m2 = (att + sigb) * (att - sigb)
+    m2[m2 <= 0] = 0
+    m = np.sqrt(m2)
+
+    sb = sdb * rho + sdf * tau
+    sf = sdf * rho + sdb * tau
+    vb = dob * rho + dof * tau
+    vf = dof * rho + dob * tau
+    w = sob * rho + sof * tau
+
+    if lai <= 0:
+        return rsoil, rsoil, rsoil, rsoil
+
+    e1 = np.exp(-m * lai)
+    e2 = e1**2
+    rinf = (att - m) / sigb
+    rinf2 = rinf**2
+    re = rinf * e1
+    denom = 1 - rinf2 * e2
+
+    J1ks = Jfunc1(ks, m, lai)
+    J2ks = Jfunc2(ks, m, lai)
+    J1ko = Jfunc1(ko, m, lai)
+    J2ko = Jfunc2(ko, m, lai)
+
+    Ps = (sf + sb * rinf) * J1ks
+    Qs = (sf * rinf + sb) * J2ks
+    Pv = (vf + vb * rinf) * J1ko
+    Qv = (vf * rinf + vb) * J2ko
+
+    rdd = rinf * (1 - e2) / denom
+    tdd = (1 - rinf2) * e1 / denom
+    tsd = (Ps - re * Qs) / denom
+    rsd = (Qs - re * Ps) / denom
+    tdo = (Pv - re * Qv) / denom
+    rdo = (Qv - re * Pv) / denom
+
+    tss = np.exp(-ks * lai)
+    too = np.exp(-ko * lai)
+    z = Jfunc3(ks, ko, lai)
+    g1 = (z - J1ks * too) / (ko + m)
+    g2 = (z - J1ko * tss) / (ks + m)
+
+    Tv1 = (vf * rinf + vb) * g1
+    Tv2 = (vf + vb * rinf) * g2
+    T1 = Tv1 * (sf + sb * rinf)
+    T2 = Tv2 * (sf * rinf + sb)
+    T3 = (rdo * Qs + tdo * Ps) * rinf
+    rsod = (T1 + T2 - T3) / (1 - rinf2)
+
+    # Hotspot effect
+    if q <= 0:
+        alf = 1e6
+    else:
+        alf = min((dso / q) * 2 / (ks + ko), 200)
+
+    if alf == 0:
+        tsstoo = tss
+        sumint = (1 - tss) / (ks * lai)
+    else:
+        fhot = lai * np.sqrt(ko * ks)
+        fint = (1 - np.exp(-alf)) * 0.05
+        sumint = 0
+        x1 = y1 = 0
+        f1 = 1
+        for i in range(20):
+            x2 = -np.log(1 - i * fint) / alf if i < 19 else 1
+            y2 = -(ko + ks) * lai * x2 + fhot * (1 - np.exp(-alf * x2)) / alf
+            f2 = np.exp(y2)
+            if y2 != y1:
+                sumint += (f2 - f1) * (x2 - x1) / (y2 - y1)
+            x1, y1, f1 = x2, y2, f2
+        tsstoo = f1
+
+    rsos = w * lai * sumint
+    rso = rsos + rsod
+
+    # Canopy–soil interaction
+    dn = 1 - rsoil * rdd
+    rddt = rdd + tdd * rsoil * tdd / dn
+    rsdt = rsd + (tsd + tss) * rsoil * tdd / dn
+    rdot = rdo + tdd * rsoil * (tdo + too) / dn
+    rsodt = rsod + ((tss + tsd) * tdo + (tsd + tss * rsoil * rdd) * too) * rsoil / dn
+    rsost = rsos + tsstoo * rsoil
+    rsot = rsost + rsodt
+
+    return np.column_stack((wl, rdot, rsot, rddt, rsdt))

prosail_method/modules/tmp7sub_ja4

@@ -0,0 +1,141 @@
+import numpy as np
+from modules.volscatt import volscatt
+from modules.dladgen import dladgen
+from modules.campbell import campbell
+from modules.Jfunc import Jfunc1, Jfunc2, Jfunc3
+
+def sail(rho, tau, LIDFa, LIDFb, TypeLidf, lai, q, tts, tto, psi, rsoil):
+    rd = np.pi / 180
+    cts = np.cos(rd * tts)
+    cto = np.cos(rd * tto)
+    ctscto = cts * cto
+    tants = np.tan(rd * tts)
+    tanto = np.tan(rd * tto)
+    cospsi = np.cos(rd * psi)
+    dso = np.sqrt(tants**2 + tanto**2 - 2 * tants * tanto * cospsi)
+
+    # Generate LIDF
+    if TypeLidf == 1:
+        lidf, litab = dladgen(LIDFa, LIDFb)
+    elif TypeLidf == 2:
+        lidf, litab = campbell(LIDFa)
+
+    # Initialize geometric coefficients
+    ks = ko = bf = sob = sof = 0
+    na = len(litab)
+
+    for i in range(na):
+        ttl = litab[i]
+        ctl = np.cos(rd * ttl)
+        chi_s, chi_o, frho, ftau = volscatt(tts, tto, psi, ttl)
+
+        ksli = chi_s / cts
+        koli = chi_o / cto
+        sobli = frho * np.pi / ctscto
+        sofli = ftau * np.pi / ctscto
+        bfli = ctl ** 2
+
+        ks += ksli * lidf[i]
+        ko += koli * lidf[i]
+        bf += bfli * lidf[i]
+        sob += sobli * lidf[i]
+        sof += sofli * lidf[i]
+
+    # Structure factors
+    sdb = 0.5 * (ks + bf)
+    sdf = 0.5 * (ks - bf)
+    dob = 0.5 * (ko + bf)
+    dof = 0.5 * (ko - bf)
+    ddb = 0.5 * (1 + bf)
+    ddf = 0.5 * (1 - bf)
+
+    sigb = ddb * rho + ddf * tau
+    sigf = ddf * rho + ddb * tau
+    att = 1 - sigf
+    m2 = (att + sigb) * (att - sigb)
+    m2[m2 <= 0] = 0
+    m = np.sqrt(m2)
+
+    sb = sdb * rho + sdf * tau
+    sf = sdf * rho + sdb * tau
+    vb = dob * rho + dof * tau
+    vf = dof * rho + dob * tau
+    w = sob * rho + sof * tau
+
+    if lai <= 0:
+        return rsoil, rsoil, rsoil, rsoil
+
+    e1 = np.exp(-m * lai)
+    e2 = e1**2
+    rinf = (att - m) / sigb
+    rinf2 = rinf**2
+    re = rinf * e1
+    denom = 1 - rinf2 * e2
+
+    J1ks = Jfunc1(ks, m, lai)
+    J2ks = Jfunc2(ks, m, lai)
+    J1ko = Jfunc1(ko, m, lai)
+    J2ko = Jfunc2(ko, m, lai)
+
+    Ps = (sf + sb * rinf) * J1ks
+    Qs = (sf * rinf + sb) * J2ks
+    Pv = (vf + vb * rinf) * J1ko
+    Qv = (vf * rinf + vb) * J2ko
+
+    rdd = rinf * (1 - e2) / denom
+    tdd = (1 - rinf2) * e1 / denom
+    tsd = (Ps - re * Qs) / denom
+    rsd = (Qs - re * Ps) / denom
+    tdo = (Pv - re * Qv) / denom
+    rdo = (Qv - re * Pv) / denom
+
+    tss = np.exp(-ks * lai)
+    too = np.exp(-ko * lai)
+    z = Jfunc3(ks, ko, lai)
+    g1 = (z - J1ks * too) / (ko + m)
+    g2 = (z - J1ko * tss) / (ks + m)
+
+    Tv1 = (vf * rinf + vb) * g1
+    Tv2 = (vf + vb * rinf) * g2
+    T1 = Tv1 * (sf + sb * rinf)
+    T2 = Tv2 * (sf * rinf + sb)
+    T3 = (rdo * Qs + tdo * Ps) * rinf
+    rsod = (T1 + T2 - T3) / (1 - rinf2)
+
+    # Hotspot effect
+    if q <= 0:
+        alf = 1e6
+    else:
+        alf = min((dso / q) * 2 / (ks + ko), 200)
+
+    if alf == 0:
+        tsstoo = tss
+        sumint = (1 - tss) / (ks * lai)
+    else:
+        fhot = lai * np.sqrt(ko * ks)
+        fint = (1 - np.exp(-alf)) * 0.05
+        sumint = 0
+        x1 = y1 = 0
+        f1 = 1
+        for i in range(20):
+            x2 = -np.log(1 - i * fint) / alf if i < 19 else 1
+            y2 = -(ko + ks) * lai * x2 + fhot * (1 - np.exp(-alf * x2)) / alf
+            f2 = np.exp(y2)
+            if y2 != y1:
+                sumint += (f2 - f1) * (x2 - x1) / (y2 - y1)
+            x1, y1, f1 = x2, y2, f2
+        tsstoo = f1
+
+    rsos = w * lai * sumint
+    rso = rsos + rsod
+
+    # Canopy–soil interaction
+    dn = 1 - rsoil * rdd
+    rddt = rdd + tdd * rsoil * tdd / dn
+    rsdt = rsd + (tsd + tss) * rsoil * tdd / dn
+    rdot = rdo + tdd * rsoil * (tdo + too) / dn
+    rsodt = rsod + ((tss + tsd) * tdo + (tsd + tss * rsoil * rdd) * too) * rsoil / dn
+    rsost = rsos + tsstoo * rsoil
+    rsot = rsost + rsodt
+
+    return np.column_stack((rdot, rsot, rddt, rsdt))

prosail_method/modules/tmpjnjd71aa

@@ -0,0 +1,23 @@
+from modules.prospect_d import prospect_d
+from modules.sail import sail
+
+def prosail(N, Cab, Car, Ant, Cbrown, Cw, Cm,
+            LIDFa, LIDFb, TypeLidf, lai, hspot,
+            tts, tto, psi, rsoil,
+            wl, nr, Kab, Kcar, Kant, Kbrown, Kw, Km):
+    """
+    PROSAIL simulator using PROSPECT-D + 4SAIL.
+
+    Returns:
+        refl: canopy directional reflectance
+        LRT : leaf reflectance/transmittance spectra
+    """
+    # Step 1: simulate leaf optical properties using PROSPECT-D
+    LRT = prospect_d(N, Cab, Car, Ant, Cbrown, Cw, Cm, wl, nr, Kab, Kcar, Kant, Kbrown, Kw, Km)
+    rho = LRT[:, 1]
+    tau = LRT[:, 2]
+
+    # Step 2: simulate canopy reflectance using 4SAIL
+    refl = sail(rho, tau, LIDFa, LIDFb, TypeLidf, lai, hspot, tts, tto, psi, rsoil, wl)
+
+    return wl, refl, LRT

prosail_method/modules/tmpkhzqjdcu

@@ -0,0 +1,141 @@
+import numpy as np
+from modules.volscatt import volscatt
+from modules.dladgen import dladgen
+from modules.campbell import campbell
+from modules.Jfunc import Jfunc1, Jfunc2, Jfunc3
+
+def sail(rho, tau, LIDFa, LIDFb, TypeLidf, lai, q, tts, tto, psi, rsoil,wl):
+    rd = np.pi / 180
+    cts = np.cos(rd * tts)
+    cto = np.cos(rd * tto)
+    ctscto = cts * cto
+    tants = np.tan(rd * tts)
+    tanto = np.tan(rd * tto)
+    cospsi = np.cos(rd * psi)
+    dso = np.sqrt(tants**2 + tanto**2 - 2 * tants * tanto * cospsi)
+
+    # Generate LIDF
+    if TypeLidf == 1:
+        lidf, litab = dladgen(LIDFa, LIDFb)
+    elif TypeLidf == 2:
+        lidf, litab = campbell(LIDFa)
+
+    # Initialize geometric coefficients
+    ks = ko = bf = sob = sof = 0
+    na = len(litab)
+
+    for i in range(na):
+        ttl = litab[i]
+        ctl = np.cos(rd * ttl)
+        chi_s, chi_o, frho, ftau = volscatt(tts, tto, psi, ttl)
+
+        ksli = chi_s / cts
+        koli = chi_o / cto
+        sobli = frho * np.pi / ctscto
+        sofli = ftau * np.pi / ctscto
+        bfli = ctl ** 2
+
+        ks += ksli * lidf[i]
+        ko += koli * lidf[i]
+        bf += bfli * lidf[i]
+        sob += sobli * lidf[i]
+        sof += sofli * lidf[i]
+
+    # Structure factors
+    sdb = 0.5 * (ks + bf)
+    sdf = 0.5 * (ks - bf)
+    dob = 0.5 * (ko + bf)
+    dof = 0.5 * (ko - bf)
+    ddb = 0.5 * (1 + bf)
+    ddf = 0.5 * (1 - bf)
+
+    sigb = ddb * rho + ddf * tau
+    sigf = ddf * rho + ddb * tau
+    att = 1 - sigf
+    m2 = (att + sigb) * (att - sigb)
+    m2[m2 <= 0] = 0
+    m = np.sqrt(m2)
+
+    sb = sdb * rho + sdf * tau
+    sf = sdf * rho + sdb * tau
+    vb = dob * rho + dof * tau
+    vf = dof * rho + dob * tau
+    w = sob * rho + sof * tau
+
+    if lai <= 0:
+        return rsoil, rsoil, rsoil, rsoil
+
+    e1 = np.exp(-m * lai)
+    e2 = e1**2
+    rinf = (att - m) / sigb
+    rinf2 = rinf**2
+    re = rinf * e1
+    denom = 1 - rinf2 * e2
+
+    J1ks = Jfunc1(ks, m, lai)
+    J2ks = Jfunc2(ks, m, lai)
+    J1ko = Jfunc1(ko, m, lai)
+    J2ko = Jfunc2(ko, m, lai)
+
+    Ps = (sf + sb * rinf) * J1ks
+    Qs = (sf * rinf + sb) * J2ks
+    Pv = (vf + vb * rinf) * J1ko
+    Qv = (vf * rinf + vb) * J2ko
+
+    rdd = rinf * (1 - e2) / denom
+    tdd = (1 - rinf2) * e1 / denom
+    tsd = (Ps - re * Qs) / denom
+    rsd = (Qs - re * Ps) / denom
+    tdo = (Pv - re * Qv) / denom
+    rdo = (Qv - re * Pv) / denom
+
+    tss = np.exp(-ks * lai)
+    too = np.exp(-ko * lai)
+    z = Jfunc3(ks, ko, lai)
+    g1 = (z - J1ks * too) / (ko + m)
+    g2 = (z - J1ko * tss) / (ks + m)
+
+    Tv1 = (vf * rinf + vb) * g1
+    Tv2 = (vf + vb * rinf) * g2
+    T1 = Tv1 * (sf + sb * rinf)
+    T2 = Tv2 * (sf * rinf + sb)
+    T3 = (rdo * Qs + tdo * Ps) * rinf
+    rsod = (T1 + T2 - T3) / (1 - rinf2)
+
+    # Hotspot effect
+    if q <= 0:
+        alf = 1e6
+    else:
+        alf = min((dso / q) * 2 / (ks + ko), 200)
+
+    if alf == 0:
+        tsstoo = tss
+        sumint = (1 - tss) / (ks * lai)
+    else:
+        fhot = lai * np.sqrt(ko * ks)
+        fint = (1 - np.exp(-alf)) * 0.05
+        sumint = 0
+        x1 = y1 = 0
+        f1 = 1
+        for i in range(20):
+            x2 = -np.log(1 - i * fint) / alf if i < 19 else 1
+            y2 = -(ko + ks) * lai * x2 + fhot * (1 - np.exp(-alf * x2)) / alf
+            f2 = np.exp(y2)
+            if y2 != y1:
+                sumint += (f2 - f1) * (x2 - x1) / (y2 - y1)
+            x1, y1, f1 = x2, y2, f2
+        tsstoo = f1
+
+    rsos = w * lai * sumint
+    rso = rsos + rsod
+
+    # Canopy–soil interaction
+    dn = 1 - rsoil * rdd
+    rddt = rdd + tdd * rsoil * tdd / dn
+    rsdt = rsd + (tsd + tss) * rsoil * tdd / dn
+    rdot = rdo + tdd * rsoil * (tdo + too) / dn
+    rsodt = rsod + ((tss + tsd) * tdo + (tsd + tss * rsoil * rdd) * too) * rsoil / dn
+    rsost = rsos + tsstoo * rsoil
+    rsot = rsost + rsodt
+
+    return np.column_stack((wl,rdot, rsot, rddt, rsdt))

prosail_method/modules/tmpsi4_5sdo

@@ -0,0 +1,141 @@
+import numpy as np
+from modules.volscatt import volscatt
+from modules.dladgen import dladgen
+from modules.campbell import campbell
+from modules.Jfunc import Jfunc1, Jfunc2, Jfunc3
+
+def sail(rho, tau, LIDFa, LIDFb, TypeLidf, lai, q, tts, tto, psi, rsoil,wl):
+    rd = np.pi / 180
+    cts = np.cos(rd * tts)
+    cto = np.cos(rd * tto)
+    ctscto = cts * cto
+    tants = np.tan(rd * tts)
+    tanto = np.tan(rd * tto)
+    cospsi = np.cos(rd * psi)
+    dso = np.sqrt(tants**2 + tanto**2 - 2 * tants * tanto * cospsi)
+
+    # Generate LIDF
+    if TypeLidf == 1:
+        lidf, litab = dladgen(LIDFa, LIDFb)
+    elif TypeLidf == 2:
+        lidf, litab = campbell(LIDFa)
+
+    # Initialize geometric coefficients
+    ks = ko = bf = sob = sof = 0
+    na = len(litab)
+
+    for i in range(na):
+        ttl = litab[i]
+        ctl = np.cos(rd * ttl)
+        chi_s, chi_o, frho, ftau = volscatt(tts, tto, psi, ttl)
+
+        ksli = chi_s / cts
+        koli = chi_o / cto
+        sobli = frho * np.pi / ctscto
+        sofli = ftau * np.pi / ctscto
+        bfli = ctl ** 2
+
+        ks += ksli * lidf[i]
+        ko += koli * lidf[i]
+        bf += bfli * lidf[i]
+        sob += sobli * lidf[i]
+        sof += sofli * lidf[i]
+
+    # Structure factors
+    sdb = 0.5 * (ks + bf)
+    sdf = 0.5 * (ks - bf)
+    dob = 0.5 * (ko + bf)
+    dof = 0.5 * (ko - bf)
+    ddb = 0.5 * (1 + bf)
+    ddf = 0.5 * (1 - bf)
+
+    sigb = ddb * rho + ddf * tau
+    sigf = ddf * rho + ddb * tau
+    att = 1 - sigf
+    m2 = (att + sigb) * (att - sigb)
+    m2[m2 <= 0] = 0
+    m = np.sqrt(m2)
+
+    sb = sdb * rho + sdf * tau
+    sf = sdf * rho + sdb * tau
+    vb = dob * rho + dof * tau
+    vf = dof * rho + dob * tau
+    w = sob * rho + sof * tau
+
+    if lai <= 0:
+        return rsoil, rsoil, rsoil, rsoil
+
+    e1 = np.exp(-m * lai)
+    e2 = e1**2
+    rinf = (att - m) / sigb
+    rinf2 = rinf**2
+    re = rinf * e1
+    denom = 1 - rinf2 * e2
+
+    J1ks = Jfunc1(ks, m, lai)
+    J2ks = Jfunc2(ks, m, lai)
+    J1ko = Jfunc1(ko, m, lai)
+    J2ko = Jfunc2(ko, m, lai)
+
+    Ps = (sf + sb * rinf) * J1ks
+    Qs = (sf * rinf + sb) * J2ks
+    Pv = (vf + vb * rinf) * J1ko
+    Qv = (vf * rinf + vb) * J2ko
+
+    rdd = rinf * (1 - e2) / denom
+    tdd = (1 - rinf2) * e1 / denom
+    tsd = (Ps - re * Qs) / denom
+    rsd = (Qs - re * Ps) / denom
+    tdo = (Pv - re * Qv) / denom
+    rdo = (Qv - re * Pv) / denom
+
+    tss = np.exp(-ks * lai)
+    too = np.exp(-ko * lai)
+    z = Jfunc3(ks, ko, lai)
+    g1 = (z - J1ks * too) / (ko + m)
+    g2 = (z - J1ko * tss) / (ks + m)
+
+    Tv1 = (vf * rinf + vb) * g1
+    Tv2 = (vf + vb * rinf) * g2
+    T1 = Tv1 * (sf + sb * rinf)
+    T2 = Tv2 * (sf * rinf + sb)
+    T3 = (rdo * Qs + tdo * Ps) * rinf
+    rsod = (T1 + T2 - T3) / (1 - rinf2)
+
+    # Hotspot effect
+    if q <= 0:
+        alf = 1e6
+    else:
+        alf = min((dso / q) * 2 / (ks + ko), 200)
+
+    if alf == 0:
+        tsstoo = tss
+        sumint = (1 - tss) / (ks * lai)
+    else:
+        fhot = lai * np.sqrt(ko * ks)
+        fint = (1 - np.exp(-alf)) * 0.05
+        sumint = 0
+        x1 = y1 = 0
+        f1 = 1
+        for i in range(20):
+            x2 = -np.log(1 - i * fint) / alf if i < 19 else 1
+            y2 = -(ko + ks) * lai * x2 + fhot * (1 - np.exp(-alf * x2)) / alf
+            f2 = np.exp(y2)
+            if y2 != y1:
+                sumint += (f2 - f1) * (x2 - x1) / (y2 - y1)
+            x1, y1, f1 = x2, y2, f2
+        tsstoo = f1
+
+    rsos = w * lai * sumint
+    rso = rsos + rsod
+
+    # Canopy–soil interaction
+    dn = 1 - rsoil * rdd
+    rddt = rdd + tdd * rsoil * tdd / dn
+    rsdt = rsd + (tsd + tss) * rsoil * tdd / dn
+    rdot = rdo + tdd * rsoil * (tdo + too) / dn
+    rsodt = rsod + ((tss + tsd) * tdo + (tsd + tss * rsoil * rdd) * too) * rsoil / dn
+    rsost = rsos + tsstoo * rsoil
+    rsot = rsost + rsodt
+
+    return np.column_stack((wl, rdot, rsot, rddt, rsdt))

prosail_method/modules/tmpt93j6gvk

@@ -0,0 +1,23 @@
+from modules.prospect_d import prospect_d
+from modules.sail import sail
+
+def prosail(N, Cab, Car, Ant, Cbrown, Cw, Cm,
+            LIDFa, LIDFb, TypeLidf, lai, hspot,
+            tts, tto, psi, rsoil,
+            wl, nr, Kab, Kcar, Kant, Kbrown, Kw, Km):
+    """
+    PROSAIL simulator using PROSPECT-D + 4SAIL.
+
+    Returns:
+        refl: canopy directional reflectance
+        LRT : leaf reflectance/transmittance spectra
+    """
+    # Step 1: simulate leaf optical properties using PROSPECT-D
+    LRT = prospect_d(N, Cab, Car, Ant, Cbrown, Cw, Cm, wl, nr, Kab, Kcar, Kant, Kbrown, Kw, Km)
+    rho = LRT[:, 1]
+    tau = LRT[:, 2]
+
+    # Step 2: simulate canopy reflectance using 4SAIL
+    refl = sail(rho, tau, LIDFa, LIDFb, TypeLidf, lai, hspot, tts, tto, psi, rsoil,wl)
+
+    return wl, refl, LRT

prosail_method/modules/volscatt.py

@@ -0,0 +1,96 @@
+import numpy as np
+
+def volscatt(tts, tto, psi, ttl):
+    """
+    Volume scattering and interception functions for a given geometry.
+
+    Inputs:
+        tts  : solar zenith angle (deg)
+        tto  : observer zenith angle (deg)
+        psi  : relative azimuth angle (deg)
+        ttl  : leaf inclination angle (deg)
+
+    Returns:
+        chi_s : solar interception function
+        chi_o : observer interception function
+        frho  : function to multiply with leaf reflectance
+        ftau  : function to multiply with leaf transmittance
+    """
+    rd = np.pi / 180.0
+
+    costs = np.cos(rd * tts)
+    costo = np.cos(rd * tto)
+    sints = np.sin(rd * tts)
+    sinto = np.sin(rd * tto)
+    cospsi = np.cos(rd * psi)
+    psir = rd * psi
+    costl = np.cos(rd * ttl)
+    sintl = np.sin(rd * ttl)
+
+    cs = costl * costs
+    co = costl * costo
+    ss = sintl * sints
+    so = sintl * sinto
+
+    # Solar side
+    if abs(ss) > 1e-6:
+        cosbts = -cs / ss
+    else:
+        cosbts = 5
+
+    if abs(so) > 1e-6:
+        cosbto = -co / so
+    else:
+        cosbto = 5
+
+    if abs(cosbts) < 1:
+        bts = np.arccos(cosbts)
+        ds = ss
+    else:
+        bts = np.pi
+        ds = cs
+
+    chi_s = 2 / np.pi * ((bts - 0.5 * np.pi) * cs + np.sin(bts) * ss)
+
+    if abs(cosbto) < 1:
+        bto = np.arccos(cosbto)
+        doo = so
+    elif tto < 90:
+        bto = np.pi
+        doo = co
+    else:
+        bto = 0
+        doo = -co
+
+    chi_o = 2 / np.pi * ((bto - 0.5 * np.pi) * co + np.sin(bto) * so)
+
+    # Bidirectional scattering coefficient
+    btran1 = abs(bts - bto)
+    btran2 = np.pi - abs(bts + bto - np.pi)
+
+    if psir <= btran1:
+        bt1 = psir
+        bt2 = btran1
+        bt3 = btran2
+    else:
+        bt1 = btran1
+        if psir <= btran2:
+            bt2 = psir
+            bt3 = btran2
+        else:
+            bt2 = btran2
+            bt3 = psir
+
+    t1 = 2 * cs * co + ss * so * cospsi
+    t2 = 0
+    if bt2 > 0:
+        t2 = np.sin(bt2) * (2 * ds * doo + ss * so * np.cos(bt1) * np.cos(bt3))
+
+    denom = 2 * np.pi * np.pi
+    frho = ((np.pi - bt2) * t1 + t2) / denom
+    ftau = (-bt2 * t1 + t2) / denom
+
+    frho = max(frho, 0)
+    ftau = max(ftau, 0)
+
+    return chi_s, chi_o, frho, ftau