import proj4 from 'proj4';


// Development tools

// Map application with coordinates projection
// https://dominoc925-pages.appspot.com/mapplets/cs_ecef.html

// The World coordinates converter
// https://twcc.fr/#

export const Projections = [
  {
    title: "RGF93 (Lambert-93)",
    url: "http://spatialreference.org/ref/epsg/rgf93-lambert-93/",
    code: "EPSG:2154",
    proj4: "+proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
  },
  {
    title: "RGF93 (CC50)",
    code: "EPSG:3950",
    url: "http://spatialreference.org/ref/epsg/3950/",
    proj4: "+proj=lcc +lat_1=49.25 +lat_2=50.75 +lat_0=50 +lon_0=3 +x_0=1700000 +y_0=9200000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
  },
  {
    title: "RGF93 (CC49)",
    code: "EPSG:3949",
    url: "http://spatialreference.org/ref/epsg/3949/",
    proj4: "+proj=lcc +lat_1=48.25 +lat_2=49.75 +lat_0=49 +lon_0=3 +x_0=1700000 +y_0=8200000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
  },
  {
    title: "RGF93 (CC48)",
    code: "EPSG:3948",
    url: "http://spatialreference.org/ref/epsg/3948/",
    proj4: "+proj=lcc +lat_1=47.25 +lat_2=48.75 +lat_0=48 +lon_0=3 +x_0=1700000 +y_0=7200000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
  },
  {
    title: "RGF93 (CC47)",
    code: "EPSG:3947",
    url: "http://spatialreference.org/ref/epsg/3947/",
    proj4: "+proj=lcc +lat_1=46.25 +lat_2=47.75 +lat_0=47 +lon_0=3 +x_0=1700000 +y_0=6200000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
  },
  {
    title: "RGF93 (CC46)",
    code: "EPSG:3946",
    url: "http://spatialreference.org/ref/epsg/3946/",
    proj4: "+proj=lcc +lat_1=45.25 +lat_2=46.75 +lat_0=46 +lon_0=3 +x_0=1700000 +y_0=5200000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
  },
  {
    title: "RGF93 (CC45)",
    code: "EPSG:3945",
    url: "http://spatialreference.org/ref/epsg/3945/",
    proj4: "+proj=lcc +lat_1=44.25 +lat_2=45.75 +lat_0=45 +lon_0=3 +x_0=1700000 +y_0=4200000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
  },
  {
    title: "RGF93 (CC44)",
    code: "EPSG:3944",
    url: "http://spatialreference.org/ref/epsg/3944/",
    proj4: "+proj=lcc +lat_1=43.25 +lat_2=44.75 +lat_0=44 +lon_0=3 +x_0=1700000 +y_0=3200000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
  },
  {
    title: "RGF93 (CC43)",
    code: "EPSG:3943",
    url: "http://spatialreference.org/ref/epsg/3943/",
    proj4: "+proj=lcc +lat_1=42.25 +lat_2=43.75 +lat_0=43 +lon_0=3 +x_0=1700000 +y_0=2200000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
  },
  {
    title: "RGF93 (CC42)",
    code: "EPSG:3942",
    url: "http://spatialreference.org/ref/epsg/3942/",
    proj4: "+proj=lcc +lat_1=41.25 +lat_2=42.75 +lat_0=42 +lon_0=3 +x_0=1700000 +y_0=1200000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
  },
  {
    title: "NTF (Lambert II)",
    code: "EPSG:27572",
    url: "http://spatialreference.org/ref/epsg/27572/",
    proj4: "+proj=lcc +lat_1=46.8 +lat_0=46.8 +lon_0=0 +k_0=0.99987742 +x_0=600000 +y_0=2200000 +a=6378249.2 +b=6356515 +towgs84=-168,-60,320,0,0,0,0 +pm=paris +units=m +no_defs"
  },
  {
    title: "NTF (Lambert III)",
    code: "EPSG:27573",
    url: "http://spatialreference.org/ref/epsg/27573/",
    proj4: "+proj=lcc +lat_1=44.10000000000001 +lat_0=44.10000000000001 +lon_0=0 +k_0=0.999877499 +x_0=600000 +y_0=3200000 +a=6378249.2 +b=6356515 +towgs84=-168,-60,320,0,0,0,0 +pm=paris +units=m +no_defs"
  },
  {
    title: "NTF",
    code: "EPSG:4275",
    url: "http://spatialreference.org/ref/epsg/4275/",
    proj4: "+proj=longlat +a=6378249.2 +b=6356515 +towgs84=-168,-60,320,0,0,0,0 +no_defs"
  },
  {
    title: "WGS84",
    code: "EPSG:4326",
    url: "https://epsg.io/4326",
    proj4: "+proj=longlat +datum=WGS84 +no_defs"
  }
];

