// Schaff Trend Cycle
//
// schaff_trend_cycle
//
// Schaff Trend Cycle
//
// 
//
// Inputs:    1
// Outputs:   2
// Constants: 1


#include <stdlib.h>
#include <math.h>
#include "ts.h"

typedef struct {
    uint_t period;
} ts_schaff_trend_cycle_params_t;

ts_parameter_t const ts_schaff_trend_cycle_p = { .n_inputs = 1, .n_outputs = 2, .n_constants = 1, .params_size = sizeof(ts_schaff_trend_cycle_params_t) };

uint_t ts_schaff_trend_cycle_offset( ts_constant_t const *c )
{
	( void )c;
	return ( 52 );
}

int ts_schaff_trend_cycle(
	fp_t * const *y, uint_t y_n,
	fp_t const * const *x, uint_t x_n,
	uint_t n, ts_constant_t const *c
)
{
	TS_PROLOGUE( ts_schaff_trend_cycle_p );
	if ( !c ) return -1;
    uint_t period = ((ts_schaff_trend_cycle_params_t *)c)->period;
    if (period == 0) period = 3;
    if (n == 0) return -1;
    if (period > n) { for (uint_t j = 0; j < n; ++j) { y[0][j] = 0; y[1][j] = 0; } return 0; }

    uint_t fast = 10;
    uint_t slow = 21;
    uint_t rsi_period = 13;

    fp_t ema_fast = x[0][0];
    fp_t ema_slow = x[0][0];
    fp_t a_fast = 2.0 / (fast + 1);
    fp_t a_slow = 2.0 / (slow + 1);

    fp_t *macd_line = calloc(n, sizeof(fp_t));
    if (!macd_line) return -1;

    for( uint_t i = 1; i < n; ++i )
    {
        ema_fast = ema_fast + a_fast * (x[0][i] - ema_fast);
        ema_slow = ema_slow + a_slow * (x[0][i] - ema_slow);
        macd_line[i] = ema_fast - ema_slow;
    }

    fp_t *rsi1 = calloc(n, sizeof(fp_t));
    fp_t *rsi2 = calloc(n, sizeof(fp_t));
    if (!rsi1 || !rsi2) { free(macd_line); free(rsi1); free(rsi2); return -1; }

    {
        fp_t avg_up = 0, avg_down = 0;
        for( uint_t i = 1; i <= rsi_period; ++i )
        {
            fp_t change = macd_line[i] - macd_line[i - 1];
            if( change > 0 ) avg_up += change;
            else avg_down += (-change);
        }
        avg_up /= rsi_period;
        avg_down /= rsi_period;
        rsi1[rsi_period] = (avg_down == 0) ? 100 : 100 - 100 / (1 + avg_up / avg_down);

        for( uint_t i = rsi_period + 1; i < n; ++i )
        {
            fp_t change = macd_line[i] - macd_line[i - 1];
            fp_t gain = change > 0 ? change : 0;
            fp_t loss = change < 0 ? (-change) : 0;
            avg_up = (avg_up * (rsi_period - 1) + gain) / rsi_period;
            avg_down = (avg_down * (rsi_period - 1) + loss) / rsi_period;
            rsi1[i] = (avg_down == 0) ? 100 : 100 - 100 / (1 + avg_up / avg_down);
        }
    }

    {
        fp_t avg_up = 0, avg_down = 0;
        for( uint_t i = rsi_period + 1; i <= rsi_period * 2; ++i )
        {
            fp_t change = rsi1[i] - rsi1[i - 1];
            if( change > 0 ) avg_up += change;
            else avg_down += (-change);
        }
        avg_up /= rsi_period;
        avg_down /= rsi_period;
        rsi2[rsi_period * 2] = (avg_down == 0) ? 100 : 100 - 100 / (1 + avg_up / avg_down);

        for( uint_t i = rsi_period * 2 + 1; i < n; ++i )
        {
            fp_t change = rsi1[i] - rsi1[i - 1];
            fp_t gain = change > 0 ? change : 0;
            fp_t loss = change < 0 ? (-change) : 0;
            avg_up = (avg_up * (rsi_period - 1) + gain) / rsi_period;
            avg_down = (avg_down * (rsi_period - 1) + loss) / rsi_period;
            rsi2[i] = (avg_down == 0) ? 100 : 100 - 100 / (1 + avg_up / avg_down);
        }
    }

    uint_t offset = slow + rsi_period * 2;
    fp_t sig_a = 1.0 / (fp_t)period;
    fp_t signal_ema = 0;
    uint_t active = 0;
    uint_t count = 0;

    for (uint_t j = 0; j < n && j < offset; ++j) { y[0][j] = 0; y[1][j] = 0; }
    for( uint_t i = offset; i < n; ++i )
    {
        fp_t hist = macd_line[i] - macd_line[i - 1];
        fp_t hist_prev = (i > 1) ? macd_line[i - 1] - macd_line[i - 2] : 0;
        int prev_above = hist_prev > 0.5;
        int curr_above = hist > 0.5;
        fp_t stc = 0;

        if( !prev_above && curr_above )
        {
            active = 1;
            count = 0;
        }
        if( active )
        {
            stc = 100;
            count++;
            if( count >= period )
                active = 0;
        }

        y[0][i] = stc;

        signal_ema = signal_ema + sig_a * (stc - signal_ema);
        y[1][i] = signal_ema;
    }

    free(macd_line);
    free(rsi1);
    free(rsi2);
	return 0;
}