/* vim: set expandtab shiftwidth=2 softtabstop=2 tw=70: */

#include <Rcpp.h>
using namespace Rcpp;

// Cross-reference work:
// 1. update registerDynamicSymbol.c with an item for this
// 2. main code should use the autogenerated wrapper in ../R/RcppExports.R

//#define DEBUG 1

// [[Rcpp::export]]
List do_curl1(NumericMatrix u, NumericMatrix v, NumericVector x, NumericVector y, NumericVector geographical)
{
  double R = 6371.0e3; // FIXME: could be an argument but there seems little need
  int ncol = u.ncol(), nrow = u.nrow();
  if (nrow != v.nrow())
    ::Rf_error("matrices u and v must have nrow values but they are %d and %d", nrow, v.nrow());
  if (ncol != v.ncol())
    ::Rf_error("matrices u and v must have ncol values but they are %d and %d", ncol, v.ncol());
  int nx = x.size(), ny = y.size();
  if (nx != nrow)
    ::Rf_error("length(x)=%d does not match ncol(u)=%d\n", nx, ncol);
  if (ny != ncol)
    ::Rf_error("length(y)=%d does not match nrow(u)=%d\n", ny, nrow);
  int isGeographical = 0.0 != geographical[0];
  NumericMatrix curl(nrow, ncol);
  // SCHEME: x or longitude is indexed by 0 <= i <= (nrow-1), while
  //         y or latitude  is indexed by 0 <= j <= (ncol-1).
  // set to NA so we can see if we are failing to fill grid correctly
  for (int i = 0; i < nrow; i++)
    for (int j = 0; j < ncol; j++)
      curl(i, j) = NA_REAL;
  double xfac=1.0, yfac = 1.0;
  if (isGeographical)
    yfac = R * M_PI / 180.0;
  for (int j = 1; j < ncol-1; j++) {
    if (isGeographical)
      xfac = yfac * cos(y[j]*M_PI/180.0);
    for (int i = 1; i < nrow-1; i++) {
      // Calculate first difference with a 5-point stencil, e.g. infer d/dx by subtracting
      // the value at i-1 from the value at i+1 etc.
      double du = u(i,j+1) - u(i,j-1);
      double dv = v(i+1,j) - v(i-1,j);
      double dx = xfac * (x[i+1] - x[i-1]);
      double dy = yfac * (y[j+1] - y[j-1]);
      curl(i, j) = dv/dx - du/dy;
    }
  }
  // bottom and top: copy neighbours above and below
  for (int i = 1; i < nrow-1; i++) {
      curl(i,      0) = curl(i,      1);
      curl(i, ncol-1) = curl(i, ncol-2);
  }
  // left and right: copy neighbors to right and left
  for (int j = 1; j < ncol-1; j++) {
      curl(     0, j) = curl(     1, j);
      curl(nrow-1, j) = curl(nrow-2, j);
  }
  // corners: use diagonal neighbour
  curl(     0,      0) = curl(     1,      1);
  curl(     0, ncol-1) = curl(     1, ncol-2);
  curl(nrow-1, 0     ) = curl(nrow-2,      1);
  curl(nrow-1, ncol-1) = curl(nrow-2, ncol-2);
  return(List::create(Named("x")=x, Named("y")=y, Named("curl")=curl));
}

