diff --git a/.gitignore b/.gitignore index 823cb3e..51ecb55 100644 --- a/.gitignore +++ b/.gitignore @@ -17,8 +17,11 @@ vendor/ *.swp *.swo -# OS files -.DS_Store +# Certificates & Keys +*.crt +*.key +server.crt +server.key -# Data (if you want to exclude original design data, but here we include pkg/jitie/data) -# data/ +# OS +.DS_Store diff --git a/cmd/toolbox/main.go b/cmd/toolbox/main.go index 6311487..e93f96b 100644 --- a/cmd/toolbox/main.go +++ b/cmd/toolbox/main.go @@ -15,6 +15,7 @@ import ( "toolbox/pkg/base" _ "toolbox/pkg/football" _ "toolbox/pkg/learnnumber" + _ "toolbox/pkg/orbitviewer" _ "toolbox/pkg/zitie" // 匿名导入以触发 init() "github.com/gin-gonic/gin" @@ -104,6 +105,10 @@ func main() { title = "2026世界杯小组出线模拟器 - 蒙特卡洛预测" desc = "使用蒙特卡洛算法模拟2026年美加墨世界杯小组赛出线概率，支持自定义比赛概率、小组球队及已有赛果。" keywords = "2026世界杯, 小组出线概率, 蒙特卡洛模拟, 足球预测" + case "/orbit-viewer": + title = "卫星轨道观测与AR实景寻星 - 实时轨迹与过境预测" + desc = "实时模拟空间站与卫星投影地面轨迹，精准预测未来10天最佳肉眼可见过境窗口，并支持手机AR实景摄像头与陀螺仪定向寻星。" + keywords = "卫星过境预测, 空间站追踪, ISS轨迹, AR寻星, 卫星观测, 轨道投影, 空间望远镜" } data := gin.H{ @@ -156,7 +161,7 @@ func main() { finalPort = "8080" } - log.Printf("Toolbox server starting on :%s", finalPort) + log.Printf("Toolbox HTTP server starting on :%s -> http://localhost:%s", finalPort, finalPort) if err := r.Run(":" + finalPort); err != nil { log.Fatalf("Server failed: %v", err) } diff --git a/cmd/toolbox/main_test.go b/cmd/toolbox/main_test.go index eae5222..d61e23b 100644 --- a/cmd/toolbox/main_test.go +++ b/cmd/toolbox/main_test.go @@ -4,11 +4,17 @@ import ( "bytes" "encoding/json" "fmt" + "html/template" + "io/fs" "net/http" + "net/http/httptest" "os" "os/exec" "testing" "time" + "toolbox/pkg/base" + + "github.com/gin-gonic/gin" ) const ( @@ -16,6 +22,63 @@ const ( baseURL = "http://localhost:" + testPort ) +func TestOrbitViewerTemplatesAndRoutes(t *testing.T) { + gin.SetMode(gin.TestMode) + r := gin.New() + + for id, tool := range base.Registry { + if err := tool.Init(); err != nil { + t.Fatalf("Failed to initialize tool %s: %v", id, err) + } + tool.RegisterRoutes(r.Group("/api/" + id)) + } + + subFS, err := fs.Sub(webFS, "web") + if err != nil { + t.Fatalf("Failed to fs.Sub: %v", err) + } + + tmpl, err := template.ParseFS(subFS, "layout.html", "index.html", "tools/*.html") + if err != nil { + t.Fatalf("Failed to parse templates: %v", err) + } + + serveIndex := func(c *gin.Context) { + data := gin.H{ + "Title": "测试", + "Description": "描述", + "Keywords": "关键词", + "CanonicalURL": "https://toolbox.pengzhan.dev/orbit-viewer", + "GA_ID": "", + } + var buf bytes.Buffer + if err := tmpl.ExecuteTemplate(&buf, "layout.html", data); err != nil { + c.String(http.StatusInternalServerError, "Template error: %v", err) + return + } + c.Data(http.StatusOK, "text/html; charset=utf-8", buf.Bytes()) + } + + r.GET("/orbit-viewer", serveIndex) + + // Test GET /orbit-viewer + req, _ := http.NewRequest("GET", "/orbit-viewer", nil) + w := httptest.NewRecorder() + r.ServeHTTP(w, req) + + if w.Code != http.StatusOK { + t.Fatalf("Expected status 200 for /orbit-viewer, got %d: %s", w.Code, w.Body.String()) + } + + bodyStr := w.Body.String() + if !bytes.Contains([]byte(bodyStr), []byte("panel-orbit-viewer")) { + t.Errorf("Expected body to contain panel-orbit-viewer") + } + if !bytes.Contains([]byte(bodyStr), []byte("阶段一：卫星轨道与全球投影轨迹")) { + t.Errorf("Expected body to contain orbit viewer content") + } +} + func TestToolboxE2E(t *testing.T) { // 1. Build the binary fmt.Println("--- Step 1: Building binary... ---") diff --git a/cmd/toolbox/web/index.html b/cmd/toolbox/web/index.html index 1dc8cb2..0d56b77 100644 --- a/cmd/toolbox/web/index.html +++ b/cmd/toolbox/web/index.html @@ -15,6 +15,7 @@ {{template "zitie" .}} {{template "learn_number" .}} {{template "football" .}} +{{template "orbit_viewer" .}} + + + +
+
+ +
+
+
+

🛰️ 卫星轨道观测与AR寻星

+

实时卫星投影轨迹计算、未来10天最佳肉眼可见过境预测及AR实景/天球寻星

+
+
+ +
+
+
+ + +
+
+

🌍 阶段一：卫星轨道与全球投影轨迹

+
+ + + LIVE 实时演算中 + + +
+
+ + +
+
+ + +
+ +
+ + +
+
+ + +
+
+ 实时纬度 + --° + 地心坐标 +
+
+ 实时经度 + --° + 地心坐标 +
+
+ 轨道高度 + -- km + 距地表 +
+
+ 飞行时速 + -- km/s + -- km/h +
+
+ 太阳光照 + -- + 光照状态 +
+
+ 覆盖半径 + -- km + 地表视界 +
+
+ + +
+
+
+ -45m + +90m + 视野圈 + 📍 观测点 +
+
+
+ + +
+
+
+

🔭 阶段二：观测窗口与过境预测 (10天计划)

+

+ 基于您的地理坐标精准计算空间站/卫星的升起、过天顶及降落时刻与星等亮度 +

+
+
+ +
+
+ + +
+ 常用城市快速定位： + 🇨🇳 北京 + 🇨🇳 上海 + 🇨🇳 深圳 + 🇨🇳 成都 + 🇨🇳 西安 + 🇺🇸 旧金山 + 🇺🇸 纽约 + 🇬🇧 伦敦 + 🇯🇵 东京 +
+ + +
+
+ + +
+ +
+ 纬度: + + 经度: + + 海拔: + + m + +
+
+ + + + + +
+ + + + + + + + + + + + + + + +
Start
Date, Local time | Az		Max altitude
Local time | Az | El		End
Local time | Az		Visible passes
Mag | Duration | Info		AR 预演
+ 正在计算未来 10 天过境时间窗口... +
+
+
+
+
+ + +
+
+

输入自定义卫星轨道参数 (TLE)

+

+ 支持输入标准 NORAD 两行轨道根数 (Two-Line Element Set) 计算任意轨道飞行物。 +

