diff --git a/html/geomag2020.js b/html/geomag2020.js new file mode 100644 index 0000000..8093d65 --- /dev/null +++ b/html/geomag2020.js @@ -0,0 +1,361 @@ +/* 2012-03-26 + Copyright 2012 Christopher Weiss (cmweiss@gmail.com) + + Suggestions for improvements are appreciated. + + Adapted from the geomagc software and World Magnetic Model of the NOAA + Satellite and Information Service, National Geophysical Data Center + http://www.ngdc.noaa.gov/geomag/WMM/DoDWMM.shtml + + geoMagFactory() requires a world magnetic model (WMM) object. The helper + function cof2Obj(), available in cof2Obj.js, takes the text of WMM.COF and + returns an object suitable for geoMagFactory(). A syncronous XMLHttpRequest + to fetch the WMM.COF is recommended in a web environment. The helper + function syncXHR(), available in syncXHR.js, takes the url of the WMM.COF + file and returns the WMM.COF file as text. + + Usage: + geoMagFactory(wmm) returns a function which can compute the Earth's + magnetic field. + The returned function requires two arguments, latitude and longitude (in + decimal degrees), and, optionally, altitude in feet (default is 0), and + a date object (default is the current system time). + + let cof = syncXHR('http://host/path/WMM.COF'), + wmm = cof2Obj(cof), + geoMag = geoMagFactory(wmm), + latitude = 40.0, // decimal degrees (north is positive) + longitude = -80.0, // decimal degrees (east is positive) + altitude = 0, // feet (optional, default is 0) + time = new Date(2012, 4, 20), // (optional, default is the current + // system time) + myGeoMag = geoMag(latitude, longitude, altitude, time), + magneticletiation = myGeoMag.dec, // Geomagnetic declination + // (letiation) in decimal degrees + // -- east is positive + magneticDip = myGeoMag.dip, // Geomagnetic dip in decimal degrees + // (down is positive) + magneticFieldIntensity = myGeoMag.ti, // Total intensity of the + // geomagnetic field in + // nanoteslas + magneticBH = myGeoMag.bh, // Horizontal intensity of the geomagnetic + // field in nT + magneticBX = myGeoMag.bx, // North component of the geomagnetic field + // in nT + magneticBY = myGeoMag.by, // East component of the geomagnetic field + // in nT + magneticBZ = myGeoMag.bz, // Vertical component of the geomagnetic + // field (down is positive) + lat = myGeoMag.lat, // input latitude + lon = myGeoMag.lon; // input longitude +*/ + +/*jslint plusplus: true */ +function geoMagFactory(wmm) { + 'use strict'; + function rad2deg(rad) { + return rad * (180 / Math.PI); + } + function deg2rad(deg) { + return deg * (Math.PI / 180); + } + + let i, model, epoch = wmm.epoch, + z = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + maxord = 12, + tc = [z.slice(), z.slice(), z.slice(), z.slice(), z.slice(), z.slice(), + z.slice(), z.slice(), z.slice(), z.slice(), z.slice(), z.slice(), + z.slice()], + sp = z.slice(), + cp = z.slice(), + pp = z.slice(), + p = [z.slice(), z.slice(), z.slice(), z.slice(), z.slice(), z.slice(), + z.slice(), z.slice(), z.slice(), z.slice(), z.slice(), z.slice(), + z.slice()], + dp = [z.slice(), z.slice(), z.slice(), z.slice(), z.slice(), z.slice(), + z.slice(), z.slice(), z.slice(), z.slice(), z.slice(), z.slice(), + z.slice()], + a = 6378.137, + b = 6356.7523142, + re = 6371.2, + a2 = a * a, + b2 = b * b, + c2 = a2 - b2, + a4 = a2 * a2, + b4 = b2 * b2, + c4 = a4 - b4, + c = [z.slice(), z.slice(), z.slice(), z.slice(), z.slice(), z.slice(), + z.slice(), z.slice(), z.slice(), z.slice(), z.slice(), z.slice(), + z.slice()], + cd = [z.slice(), z.slice(), z.slice(), z.slice(), z.slice(), z.slice(), + z.slice(), z.slice(), z.slice(), z.slice(), z.slice(), z.slice(), + z.slice()], + n, m, + snorm = [z.slice(), z.slice(), z.slice(), z.slice(), z.slice(), + z.slice(), z.slice(), z.slice(), z.slice(), z.slice(), z.slice(), + z.slice(), z.slice()], + j, + k = [z.slice(), z.slice(), z.slice(), z.slice(), z.slice(), z.slice(), + z.slice(), z.slice(), z.slice(), z.slice(), z.slice(), z.slice(), + z.slice()], + flnmj, + fn = [0, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13], + fm = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12], + D2; + + tc[0][0] = 0; + sp[0] = 0.0; + cp[0] = 1.0; + pp[0] = 1.0; + p[0][0] = 1; + + model = wmm.wmm; + for (i in model) { + if (model.hasOwnProperty(i)) { + if (model[i].m <= model[i].n) { + c[model[i].m][model[i].n] = model[i].gnm; + cd[model[i].m][model[i].n] = model[i].dgnm; + if (model[i].m !== 0) { + c[model[i].n][model[i].m - 1] = model[i].hnm; + cd[model[i].n][model[i].m - 1] = model[i].dhnm; + } + } + } + } + wmm = null; + model = null; + + /* CONVERT SCHMIDT NORMALIZED GAUSS COEFFICIENTS TO UNNORMALIZED */ + snorm[0][0] = 1; + + for (n = 1; n <= maxord; n++) { + snorm[0][n] = snorm[0][n - 1] * (2 * n - 1) / n; + j = 2; + + for (m = 0, D2 = (n - m + 1); D2 > 0; D2--, m++) { + k[m][n] = (((n - 1) * (n - 1)) - (m * m)) / + ((2 * n - 1) * (2 * n - 3)); + if (m > 0) { + flnmj = ((n - m + 1) * j) / (n + m); + snorm[m][n] = snorm[m - 1][n] * Math.sqrt(flnmj); + j = 1; + c[n][m - 1] = snorm[m][n] * c[n][m - 1]; + cd[n][m - 1] = snorm[m][n] * cd[n][m - 1]; + } + c[m][n] = snorm[m][n] * c[m][n]; + cd[m][n] = snorm[m][n] * cd[m][n]; + } + } + k[1][1] = 0.0; + + return function (glat, glon, h, date) { + function decimalDate(date) { + date = date || new Date(); + let year = date.getUTCFullYear(), + daysInYear = 365 + + (((year % 400 === 0) || (year % 4 === 0 && (year % 100 > 0))) ? 1 : 0), + msInYear = daysInYear * 24 * 60 * 60 * 1000; + + return date.getUTCFullYear() + (date.valueOf() - Date.UTC(year, 0)) / msInYear; + } + + let alt = (h / 3280.8399) || 0, // convert h (in feet) to kilometers or set default of 0 + time = decimalDate(date), + dt = time - epoch, + rlat = deg2rad(glat), + rlon = deg2rad(glon), + srlon = Math.sin(rlon), + srlat = Math.sin(rlat), + crlon = Math.cos(rlon), + crlat = Math.cos(rlat), + srlat2 = srlat * srlat, + crlat2 = crlat * crlat, + q, + q1, + q2, + ct, + st, + r2, + r, + d, + ca, + sa, + aor, + ar, + br = 0.0, + bt = 0.0, + bp = 0.0, + bpp = 0.0, + par, + temp1, + temp2, + parp, + D4, + bx, + by, + bz, + bh, + ti, + dec, + dip, + gv; + sp[1] = srlon; + cp[1] = crlon; + + /* CONVERT FROM GEODETIC COORDS. TO SPHERICAL COORDS. */ + q = Math.sqrt(a2 - c2 * srlat2); + q1 = alt * q; + q2 = ((q1 + a2) / (q1 + b2)) * ((q1 + a2) / (q1 + b2)); + ct = srlat / Math.sqrt(q2 * crlat2 + srlat2); + st = Math.sqrt(1.0 - (ct * ct)); + r2 = (alt * alt) + 2.0 * q1 + (a4 - c4 * srlat2) / (q * q); + r = Math.sqrt(r2); + d = Math.sqrt(a2 * crlat2 + b2 * srlat2); + ca = (alt + d) / r; + sa = c2 * crlat * srlat / (r * d); + + for (m = 2; m <= maxord; m++) { + sp[m] = sp[1] * cp[m - 1] + cp[1] * sp[m - 1]; + cp[m] = cp[1] * cp[m - 1] - sp[1] * sp[m - 1]; + } + + aor = re / r; + ar = aor * aor; + + for (n = 1; n <= maxord; n++) { + ar = ar * aor; + for (m = 0, D4 = (n + m + 1); D4 > 0; D4--, m++) { + + /* + COMPUTE UNNORMALIZED ASSOCIATED LEGENDRE POLYNOMIALS + AND DERIVATIVES VIA RECURSION RELATIONS + */ + if (n === m) { + p[m][n] = st * p[m - 1][n - 1]; + dp[m][n] = st * dp[m - 1][n - 1] + ct * + p[m - 1][n - 1]; + } else if (n === 1 && m === 0) { + p[m][n] = ct * p[m][n - 1]; + dp[m][n] = ct * dp[m][n - 1] - st * p[m][n - 1]; + } else if (n > 1 && n !== m) { + if (m > n - 2) { p[m][n - 2] = 0; } + if (m > n - 2) { dp[m][n - 2] = 0.0; } + p[m][n] = ct * p[m][n - 1] - k[m][n] * p[m][n - 2]; + dp[m][n] = ct * dp[m][n - 1] - st * p[m][n - 1] - + k[m][n] * dp[m][n - 2]; + } + + /* + TIME ADJUST THE GAUSS COEFFICIENTS + */ + + tc[m][n] = c[m][n] + dt * cd[m][n]; + if (m !== 0) { + tc[n][m - 1] = c[n][m - 1] + dt * cd[n][m - 1]; + } + + /* + ACCUMULATE TERMS OF THE SPHERICAL HARMONIC EXPANSIONS + */ + par = ar * p[m][n]; + if (m === 0) { + temp1 = tc[m][n] * cp[m]; + temp2 = tc[m][n] * sp[m]; + } else { + temp1 = tc[m][n] * cp[m] + tc[n][m - 1] * sp[m]; + temp2 = tc[m][n] * sp[m] - tc[n][m - 1] * cp[m]; + } + bt = bt - ar * temp1 * dp[m][n]; + bp += (fm[m] * temp2 * par); + br += (fn[n] * temp1 * par); + /* + SPECIAL CASE: NORTH/SOUTH GEOGRAPHIC POLES + */ + if (st === 0.0 && m === 1) { + if (n === 1) { + pp[n] = pp[n - 1]; + } else { + pp[n] = ct * pp[n - 1] - k[m][n] * pp[n - 2]; + } + parp = ar * pp[n]; + bpp += (fm[m] * temp2 * parp); + } + } + } + + bp = (st === 0.0 ? bpp : bp / st); + /* + ROTATE MAGNETIC VECTOR COMPONENTS FROM SPHERICAL TO + GEODETIC COORDINATES + */ + bx = -bt * ca - br * sa; + by = bp; + bz = bt * sa - br * ca; + + /* + COMPUTE DECLINATION (DEC), INCLINATION (DIP) AND + TOTAL INTENSITY (TI) + */ + bh = Math.sqrt((bx * bx) + (by * by)); + ti = Math.sqrt((bh * bh) + (bz * bz)); + dec = rad2deg(Math.atan2(by, bx)); + dip = rad2deg(Math.atan2(bz, bh)); + + /* + COMPUTE MAGNETIC GRID letIATION IF THE CURRENT + GEODETIC POSITION IS IN THE ARCTIC OR ANTARCTIC + (I.E. GLAT > +55 DEGREES OR GLAT < -55 DEGREES) + OTHERWISE, SET MAGNETIC GRID letIATION TO -999.0 + */ + + if (Math.abs(glat) >= 55.0) { + if (glat > 0.0 && glon >= 0.0) { + gv = dec - glon; + } else if (glat > 0.0 && glon < 0.0) { + gv = dec + Math.abs(glon); + } else if (glat < 0.0 && glon >= 0.0) { + gv = dec + glon; + } else if (glat < 0.0 && glon < 0.0) { + gv = dec - Math.abs(glon); + } + if (gv > 180.0) { + gv -= 360.0; + } else if (gv < -180.0) { gv += 360.0; } + } + + return {dec: dec, dip: dip, ti: ti, bh: bh, bx: bx, by: by, bz: bz, lat: glat, lon: glon, gv: gv, epoch: epoch}; + }; +} +/* + cof2Obj.js + Converts the WMM.COF text to a JSON object usable by geoMagFactory(). +*/ + +function cof2Obj() { +let cof=" 2020.0 WMM-2020 12/10/2019\n 1 0 -29404.5 0.0 6.7 0.0\n 1 1 -1450.7 4652.9 7.7 -25.1\n 2 0 -2500.0 0.0 -11.5 0.0\n 2 1 2982.0 -2991.6 -7.1 -30.2\n 2 2 1676.8 -734.8 -2.2 -23.9\n 3 0 1363.9 0.0 2.8 0.0\n 3 1 -2381.0 -82.2 -6.2 5.7\n 3 2 1236.2 241.8 3.4 -1.0\n 3 3 525.7 -542.9 -12.2 1.1\n 4 0 903.1 0.0 -1.1 0.0\n 4 1 809.4 282.0 -1.6 0.2\n 4 2 86.2 -158.4 -6.0 6.9\n 4 3 -309.4 199.8 5.4 3.7\n 4 4 47.9 -350.1 -5.5 -5.6\n 5 0 -234.4 0.0 -0.3 0.0\n 5 1 363.1 47.7 0.6 0.1\n 5 2 187.8 208.4 -0.7 2.5\n 5 3 -140.7 -121.3 0.1 -0.9\n 5 4 -151.2 32.2 1.2 3.0\n 5 5 13.7 99.1 1.0 0.5\n 6 0 65.9 0.0 -0.6 0.0\n 6 1 65.6 -19.1 -0.4 0.1\n 6 2 73.0 25.0 0.5 -1.8\n 6 3 -121.5 52.7 1.4 -1.4\n 6 4 -36.2 -64.4 -1.4 0.9\n 6 5 13.5 9.0 -0.0 0.1\n 6 6 -64.7 68.1 0.8 1.0\n 7 0 80.6 0.0 -0.1 0.0\n 7 1 -76.8 -51.4 -0.3 0.5\n 7 2 -8.3 -16.8 -0.1 0.6\n 7 3 56.5 2.3 0.7 -0.7\n 7 4 15.8 23.5 0.2 -0.2\n 7 5 6.4 -2.2 -0.5 -1.2\n 7 6 -7.2 -27.2 -0.8 0.2\n 7 7 9.8 -1.9 1.0 0.3\n 8 0 23.6 0.0 -0.1 0.0\n 8 1 9.8 8.4 0.1 -0.3\n 8 2 -17.5 -15.3 -0.1 0.7\n 8 3 -0.4 12.8 0.5 -0.2\n 8 4 -21.1 -11.8 -0.1 0.5\n 8 5 15.3 14.9 0.4 -0.3\n 8 6 13.7 3.6 0.5 -0.5\n 8 7 -16.5 -6.9 0.0 0.4\n 8 8 -0.3 2.8 0.4 0.1\n 9 0 5.0 0.0 -0.1 0.0\n 9 1 8.2 -23.3 -0.2 -0.3\n 9 2 2.9 11.1 -0.0 0.2\n 9 3 -1.4 9.8 0.4 -0.4\n 9 4 -1.1 -5.1 -0.3 0.4\n 9 5 -13.3 -6.2 -0.0 0.1\n 9 6 1.1 7.8 0.3 -0.0\n 9 7 8.9 0.4 -0.0 -0.2\n 9 8 -9.3 -1.5 -0.0 0.5\n 9 9 -11.9 9.7 -0.4 0.2\n 10 0 -1.9 0.0 0.0 0.0\n 10 1 -6.2 3.4 -0.0 -0.0\n 10 2 -0.1 -0.2 -0.0 0.1\n 10 3 1.7 3.5 0.2 -0.3\n 10 4 -0.9 4.8 -0.1 0.1\n 10 5 0.6 -8.6 -0.2 -0.2\n 10 6 -0.9 -0.1 -0.0 0.1\n 10 7 1.9 -4.2 -0.1 -0.0\n 10 8 1.4 -3.4 -0.2 -0.1\n 10 9 -2.4 -0.1 -0.1 0.2\n 10 10 -3.9 -8.8 -0.0 -0.0\n 11 0 3.0 0.0 -0.0 0.0\n 11 1 -1.4 -0.0 -0.1 -0.0\n 11 2 -2.5 2.6 -0.0 0.1\n 11 3 2.4 -0.5 0.0 0.0\n 11 4 -0.9 -0.4 -0.0 0.2\n 11 5 0.3 0.6 -0.1 -0.0\n 11 6 -0.7 -0.2 0.0 0.0\n 11 7 -0.1 -1.7 -0.0 0.1\n 11 8 1.4 -1.6 -0.1 -0.0\n 11 9 -0.6 -3.0 -0.1 -0.1\n 11 10 0.2 -2.0 -0.1 0.0\n 11 11 3.1 -2.6 -0.1 -0.0\n 12 0 -2.0 0.0 0.0 0.0\n 12 1 -0.1 -1.2 -0.0 -0.0\n 12 2 0.5 0.5 -0.0 0.0\n 12 3 1.3 1.3 0.0 -0.1\n 12 4 -1.2 -1.8 -0.0 0.1\n 12 5 0.7 0.1 -0.0 -0.0\n 12 6 0.3 0.7 0.0 0.0\n 12 7 0.5 -0.1 -0.0 -0.0\n 12 8 -0.2 0.6 0.0 0.1\n 12 9 -0.5 0.2 -0.0 -0.0\n 12 10 0.1 -0.9 -0.0 -0.0\n 12 11 -1.1 -0.0 -0.0 0.0\n 12 12 -0.3 0.5 -0.1 -0.1\n999999999999999999999999999999999999999999999999\n999999999999999999999999999999999999999999999999\n" + 'use strict'; + let modelLines = cof.split('\n'), + wmm = [], + i, vals, epoch, model, modelDate; + for (i in modelLines) { + if (modelLines.hasOwnProperty(i)) { + vals = modelLines[i].replace(/^\s+|\s+$/g, "").split(/\s+/); + if (vals.length === 3) { + epoch = parseFloat(vals[0]); + model = vals[1]; + modelDate = vals[2]; + } else if (vals.length === 6) { + wmm.push({ + n: parseInt(vals[0], 10), + m: parseInt(vals[1], 10), + gnm: parseFloat(vals[2]), + hnm: parseFloat(vals[3]), + dgnm: parseFloat(vals[4]), + dhnm: parseFloat(vals[5]) + }); + } + } + } + + return {epoch: epoch, model: model, modelDate: modelDate, wmm: wmm}; +} diff --git a/html/index.html b/html/index.html index 785b82c..ffed3ae 100644 --- a/html/index.html +++ b/html/index.html @@ -31,6 +31,7 @@ + tar1090 @@ -314,6 +315,38 @@ + +
+
+ Wind +
+
+
+
+
+ Speed: +
+
+ n/a +
+
+
+
+ Direction: +
+
+ n/a +
+
+
+
+ Crab angle: +
+
+ n/a +
+
+
Speed @@ -418,10 +451,10 @@
- Track Rate: + True Heading:
- n/a + n/a
@@ -434,18 +467,26 @@
- Roll: + Magnetic Decl.:
- n/a + n/a
- True Heading: + Track Rate:
- n/a + n/a +
+
+
+
+ Roll: +
+
+ n/a
diff --git a/html/script.js b/html/script.js index 2618d27..88e7efa 100644 --- a/html/script.js +++ b/html/script.js @@ -87,6 +87,7 @@ var labelFill = null; var blackFill = null; var labelStroke = null; var legSel = -1; +let geoMag = null; var shareLink = ''; @@ -1139,6 +1140,8 @@ function parse_history() { if (tempTrails) selectAllPlanes(); + + geoMag = geoMagFactory(cof2Obj()); } // Make a LineString with 'points'-number points @@ -1804,6 +1807,54 @@ function refreshSelected() { $('#selected_squawk2').text(selected.squawk); } + let magResult = null; + + if (geoMag && selected.position != null) { + let lon = selected.position[0]; + let lat = selected.position[1]; + let alt = selected.altitude == "ground" ? 0 : selected.altitude; + magResult = geoMag(lat, lon, alt); + $('#selected_mag_declination').text(format_track_brief(magResult.dec)); + } else { + $('#selected_mag_declination').text('n/a'); + } + + if (magResult && selected.gs != null && selected.tas != null && selected.track != null && selected.mag_heading != null) { + + const trk = (Math.PI / 180) * selected.track; + const hdg = (Math.PI / 180) * (selected.mag_heading + magResult.dec); + const tas = selected.tas; + const gs = selected.gs; + const ws = Math.round(Math.sqrt(Math.pow(tas - gs, 2) + 4 * tas * gs * Math.pow(Math.sin((hdg - trk) / 2), 2))); + let wd = trk + Math.atan2(tas * Math.sin(hdg - trk), tas * Math.cos(hdg - trk) - gs); + if (wd < 0) { + wd = wd + 2 * Math.PI; + } + if (wd > 2 * Math.PI) { + wd = wd - 2 * Math.PI; + } + wd = Math.round((180 / Math.PI) * wd); + $('#selected_wd').text(format_track_long(wd, DisplayUnits)); + $('#selected_ws').text(format_speed_long(ws, DisplayUnits)); + } else { + $('#selected_wd').text('n/a'); + $('#selected_ws').text('n/a'); + } + + if (magResult && selected.mag_heading != null && selected.track != null) + $('#selected_crab').text(format_track_brief(selected.mag_heading + magResult.dec - selected.track)); + else if (selected.true_heading != null && selected.track != null) + $('#selected_crab').text(format_track_brief(selected.true_heading - selected.track)); + else + $('#selected_crab').text('n/a'); + + $('#selected_mag_heading').text(format_track_brief(selected.mag_heading)); + + if (selected.true_heading == null && selected.mag_heading != null && magResult) + $('#selected_true_heading').text(format_track_brief(selected.mag_heading + magResult.dec)); + else + $('#selected_true_heading').text(format_track_brief(selected.true_heading)); + $('#selected_speed1').text(format_speed_long(selected.gs, DisplayUnits)); $('#selected_speed2').text(format_speed_long(selected.gs, DisplayUnits)); $('#selected_ias').text(format_speed_long(selected.ias, DisplayUnits)); @@ -1885,8 +1936,6 @@ function refreshSelected() { } $('#selected_altitude_geom').text(format_altitude_long(selected.alt_geom, selected.geom_rate, DisplayUnits)); - $('#selected_mag_heading').text(format_track_brief(selected.mag_heading)); - $('#selected_true_heading').text(format_track_brief(selected.true_heading)); $('#selected_ias').text(format_speed_long(selected.ias, DisplayUnits)); $('#selected_tas').text(format_speed_long(selected.tas, DisplayUnits)); if (selected.mach == null) {