// [[Rcpp::export]]
List do_curl2(NumericMatrix u, NumericMatrix v, NumericVector x, NumericVector y, NumericVector geographical)
{
  double R = 6371.0e3; // FIXME: could be an argument but there seems little need
  int nrow = u.nrow(), ncol = u.ncol();
  if (nrow != v.nrow())
    ::Rf_error("matrices u and v must have equal nrow values but they are %d and %d", nrow, v.nrow());
  if (ncol != v.ncol())
    ::Rf_error("matrices u and v must have equal ncol values but they are %d and %d", ncol, v.ncol());
  int nx = x.size(), ny = y.size();
  if (nx != nrow)
    ::Rf_error("length(x)=%d does not match ncol(u)=%d", nx, ncol);
  if (ny != ncol)
    ::Rf_error("length(y)=%d does not match nrow(u)=%d", ny, nrow);
  int isGeographical = 0.0 != geographical[0];
  NumericMatrix curl(nrow-1, ncol-1);
  // SCHEME: x or longitude is indexed by 0 <= i <= (nrow-1), while
  //         y or latitude  is indexed by 0 <= j <= (ncol-1).
  // set to NA so we can see if we are failing to fill grid correctly
  for (int i = 0; i < nrow-1; i++)
    for (int j = 0; j < ncol-1; j++)
      curl(i,j) = NA_REAL;
  double xfac=1.0, yfac = 1.0;
  if (isGeographical)
    yfac = R * M_PI / 180.0;
#ifdef DEBUG
  Rprintf("nrow= %d\n", nrow);
  Rprintf("ncol= %d\n", ncol);
  Rprintf("2.5*yfac=%.5e\n", 2.5*yfac);
#endif
  for (int j = 0; j < ncol-1; j++) {
    // For xfac, use cosine factor as the average of values at the north
    // and south sides of the grid box.
    if (isGeographical)
      xfac = yfac * 0.5*(cos(y[j]*M_PI/180.0)+cos(y[j+1]*M_PI/180.0));
#ifdef DEBUG
    Rprintf("at j=%d, have yp[%d]=%.5f and 2.5*xfac=%.5e\n", j  , j  , y[j  ], 2.5*yfac * cos(y[j  ]*M_PI/180.0));
    Rprintf("at j=%d, have yp[%d]=%.5f and 2.5*xfac=%.5e\n", j+1, j+1, y[j+1], 2.5*yfac * cos(y[j+1]*M_PI/180.0));
#endif
    for (int i = 0; i < nrow-1; i++) {
      // Calculate derivatives centred within each grid box. Thus, for
      // du, we must avg across x (i.e i) at two grid points, and then
      // difference across y (i.e. j).
      double du = 0.5*(u(i  , j+1) + u(i+1, j+1)) - 0.5*(u(i, j) + u(i+1, j  ));
      double dv = 0.5*(v(i+1, j  ) + v(i+1, j+1)) - 0.5*(v(i, j) + v(i  , j+1));
      double dx = xfac * (x[i+1] - x[i]);
      double dy = yfac * (y[j+1] - y[j]);
      curl(i, j) = dv/dx - du/dy;
#ifdef DEBUG
      if (i == 0 && j == 0) {
        Rprintf("du = 0.5*(%g + %g) - 0.5*(%g + %g)\n", u(i,j+1),u(i+1,j+1),u(i,j),u(i+1,j));
        Rprintf("dv = 0.5*(%g + %g) - 0.5*(%g + %g)\n", v(i+1,j),v(i+1,j+1),v(i,j),v(i,j+1));
        Rprintf("x[i=%d,%d]=(%.1f,%.1f), y[j=%d,%d]=(%.1f,%.1f)\n", i, i+1, x[i], x[i+1], j,j+1,y[j], y[j+1]);
        Rprintf("du=%.4f dv=%.4f dx=%g dy=%g dv/dx=%.3e du/dy=%.3e curl=%.3e\n",
            du,dv,dx,dy,dv/dx,du/dy,curl(i,j));
      }
#endif
    }
  }
  // Construct xnew and ynew, which we will return.
  NumericVector xnew(nrow-1), ynew(ncol-1);
  for (int i = 0; i < nrow - 1; i++)
    xnew[i] = 0.5 * (x[i] + x[i+1]);
  for (int j = 0; j < ncol - 1; j++)
    ynew[j] = 0.5 * (y[j] + y[j+1]);
  return(List::create(Named("x")=xnew, Named("y")=ynew, Named("curl")=curl));
}

