#include <iostream>

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

extern "C" {
  
__declspec(dllexport) 
void vc_avar(
    double y[],
    int n,
    double tau0, 
    double tau[],
    double ADEV[],
    double MDEV[],
    double TDEV[],
    int *N){
      
// function vc_avar : computes the ALLAN VARIANCE
// INPUTS : 
// y = pointer to data, already allocated, length = n
// n = length of data 
// tau0 = sampling periode in seconds
// OUTPUTS : 
// tau = pointer to tau , already allocated, length >= N
// ADEV = pointer to overlapped Allan STD DEV, already allocated , length >= N
// MDEV = pointer to modified Allan STD DEV, already allocated , length >=N
// TDEV = pointer to Time Allan STD DEV, already allocated , length >= N
// N = length of results tau, ADEV, MDEV and TDEV. 
//
// Copyright Alaa MAKDISSI 2003
// free for personal use only.
   int k;
   double  *D;
   int m;
   double s,v,u;
   *N=floor(1+ log((n-1.0)/3.0)/log(2.0) ); // nb of points for the results
   //*N = 2;
   m=1; // m = 2^j ; j=0 ==> tau=tau0
   D= new double [n];
   for (int j=0; j<*N;j++)
   {
   tau[j]=m*tau0;
   s=0; v=0; //s =Y(0) and v =Y(m)
   for (k=0;k<m;k++)
   {
   s=s+y[k]; 
   v=v+y[k+m]; 
   }
   D[0]=(v-s)/m;
//D(k)
   for (k=1;k<n+1-2*m;k++) D[k]=D[k-1]+(y[k-1]-2*y[k+m-1]+y[k+2*m-1])/m;
   //for (k=0;k<n+1-2*m;k++) std::cout << D[k] << std::endl;
//ADEV
   s=0;
   for (k=0;k<n+1-2*m;k++)
   {
   v=D[k]*D[k];
   s=s+v;
   }
   ADEV[j]=sqrt(s/(n+1-2*m)/2);
//MDEV
   v=0 ;
   for(int j=0;j<n+2-(3*m);j++)
   {
   s=0;
   for( k=0;k<m;k++)
   {
   s=s+D[j+k];
   }
   v=v+s*s;
   }
   v=v/(n+2-3*m);
   MDEV[j]=sqrt(v/2)/m;
//TDEV
   TDEV[j]=MDEV[j]*m*tau0/sqrt(3.0);
   m=m*2; // m = 2^j
   }
   delete [] D;
} //end vc_avar
}

int main(){
  double test[] = {0.5,-0.5,0.25,-0.25,0.5,-0.5};
  double tau[16],ADEV[16],MDEV[16],TDEV[16];
  int rlen;
  vc_avar(test,6,0.01,tau,ADEV,MDEV,TDEV,&rlen);
  
  for(int i = 0; i < rlen; i++){
    std::cout << tau[i]
              << ", ADEV:" << ADEV[i]
              << ", MDEV:" << MDEV[i]
              << ", TDEV:" << TDEV[i]
              << std::endl; 
  }
  
  return 0; 
}