Class: Astronoby::Precession

Inherits:
Object
  • Object
show all
Defined in:
lib/astronoby/precession.rb

Overview

Computes precession matrices using the IAU 2006 Fukushima-Williams angles with ICRS frame bias. Also provides coordinate precession using the classical method.

Class Method Summary collapse

Instance Method Summary collapse

Constructor Details

#initialize(instant:) ⇒ Precession

Returns a new instance of Precession.

Parameters:



17
18
19
 
# File 'lib/astronoby/precession.rb', line 17

def initialize(instant:)
  @instant = instant
end

Class Method Details

.for_equatorial_coordinates(coordinates:, epoch:) ⇒ Astronoby::Coordinates::Equatorial

Precesses equatorial coordinates from their current epoch to a new epoch.

Parameters:

Returns:



79
80
81
 
# File 'lib/astronoby/precession.rb', line 79

def self.for_equatorial_coordinates(coordinates:, epoch:)
  precess(coordinates, epoch)
end

.matrix_for(instant) ⇒ Matrix

Computes the combined precession-bias matrix for a given instant.

Parameters:

Returns:

  • (Matrix) —

    the 3x3 precession-bias matrix PB(t)



12
13
14
 
# File 'lib/astronoby/precession.rb', line 12

def self.matrix_for(instant)
  new(instant: instant).matrix
end

.matrix_for_epoch(epoch) ⇒ Matrix

Computes the classical precession matrix for a given epoch.

Parameters:

  • epoch (Numeric) —

    a Julian Date epoch

Returns:

  • (Matrix) —

    3x3 precession matrix



115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
 
# File 'lib/astronoby/precession.rb', line 115

def self.matrix_for_epoch(epoch)
  t = (epoch - JulianDate::DEFAULT_EPOCH) / Constants::DAYS_PER_JULIAN_CENTURY

  zeta = Angle.from_degrees(
    0.6406161 * t + 0.0000839 * t * t + 0.000005 * t * t * t
  )
  z = Angle.from_degrees(
    0.6406161 * t + 0.0003041 * t * t + 0.0000051 * t * t * t
  )
  theta = Angle.from_degrees(
    0.5567530 * t - 0.0001185 * t * t - 0.0000116 * t * t * t
  )

  cx = zeta.cos
  sx = zeta.sin
  cz = z.cos
  sz = z.sin
  ct = theta.cos
  st = theta.sin

  Matrix[
    [
      cx * ct * cz - sx * sz,
      cx * ct * sz + sx * cz,
      cx * st
    ],
    [
      -sx * ct * cz - cx * sz,
      -sx * ct * sz + cx * cz,
      -sx * st
    ],
    [
      -st * cz,
      -st * sz,
      ct
    ]
  ]
end

.precess(coordinates, epoch) ⇒ Astronoby::Coordinates::Equatorial

Returns precessed coordinates.

Parameters:

Returns:



87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
 
# File 'lib/astronoby/precession.rb', line 87

def self.precess(coordinates, epoch)
  matrix_a = matrix_for_epoch(coordinates.epoch)
  matrix_b = matrix_for_epoch(epoch).transpose

  vector = Vector[
    coordinates.right_ascension.cos * coordinates.declination.cos,
    coordinates.right_ascension.sin * coordinates.declination.cos,
    coordinates.declination.sin
  ]

  s = matrix_a * vector
  w = matrix_b * s

  Coordinates::Equatorial.new(
    right_ascension: Util::Trigonometry.adjustement_for_arctangent(
      Angle.from_radians(w[1]),
      Angle.from_radians(w[0]),
      Angle.atan(w[1] / w[0])
    ),
    declination: Angle.asin(w[2]),
    epoch: epoch
  )
end

Instance Method Details

#matrix ⇒ Matrix

Computes the combined precession-bias matrix PB(t) = R1(-εA) R3(-ψ̄) R1(φ̄) R3(γ̄)

Returns:

  • (Matrix) —

    the 3x3 precession-bias matrix



24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
 
# File 'lib/astronoby/precession.rb', line 24

def matrix
  # Fukushima-Williams precession angles including frame bias.
  # Source:
  #  IERS Conventions 2010, Section 5.6.4
  #  Wallace & Capitaine (2006), A&A 459, 981
  #  ERFA eraPfw06 / eraFw2m
  # PB(t) = R1(−εA) R3(−ψ̄) R1(φ̄) R3(γ̄)

  cache.fetch(cache_key) do
    gamma_bar = ((((
      +0.0000000260 * t +
      -0.000002788) * t +
      -0.00031238) * t +
      +0.4932044) * t +
      +10.556378) * t +
      -0.052928

    phi_bar = ((((
      -0.0000000176 * t +
      -0.000000440) * t +
      +0.00053289) * t +
      +0.0511268) * t +
      -46.811016) * t +
      +84381.412819

    psi_bar = ((((
      -0.0000000148 * t +
      -0.000026452) * t +
      -0.00018522) * t +
      +1.5584175) * t +
      +5038.481484) * t +
      -0.041775

    gamma_bar = Angle.from_degree_arcseconds(gamma_bar)
    phi_bar = Angle.from_degree_arcseconds(phi_bar)
    psi_bar = Angle.from_degree_arcseconds(psi_bar)
    eps_a = MeanObliquity.at(@instant)

    Rotation.about_x(-eps_a) * Rotation.about_z(-psi_bar) *
      Rotation.about_x(phi_bar) * Rotation.about_z(gamma_bar)
  end
end