2
$\begingroup$

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.

$\endgroup$
11
  • 1
    $\begingroup$ 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. $\endgroup$ Sep 7, 2023 at 2:36
  • 5
    $\begingroup$ Are you looking for the subsolar point? $\endgroup$
    – Mike G
    Sep 7, 2023 at 3:34
  • 4
    $\begingroup$ 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. celestialprogramming.com/snippets/geographicPosition.html $\endgroup$ Sep 7, 2023 at 3:58
  • 1
    $\begingroup$ "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, which is essentially aligned to the J2000.0 equatorial frame. $\endgroup$
    – PM 2Ring
    Sep 7, 2023 at 8:07
  • 2
    $\begingroup$ Solution removed from post and put into an answer. $\endgroup$ Sep 8, 2023 at 15:54

1 Answer 1

4
$\begingroup$

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

$\endgroup$

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .