//------------------------------------------------------------------
//
//------------------------------------------------------------------
#property indicator_separate_window
#property indicator_buffers 7
#property indicator_color1  clrDimGray
#property indicator_color2  clrDarkGray
#property indicator_color3  clrRed
#property indicator_color4  clrLimeGreen
#property indicator_color5  clrSandyBrown
#property indicator_color6  clrLimeGreen
#property indicator_color7  clrSandyBrown
#property indicator_width2  2
#property indicator_width3  2
#property strict



extern int                period   = 20;          // RSI period
extern ENUM_APPLIED_PRICE Price    = PRICE_CLOSE; // RSI price
extern double             Factor   = 1;           // Bands deviations
extern int                LastBar  = 0;           // Last bar in the past data
extern int                PastBars = 200;         // Number of past bars
extern int                FutBars  = 50;          // Number of bars to predict 
extern int                HarmNo   = 8;           // Number of frequencies; HarmNo=0 computes PastBars harmonics
extern double             FreqTOL  = 0.0001;      // Tolerance of frequency calculation
extern bool               EveryTick = false;       // Calculate on every tick?






double in[],pv[],fv[],upp[],dnp[],upf[],dnf[];

//------------------------------------------------------------------
//
//------------------------------------------------------------------
//
//
//
//
//

int init()
{
   if(HarmNo==0 || HarmNo>PastBars) HarmNo=PastBars; LastBar=MathMax(LastBar,0);
      SetIndexBuffer(0,in);
      SetIndexBuffer(1,pv);  SetIndexShift(1,       -LastBar);                        // past data vector 0..np-1; 0 corresponds to bar=LastBar
      SetIndexBuffer(2,fv);  SetIndexShift(2,FutBars-LastBar);                        // future data vector i=0..nf; nf corresponds to bar=LastBar
      SetIndexBuffer(3,upp); SetIndexShift(3,       -LastBar); SetIndexLabel(3,NULL); // past data vector
      SetIndexBuffer(4,dnp); SetIndexShift(4,       -LastBar); SetIndexLabel(4,NULL); // past data vector
      SetIndexBuffer(5,upf); SetIndexShift(5,FutBars-LastBar); SetIndexLabel(5,NULL); // future data vector
      SetIndexBuffer(6,dnf); SetIndexShift(6,FutBars-LastBar); SetIndexLabel(6,NULL); // future data vector
   IndicatorShortName("Cycle Extrapolator ("+(string)period+","+(string)PastBars+","+(string)HarmNo+")");  
   return(0);
}
int deinit(){ return(0); }

//------------------------------------------------------------------
//
//------------------------------------------------------------------
//
//
//
//
//

