explore longitude array problem

[mdsumner]

f <- 360
while(TRUE) {

lon <- seq(-180, 180, length.out = f)
f <- f * 5
if (length(unique(diff(lon)) > 1)) stop() 
}

range(diff(lon))

src <-
"#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
NumericVector rcpp_double(IntegerVector n){
  NumericVector v(n[0]);
  double start = -180.0; 
  double end = 180.0;
  double dv = (end - start)/n[0]; 
  for(int i=0; i<v.length(); ++i){
    v[i] = start + i * dv; 
  }
  return(v);
}

// [[Rcpp::export]]
NumericVector rcpp_float(IntegerVector n){
  NumericVector v(n[0]);
  double start = -180.0; 
  double end = 180.0;
  float dv = (float)(end - start)/(double)n[0]; 
  for(int i=0; i<v.length(); ++i){
   float val = start + i * dv; 
    v[i] = val; 
  }
  return(v);
}

"
library(Rcpp)
sourceCpp(code = src)
eps <- 4 * (1/1e5) #delta / 500

n_vals <- c(360 * 10^(seq(5, 0, by = -0.5)))
df <- numeric(length(n_vals))
plt <- F
if (plt) par(mfrow = c(length(n_vals), 1L), mar = rep(0, 4))

for (i in seq_along(n_vals)) {
  n <- n_vals[i]
  delta <- 360/n
  lon <- rcpp_double(n)
  lonf <- rcpp_float(n)

  if (plt) {
plot(NA, type = "n", xlab = "", ylab = "", xlim = c(-180, 180), ylim = c(-1, 1) * eps +  delta, main = n, sub = max(range(diff(lonf))))
points(lon[-1], diff(lon), pch = 19, cex = 0.8)
points(lonf[-1], diff(lonf), pch = 19, cex = 0.2, col = "red")
  }
  df[i ] <- max(abs(lon - lonf))
print(sprintf("n: %i, diff: %e, dist: %f", 
              as.integer(n), 
              df[i], geodist::geodist(cbind(lon = c(0, df[i]), lat = 0), measure = "geodesic", sequential = T)))
}