//+------------------------------------------------------------------+
//|                                              Wave-Pm Smooth.mq4  |
//|                                                                  |
//|                                                                  |
//+------------------------------------------------------------------+
#property copyright ""
#property link      ""


#property indicator_separate_window
#property indicator_buffers    2
#property indicator_maximum    1
#property indicator_minimum    0.2
#property indicator_color1     LimeGreen
#property indicator_color2     Red
#property indicator_width1     1
#property indicator_width2     2
#property indicator_level1     0.30
#property indicator_level2     0.50
#property indicator_level3     0.70
#property indicator_level4     0.90
#property indicator_levelcolor MediumOrchid

//
//
//
//
//

extern string note_TimeFrames        = "M1;5,15,30,60H1;240H4;1440D1;10080W1;43200MN|0-CurrentTF";
extern string TimeFrame              = "Current time frame";
extern int    LongBandsPeriod        = 55;
extern int    LongBandsMode          = 0;
extern int    LongBandsPrice         = 0;
extern double LongBandsDeviations    = 2.2;
extern double LongSmoothLength       = 10.0;
extern double LongSmoothPhase        = 0.0;
extern bool   LongSmoothDouble       = false;
extern int    ShortBandsPeriod       = 14;
extern int    ShortBandsMode         = 0;
extern int    ShortBandsPrice        = 0;
extern double ShortBandsDeviations   = 2.2;
extern double ShortSmoothLength      = 10.0;
extern double ShortSmoothPhase       = 0.0;
extern bool   ShortSmoothDouble      = false;
extern int    periods_characteristic = 100;
extern bool   Interpolate            = true;

extern bool   alertsOn               = true;
extern bool   alertsOnSlowSlope      = true;
extern bool   alertsOnSlow50Cross    = true;
extern bool   alertsOnSlow90Cross    = true;
extern bool   alertsOnFastSlope      = true;
extern bool   alertsOnFast50Cross    = true;
extern bool   alertsOnFast90Cross    = true;
extern bool   alertsOnFastSlowSlope  = true;
extern bool   alertsOnCurrent        = false;
extern bool   alertsMessage          = true;
extern bool   alertsSound            = true;
extern bool   alertsEmail            = false;

extern string  __                    = "arrows settings";
extern bool   ShowArrows             = true;
extern string arrowsIdentifier       = "wave-pm arrows";
extern color  arrowsUpColor          = Lime;
extern color  arrowsDnColor          = Red;

extern bool   verticalLinesVisible   = true;
extern string verticalLinesID        = "wave-pm lines";
extern color  verticalLinesUpColor   = DeepSkyBlue;
extern color  verticalLinesDownColor = PaleVioletRed;
extern int    verticalLinesStyle     = STYLE_DOT;
extern int    verticalLinesWidth     = 0;

//
//
//
//
//

double Shortoscillator[];
double Longoscillator[];
double ShortDev[];
double LongDev[];
double LongBollMd[];
double ShortBollMd[];
double slope[];
double trends[][7];

//
//
//
//
//

#define _tup1 0
#define _tup2 1
#define _tup3 2
#define _tdn1 3
#define _tdn2 4
#define _tdn3 5
#define _tmi  6

//
//
//
//
//

string indicatorFileName;
bool   calculateValue;
bool   returnBars;
int    timeFrame;

//
//
//
//
//

int init()
  {
   IndicatorBuffers(7);   
   SetIndexBuffer(0,Shortoscillator);
   SetIndexBuffer(1,Longoscillator);
   SetIndexBuffer(2,ShortDev); 
   SetIndexBuffer(3,LongDev); 
   SetIndexBuffer(4,LongBollMd); 
   SetIndexBuffer(5,ShortBollMd);
   SetIndexBuffer(6,slope);
   
   //
   //
   //
   //
   //
   
   indicatorFileName = WindowExpertName();
   calculateValue    = (TimeFrame=="calculateValue"); if (calculateValue) return(0);
   returnBars        = (TimeFrame=="returnBars");     if (returnBars)     return(0);
   timeFrame         = stringToTimeFrame(TimeFrame);
   
   //
   //
   //
   //
   //
   
   IndicatorShortName(timeFrameToString(timeFrame)+  "  wave-pm smooth ");
   
  return(0);
  }
  
//
//
//
//
//

int deinit()
{

   if (!calculateValue && ShowArrows) deleteArrows();
   
   string lookFor       = verticalLinesID+":";
   int    lookForLength = StringLen(lookFor);
   for (int i=ObjectsTotal()-1; i>=0; i--)
   {
      string objectName = ObjectName(i);
         if (StringSubstr(objectName,0,lookForLength) == lookFor) ObjectDelete(objectName);
   }
   return(0);
}

//+------------------------------------------------------------------+
//|                                                                  |
//+------------------------------------------------------------------+
//
//
//
//

int start()
  {
   int counted_bars=IndicatorCounted();
   int i,r,j,limit;

   if(counted_bars<0) return(-1);
   if(counted_bars>0) counted_bars--;
         limit=MathMin(Bars-1,Bars-counted_bars-1);
         if (returnBars) { Shortoscillator[0] = limit+1; return(0); }
           
    //
    //
    //
    //
    //
           
    if (calculateValue || timeFrame == Period())
    {
    
    
    if (ArrayRange(trends,0)!=Bars) {  ArrayResize(trends,Bars); }
   
   
    for(i = limit, r=Bars-i-1; i>=0; i--,r++)
     {
      
      double Longsum  = 0.0;
      double Shortsum = 0.0;
      LongBollMd[i]   = iMA(NULL,0,LongBandsPeriod, 0,LongBandsMode, LongBandsPrice, i);
      ShortBollMd[i]  = iMA(NULL,0,ShortBandsPeriod,0,ShortBandsMode,ShortBandsPrice,i);
      
      //
      //
      //
      //
      //
      
      for(j=0; j<LongBandsPeriod; j++)
      {
      double Longpr  = iMA(NULL,0,1,0,MODE_SMA,LongBandsPrice,i+j);
             Longsum+=(Longpr-LongBollMd[i])*(Longpr-LongBollMd[i]);
      }
      
      //
      //
      //
      //
      //
      
      for(j=0; j<ShortBandsPeriod; j++)
      {
      double Shortpr  = iMA(NULL,0,1,0,MODE_SMA,ShortBandsPrice,i+j);
             Shortsum+=(Shortpr-ShortBollMd[i])*(Shortpr-ShortBollMd[i]);
      }
      
      //
      //
      //
      //
      //
      
      LongDev[i]         = MathSqrt(Longsum/LongBandsPeriod);
      Longoscillator[i]  = iDSmooth(OscillatorLine(LongDev,i), LongSmoothLength, LongSmoothPhase, LongSmoothDouble,i,  0);
      ShortDev[i]        = MathSqrt(Shortsum/ShortBandsPeriod);
      Shortoscillator[i] = iDSmooth(OscillatorLine(ShortDev,i),ShortSmoothLength,ShortSmoothPhase,ShortSmoothDouble,i,20);
                slope[i] = slope[i+1];
                
      if(Longoscillator[i] > Longoscillator[i+1] && Shortoscillator[i] > Shortoscillator[i+1]) slope[i] =  1;  
      if(Longoscillator[i] < Longoscillator[i+1] && Shortoscillator[i] < Shortoscillator[i+1]) slope[i] = -1;  
                
      if (!calculateValue) manageArrow(i);          
                
      
      manageLines(i);
           
      setTrends(i,r);
       
     }
     
   manageAlerts();
   return(0);
  }
  
  //
  //
  //
  //
  //

  limit = MathMax(limit,MathMin(Bars,iCustom(NULL,timeFrame,indicatorFileName,"returnBars",0,0)*timeFrame/Period()));
   
   for(i = limit, r=Bars-i-1; i>=0; i--,r++)
   {
   int y = iBarShift(NULL,timeFrame,Time[i]);
   Shortoscillator[i] = iCustom(NULL,timeFrame,indicatorFileName,"","calculateValue",LongBandsPeriod,LongBandsMode,LongBandsPrice,LongBandsDeviations,LongSmoothLength,LongSmoothPhase,LongSmoothDouble,ShortBandsPeriod,ShortBandsMode,ShortBandsPrice,ShortBandsDeviations,ShortSmoothLength,ShortSmoothPhase,ShortSmoothDouble,periods_characteristic,0,y);
   Longoscillator[i]  = iCustom(NULL,timeFrame,indicatorFileName,"","calculateValue",LongBandsPeriod,LongBandsMode,LongBandsPrice,LongBandsDeviations,LongSmoothLength,LongSmoothPhase,LongSmoothDouble,ShortBandsPeriod,ShortBandsMode,ShortBandsPrice,ShortBandsDeviations,ShortSmoothLength,ShortSmoothPhase,ShortSmoothDouble,periods_characteristic,1,y);
   slope[i]           = iCustom(NULL,timeFrame,indicatorFileName,"","calculateValue",LongBandsPeriod,LongBandsMode,LongBandsPrice,LongBandsDeviations,LongSmoothLength,LongSmoothPhase,LongSmoothDouble,ShortBandsPeriod,ShortBandsMode,ShortBandsPrice,ShortBandsDeviations,ShortSmoothLength,ShortSmoothPhase,ShortSmoothDouble,periods_characteristic,6,y);
   
   manageLines(i);
   
   manageArrow(i);
   
   setTrends(i,r);
   
   //
   //
   //
   //
   //
   
   if (timeFrame <= Period() || y==iBarShift(NULL,timeFrame,Time[i-1])) continue;
   if (!Interpolate) continue;
   
   //
   //
   //
   //
   //
         
   datetime time = iTime(NULL,timeFrame,y);
      for(int n = 1; i+n < Bars && Time[i+n] >= time; n++) continue;	
      double factor = 1.0 / n;
      for(j = 1; j < n; j++)
      {
  	   Shortoscillator[i+j] = j*factor*Shortoscillator[i+n] + (1.0-j*factor)* Shortoscillator[i];
  	   Longoscillator[i+j]  = j*factor*Longoscillator[i+n]  + (1.0-j*factor)* Longoscillator[i];

      }
   }

   //
   //
   //
   //
   //
   manageAlerts(); 
   return(0);
}

