#include "sierrachart.h"
float four1(SCStudyGraphRef sc, SCStudyGraphRef data, SCStudyGraphRef Out, int nn);

SCDLLName("FFT2")
SCSFExport scsf_FFT2(SCStudyGraphRef sc)

{
	SCSubgraphRef FFT		 	= sc.Subgraph[0];
	SCInputRef INperiod 		= sc.Input[0];
	
	
	if (sc.SetDefaults)
	{
		sc.GraphName = "Fast Fourier Transform";

		sc.AutoLoop = true;

		sc.GraphRegion = 0;

		FFT.Name = "FFT";
		FFT.PrimaryColor = COLOR_CYAN;
		FFT.SecondaryColor = COLOR_MAGENTA;
		FFT.SecondaryColorUsed = 3;

        INperiod.Name="Period Length";
        INperiod.SetInt(55);
        INperiod.SetIntLimits(1,10000);			

		// During development set this flag to 1, so the DLL can be modified. When development is done, set it to 0 to improve performance.
		sc.FreeDLL = 1;

		//sc.CalculationPrecedence = LOW_PREC_LEVEL;
		return;
	}
	
	long unsigned nn 				= INperiod.GetInt();

	SCFloatArrayRef Out				= sc.Subgraph[0].Arrays[1];	
	SCFloatArrayRef data			= sc.Subgraph[0].Arrays[2];		

	
	if(sc.GetBarHasClosedStatus() == BHCS_BAR_HAS_CLOSED)
	{
		data[sc.Index] = sc.Close[sc.Index];
		FFT[sc.Index]=four1(sc, data, Out, nn);
	}	
}
float four1(SCStudyGraphRef sc, SCFloatArrayRef data, SCFloatArrayRef Out, long unsigned nn)
//void four1(SCStudyGraphRef sc, SCFloatArrayRef In, SCFloatArrayRef Out, double* data, unsigned long nn)
{
    unsigned long n, mmax, m, j, istep, i;
    double wtemp, wr, wpr, wpi, wi, theta;
    double tempr, tempi;

    // reverse-binary reindexing
    n = nn<<1;
    j=1;
    for (i=1; i<n; i+=2) {
        if (j>i) {
            swap(data[j-1], data[i-1]);
            swap(data[j], data[i]);
        }
        m = nn;
        while (m>=2 && j>m) {
            j -= m;
            m >>= 1;
        }
        j += m;
    };

    // here begins the Danielson-Lanczos section
    mmax=2;
    while (n>mmax) {
        istep = mmax<<1;
        theta = -(2*M_PI/mmax);
        wtemp = sin(0.5*theta);
        wpr = -2.0*wtemp*wtemp;
        wpi = sin(theta);
        wr = 1.0;
        wi = 0.0;
        for (m=1; m < mmax; m += 2) {
            for (i=m; i <= n; i += istep) {
                j=i+mmax;
                tempr = wr*data[j-1] - wi*data[j];
                tempi = wr * data[j] + wi*data[j-1];

                data[j-1] = data[i-1] - tempr;
                data[j] = data[i] - tempi;
                data[i-1] += tempr;
                data[i] += tempi;
            }
            wtemp=wr;
            wr += wr*wpr - wi*wpi;
            wi += wi*wpr + wtemp*wpi;
        }
        mmax=istep;
    }
	Out[i]=data[i];
	return (Out[i]);
}