City-Lifeline-Cockpit/tools/city-massing/verify-roof-geodesy.mjs

import assert from 'node:assert/strict'
import {readFileSync} from 'node:fs'
import {enuToGeodetic} from './roof-geodesy.mjs'
const identity=[1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1]
// Analytic WGS84 forward coordinates are independent of the Bowring inverse.
let maxAngleError=0,maxHeightError=0
for(const lon of [110.32875,110.395,110.46578])for(const lat of [28.39172,28.454,28.4913])for(const height of [0,150,450]){
 const L=lon*Math.PI/180,P=lat*Math.PI/180,a=6378137,e2=6.6943799901413165e-3,N=a/Math.sqrt(1-e2*Math.sin(P)**2)
 const X=(N+height)*Math.cos(P)*Math.cos(L),Y=(N+height)*Math.cos(P)*Math.sin(L),Z=(N*(1-e2)+height)*Math.sin(P)
 const actual=enuToGeodetic(identity,X,Y,Z)
 maxAngleError=Math.max(maxAngleError,Math.abs(actual.longitude-L),Math.abs(actual.latitude-P));maxHeightError=Math.max(maxHeightError,Math.abs(actual.height-height))
 assert.ok(Math.abs(actual.latitude-P)<1e-10&&Math.abs(actual.longitude-L)<1e-10)
 assert.ok(Math.abs(actual.height-height)<.001)
}
const input=JSON.parse(readFileSync('tools/city-massing/city/cad-bim/instanced-active-input.json'))
const origin=enuToGeodetic(input.transform,0,0,0)
assert.ok(Math.abs(origin.longitude*180/Math.PI-input.origin.lon)<1e-8&&Math.abs(origin.latitude*180/Math.PI-input.origin.lat)<1e-8)
for(const point of [[1400,-1200,150],[3000,-6000,200],[-3000,3000,100]]){
 const m=input.transform,[x,y,z]=point,X=m[0]*x+m[4]*y+m[8]*z+m[12],Y=m[1]*x+m[5]*y+m[9]*z+m[13],Z=m[2]*x+m[6]*y+m[10]*z+m[14]
 assert.deepEqual(enuToGeodetic(m,x,y,z),enuToGeodetic(identity,X,Y,Z))
}
console.log(JSON.stringify({passed:true,analyticWgs84RoundTrips:27,maxAngleErrorRad:maxAngleError,maxHeightErrorM:maxHeightError,transformChecked:true}))