//+------------------------------------------------------------------+
//|                                                                  |
//+------------------------------------------------------------------+


void setTrends(int i, int r)
{
   
   trends[r][_tup1] = trends[r-1][_tup1];
   trends[r][_tup2] = trends[r-1][_tup2];
   trends[r][_tup3] = trends[r-1][_tup3];
   trends[r][_tdn1] = trends[r-1][_tdn1];
   trends[r][_tdn2] = trends[r-1][_tdn2];
   trends[r][_tdn3] = trends[r-1][_tdn3];
   trends[r][_tmi]  = trends[r-1][_tmi];

      if (Longoscillator[i]  > Longoscillator[i+1])                                               trends[r][_tup1] =  1;
      if (Longoscillator[i]  < Longoscillator[i+1])                                               trends[r][_tup1] = -1;
      if (Longoscillator[i]  > Longoscillator[i+1]  && Longoscillator[i]  < 0.50)                 trends[r][_tup2] =  1;
      if (Longoscillator[i]  < Longoscillator[i+1]  && Longoscillator[i]  < 0.50)                 trends[r][_tup2] = -1;
      if (Longoscillator[i]  > Longoscillator[i+1]  && Longoscillator[i]  > 0.90)                 trends[r][_tup3] =  1;
      if (Longoscillator[i]  < Longoscillator[i+1]  && Longoscillator[i]  > 0.90)                 trends[r][_tup3] = -1;
      if (Shortoscillator[i] > Shortoscillator[i+1])                                              trends[r][_tdn1] =  1;
      if (Shortoscillator[i] < Shortoscillator[i+1])                                              trends[r][_tdn1] = -1;
      if (Shortoscillator[i] > Shortoscillator[i+1] && Shortoscillator[i] < 0.50)                 trends[r][_tdn2] =  1;
      if (Shortoscillator[i] < Shortoscillator[i+1] && Shortoscillator[i] < 0.50)                 trends[r][_tdn2] = -1;
      if (Shortoscillator[i] > Shortoscillator[i+1] && Shortoscillator[i] > 0.90)                 trends[r][_tdn3] =  1;
      if (Shortoscillator[i] < Shortoscillator[i+1] && Shortoscillator[i] > 0.90)                 trends[r][_tdn3] = -1;
      if (Longoscillator[i]  > Longoscillator[i+1]  && Shortoscillator[i] > Shortoscillator[i+1]) trends[r][_tmi]  =  1;
      if (Longoscillator[i]  < Longoscillator[i+1]  && Shortoscillator[i] < Shortoscillator[i+1]) trends[r][_tmi]  = -1;
                           
}

//
//
//
//
//

