Calculating the ascendant and midheaven

[KurtPreston]

I am trying to calculate the ascendant and midheaven astrological points. My understanding is that the ascendant is calculated by determining the ecliptic longitude of due east on the horizon at a specific place and time. The midheaven is the ecliptic longitude of the zenith.

I have tried implementing this using astronomy-engine, but my implementation seems to produce incorrect results. I am unsure whether it is a mistake in my understanding of the astrology, or if it's a bug in the code, but any help would be appreciated.

export function calculateAscendant(params: {
  latitude: number;
  longitude: number;
  date: Date;
}) {
  const {latitude, longitude, date} = params;
  const observer = new Astronomy.Observer(latitude, longitude, 0);
  const horizVector: Astronomy.Vector = new Astronomy.Vector(
    0, // north-south
    1, // west-east
    0, // up-down
    new Astronomy.AstroTime(date)
  );
  const rotation = Astronomy.Rotation_HOR_EQJ(date, observer);
  const eqjVector = Astronomy.RotateVector(rotation, horizVector);
  const ecliptic = Astronomy.Ecliptic(eqjVector);
  return ecliptic.elon;
}

To get the midheaven, I instead initialize the horizVector with this:

  const horizVector: Astronomy.Vector = new Astronomy.Vector(
    0, // north-south
    0, // west-east
    1, // up-down
    new Astronomy.AstroTime(date)
  );

Note, a similar question was asked before here.

Also, thank you for this library! I love this project.

[cosinekitty]

So far your approach looks mostly correct.

I spotted one problem with the vector you intended to point due east. It actually points west: x = north, y = west, z = up. So the first vector should be (0, -1, 0).

Your midheaven calculation looks correct. Did it give you the expected answer?

[gsvjk]

@KurtPreston hi kurt, using your code i get correct values. i just verified with jhora.

[cosinekitty]

@KurtPreston, I'm looking at the formula on the Wikipedia page you mentioned, specifically the one from Meeus 1991:

It may take me a day or two, but I will study this formula so that I can understand it and figure out why it's giving a different answer than the rotation matrix approach. It is possible the Meeus formula answers a different question. It is also possible we are all looking at this wrong somehow...

[KurtPreston]

@cosinekitty Thank you for offering to look into it, but after scratching my head over this for a few days, I realize I had the definition of ascendant slightly wrong.

The horizontal vector should not be pointing due east, but rather aimed at the point of the Eastern horizon intersecting with the ecliptic. I'm not totally certain how I should calculate that angle yet, but it would explain the discrepancy.

If I can fix the matrix implementation, I'll do a comparison of values to see how close it is to the angle returned from the Meeus formula.

Any advice on how to calculate the horizontal vector with the correct azimuth would be appreciated. Again, thanks!

[cosinekitty]

Yes, this is an interesting problem! The new definition makes sense, but it's not yet obvious how to calculate it. Let me think about it...

[cosinekitty]

So far, the best approach I can figure is to implement a search. It would scan around the eastern horizon, with the azimuth as a free variable in the range (0, +180). From this azimuth, produce a horizontal vector: (variable azimuth, altitude = 0, distance = 1). Convert that horizontal vector to an ecliptic vector.

Here you have options. If you convert to a cartesian ecliptic vector (ex, ey, ez), you can search for whatever azimuth makes ez closest to zero.

If you convert to spherical coordinates (elon, elat, distance), you can search for the azimuth that makes elat closest to zero.

Astronomy Engine provides a customizable search function Search, but unfortunately, it is geared for time as the free variable, not azimuth angles.

Overall, I think the Meeus formula is the best solution for this particular application because it is a direct solution rather than a search.

[KurtPreston]

Makes sense. In case it proves useful to anyone else, here's an implementation of the Meeus formulas for ascendant and midheaven using astronomy-engine to get the ecliptic obliquity and sidereal time:

import * as Astronomy from 'astronomy-engine';

