//
//  MarkVariogram.cpp
//  
//
//  Created by Arne Pommerening on 1/25/13.
//  Copyright 2013 Philodendron International. All rights reserved.
//

#include <Rcpp.h>
#include <cmath>
#include <ctime>
using namespace Rcpp;
using namespace std;


/**
 * Estimates the mark variogram and the corresponding marked L function.
 */
// [[Rcpp::export]]
DataFrame calcMarkVariogram(DataFrame TreeList, double xmax, double ymax, int m, double dd, double h, double mvar, int kernel) {
	vector<double> x = as< vector<double> > (TreeList["x"]);
    vector<double> y = as< vector<double> > (TreeList["y"]);
    vector<double> dbh = as< vector<double> > (TreeList["dbh"]);
    const double pi = 4.0 * atan(1.0);
    int n = dbh.size();
    double la = n / (xmax * ymax);
    double ql = la * la;
	NumericVector g(m), xv(m), r(m), xl(m), xk(m);
    for(int i = 0; i < m; i++) {
        g[i] = 0;
        xk[i] = 0;
        xv[i] = 0;
        xl[i] = 0;
		r[i] = 0;
    }
    double rz = m * dd + h;
    for(int i = 1; i < n; i++) {
        double xx = x[i];
        double yy = y[i];
        double zz = dbh[i];
        for(int j = 0; j < i; j++) {
            double cx = fabs(xx - x[j]);
            if(cx < rz) {
                double cy = fabs(yy - y[j]);
                if(cy < rz) {
                    double tz = sqrt(cx * cx + cy * cy);
                    if(tz < rz) {
                        double cc = 1 / ((xmax - cx) * (ymax - cy));
                        for(int l = 0; l < m; l++) {
                            double xr = (l + 1) * dd;
                            if(tz <= xr) {
								xk[l] += 2 * cc; 
								xl[l] += 2 * cc * pow(zz - dbh[j], 2) * 0.5; 
							}
							double t = fabs(xr - tz);
                            if(t <= h) {
                                double rh = 0;
                                if(kernel == 0)
                                  rh = (1 - pow(t / h, 2)) * 0.75 / h;
                                else
                                  rh = 0.5 / h;
                                double cr = pi * xr;
                                rh *= cc / cr;
                                g[l] += rh;
                                xv[l] += rh * pow(zz - dbh[j], 2) * 0.5;
                            }
                        }
                    }
                }
            }    
        }
    }
    for(int l = 0; l < m; l++) {
		xk[l] = sqrt(xk[l] / (pi * ql));
		r[l] = (l + 1) * dd; 
        xl[l] = sqrt(xl[l] / (pi * ql * mvar));
		if(g[l] < 1E-20)
            xv[l] = -1;
        else
            xv[l] /= g[l];
    }
    return DataFrame::create(Named("r") = r, Named("L") = xk, Named("mv") = xv, Named("Lm") = xl);
}

/**
 * Calculates the minimum distance between points in an observation window.
 */
// [[Rcpp::export]]
double calcMinDistance(NumericVector x, NumericVector y) {
	double dm = 1E50;
	double xx = 0;
    int n = x.size();
    for(int i = 1; i < n; i++) {
        for(int j = 0; j < i; j++) {
			xx = pow(x[i] - x[j], 2) + pow(y[i] - y[j], 2);
			dm = min(xx, dm);
		}
	}
	return sqrt(dm);
}