predictor/cmd/compare_prediction/main.go

package main

import (
	"context"
	"encoding/json"
	"fmt"
	"io"
	"math"
	"net/http"
	"os"
	"time"

	"predictor-refactored/internal/dataset"
	"predictor-refactored/internal/downloader"
	"predictor-refactored/internal/prediction"

	"go.uber.org/zap"
)

// Downloads a few forecast steps and runs a prediction, then compares
// against the public Tawhiri API.
func main() {
	log, _ := zap.NewDevelopment()

	cfg := &downloader.Config{
		DataDir:  os.TempDir(),
		Parallel: 4,
	}
	dl := downloader.NewDownloader(cfg, log)

	ctx := context.Background()

	// Find latest run
	run, err := dl.FindLatestRun(ctx)
	if err != nil {
		fmt.Fprintf(os.Stderr, "FindLatestRun: %v\n", err)
		os.Exit(1)
	}
	fmt.Printf("Using run: %s\n", run.Format("2006010215"))

	// Create dataset and download first 10 steps (0-27 hours, enough for a prediction)
	dsPath := fmt.Sprintf("/tmp/pred_test_%s.bin", run.Format("2006010215"))
	defer os.Remove(dsPath)

	ds, err := dataset.Create(dsPath)
	if err != nil {
		fmt.Fprintf(os.Stderr, "Create: %v\n", err)
		os.Exit(1)
	}

	date := run.Format("20060102")
	runHour := run.Hour()
	stepsToDownload := []int{0, 3, 6, 9, 12, 15, 18, 21, 24, 27}

	fmt.Printf("Downloading %d steps...\n", len(stepsToDownload))
	for _, step := range stepsToDownload {
		hourIdx := dataset.HourIndex(step)
		fmt.Printf("  step %d (hour idx %d)...\n", step, hourIdx)

		urlA := dataset.GribURL(date, runHour, step)
		if err := dl.DownloadAndBlit(ctx, ds, urlA, hourIdx, dataset.LevelSetA); err != nil {
			fmt.Fprintf(os.Stderr, "  pgrb2 step %d: %v\n", step, err)
			os.Exit(1)
		}

		urlB := dataset.GribURLB(date, runHour, step)
		if err := dl.DownloadAndBlit(ctx, ds, urlB, hourIdx, dataset.LevelSetB); err != nil {
			fmt.Fprintf(os.Stderr, "  pgrb2b step %d: %v\n", step, err)
			os.Exit(1)
		}
	}
	ds.Flush()
	fmt.Println("Download complete")

	// Set dataset time
	ds.DSTime = run

	// Run our prediction
	launchLat := 52.2135
	launchLon := 0.0964 // already in [0, 360)
	launchAlt := 0.0
	ascentRate := 5.0
	burstAlt := 30000.0
	descentRate := 5.0

	// Launch 3 hours into the forecast
	launchTime := run.Add(3 * time.Hour)
	launchTimestamp := float64(launchTime.Unix())
	dsEpoch := float64(run.Unix())

	warnings := &prediction.Warnings{}
	stages := prediction.StandardProfile(ascentRate, burstAlt, descentRate, ds, dsEpoch, warnings, nil)
	results := prediction.RunPrediction(launchTimestamp, launchLat, launchLon, launchAlt, stages)

	fmt.Printf("\n=== Our prediction ===\n")
	for i, sr := range results {
		stage := "ascent"
		if i == 1 {
			stage = "descent"
		}
		first := sr.Points[0]
		last := sr.Points[len(sr.Points)-1]
		fmt.Printf("  %s: %d points, start=(%.4f, %.4f, %.0fm) end=(%.4f, %.4f, %.0fm)\n",
			stage, len(sr.Points),
			first.Lat, first.Lng, first.Alt,
			last.Lat, last.Lng, last.Alt)
	}

	// Get landing point
	var ourLandLat, ourLandLon float64
	if len(results) >= 2 {
		last := results[1].Points[len(results[1].Points)-1]
		ourLandLat = last.Lat
		ourLandLon = last.Lng
		if ourLandLon > 180 {
			ourLandLon -= 360
		}
	}
	fmt.Printf("  Landing: lat=%.4f, lon=%.4f\n", ourLandLat, ourLandLon)

	// Compare against public Tawhiri API
	fmt.Printf("\n=== Tawhiri API comparison ===\n")
	tawhiriLandLat, tawhiriLandLon, err := queryTawhiri(launchLat, launchLon, launchAlt, launchTime, ascentRate, burstAlt, descentRate)
	if err != nil {
		fmt.Printf("  Tawhiri API error: %v\n", err)
		fmt.Println("  (Cannot compare — Tawhiri may use a different dataset)")
		ds.Close()
		return
	}
	fmt.Printf("  Tawhiri landing: lat=%.4f, lon=%.4f\n", tawhiriLandLat, tawhiriLandLon)

	dist := haversine(ourLandLat, ourLandLon, tawhiriLandLat, tawhiriLandLon)
	fmt.Printf("  Distance between landing points: %.2f km\n", dist/1000)

	if dist < 1000 {
		fmt.Println("  CLOSE MATCH (< 1 km)")
	} else if dist < 50000 {
		fmt.Printf("  MODERATE DIFFERENCE (%.1f km) — likely different datasets\n", dist/1000)
	} else {
		fmt.Printf("  LARGE DIFFERENCE (%.1f km) — possible bug\n", dist/1000)
	}

	ds.Close()
}

func queryTawhiri(lat, lon, alt float64, launchTime time.Time, ascentRate, burstAlt, descentRate float64) (landLat, landLon float64, err error) {
	url := fmt.Sprintf(
		"https://api.v2.sondehub.org/tawhiri?launch_latitude=%.4f&launch_longitude=%.4f&launch_altitude=%.0f&launch_datetime=%s&ascent_rate=%.1f&burst_altitude=%.0f&descent_rate=%.1f",
		lat, lon, alt, launchTime.Format(time.RFC3339), ascentRate, burstAlt, descentRate)

	resp, err := http.Get(url)
	if err != nil {
		return 0, 0, err
	}
	defer resp.Body.Close()

	body, _ := io.ReadAll(resp.Body)
	if resp.StatusCode != 200 {
		return 0, 0, fmt.Errorf("HTTP %d: %s", resp.StatusCode, string(body))
	}

	var result struct {
		Prediction []struct {
			Stage      string `json:"stage"`
			Trajectory []struct {
				Latitude  float64 `json:"latitude"`
				Longitude float64 `json:"longitude"`
				Altitude  float64 `json:"altitude"`
			} `json:"trajectory"`
		} `json:"prediction"`
	}

	if err := json.Unmarshal(body, &result); err != nil {
		return 0, 0, err
	}

	for _, stage := range result.Prediction {
		if stage.Stage == "descent" && len(stage.Trajectory) > 0 {
			last := stage.Trajectory[len(stage.Trajectory)-1]
			return last.Latitude, last.Longitude, nil
		}
	}

	return 0, 0, fmt.Errorf("no descent stage found")
}

func haversine(lat1, lon1, lat2, lon2 float64) float64 {
	const R = 6371000.0
	phi1 := lat1 * math.Pi / 180
	phi2 := lat2 * math.Pi / 180
	dphi := (lat2 - lat1) * math.Pi / 180
	dlam := (lon2 - lon1) * math.Pi / 180
	a := math.Sin(dphi/2)*math.Sin(dphi/2) + math.Cos(phi1)*math.Cos(phi2)*math.Sin(dlam/2)*math.Sin(dlam/2)
	return R * 2 * math.Atan2(math.Sqrt(a), math.Sqrt(1-a))
}