Skip to content

Commit

Permalink
[api]: Improve TPPA tests
Browse files Browse the repository at this point in the history
  • Loading branch information
@tiagohm
tiagohm committed Oct 20, 2024
1 parent a1502fd commit 05db889
Show file tree
Hide file tree
Showing 3 changed files with 216 additions and 46 deletions.
1 change: 1 addition & 0 deletions nebulosa-alignment/build.gradle.kts
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@ dependencies {
api(project(":nebulosa-time"))
api(project(":nebulosa-platesolver"))
implementation(project(":nebulosa-log"))
testImplementation(project(":nebulosa-nova"))
testImplementation(project(":nebulosa-test"))
}

Expand Down
35 changes: 35 additions & 0 deletions ...lignment/src/main/kotlin/nebulosa/alignment/polar/point/three/ThreePointPolarAlignment.kt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,9 @@ import nebulosa.platesolver.PlateSolverException
import nebulosa.time.UTC
import nebulosa.util.Resettable
import java.nio.file.Path
import kotlin.math.cos
import kotlin.math.sin
import kotlin.math.tan

/**
* Three Point Polar Alignment almost anywhere in the sky.
Expand Down Expand Up @@ -107,5 +110,37 @@ data class ThreePointPolarAlignment(
companion object {

const val DEFAULT_RADIUS: Angle = 5 * DEG2RAD

/**
* Returns the RA/DEC polar alignment error in radians, given
* the [hourAngle] and [declination] of current target,
* the [latitude] of current observation site and
* the [azimuthError], [altitudeError] alignment errors.
*/
fun computePAE(
hourAngle: Angle, declination: Angle, latitude: Angle,
azimuthError: Angle, altitudeError: Angle
): DoubleArray {
// Source: https://sourceforge.net/p/sky-simulator/code/ci/default/tree/sky_annotation.pas
// Polar error calculation based on two celestial reference points and the error of the telescope mount at these point(s).
// Based on formulas from Ralph Pass documented at https://rppass.com/align.pdf.
// They are based on the book “Telescope Control’ by Trueblood and Genet, p.111
// Ralph added sin(latitude) term in the equation for the error in RA.

// For one reference image the difference in RA and DEC caused by the misalignment of the polar axis, formula (3):
// delta_ra:= de * TAN(dec)*SIN(h) + da * (sin(lat)- COS(lat)*(TAN(dec1)*COS(h_1))
// delta_dec:=de * COS(h) + da * COS(lat)*SIN(h))}
// raJnow0:=raJnow0-dRa;
// decJnow0:=decJnow0+dDec;

val cosHA = cos(hourAngle)
val sinHA = sin(hourAngle)
val tanDEC = tan(declination)
val cosLat = cos(latitude)
val sinLat = sin(latitude)
val dRA = altitudeError * tanDEC * sinHA + azimuthError * (sinLat - cosLat * tanDEC * cosHA)
val dDEC = -altitudeError * cosHA + azimuthError * cosLat * sinHA
return doubleArrayOf(dRA, dDEC)
}
}
}
226 changes: 180 additions & 46 deletions nebulosa-alignment/src/test/kotlin/ThreePointPolarAlignmentTest.kt
View file
Original file line number Diff line number Diff line change
@@ -1,17 +1,19 @@
import io.kotest.matchers.doubles.plusOrMinus
import io.kotest.matchers.doubles.shouldBeExactly
import io.kotest.matchers.shouldBe
import nebulosa.alignment.polar.point.three.PolarErrorDetermination
import nebulosa.alignment.polar.point.three.PolarErrorDetermination.Companion.stenographicProjection
import nebulosa.alignment.polar.point.three.Position
import nebulosa.math.arcsec
import nebulosa.math.deg
import nebulosa.math.hours
import nebulosa.math.toArcsec
import nebulosa.alignment.polar.point.three.ThreePointPolarAlignment
import nebulosa.erfa.SphericalCoordinate
import nebulosa.math.*
import nebulosa.nova.position.Geoid
import nebulosa.nova.position.ICRF
import nebulosa.platesolver.PlateSolution
import nebulosa.time.TimeYMDHMS
import nebulosa.time.UTC
import org.junit.jupiter.api.Disabled
import org.junit.jupiter.api.Test
import kotlin.random.Random