void manageAlerts()
{
   if (!calculateValue && alertsOn)
   {
      if (alertsOnCurrent)
           int whichBar = 0;
      else     whichBar = 1; whichBar = iBarShift(NULL,0,iTime(NULL,timeFrame,whichBar)); 
                             whichBar = Bars-whichBar-1;

      //
      //
      //
      //
      //
            
      static datetime time1 = 0;
      static string   mess1 = "";
      if (alertsOnSlowSlope && trends[whichBar][_tup1] != trends[whichBar-1][_tup1])
      {
         if (trends[whichBar][_tup1] ==  1) doAlert(time1,mess1,whichBar,"slow says trend");
         if (trends[whichBar][_tup1] == -1) doAlert(time1,mess1,whichBar,"slow says no trend");
      }
      static datetime time2 = 0;
      static string   mess2 = "";
      if (alertsOnSlow50Cross && trends[whichBar][_tup2] != trends[whichBar-1][_tup2])
      {
         if (trends[whichBar][_tup2] ==  1) doAlert(time2,mess2,whichBar,"slow says trend and less than 50");
         if (trends[whichBar][_tup2] == -1) doAlert(time2,mess2,whichBar,"slow says no trend but less than 50");
      }
      static datetime time3 = 0;
      static string   mess3 = "";
      if (alertsOnSlow90Cross && trends[whichBar][_tup3] != trends[whichBar-1][_tup3])
      {
         if (trends[whichBar][_tup3] ==  1) doAlert(time3,mess3,whichBar,"slow says trend but more than 90");
         if (trends[whichBar][_tup3] == -1) doAlert(time3,mess3,whichBar,"slow says no trend and more than 90");
      }
      static datetime time4 = 0;
      static string   mess4 = "";
      if (alertsOnFastSlope && trends[whichBar][_tdn1] != trends[whichBar-1][_tdn1])
      {
         if (trends[whichBar][_tdn1] ==  1) doAlert(time4,mess4,whichBar,"fast says trend");
         if (trends[whichBar][_tdn1] == -1) doAlert(time4,mess4,whichBar,"fast says no trend");
      }
      static datetime time5 = 0;
      static string   mess5 = "";
      if (alertsOnFast50Cross && trends[whichBar][_tdn2] != trends[whichBar-1][_tdn2])
      {
         if (trends[whichBar][_tdn2] ==  1) doAlert(time5,mess5,whichBar,"fast says trend and less than 50");
         if (trends[whichBar][_tdn2] == -1) doAlert(time5,mess5,whichBar,"fast says no trend but less than 50");
      }
      static datetime time6 = 0;
      static string   mess6 = "";
      if (alertsOnFast90Cross && trends[whichBar][_tdn3] != trends[whichBar-1][_tdn3])
      {
         if (trends[whichBar][_tdn3] ==  1) doAlert(time6,mess6,whichBar,"fast says trend but more than 90");
         if (trends[whichBar][_tdn3] == -1) doAlert(time6,mess6,whichBar,"fast says no trend but more than 90");
      }
      static datetime time7 = 0;
      static string   mess7 = "";
      if (alertsOnFastSlowSlope && trends[whichBar][_tmi] != trends[whichBar-1][_tmi])
      {
         if (trends[whichBar][_tmi] ==  1) doAlert(time7,mess7,whichBar,"fast and slow both say trend");
         if (trends[whichBar][_tmi] == -1) doAlert(time7,mess7,whichBar,"fast and slow both say no trend");
      }
   }
}

//
//
//
//
//

void doAlert(datetime& previousTime, string& previousAlert, int forBar, string doWhat)
{
   string message;
   
   if (previousAlert != doWhat || previousTime != Time[forBar]) {
       previousAlert  = doWhat;
       previousTime   = Time[forBar];

       //
       //
       //
       //
       //

       message =  StringConcatenate(Symbol()," ",timeFrameToString(timeFrame)," at ",TimeToStr(TimeLocal(),TIME_SECONDS)," Wave-PM ",doWhat);
          if (alertsMessage) Alert(message);
          if (alertsEmail)   SendMail(StringConcatenate(Symbol()," Wave-PM "),message);
          if (alertsSound)   PlaySound("alert2.wav");
   }
}

//
//
//
//
//

//+-------------------------------------------------------------------
//|                                                                  
//+-------------------------------------------------------------------
//
//
//
//
//

string sTfTable[] = {"M1","M5","M15","M30","H1","H4","D1","W1","MN"};
int    iTfTable[] = {1,5,15,30,60,240,1440,10080,43200};

