@MikeG's comment pointed me in exactly the right direction, that of subsolar points. The Wikipedia page listed equations that came from a paper by NASA scientists. Following their equation as exactly as possible and adapting their FORTRAN code yields the following working code in TypeScript:
// Type for returning the apparent coordinates of the Sun's location in the sky
type SunLatLong = {
/** Longitude of the sun's apparent position in degrees */
latitude: number;
/** Latitude of the sun's apparent position in degrees */
longitude: number;
}
// Type for returning the unit vector pointing from the observer to the center of the Sun
type ObserverSunUnitVector = {
s_x: number;
s_y: number;
s_z: number;
}
// Type for returning the altitude and azimuth of the Sun for an observer
type SunAltAzi = {
/** Altitude of the Sun from the horizon in degrees from the observer's perspective */
altitude: number;
/** Azimuth of the Sun in degrees (south-clockwise convention) from the observer's perspective */
azimuth: number;
}
/**
* Returns the sun coordinates given a date.
*
* @remarks
* This method is adapted from the FORTRAN code in a paper by NASA workers Zhang, Stackhouse Jr., Macpherson, and Mikovitz, which can be found at https://doi.org/10.1016/j.renene.2021.03.047
*
* @param date - the UTC date for which to get relevant values. Uses JavaScript's Date object as type, but is interpreted as the count of milliseconds since the Unix epoch.
* @returns The coordinates (in degrees) of the sun in the sky using 0 N, 0 E as origin
*
* @beta
*/
export function date_to_sun_lat_long(date: Date) {
const degs = 180 / Math.PI
const rads = Math.PI / 180
const sundate = date || new Date(Date.now()) // UTC in milliseconds from Jan 1 1970
const jd = +sundate/86400000 + 2440587.5
const n = jd - 2451545.0 // fractional days since J2000
let L: number = (280.466 + 0.9856474 * n) % 360
let g: number = (357.528 + 0.9856003 * n) % 360
let lambda: number = (L + (1.915 * Math.sin(g * rads)) + 0.020 * Math.sin(2 * g * rads)) % 360
let epsilon: number = 23.440 - 0.0000004 * n
let alpha: number = (Math.atan2(Math.cos(epsilon * rads) * Math.sin(lambda * rads), Math.cos(lambda * rads)) * degs) % 360 // alpha in same quadrant as lambda
let delta: number = Math.asin(Math.sin(epsilon * rads) * Math.sin(lambda * rads)) * degs
let e_min: number = (((L - alpha) + 180) % 360) - 180 // in degrees
let T_gmt: number = sundate.getUTCHours() + (sundate.getUTCMinutes() / 60) + (sundate.getUTCSeconds() / 3600) // decimal hours since beginning of UTC day in question
let subsolar_point_lat: number = delta // Degrees
let subsolar_point_lon: number = (-15 * (T_gmt - 12 + (e_min * 4 / 60))) % 360 // Degrees
let SunLatLong: SunLatLong = {
latitude: subsolar_point_lat,
longitude: subsolar_point_lon,
}
return SunLatLong
}
/**
* Returns the sun coordinates, unit bector from observer to sun, and alitude and azimuth of the sun given a date and observer coordinates.
*
* @remarks
* This method is adapted from the FORTRAN code in a paper by NASA workers Zhang, Stackhouse Jr., Macpherson, and Mikovitz, which can be found at https://doi.org/10.1016/j.renene.2021.03.047
*
* @param date - the UTC date for which to get relevant values. Uses JavaScript's Date object as type, but is interpreted as the count of milliseconds since the Unix epoch.
* @param observer_lat - The latitude of the observer on Earth, measured in degrees
* @param observer_lon - The longitude of the observer on Earth, measured in degrees
* @returns (1) The coordinates (in degrees) of the sun in the sky using 0 N, 0 E as origin
* (2) The unit vector pointing from the observer to the center of the Sun
* (3) The altitude and azimuth of the sun from the observer's point of view
*
* @beta
*/
export function get_SunCoords_SunUnitVector_SunAziElev(date: Date, observer_lat: number, observer_lon: number) {
const degs = 180 / Math.PI
const rads = Math.PI / 180
const sundate = date || new Date(Date.now())
const sun_lat_long = date_to_sun_lat_long(sundate)
const phi_s = sun_lat_long.latitude
const lambda_s = sun_lat_long.longitude
const phi_o = observer_lat
const lambda_o = observer_lon
type AllValues = {
SunLatLong: SunLatLong,
ObserverSunUnitVector: ObserverSunUnitVector,
SunAltAzi: SunAltAzi
}
let ObserverSunUnitVector: ObserverSunUnitVector = {
s_x: Math.cos(phi_s * rads) * Math.sin((lambda_s * rads) - (lambda_o * rads)),
s_y: (Math.cos(phi_o * rads) * Math.sin(phi_s * rads)) - (Math.sin(phi_o * rads) * Math.cos(phi_s * rads) * Math.cos((lambda_s * rads) - (lambda_o * rads))),
s_z: (Math.sin(phi_o * rads) * Math.sin(phi_s * rads)) + (Math.cos(phi_o * rads) * Math.cos(phi_s * rads) * Math.cos((lambda_s * rads) - (lambda_o * rads)))
}
let SunAltAzi: SunAltAzi = {
altitude: Math.acos(ObserverSunUnitVector.s_z * rads) / rads,
azimuth: Math.atan2(-ObserverSunUnitVector.s_x * rads, -ObserverSunUnitVector.s_y * rads) / rads,
}
let AllVals: AllValues = {
SunLatLong: sun_lat_long,
ObserverSunUnitVector: ObserverSunUnitVector,
SunAltAzi: SunAltAzi,
}
return AllVals
}
The above gets me within 0.4 degrees longitude and exact latitude, which is accurate enough for my application. I suspect the longitude error may be due to a lack of specific leap year/leap second accounting, which the authors do in their code.
EDIT: After further research, I found the true solution to my original question of accounting for vernal equinox offset. The answer lies in using Sidereal Time. Many thanks to Greg Miller for his site's source code and his correspondence.
/**
* Return type with latitude and longitude coordinates (in both degrees and radians) of the subsolar point for a given time.
*
* @param latitude_rads Latitude (in radians) of the subsolar point
* @param longitude_rads Longitude (in radians) of the subsolar point
* @param latitude_degs Latitude (in degrees) of the subsolar point
* @param longitude_degs Longitude (in degrees) of the subsolar point
*
* @beta
*/
type LatLongFromRaDec = {
latitude_rads: number;
longitude_rads: number;
latitude_degs: number;
longitude_degs: number;
}
/**
* Return type with right ascension and declination of the Sun for a given time.
*
* @param right_ascension Right Ascension of the Sun for a given time
* @param declination Declination of the Sun for a given time
*
* @beta
*/
type RaDec = {
right_ascension: number;
declination: number;
}
/**
* Returns the right ascension and declination of the Sun given a date.
*
* @remarks
* This method borrows heavily/outright copies from Greg Miller's source code at https://www.celestialprogramming.com/snippets/geographicPosition.html
*
* @param date - the UTC date for which to get relevant values. Uses JavaScript's Date object as type, but is interpreted as the count of milliseconds since the Unix epoch.
* @returns {RaDec} The right ascension and declination of the Sun
*
* @beta
*/
function sunPosition(date: Date) {
const sundate = date || new Date(Date.now()) // UTC in milliseconds from Jan 1 1970
const jd = +sundate/86400000 + 2440587.5
let n: number = jd - 2451545.0;
let L: number = (280.466 + 0.9856474*n) % 360;
let g: number = ((357.528+ 0.9856003*n) % 360) * rads;
if( L < 0 ) { L += 360; }
if( g < 0 ) { g += Math.PI * 2.0; }
let lambda: number = (L + 1.915 * Math.sin(g) + 0.020 * Math.sin(2 * g)) * rads;
let beta: number = 0.0;
let epsilon: number = (23.440-0.0000004*n)*rads;
let right_ascension: number = Math.atan2(Math.cos(epsilon) * Math.sin(lambda), Math.cos(lambda));
let declination: number = Math.asin(Math.sin(epsilon) * Math.sin(lambda));
if( right_ascension < 0 ) { right_ascension += Math.PI * 2; }
let SunPosition: RaDec = {
right_ascension: right_ascension,
declination: declination,
}
return SunPosition;
}
/**
* Returns the geographic coordinates in degrees and radians of the subsolar point given a GMST date.
*
* @remarks
* This method borrows heavily/outright copies from Greg Miller's source code at https://www.celestialprogramming.com/snippets/geographicPosition.html
*
* @param right_ascension - The right ascension of the Sun.
* @param declination - The declination of the Sun.
* @param gmst - The Greenwich Mean Sidereal Time for the date in question.
* @returns {LatLongFromRaDec} The latitude and longitude in both degrees and radians of the subsolar point.
*
* @beta
*/
function getGeographicPosition(right_ascension: number, declination: number, gmst: number){
const latitude_rads = declination;
let longitude_rads = right_ascension - gmst;
if (longitude_rads > 2 * Math.PI){
longitude_rads -= 2 * Math.PI;
}
if (longitude_rads > Math.PI){
longitude_rads -= 2 * Math.PI;
}
if(longitude_rads < -Math.PI){
longitude_rads += 2 * Math.PI;
}
let LatLong: LatLongFromRaDec = {
latitude_rads: latitude_rads,
longitude_rads: longitude_rads,
latitude_degs: latitude_rads * degs,
longitude_degs: longitude_rads * degs,
}
return LatLong;
}
/**
* Returns the Earth Rotation Angle in radians given a GMST date.
*
* @remarks
* This method borrows heavily/outright copies from Greg Miller's source code at https://www.celestialprogramming.com/snippets/geographicPosition.html
*
* @param date - the UTC date for which to get relevant values. Uses JavaScript's Date object as type, but is interpreted as the count of milliseconds since the Unix epoch.
* @returns {LatLongFromRaDec} The latitude and longitude in both degrees and radians of the subsolar point.
*
* @beta
*/
function earthRotationAngle(date: Date){
const sundate = date || new Date(Date.now()) // UTC in milliseconds from Jan 1 1970
const jd = +sundate/86400000 + 2440587.5
const t = jd - 2451545.0;
const frac = jd % 1.0;
let era: number = ((Math.PI * 2) * (0.7790572732640 + 0.00273781191135448 * t + frac)) % (2 * Math.PI);
if( era < 0 ) { era += Math.PI * 2 };
return era;
}
/**
* Returns the Greenwich Mean Sidereal Time for a given UTC timestamp.
*
* @remarks
* This method borrows heavily/outright copies from Greg Miller's source code at https://www.celestialprogramming.com/snippets/geographicPosition.html
*
* @param date - the UTC date for which to get relevant values. Uses JavaScript's Date object as type, but is interpreted as the count of milliseconds since the Unix epoch.
* @returns gmst - The Greenwich Mean Sidereal Time for the given UTC timestamp.
*
* @beta
*/
function greenwichMeanSiderealTime(date: Date){
//The IAU Resolutions on Astronomical Reference Systems, Time Scales, and Earth Rotation Models Explanation and Implementation (George H. Kaplan)
//https://arxiv.org/pdf/astro-ph/0602086.pdf
const sundate = date || new Date(Date.now()) // UTC in milliseconds from Jan 1 1970
const jd = +sundate/86400000 + 2440587.5
const t = (+jd - 2451545.0) / 36525.0;
const era = earthRotationAngle(date);
// EQ 2.12
let gmst: number = (era + (0.014506 + (4612.15739966 * t) + (1.39667721 * t*t) + (-0.00009344 * t*t*t) + (0.00001882 * t*t*t*t)) / 60 / 60 * (Math.PI / 180)) % (2 * Math.PI);
if ( gmst < 0 ) { gmst += 2 * Math.PI };
return gmst;
}
/**
* Returns returns the latitude and longitude of the subsolar point at given a timestamp, in both degrees and radians..
*
* @remarks
* This method borrows heavily/outright copies from Greg Miller's source code at https://www.celestialprogramming.com/snippets/geographicPosition.html
*
* @param sundate - the UTC date for which to get relevant values. Uses JavaScript's Date object as type, but is interpreted as the count of milliseconds since the Unix epoch.
* @returns {LatLongFromRaDec} Latitude (in both degrees and radians) and longitude (in both degrees and radians) of the subsolar point.
*
* @beta
*/
export function computeSubsolarPointCoordinatesWithGMST(sundate: Date) {
const sun = sunPosition(sundate);
const gmst = greenwichMeanSiderealTime(sundate);
let gp = getGeographicPosition(sun.right_ascension, sun.declination, gmst);
return gp;
}
I hope this can be of help to anybody who finds themselves in the same spot I was in.
