import adjust_lon from '../common/adjust_lon'; import {HALF_PI, EPSLN} from '../constants/values'; import mlfn from '../common/mlfn'; import e0fn from '../common/e0fn'; import e1fn from '../common/e1fn'; import e2fn from '../common/e2fn'; import e3fn from '../common/e3fn'; import gN from '../common/gN'; import asinz from '../common/asinz'; import imlfn from '../common/imlfn'; export function init() { this.sin_p12 = Math.sin(this.lat0); this.cos_p12 = Math.cos(this.lat0); } export function forward(p) { var lon = p.x; var lat = p.y; var sinphi = Math.sin(p.y); var cosphi = Math.cos(p.y); var dlon = adjust_lon(lon - this.long0); var e0, e1, e2, e3, Mlp, Ml, tanphi, Nl1, Nl, psi, Az, G, H, GH, Hs, c, kp, cos_c, s, s2, s3, s4, s5; if (this.sphere) { if (Math.abs(this.sin_p12 - 1) <= EPSLN) { //North Pole case p.x = this.x0 + this.a * (HALF_PI - lat) * Math.sin(dlon); p.y = this.y0 - this.a * (HALF_PI - lat) * Math.cos(dlon); return p; } else if (Math.abs(this.sin_p12 + 1) <= EPSLN) { //South Pole case p.x = this.x0 + this.a * (HALF_PI + lat) * Math.sin(dlon); p.y = this.y0 + this.a * (HALF_PI + lat) * Math.cos(dlon); return p; } else { //default case cos_c = this.sin_p12 * sinphi + this.cos_p12 * cosphi * Math.cos(dlon); c = Math.acos(cos_c); kp = c ? c / Math.sin(c) : 1; p.x = this.x0 + this.a * kp * cosphi * Math.sin(dlon); p.y = this.y0 + this.a * kp * (this.cos_p12 * sinphi - this.sin_p12 * cosphi * Math.cos(dlon)); return p; } } else { e0 = e0fn(this.es); e1 = e1fn(this.es); e2 = e2fn(this.es); e3 = e3fn(this.es); if (Math.abs(this.sin_p12 - 1) <= EPSLN) { //North Pole case Mlp = this.a * mlfn(e0, e1, e2, e3, HALF_PI); Ml = this.a * mlfn(e0, e1, e2, e3, lat); p.x = this.x0 + (Mlp - Ml) * Math.sin(dlon); p.y = this.y0 - (Mlp - Ml) * Math.cos(dlon); return p; } else if (Math.abs(this.sin_p12 + 1) <= EPSLN) { //South Pole case Mlp = this.a * mlfn(e0, e1, e2, e3, HALF_PI); Ml = this.a * mlfn(e0, e1, e2, e3, lat); p.x = this.x0 + (Mlp + Ml) * Math.sin(dlon); p.y = this.y0 + (Mlp + Ml) * Math.cos(dlon); return p; } else { //Default case tanphi = sinphi / cosphi; Nl1 = gN(this.a, this.e, this.sin_p12); Nl = gN(this.a, this.e, sinphi); psi = Math.atan((1 - this.es) * tanphi + this.es * Nl1 * this.sin_p12 / (Nl * cosphi)); Az = Math.atan2(Math.sin(dlon), this.cos_p12 * Math.tan(psi) - this.sin_p12 * Math.cos(dlon)); if (Az === 0) { s = Math.asin(this.cos_p12 * Math.sin(psi) - this.sin_p12 * Math.cos(psi)); } else if (Math.abs(Math.abs(Az) - Math.PI) <= EPSLN) { s = -Math.asin(this.cos_p12 * Math.sin(psi) - this.sin_p12 * Math.cos(psi)); } else { s = Math.asin(Math.sin(dlon) * Math.cos(psi) / Math.sin(Az)); } G = this.e * this.sin_p12 / Math.sqrt(1 - this.es); H = this.e * this.cos_p12 * Math.cos(Az) / Math.sqrt(1 - this.es); GH = G * H; Hs = H * H; s2 = s * s; s3 = s2 * s; s4 = s3 * s; s5 = s4 * s; c = Nl1 * s * (1 - s2 * Hs * (1 - Hs) / 6 + s3 / 8 * GH * (1 - 2 * Hs) + s4 / 120 * (Hs * (4 - 7 * Hs) - 3 * G * G * (1 - 7 * Hs)) - s5 / 48 * GH); p.x = this.x0 + c * Math.sin(Az); p.y = this.y0 + c * Math.cos(Az); return p; } } } export function inverse(p) { p.x -= this.x0; p.y -= this.y0; var rh, z, sinz, cosz, lon, lat, con, e0, e1, e2, e3, Mlp, M, N1, psi, Az, cosAz, tmp, A, B, D, Ee, F, sinpsi; if (this.sphere) { rh = Math.sqrt(p.x * p.x + p.y * p.y); if (rh > (2 * HALF_PI * this.a)) { return; } z = rh / this.a; sinz = Math.sin(z); cosz = Math.cos(z); lon = this.long0; if (Math.abs(rh) <= EPSLN) { lat = this.lat0; } else { lat = asinz(cosz * this.sin_p12 + (p.y * sinz * this.cos_p12) / rh); con = Math.abs(this.lat0) - HALF_PI; if (Math.abs(con) <= EPSLN) { if (this.lat0 >= 0) { lon = adjust_lon(this.long0 + Math.atan2(p.x, - p.y)); } else { lon = adjust_lon(this.long0 - Math.atan2(-p.x, p.y)); } } else { /*con = cosz - this.sin_p12 * Math.sin(lat); if ((Math.abs(con) < EPSLN) && (Math.abs(p.x) < EPSLN)) { //no-op, just keep the lon value as is } else { var temp = Math.atan2((p.x * sinz * this.cos_p12), (con * rh)); lon = adjust_lon(this.long0 + Math.atan2((p.x * sinz * this.cos_p12), (con * rh))); }*/ lon = adjust_lon(this.long0 + Math.atan2(p.x * sinz, rh * this.cos_p12 * cosz - p.y * this.sin_p12 * sinz)); } } p.x = lon; p.y = lat; return p; } else { e0 = e0fn(this.es); e1 = e1fn(this.es); e2 = e2fn(this.es); e3 = e3fn(this.es); if (Math.abs(this.sin_p12 - 1) <= EPSLN) { //North pole case Mlp = this.a * mlfn(e0, e1, e2, e3, HALF_PI); rh = Math.sqrt(p.x * p.x + p.y * p.y); M = Mlp - rh; lat = imlfn(M / this.a, e0, e1, e2, e3); lon = adjust_lon(this.long0 + Math.atan2(p.x, - 1 * p.y)); p.x = lon; p.y = lat; return p; } else if (Math.abs(this.sin_p12 + 1) <= EPSLN) { //South pole case Mlp = this.a * mlfn(e0, e1, e2, e3, HALF_PI); rh = Math.sqrt(p.x * p.x + p.y * p.y); M = rh - Mlp; lat = imlfn(M / this.a, e0, e1, e2, e3); lon = adjust_lon(this.long0 + Math.atan2(p.x, p.y)); p.x = lon; p.y = lat; return p; } else { //default case rh = Math.sqrt(p.x * p.x + p.y * p.y); Az = Math.atan2(p.x, p.y); N1 = gN(this.a, this.e, this.sin_p12); cosAz = Math.cos(Az); tmp = this.e * this.cos_p12 * cosAz; A = -tmp * tmp / (1 - this.es); B = 3 * this.es * (1 - A) * this.sin_p12 * this.cos_p12 * cosAz / (1 - this.es); D = rh / N1; Ee = D - A * (1 + A) * Math.pow(D, 3) / 6 - B * (1 + 3 * A) * Math.pow(D, 4) / 24; F = 1 - A * Ee * Ee / 2 - D * Ee * Ee * Ee / 6; psi = Math.asin(this.sin_p12 * Math.cos(Ee) + this.cos_p12 * Math.sin(Ee) * cosAz); lon = adjust_lon(this.long0 + Math.asin(Math.sin(Az) * Math.sin(Ee) / Math.cos(psi))); sinpsi = Math.sin(psi); lat = Math.atan2((sinpsi - this.es * F * this.sin_p12) * Math.tan(psi), sinpsi * (1 - this.es)); p.x = lon; p.y = lat; return p; } } } export var names = ["Azimuthal_Equidistant", "aeqd"]; export default { init: init, forward: forward, inverse: inverse, names: names };