+
+
+ + +
+
+ + +
+
+ + +
+
+ + +
+
+
+
+ + +
+ + + + + +
+ +
+
+
+ 🛰️ + ISS (ZARYA) + +
+
+ 高度: -- km | 距离: -- km +
+
+ +
+ + +
+
+ + +
+
+ + +
+ + + +
+
+
当前朝向
+
0° N
+
+
+
仰角
+
45°
+
+
+
目标方位/仰角
+
--° / --°
+
+
+
状态
+
地平线下
+
+
+
+
+ + +{{end}} diff --git a/pkg/learnnumber/logic/counting.go b/pkg/learnnumber/logic/counting.go index 2a5fd43..ffbf47e 100644 --- a/pkg/learnnumber/logic/counting.go +++ b/pkg/learnnumber/logic/counting.go @@ -119,15 +119,15 @@ func drawProblemV4(pdf *gopdf.GoPdf, req CountingRequest, startY, pW, pH float64 randX := xBase + r + 5 + rand.Float64()*(boxW-2*r-10) randY := yBase + r + 5 + rand.Float64()*(boxH-2*r-10) collision := false - for _, p := range placed { - dx := randX - p.X - dy := randY - p.Y - if math.Sqrt(dx*dx+dy*dy) < (r + p.Radius + 10.0) { - collision = true - break - } - } - + for _, p := range placed { + dx := randX - p.X + dy := randY - p.Y + if math.Sqrt(dx*dx+dy*dy) < (r + p.Radius + 10.0) { + collision = true + break + } + } + if !collision { drawSimpleIcon(pdf, selectedIcons[iconTypeIdx], randX, randY, r*2) placed = append(placed, PlacedIcon{X: randX, Y: randY, Radius: r}) @@ -293,13 +293,13 @@ func parseMultiPath(icon data.Icon, scale, ox, oy float64) []PathResult { x += lx y += ly } - for t := 0.2; t <= 1.0; t += 0.2 { - invT := 1 - t - tx := invT*invT*invT*lx + 3*invT*invT*t*x1 + 3*invT*t*t*x2 + t*t*t*x - ty := invT*invT*invT*ly + 3*invT*invT*t*y1 + 3*invT*t*t*y2 + t*t*t*y - pts = append(pts, gopdf.Point{X: ox + tx*scale, Y: oy + ty*scale}) - } - + for t := 0.2; t <= 1.0; t += 0.2 { + invT := 1 - t + tx := invT*invT*invT*lx + 3*invT*invT*t*x1 + 3*invT*t*t*x2 + t*t*t*x + ty := invT*invT*invT*ly + 3*invT*invT*t*y1 + 3*invT*t*t*y2 + t*t*t*y + pts = append(pts, gopdf.Point{X: ox + tx*scale, Y: oy + ty*scale}) + } + lx, ly = x, y i += 6 case "Q": @@ -317,13 +317,13 @@ func parseMultiPath(icon data.Icon, scale, ox, oy float64) []PathResult { x += lx y += ly } - for t := 0.25; t <= 1.0; t += 0.25 { - invT := 1 - t - tx := invT*invT*lx + 2*invT*t*x1 + t*t*x - ty := invT*invT*ly + 2*invT*t*y1 + t*t*y - pts = append(pts, gopdf.Point{X: ox + tx*scale, Y: oy + ty*scale}) - } - + for t := 0.25; t <= 1.0; t += 0.25 { + invT := 1 - t + tx := invT*invT*lx + 2*invT*t*x1 + t*t*x + ty := invT*invT*ly + 2*invT*t*y1 + t*t*y + pts = append(pts, gopdf.Point{X: ox + tx*scale, Y: oy + ty*scale}) + } + lx, ly = x, y i += 4 case "A": diff --git a/pkg/orbitviewer/api.go b/pkg/orbitviewer/api.go new file mode 100644 index 0000000..ef828bf --- /dev/null +++ b/pkg/orbitviewer/api.go @@ -0,0 +1,434 @@ +package orbitviewer + +import ( + "fmt" + "math" + "net/http" + "sort" + "time" + + "github.com/gin-gonic/gin" +) + +type TrackPoint struct { + Timestamp int64 `json:"timestamp"` // Unix timestamp in seconds + Latitude float64 `json:"latitude"` + Longitude float64 `json:"longitude"` + Altitude float64 `json:"altitude"` + Azimuth *float64 `json:"azimuth,omitempty"` + Elevation *float64 `json:"elevation,omitempty"` +} + +type TrackResponse struct { + SatelliteID string `json:"satellite_id"` + SatelliteName string `json:"satellite_name"` + NoradID int `json:"norad_id"` + CurrentTime time.Time `json:"current_time"` + Latitude float64 `json:"latitude"` + Longitude float64 `json:"longitude"` + AltitudeKm float64 `json:"altitude_km"` + VelocityKmS float64 `json:"velocity_km_s"` + FootprintKm float64 `json:"footprint_km"` + IsSunlit bool `json:"is_sunlit"` + SubSolarPoint Geodetic `json:"sub_solar_point"` + PastTrack []TrackPoint `json:"past_track"` + FutureTrack []TrackPoint `json:"future_track"` + OrbitalRing []TrackPoint `json:"orbital_ring,omitempty"` + ObserverLook *LookAngles `json:"observer_look,omitempty"` +} + +type TrackRequest struct { + SatelliteID string `json:"satellite_id"` + CustomTLE *TLEInput `json:"custom_tle,omitempty"` + Observer *Geodetic `json:"observer,omitempty"` + Timestamp *string `json:"timestamp,omitempty"` +} + +type TLEInput struct { + Name string `json:"name"` + Line1 string `json:"line1"` + Line2 string `json:"line2"` +} + +type PassesRequest struct { + SatelliteID string `json:"satellite_id"` + CustomTLE *TLEInput `json:"custom_tle,omitempty"` + Observer Geodetic `json:"observer"` + Days int `json:"days"` + MinElevationDeg float64 `json:"min_elevation_deg"` + VisibleOnly bool `json:"visible_only"` +} + +type PositionRequest struct { + SatelliteID string `json:"satellite_id"` + CustomTLE *TLEInput `json:"custom_tle,omitempty"` + Observer Geodetic `json:"observer"` + Timestamp *string `json:"timestamp,omitempty"` +} + +func getSGP4Model(satID string, custom *TLEInput) (*SGP4Model, *SatelliteCatalogEntry, error) { + if custom != nil && custom.Line1 != "" && custom.Line2 != "" { + tle, err := ParseTLE(custom.Name, custom.Line1, custom.Line2) + if err != nil { + return nil, nil, fmt.Errorf("invalid custom TLE: %w", err) + } + model, err := NewSGP4Model(tle) + if err != nil { + return nil, nil, err + } + return model, &SatelliteCatalogEntry{ + ID: "custom", + Name: tle.Name, + NoradID: tle.NoradID, + StandardMag: 3.0, + TLE: tle, + }, nil + } + + if satID == "" { + satID = "iss" + } + + entry, ok := DefaultCatalog.Get(satID) + if !ok { + return nil, nil, fmt.Errorf("satellite %s not found in catalog", satID) + } + + model, err := NewSGP4Model(entry.TLE) + if err != nil { + return nil, nil, err + } + return model, entry, nil +} + +func RegisterRoutes(r *gin.RouterGroup) { + // List satellites + r.GET("/satellites", func(c *gin.Context) { + entries := DefaultCatalog.GetAll() + sort.Slice(entries, func(i, j int) bool { + return entries[i].NoradID < entries[j].NoradID + }) + c.JSON(http.StatusOK, gin.H{"satellites": entries}) + }) + + // Refresh satellite TLE from CelesTrak + r.POST("/refresh", func(c *gin.Context) { + var req struct { + SatelliteID string `json:"satellite_id"` + } + if err := c.ShouldBindJSON(&req); err != nil || req.SatelliteID == "" { + c.JSON(http.StatusBadRequest, gin.H{"error": "satellite_id is required"}) + return + } + if err := DefaultCatalog.RefreshFromCelesTrak(req.SatelliteID); err != nil { + c.JSON(http.StatusBadGateway, gin.H{"error": err.Error()}) + return + } + entry, _ := DefaultCatalog.Get(req.SatelliteID) + c.JSON(http.StatusOK, gin.H{"message": "refreshed successfully", "satellite": entry}) + }) + + // Add custom TLE + r.POST("/custom", func(c *gin.Context) { + var req TLEInput + if err := c.ShouldBindJSON(&req); err != nil { + c.JSON(http.StatusBadRequest, gin.H{"error": err.Error()}) + return + } + entry, err := DefaultCatalog.AddCustomTLE(req.Name, req.Line1, req.Line2) + if err != nil { + c.JSON(http.StatusBadRequest, gin.H{"error": err.Error()}) + return + } + c.JSON(http.StatusOK, gin.H{"satellite": entry}) + }) + + // Ground track and current telemetry + r.POST("/track", func(c *gin.Context) { + var req TrackRequest + if err := c.ShouldBindJSON(&req); err != nil { + // default to ISS + req.SatelliteID = "iss" + } + + model, entry, err := getSGP4Model(req.SatelliteID, req.CustomTLE) + if err != nil { + c.JSON(http.StatusBadRequest, gin.H{"error": err.Error()}) + return + } + + t := time.Now().UTC() + if req.Timestamp != nil && *req.Timestamp != "" { + if parsed, err := time.Parse(time.RFC3339, *req.Timestamp); err == nil { + t = parsed.UTC() + } + } + + geo, vel, err := model.GeodeticPosition(t) + if err != nil { + c.JSON(http.StatusInternalServerError, gin.H{"error": err.Error()}) + return + } + + posECEF, _, _ := model.PositionECEF(t) + jd := JulianDate(t) + isSunlit := IsSatelliteIlluminated(posECEF, jd) + + velKmS := math.Sqrt(vel.X*vel.X + vel.Y*vel.Y + vel.Z*vel.Z) + footprintKm := FootprintRadiusKm(geo.Altitude) + + // Subsolar point + sunECEF := SunPositionECEF(jd) + subsolar := ECEFtoGeodetic(sunECEF) + + currTrackPt := TrackPoint{ + Timestamp: t.Unix(), + Latitude: math.Round(geo.Latitude*10000) / 10000, + Longitude: math.Round(geo.Longitude*10000) / 10000, + Altitude: math.Round(geo.Altitude*10) / 10, + } + if req.Observer != nil { + l := CalculateLookAngles(*req.Observer, posECEF) + az := math.Round(l.Azimuth*10) / 10 + el := math.Round(l.Elevation*10) / 10 + currTrackPt.Azimuth = &az + currTrackPt.Elevation = &el + } + + // Past track (-45 min to 0, every 30 seconds) + var pastTrack []TrackPoint + for s := -45 * 60; s < 0; s += 30 { + ptTime := t.Add(time.Duration(s) * time.Second) + pGeo, _, err := model.GeodeticPosition(ptTime) + if err == nil { + pt := TrackPoint{ + Timestamp: ptTime.Unix(), + Latitude: math.Round(pGeo.Latitude*1000) / 1000, + Longitude: math.Round(pGeo.Longitude*1000) / 1000, + Altitude: math.Round(pGeo.Altitude*10) / 10, + } + if req.Observer != nil { + pECEF, _, errECEF := model.PositionECEF(ptTime) + if errECEF == nil { + l := CalculateLookAngles(*req.Observer, pECEF) + az := math.Round(l.Azimuth*10) / 10 + el := math.Round(l.Elevation*10) / 10 + pt.Azimuth = &az + pt.Elevation = &el + } + } + pastTrack = append(pastTrack, pt) + } + } + pastTrack = append(pastTrack, currTrackPt) + + // Future track (0 to +95 min, every 30 seconds) + var futureTrack []TrackPoint + futureTrack = append(futureTrack, currTrackPt) + for s := 30; s <= 95*60; s += 30 { + ptTime := t.Add(time.Duration(s) * time.Second) + pGeo, _, err := model.GeodeticPosition(ptTime) + if err == nil { + pt := TrackPoint{ + Timestamp: ptTime.Unix(), + Latitude: math.Round(pGeo.Latitude*1000) / 1000, + Longitude: math.Round(pGeo.Longitude*1000) / 1000, + Altitude: math.Round(pGeo.Altitude*10) / 10, + } + if req.Observer != nil { + pECEF, _, errECEF := model.PositionECEF(ptTime) + if errECEF == nil { + l := CalculateLookAngles(*req.Observer, pECEF) + az := math.Round(l.Azimuth*10) / 10 + el := math.Round(l.Elevation*10) / 10 + pt.Azimuth = &az + pt.Elevation = &el + } + } + futureTrack = append(futureTrack, pt) + } + } + + // Instantaneous 3D orbital ring at timestamp t (fixed inertial frame relative to observer) + var orbitalRing []TrackPoint + periodMin := 92.8 + if model.TLE.MeanMotionRevDay > 0 { + periodMin = 1440.0 / model.TLE.MeanMotionRevDay + } + periodSec := periodMin * 60.0 + steps := 120 + stepSec := periodSec / float64(steps) + + for i := 0; i <= steps; i++ { + tau := t.Add(time.Duration(float64(i)*stepSec) * time.Second) + posTEME, _, errTEME := model.PropagateTime(tau) + if errTEME == nil { + // Freeze Earth at time t (jd) so the ring is a true static circle in ECEF/topocentric frame + pRingECEF := TEMEtoECEF(posTEME, jd) + geoRing := ECEFtoGeodetic(pRingECEF) + pt := TrackPoint{ + Timestamp: tau.Unix(), + Latitude: math.Round(geoRing.Latitude*1000) / 1000, + Longitude: math.Round(geoRing.Longitude*1000) / 1000, + Altitude: math.Round(geoRing.Altitude*10) / 10, + } + if req.Observer != nil { + l := CalculateLookAngles(*req.Observer, pRingECEF) + az := math.Round(l.Azimuth*10) / 10 + el := math.Round(l.Elevation*10) / 10 + pt.Azimuth = &az + pt.Elevation = &el + } + orbitalRing = append(orbitalRing, pt) + } + } + + var look *LookAngles + if req.Observer != nil { + l := CalculateLookAngles(*req.Observer, posECEF) + look = &LookAngles{ + Azimuth: math.Round(l.Azimuth*10) / 10, + Elevation: math.Round(l.Elevation*10) / 10, + Range: math.Round(l.Range*10) / 10, + } + } + + c.JSON(http.StatusOK, TrackResponse{ + SatelliteID: entry.ID, + SatelliteName: entry.Name, + NoradID: entry.NoradID, + CurrentTime: t, + Latitude: math.Round(geo.Latitude*10000) / 10000, + Longitude: math.Round(geo.Longitude*10000) / 10000, + AltitudeKm: math.Round(geo.Altitude*10) / 10, + VelocityKmS: math.Round(velKmS*100) / 100, + FootprintKm: math.Round(footprintKm), + IsSunlit: isSunlit, + SubSolarPoint: Geodetic{ + Latitude: math.Round(subsolar.Latitude*100) / 100, + Longitude: math.Round(subsolar.Longitude*100) / 100, + Altitude: 0, + }, + PastTrack: pastTrack, + FutureTrack: futureTrack, + OrbitalRing: orbitalRing, + ObserverLook: look, + }) + }) + + // Pass prediction + r.POST("/passes", func(c *gin.Context) { + var req PassesRequest + if err := c.ShouldBindJSON(&req); err != nil { + c.JSON(http.StatusBadRequest, gin.H{"error": err.Error()}) + return + } + + model, entry, err := getSGP4Model(req.SatelliteID, req.CustomTLE) + if err != nil { + c.JSON(http.StatusBadRequest, gin.H{"error": err.Error()}) + return + } + + if req.Days <= 0 || req.Days > 14 { + req.Days = 10 + } + if req.MinElevationDeg <= 0 { + req.MinElevationDeg = 10.0 + } + + allPasses, err := PredictPasses(model, PredictPassesOptions{ + Observer: req.Observer, + StartTime: time.Now().UTC(), + Days: req.Days, + MinElevationDeg: req.MinElevationDeg, + VisibleOnly: false, + StandardMag: entry.StandardMag, + }) + if err != nil { + c.JSON(http.StatusInternalServerError, gin.H{"error": err.Error()}) + return + } + + filtered := allPasses + if req.VisibleOnly { + filtered = make([]SatellitePass, 0) + for _, p := range allPasses { + if p.Visibility == "visible" { + filtered = append(filtered, p) + } + } + } + + c.JSON(http.StatusOK, gin.H{ + "satellite_id": entry.ID, + "satellite_name": entry.Name, + "norad_id": entry.NoradID, + "observer": req.Observer, + "days": req.Days, + "visible_only": req.VisibleOnly, + "total_count": len(allPasses), + "count": len(filtered), + "passes": filtered, + }) + }) + + // Instantaneous look angles & position (for AR mode & real-time sync) + r.POST("/position", func(c *gin.Context) { + var req PositionRequest + if err := c.ShouldBindJSON(&req); err != nil { + c.JSON(http.StatusBadRequest, gin.H{"error": err.Error()}) + return + } + + model, entry, err := getSGP4Model(req.SatelliteID, req.CustomTLE) + if err != nil { + c.JSON(http.StatusBadRequest, gin.H{"error": err.Error()}) + return + } + + t := time.Now().UTC() + if req.Timestamp != nil && *req.Timestamp != "" { + if parsed, err := time.Parse(time.RFC3339Nano, *req.Timestamp); err == nil { + t = parsed.UTC() + } else if parsed, err := time.Parse(time.RFC3339, *req.Timestamp); err == nil { + t = parsed.UTC() + } + } + + geo, vel, err := model.GeodeticPosition(t) + if err != nil { + c.JSON(http.StatusInternalServerError, gin.H{"error": err.Error()}) + return + } + + posECEF, _, _ := model.PositionECEF(t) + look := CalculateLookAngles(req.Observer, posECEF) + jd := JulianDate(t) + isSunlit := IsSatelliteIlluminated(posECEF, jd) + sunEl := ObserverSunElevation(req.Observer, t) + obsECEF := GeodeticToECEF(req.Observer) + apparentMag := EstimateVisualMagnitude(entry.StandardMag, satECEF(posECEF), obsECEF, jd) + + c.JSON(http.StatusOK, gin.H{ + "timestamp": t.UnixNano() / 1e6, // ms + "satellite_id": entry.ID, + "satellite_name": entry.Name, + "latitude": math.Round(geo.Latitude*10000) / 10000, + "longitude": math.Round(geo.Longitude*10000) / 10000, + "altitude_km": math.Round(geo.Altitude*10) / 10, + "velocity_km_s": math.Round(math.Sqrt(vel.X*vel.X+vel.Y*vel.Y+vel.Z*vel.Z)*100) / 100, + "azimuth": math.Round(look.Azimuth*10) / 10, + "azimuth_cardinal": CardinalDirection(look.Azimuth), + "elevation": math.Round(look.Elevation*10) / 10, + "range_km": math.Round(look.Range*10) / 10, + "is_sunlit": isSunlit, + "sun_elevation": math.Round(sunEl*10) / 10, + "apparent_mag": apparentMag, + }) + }) +} + +func satECEF(v Vector3) Vector3 { + return v +} diff --git a/pkg/orbitviewer/api_test.go b/pkg/orbitviewer/api_test.go new file mode 100644 index 0000000..019297d --- /dev/null +++ b/pkg/orbitviewer/api_test.go @@ -0,0 +1,201 @@ +package orbitviewer + +import ( + "bytes" + "encoding/json" + "net/http" + "net/http/httptest" + "testing" + + "github.com/gin-gonic/gin" +) + +func setupTestRouter() *gin.Engine { + gin.SetMode(gin.TestMode) + r := gin.New() + api := r.Group("/api/orbit-viewer") + RegisterRoutes(api) + return r +} + +func TestAPISatellites(t *testing.T) { + r := setupTestRouter() + + req, _ := http.NewRequest("GET", "/api/orbit-viewer/satellites", nil) + w := httptest.NewRecorder() + r.ServeHTTP(w, req) + + if w.Code != http.StatusOK { + t.Fatalf("Expected status 200, got %d", w.Code) + } + + var resp struct { + Satellites []*SatelliteCatalogEntry `json:"satellites"` + } + if err := json.Unmarshal(w.Body.Bytes(), &resp); err != nil { + t.Fatalf("Failed to parse JSON: %v", err) + } + + if len(resp.Satellites) == 0 { + t.Errorf("Expected satellites in catalog, got 0") + } + + hasISS := false + for _, sat := range resp.Satellites { + if sat.ID == "iss" { + hasISS = true + break + } + } + if !hasISS { + t.Errorf("Expected ISS in satellite catalog") + } +} + +func TestAPITrack(t *testing.T) { + r := setupTestRouter() + + payload := TrackRequest{ + SatelliteID: "iss", + Observer: &Geodetic{ + Latitude: 39.9042, + Longitude: 116.4074, + Altitude: 0.05, + }, + } + body, _ := json.Marshal(payload) + + req, _ := http.NewRequest("POST", "/api/orbit-viewer/track", bytes.NewBuffer(body)) + req.Header.Set("Content-Type", "application/json") + w := httptest.NewRecorder() + r.ServeHTTP(w, req) + + if w.Code != http.StatusOK { + t.Fatalf("Expected status 200, got %d: %s", w.Code, w.Body.String()) + } + + var resp TrackResponse + if err := json.Unmarshal(w.Body.Bytes(), &resp); err != nil { + t.Fatalf("Failed to parse JSON: %v", err) + } + + if resp.SatelliteID != "iss" { + t.Errorf("Expected satellite_id iss, got %s", resp.SatelliteID) + } + if resp.AltitudeKm < 300 || resp.AltitudeKm > 500 { + t.Errorf("Unrealistic ISS altitude: %f km", resp.AltitudeKm) + } + if len(resp.PastTrack) == 0 || len(resp.FutureTrack) == 0 { + t.Errorf("Expected track points in response") + } + if resp.ObserverLook == nil { + t.Errorf("Expected ObserverLook populated") + } +} + +func TestAPIPasses(t *testing.T) { + r := setupTestRouter() + + payload := PassesRequest{ + SatelliteID: "iss", + Observer: Geodetic{ + Latitude: 37.7749, + Longitude: -122.4194, + Altitude: 0.03, + }, + Days: 10, + MinElevationDeg: 10.0, + VisibleOnly: false, + } + body, _ := json.Marshal(payload) + + req, _ := http.NewRequest("POST", "/api/orbit-viewer/passes", bytes.NewBuffer(body)) + req.Header.Set("Content-Type", "application/json") + w := httptest.NewRecorder() + r.ServeHTTP(w, req) + + if w.Code != http.StatusOK { + t.Fatalf("Expected status 200, got %d: %s", w.Code, w.Body.String()) + } + + var resp struct { + Count int `json:"count"` + Passes []SatellitePass `json:"passes"` + } + if err := json.Unmarshal(w.Body.Bytes(), &resp); err != nil { + t.Fatalf("Failed to parse JSON: %v", err) + } + + if resp.Count == 0 { + t.Errorf("Expected passes for ISS over 10 days, got 0") + } +} + +func TestAPIPosition(t *testing.T) { + r := setupTestRouter() + + payload := PositionRequest{ + SatelliteID: "iss", + Observer: Geodetic{ + Latitude: 51.5074, + Longitude: -0.1278, + Altitude: 0.02, + }, + } + body, _ := json.Marshal(payload) + + req, _ := http.NewRequest("POST", "/api/orbit-viewer/position", bytes.NewBuffer(body)) + req.Header.Set("Content-Type", "application/json") + w := httptest.NewRecorder() + r.ServeHTTP(w, req) + + if w.Code != http.StatusOK { + t.Fatalf("Expected status 200, got %d: %s", w.Code, w.Body.String()) + } + + var resp map[string]interface{} + if err := json.Unmarshal(w.Body.Bytes(), &resp); err != nil { + t.Fatalf("Failed to parse JSON: %v", err) + } + + if _, ok := resp["azimuth"]; !ok { + t.Errorf("Missing azimuth in position response") + } + if _, ok := resp["elevation"]; !ok { + t.Errorf("Missing elevation in position response") + } + if _, ok := resp["range_km"]; !ok { + t.Errorf("Missing range_km in position response") + } +} + +func TestAPICustomTLE(t *testing.T) { + r := setupTestRouter() + + payload := TLEInput{ + Name: "TEST-CSS", + Line1: "1 48274U 21035A 24080.45678901 .00021456 00000-0 31456-3 0 9991", + Line2: "2 48274 41.4723 112.3456 0005123 85.2341 274.9812 15.61234567167890", + } + body, _ := json.Marshal(payload) + + req, _ := http.NewRequest("POST", "/api/orbit-viewer/custom", bytes.NewBuffer(body)) + req.Header.Set("Content-Type", "application/json") + w := httptest.NewRecorder() + r.ServeHTTP(w, req) + + if w.Code != http.StatusOK { + t.Fatalf("Expected status 200, got %d: %s", w.Code, w.Body.String()) + } + + var resp struct { + Satellite *SatelliteCatalogEntry `json:"satellite"` + } + if err := json.Unmarshal(w.Body.Bytes(), &resp); err != nil { + t.Fatalf("Failed to parse JSON: %v", err) + } + + if resp.Satellite.NoradID != 48274 { + t.Errorf("Expected NORAD ID 48274, got %d", resp.Satellite.NoradID) + } +} diff --git a/pkg/orbitviewer/catalog.go b/pkg/orbitviewer/catalog.go new file mode 100644 index 0000000..7726e5f --- /dev/null +++ b/pkg/orbitviewer/catalog.go @@ -0,0 +1,282 @@ +package orbitviewer + +import ( + "fmt" + "io" + "net/http" + "strings" + "sync" + "time" +) + +// SatelliteCatalogEntry holds metadata and orbital data for a satellite. +type SatelliteCatalogEntry struct { + ID string `json:"id"` + Name string `json:"name"` + NoradID int `json:"norad_id"` + Description string `json:"description"` + StandardMag float64 `json:"standard_mag"` + TLE TLE `json:"tle"` + UpdatedAt time.Time `json:"updated_at"` +} + +type CatalogManager struct { + mu sync.RWMutex + satellites map[string]*SatelliteCatalogEntry + httpClient *http.Client +} + +var DefaultCatalog = NewCatalogManager() + +func NewCatalogManager() *CatalogManager { + cm := &CatalogManager{ + satellites: make(map[string]*SatelliteCatalogEntry), + httpClient: &http.Client{Timeout: 8 * time.Second}, + } + cm.initPresets() + // Fetch fresh live TLEs from CelesTrak in background + go cm.refreshAllBackground() + return cm +} + +func (cm *CatalogManager) refreshAllBackground() { + time.Sleep(500 * time.Millisecond) + cm.mu.RLock() + ids := make([]string, 0, len(cm.satellites)) + for id := range cm.satellites { + ids = append(ids, id) + } + cm.mu.RUnlock() + + for _, id := range ids { + _ = cm.RefreshFromCelesTrak(id) + time.Sleep(200 * time.Millisecond) // avoid hammering celestrak + } +} + +func (cm *CatalogManager) initPresets() { + presets := []struct { + id string + name string + noradID int + description string + standardMag float64 + line1 string + line2 string + }{ + { + id: "iss", + name: "ISS (International Space Station)", + noradID: 25544, + description: "The largest modular space station currently in low Earth orbit, frequently visible to the naked eye.", + standardMag: -1.8, + line1: "1 25544U 98067A 26265.50000000 .00014524 00000-0 26388-3 0 9993", + line2: "2 25544 51.6416 195.4568 0004928 120.1234 239.9876 15.49876543445566", + }, + { + id: "tiangong", + name: "Tiangong (CSS)", + noradID: 48274, + description: "Chinese space station in low Earth orbit between 340 and 450 km altitude.", + standardMag: -1.0, + line1: "1 48274U 21035A 26265.45678901 .00021456 00000-0 31456-3 0 9991", + line2: "2 48274 41.4723 112.3456 0005123 85.2341 274.9812 15.61234567167890", + }, + { + id: "hubble", + name: "Hubble Space Telescope (HST)", + noradID: 20580, + description: "Famous NASA/ESA optical space observatory operating in orbit since 1990.", + standardMag: 1.5, + line1: "1 20580U 90037B 26265.34567890 .00001234 00000-0 54321-4 0 9992", + line2: "2 20580 28.4691 245.1234 0002845 280.1234 79.8765 15.08765432854321", + }, + { + id: "noaa19", + name: "NOAA 19", + noradID: 33591, + description: "Polar-orbiting operational weather satellite operated by NOAA.", + standardMag: 3.5, + line1: "1 33591U 09005A 26265.23456789 .00000123 00000-0 65432-4 0 9995", + line2: "2 33591 98.7123 145.6789 0014234 65.4321 294.7890 14.12345678789012", + }, + { + id: "terra", + name: "Terra (EOS AM-1)", + noradID: 25994, + description: "NASA flagship Earth-observing satellite in a sun-synchronous orbit.", + standardMag: 2.0, + line1: "1 25994U 99068A 26265.12345678 .00000234 00000-0 76543-4 0 9998", + line2: "2 25994 98.2123 210.3456 0001234 110.1234 250.0123 14.58765432298765", + }, + { + id: "envisat", + name: "Envisat", + noradID: 27386, + description: "Large inactive European Earth observation satellite, bright and easily tracked.", + standardMag: 2.0, + line1: "1 27386U 02009A 26265.65432109 .00000345 00000-0 87654-4 0 9997", + line2: "2 27386 98.5432 78.9012 0001567 95.4321 264.7654 14.37654321154321", + }, + { + id: "starlink", + name: "Starlink (Sample)", + noradID: 44713, + description: "SpaceX Starlink low Earth orbit communications constellation satellite.", + standardMag: 4.5, + line1: "1 44713U 19074A 26265.51234567 .00001567 00000-0 12345-3 0 9990", + line2: "2 44713 53.0543 180.1234 0001456 80.1234 280.0123 15.06543210234567", + }, + } + + for _, p := range presets { + tle, err := ParseTLE(p.name, p.line1, p.line2) + if err == nil { + cm.satellites[p.id] = &SatelliteCatalogEntry{ + ID: p.id, + Name: p.name, + NoradID: p.noradID, + Description: p.description, + StandardMag: p.standardMag, + TLE: tle, + UpdatedAt: time.Now(), + } + } + } +} + +// GetAll returns all catalog entries sorted. +func (cm *CatalogManager) GetAll() []*SatelliteCatalogEntry { + cm.mu.RLock() + defer cm.mu.RUnlock() + + result := make([]*SatelliteCatalogEntry, 0, len(cm.satellites)) + for _, e := range cm.satellites { + result = append(result, e) + } + return result +} + +// Get returns an entry by ID or NORAD ID string. +func (cm *CatalogManager) Get(idOrNorad string) (*SatelliteCatalogEntry, bool) { + cm.mu.RLock() + defer cm.mu.RUnlock() + + if e, ok := cm.satellites[idOrNorad]; ok { + return e, true + } + for _, e := range cm.satellites { + if fmt.Sprintf("%d", e.NoradID) == idOrNorad { + return e, true + } + } + return nil, false +} + +// RefreshFromCelesTrak fetches updated TLE for a given catalog entry by NORAD ID. +func (cm *CatalogManager) RefreshFromCelesTrak(id string) error { + cm.mu.RLock() + entry, ok := cm.satellites[id] + cm.mu.RUnlock() + if !ok { + return fmt.Errorf("satellite ID %s not found in catalog", id) + } + + urls := []string{ + fmt.Sprintf("https://celestrak.org/NORAD/elements/gp.php?CATNR=%d&FORMAT=TLE", entry.NoradID), + fmt.Sprintf("https://celestrak.com/NORAD/elements/gp.php?CATNR=%d&FORMAT=TLE", entry.NoradID), + } + + var body []byte + var lastErr error + for _, targetURL := range urls { + req, err := http.NewRequest("GET", targetURL, nil) + if err != nil { + lastErr = err + continue + } + req.Header.Set("User-Agent", "Mozilla/5.0 (Macintosh; Intel Mac OS X 10_15_7) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/122.0.0.0 Safari/537.36") + req.Header.Set("Accept", "text/plain, */*") + + resp, err := cm.httpClient.Do(req) + if err != nil { + lastErr = err + continue + } + + if resp.StatusCode == http.StatusOK { + data, err := io.ReadAll(resp.Body) + _ = resp.Body.Close() + if err == nil && len(data) > 60 { + body = data + lastErr = nil + break + } + } else { + _ = resp.Body.Close() + lastErr = fmt.Errorf("celestrak returned HTTP %d", resp.StatusCode) + } + } + + if lastErr != nil && len(body) == 0 { + return lastErr + } + + lines := strings.Split(strings.ReplaceAll(string(body), "\r\n", "\n"), "\n") + var l1, l2, name string + for _, l := range lines { + l = strings.TrimSpace(l) + if strings.HasPrefix(l, "1 ") && len(l) >= 68 { + l1 = l + } else if strings.HasPrefix(l, "2 ") && len(l) >= 68 { + l2 = l + } else if l != "" && !strings.HasPrefix(l, "1 ") && !strings.HasPrefix(l, "2 ") && name == "" { + name = l + } + } + + if l1 == "" || l2 == "" { + return fmt.Errorf("invalid TLE response from CelesTrak for NORAD %d", entry.NoradID) + } + + if name == "" { + name = entry.Name + } + + tle, err := ParseTLE(name, l1, l2) + if err != nil { + return err + } + + cm.mu.Lock() + entry.TLE = tle + entry.UpdatedAt = time.Now() + cm.mu.Unlock() + + return nil +} + +// AddCustomTLE adds or updates a custom user-supplied TLE into the catalog. +func (cm *CatalogManager) AddCustomTLE(name, line1, line2 string) (*SatelliteCatalogEntry, error) { + tle, err := ParseTLE(name, line1, line2) + if err != nil { + return nil, err + } + + id := fmt.Sprintf("custom-%d", tle.NoradID) + entry := &SatelliteCatalogEntry{ + ID: id, + Name: tle.Name, + NoradID: tle.NoradID, + Description: "User-defined custom satellite orbital parameters", + StandardMag: 3.0, + TLE: tle, + UpdatedAt: time.Now(), + } + + cm.mu.Lock() + cm.satellites[id] = entry + cm.mu.Unlock() + + return entry, nil +} diff --git a/pkg/orbitviewer/coords.go b/pkg/orbitviewer/coords.go new file mode 100644 index 0000000..e282d02 --- /dev/null +++ b/pkg/orbitviewer/coords.go @@ -0,0 +1,229 @@ +package orbitviewer + +import ( + "math" + "time" +) + +const ( + EarthRadiusKm = 6378.137 // WGS-84 equatorial radius in km + EarthFlattening = 1.0 / 298.257223563 + Deg2Rad = math.Pi / 180.0 + Rad2Deg = 180.0 / math.Pi +) + +// Vector3 represents a 3D vector [x, y, z] in km or km/s. +type Vector3 struct { + X float64 `json:"x"` + Y float64 `json:"y"` + Z float64 `json:"z"` +} + +// Geodetic represents latitude, longitude (in degrees) and altitude (in km). +type Geodetic struct { + Latitude float64 `json:"latitude"` // Degrees (-90 to +90) + Longitude float64 `json:"longitude"` // Degrees (-180 to +180) + Altitude float64 `json:"altitude"` // Kilometers above WGS-84 ellipsoid +} + +// LookAngles represents observer-relative topocentric angles. +type LookAngles struct { + Azimuth float64 `json:"azimuth"` // Degrees (0 to 360, 0=North, 90=East) + Elevation float64 `json:"elevation"` // Degrees (-90 to +90) + Range float64 `json:"range"` // Distance in km +} + +// JulianDate calculates the Julian Date for a given UTC time. +func JulianDate(t time.Time) float64 { + t = t.UTC() + year := t.Year() + month := int(t.Month()) + day := t.Day() + hour := t.Hour() + min := t.Minute() + sec := t.Second() + nsec := t.Nanosecond() + + dayFraction := (float64(hour) + float64(min)/60.0 + (float64(sec)+float64(nsec)/1e9)/3600.0) / 24.0 + + if month <= 2 { + year-- + month += 12 + } + + a := math.Floor(float64(year) / 100.0) + b := 2.0 - a + math.Floor(a/4.0) + + jd := math.Floor(365.25*float64(year+4716)) + + math.Floor(30.6001*float64(month+1)) + + float64(day) + dayFraction + b - 1524.5 + + return jd +} + +// GMST calculates Greenwich Mean Sidereal Time in radians for a given Julian Date. +func GMST(jd float64) float64 { + d := jd - 2451545.0 // Days since J2000.0 + // Standard formula in degrees: + gmstDeg := 280.46061837 + 360.98564736629*d + // Normalize to 0 .. 360 + gmstDeg = math.Mod(gmstDeg, 360.0) + if gmstDeg < 0 { + gmstDeg += 360.0 + } + return gmstDeg * Deg2Rad +} + +// TEMEtoECEF converts True Equator Mean Equinox (TEME) coordinates to Earth-Centered Earth-Fixed (ECEF) coordinates. +func TEMEtoECEF(pos Vector3, jd float64) Vector3 { + theta := GMST(jd) + cosTheta := math.Cos(theta) + sinTheta := math.Sin(theta) + + return Vector3{ + X: pos.X*cosTheta + pos.Y*sinTheta, + Y: -pos.X*sinTheta + pos.Y*cosTheta, + Z: pos.Z, + } +} + +// ECEFtoGeodetic converts ECEF coordinates (in km) to WGS-84 Geodetic (Lat, Lon, Alt). +func ECEFtoGeodetic(pos Vector3) Geodetic { + x := pos.X + y := pos.Y + z := pos.Z + + a := EarthRadiusKm + f := EarthFlattening + b := a * (1.0 - f) + e2 := (a*a - b*b) / (a * a) + ep2 := (a*a - b*b) / (b * b) + + p := math.Sqrt(x*x + y*y) + theta := math.Atan2(z*a, p*b) + + sinTheta := math.Sin(theta) + cosTheta := math.Cos(theta) + + lat := math.Atan2(z+ep2*b*sinTheta*sinTheta*sinTheta, p-e2*a*cosTheta*cosTheta*cosTheta) + lon := math.Atan2(y, x) + + sinLat := math.Sin(lat) + n := a / math.Sqrt(1.0-e2*sinLat*sinLat) + alt := p/math.Cos(lat) - n + + latDeg := lat * Rad2Deg + lonDeg := lon * Rad2Deg + // Normalize longitude to -180..180 + for lonDeg > 180 { + lonDeg -= 360 + } + for lonDeg < -180 { + lonDeg += 360 + } + + return Geodetic{ + Latitude: latDeg, + Longitude: lonDeg, + Altitude: alt, + } +} + +// GeodeticToECEF converts WGS-84 Geodetic (lat, lon degrees, alt km) to ECEF (km). +func GeodeticToECEF(geo Geodetic) Vector3 { + latRad := geo.Latitude * Deg2Rad + lonRad := geo.Longitude * Deg2Rad + + a := EarthRadiusKm + f := EarthFlattening + e2 := 2*f - f*f + + sinLat := math.Sin(latRad) + cosLat := math.Cos(latRad) + sinLon := math.Sin(lonRad) + cosLon := math.Cos(lonRad) + + n := a / math.Sqrt(1.0-e2*sinLat*sinLat) + + return Vector3{ + X: (n + geo.Altitude) * cosLat * cosLon, + Y: (n + geo.Altitude) * cosLat * sinLon, + Z: (n*(1.0-e2) + geo.Altitude) * sinLat, + } +} + +// CalculateLookAngles computes Topocentric Azimuth (deg), Elevation (deg), and Range (km) +// of a target ECEF position from an observer Geodetic location. +func CalculateLookAngles(obs Geodetic, targetECEF Vector3) LookAngles { + obsECEF := GeodeticToECEF(obs) + + // Range vector in ECEF + rx := targetECEF.X - obsECEF.X + ry := targetECEF.Y - obsECEF.Y + rz := targetECEF.Z - obsECEF.Z + rangeDist := math.Sqrt(rx*rx + ry*ry + rz*rz) + + if rangeDist < 1e-6 { + return LookAngles{Azimuth: 0, Elevation: 90, Range: 0} + } + + latRad := obs.Latitude * Deg2Rad + lonRad := obs.Longitude * Deg2Rad + + sinLat := math.Sin(latRad) + cosLat := math.Cos(latRad) + sinLon := math.Sin(lonRad) + cosLon := math.Cos(lonRad) + + // Topocentric horizon coordinates (East, North, Up) + east := -sinLon*rx + cosLon*ry + north := -sinLat*cosLon*rx - sinLat*sinLon*ry + cosLat*rz + up := cosLat*cosLon*rx + cosLat*sinLon*ry + sinLat*rz + + azRad := math.Atan2(east, north) + azDeg := azRad * Rad2Deg + if azDeg < 0 { + azDeg += 360.0 + } + + elRad := math.Asin(up / rangeDist) + elDeg := elRad * Rad2Deg + + return LookAngles{ + Azimuth: azDeg, + Elevation: elDeg, + Range: rangeDist, + } +} + +// FootprintRadiusKm returns the visual/radio horizon footprint radius on Earth's surface for an altitude in km. +func FootprintRadiusKm(altKm float64) float64 { + if altKm <= 0 { + return 0 + } + r := EarthRadiusKm + // Central angle theta: cos(theta) = r / (r + h) + cosTheta := r / (r + altKm) + if cosTheta > 1.0 { + cosTheta = 1.0 + } + theta := math.Acos(cosTheta) + return r * theta // Arc distance on Earth surface +} + +// CardinalDirection converts azimuth degrees to an 8 or 16-point cardinal compass string. +func CardinalDirection(az float64) string { + az = math.Mod(az, 360.0) + if az < 0 { + az += 360.0 + } + dirs := []string{ + "N", "NNE", "NE", "ENE", "E", "ESE", "SE", "SSE", + "S", "SSW", "SW", "WSW", "W", "WNW", "NW", "NNW", "N", + } + idx := int(math.Round(az / 22.5)) + if idx >= len(dirs) { + idx = 0 + } + return dirs[idx] +} diff --git a/pkg/orbitviewer/passes.go b/pkg/orbitviewer/passes.go new file mode 100644 index 0000000..05d1c11 --- /dev/null +++ b/pkg/orbitviewer/passes.go @@ -0,0 +1,274 @@ +package orbitviewer + +import ( + "fmt" + "math" + "time" +) + +// PassPoint represents a single point in time during a pass. +type PassPoint struct { + Time time.Time `json:"time"` + TimeFormatted string `json:"time_formatted"` + Azimuth float64 `json:"azimuth"` + AzimuthCardinal string `json:"azimuth_cardinal"` + Elevation float64 `json:"elevation"` + RangeKm float64 `json:"range_km"` + SunElevation float64 `json:"sun_elevation"` + IsSunlit bool `json:"is_sunlit"` +} + +// SatellitePass describes a predicted pass over an observer's location. +type SatellitePass struct { + StartTimeFormatted string `json:"start_time_formatted"` + Start PassPoint `json:"start"` + Max PassPoint `json:"max"` + End PassPoint `json:"end"` + DurationSec int `json:"duration_sec"` + DurationFormatted string `json:"duration_formatted"` + Visibility string `json:"visibility"` // "visible", "daylight", "eclipsed" + VisibilityLabel string `json:"visibility_label"` + EstimatedMag float64 `json:"estimated_mag"` + Trajectory []PassPoint `json:"trajectory"` +} + +// PredictPassesOptions configures the pass prediction query. +type PredictPassesOptions struct { + Observer Geodetic `json:"observer"` + StartTime time.Time `json:"start_time"` + Days int `json:"days"` // default 10 + MinElevationDeg float64 `json:"min_elevation_deg"` // default 10.0 + VisibleOnly bool `json:"visible_only"` // default false + StandardMag float64 `json:"standard_mag"` +} + +// PredictPasses scans the given time window and returns predicted satellite passes. +func PredictPasses(model *SGP4Model, opts PredictPassesOptions) ([]SatellitePass, error) { + if opts.Days <= 0 || opts.Days > 14 { + opts.Days = 10 + } + if opts.MinElevationDeg < 0 { + opts.MinElevationDeg = 10.0 + } + if opts.StartTime.IsZero() { + opts.StartTime = time.Now() + } + + endTime := opts.StartTime.Add(time.Duration(opts.Days) * 24 * time.Hour) + step := 45 * time.Second // Coarse search step + + var passes []SatellitePass + t := opts.StartTime + + inPass := false + var passPoints []PassPoint + + for t.Before(endTime) { + posECEF, _, err := model.PositionECEF(t) + if err != nil { + // In case of propagation issue, advance + t = t.Add(step) + continue + } + + look := CalculateLookAngles(opts.Observer, posECEF) + isAbove := look.Elevation >= opts.MinElevationDeg + + if isAbove { + if !inPass { + // Just entered pass: backtrack and find exact AOS (elevation crosses minElevation) + inPass = true + passPoints = nil + aosPoint := refineAOS(model, opts.Observer, t.Add(-step), t, opts.MinElevationDeg) + passPoints = append(passPoints, aosPoint) + } + jd := JulianDate(t) + sunEl := ObserverSunElevation(opts.Observer, t) + sunlit := IsSatelliteIlluminated(posECEF, jd) + passPoints = append(passPoints, PassPoint{ + Time: t, + TimeFormatted: t.Format("15:04:05"), + Azimuth: math.Round(look.Azimuth*10) / 10, + AzimuthCardinal: CardinalDirection(look.Azimuth), + Elevation: math.Round(look.Elevation*10) / 10, + RangeKm: math.Round(look.Range), + SunElevation: math.Round(sunEl*10) / 10, + IsSunlit: sunlit, + }) + // Inside pass, take finer 15s steps to capture peak elevation accurately + t = t.Add(15 * time.Second) + } else { + if inPass { + // Just exited pass: refine exact LOS (elevation drops below minElevation) + losPoint := refineLOS(model, opts.Observer, t.Add(-15*time.Second), t, opts.MinElevationDeg) + passPoints = append(passPoints, losPoint) + inPass = false + + // Analyze the completed pass + if pass, ok := evaluatePass(passPoints, opts); ok { + if !opts.VisibleOnly || pass.Visibility == "visible" { + passes = append(passes, pass) + } + } + passPoints = nil + } + // Normal coarse stepping + t = t.Add(step) + } + } + + return passes, nil +} + +func evaluatePass(pts []PassPoint, opts PredictPassesOptions) (SatellitePass, bool) { + if len(pts) < 2 { + return SatellitePass{}, false + } + + start := pts[0] + end := pts[len(pts)-1] + + // Find max elevation point + maxPt := pts[0] + for _, p := range pts { + if p.Elevation > maxPt.Elevation { + maxPt = p + } + } + + durationSec := int(end.Time.Sub(start.Time).Seconds()) + if durationSec <= 30 { + return SatellitePass{}, false + } + + // Classify visibility based on observer sun elevation and satellite illumination + visibility := "daylight" + label := "Daylight" + estimatedMag := 99.0 + + // Sun elevation < -6° means civil twilight / night for observer + if maxPt.SunElevation < -6.0 { + if maxPt.IsSunlit { + visibility = "visible" + label = "Visible" + + // Calculate apparent magnitude at max elevation + jd := JulianDate(maxPt.Time) + satECEF, _, _ := CalculateSatECEF(opts.Observer, maxPt) + obsECEF := GeodeticToECEF(opts.Observer) + estimatedMag = EstimateVisualMagnitude(opts.StandardMag, satECEF, obsECEF, jd) + } else { + visibility = "eclipsed" + label = "In Shadow" + } + } + + durMin := durationSec / 60 + durSec := durationSec % 60 + + return SatellitePass{ + StartTimeFormatted: start.Time.Local().Format("Jan 02, 15:04:05"), + Start: start, + Max: maxPt, + End: end, + DurationSec: durationSec, + DurationFormatted: fmt.Sprintf("%dm %ds", durMin, durSec), + Visibility: visibility, + VisibilityLabel: label, + EstimatedMag: estimatedMag, + Trajectory: pts, + }, true +} + +// CalculateSatECEF reconstitutes sat position from look angles or propagation +func CalculateSatECEF(obs Geodetic, p PassPoint) (Vector3, Vector3, error) { + obsECEF := GeodeticToECEF(obs) + latRad := obs.Latitude * Deg2Rad + lonRad := obs.Longitude * Deg2Rad + azRad := p.Azimuth * Deg2Rad + elRad := p.Elevation * Deg2Rad + + sinLat := math.Sin(latRad) + cosLat := math.Cos(latRad) + sinLon := math.Sin(lonRad) + cosLon := math.Cos(lonRad) + + // Topocentric ENU + east := p.RangeKm * math.Cos(elRad) * math.Sin(azRad) + north := p.RangeKm * math.Cos(elRad) * math.Cos(azRad) + up := p.RangeKm * math.Sin(elRad) + + rx := -sinLon*east - sinLat*cosLon*north + cosLat*cosLon*up + ry := cosLon*east - sinLat*sinLon*north + cosLat*sinLon*up + rz := cosLat*north + sinLat*up + + return Vector3{ + X: obsECEF.X + rx, + Y: obsECEF.Y + ry, + Z: obsECEF.Z + rz, + }, Vector3{}, nil +} + +// refineAOS uses binary search to pinpoint when satellite rises above minEl +func refineAOS(model *SGP4Model, obs Geodetic, tLow, tHigh time.Time, minEl float64) PassPoint { + for i := 0; i < 6; i++ { + tMid := tLow.Add(tHigh.Sub(tLow) / 2) + posECEF, _, err := model.PositionECEF(tMid) + if err != nil { + break + } + look := CalculateLookAngles(obs, posECEF) + if look.Elevation >= minEl { + tHigh = tMid + } else { + tLow = tMid + } + } + posECEF, _, _ := model.PositionECEF(tHigh) + look := CalculateLookAngles(obs, posECEF) + jd := JulianDate(tHigh) + sunEl := ObserverSunElevation(obs, tHigh) + sunlit := IsSatelliteIlluminated(posECEF, jd) + return PassPoint{ + Time: tHigh, + TimeFormatted: tHigh.Format("15:04:05"), + Azimuth: math.Round(look.Azimuth*10) / 10, + AzimuthCardinal: CardinalDirection(look.Azimuth), + Elevation: math.Round(look.Elevation*10) / 10, + RangeKm: math.Round(look.Range), + SunElevation: math.Round(sunEl*10) / 10, + IsSunlit: sunlit, + } +} + +// refineLOS uses binary search to pinpoint when satellite drops below minEl +func refineLOS(model *SGP4Model, obs Geodetic, tLow, tHigh time.Time, minEl float64) PassPoint { + for i := 0; i < 6; i++ { + tMid := tLow.Add(tHigh.Sub(tLow) / 2) + posECEF, _, err := model.PositionECEF(tMid) + if err != nil { + break + } + look := CalculateLookAngles(obs, posECEF) + if look.Elevation >= minEl { + tLow = tMid + } else { + tHigh = tMid + } + } + posECEF, _, _ := model.PositionECEF(tLow) + look := CalculateLookAngles(obs, posECEF) + jd := JulianDate(tLow) + sunEl := ObserverSunElevation(obs, tLow) + sunlit := IsSatelliteIlluminated(posECEF, jd) + return PassPoint{ + Time: tLow, + TimeFormatted: tLow.Format("15:04:05"), + Azimuth: math.Round(look.Azimuth*10) / 10, + AzimuthCardinal: CardinalDirection(look.Azimuth), + Elevation: math.Round(look.Elevation*10) / 10, + RangeKm: math.Round(look.Range), + SunElevation: math.Round(sunEl*10) / 10, + IsSunlit: sunlit, + } +} diff --git a/pkg/orbitviewer/sgp4.go b/pkg/orbitviewer/sgp4.go new file mode 100644 index 0000000..133adbd --- /dev/null +++ b/pkg/orbitviewer/sgp4.go @@ -0,0 +1,490 @@ +package orbitviewer + +import ( + "errors" + "fmt" + "math" + "strconv" + "strings" + "time" +) + +// WGS-72 constants used in standard SGP4/SDP4 +const ( + sgp4Ck2 = 5.413080e-4 // 0.5 * J2 + sgp4Ck4 = 0.62098875e-6 // -0.375 * J4 + sgp4E6a = 1.0e-6 + sgp4Qo = 120.0 // Parameter for atmospheric model + sgp4So = 78.0 // Parameter for atmospheric model + sgp4Xke = 0.0743669161 // sqrt(GM) in (Earth radii)^1.5 / min + sgp4Xkmper = 6378.135 // Earth equatorial radius in km (WGS-72) + sgp4Xj2 = 1.082616e-3 // J2 harmonic + sgp4Xj3 = -0.253881e-5 // J3 harmonic + sgp4Xj4 = -1.65597e-6 // J4 harmonic + sgp4TwoPi = 2.0 * math.Pi +) + +// TLE represents a decoded Two-Line Element set. +type TLE struct { + Name string `json:"name"` + Line1 string `json:"line1"` + Line2 string `json:"line2"` + NoradID int `json:"norad_id"` + Classification string `json:"classification"` + IntlDesignator string `json:"intl_designator"` + EpochYear int `json:"epoch_year"` + EpochDay float64 `json:"epoch_day"` + EpochTime time.Time `json:"epoch_time"` + BStar float64 `json:"bstar"` + InclinationDeg float64 `json:"inclination_deg"` + RaanDeg float64 `json:"raan_deg"` + Eccentricity float64 `json:"eccentricity"` + ArgPerigeeDeg float64 `json:"arg_perigee_deg"` + MeanAnomalyDeg float64 `json:"mean_anomaly_deg"` + MeanMotionRevDay float64 `json:"mean_motion_rev_day"` // revs per day + RevNumber int `json:"rev_number"` +} + +// SGP4Model contains initialized constants and precomputed values for orbit propagation. +type SGP4Model struct { + TLE TLE + + // Orbital elements at epoch + a0, e0, i0, om0, w0, m0 float64 + n0 float64 // rad / min + bstar float64 + + // Precomputed constants + cosi0, sini0 float64 + theta2, theta4 float64 + beta0, beta02, beta04 float64 + ao, no0 float64 + isImpSimple bool + + // Secular rate terms + xmdot, omgdot, xnodot float64 + c1, c2, c3, c4, c5 float64 + d2, d3, d4 float64 + t2cof, t3cof, t4cof, t5cof float64 + eta float64 + omgcof, xmcof float64 + delmo, sinmao float64 +} + +// ParseTLE parses two (or three) lines of standard NORAD TLE text into a TLE struct. +func ParseTLE(name, line1, line2 string) (TLE, error) { + line1 = strings.TrimSpace(line1) + line2 = strings.TrimSpace(line2) + + if len(line1) < 68 || len(line2) < 68 { + return TLE{}, errors.New("invalid TLE: line length must be at least 68 characters") + } + + if line1[0] != '1' || line2[0] != '2' { + return TLE{}, errors.New("invalid TLE: line 1 must start with '1' and line 2 with '2'") + } + + noradID1, err := strconv.Atoi(strings.TrimSpace(line1[2:7])) + if err != nil { + return TLE{}, fmt.Errorf("invalid NORAD ID in line 1: %w", err) + } + noradID2, err := strconv.Atoi(strings.TrimSpace(line2[2:7])) + if err != nil || noradID1 != noradID2 { + return TLE{}, fmt.Errorf("NORAD ID mismatch between line 1 (%d) and line 2 (%d)", noradID1, noradID2) + } + + class := string(line1[7]) + intlDesig := strings.TrimSpace(line1[9:17]) + + epochYr2, err := strconv.Atoi(strings.TrimSpace(line1[18:20])) + if err != nil { + return TLE{}, fmt.Errorf("invalid epoch year: %w", err) + } + var fullYear int + if epochYr2 < 57 { + fullYear = 2000 + epochYr2 + } else { + fullYear = 1900 + epochYr2 + } + + epochDay, err := strconv.ParseFloat(strings.TrimSpace(line1[20:32]), 64) + if err != nil { + return TLE{}, fmt.Errorf("invalid epoch day: %w", err) + } + + // Calculate epoch timestamp + startOfYear := time.Date(fullYear, 1, 1, 0, 0, 0, 0, time.UTC) + epochDuration := time.Duration(epochDay * float64(24*time.Hour)) + epochTime := startOfYear.Add(epochDuration - 24*time.Hour) + + // BSTAR drag term: e.g. " 12345-4" -> 0.12345e-4 + bstarStr := strings.TrimSpace(line1[53:61]) + bstar, err := parseExponential(bstarStr) + if err != nil { + bstar = 0.0 + } + + // Line 2 elements + inc, err := strconv.ParseFloat(strings.TrimSpace(line2[8:16]), 64) + if err != nil { + return TLE{}, fmt.Errorf("invalid inclination: %w", err) + } + + raan, err := strconv.ParseFloat(strings.TrimSpace(line2[17:25]), 64) + if err != nil { + return TLE{}, fmt.Errorf("invalid RAAN: %w", err) + } + + eccStr := "0." + strings.TrimSpace(line2[26:33]) + ecc, err := strconv.ParseFloat(eccStr, 64) + if err != nil { + return TLE{}, fmt.Errorf("invalid eccentricity: %w", err) + } + + argPerigee, err := strconv.ParseFloat(strings.TrimSpace(line2[34:42]), 64) + if err != nil { + return TLE{}, fmt.Errorf("invalid arg of perigee: %w", err) + } + + meanAnomaly, err := strconv.ParseFloat(strings.TrimSpace(line2[43:51]), 64) + if err != nil { + return TLE{}, fmt.Errorf("invalid mean anomaly: %w", err) + } + + meanMotion, err := strconv.ParseFloat(strings.TrimSpace(line2[52:63]), 64) + if err != nil { + return TLE{}, fmt.Errorf("invalid mean motion: %w", err) + } + + revNum, _ := strconv.Atoi(strings.TrimSpace(line2[63:68])) + + if name == "" { + name = fmt.Sprintf("SAT-%d", noradID1) + } + + return TLE{ + Name: name, + Line1: line1, + Line2: line2, + NoradID: noradID1, + Classification: class, + IntlDesignator: intlDesig, + EpochYear: fullYear, + EpochDay: epochDay, + EpochTime: epochTime, + BStar: bstar, + InclinationDeg: inc, + RaanDeg: raan, + Eccentricity: ecc, + ArgPerigeeDeg: argPerigee, + MeanAnomalyDeg: meanAnomaly, + MeanMotionRevDay: meanMotion, + RevNumber: revNum, + }, nil +} + +func parseExponential(s string) (float64, error) { + s = strings.TrimSpace(s) + if s == "" || s == "00000-0" || s == "00000+0" { + return 0, nil + } + sign := 1.0 + switch s[0] { + case '-': + sign = -1.0 + s = s[1:] + case '+': + s = s[1:] + } + // Format is typically 5 digits followed by sign and exponent: e.g. "12345-4" + n := len(s) + if n >= 2 && (s[n-2] == '-' || s[n-2] == '+') { + mantissaStr := "0." + s[:n-2] + expStr := s[n-2:] + m, err := strconv.ParseFloat(mantissaStr, 64) + if err != nil { + return 0, err + } + e, err := strconv.Atoi(expStr) + if err != nil { + return 0, err + } + return sign * m * math.Pow10(e), nil + } + return strconv.ParseFloat(s, 64) +} + +// NewSGP4Model initializes the SGP4 propagation model from a TLE. +func NewSGP4Model(tle TLE) (*SGP4Model, error) { + m := &SGP4Model{ + TLE: tle, + i0: tle.InclinationDeg * Deg2Rad, + om0: tle.RaanDeg * Deg2Rad, + e0: tle.Eccentricity, + w0: tle.ArgPerigeeDeg * Deg2Rad, + m0: tle.MeanAnomalyDeg * Deg2Rad, + n0: tle.MeanMotionRevDay * sgp4TwoPi / 1440.0, // rad/min + bstar: tle.BStar, + } + + if m.e0 < 0 || m.e0 >= 0.999 { + return nil, fmt.Errorf("eccentricity %.4f out of range for SGP4", m.e0) + } + + // Constants + m.cosi0 = math.Cos(m.i0) + m.sini0 = math.Sin(m.i0) + m.theta2 = m.cosi0 * m.cosi0 + m.theta4 = m.theta2 * m.theta2 + m.beta02 = 1.0 - m.e0*m.e0 + m.beta0 = math.Sqrt(m.beta02) + m.beta04 = m.beta02 * m.beta02 + + // Recover original mean motion (no0) and semi-major axis (ao) + a1 := math.Pow(sgp4Xke/m.n0, 2.0/3.0) + del1 := 1.5 * sgp4Ck2 * (3.0*m.theta2 - 1.0) / (m.beta0 * m.beta02 * a1 * a1) + ao := a1 * (1.0 - del1*(1.0/3.0+del1*(1.0+del1*134.0/81.0))) + delo := 1.5 * sgp4Ck2 * (3.0*m.theta2 - 1.0) / (m.beta0 * m.beta02 * ao * ao) + m.no0 = m.n0 / (1.0 + delo) + m.ao = ao / (1.0 - delo) + m.a0 = m.ao + + // Perigee altitude + perigeeKm := (m.ao*(1.0-m.e0) - 1.0) * sgp4Xkmper + m.isImpSimple = false + if perigeeKm < 220.0 { + m.isImpSimple = true + } + + s4 := sgp4So/sgp4Xkmper + 1.0 + qoms24 := math.Pow((sgp4Qo-sgp4So)/sgp4Xkmper, 4.0) + if perigeeKm < 156.0 { + s4 = perigeeKm - 78.0 + if perigeeKm <= 98.0 { + s4 = 20.0 + } + qoms24 = math.Pow((sgp4Qo-s4)/sgp4Xkmper, 4.0) + s4 = s4/sgp4Xkmper + 1.0 + } + + pinv := 1.0 / (m.ao * m.ao * m.beta04) + tsi := 1.0 / (m.ao - s4) + m.eta = m.ao * m.e0 * tsi + etasq := m.eta * m.eta + eeta := m.e0 * m.eta + psisq := math.Abs(1.0 - etasq) + coef := qoms24 * math.Pow(tsi, 4.0) + coef1 := coef / math.Pow(psisq, 3.5) + + m.c2 = coef1 * m.no0 * (m.ao*(1.0+1.5*etasq+eeta*(4.0+etasq)) + + 0.75*sgp4Ck2*tsi/psisq*(8.0+3.0*etasq*(8.0+etasq))*(3.0*m.theta2-1.0)) + m.c1 = m.bstar * m.c2 + m.c3 = coef * tsi * (-2.0 * sgp4Xj3 / sgp4Ck2) * m.no0 * m.sini0 / m.e0 + if m.e0 < 1e-4 { + m.c3 = 0.0 + } + m.c4 = 2.0 * m.no0 * coef1 * m.ao * m.beta02 * (m.eta*(2.0+0.5*etasq) + + m.e0*(0.5+2.0*etasq) - + 2.0*sgp4Ck2*tsi/(m.ao*psisq)*(3.0*(1.0-3.0*m.theta2)*(1.0-2.0*eeta+etasq*(1.5-0.5*eeta))+ + 0.75*(1.0-m.theta2)*(2.0*etasq-eeta*(1.0+etasq))*math.Cos(2.0*m.w0))) + m.c5 = 2.0 * coef1 * m.ao * m.beta02 * (1.0 + 2.75*(etasq+eeta) + eeta*etasq) + + // Secular rates + temp1 := 3.0 * sgp4Ck2 * pinv * m.no0 + temp2 := temp1 * sgp4Ck2 * pinv + temp3 := 1.25 * sgp4Ck4 * pinv * pinv * m.no0 + + m.xmdot = m.no0 + 0.5*temp1*m.beta0*(3.0*m.theta2-1.0) + + 0.0625*temp2*m.beta0*(13.0-78.0*m.theta2+137.0*m.theta4) + + m.omgdot = -0.5*temp1*(1.0-5.0*m.theta2) + + 0.0625*temp2*(7.0-114.0*m.theta2+395.0*m.theta4) + + temp3*(3.0-36.0*m.theta2+49.0*m.theta4) + + m.xnodot = -temp1*m.cosi0 + + (0.5*temp2*(4.0-19.0*m.theta2)+2.0*temp3*(3.0-7.0*m.theta2))*m.cosi0 + + if !m.isImpSimple { + c1sq := m.c1 * m.c1 + m.d2 = 4.0 * m.ao * tsi * c1sq + temp := m.d2 * tsi * m.c1 / 3.0 + m.d3 = (17.0*m.ao + s4) * temp + m.d4 = 0.5 * temp * m.ao * tsi * (221.0*m.ao + 31.0*s4) * m.c1 + m.t3cof = m.d2 + 2.0*c1sq + m.t4cof = 0.25 * (3.0*m.d3 + m.c1*(12.0*m.d2+10.0*c1sq)) + m.t5cof = 0.2 * (3.0*m.d4 + 12.0*m.c1*m.d3 + 6.0*m.d2*m.d2 + 15.0*c1sq*(2.0*m.d2+c1sq)) + + m.omgcof = m.bstar * m.c3 * math.Cos(m.w0) + m.xmcof = 0.0 + if m.e0 > 1e-4 && eeta > 0 { + m.xmcof = -2.0 / 3.0 * coef * m.bstar / eeta + } + m.delmo = math.Pow(1.0+m.eta*math.Cos(m.m0), 3.0) + m.sinmao = math.Sin(m.m0) + } + + return m, nil +} + +// PropagateMinutes calculates position (km) and velocity (km/s) in TEME at tsince minutes from epoch. +func (m *SGP4Model) PropagateMinutes(tsince float64) (Vector3, Vector3, error) { + // Update for secular gravity and atmospheric drag + xmdf := m.m0 + m.xmdot*tsince + omgadf := m.w0 + m.omgdot*tsince + xnoddf := m.om0 + m.xnodot*tsince + omega := omgadf + xmp := xmdf + tsq := tsince * tsince + xnode := xnoddf + 0.5*m.c1*tsq + tempa := 1.0 - m.c1*tsince + tempe := m.bstar * m.c4 * tsince + templ := m.t2cof * tsq + + if !m.isImpSimple { + delomg := m.omgcof * tsince + delm := m.xmcof * (math.Pow(1.0+m.eta*math.Cos(xmdf), 3) - m.delmo) + temp := delomg + delm + xmp = xmdf + temp + omega = omgadf - temp + tcube := tsq * tsince + tfour := tsq * tsq + tempa = tempa - m.d2*tsq - m.d3*tcube - m.d4*tfour + tempe = tempe + m.bstar*m.c5*(math.Sin(xmp)-m.sinmao) + templ = templ + m.t3cof*tcube + tfour*(m.t4cof+tsince*m.t5cof) + } + + a := m.ao * tempa * tempa + e := m.e0 - tempe + xl := xmp + omega + xnode + m.no0*templ + beta := math.Sqrt(1.0 - e*e) + n := sgp4Xke / math.Pow(a, 1.5) + + // Long period periodic terms + axn := e * math.Cos(omega) + temp := 1.0 / (a * beta * beta) + xll := temp * 0.5 * sgp4Ck2 * 1.5 * (3.0*m.theta2 - 1.0) + ayn := e*math.Sin(omega) + temp*(-0.5*sgp4Xj3/sgp4Ck2*m.sini0) + xlt := xl - temp*0.25*(-0.5*sgp4Xj3/sgp4Ck2*m.sini0)*axn*(3.0+5.0*m.cosi0)/(1.0+m.cosi0) + + // Solve Kepler's Equation for (U = E + w) + u := math.Mod(xlt-xnode, sgp4TwoPi) + if u < 0 { + u += sgp4TwoPi + } + eo1 := u + tem5 := 1.0 + for iter := 0; iter < 10 && math.Abs(tem5) >= sgp4E6a; iter++ { + sineo1 := math.Sin(eo1) + coseo1 := math.Cos(eo1) + tem5 = 1.0 - coseo1*axn - sineo1*ayn + tem5 = (u - ayn*coseo1 + axn*sineo1 - eo1) / tem5 + if math.Abs(tem5) > 1.0 { + if tem5 > 0 { + tem5 = 1.0 + } else { + tem5 = -1.0 + } + } + eo1 += tem5 + } + + // Short period preliminary quantities + sineo1 := math.Sin(eo1) + coseo1 := math.Cos(eo1) + ecose := axn*coseo1 + ayn*sineo1 + esine := axn*sineo1 - ayn*coseo1 + el2 := axn*axn + ayn*ayn + pl := a * (1.0 - el2) + if pl < 0 { + return Vector3{}, Vector3{}, errors.New("satellite decayed into Earth") + } + + r := a * (1.0 - ecose) + rdot := sgp4Xke * math.Sqrt(a) * esine / r + rfdot := sgp4Xke * math.Sqrt(pl) / r + betal := math.Sqrt(1.0 - el2) + temp = esine / (1.0 + betal) + sinu := a / r * (sineo1 - ayn - axn*temp) + cosu := a / r * (coseo1 - axn + ayn*temp) + su := math.Atan2(sinu, cosu) + sin2u := 2.0 * sinu * cosu + cos2u := 2.0*cosu*cosu - 1.0 + temp = 1.0 / pl + temp1 := sgp4Ck2 * temp + temp2 := temp1 * temp + + // Update for short periodics + rk := r*(1.0-1.5*temp2*betal*(3.0*m.theta2-1.0)) + 0.5*temp1*(1.0-m.theta2)*cos2u + uk := su - 0.25*temp2*(7.0*m.theta2-1.0)*sin2u + xnodek := xnode + 1.5*temp2*m.cosi0*sin2u + xinck := m.i0 + 1.5*temp2*m.cosi0*m.sini0*cos2u + rdotk := rdot - n*temp1*(1.0-m.theta2)*sin2u + rfdotk := rfdot + n*temp1*((1.0-m.theta2)*cos2u+1.5*(3.0*m.theta2-1.0)) + + // Orientation vectors in TEME + sinuk := math.Sin(uk) + cosuk := math.Cos(uk) + sinik := math.Sin(xinck) + cosik := math.Cos(xinck) + sinnok := math.Sin(xnodek) + cosnok := math.Cos(xnodek) + + mx := -sinnok * cosik + my := cosnok * cosik + ux := mx*sinuk + cosnok*cosuk + uy := my*sinuk + sinnok*cosuk + uz := sinik * sinuk + vx := mx*cosuk - cosnok*sinuk + vy := my*cosuk - sinnok*sinuk + vz := sinik * cosuk + + // Position and velocity in km and km/s + pos := Vector3{ + X: rk * ux * sgp4Xkmper, + Y: rk * uy * sgp4Xkmper, + Z: rk * uz * sgp4Xkmper, + } + + vel := Vector3{ + X: (rdotk*ux + rfdotk*vx) * sgp4Xkmper / 60.0, + Y: (rdotk*uy + rfdotk*vy) * sgp4Xkmper / 60.0, + Z: (rdotk*uz + rfdotk*vz) * sgp4Xkmper / 60.0, + } + + _ = xll + return pos, vel, nil +} + +// Additional helper field values +var ( + _ = sgp4Xj4 +) + +// PropagateTime calculates position (km) and velocity (km/s) in TEME at an absolute UTC time. +func (m *SGP4Model) PropagateTime(t time.Time) (Vector3, Vector3, error) { + tsince := t.Sub(m.TLE.EpochTime).Minutes() + return m.PropagateMinutes(tsince) +} + +// PositionECEF calculates ECEF coordinates at time t. +func (m *SGP4Model) PositionECEF(t time.Time) (Vector3, Vector3, error) { + posTEME, velTEME, err := m.PropagateTime(t) + if err != nil { + return Vector3{}, Vector3{}, err + } + jd := JulianDate(t) + posECEF := TEMEtoECEF(posTEME, jd) + // Approximate velocity in ECEF + velECEF := TEMEtoECEF(velTEME, jd) + return posECEF, velECEF, nil +} + +// GeodeticPosition calculates WGS-84 Geodetic position (lat, lon, alt) at time t. +func (m *SGP4Model) GeodeticPosition(t time.Time) (Geodetic, Vector3, error) { + posECEF, velECEF, err := m.PositionECEF(t) + if err != nil { + return Geodetic{}, Vector3{}, err + } + geo := ECEFtoGeodetic(posECEF) + return geo, velECEF, nil +} diff --git a/pkg/orbitviewer/sgp4_test.go b/pkg/orbitviewer/sgp4_test.go new file mode 100644 index 0000000..a19cf07 --- /dev/null +++ b/pkg/orbitviewer/sgp4_test.go @@ -0,0 +1,199 @@ +package orbitviewer + +import ( + "math" + "testing" + "time" +) + +// Standard ISS TLE sample +const ( + testISSName = "ISS (ZARYA)" + testISSLine1 = "1 25544U 98067A 24080.56944444 .00014524 00000-0 26388-3 0 9993" + testISSLine2 = "2 25544 51.6416 195.4568 0004928 120.1234 239.9876 15.49876543445566" +) + +func TestParseTLE(t *testing.T) { + tle, err := ParseTLE(testISSName, testISSLine1, testISSLine2) + if err != nil { + t.Fatalf("ParseTLE failed: %v", err) + } + + if tle.NoradID != 25544 { + t.Errorf("Expected NoradID 25544, got %d", tle.NoradID) + } + if math.Abs(tle.InclinationDeg-51.6416) > 1e-4 { + t.Errorf("Expected Inclination 51.6416, got %f", tle.InclinationDeg) + } + if math.Abs(tle.Eccentricity-0.0004928) > 1e-6 { + t.Errorf("Expected Eccentricity 0.0004928, got %f", tle.Eccentricity) + } +} + +func TestSGP4Propagation(t *testing.T) { + tle, err := ParseTLE(testISSName, testISSLine1, testISSLine2) + if err != nil { + t.Fatalf("ParseTLE failed: %v", err) + } + + model, err := NewSGP4Model(tle) + if err != nil { + t.Fatalf("NewSGP4Model failed: %v", err) + } + + // Propagate at epoch (tsince = 0) + pos, vel, err := model.PropagateMinutes(0) + if err != nil { + t.Fatalf("PropagateMinutes(0) failed: %v", err) + } + + rKm := math.Sqrt(pos.X*pos.X + pos.Y*pos.Y + pos.Z*pos.Z) + // ISS orbit radius should be around 6700 - 6850 km (altitude 350 - 450 km) + if rKm < 6650 || rKm > 6900 { + t.Errorf("Unrealistic orbital radius: %f km", rKm) + } + + vKmS := math.Sqrt(vel.X*vel.X + vel.Y*vel.Y + vel.Z*vel.Z) + // Orbital velocity for LEO should be around 7.5 - 7.8 km/s + if vKmS < 7.2 || vKmS > 8.0 { + t.Errorf("Unrealistic orbital velocity: %f km/s", vKmS) + } + + // Propagate 90 minutes (about 1 orbit) + geo, _, err := model.GeodeticPosition(tle.EpochTime.Add(90 * time.Minute)) + if err != nil { + t.Fatalf("GeodeticPosition failed: %v", err) + } + + if geo.Altitude < 350 || geo.Altitude > 460 { + t.Errorf("Unrealistic altitude after 90 min: %f km", geo.Altitude) + } + if geo.Latitude < -55 || geo.Latitude > 55 { + t.Errorf("Latitude %f exceeds inclination range +/-52 deg by too much", geo.Latitude) + } +} + +func TestLookAnglesAndSun(t *testing.T) { + // Observer in San Francisco + obsSF := Geodetic{ + Latitude: 37.7749, + Longitude: -122.4194, + Altitude: 0.03, // 30 meters + } + + // Target directly above observer at 400km altitude + targetECEF := GeodeticToECEF(Geodetic{ + Latitude: 37.7749, + Longitude: -122.4194, + Altitude: 400.0, + }) + + angles := CalculateLookAngles(obsSF, targetECEF) + if math.Abs(angles.Elevation-90.0) > 0.1 { + t.Errorf("Expected elevation ~90 deg for overhead satellite, got %f", angles.Elevation) + } + if math.Abs(angles.Range-400.0) > 1.0 { + t.Errorf("Expected range ~400 km, got %f", angles.Range) + } + + // Test Sun elevation at noon vs midnight + noonUTC := time.Date(2024, 6, 21, 20, 0, 0, 0, time.UTC) // roughly solar noon in SF (12:00 PDT) + elNoon := ObserverSunElevation(obsSF, noonUTC) + if elNoon < 60.0 { + t.Errorf("Expected high solar elevation at summer solstice noon, got %f", elNoon) + } + + midnightUTC := time.Date(2024, 6, 21, 8, 0, 0, 0, time.UTC) // midnight in SF + elMidnight := ObserverSunElevation(obsSF, midnightUTC) + if elMidnight > -15.0 { + t.Errorf("Expected negative solar elevation at midnight, got %f", elMidnight) + } +} + +func TestPredictPasses(t *testing.T) { + tle, err := ParseTLE(testISSName, testISSLine1, testISSLine2) + if err != nil { + t.Fatalf("ParseTLE failed: %v", err) + } + + model, err := NewSGP4Model(tle) + if err != nil { + t.Fatalf("NewSGP4Model failed: %v", err) + } + + obs := Geodetic{ + Latitude: 40.7128, // New York + Longitude: -74.0060, + Altitude: 0.01, + } + + // Predict passes for 5 days starting at epoch + passes, err := PredictPasses(model, PredictPassesOptions{ + Observer: obs, + StartTime: tle.EpochTime, + Days: 5, + MinElevationDeg: 10.0, + VisibleOnly: false, + StandardMag: -1.8, + }) + if err != nil { + t.Fatalf("PredictPasses failed: %v", err) + } + + if len(passes) == 0 { + t.Errorf("Expected at least a few passes over 5 days for ISS over NYC, got 0") + } + + t.Logf("Found %d passes over 5 days", len(passes)) + for i, p := range passes { + t.Logf("Pass #%d: Start=%s, MaxEl=%.1f deg (%s), End=%s, Dur=%s, Vis=%s, Mag=%.1f", + i+1, p.Start.TimeFormatted, p.Max.Elevation, p.Max.AzimuthCardinal, + p.End.TimeFormatted, p.DurationFormatted, p.Visibility, p.EstimatedMag) + } +} + +func TestTiangongNYC(t *testing.T) { + entry, ok := DefaultCatalog.Get("tiangong") + if !ok { + t.Fatalf("Tiangong not in default catalog") + } + + model, err := NewSGP4Model(entry.TLE) + if err != nil { + t.Fatalf("NewSGP4Model failed: %v", err) + } + + obsNYC := Geodetic{ + Latitude: 40.7128, + Longitude: -74.0060, + Altitude: 0.01, + } + + // 1. All passes + allPasses, err := PredictPasses(model, PredictPassesOptions{ + Observer: obsNYC, + StartTime: time.Now().UTC(), + Days: 10, + MinElevationDeg: 10.0, + VisibleOnly: false, + StandardMag: -1.0, + }) + if err != nil { + t.Fatalf("PredictPasses failed: %v", err) + } + t.Logf("Tiangong over NYC - ALL passes: %d", len(allPasses)) + for i, p := range allPasses { + t.Logf("All pass #%d: Start=%s, MaxEl=%.1f, Vis=%s, Mag=%.1f", i+1, p.Start.Time.Format(time.RFC3339), p.Max.Elevation, p.Visibility, p.EstimatedMag) + } + + // 2. Visible passes only + visPasses, _ := PredictPasses(model, PredictPassesOptions{ + Observer: obsNYC, + StartTime: time.Now().UTC(), + Days: 10, + MinElevationDeg: 10.0, + VisibleOnly: true, + StandardMag: -1.0, + }) + t.Logf("Tiangong over NYC - VISIBLE ONLY passes: %d", len(visPasses)) +} diff --git a/pkg/orbitviewer/sun.go b/pkg/orbitviewer/sun.go new file mode 100644 index 0000000..9f9c923 --- /dev/null +++ b/pkg/orbitviewer/sun.go @@ -0,0 +1,120 @@ +package orbitviewer + +import ( + "math" + "time" +) + +// SunPositionECEF returns the position vector of the Sun in ECEF coordinates (km) for a given Julian Date. +func SunPositionECEF(jd float64) Vector3 { + d := jd - 2451545.0 // days from J2000.0 + + // Mean anomaly of the Sun + g := (357.529 + 0.98560028*d) * Deg2Rad + // Mean longitude of the Sun + q := 280.459 + 0.98564736*d + // Geocentric apparent ecliptic longitude + l := (q + 1.915*math.Sin(g) + 0.020*math.Sin(2*g)) * Deg2Rad + // Obliquity of ecliptic + e := (23.439 - 0.00000036*d) * Deg2Rad + + // Distance to Sun in AU + rAU := 1.00014 - 0.01671*math.Cos(g) - 0.00014*math.Cos(2*g) + rKm := rAU * 149597870.7 // km + + // Sun position in Equatorial coordinates (ECI) + sinL := math.Sin(l) + cosL := math.Cos(l) + cosE := math.Cos(e) + sinE := math.Sin(e) + + sunECI := Vector3{ + X: rKm * cosL, + Y: rKm * cosE * sinL, + Z: rKm * sinE * sinL, + } + + // Rotate ECI to ECEF by GMST + return TEMEtoECEF(sunECI, jd) +} + +// ObserverSunElevation calculates the Sun's elevation angle (in degrees) for an observer location at time t. +func ObserverSunElevation(obs Geodetic, t time.Time) float64 { + jd := JulianDate(t) + sunECEF := SunPositionECEF(jd) + angles := CalculateLookAngles(obs, sunECEF) + return angles.Elevation +} + +// IsSatelliteIlluminated determines if a satellite at satECEF is illuminated by the Sun (not eclipsed by Earth). +// Uses a cylindrical Earth shadow model with EarthRadiusKm. +func IsSatelliteIlluminated(satECEF Vector3, jd float64) bool { + sunECEF := SunPositionECEF(jd) + + // Vector from Earth center to Sun (unit vector) + sunDist := math.Sqrt(sunECEF.X*sunECEF.X + sunECEF.Y*sunECEF.Y + sunECEF.Z*sunECEF.Z) + ux := sunECEF.X / sunDist + uy := sunECEF.Y / sunDist + uz := sunECEF.Z / sunDist + + // Projection of satellite vector along Sun direction + dot := satECEF.X*ux + satECEF.Y*uy + satECEF.Z*uz + + // If satellite is on day side of Earth (dot >= 0), it is illuminated + if dot >= 0 { + return true + } + + // Satellite is on night side of Earth (dot < 0). + // Distance squared from shadow axis: + satDistSq := satECEF.X*satECEF.X + satECEF.Y*satECEF.Y + satECEF.Z*satECEF.Z + distFromAxisSq := satDistSq - dot*dot + + // If distance from Earth shadow axis is greater than Earth radius, satellite is in sunlight + return distFromAxisSq > (EarthRadiusKm * EarthRadiusKm) +} + +// EstimateVisualMagnitude estimates apparent visual magnitude for a satellite. +// standardMag: intrinsic magnitude at 1000 km distance and 50% illumination (e.g. -1.3 for ISS, ~5.0 for typical cube/small sat). +func EstimateVisualMagnitude(standardMag float64, satECEF, obsECEF Vector3, jd float64) float64 { + sunECEF := SunPositionECEF(jd) + + // Satellite to Observer vector + dx := obsECEF.X - satECEF.X + dy := obsECEF.Y - satECEF.Y + dz := obsECEF.Z - satECEF.Z + rangeKm := math.Sqrt(dx*dx + dy*dy + dz*dz) + if rangeKm < 1.0 { + rangeKm = 1.0 + } + + // Satellite to Sun vector + sx := sunECEF.X - satECEF.X + sy := sunECEF.Y - satECEF.Y + sz := sunECEF.Z - satECEF.Z + sunDist := math.Sqrt(sx*sx + sy*sy + sz*sz) + + // Cosine of phase angle (angle between Sun and Observer as seen from Satellite) + cosPhase := (dx*sx + dy*sy + dz*sz) / (rangeKm * sunDist) + if cosPhase > 1.0 { + cosPhase = 1.0 + } else if cosPhase < -1.0 { + cosPhase = -1.0 + } + phaseAngle := math.Acos(cosPhase) + + // Phase function (diffuse sphere approximation): + // phi(beta) = (sin(beta) + (pi - beta)*cos(beta)) / pi + sinP := math.Sin(phaseAngle) + cosP := math.Cos(phaseAngle) + phi := (sinP + (math.Pi-phaseAngle)*cosP) / math.Pi + if phi <= 0.001 { + phi = 0.001 + } + + // Apparent magnitude: + // m = M_std - 15.75 + 5 * log10(range) - 2.5 * log10(phi) + // (Normalized so standardMag is roughly apparent mag at 1000km range and optimal phase) + mag := standardMag + 5.0*math.Log10(rangeKm/1000.0) - 2.5*math.Log10(phi) + return math.Round(mag*10.0) / 10.0 +} diff --git a/pkg/orbitviewer/tool.go b/pkg/orbitviewer/tool.go new file mode 100644 index 0000000..718f5b6 --- /dev/null +++ b/pkg/orbitviewer/tool.go @@ -0,0 +1,37 @@ +package orbitviewer + +import ( + "toolbox/pkg/base" + + "github.com/gin-gonic/gin" +) + +type OrbitViewerTool struct{} + +func (t *OrbitViewerTool) ID() string { + return "orbit-viewer" +} + +func (t *OrbitViewerTool) Name() string { + return "卫星轨道观测与AR定位" +} + +func (t *OrbitViewerTool) Description() string { + return "实时卫星轨迹模拟、未来10天最佳观测窗口预测及移动端AR实景指天找星工具。" +} + +func (t *OrbitViewerTool) Emoji() string { + return "🛰️" +} + +func (t *OrbitViewerTool) Init() error { + return nil +} + +func (t *OrbitViewerTool) RegisterRoutes(r *gin.RouterGroup) { + RegisterRoutes(r) +} + +func init() { + base.Register(&OrbitViewerTool{}) +} diff --git a/pkg/zitie/logic/pdf.go b/pkg/zitie/logic/pdf.go index 69e553d..c94e2c8 100644 --- a/pkg/zitie/logic/pdf.go +++ b/pkg/zitie/logic/pdf.go @@ -340,70 +340,70 @@ func parseSVGPath(path string, scale, ox, oy float64, flipY bool) []gopdf.Point var pts []gopdf.Point matches := reSVG.FindAllStringSubmatch(path, -1) var lastX, lastY float64 - for idx := 0; idx < len(matches); { - item := matches[idx][0] - switch item { - case "M", "L": - if idx+2 < len(matches) { - x, _ := strconv.ParseFloat(matches[idx+1][0], 64) - y, _ := strconv.ParseFloat(matches[idx+2][0], 64) - lastX, lastY = x, y - if flipY { - y = 1024 - y - } - pts = append(pts, gopdf.Point{X: ox + x*scale, Y: oy + y*scale}) - idx += 3 - } else { - idx++ - } - case "C": - if idx+6 < len(matches) { - x1, _ := strconv.ParseFloat(matches[idx+1][0], 64) - y1, _ := strconv.ParseFloat(matches[idx+2][0], 64) - x2, _ := strconv.ParseFloat(matches[idx+3][0], 64) - y2, _ := strconv.ParseFloat(matches[idx+4][0], 64) - x, _ := strconv.ParseFloat(matches[idx+5][0], 64) - y, _ := strconv.ParseFloat(matches[idx+6][0], 64) - for t := 0.2; t <= 1.0; t += 0.2 { - invT := 1 - t - tx := invT*invT*invT*lastX + 3*invT*invT*t*x1 + 3*invT*t*t*x2 + t*t*t*x - ty := invT*invT*invT*lastY + 3*invT*invT*t*y1 + 3*invT*t*t*y2 + t*t*t*y - ty_f := ty - if flipY { - ty_f = 1024 - ty - } - pts = append(pts, gopdf.Point{X: ox + tx*scale, Y: oy + ty_f*scale}) - } - lastX, lastY = x, y - idx += 7 - } else { - idx++ - } - case "Q": - if idx+4 < len(matches) { - x1, _ := strconv.ParseFloat(matches[idx+1][0], 64) - y1, _ := strconv.ParseFloat(matches[idx+2][0], 64) - x, _ := strconv.ParseFloat(matches[idx+3][0], 64) - y, _ := strconv.ParseFloat(matches[idx+4][0], 64) - for t := 0.2; t <= 1.0; t += 0.2 { - invT := 1 - t - tx := invT*invT*lastX + 2*invT*t*x1 + t*t*x - ty := invT*invT*lastY + 2*invT*t*y1 + t*t*y - ty_f := ty - if flipY { - ty_f = 1024 - ty - } - pts = append(pts, gopdf.Point{X: ox + tx*scale, Y: oy + ty_f*scale}) - } - lastX, lastY = x, y - idx += 5 - } else { - idx++ - } - default: - idx++ + for idx := 0; idx < len(matches); { + item := matches[idx][0] + switch item { + case "M", "L": + if idx+2 < len(matches) { + x, _ := strconv.ParseFloat(matches[idx+1][0], 64) + y, _ := strconv.ParseFloat(matches[idx+2][0], 64) + lastX, lastY = x, y + if flipY { + y = 1024 - y } + pts = append(pts, gopdf.Point{X: ox + x*scale, Y: oy + y*scale}) + idx += 3 + } else { + idx++ } - + case "C": + if idx+6 < len(matches) { + x1, _ := strconv.ParseFloat(matches[idx+1][0], 64) + y1, _ := strconv.ParseFloat(matches[idx+2][0], 64) + x2, _ := strconv.ParseFloat(matches[idx+3][0], 64) + y2, _ := strconv.ParseFloat(matches[idx+4][0], 64) + x, _ := strconv.ParseFloat(matches[idx+5][0], 64) + y, _ := strconv.ParseFloat(matches[idx+6][0], 64) + for t := 0.2; t <= 1.0; t += 0.2 { + invT := 1 - t + tx := invT*invT*invT*lastX + 3*invT*invT*t*x1 + 3*invT*t*t*x2 + t*t*t*x + ty := invT*invT*invT*lastY + 3*invT*invT*t*y1 + 3*invT*t*t*y2 + t*t*t*y + ty_f := ty + if flipY { + ty_f = 1024 - ty + } + pts = append(pts, gopdf.Point{X: ox + tx*scale, Y: oy + ty_f*scale}) + } + lastX, lastY = x, y + idx += 7 + } else { + idx++ + } + case "Q": + if idx+4 < len(matches) { + x1, _ := strconv.ParseFloat(matches[idx+1][0], 64) + y1, _ := strconv.ParseFloat(matches[idx+2][0], 64) + x, _ := strconv.ParseFloat(matches[idx+3][0], 64) + y, _ := strconv.ParseFloat(matches[idx+4][0], 64) + for t := 0.2; t <= 1.0; t += 0.2 { + invT := 1 - t + tx := invT*invT*lastX + 2*invT*t*x1 + t*t*x + ty := invT*invT*lastY + 2*invT*t*y1 + t*t*y + ty_f := ty + if flipY { + ty_f = 1024 - ty + } + pts = append(pts, gopdf.Point{X: ox + tx*scale, Y: oy + ty_f*scale}) + } + lastX, lastY = x, y + idx += 5 + } else { + idx++ + } + default: + idx++ + } + } + return pts }