This commit is contained in:
2026-10-21 19:32:02 -04:00
commit b74d278c8e
96 changed files with 54476 additions and 0 deletions
+46
View File
@@ -0,0 +1,46 @@
var GeodeticCoord = require('./geodetic.js');
module.exports = CartesianCoord;
function CartesianCoord(options){
this.x = options.x;
this.y = options.y;
this.z = options.z;
}
CartesianCoord.prototype.toGeodetic = function(ellipsoid){
var modified_b = this.z < 0 ? -ellipsoid.b : ellipsoid.b;
var x = this.x;
var y = this.y;
var z = this.z;
var r = Math.sqr(x*x + y*y);
var e = (modified_b*z - (ellipsoid.a*ellipsoid.a - modified_b*modified_b)) / (ellipsoid.a*r);
var f = (modified_b*z + (ellipsoid.a*ellipsoid.a - modified_b*modified_b)) / (ellipsoid.a*r);
var p = (4/3) * (e*f+1);
var q = 2 * (e*e + f*f);
var d = p*p*p + q*q;
var v;
if(d > 0) {
v = Math.pow(Math.sqrt(d) - q, 1/3) - Math.pow(Math.sqrt(d) + q, 1/3);
} else {
v = 2 * Math.sqrt(-p) * Math.cos(Math.acos(q / (p * Math.sqrt(-p)) ) / 3);
}
if(v*v < Math.abs(p)) {
v = -(v*v*v + 2*q) / 3*p;
}
var g = (Math.sqrt(e*e + v) + e) / 2;
var t = Math.sqrt(g*g + (f - v*g)/(2*g - e)) - g;
var rlat = Math.atan( (ellipsoid.a*(1 - t*t)) / (2*modified_b*t) );
var zlong = Math.atan2(y, x);
if(zlong < 0) zlong += 2*Math.PI;
var height_above_ellipsoid = (r - ellipsoid.a*t) * Math.cos(rlat) + (z - modified_b) * Math.sin(rlat);
var lambda = zlong * Math.PI/180;
if(lambda > 180) lambda -= 360;
return new GeodeticCoord({
lambda: lambda,
phi: rlat * Math.PI/180,
height_above_ellipsoid: height_above_ellipsoid,
ellipsoid: ellipsoid
});
};
+31
View File
@@ -0,0 +1,31 @@
var SphericalCoord = require('./spherical.js');
var Ellipsoid = require('../ellipsoid.js');
module.exports = GeodeticCoord;
function GeodeticCoord(options){
this.options = options;
this.lambda = options.lambda || options.lon; // longitude
this.phi = options.phi || options.lat; // geodetic latitude
this.height_above_ellipsoid = options.height_above_ellipsoid || 0;
this.ellipsoid = options.ellipsoid || new Ellipsoid();
}
GeodeticCoord.prototype.toSpherical = function() {
var ellipsoid = this.ellipsoid;
var coslat = Math.cos(this.phi*Math.PI/180);
var sinlat = Math.sin(this.phi*Math.PI/180);
var rc = ellipsoid.a / Math.sqrt(1 - ellipsoid.epssq * sinlat * sinlat);
var xp = (rc + this.height_above_ellipsoid) * coslat;
var zp = (rc * (1 - ellipsoid.epssq) + this.height_above_ellipsoid) * sinlat;
var r = Math.sqrt(xp*xp + zp*zp);
return new SphericalCoord({
r: r,
phig: 180/Math.PI*Math.asin(zp / r),
lambda: this.lambda
});
};
GeodeticCoord.prototype.clone = function(){
return new GeodeticCoord(this.options);
};
+5
View File
@@ -0,0 +1,5 @@
module.exports = {
Cartesian: require('./cartesian.js'),
Spherical: require('./spherical.js'),
Geodetic: require('./geodetic.js')
};
+54
View File
@@ -0,0 +1,54 @@
var CartesianCoord = require('./cartesian.js');
var LegendreFunction = require('../legendre.js');
module.exports = SphericalCoord;
function SphericalCoord(options){
this.lambda = options.lambda; // longitude
this.phig = options.phig; // geocentric latitude
this.r = options.r; // distance from center of ellipsoid
}
SphericalCoord.prototype.toCartesian = function(){
var radphi = this.phig * Math.PI/180;
var radlambda = this.lambda * Math.PI/180;
return new CartesianCoord({
x: this.r * cos(radphi) * cos(radlambda),
y: this.r * cos(radphi) * sin(radlambda),
z: this.r * sin(radphi)
});
};
SphericalCoord.prototype.toGeodetic = function(ellipsoid){
return this.toCartesian().toGeodetic(ellipsoid);
};
SphericalCoord.prototype.getHarmonicVariables = function(ellipsoid, n_max){
var m, n;
var cos_lambda = Math.cos(Math.PI/180 * this.lambda);
var sin_lambda = Math.sin(Math.PI/180 * this.lambda);
var cos_mlambda = [1.0, cos_lambda];
var sin_mlambda = [0.0, sin_lambda];
var relative_radius_power = [
(ellipsoid.re/this.r) * (ellipsoid.re/this.r)
];
for(n = 1; n <= n_max; n++){
relative_radius_power[n] = relative_radius_power[n-1] * (ellipsoid.re/this.r);
}
for(m = 2; m <= n_max; m++){
cos_mlambda[m] = cos_mlambda[m-1]*cos_lambda - sin_mlambda[m-1]*sin_lambda;
sin_mlambda[m] = cos_mlambda[m-1]*sin_lambda + sin_mlambda[m-1]*cos_lambda;
}
return {
relative_radius_power: relative_radius_power,
cos_mlambda: cos_mlambda,
sin_mlambda: sin_mlambda
};
};
SphericalCoord.prototype.getLegendreFunction = function(n_max) {
return new LegendreFunction(this, n_max);
};
+17
View File
@@ -0,0 +1,17 @@
module.exports = Ellipsoid;
function Ellipsoid(options) {
options = options || {
a: 6378.137,
b: 6356.7523142,
fla: 1 / 298.257223563,
re: 6371.2
};
this.a = options.a; // semi-major axis
this.b = options.b; // semi-minor axis
this.fla = options.fla; // flattening
this.re = options.re; // mean radius
this.eps = Math.sqrt(1-(this.b*this.b)/(this.a*this.a)); // first eccentricity
this.epssq = this.eps*this.eps; // first eccentricity squared
}
+168
View File
@@ -0,0 +1,168 @@
module.exports = LegendreFunction;
function LegendreFunction(coord_spherical, n_max){
var sin_phi = Math.sin(Math.PI/180 * coord_spherical.phig);
var result;
if(n_max <= 16 || (1 - Math.abs(sin_phi)) < 1e-10){
result = PcupLow(sin_phi, n_max);
} else {
result = PcupHigh(sin_phi, n_max);
}
this.pcup = result.pcup;
this.dpcup = result.dpcup;
}
function PcupLow(x, n_max){
var k, z, n, m, i, i1, i2, num_terms;
var schmidt_quasi_norm = [1.0];
var pcup = [1.0];
var dpcup = [0.0];
z = Math.sqrt((1-x)*(1+x));
num_terms = (n_max+1) * (n_max+2) / 2;
for(n = 1; n <= n_max; n++) {
for(m = 0; m <= n; m++) {
i = n*(n+1)/2 + m;
if(n == m){
i1 = (n-1)*n/2 + m - 1;
pcup[i] = z*pcup[i1];
dpcup[i] = z*dpcup[i1] + x*pcup[i1];
} else if (n == 1 && m == 0) {
i1 = (n-1)*n/2 + m;
pcup[i] = x*pcup[i1];
dpcup[i] = x*dpcup[i1] - z*pcup[i1];
} else if (n > 1 && n != m) {
i1 = (n-2)*(n-1)/2 + m;
i2 = (n-1)*n/2 + m;
if(m > n - 2){
pcup[i] = x*pcup[i2];
dpcup[i] = x*dpcup[i2] - z*pcup[i2];
} else {
k = ((n-1)*(n-1) - m*m) / ((2*n-1)*(2*n-3));
pcup[i] = x*pcup[i2] - k*pcup[i1];
dpcup[i] = x*dpcup[i2] - z*pcup[i2] - k*dpcup[i1];
}
}
}
}
for(n = 1; n <= n_max; n++) {
i = n*(n+1)/2;
i1 = (n-1)*n/2;
schmidt_quasi_norm[i] = schmidt_quasi_norm[i1] * (2*n-1) / n;
for(m = 1; m <= n; m++){
i = n*(n+1)/2 + m;
i1 = n*(n+1)/2 + m - 1;
schmidt_quasi_norm[i] = schmidt_quasi_norm[i1] * Math.sqrt(((n - m + 1) * (m == 1 ? 2 : 1)) / (n + m));
}
}
for(n = 1; n <= n_max; n++) {
for(m = 0; m <= n; m++) {
i = n*(n+1)/2+m;
pcup[i] *= schmidt_quasi_norm[i];
dpcup[i] *= -schmidt_quasi_norm[i];
}
}
return {
pcup: pcup,
dpcup: dpcup
};
}
function PcupHigh(x, n_max){
if(Math.abs(x) == 1.0){
throw new Error("Error in PcupHigh: derivative cannot be calculated at poles");
}
var n, m, k;
var num_terms = (n_max + 1) * (n_max + 2)/2;
var f1 = [];
var f2 = [];
var pre_sqr = [];
var scalef = 1.0e-280;
for(n = 0; n <= 2*n_max + 1; ++n){
pre_sqr[n] = Math.sqrt(n);
}
k = 2;
for(n = 0; n <= n_max; n++){
k++;
f1[k] = (2*n - 1) / n;
f2[k] = (n - 1) / n;
for(m = 1; m <= n - 2; m++){
k++;
f1[k] = (2*n - 1) / pre_sqr[n+m] / pre_sqr[n-m];
f2[k] = pre_sqr[n-m-1] * pre_sqr[n+m-1] / pre_sqr[n+m] / pre_sqr[n-m];
}
k += 2;
}
var z = Math.sqrt((1-x)*(1+x));
var plm;
var pm1 = x;
var pm2 = 1;
var pcup = [1.0, pm1];
var dpcup = [0.0, z];
if(n_max == 0){
throw new Error("Error in PcupHigh: n_max must be greater than 0");
}
k = 1;
for(n = 2; n <= n_max; n++)
{
k = k + n;
plm = f1[k] * x * pm1 - f2[k] * pm2;
pcup[k] = plm;
dpcup[k] = n * (pm1 - x * plm) / z;
pm2 = pm1;
pm1 = plm;
}
var pmm = pre_sqr[2] * scalef;
var rescalem = 1/scalef;
var kstart = 0;
for(var m = 1; m <= n_max - 1; ++m){
rescalem *= z;
//calculate pcup(m,m)
kstart = kstart + m + 1;
pmm = pmm * pre_sqr[2*m+1] / pre_sqr[2*m];
pcup[kstart] = pmm * rescalem / pre_sqr[2*m+1];
dpcup[kstart] = -( m*x*pcup[kstart] /z );
pm2 = pmm / pre_sqr[2*m + 1];
//calculate pcup(m+1,m)
k = kstart + m + 1;
pm1 = x * pre_sqr[2*m+1] * pm2;
pcup[k] = pm1*rescalem;
dpcup[k] = (pm2*rescalem *pre_sqr[2*m+1] - x*(m+1)*pcup[k]) / z;
//calculate pcup(n,m)
for(n = m + 2; n <= n_max; ++n){
k = k + n;
plm = x*f1[k]*pm1 - f2[k]*pm2;
pcup[k] = plm*rescalem;
dpcup[k] = (pre_sqr[n+m]*pre_sqr[n-m]*pm1*rescalem - n*x*pcup[k]) / z;
pm2 = pm1;
pm1 = plm;
}
}
//calculate pcup(n_max,n_max)
rescalem = rescalem*z;
kstart = kstart + m + 1;
pmm = pmm / pre_sqr[2*n_max];
pcup[kstart] = pmm * rescalem;
dpcup[kstart] = -n_max * x * pcup[kstart] / z;
return {
pcup: pcup,
dpcup: dpcup
};
}
+98
View File
@@ -0,0 +1,98 @@
module.exports = MagneticElements;
function MagneticElements(options){
this.x = options.x;
this.y = options.y;
this.z = options.z;
this.h = options.h;
this.f = options.f;
this.incl = options.incl;
this.decl = options.decl;
this.gv = options.gv;
this.xdot = options.xdot;
this.ydot = options.ydot;
this.zdot = options.zdot;
this.hdot = options.hdot;
this.fdot = options.fdot;
this.incldot = options.incldot;
this.decldot = options.decldot;
this.gvdot = options.gvdot;
}
MagneticElements.prototype.clone = function(){
return new MagneticElements(this);
};
MagneticElements.prototype.scale = function(factor){
return new MagneticElements({
x : this.x * factor,
y : this.y * factor,
z : this.z * factor,
h : this.h * factor,
f : this.f * factor,
incl : this.incl * factor,
decl : this.decl * factor,
gv : this.gv * factor,
xdot : this.xdot * factor,
ydot : this.ydot * factor,
zdot : this.zdot * factor,
hdot : this.hdot * factor,
fdot : this.fdot * factor,
incldot : this.incldot * factor,
decldot : this.decldot * factor,
gvdot : this.gvdot * factor
});
};
MagneticElements.prototype.subtract = function(subtrahend){
return new MagneticElements({
x : this.x - subtrahend.x,
y : this.y - subtrahend.y,
z : this.z - subtrahend.z,
h : this.h - subtrahend.h,
f : this.f - subtrahend.f,
incl : this.incl - subtrahend.incl,
decl : this.decl - subtrahend.decl,
gv : this.gv - subtrahend.gv,
xdot : this.xdot - subtrahend.xdot,
ydot : this.ydot - subtrahend.ydot,
zdot : this.zdot - subtrahend.zdot,
hdot : this.hdot - subtrahend.hdot,
fdot : this.fdot - subtrahend.fdot,
incldot : this.incldot - subtrahend.incldot,
decldot : this.decldot - subtrahend.decldot,
gvdot : this.gvdot - subtrahend.gvdot
});
};
MagneticElements.fromGeoMagneticVector = function(magnetic_vector){
var bx = magnetic_vector.bx;
var by = magnetic_vector.by;
var bz = magnetic_vector.bz;
var h = Math.sqrt(bx*bx + by*by);
return new MagneticElements({
x: bx,
y: by,
z: bz,
h: h,
f: Math.sqrt(h*h + bz*bz),
decl: 180/Math.PI * Math.atan2(by, bx),
incl: 180/Math.PI * Math.atan2(bz, h)
});
};
MagneticElements.fromSecularVariationVector = function(magnetic_variation){
var bx = magnetic_variation.bx;
var by = magnetic_variation.by;
var bz = magnetic_variation.bz;
return new MagneticElements({
xdot: bx,
ydot: by,
zdot: bz,
hdot: (this.x*bx + this.y*by) / this.h,
fdot: (this.x*bx + this.y*by + this.z*bz) / this.f,
decldot: 180/Math.PI * (this.x*this.ydot - this.y*this.xdot) / (this.h*this.h),
incldot: 180/Math.PI * (this.h*this.zdot - this.z*this.hdot) / (this.f*this.f),
gvdot: this.decldot
});
};
+16
View File
@@ -0,0 +1,16 @@
module.exports = MagneticVector;
function MagneticVector(options){
this.bx = options.bx;
this.by = options.by;
this.bz = options.bz;
}
MagneticVector.prototype.rotate = function(coord_spherical, coord_geodetic) {
var psi = Math.PI/180 * (coord_spherical.phig - coord_geodetic.phi);
return new MagneticVector({
bz: this.bx*Math.sin(psi) + this.bz*Math.cos(psi),
bx: this.bx*Math.cos(psi) - this.bz*Math.sin(psi),
by: this.by
});
};
+140
View File
@@ -0,0 +1,140 @@
var GeodeticCoord = require('./coords/geodetic.js');
var MagneticElements = require('./magnetic_elements.js');
var MagneticVector = require('./magnetic_vector.js');
module.exports = Model;
function Model(options){
this.epoch = options.epoch;
this.start_date = new Date(options.start_date);
this.end_date = new Date(options.end_date);
this.name = options.name || "";
this.main_field_coeff_g = options.main_field_coeff_g || [0];
this.main_field_coeff_h = options.main_field_coeff_h || [0];
this.secular_var_coeff_g = options.secular_var_coeff_g || [0];
this.secular_var_coeff_h = options.secular_var_coeff_h || [0];
this.n_max = options.n_max || 0;
this.n_max_sec_var = options.n_max_sec_var || 0;
this.num_terms = this.n_max*(this.n_max+1)/2 + this.n_max;
}
Model.prototype.getTimedModel = function(date, allowOutOfBoundsModel = false){
var year_int = date.getUTCFullYear();
var fractional_year = (date.valueOf() - Date.UTC(year_int)) / (1000*3600*24*365);
var year = year_int + fractional_year;
var dyear = year - this.epoch;
if(date < this.start_date || date > this.end_date){
if(allowOutOfBoundsModel) {
console.error("Model is only valid from "+this.start_date.toDateString()+" to "+this.end_date.toDateString());
} else {
throw new Error("Model is only valid from "+this.start_date.toDateString()+" to "+this.end_date.toDateString());
}
}
var model = new Model({
epoch: this.epoch,
n_max: this.n_max,
n_max_sec_var: this.n_max_sec_var,
name: this.name,
start_date: this.start_date,
end_date: this.end_date,
});
var a = model.n_max_sec_var;
var b = a*(a + 1)/2 + a;
for(var n = 1; n <= this.n_max; n++){
for(var m = 0; m <= n; m++){
var i = n * (n + 1)/2 + m;
var hnm = this.main_field_coeff_h[i];
var gnm = this.main_field_coeff_g[i];
var dhnm = this.secular_var_coeff_h[i];
var dgnm = this.secular_var_coeff_g[i];
if(i <= b){
model.main_field_coeff_h[i] = hnm + dyear*dhnm;
model.main_field_coeff_g[i] = gnm + dyear*dgnm;
model.secular_var_coeff_h[i] = dhnm;
model.secular_var_coeff_g[i] = dgnm;
} else {
model.main_field_coeff_h[i] = hnm;
model.main_field_coeff_g[i] = gnm;
}
}
}
return model;
};
Model.prototype.getSummation = function(legendre, sph_variables, coord_spherical){
var bx = 0;
var by = 0;
var bz = 0;
var n, m, i, k;
var relative_radius_power = sph_variables.relative_radius_power;
var cos_mlambda = sph_variables.cos_mlambda;
var sin_mlambda = sph_variables.sin_mlambda;
var g = this.main_field_coeff_g;
var h = this.main_field_coeff_h;
for(n = 1; n <= this.n_max; n++){
for(m = 0; m <= n; m++){
i = n*(n+1)/2 + m;
bz -= relative_radius_power[n] *
(g[i]*cos_mlambda[m] + h[i]*sin_mlambda[m]) *
(n+1) * legendre.pcup[i];
by += relative_radius_power[n] *
(g[i]*sin_mlambda[m] - h[i]*cos_mlambda[m]) *
(m) * legendre.pcup[i];
bx -= relative_radius_power[n] *
(g[i]*cos_mlambda[m] + h[i]*sin_mlambda[m]) *
legendre.dpcup[i];
}
}
var cos_phi = Math.cos(Math.PI/180 * coord_spherical.phig);
if(Math.abs(cos_phi) > 1e-10){
by = by / cos_phi;
} else {
//special calculation around poles
by = 0;
var schmidt_quasi_norm1 = 1.0, schmidt_quasi_norm2, schmidt_quasi_norm3;
var pcup_s = [1];
var sin_phi = Math.sin(Math.PI/180 * coord_spherical.phig);
for(n = 1; n <= this.n_max; n++){
i = n*(n+1)/2 + 1;
schmidt_quasi_norm2 = schmidt_quasi_norm1 * (2*n-1) / n;
schmidt_quasi_norm3 = schmidt_quasi_norm2 * Math.sqrt(2*n/(n+1));
schmidt_quasi_norm1 = schmidt_quasi_norm2;
if(n == 1){
pcup_s[n] = pcup_s[n-1];
} else {
k = ((n-1)*(n-1) - 1) / ((2*n-1)*(2*n-3));
pcup_s[n] = sin_phi * pcup_s[n-1] - k*pcup_s[n-2];
}
by += relative_radius_power[n] *
(g[i]*sin_mlambda[1] - h[i]*cos_mlambda[1]) *
pcup_s[n] * schmidt_quasi_norm3;
}
}
return new MagneticVector({
bx: bx, by: by, bz: bz
});
};
Model.prototype.point = function(coords) {
// Extract altitude if provided, otherwise default to 0
var altitude = (coords.length > 2 && typeof coords[2] === 'number') ? coords[2] : 0;
var coord_geodetic = new GeodeticCoord({
lat: coords[0],
lon: coords[1],
height_above_ellipsoid: altitude
});
var coord_spherical = coord_geodetic.toSpherical();
var legendre = coord_spherical.getLegendreFunction(this.n_max);
var harmonic_variables = coord_spherical.getHarmonicVariables(coord_geodetic.ellipsoid, this.n_max);
var magnetic_vector_sph = this.getSummation(legendre, harmonic_variables, coord_spherical);
var magnetic_vector_geo = magnetic_vector_sph.rotate(coord_spherical, coord_geodetic);
return MagneticElements.fromGeoMagneticVector(magnetic_vector_geo);
};
+40
View File
@@ -0,0 +1,40 @@
module.exports = function(str){
var result = {
main_field_coeff_g: [0],
main_field_coeff_h: [0],
secular_var_coeff_g: [0],
secular_var_coeff_h: [0],
n_max: 0,
n_max_sec_var: 0
};
var lines = str.split(/\r?\n/);
for(var i = 0; i < lines.length; i++) {
var line = lines[i];
var vals = line.match(/\S+/g);
if(!vals) continue;
if(vals.length === 3){
result.epoch = parseFloat(vals[0]);
result.name = vals[1];
result.start_date = new Date(vals[2]);
result.end_date = new Date(vals[2]);
var year = result.start_date.getFullYear();
result.end_date.setFullYear(year + 5);
}
if(vals.length === 6){
var n = parseInt(vals[0]);
var m = parseInt(vals[1]);
if(m <= n) {
var i = (n*(n+1)/2+m);
result.main_field_coeff_g[i] = parseFloat(vals[2]);
result.main_field_coeff_h[i] = parseFloat(vals[3]);
result.secular_var_coeff_g[i] = parseFloat(vals[4]);
result.secular_var_coeff_h[i] = parseFloat(vals[5]);
}
if(n > result.n_max){
result.n_max = n;
result.n_max_sec_var = n;
}
}
}
return result;
};