class ThreePointPolarAlignmentTest {

Expand All @@ -20,8 +22,8 @@ class ThreePointPolarAlignmentTest {
val a = Position("04:14:08".hours, "-05 26 10".deg, LNG, SLAT, UTC(TimeYMDHMS(2024, 2, 10, 22, 59, 13.3739)))
val b = Position(a.vector, LNG, SLAT)

a.topocentric.azimuth shouldBeExactly b.topocentric.azimuth
a.topocentric.altitude shouldBeExactly b.topocentric.altitude
a.topocentric.azimuth shouldBe (b.topocentric.azimuth plusOrMinus 1e-14)
a.topocentric.altitude shouldBe (b.topocentric.altitude plusOrMinus 1e-14)
}

@Test
Expand All @@ -34,65 +36,197 @@ class ThreePointPolarAlignmentTest {

@Test
fun destinationCoordinates() {
val position1 = Position("05 35 18".hours, "-05 23 26".deg, LNG, SLAT, UTC(TimeYMDHMS(2024, 2, 10, 22, 58, 42.4979)))
val position2 = Position("04 54 45".hours, "-05 24 50".deg, LNG, SLAT, UTC(TimeYMDHMS(2024, 2, 10, 22, 58, 58.1655)))
val position3 = Position("04 14 08".hours, "-05 26 10".deg, LNG, SLAT, UTC(TimeYMDHMS(2024, 2, 10, 22, 59, 13.3739)))
val now = UTC(TimeYMDHMS(2024, 2, 10, 22, 58, 42.4979))
val position1 = Position("05 35 18".hours, "-05 23 26".deg, LNG, SLAT, now)
val position2 = Position("04 54 45".hours, "-05 24 50".deg, LNG, SLAT, now)
val position3 = Position("04 14 08".hours, "-05 26 10".deg, LNG, SLAT, now)
val initialFrame = PlateSolution(true, 0.0, 1.0.arcsec, "04 14 08".hours, "-05 26 10".deg, 1280.0, 1024.0)
val pe = PolarErrorDetermination(initialFrame, position1, position2, position3, LNG, SLAT)

with(pe.destinationCoordinates(0.0, 0.0)) {
this[0] shouldBe ("04 14 08".hours plusOrMinus 1e-8)
this[1] shouldBe ("-05 26 10".deg plusOrMinus 1e-8)
with(pe.destinationCoordinates(0.0, 0.0, now)) {
this[0] shouldBe ("04 14 08".hours plusOrMinus 1e-11)
this[1] shouldBe ("-05 26 10".deg plusOrMinus 1e-11)
}
}

@Test
fun perfectlyAligned() {
val now = UTC(TimeYMDHMS(2024, 10, 20, 0, 34, 47.15))
val position1 = Position("23 24 57.1".hours, "-022 38 11.6".deg, LNG, SLAT, now)
val position2 = Position("22 24 53.6".hours, "-022 37 35.5".deg, LNG, SLAT, now)
val position3 = Position("21 24 50.6".hours, "-022 36 31.5".deg, LNG, SLAT, now)

val initialFrame = PlateSolution(true, 0.0, 1.0.arcsec, "21 24 50.6".hours, "-022 36 31.5".deg, 1280.0, 1024.0)
val pe = PolarErrorDetermination(initialFrame, position1, position2, position3, LNG, SLAT)
val (az, alt) = pe.compute()

az.toArcsec shouldBe (0.0 plusOrMinus 10.0)
alt.toArcsec shouldBe (0.0 plusOrMinus 65.0)
with(computePAE(LNG, SLAT, 0.arcmin, 0.arcmin).first) {
this[0].toArcsec shouldBe (0.0 plusOrMinus 6.0)
this[1].toArcsec shouldBe (0.0 plusOrMinus 30.0) // refraction error
}
}

@Test
fun badSouthernPolarAligned() {
val now = UTC(TimeYMDHMS(2024, 10, 20, 0, 49, 0.0))
val position1 = Position("23 24 55.4".hours, "-022°52'13.7".deg, LNG, SLAT, now)
val position2 = Position("22 24 46.7".hours, "-022°47'14.9".deg, LNG, SLAT, now)
val position3 = Position("21 24 40.5".hours, "-022°41'05.4".deg, LNG, SLAT, now)
val initialFrame = PlateSolution(true, 0.0, 1.0.arcsec, "21 24 40.5".hours, "-022°41'05.4".deg, 1280.0, 1024.0)
val pe = PolarErrorDetermination(initialFrame, position1, position2, position3, LNG, SLAT)
val (az, alt) = pe.compute()
with(computePAE(LNG, SLAT, 8.arcmin, 8.arcmin).first) {
this[0].toArcsec shouldBe (480.0 plusOrMinus 6.0)
this[1].toArcsec shouldBe (480.0 plusOrMinus 30.0)
}

with(computePAE(LNG, SLAT, (-8).arcmin, 8.arcmin).first) {
this[0].toArcsec shouldBe (-480.0 plusOrMinus 6.0)
this[1].toArcsec shouldBe (480.0 plusOrMinus 30.0)
}

with(computePAE(LNG, SLAT, (-8).arcmin, (-8).arcmin).first) {
this[0].toArcsec shouldBe (-480.0 plusOrMinus 6.0)
this[1].toArcsec shouldBe (-480.0 plusOrMinus 30.0)
}

with(computePAE(LNG, SLAT, 8.arcmin, (-8).arcmin).first) {
this[0].toArcsec shouldBe (480.0 plusOrMinus 6.0)
this[1].toArcsec shouldBe (-480.0 plusOrMinus 30.0)
}

with(computePAE(LNG, SLAT, 8.arcmin, 0.arcmin).first) {
this[0].toArcsec shouldBe (480.0 plusOrMinus 6.0)
this[1].toArcsec shouldBe (0.0 plusOrMinus 30.0)
}

with(computePAE(LNG, SLAT, (-8).arcmin, 0.arcmin).first) {
this[0].toArcsec shouldBe (-480.0 plusOrMinus 6.0)
this[1].toArcsec shouldBe (0.0 plusOrMinus 30.0)
}

az.toArcsec shouldBe (900.0 plusOrMinus 10.0)
alt.toArcsec shouldBe (900.0 plusOrMinus 60.0)
with(computePAE(LNG, SLAT, 0.arcmin, (-8).arcmin).first) {
this[0].toArcsec shouldBe (0.0 plusOrMinus 6.0)
this[1].toArcsec shouldBe (-480.0 plusOrMinus 30.0)
}

with(computePAE(LNG, SLAT, 0.arcmin, (-8).arcmin).first) {
this[0].toArcsec shouldBe (0.0 plusOrMinus 6.0)
this[1].toArcsec shouldBe (-480.0 plusOrMinus 30.0)
}
}

@Test
fun badNorthernPolarAligned() {
val now = UTC(TimeYMDHMS(2024, 10, 20, 0, 53, 0.0))
val position1 = Position("23 24 09.0".hours, "-022 51 56.0".deg, LNG, NLAT, now)
val position2 = Position("22 24 00.5".hours, "-022 46 54.1".deg, LNG, NLAT, now)
val position3 = Position("21 23 54.4".hours, "-022 40 42.4".deg, LNG, NLAT, now)
val initialFrame = PlateSolution(true, 0.0, 1.0.arcsec, "21 23 54.4".hours, "-022 40 42.4".deg, 1280.0, 1024.0)
val pe = PolarErrorDetermination(initialFrame, position1, position2, position3, LNG, NLAT)
val (az, alt) = pe.compute()

az.toArcsec shouldBe (900.0 plusOrMinus 25.0)
alt.toArcsec shouldBe (900.0 plusOrMinus 90.0)
with(computePAE(LNG, NLAT, 8.arcmin, 8.arcmin).first) {
this[0].toArcsec shouldBe (480.0 plusOrMinus 6.0)
this[1].toArcsec shouldBe (480.0 plusOrMinus 30.0)
}

with(computePAE(LNG, NLAT, (-8).arcmin, 8.arcmin).first) {
this[0].toArcsec shouldBe (-480.0 plusOrMinus 6.0)
this[1].toArcsec shouldBe (480.0 plusOrMinus 30.0)
}

with(computePAE(LNG, NLAT, (-8).arcmin, (-8).arcmin).first) {
this[0].toArcsec shouldBe (-480.0 plusOrMinus 6.0)
this[1].toArcsec shouldBe (-480.0 plusOrMinus 30.0)
}

with(computePAE(LNG, NLAT, 8.arcmin, (-8).arcmin).first) {
this[0].toArcsec shouldBe (480.0 plusOrMinus 6.0)
this[1].toArcsec shouldBe (-480.0 plusOrMinus 30.0)
}

with(computePAE(LNG, NLAT, 8.arcmin, 0.arcmin).first) {
this[0].toArcsec shouldBe (480.0 plusOrMinus 6.0)
this[1].toArcsec shouldBe (0.0 plusOrMinus 30.0)
}

with(computePAE(LNG, NLAT, (-8).arcmin, 0.arcmin).first) {
this[0].toArcsec shouldBe (-480.0 plusOrMinus 6.0)
this[1].toArcsec shouldBe (0.0 plusOrMinus 30.0)
}

with(computePAE(LNG, NLAT, 0.arcmin, (-8).arcmin).first) {
this[0].toArcsec shouldBe (0.0 plusOrMinus 6.0)
this[1].toArcsec shouldBe (-480.0 plusOrMinus 30.0)
}

with(computePAE(LNG, NLAT, 0.arcmin, (-8).arcmin).first) {
this[0].toArcsec shouldBe (0.0 plusOrMinus 6.0)
this[1].toArcsec shouldBe (-480.0 plusOrMinus 30.0)
}
}

@Test
fun tooBadPolarAlignment() {
with(computePAE(LNG, SLAT, 1.0.deg, 1.0.deg).first) {
this[0].toArcsec shouldBe (3600.0 plusOrMinus 30.0)
this[1].toArcsec shouldBe (3600.0 plusOrMinus 90.0)
}
}

@Test
@Disabled
fun updateWithoutSeeing() {
val (pae, ped) = computePAE(LNG, NLAT, 0.arcmin, (-8).arcmin)

with(pae) {
this[0].toArcsec shouldBe (0.0 plusOrMinus 6.0)
this[1].toArcsec shouldBe (-480.0 plusOrMinus 30.0)
}

val rightAscension = ped.currentReferenceFrame.rightAscension
val declination = ped.currentReferenceFrame.declination
val referenceFrame = PlateSolution(true, 0.0, 1.0.arcsec, rightAscension, declination, 1280.0, 1024.0)

repeat(500) {
with(ped.update(PAE_TIME, pae[0], pae[1], referenceFrame)) {
this[0].toArcsec shouldBe (0.0 plusOrMinus 6.0)
this[1].toArcsec shouldBe (-480.0 plusOrMinus 30.0)
}
}
}

@Test
@Disabled
fun updateWithSeeing() {
val (pae, ped) = computePAE(LNG, NLAT, 0.arcmin, (-8).arcmin)

with(pae) {
this[0].toArcsec shouldBe (0.0 plusOrMinus 6.0)
this[1].toArcsec shouldBe (-480.0 plusOrMinus 30.0)
}

val currentReferenceFrame = ped.currentReferenceFrame
val (initialAzimuthError, initialAltitudeError) = pae

repeat(500) {
val seeing = Random.nextDouble(-1.0, 1.0).arcsec
val rightAscension = currentReferenceFrame.rightAscension + seeing
val declination = currentReferenceFrame.declination + seeing
val referenceFrame = PlateSolution(true, 0.0, 1.0.arcsec, rightAscension, declination, 1280.0, 1024.0)

with(ped.update(PAE_TIME, initialAzimuthError, initialAltitudeError, referenceFrame)) {
this[0].toArcsec shouldBe (0.0 plusOrMinus 65.0)
this[1].toArcsec shouldBe (-480.0 plusOrMinus 65.0)
}
}
}

companion object {

@JvmStatic private val SLAT = "-022 30 00".deg
@JvmStatic private val NLAT = "+022 30 00".deg
@JvmStatic private val LNG = "-045 30 00".deg
private val SLAT = (-22.0).deg
private val NLAT = 22.0.deg
private val LNG = (-45.5).deg

private val PAE_TIME = UTC(TimeYMDHMS(2024, 10, 20, 12, 0, 0.0))

private fun computePAE(longitude: Angle, latitude: Angle, azimuthError: Angle, altitudeError: Angle): Pair<DoubleArray, PolarErrorDetermination> {
val location = Geoid.IERS2010.lonLat(longitude, latitude)
val rightAscension = location.lstAt(PAE_TIME) // zenith
val declination = latitude

fun computePosition(raOffset: Angle): Pair<Position, SphericalCoordinate> {
val (dRA, dDEC) = ThreePointPolarAlignment.computePAE(-raOffset, declination, latitude, azimuthError, altitudeError)
val solver = ICRF.equatorial(rightAscension + raOffset + dRA, declination + dDEC, time = PAE_TIME, epoch = PAE_TIME, center = location).equatorial()
val position = Position(solver.longitude, solver.latitude, longitude, latitude, PAE_TIME)
return position to solver
}

val (position1) = computePosition((-1.0).hours)
val (position2) = computePosition(0.0.hours)
val (position3, solver3) = computePosition(1.0.hours)

val initialFrame = PlateSolution(true, 0.0, 1.0.arcsec, solver3.longitude, solver3.latitude, 1280.0, 1024.0)
val pe = PolarErrorDetermination(initialFrame, position1, position2, position3, longitude, latitude)

return pe.compute() to pe
}
}
}

0 comments on commit 05db889

Please sign in to comment.