Satellite.hs

-- Reaktor orbital challenge solution
-- Will Yager
import Data.Map (Map, keys, (!), fromList)
import Data.Set (Set, empty, notMember, insert)
import Data.List (find)
import Data.List.Split (splitOn)
data Cartesian = Cartesian Double Double Double deriving (Eq, Show)
data Polar = Polar {lat :: Double, lon :: Double, alt :: Double} deriving (Show)
data Line = Line {a :: Cartesian, b :: Cartesian}
data ID = Start | ID Int | End deriving (Eq, Ord, Show)
data Tree = Tree ID [Tree] deriving (Show)
data Path = Path [ID] deriving (Show)

-- Finds a shortest path (number of hops) to id, if it exists.
shortestPathTo :: ID -> Tree -> Maybe Path
shortestPathTo id tree = find (endsWith id) (pathsFrom tree)
    where
    endsWith id (Path (sat:_)) | sat == id = True
    endsWith id _                          = False

-- Does a breadth-first traversal of all paths
pathsFrom :: Tree -> [Path]
pathsFrom tree = map snd . concat . takeWhile (not . null) $ allPaths 
    where allPaths = iterate expand [([tree],Path [])]

-- Calculate the next layer of paths in the path tree
expand :: [([Tree], Path)] -> [([Tree], Path)]
expand paths = do
    (trees, Path steps) <- paths
    Tree new trees' <- trees
    return (trees', Path (new : steps))

-- Generate a map from satellite to its visible neighbors
visibleFrom :: Map ID Cartesian -> Map ID [ID]
visibleFrom satellites = fmap visible satellites
    where
    visible origin = filter (not . earthIntersects . Line origin . toSatellite) sats
    sats = keys satellites
    toSatellite sat = satellites ! sat

-- Generate a tree of all paths
treeFrom :: Map ID [ID] -> ID -> Tree
treeFrom = treeFrom' empty

-- Helper function for treeFrom. Prevents cycles in the graph.
treeFrom' :: Set ID -> Map ID [ID] -> ID -> Tree
treeFrom' visited visible start = Tree start branches
    where
    visited' = insert start visited 
    new      = filter (`notMember` visited') (visible ! start)
    branches = map (treeFrom' visited' visible) new

p2c :: Polar -> Cartesian
p2c (Polar lat lon r) = Cartesian x y z
    where θ = lat * pi / 180
          ϕ = lon * pi / 180
          x = r * cos θ * cos ϕ
          y = r * cos θ * sin ϕ
          z = r * sin θ

(Cartesian a b c) • (Cartesian x y z) = a*x + b*y + c*z
(Cartesian a b c) `sub` (Cartesian x y z) = Cartesian (a-x) (b-y) (c-z)
(Cartesian a b c) `plus` (Cartesian x y z) = Cartesian (a+x) (b+y) (c+z)
scale v (Cartesian a b c) = Cartesian (a*v) (b*v) (c*v)

origin :: Cartesian
origin = Cartesian 0 0 0

closestToOrigin :: Line -> Double
closestToOrigin (Line a b)
    | c1 <= 0   = norm origin a
    | c2 <= c1  = norm origin b
    | otherwise = norm origin p
    where
    v = b `sub` a
    w = origin `sub` a
    c1 = w • v
    c2 = v • v
    norm a b = sqrt (square (a `sub` b))
    square v = v • v
    bs = c1 / c2
    p = a `plus` (scale bs v)

-- Assumes sphere origin is at (0,0,0) with radius of 6371
earthIntersects :: Line -> Bool
earthIntersects line = closestToOrigin line < 6371

parseSat :: [String] -> (ID, Cartesian)
parseSat ['S':'A':'T':id,lat,lon,alt] = (ID (read id), p2c polar)
    where polar = Polar (read lat) (read lon) (read alt + 6371)

parseEnd :: [String] -> (Cartesian, Cartesian)
parseEnd ["ROUTE", lat1, lon1, lat2, lon2] = (p2c polar1, p2c polar2)
    where polar1 = Polar (read lat1) (read lon1) 6371.01 -- Assume 10 meters agl
          polar2 = Polar (read lat2) (read lon2) 6371.01 -- Assume 10 meters agl

main = do
    let process = tail . map (splitOn ",") . lines
    input <- process <$> readFile "data.csv"
    let sats = map parseSat $ init input
    let (start, end) = parseEnd $ last input
    let locs = fromList $ (Start,start):(End,end):sats
    let visibility = visibleFrom locs
    let startTree = treeFrom visibility Start
    let shortest = shortestPathTo End startTree
    case shortest of
        Nothing -> print "There is no path."
        Just (Path steps) -> do
            let format (ID id) = "SAT" ++ show id ++ ","
            putStrLn . init . concatMap format . tail . reverse . tail $ steps