function calculateAscendant(params: {
  latitude: number;
  longitude: number;
  date: Date;
}) {
  const {date, latitude} = params;
  const localSiderealRadians = localSiderealTimeRadians(params);
  const eclipticObliquity = degreesToRadians(
    Astronomy.e_tilt(new Astronomy.AstroTime(date)).tobl
  );
  const x =
    Math.sin(localSiderealRadians) * Math.cos(eclipticObliquity) +
    Math.tan(degreesToRadians(latitude)) * Math.sin(eclipticObliquity);
  const y = -1 * Math.cos(localSiderealRadians);
  const celestialLongitudeRadians = Math.atan(y / x);
  let ascendantDegrees = radiansToDegrees(celestialLongitudeRadians);

  // Correcting the quadrant
  if (x < 0) {
    ascendantDegrees += 180;
  } else {
    ascendantDegrees += 360;
  }
  if (ascendantDegrees < 180) {
    ascendantDegrees += 180;
  } else {
    ascendantDegrees -= 180;
  }
  return ascendantDegrees;
}

function calculateMidheaven(params: {
  latitude: number;
  longitude: number;
  date: Date;
}) {
  const {date} = params;
  const localSiderealRadians = localSiderealTimeRadians(params);
  const eclipticObliquity = degreesToRadians(
    Astronomy.e_tilt(new Astronomy.AstroTime(date)).tobl
  );
  const numerator = Math.tan(localSiderealRadians);
  const denominator = Math.cos(degreesToRadians(eclipticObliquity));
  let midheavenDegrees = radiansToDegrees(Math.atan(numerator / denominator));

  // Correcting the quadrant
  if (midheavenDegrees < 0) {
    midheavenDegrees += 360;
  }
  const localSiderealDegrees = radiansToDegrees(localSiderealRadians);
  if (midheavenDegrees > localSiderealDegrees) {
    midheavenDegrees -= 180;
  }
  if (midheavenDegrees < 0) {
    midheavenDegrees += 180;
  }
  if (midheavenDegrees < 180 && localSiderealDegrees >= 180) {
    midheavenDegrees += 180;
  }
  midheavenDegrees = (midheavenDegrees + 360) % 360;
  return midheavenDegrees;
}

export function localSiderealTimeRadians(params: {
  longitude: number;
  date: Date;
}): number {
  const {longitude, date} = params;
  const grenichSiderealTime = Astronomy.SiderealTime(date);
  const localSiderealTime = grenichSiderealTime + degreesToHours(longitude);
  let localSiderealDegrees = hoursToDegrees(localSiderealTime);
  localSiderealDegrees = (localSiderealDegrees + 360) % 360;
  const localSiderealRadians = degreesToRadians(localSiderealDegrees);
  return localSiderealRadians;
}

function hoursToDegrees(hours: number): number {
  return hours * 15;
}

function degreesToHours(degrees: number): number {
  return degrees / 15;
}

function degreesToRadians(degrees: number): number {
  return degrees / 180 * Math.PI
}

function radiansToDegrees(radians: number): number {
  return radians * 180 / Math.PI;
}

[ebraminio]

Thanks to your code @KurtPreston but as e_tilt wasn't exposed in the Kotlin port my code turned it like this,

// Extracted from not exposed calculations of astronomy library
private fun eclipticObliquity(time: Time): Double {
    val t = time.tt / 36525.0
    val asec =
        ((((-0.0000000434 * t
                - 0.000000576) * t
                + 0.00200340) * t
                - 0.0001831) * t
                - 46.836769) * t + 84381.406
    val mobl = asec / 3600
    val tobl = mobl + (time.nutationEps() / 3600)
    return Math.toRadians(tobl)
}