export const Proj4Defs = Projections.map(p => [p.code, p.proj4]);

proj4.defs(Proj4Defs);

export function project(from, to, y, x) {
  console.log(`Projecting coordinates '${x},${y}' from '${from}' to '${to}'`);
  const coords = proj4(from, to, [x, y]);
  return { x: coords[0], y: coords[1] };
};

// Voir https://geodesie.ign.fr/contenu/fichiers/documentation/pedagogiques/TransformationsCoordonneesGeodesiques.pdf
export const WGS84Ellipsoid = createEllipsoid(6378137, 1/298.257223563);

function createEllipsoid(a, f) {
  const b = a * (1 - f);
  const asqr = a * a;
  const bsqr = b * b;
  const e = Math.sqrt((asqr-bsqr)/asqr);
  const eprime = Math.sqrt((asqr-bsqr)/bsqr);
  return { a, f, b, e, bsqr, asqr, eprime };
}

function degrees(angle) {
  return angle * (180 / Math.PI);
}

export function XYZtoLLA(x, y, z, ellipsoid) {

  const p = Math.sqrt(x*x + y*y);
  const theta = Math.atan((z*ellipsoid.a)/(p*ellipsoid.b));

  const sintheta = Math.sin(theta);
  const costheta = Math.cos(theta);

  const num = z + ellipsoid.eprime * ellipsoid.eprime * ellipsoid.b * sintheta * sintheta * sintheta;
  const denom = p - ellipsoid.e * ellipsoid.e * ellipsoid.a * costheta * costheta * costheta;

  const lat  = Math.atan(num/denom);
  let lon = Math.atan(y/x);
  const n = getN(lat, ellipsoid);
  const alt  = (p / Math.cos(lat)) - n;

  if (x < 0 && y < 0) {
    lon = lon - Math.PI;
  }

  if (x < 0 && y > 0) {
    lon = lon + Math.PI;
  }

  return { latitude: degrees(lat), longitude: degrees(lon), altitude: alt };

}

function getN(lat, ellipsoid) {
  const sinlat = Math.sin(lat);
  const denom = Math.sqrt(1-ellipsoid.e*ellipsoid.e*sinlat*sinlat);
  const n = ellipsoid.a / denom;
  return n;
}

export function LLAtoXYZ(latitude, longitude, altitude, ellipsoid) {

  const cosLat = Math.cos(latitude * Math.PI / 180.0);
  const sinLat = Math.sin(latitude * Math.PI / 180.0);
  const cosLon = Math.cos(longitude * Math.PI / 180.0);
  const sinLon = Math.sin(longitude * Math.PI / 180.0);
  
  const rad = ellipsoid.a;
  const f = ellipsoid.f;

  const C = 1.0 / Math.sqrt(cosLat * cosLat + (1 - f) * (1 - f) * sinLat * sinLat);
  const S = (1.0 - f) * (1.0 - f) * C;
  const h = altitude;

  const x = (rad * C + h) * cosLat * cosLon;
  const y = (rad * C + h) * cosLat * sinLon;
  const z = (rad * S + h) * sinLat;

  return { x, y, z };

}