// Licensed under the Apache License, Version 2.0
// © mihakralj
//@version=6
indicator("Ehlers Hilbert Transform Phasor Components (HT_PHASOR)", "HT_PHASOR", overlay=false)

//@function Calculates Hilbert Transform Phasor Components using TA-Lib algorithm
//@param source Series to analyze for phasor components
//@returns Tuple [inphase, quadrature] - raw I1[3] and Q1 components
ht_phasor(series float source) =>
    var float detrender = 0.0
    var float i1 = 0.0
    var float q1 = 0.0
    var float ji = 0.0
    var float jq = 0.0
    var float i2 = 0.0
    var float q2 = 0.0
    var float re = 0.0
    var float im = 0.0
    var float period = 0.0
    var float smooth_period = 0.0

    float price = nz(source)

    // Step 1: WMA smoothing (4-tap: [4,3,2,1]/10)
    float smooth_price = (4.0 * price + 3.0 * nz(price[1]) + 2.0 * nz(price[2]) + nz(price[3])) / 10.0

    // Bandwidth uses period (not smooth_period) per TA-Lib
    float bandwidth = 0.075 * period + 0.54

    // Step 2: Hilbert FIR detrender
    detrender := (0.0962 * smooth_price + 0.5769 * nz(smooth_price[2]) - 0.5769 * nz(smooth_price[4]) - 0.0962 * nz(smooth_price[6])) * bandwidth

    // Step 3: Q1 computation (Hilbert FIR on detrender)
    q1 := (0.0962 * detrender + 0.5769 * nz(detrender[2]) - 0.5769 * nz(detrender[4]) - 0.0962 * nz(detrender[6])) * bandwidth

    // Step 4: I1 = detrender delayed 3 bars
    i1 := nz(detrender[3])

    // Step 5: Advance phase via JI/JQ
    ji := (0.0962 * i1 + 0.5769 * nz(i1[2]) - 0.5769 * nz(i1[4]) - 0.0962 * nz(i1[6])) * bandwidth
    jq := (0.0962 * q1 + 0.5769 * nz(q1[2]) - 0.5769 * nz(q1[4]) - 0.0962 * nz(q1[6])) * bandwidth

    // Step 6: Smooth I2/Q2 with 2-bar EMA (used internally for period calc)
    i2 := 0.2 * (i1 - jq) + 0.8 * nz(i2[1])
    q2 := 0.2 * (q1 + ji) + 0.8 * nz(q2[1])

    // Step 7: Homodyne discriminator
    re := 0.2 * (i2 * nz(i2[1]) + q2 * nz(q2[1])) + 0.8 * nz(re[1])
    im := 0.2 * (i2 * nz(q2[1]) - q2 * nz(i2[1])) + 0.8 * nz(im[1])

    // Step 8: Period from atan (NOT atan2) + rate limiting + clamping
    float prev_period = period
    if math.abs(im) > 1e-12 and math.abs(re) > 1e-12
        float angle = math.atan(im / re)
        if math.abs(angle) > 1e-12
            period := 2.0 * math.pi / angle

    // Rate limit: ±50% bar-to-bar
    if prev_period > 0
        period := math.min(period, 1.5 * prev_period)
        period := math.max(period, 0.67 * prev_period)

    // Clamp to valid range
    period := math.max(6.0, math.min(50.0, period))

    // Step 9: Smooth period with 0.2/0.8 EMA
    period := 0.2 * period + 0.8 * prev_period

    // Step 10: Smooth smoothPeriod with 0.33/0.67 EMA
    smooth_period := 0.33 * period + 0.67 * smooth_period

    // Step 11: Output raw I1[3] (inPhase) and Q1 (quadrature) per TA-Lib HT_PHASOR
    // TA-Lib outputs the detrender delayed by 3 bars as InPhase, and the raw Q1 as Quadrature
    float inphase_out = nz(i1[3])
    float quadrature_out = q1
    [inphase_out, quadrature_out]

// ---------- Main loop ----------

// Inputs
i_source = input.source(hlc3, "Source")

// Calculation
[inphase, quadrature] = ht_phasor(i_source)

// Plot
plot(inphase, "InPhase", color=color.yellow, linewidth=2)
plot(quadrature, "Quadrature", color=color.blue, linewidth=2)
hline(0, "Zero", color=color.gray, linestyle=hline.style_solid)