/*
 Just a test to see what this mex stuff can do
 */

#include <iostream>
#include <cstring>
#include <sstream>
#include <string>

#include <math.h>
#include "mex.h"

using namespace std;

 void mexFunction( int nlhs, mxArray *plhs[], 
		  int nrhs, const mxArray*prhs[] )
{
 
 double*X = mxGetPr(prhs[0]);
 double*T = mxGetPr(prhs[1]);
 
 double test;
 
 int index_X,index_T;
 double first_comp,second_comp,third_comp;
     
 int m_x,n_x;
 m_x = (int) mxGetM(prhs[0]); 
 n_x = (int) mxGetN(prhs[0]);
 
 double std_gaussian = 5;
 
//  plhs[0] = mxCreateDouble(1,1,mxREAL);
 plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL);
 double * data = (double *) mxGetData(plhs[0]);

 int m_t,n_t;
 m_t = (int) mxGetM(prhs[1]); 
 n_t = (int) mxGetN(prhs[1]);
 
//  data[0] = 9;
// //  
//  mexPrintf("%d",m_x);
//  mexPrintf("%d",n_x);
//  mexPrintf("%d",m_t);
//  mexPrintf("%d",n_t);
 
//  mexPrintf("%lf",126.46);
 
 data[0] = 0;
 

 
 for (int i = 0;i < m_x; i++)
{
    for (int j = 0;j < m_t; j++)
    {
        
        index_X = m_x*(0) + i;
        index_T = m_t*(0) + j;
        
        first_comp = *(X + index_X) - *(T + index_T);
        
        index_X = m_x*(1) + i;
        index_T = m_t*(1) + j;
         
        second_comp = *(X + index_X) - *(T + index_T);
    
        index_X = m_x*(2) + i;
        index_T = m_t*(2) + j;
        
        third_comp = *(X + index_X) - *(T + index_T);
//         
//         mexPrintf("%lf",first_comp);
//         mexPrintf("\n");
//         mexPrintf("%lf",second_comp);
//                 mexPrintf("\n");
// 
//         mexPrintf("%lf",third_comp);
//                 mexPrintf("\n");

        test = (pow(first_comp,2) + pow(second_comp,2) + pow(third_comp,2));
//         mexPrintf("%lf",test);
//                 mexPrintf("\n");

        if(test <= 36)
        {
            data[0] = data[0] + 1;
            break;
        }
        
        
//         data[0] = data[0] + exp(-(pow(first_comp,2) + pow(second_comp,2) + pow(third_comp,2))/pow(std_gaussian,2));
    }
 }
 
//  *mxGetPr(plhs[0]) = 8;
 
//      double*answer = mxGetPr(plhs[0]);
//       mexPrintf("%d",m);
//       mexPrintf("%d",n);

//      answer = 8.0;
//      double*answer = mxGetPr(plhs[0]);
//      double intermediate = (-4^2)*-1;
//      double answer_zero = exp(intermediate);
//      answer = *answer_zero;
 }