Ehlers Phasor Analysis (PHASOR)

// The MIT License (MIT)
// © mihakralj
//@version=6
indicator("Ehlers Phasor Analysis (PHASOR)", shorttitle="PHASOR", overlay=false)

//@function Calculates the Ehlers Phasor Angle, Derived Period, and Trend State.
//@param src The source series to analyze.
//@param period The fixed cycle period to correlate against. Default is 28.
//@returns A tuple: `[float finalPhasorAngle, float derivedPeriod, int trendState]`.
phasor(series float src, simple int period = 28) =>
    float sx_corr = 0.0
    float sy_cos_corr = 0.0
    float sxx_corr = 0.0
    float sxy_cos_corr = 0.0
    float syy_cos_corr = 0.0
    for i = 0 to period - 1
        float x_val = nz(src[i])
        float y_val_cos = math.cos(2 * math.pi * i / period)
        sx_corr += x_val
        sy_cos_corr += y_val_cos
        sxx_corr += x_val * x_val
        sxy_cos_corr += x_val * y_val_cos
        syy_cos_corr += y_val_cos * y_val_cos
    float real_part = 0.0
    float den_cos = (period * sxx_corr - sx_corr * sx_corr) * (period * syy_cos_corr - sy_cos_corr * sy_cos_corr)
    if den_cos > 0
        real_part := (period * sxy_cos_corr - sx_corr * sy_cos_corr) / math.sqrt(den_cos)
    sx_corr := 0.0
    sxx_corr := 0.0
    float sy_sin_corr = 0.0
    float sxy_sin_corr = 0.0
    float syy_sin_corr = 0.0
    for i = 0 to period - 1
        float x_val = nz(src[i])
        float y_val_sin = -math.sin(2 * math.pi * i / period) // Negative sine as per Ehlers
        sx_corr += x_val
        sxx_corr += x_val * x_val
        sy_sin_corr += y_val_sin
        sxy_sin_corr += x_val * y_val_sin
        syy_sin_corr += y_val_sin * y_val_sin
    float imag_part = 0.0
    float den_sin = (period * sxx_corr - sx_corr * sx_corr) * (period * syy_sin_corr - sy_sin_corr * sy_sin_corr)
    if den_sin > 0
        imag_part := (period * sxy_sin_corr - sx_corr * sy_sin_corr) / math.sqrt(den_sin)
    float current_raw_phase = 0.0
    if real_part != 0.0
        current_raw_phase := 90.0 - math.atan(imag_part / real_part) * 180.0 / math.pi
        if real_part < 0.0
            current_raw_phase -= 180.0
    else if imag_part != 0.0
        current_raw_phase := imag_part > 0.0 ? 0.0 : 180.0
    var float core_Phasor_unwrapped_state = na
    if not na(core_Phasor_unwrapped_state[1])
        float diff = current_raw_phase - core_Phasor_unwrapped_state[1]
        if diff > 180.0
            current_raw_phase -= 360.0
        else if diff < -180.0
            current_raw_phase += 360.0
    core_Phasor_unwrapped_state := na(core_Phasor_unwrapped_state[1]) ? current_raw_phase : core_Phasor_unwrapped_state[1] + (current_raw_phase - core_Phasor_unwrapped_state[1])
    float calculated_Phasor_val = core_Phasor_unwrapped_state
    var float final_Phasor_state = na
    if na(final_Phasor_state[1])
        final_Phasor_state := calculated_Phasor_val
    else
        if calculated_Phasor_val < final_Phasor_state[1] and ((calculated_Phasor_val > -135 and final_Phasor_state[1] < 135) or (calculated_Phasor_val < -90 and final_Phasor_state[1] < -90))
            final_Phasor_state := final_Phasor_state[1]
        else
            final_Phasor_state := calculated_Phasor_val
    var float derivedPeriod_calc_state = na 
    float angle_Change_For_Period = final_Phasor_state - nz(final_Phasor_state[1], final_Phasor_state)
    if nz(angle_Change_For_Period) == 0 and not na(derivedPeriod_calc_state[1])
        if derivedPeriod_calc_state[1] != 0
            angle_Change_For_Period := 360.0 / derivedPeriod_calc_state[1]
        else
            angle_Change_For_Period := 0.0
    if nz(angle_Change_For_Period) <= 0 and not na(derivedPeriod_calc_state[1])
        if derivedPeriod_calc_state[1] != 0
            angle_Change_For_Period := 360.0 / derivedPeriod_calc_state[1]
        else
            angle_Change_For_Period := 0.0
    if nz(angle_Change_For_Period) != 0.0
        derivedPeriod_calc_state := 360.0 / angle_Change_For_Period
    else if not na(derivedPeriod_calc_state[1])
        derivedPeriod_calc_state := derivedPeriod_calc_state[1]
    else
        derivedPeriod_calc_state := 60.0 
    derivedPeriod_calc_state := math.max(1.0, math.min(derivedPeriod_calc_state, 60.0))
    var int trendState_calc_state = 0 
    float angle_Change_For_State = final_Phasor_state - nz(final_Phasor_state[1], final_Phasor_state)
    int currentTrendState_calc = 0 
    if angle_Change_For_State <= 6.0 
        if final_Phasor_state >= 90.0 or final_Phasor_state <= -90.0
            currentTrendState_calc := 1 
        else if final_Phasor_state > -90.0 and final_Phasor_state < 90.0
            currentTrendState_calc := -1 
    trendState_calc_state := currentTrendState_calc
    [final_Phasor_state, derivedPeriod_calc_state, trendState_calc_state]

// ---------- Inputs ----------
i_period = input.int(28, "Period", minval=1, group="Phasor Settings")
i_source = input.source(close, "Source", group="Phasor Settings")
showDerivedPeriod = input.bool(false, "Show Derived Period", group="Optional Plots", inline="derived_period")
showTrendState = input.bool(false, "Show Trend State Variable", group="Optional Plots", inline="trend_state")

// ---------- Calculations ----------
// Call the main function to get all values
[phasorAngle, derivedPeriodValue, trendStateValue] = phasor(i_source, i_period)

// ---------- Plotting Phasor Angle ----------
plot(phasorAngle, "Phasor Angle", color=color.yellow, linewidth=2)


// ---------- Optional Plots ----------
// Plot for Derived Period
plot(showDerivedPeriod ? derivedPeriodValue : na, "Derived Period", color=color.yellow, linewidth=2)

// Plot for Trend State
plot(showTrendState ? trendStateValue : na, "Trend State", color=color.yellow, linewidth=2, style=plot.style_histogram)