267 lines
8.8 KiB
Go
267 lines
8.8 KiB
Go
package calculator
|
|
|
|
import (
|
|
"bufio"
|
|
"bytes"
|
|
_ "embed"
|
|
"fmt"
|
|
"io"
|
|
"math"
|
|
"os"
|
|
"strconv"
|
|
"strings"
|
|
"time"
|
|
|
|
"github.com/hebl/gofa"
|
|
)
|
|
|
|
const (
|
|
ARC_SECONDS_PER_RADIAN = 206264.80624709636
|
|
)
|
|
|
|
// https://www.iers.org/IERS/EN/DataProducts/EarthOrientationData/eop.html
|
|
// https://maia.usno.navy.mil/ser7/finals2000A.all
|
|
type EOPEntry struct {
|
|
Year int
|
|
Month int
|
|
Day int
|
|
MJD float64
|
|
Xp float64 // 弧度 rad
|
|
Yp float64 // 弧度 rad
|
|
Dut1 float64
|
|
}
|
|
|
|
// https://earth-info.nga.mil/#eopp
|
|
type EOPP5Line struct {
|
|
Ta float64
|
|
A, B float64
|
|
C, D, P [2]float64
|
|
E, F float64
|
|
G, H, Q [2]float64
|
|
Tb float64
|
|
I, J float64
|
|
K, L, R [4]float64
|
|
}
|
|
|
|
type EOPTable struct {
|
|
URL string
|
|
Local2000A string
|
|
EOPs map[float64]*EOPEntry
|
|
EOPP5Line EOPP5Line
|
|
UseEOPP5Line bool
|
|
}
|
|
|
|
var EOP *EOPTable
|
|
|
|
func NewEOPTable() *EOPTable {
|
|
eop := &EOPTable{
|
|
URL: "https://maia.usno.navy.mil/ser7/finals2000A.all",
|
|
Local2000A: "data/EOP/finals2000A.all",
|
|
EOPs: make(map[float64]*EOPEntry),
|
|
UseEOPP5Line: false, // 当天使用 EOPP5Line 预测结果无法通过测试
|
|
}
|
|
return eop
|
|
}
|
|
|
|
func (eop *EOPTable) Load(eopData []byte, eopp5Line []byte) error {
|
|
f, err := os.Open(eop.Local2000A)
|
|
if err != nil {
|
|
return eop.loadEmbed(eopData, eopp5Line)
|
|
}
|
|
defer f.Close()
|
|
|
|
data, err := io.ReadAll(f)
|
|
if err != nil {
|
|
return eop.loadEmbed(eopData, eopp5Line)
|
|
}
|
|
|
|
return eop.loadEmbed(data, eopp5Line)
|
|
}
|
|
|
|
func (eop *EOPTable) loadEmbed(eopData []byte, eopp5Line []byte) error {
|
|
reader := bufio.NewReader(bytes.NewReader(eopData))
|
|
for {
|
|
line, _, err := reader.ReadLine()
|
|
if err == io.EOF {
|
|
break
|
|
}
|
|
if err != nil {
|
|
return err
|
|
}
|
|
|
|
var e EOPEntry
|
|
if err := e.line2EOPEntry(line); err != nil {
|
|
break
|
|
}
|
|
|
|
eop.EOPs[e.MJD] = &e
|
|
}
|
|
|
|
eop.EOPP5Line.parse(eopp5Line)
|
|
|
|
return nil
|
|
}
|
|
|
|
func (eop EOPTable) Get(mjd float64) (*EOPEntry, bool) {
|
|
e, ok := eop.EOPs[mjd]
|
|
if !ok && eop.UseEOPP5Line {
|
|
xp, yp, dut1 := eop.EOPP5Line.Predict(mjd)
|
|
e = &EOPEntry{
|
|
MJD: mjd,
|
|
Xp: xp / ARC_SECONDS_PER_RADIAN,
|
|
Yp: yp / ARC_SECONDS_PER_RADIAN,
|
|
Dut1: dut1,
|
|
}
|
|
ok = true
|
|
}
|
|
return e, ok
|
|
}
|
|
|
|
func (eop EOPTable) GetByUTCTime(utctime time.Time) (*EOPEntry, bool) {
|
|
var djm0, mjd float64
|
|
gofa.Cal2jd(int(utctime.Year()), int(utctime.Month()), int(utctime.Day()), &djm0, &mjd)
|
|
return eop.Get(mjd)
|
|
}
|
|
|
|
// https://maia.usno.navy.mil/ser7/readme.finals2000A
|
|
// The format of the finals2000A.data, finals2000A.daily, and finals2000A.all files is:
|
|
|
|
// Col.# Format Quantity
|
|
// ------- ------ -------------------------------------------------------------
|
|
// 1-2 I2 year (to get true calendar year, add 1900 for MJD<=51543 or add 2000 for MJD>=51544)
|
|
// 3-4 I2 month number
|
|
// 5-6 I2 day of month
|
|
// 7 X [blank]
|
|
// 8-15 F8.2 fractional Modified Julian Date (MJD UTC)
|
|
// 16 X [blank]
|
|
// 17 A1 IERS (I) or Prediction (P) flag for Bull. A polar motion values
|
|
// 18 X [blank]
|
|
// 19-27 F9.6 Bull. A PM-x (sec. of arc)
|
|
// 28-36 F9.6 error in PM-x (sec. of arc)
|
|
// 37 X [blank]
|
|
// 38-46 F9.6 Bull. A PM-y (sec. of arc)
|
|
// 47-55 F9.6 error in PM-y (sec. of arc)
|
|
// 56-57 2X [blanks]
|
|
// 58 A1 IERS (I) or Prediction (P) flag for Bull. A UT1-UTC values
|
|
// 59-68 F10.7 Bull. A UT1-UTC (sec. of time)
|
|
// 69-78 F10.7 error in UT1-UTC (sec. of time)
|
|
// 79 X [blank]
|
|
// 80-86 F7.4 Bull. A LOD (msec. of time) -- NOT ALWAYS FILLED
|
|
// 87-93 F7.4 error in LOD (msec. of time) -- NOT ALWAYS FILLED
|
|
// 94-95 2X [blanks]
|
|
// 96 A1 IERS (I) or Prediction (P) flag for Bull. A nutation values
|
|
// 97 X [blank]
|
|
// 98-106 F9.3 Bull. A dX wrt IAU2000A Nutation (msec. of arc), Free Core Nutation NOT Removed
|
|
// 107-115 F9.3 error in dX (msec. of arc)
|
|
// 116 X [blank]
|
|
// 117-125 F9.3 Bull. A dY wrt IAU2000A Nutation (msec. of arc), Free Core Nutation NOT Removed
|
|
// 126-134 F9.3 error in dY (msec. of arc)
|
|
// 135-144 F10.6 Bull. B PM-x (sec. of arc)
|
|
// 145-154 F10.6 Bull. B PM-y (sec. of arc)
|
|
// 155-165 F11.7 Bull. B UT1-UTC (sec. of time)
|
|
// 166-175 F10.3 Bull. B dX wrt IAU2000A Nutation (msec. of arc)
|
|
// 176-185 F10.3 Bull. B dY wrt IAU2000A Nutation (msec. of arc)
|
|
func (e *EOPEntry) line2EOPEntry(line []byte) error {
|
|
if len(line) < 68 {
|
|
return fmt.Errorf("invalid eop line")
|
|
}
|
|
|
|
e.Year, _ = strconv.Atoi(string(line[0:2]))
|
|
e.Month, _ = strconv.Atoi(strings.TrimSpace(string(line[2:4])))
|
|
e.Day, _ = strconv.Atoi(strings.TrimSpace(string(line[4:6])))
|
|
e.MJD, _ = strconv.ParseFloat(string(line[7:15]), 64)
|
|
|
|
e.Xp, _ = strconv.ParseFloat(strings.TrimSpace(string(line[18:27])), 64)
|
|
e.Yp, _ = strconv.ParseFloat(strings.TrimSpace(string(line[37:46])), 64)
|
|
|
|
e.Xp = e.Xp / ARC_SECONDS_PER_RADIAN
|
|
e.Yp = e.Yp / ARC_SECONDS_PER_RADIAN
|
|
|
|
e.Dut1, _ = strconv.ParseFloat(strings.TrimSpace(string(line[58:68])), 64)
|
|
|
|
return nil
|
|
}
|
|
|
|
func (e *EOPP5Line) parse(eopp5Line []byte) error {
|
|
e.P[0] = 1.0
|
|
e.P[1] = 1.0
|
|
for i := 0; i < 4; i++ {
|
|
e.R[i] = 1.0
|
|
}
|
|
|
|
reader := bufio.NewReader(bytes.NewReader(eopp5Line))
|
|
var lines [][]byte
|
|
for {
|
|
line, _, err := reader.ReadLine()
|
|
if err == io.EOF {
|
|
break
|
|
}
|
|
if err != nil {
|
|
return err
|
|
}
|
|
|
|
lines = append(lines, line)
|
|
}
|
|
|
|
if len(lines) < 5 {
|
|
return fmt.Errorf("invalid eopp5 lines")
|
|
}
|
|
|
|
e.Ta, _ = strconv.ParseFloat(strings.TrimSpace(string(lines[0][0:10])), 64)
|
|
e.A, _ = strconv.ParseFloat(strings.TrimSpace(string(lines[0][10:20])), 64)
|
|
e.B, _ = strconv.ParseFloat(strings.TrimSpace(string(lines[0][20:30])), 64)
|
|
e.C[0], _ = strconv.ParseFloat(strings.TrimSpace(string(lines[0][30:40])), 64)
|
|
e.C[1], _ = strconv.ParseFloat(strings.TrimSpace(string(lines[0][40:50])), 64)
|
|
e.D[0], _ = strconv.ParseFloat(strings.TrimSpace(string(lines[0][50:60])), 64)
|
|
e.D[1], _ = strconv.ParseFloat(strings.TrimSpace(string(lines[0][60:70])), 64)
|
|
e.P[0], _ = strconv.ParseFloat(strings.TrimSpace(string(lines[0][70:76])), 64)
|
|
|
|
e.P[1], _ = strconv.ParseFloat(strings.TrimSpace(string(lines[1][0:6])), 64)
|
|
e.E, _ = strconv.ParseFloat(strings.TrimSpace(string(lines[1][6:16])), 64)
|
|
e.F, _ = strconv.ParseFloat(strings.TrimSpace(string(lines[1][16:26])), 64)
|
|
e.G[0], _ = strconv.ParseFloat(strings.TrimSpace(string(lines[1][26:36])), 64)
|
|
e.G[1], _ = strconv.ParseFloat(strings.TrimSpace(string(lines[1][36:46])), 64)
|
|
e.H[0], _ = strconv.ParseFloat(strings.TrimSpace(string(lines[1][46:56])), 64)
|
|
e.H[1], _ = strconv.ParseFloat(strings.TrimSpace(string(lines[1][56:66])), 64)
|
|
e.Q[0], _ = strconv.ParseFloat(strings.TrimSpace(string(lines[1][66:72])), 64)
|
|
e.Q[1], _ = strconv.ParseFloat(strings.TrimSpace(string(lines[1][72:78])), 64)
|
|
|
|
e.Tb, _ = strconv.ParseFloat(strings.TrimSpace(string(lines[2][0:10])), 64)
|
|
e.I, _ = strconv.ParseFloat(strings.TrimSpace(string(lines[2][10:20])), 64)
|
|
e.J, _ = strconv.ParseFloat(strings.TrimSpace(string(lines[2][20:30])), 64)
|
|
e.K[0], _ = strconv.ParseFloat(strings.TrimSpace(string(lines[2][30:40])), 64)
|
|
e.K[1], _ = strconv.ParseFloat(strings.TrimSpace(string(lines[2][40:50])), 64)
|
|
e.K[2], _ = strconv.ParseFloat(strings.TrimSpace(string(lines[2][50:60])), 64)
|
|
e.K[3], _ = strconv.ParseFloat(strings.TrimSpace(string(lines[2][60:70])), 64)
|
|
|
|
e.L[0], _ = strconv.ParseFloat(strings.TrimSpace(string(lines[3][0:10])), 64)
|
|
e.L[1], _ = strconv.ParseFloat(strings.TrimSpace(string(lines[3][10:20])), 64)
|
|
e.L[2], _ = strconv.ParseFloat(strings.TrimSpace(string(lines[3][20:30])), 64)
|
|
e.L[3], _ = strconv.ParseFloat(strings.TrimSpace(string(lines[3][30:40])), 64)
|
|
e.R[0], _ = strconv.ParseFloat(strings.TrimSpace(string(lines[3][40:49])), 64)
|
|
e.R[1], _ = strconv.ParseFloat(strings.TrimSpace(string(lines[3][49:58])), 64)
|
|
e.R[2], _ = strconv.ParseFloat(strings.TrimSpace(string(lines[3][58:67])), 64)
|
|
e.R[3], _ = strconv.ParseFloat(strings.TrimSpace(string(lines[3][67:76])), 64)
|
|
|
|
return nil
|
|
}
|
|
|
|
// https://earth-info.nga.mil/php/download.php?file=gnss-eopp
|
|
func (e *EOPP5Line) Predict(mjd float64) (xp, yp, dut1 float64) {
|
|
xp = e.A + e.B*(mjd-e.Ta) +
|
|
e.C[0]*math.Sin(2*math.Pi*(mjd-e.Ta)/e.P[0]) + e.D[0]*math.Cos(2*math.Pi*(mjd-e.Ta)/e.P[0]) +
|
|
e.C[1]*math.Sin(2*math.Pi*(mjd-e.Ta)/e.P[1]) + e.D[1]*math.Cos(2*math.Pi*(mjd-e.Ta)/e.P[1])
|
|
|
|
yp = e.E + e.F*(mjd-e.Ta) +
|
|
e.G[0]*math.Sin(2*math.Pi*(mjd-e.Ta)/e.Q[0]) + e.H[0]*math.Cos(2*math.Pi*(mjd-e.Ta)/e.Q[0]) +
|
|
e.G[1]*math.Sin(2*math.Pi*(mjd-e.Ta)/e.Q[1]) + e.H[1]*math.Cos(2*math.Pi*(mjd-e.Ta)/e.Q[1])
|
|
|
|
dut1 = e.I + e.J*(mjd-e.Tb) +
|
|
e.K[0]*math.Sin(2*math.Pi*(mjd-e.Tb)/e.R[0]) + e.L[0]*math.Cos(2*math.Pi*(mjd-e.Tb)/e.R[0]) +
|
|
e.K[1]*math.Sin(2*math.Pi*(mjd-e.Tb)/e.R[1]) + e.L[1]*math.Cos(2*math.Pi*(mjd-e.Tb)/e.R[1]) +
|
|
e.K[2]*math.Sin(2*math.Pi*(mjd-e.Tb)/e.R[2]) + e.L[2]*math.Cos(2*math.Pi*(mjd-e.Tb)/e.R[2]) +
|
|
e.K[3]*math.Sin(2*math.Pi*(mjd-e.Tb)/e.R[3]) + e.L[3]*math.Cos(2*math.Pi*(mjd-e.Tb)/e.R[3])
|
|
|
|
return
|
|
}
|