//
//
//
//
//

int stringToTimeFrame(string tfs)
{
   tfs = stringUpperCase(tfs);
   for (int i=ArraySize(iTfTable)-1; i>=0; i--)
         if (tfs==sTfTable[i] || tfs==""+iTfTable[i]) return(MathMax(iTfTable[i],Period()));
                                                      return(Period());
}
string timeFrameToString(int tf)
{
   for (int i=ArraySize(iTfTable)-1; i>=0; i--) 
         if (tf==iTfTable[i]) return(sTfTable[i]);
                              return("");
}

//
//
//
//
//

string stringUpperCase(string str)
{
   string   s = str;

   for (int length=StringLen(str)-1; length>=0; length--)
   {
      int char = StringGetChar(s, length);
         if((char > 96 && char < 123) || (char > 223 && char < 256))
                     s = StringSetChar(s, length, char - 32);
         else if(char > -33 && char < 0)
                     s = StringSetChar(s, length, char + 224);
   }
   return(s);
} 

//+-------------------------------------------------------------------
//|                                                                  
//+-------------------------------------------------------------------
//
//
//
//
//


void manageArrow(int i)
{
   if (ShowArrows)
   {
      deleteArrow(Time[i]);
      if (slope[i]!=slope[i+1])
      {
         if (slope[i] == 1) drawArrow(i,arrowsUpColor,171,false);
         if (slope[i] ==-1) drawArrow(i,arrowsDnColor,171,true);
      }
   }
}               

//
//
//
//
//

void drawArrow(int i,color theColor,int theCode, bool up)
{
   string name = arrowsIdentifier+":"+Time[i];
   double gap  = iATR(NULL,0,20,i)/2.0;   
   
      //
      //
      //
      //
      //
      
      ObjectCreate(name,OBJ_ARROW,0,Time[i],0);
         ObjectSet(name,OBJPROP_ARROWCODE,theCode);
         ObjectSet(name,OBJPROP_COLOR,theColor);
  
         if (up)
               ObjectSet(name,OBJPROP_PRICE1,High[i]+gap);
         else  ObjectSet(name,OBJPROP_PRICE1,Low[i] -gap);
}

//
//
//
//
//

void deleteArrows()
{
   string lookFor       = arrowsIdentifier+":";
   int    lookForLength = StringLen(lookFor);
   for (int i=ObjectsTotal()-1; i>=0; i--)
   {
      string objectName = ObjectName(i);
         if (StringSubstr(objectName,0,lookForLength) == lookFor) ObjectDelete(objectName);
   }
}
void deleteArrow(datetime time)
{
   string lookFor = arrowsIdentifier+":"+time; ObjectDelete(lookFor);
}

//
//
//
//
//


void manageLines(int i)
{
   if (!calculateValue && verticalLinesVisible )
   {
         deleteLine(Time[i]);
         if (slope[i]!=slope[i+1])
         {
            if (slope[i] == 1) drawLine(i,verticalLinesUpColor);
            if (slope[i] ==-1) drawLine(i,verticalLinesDownColor);
         }
   }
}               

//
//
//
//
//

void drawLine(int i,color theColor)
{
   string name = verticalLinesID+":"+Time[i];
   
      //
      //
      //
      //
      //
      
      ObjectCreate(name,OBJ_VLINE,0,Time[i],0);
         ObjectSet(name,OBJPROP_COLOR,theColor);
         ObjectSet(name,OBJPROP_STYLE,verticalLinesStyle);
         ObjectSet(name,OBJPROP_WIDTH,verticalLinesWidth);
         ObjectSet(name,OBJPROP_BACK,true);
}

//
//
//
//
//

void deleteLine(datetime time)
{
   string lookFor = verticalLinesID+":"+time; ObjectDelete(lookFor);
}

//
//
//
//
//

double OscillatorLine(double indicator[], int StartBar)
  {
      double S=0;
      double Result;
      int ArrayLong=periods_characteristic;
      
      for(int j=StartBar;j<ArrayLong+StartBar;j++)
      {
      
      S += MathPow((indicator[j]/Point),2);
      
      }
      
      S /= ArrayLong;
      S  = MathSqrt(S)*Point;
      
      if(S!=0) { Result=indicator[StartBar]/S; } 
        
      Result=MathTanh(Result);
      return (Result);
  }

//
//
//
//
//

double MathTanh(double x)
{ 
   double exp;
   double returnNum;
   
   if(x>0)
     {
       exp=MathExp(-2*x);
       returnNum= (1-exp)/(1+exp);
       return (returnNum);
     }
   else
     {
       exp=MathExp(2*x);
       returnNum=(exp-1)/(1+exp);
       return (returnNum);
     }
}

//+-------------------------------------------------------------------
//|                                                                  
//+-------------------------------------------------------------------
//
//

double wrk[][40];

#define bsmax  5
#define bsmin  6
#define volty  7
#define vsum   8
#define avolty 9

//
//
//
//
//

double iDSmooth(double price, double length, double phase, bool isDouble, int i, int s=0)
{
   if (isDouble)
         return (iSmooth(iSmooth(price,MathSqrt(length),phase,i,s),MathSqrt(length),phase,i,s+10));
   else  return (iSmooth(price,length,phase,i,s));
}

//
//
//
//
//

double iSmooth(double price, double length, double phase, int i, int s=0)
{
   if (length <=1) return(price);
   if (ArrayRange(wrk,0) != Bars) ArrayResize(wrk,Bars);
   
   int r = Bars-i-1; 
      if (r==0) { for(int k=0; k<7; k++) wrk[r][k+s]=price; for(; k<10; k++) wrk[r][k+s]=0; return(price); }

   //
   //
   //
   //
   //
   
      double len1   = MathMax(MathLog(MathSqrt(0.5*(length-1)))/MathLog(2.0)+2.0,0);
      double pow1   = MathMax(len1-2.0,0.5);
      double del1   = price - wrk[r-1][bsmax+s];
      double del2   = price - wrk[r-1][bsmin+s];
      double div    = 1.0/(10.0+10.0*(MathMin(MathMax(length-10,0),100))/100);
      int    forBar = MathMin(r,10);
	
         wrk[r][volty+s] = 0;
               if(MathAbs(del1) > MathAbs(del2)) wrk[r][volty+s] = MathAbs(del1); 
               if(MathAbs(del1) < MathAbs(del2)) wrk[r][volty+s] = MathAbs(del2); 
         wrk[r][vsum+s] =	wrk[r-1][vsum+s] + (wrk[r][volty+s]-wrk[r-forBar][volty+s])*div;
         
         //
         //
         //
         //
         //
   
         wrk[r][avolty+s] = wrk[r-1][avolty+s]+(2.0/(MathMax(4.0*length,30)+1.0))*(wrk[r][vsum+s]-wrk[r-1][avolty+s]);
            if (wrk[r][avolty+s] > 0)
               double dVolty = wrk[r][volty+s]/wrk[r][avolty+s]; else dVolty = 0;   
	               if (dVolty > MathPow(len1,1.0/pow1)) dVolty = MathPow(len1,1.0/pow1);
                  if (dVolty < 1)                      dVolty = 1.0;

      //
      //
      //
      //
      //
	        
   	double pow2 = MathPow(dVolty, pow1);
      double len2 = MathSqrt(0.5*(length-1))*len1;
      double Kv   = MathPow(len2/(len2+1), MathSqrt(pow2));

         if (del1 > 0) wrk[r][bsmax+s] = price; else wrk[r][bsmax+s] = price - Kv*del1;
         if (del2 < 0) wrk[r][bsmin+s] = price; else wrk[r][bsmin+s] = price - Kv*del2;
	
   //
   //
   //
   //
   //
      
      double R     = MathMax(MathMin(phase,100),-100)/100.0 + 1.5;
      double beta  = 0.45*(length-1)/(0.45*(length-1)+2);
      double alpha = MathPow(beta,pow2);

         wrk[r][0+s] = price + alpha*(wrk[r-1][0+s]-price);
         wrk[r][1+s] = (price - wrk[r][0+s])*(1-beta) + beta*wrk[r-1][1+s];
         wrk[r][2+s] = (wrk[r][0+s] + R*wrk[r][1+s]);
         wrk[r][3+s] = (wrk[r][2+s] - wrk[r-1][4+s])*MathPow((1-alpha),2) + MathPow(alpha,2)*wrk[r-1][3+s];
         wrk[r][4+s] = (wrk[r-1][4+s] + wrk[r][3+s]); 

   //
   //
   //
   //
   //

   return(wrk[r][4+s]);
}