<p>@MikeG's comment pointed me in exactly the right direction, that of <a href="https://en.wikipedia.org/wiki/Subsolar_point" rel="nofollow noreferrer">subsolar points</a>. The Wikipedia page listed equations that came from <a href="https://www.sciencedirect.com/science/article/pii/S0960148121004031?via%3Dihub" rel="nofollow noreferrer">a paper by NASA scientists</a>. Following their equation as exactly as possible and adapting their FORTRAN code yields the following working code in TypeScript:</p>
<pre><code>// 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
}
</code></pre>
<p>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.</p>
<p><strong>EDIT:</strong> 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 <a href="https://www.celestialprogramming.com/snippets/geographicPosition.html" rel="nofollow noreferrer">his site</a>'s source code and his correspondence.</p>
<pre><code>/** 
  * 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 &lt; 0 ) { L += 360; }
    if( g &lt; 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 &lt; 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 &gt; 2 * Math.PI){
    longitude_rads -= 2 * Math.PI;
  }
  if (longitude_rads &gt; Math.PI){
    longitude_rads -= 2 * Math.PI;
  }
  if(longitude_rads &lt; -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 &lt; 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 &lt; 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;
}
</code></pre>
<p>I hope this can be of help to anybody who finds themselves in the same spot I was in.</p>

CC BY-SA 4.0
2023-09-08T15:53:32.640Z
2023-09-12T23:02:26.623Z
2023-09-12T23:02:26.623Z
Michael Bonnet
<p>I am trying to calculate the position of the sun, in GPS coordinates terms, from a given date. As of right now, I know how to accurately calculate the equatorial Right Ascension and Declination, and even Sun Position Vectors. However, these coordinates systems are based on the vernal equinox, which I can't seem to figure out how to calculate, or even find an equation for.</p>
<p>My gameplan is to take the offset in degrees (or whatever transformation) of the equatorial system from the GPS frame, and add that offset my equatorial terms to find the position of the sun in the terms I'm looking for.</p>
<p>Astropy isn't an option for me - I need the actual pure math to then implement myself in a different language without such good scientific support.</p>

CC BY-SA 4.0
2023-09-07T02:22:23.187Z
Hi and welcome to AstronomySE. I’m also a “purist” and prefer pure math to Mathematica, Wolfram, or Python. I’m not sure I understand your question, though. You want to know how to calculate the position of the vernal point… relative to *what* exactly? It is the reference point itself, so unless you want to know where it is from a given point on the Earth (should I understand the position you give as Earth-based?), or to a set of distant celestial objects that are “fixed” (such as quasars), then the question makes no sense…
Pierre Paquette
CC BY-SA 4.0
2023-09-07T02:36:54.390Z
I suppose I want to know where it is from a given point on the Earth, with that point being Earth-based - specifically 0°00'00.0"N 0°00'00.0"E. or 0.000000, 0.000000 in GPS coordinates. This is because that point is the "origin" in the project I'm working on, and need to know the very rough position of the Sun relative to that origin point. For example, if at a certain datetime, the sun is directly overhead for an observer at 0 N, 45 W, the "coordinates" I'm looking to calculate would be 0.00, -45.00.
Michael Bonnet
CC BY-SA 4.0
2023-09-07T02:51:24.260Z
Would the [altitude and azimuth](https://en.wikipedia.org/wiki/Horizontal_coordinate_system) of the Sun meet your need? If not, why not?
Mike G
CC BY-SA 4.0
2023-09-07T02:59:27.853Z
@MikeG altitude and azimuth would not meet my need because as I understand it, it is observer-based and also reliant upon vernal equinox - I need to convert from such values to something akin to ECI/ECEF, as I need to have the "location" of the sun represented akin to how satellites might be represented in their movement around the Earth, as a point in space relative to the Earth. I am attempting to model what parts of the globe are illuminated by the sun at arbitrary points in time.
Michael Bonnet
CC BY-SA 4.0
2023-09-07T03:34:29.010Z
Are you looking for the [subsolar point](https://en.wikipedia.org/wiki/Subsolar_point)?
Mike G
CC BY-SA 4.0
2023-09-07T03:58:08.943Z
It kinda sounds like you want "Greenwich sidereal time", it relates the Earth's rotation in relation to the equinox.  This link shows how to use it to relate any RA/Dec to lat/lon. https://www.celestialprogramming.com/snippets/geographicPosition.html
Greg Miller
CC BY-SA 4.0
2023-09-07T08:07:35.080Z
*"these coordinates systems are based on the vernal equinox"*. There are a few options. Coordinates of a body at a particular time can be given in reference to the location of the equinox at that time (aka equinox of date), or more commonly in reference to a standard equinox at a particular epoch (eg, J2000.0 uses the equinox of Noon 1 Jan 2000). But the modern approach uses the [ICRF](https://en.wikipedia.org/wiki/International_Celestial_Reference_System_and_its_realizations), which is essentially aligned to the J2000.0 equatorial frame.
PM 2Ring
CC BY-SA 4.0
2023-09-07T16:06:54.147Z
Possible duplicate of https://astronomy.stackexchange.com/questions/50257/how-do-i-calculate-the-closest-point-to-the-sun-on-earths-surface
JohnHoltz
CC BY-SA 4.0
2023-09-07T17:32:09.460Z
The link posted by @MikeG was what I was looking for, the subsolar point. The equation in the Wikipedia article works, so long as you follow their exact method in the referenced paper, and note that the T(GMT) term is decimal hours in the given day.
Michael Bonnet
CC BY-SA 4.0
2023-09-08T01:16:09.883Z
I “third” Mike G and PM 2Ring. Post your solution (with as many details as possible) as an answer, please.
Pierre Paquette
CC BY-SA 4.0
2023-09-08T15:54:01.800Z
Solution removed from post and put into an answer.
Michael Bonnet
<the-sun><orbit><positional-astronomy>
Calculating Position of Sun in GPS Frame Relative to 0°00'00.0"N 0°00'00.0"E: Determining Vernal Equinox Offset
Michael Bonnet
# Calculating Position of Sun in GPS Frame Relative to 0°00'00.0"N 0°00'00.0"E: Determining Vernal Equinox Offset

(Asked by: username_0 on 2023-09-07)

I am trying to calculate the position of the sun, in GPS coordinates terms, from a given date. As of right now, I know how to accurately calculate the equatorial Right Ascension and Declination, and even Sun Position Vectors. However, these coordinates systems are based on the vernal equinox, which I can't seem to figure out how to calculate, or even find an equation for.

My gameplan is to take the offset in degrees (or whatever transformation) of the equatorial system from the GPS frame, and add that offset my equatorial terms to find the position of the sun in the terms I'm looking for.

Astropy isn't an option for me - I need the actual pure math to then implement myself in a different language without such good scientific support.

**Tags:** the-sun, orbit, positional-astronomy
> **Comments:**
> **username_1**: Hi and welcome to AstronomySE. I’m also a “purist” and prefer pure math to Mathematica, Wolfram, or Python. I’m not sure I understand your question, though. You want to know how to calculate the position of the vernal point… relative to *what* exactly? It is the reference point itself, so unless you want to know where it is from a given point on the Earth (should I understand the position you give as Earth-based?), or to a set of distant celestial objects that are “fixed” (such as quasars), then the question makes no sense…  
> **username_0**: I suppose I want to know where it is from a given point on the Earth, with that point being Earth-based - specifically 0°00'00.0"N 0°00'00.0"E. or 0.000000, 0.000000 in GPS coordinates. This is because that point is the "origin" in the project I'm working on, and need to know the very rough position of the Sun relative to that origin point. For example, if at a certain datetime, the sun is directly overhead for an observer at 0 N, 45 W, the "coordinates" I'm looking to calculate would be 0.00, -45.00.  
> **username_2**: Would the [altitude and azimuth](https://en.wikipedia.org/wiki/Horizontal_coordinate_system) of the Sun meet your need? If not, why not?  
> **username_0**: @MikeG altitude and azimuth would not meet my need because as I understand it, it is observer-based and also reliant upon vernal equinox - I need to convert from such values to something akin to ECI/ECEF, as I need to have the "location" of the sun represented akin to how satellites might be represented in their movement around the Earth, as a point in space relative to the Earth. I am attempting to model what parts of the globe are illuminated by the sun at arbitrary points in time.  
> **username_2**: Are you looking for the [subsolar point](https://en.wikipedia.org/wiki/Subsolar_point)?  
> **username_3**: It kinda sounds like you want "Greenwich sidereal time", it relates the Earth's rotation in relation to the equinox.  This link shows how to use it to relate any RA/Dec to lat/lon. https://www.celestialprogramming.com/snippets/geographicPosition.html  
> **username_4**: *"these coordinates systems are based on the vernal equinox"*. There are a few options. Coordinates of a body at a particular time can be given in reference to the location of the equinox at that time (aka equinox of date), or more commonly in reference to a standard equinox at a particular epoch (eg, J2000.0 uses the equinox of Noon 1 Jan 2000). But the modern approach uses the [ICRF](https://en.wikipedia.org/wiki/International_Celestial_Reference_System_and_its_realizations), which is essentially aligned to the J2000.0 equatorial frame.  
> **username_5**: Possible duplicate of https://astronomy.stackexchange.com/questions/50257/how-do-i-calculate-the-closest-point-to-the-sun-on-earths-surface  
> **username_0**: The link posted by @MikeG was what I was looking for, the subsolar point. The equation in the Wikipedia article works, so long as you follow their exact method in the referenced paper, and note that the T(GMT) term is decimal hours in the given day.  
> **username_1**: I “third” Mike G and PM 2Ring. Post your solution (with as many details as possible) as an answer, please.  
> **username_0**: Solution removed from post and put into an answer.  

---

## 1 answers

### Answer by username_0 on 2023-09-08. Score: 4 | Accepted answer

@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.


---