double x[],xp[],xf[];
int start()
{
   static datetime lastTime=0;
   int counted_bars=IndicatorCounted();
      if(counted_bars<0) return(-1);
      if(counted_bars>0) counted_bars--;
         int limit=MathMin(Bars-counted_bars,Bars-1);
         if (lastTime!=0 && (lastTime==Time[0] && !EveryTick)) return(0); 
                             lastTime =Time[0];

         //
         //
         //
         //
         //

         if (Bars<PastBars) return(0);
         int pastBars=MathMin(PastBars-LastBar,Bars);
         if (ArraySize(x)!=pastBars) 
         { 
             ArrayResize(x ,pastBars); 
             ArrayResize(xp,pastBars); 
         }
         int futBars=MathMin(pastBars-1,FutBars);
         if (ArraySize(xf)!=futBars+1) 
             ArrayResize(xf,futBars+1); 

         //
         //
         //
         //
         //
           
         double av=0;
             for(int i=limit;                              i>=0;       i--) in[i]=iRSI(NULL,0,period,Price,i)-50; 
             for(int i=MathMin(pastBars+LastBar-1,Bars-1); i>=LastBar; i--) av+=in[i]; av/=pastBars; av/=pastBars;
               ArrayInitialize(xp,av);
               ArrayInitialize(xf,av);
                  ArrayCopy(x,in,0,LastBar,pastBars);
   
   //
   //
   //
   //
   //
   
   double w,m,c,s,dev=0;
   for(int harm=1;harm<=HarmNo && !IsStopped();harm++)
   {
      Freq(pastBars,FreqTOL,x,xp,w,m,c,s);
      for(int i=0; i<pastBars; i++) 
      {
                           xp[i] += m+c*cos(w*i)+s*sin(w*i);
         if (i<=futBars)   xf[i] += m+c*cos(w*i)-s*sin(w*i);
         if (harm==HarmNo) dev   += MathPow(x[i]-xp[i],2);
      }       
   }
   dev = MathSqrt(dev/(pastBars+1))*Factor;
   ArrayCopy(pv,xp); ArrayCopy(fv,xf); 
   for(int i=0;i<pastBars;i++) 
   {
      upp[i] = pv[i] + dev;
      dnp[i] = pv[i] - dev;
      upf[i] = (i<=futBars) ? fv[i] + dev : EMPTY_VALUE;
      dnf[i] = (i<+futBars) ? fv[i] - dev : EMPTY_VALUE;
   }      
   #define _invert(_arr,_1,_2) tmp=_arr[_1]; _arr[_1]=_arr[_2-_1]; _arr[_2-_1]=tmp;
   for(int i=0;i<=(futBars-1)/2;i++)
   {
      double tmp;
          _invert(fv ,i,FutBars);
          _invert(upf,i,FutBars);
          _invert(dnf,i,FutBars);
   } 

   //
   //
   //
   //
   //

   SetIndexDrawBegin(1,Bars-PastBars);
   SetIndexDrawBegin(3,Bars-PastBars);
   SetIndexDrawBegin(4,Bars-PastBars);
   SetIndexDrawBegin(2,Bars-FutBars-1);
   SetIndexDrawBegin(5,Bars-FutBars-1);
   SetIndexDrawBegin(6,Bars-FutBars-1);
   return(0); 
}

//------------------------------------------------------------------
//
//    Quinn and Fernandes algorithm
//
//------------------------------------------------------------------
//
//
//
//
//
//

double _z[];
void Freq(int np, double freqTol, double& _x[], double& values[], double& w, double& m, double& c, double& s)
{
   if (ArraySize(_z)!=np) ArrayResize(_z,np);
   
   //
   //
   //
   //
   //
   
   double a=0.0, b=2.0; _z[0]=_x[0]-values[0];
   while(MathAbs(a-b)>freqTol)
   {
      a=b;
            _z[1]=_x[1]-values[1]+a*_z[0];
               double num=_z[0]*_z[1];
               double den=_z[0]*_z[0];
               for(int i=2;i<np;i++)
               {
                  _z[i] = _x[i]-values[i]+a*_z[i-1]-_z[i-2];
                     num  += _z[i-1]*(_z[i]+_z[i-2]);
                     den  += _z[i-1]*_z[i-1];
               }
               b = (den!=0) ? num/den : 0;
   }
   w=MathArccos(MathMin(MathMax(b/2.0,-1),1));

   //
   //
   //
   //
   //
      
   double Sc=0,Ss=0,Scc=0,Sss=0,Scs=0,Sx=0,Sxc=0,Sxs=0;;
      for(int i=0;i<np;i++)
      {
         double wcos = cos(w*i);
         double wsin = sin(w*i);
                Sc  += wcos;
                Ss  += wsin;
                Scc += wcos*wcos;
                Sss += wsin*wsin;
                Scs += wcos*wsin;
                Sx  += (_x[i]-values[i]);
                Sxc += (_x[i]-values[i])*wcos;
                Sxs += (_x[i]-values[i])*wsin;
      }
      Sc/=np; Ss/=np; Scc/=np; Sss/=np; Scs/=np; Sx/=np; Sxc/=np; Sxs/=np; m=Sx; c=0.0;  s=0.0;
      if(w!=0.0)
      {
         double den = MathPow(Scs-Sc*Ss,2)-(Scc-Sc*Sc)*(Sss-Ss*Ss);
                c = ((Sxs-Sx*Ss)*(Scs-Sc*Ss)-(Sxc-Sx*Sc)*(Sss-Ss*Ss))/den;
                s = ((Sxc-Sx*Sc)*(Scs-Sc*Ss)-(Sxs-Sx*Ss)*(Scc-Sc*Sc))/den;
                m = Sx-c*Sc-s*Ss;
      }
   return;
}