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