fun calculateAscendant(latitude: Double, longitude: Double, time: Time): Double {
    val localSiderealRadians = localSiderealTimeRadians(longitude, time)
    val eclipticObliquity = eclipticObliquity(time)
    val x = sin(localSiderealRadians) * cos(eclipticObliquity) +
            tan(Math.toRadians(latitude)) * sin(eclipticObliquity)
    val y = -1 * cos(localSiderealRadians)
    val celestialLongitudeRadians = atan(y / x)
    var ascendantDegrees = Math.toDegrees(celestialLongitudeRadians)

    // Correcting the quadrant
    ascendantDegrees += if (x < 0) 180.0 else 360.0
    if (ascendantDegrees < 180.0) ascendantDegrees += 180.0 else ascendantDegrees -= 180.0
    return ascendantDegrees
}

fun calculateMidheaven(longitude: Double, time: Time): Double {
    val localSiderealRadians = localSiderealTimeRadians(longitude, time);
    val eclipticObliquity = eclipticObliquity(time)
    val numerator = tan(localSiderealRadians)
    val denominator = cos(eclipticObliquity)
    var midheavenDegrees = Math.toDegrees(atan(numerator / denominator))

    // Correcting the quadrant
    if (midheavenDegrees < 0) midheavenDegrees += 360
    val localSiderealDegrees = Math.toDegrees(localSiderealRadians)
    if (midheavenDegrees > localSiderealDegrees) midheavenDegrees -= 180
    if (midheavenDegrees < 0) midheavenDegrees += 180
    if (midheavenDegrees < 180 && localSiderealDegrees >= 180) midheavenDegrees += 180
    midheavenDegrees = (midheavenDegrees + 360) % 360
    return midheavenDegrees
}

private fun hoursToDegrees(hours: Double) = hours * 15
private fun degreesToHours(degrees: Double) = degrees / 15

private fun localSiderealTimeRadians(longitude: Double, time: Time): Double {
    val greenwichSiderealTime = siderealTime(time)
    val localSiderealTime = greenwichSiderealTime + degreesToHours(longitude)
    var localSiderealDegrees = hoursToDegrees(localSiderealTime)
    localSiderealDegrees = (localSiderealDegrees + 360) % 360
    val localSiderealRadians = Math.toRadians(localSiderealDegrees)
    return localSiderealRadians
}

Maybe exposing earthTilt in the Kotlin port can be nice.

Btw, thanks @cosinekitty for providing this library and maintaining this community. I haven't forgotten your help on providing the Kotlin port which over time has helped the free and open source app I'm maintaining to become unique by it's own. Thanks again

[ebraminio]

Turned out there is a small issue in midheaven calculation and that is degreesToRadians is being applied to eclipticObliquity twice, I've fixed that in my calculation.

[gsvjk]

@KurtPreston good job kurt. I ended up implementing the same - asc, midheaven and lunar nodes. I was able to verify asc and lunar nodes but not mid heaven. Any source where i can verify the same ?

[gsvjk]

@KurtPreston verified the midheaven calculation here. midheaven is off by 2 degrees.

        latitude: 18.966666,
        longitude: 72.833333
        2024-03-30T08:03:00.000Z // 2024 30th mar, 13:33 pm +5:30
        
        expected value : 23.54
        returned value: 21.78 // off by 2 degrees
        

[KurtPreston]

@gsvjk If I change the eclipticObliquity to use the J2000 value of 23.4392911 instead of retrieving it from astronomy-engine using Astronomy.e_tilt, then the formula outputs 23.54.

My guess is that astro-seek is using a hard-coded value for the obliquity instead of looking up it's value at the specified date.

[ebraminio]

I think midheaven calculation mismatch was due to the need of a small fix here,

  const eclipticObliquity = degreesToRadians(
    Astronomy.e_tilt(new Astronomy.AstroTime(date)).tobl
  );
  const numerator = Math.tan(localSiderealRadians);
  const denominator = Math.cos(degreesToRadians(eclipticObliquity));

That turns eclipticObliquity to radians twice.

Thanks again @KurtPreston for sharing your code as it was a solid basis for my implementation in #391