#property indicator_separate_window
#property indicator_buffers 2
#property indicator_color1  clrDodgerBlue
#property indicator_style1  STYLE_SOLID
#property indicator_width1  1
#property indicator_color2  clrDodgerBlue
#property indicator_style2  STYLE_DASH
#property indicator_width2  1


#include <SWISSEPH.mqh>


input BODY HeavenlyBody=Moon;
input int ExtrapolatedBars=200;

//--- indicator buffers
double Result[], Extrapolation[];

int flag=0;
datetime Now=0;


//+------------------------------------------------------------------+
//| Custom indicator initialization function                         |
//+------------------------------------------------------------------+
int init() {

   IndicatorShortName("Magnetic field: "+HeavenlyBody);
   IndicatorBuffers(2);
   SetIndexBuffer(0,Result,DRAW_LINE);
   SetIndexBuffer(1,Extrapolation,DRAW_LINE);
   SetIndexShift(1,ExtrapolatedBars);   
   
   flag=SEFLG_TRUEPOS+SEFLG_XYZ+SEFLG_SPEED;
   return(0);
}

//+------------------------------------------------------------------+
//| Custom indicator deinitialization function                       |
//+------------------------------------------------------------------+
int deinit() {

   swe_close();
   return(0);
}


double scalar(double &a[], double &b[]){ return (a[0]*b[0]+a[1]*b[1]+a[2]*b[2]);}
double norm(double &a[]){return (MathSqrt(scalar(a,a)));}
double angle(double &a[], double &b[]){return ((180.0/3.14159)*MathArcsin(scalar(a,b)/(norm(a)*norm(b))));};
double vector(double &a[], double &b[], double &c[]){ c[0]=a[1]*b[2]-b[1]*a[2]; c[1]=a[2]*b[0]-a[0]*b[2]; c[2]=a[0]*b[1]-b[0]*a[1]; return(0);}


//+------------------------------------------------------------------+
//| Custom indicator iteration function                              |
//+------------------------------------------------------------------+
int start() {
                
   int    i, limit, counted_bars=IndicatorCounted(); 
   if (counted_bars>0) counted_bars--;
   limit=Bars-counted_bars;                
   
                
   double a[6], b[6], c[3];  


   if (Time[0]!=Now)
   {
   Now=Time[0];  
   //do this once per bar:
   
   for( i=limit-1; i>=0; i--) {   
      calculate(Time[i],HeavenlyBody,a,flag);
      c[0]=a[3]; c[1]=a[4]; c[2]=a[5]; 
      vector(c,a,b);  
      Result[i]=norm(b)/(norm(a)*norm(a)*norm(a));
   }
   

   //do the extrapolation once per bar:
   for (int j=0; j<ExtrapolatedBars; j++)
      {
      calculate(Time[0]+(ExtrapolatedBars-j)*PeriodSeconds(),HeavenlyBody,a,flag);
      c[0]=a[3]; c[1]=a[4]; c[2]=a[5]; 
      vector(c,a,b);  
      Extrapolation[j]=norm(b)/(norm(a)*norm(a)*norm(a));     
      }
   
   }
                  
   return(0);
}