diff --git a/README.md b/README.md
--- a/README.md
+++ b/README.md
@@ -2,11 +2,195 @@
 Telescope
 =========
 
-
 [![Hackage](https://img.shields.io/hackage/v/telescope.svg&color=success)](https://hackage.haskell.org/package/telescope)
 
 Haskell library to read and write Astronomical images from telescopes
 
-* [FITS](https://fits.gsfc.nasa.gov/fits_standard.html) File Format
-* [ASDF](https://asdf-standard.readthedocs.io/) File Format
+* [FITS](https://fits.gsfc.nasa.gov/fits_standard.html) (Flexible Image Transport System) Files
+* [ASDF](https://asdf-standard.readthedocs.io/) (Advanced Scientific Data Format) Files
 
+FITS Example
+-------
+
+```
+import Telescope.Fits
+import Data.ByteString qualified as BS
+
+main :: IO ()
+main = do
+  inp <- BS.readFile "hubble.fits"
+  fits <- decode inp
+```
+
+FITS files start with a primary *Header Data Unit*
+
+```
+  print fits.primaryHDU
+```
+
+```
+DataHDU
+  Header: 490 records
+  data: 0 bytes
+  dimensions:
+    format: Int8
+    axes: []
+```
+
+HDUs may contain both data and a *header*
+
+```
+  print fits.primaryHDU.header
+```
+
+```
+SIMPLE  =                    T / conforms to FITS standard
+BITPIX  =                    8 / array data type
+NAXIS   =                    0 / number of array dimensions
+EXTEND  =                    T
+DATE    = '2009-11-10'         / date this file was written (yyyy-mm-dd)
+FILETYPE= 'SCI'                / type of data found in data file
+
+TELESCOP= 'HST'                / telescope used to acquire data
+INSTRUME= 'WFPC2'              / identifier for instrument used to acquire data
+EQUINOX =               2000.0 / equinox of celestial coord. system
+...
+```
+
+The primary HDU is followed by multiple extensions in a similar format
+
+```
+  mapM_ print fits.extensions
+```
+
+```
+Image: DataHDU
+  Header: 103 records
+  data: 12960000 bytes
+  dimensions:
+    format: Float
+    axes: [1800,1800]
+
+Image: DataHDU
+  Header: 124 records
+  data: 12960000 bytes
+  dimensions:
+    format: Float
+    axes: [1800,1800]
+
+Image: DataHDU
+  Header: 123 records
+  data: 12960000 bytes
+  dimensions:
+    format: Int32
+    axes: [1800,1800]
+```
+
+We can use `Array` from [Data.Massiv](https://hackage.haskell.org/package/massiv) to parse, manipulate, and display data
+
+```
+  Image img : _ <- pure fits.extensions
+  arr <- decodeDataArray @Ix2 @Float img.dataArray
+  writeImage "hubble1.png" (heatmap arr)
+```
+
+
+![Hubble Output](https://raw.githubusercontent.com/DKISTDC/telescope.hs/main/example/output/hubble1.png)
+
+#### Parsing Headers
+
+Parse header keywords into Haskell records
+
+```
+import Telescope.Fits
+import Effectful
+
+data HubbleInfo = HubbleInfo
+  { telescop :: Text
+  , instrume :: Text
+  , equinox :: Float
+  } deriving (Generic, FromHeader, ToHeader)
+
+main = do
+  ...
+  let h = fits.primaryHDU.header
+  print $ lookupKeyword "INSTRUME" h
+
+  case runPureEff $ runParser $ parseHeader h of
+    Left e -> fail $ show e
+    Right info -> do
+      print info.telescop
+      print info.equinox
+```
+
+ASDF Example
+------------
+
+ASDF files contain a YAML tree and binary data afterwards. They are readable and the tree is editable in a text editor.
+
+```
+example :: IO ()
+example = do
+  inp <- BS.readFile "../samples/dkist.asdf"
+  dset :: Dataset <- decodeM inp
+  print dset.info
+
+data DKISTInfo = DKISTInfo
+  { instrumentName :: Text
+  , frameCount :: Int
+  }
+  deriving (Show, Generic, FromAsdf)
+
+
+data Dataset = Dataset
+  { info :: DKISTInfo
+  }
+instance FromAsdf Dataset where
+  -- parse .dataset.meta.inventory into DKISTInfo
+  parseValue val =
+    case val of
+      Object o -> do
+        dset <- o .: "dataset"
+        meta <- dset .: "meta"
+        info <- meta .: "inventory"
+        pure $ Dataset info
+      k -> expected "Object" k
+```
+
+
+```
+DKISTInfo {instrumentName = "VISP", frameCount = 1960}
+```
+
+Data can be parsed directly from the YAML tree or from binary-encoded data blocks (an NDArray)
+
+```
+data Example = Example
+  { foo :: Int32
+  , powers :: Maybe Powers
+  , sequence :: [Int64]
+  , random :: Array D Ix1 Double
+  }
+  deriving (Generic, FromAsdf)
+
+example :: IO ()
+example = do
+  inp <- BS.readFile "../samples/example.asdf"
+  ex :: Example <- decodeM inp
+  print ex.name
+  print ex.items
+  print $ take 30 ex.sequence
+  print $ take 10 $ M.toList ex.random
+```
+
+```
+"Monty"
+["one","two","three","four","five"]
+[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29]
+[0.7842310832387069,0.2279459557291822,0.9534462812074139,0.5100515929833191,0.6597920763222204,0.8778040169160855,0.8079416746447109,0.5373925949744411,0.5169365152585088,0.436101639340324]
+```
+
+
+### Collaborators
+
+* The FITS code was developed in collaboration with @krakrjak
diff --git a/example/output/hubble1.png b/example/output/hubble1.png
new file mode 100644
Binary files /dev/null and b/example/output/hubble1.png differ
diff --git a/src/Telescope/Asdf.hs b/src/Telescope/Asdf.hs
--- a/src/Telescope/Asdf.hs
+++ b/src/Telescope/Asdf.hs
@@ -1,31 +1,78 @@
+{- |
+Module:      Telescope.Asdf
+Copyright:   (c) 2024 Sean Hess
+License:     BSD3
+Maintainer:  Sean Hess <shess@nso.edu>
+Stability:   experimental
+Portability: portable
+
+Read, and Write ASDF (Advanced Scientific Data Format) files
+
+> import Data.ByteString qualified as BS
+> import Telescope.Asdf
+>
+> data Example = Example
+>   { name :: Text
+>   , items :: [Text]
+>   , sequence :: [Int64]
+>   , random :: Array D Ix1 Double
+>   }
+>   deriving (Generic, FromAsdf)
+>
+> example :: IO ()
+> example = do
+>   inp <- BS.readFile "samples/example.asdf"
+>   ex :: Example <- decodeM inp
+>   print ex.name
+>   print ex.items
+>   print $ take 30 ex.sequence
+>   print $ take 10 $ M.toList ex.random
+-}
 module Telescope.Asdf
-  ( ToAsdf (..)
-  , FromAsdf (..)
+  ( -- * Encoding
+    encodeM
+  , encode
+  , ToAsdf (..)
+
+    -- * Decoding
   , decodeM
   , decodeEither
   , decode
-  , encodeM
-  , encode
+  , FromAsdf (..)
+  , (.:)
+  , (.:?)
   , AsdfError
+  , Parser
+
+    -- * Binary Data
   , FromNDArray (..)
   , ToNDArray (..)
-  , SchemaTag
+  , NDArrayData (..)
+
+    -- * ASDF Tree
+  , Asdf (..)
   , Node (..)
   , Value (..)
   , Key
   , Object
   , fromValue
-  , NDArrayData (..)
+  , SchemaTag
+
+    -- * JSON Reference
   , jsonPointer
   , jsonReference
   , JSONReference (..)
   , JSONPointer (..)
+
+    -- * YAML Anchors
   , Anchor (..)
-  , Parser
-  , Asdf (..)
+
+    -- ** Exports
+  , Generic
   )
 where
 
+import GHC.Generics
 import Telescope.Asdf.Class
 import Telescope.Asdf.Core (Asdf (..))
 import Telescope.Asdf.Encoding
diff --git a/src/Telescope/Asdf/Class.hs b/src/Telescope/Asdf/Class.hs
--- a/src/Telescope/Asdf/Class.hs
+++ b/src/Telescope/Asdf/Class.hs
@@ -2,7 +2,6 @@
 
 module Telescope.Asdf.Class where
 
-import Data.List ((!?))
 import Data.List.NonEmpty (NonEmpty (..))
 import Data.List.NonEmpty qualified as NE
 import Data.Massiv.Array (Array, Prim)
@@ -170,6 +169,12 @@
   toValue as = toValue $ NE.toList as
 
 
+instance {-# OVERLAPS #-} FromAsdf Object where
+  parseValue = \case
+    Object o -> pure o
+    node -> expected "Object" node
+
+
 instance (ToAsdf a, ToAsdf b) => ToAsdf (a, b) where
   toValue (a, b) = Array [toNode a, toNode b]
 instance (FromAsdf a, FromAsdf b) => FromAsdf (a, b) where
@@ -344,6 +349,7 @@
 
 
 instance ToAsdf UTCTime where
+  schema _ = "!time/time-1.1.0"
   toValue t = String $ pack $ iso8601Show t
 instance FromAsdf UTCTime where
   parseValue v = do
@@ -398,10 +404,9 @@
 -}
 (!) :: (FromAsdf a, Parser :> es) => [Node] -> Int -> Eff es a
 ns ! n = do
-  case ns !? n of
-    Nothing -> parseFail $ "Index " ++ show n ++ " not found"
-    Just node ->
-      parseAt (Index n) $ parseNode node
+  if n >= length ns
+    then parseFail $ "Index " ++ show n ++ " not found"
+    else parseAt (Index n) $ parseNode $ ns !! n
 
 
 -- | Generically serialize records to an 'Object'
diff --git a/src/Telescope/Asdf/Core.hs b/src/Telescope/Asdf/Core.hs
--- a/src/Telescope/Asdf/Core.hs
+++ b/src/Telescope/Asdf/Core.hs
@@ -24,6 +24,9 @@
   | Pixel
   | Degrees
   | Nanometers
+  | Meters
+  | Kilometers
+  | Arcseconds
   | Unit Text
   deriving (Eq)
 
@@ -35,6 +38,9 @@
     Pixel -> "pixel"
     Degrees -> "deg"
     Nanometers -> "nm"
+    Kilometers -> "km"
+    Arcseconds -> "arcsec"
+    Meters -> "m"
     (Unit t) -> String t
 instance FromAsdf Unit where
   parseValue = \case
@@ -43,6 +49,9 @@
     String "pixel" -> pure Pixel
     String "pix" -> pure Pixel
     String "nm" -> pure Nanometers
+    String "km" -> pure Kilometers
+    String "m" -> pure Meters
+    String "arcsec" -> pure Arcseconds
     String t -> pure $ Unit t
     val -> expected "String" val
 
diff --git a/src/Telescope/Asdf/Encoding.hs b/src/Telescope/Asdf/Encoding.hs
--- a/src/Telescope/Asdf/Encoding.hs
+++ b/src/Telescope/Asdf/Encoding.hs
@@ -96,7 +96,7 @@
 streamAsdfFile :: (Error AsdfError :> es, IOE :> es) => Encoded Tree -> [Encoded Block] -> Eff es (Object, Anchors)
 streamAsdfFile (Encoded inp) ebks = do
   blocks <- mapM decodeBlock ebks
-  runParseError . runYamlError $ do
+  runYamlError $ do
     runReader blocks . runState @Anchors mempty . runResource . runConduit $ Yaml.decode inp .| sinkTree
 
 
@@ -118,5 +118,7 @@
 runYamlError = runErrorNoCallStackWith @YamlError (throwError . YamlError . show)
 
 
-runParseError :: (Error AsdfError :> es) => Eff (Error ParseError : es) a -> Eff es a
-runParseError = runErrorNoCallStackWith @ParseError (throwError . ParseError . show)
+runParseError :: (Error AsdfError :> es) => Eff es (Either ParseError a) -> Eff es a
+runParseError eff = do
+  res <- eff
+  either (throwError . ParseError . show) pure res
diff --git a/src/Telescope/Asdf/Encoding/Stream.hs b/src/Telescope/Asdf/Encoding/Stream.hs
--- a/src/Telescope/Asdf/Encoding/Stream.hs
+++ b/src/Telescope/Asdf/Encoding/Stream.hs
@@ -4,7 +4,6 @@
 import Data.ByteString (ByteString)
 import Data.Conduit.Combinators (peek)
 import Data.Conduit.Combinators qualified as C
-import Data.List ((!?))
 import Data.String (fromString)
 import Data.Text (pack, unpack)
 import Data.Text qualified as T
@@ -21,7 +20,7 @@
 import Telescope.Asdf.NDArray (NDArrayData (..))
 import Telescope.Asdf.Node
 import Telescope.Data.Axes
-import Telescope.Data.Parser (runPureParser)
+import Telescope.Data.Parser (runParser)
 import Text.Libyaml (Event (..), MappingStyle (..), SequenceStyle (..), Style (..), Tag (..))
 import Text.Libyaml qualified as Yaml
 import Text.Read (readMaybe)
@@ -232,14 +231,17 @@
     case val of
       Integer s -> do
         blocks <- ask
-        case blocks !? fromIntegral s of
-          Nothing -> throwError $ NDArrayMissingBlock s
-          Just (BlockData b) -> pure b
+        if fromIntegral s >= length blocks
+          then throwError $ NDArrayMissingBlock s
+          else do
+            let BlockData b = blocks !! fromIntegral s
+            pure b
       _ -> throwError $ NDArrayExpected "Source" val
 
   parseLocal :: (FromAsdf a, Error YamlError :> es) => String -> Value -> Eff es a
-  parseLocal expected val =
-    case runPureParser . runReader @Anchors mempty $ parseValue val of
+  parseLocal expected val = do
+    res <- runParser . runReader @Anchors mempty $ parseValue val
+    case res of
       Left _ -> throwError $ NDArrayExpected expected val
       Right a -> pure a
 
diff --git a/src/Telescope/Asdf/GWCS.hs b/src/Telescope/Asdf/GWCS.hs
--- a/src/Telescope/Asdf/GWCS.hs
+++ b/src/Telescope/Asdf/GWCS.hs
@@ -11,6 +11,7 @@
 import Data.String (IsString)
 import Data.Text (Text, pack)
 import Data.Text qualified as T
+import Data.Time.Clock (UTCTime)
 import GHC.Generics
 import Telescope.Asdf
 import Telescope.Asdf.Core
@@ -278,12 +279,12 @@
       ]
 
 
-data CelestialFrame = CelestialFrame
+data CelestialFrame ref = CelestialFrame
   { name :: Text
   , axes :: NonEmpty FrameAxis
-  , referenceFrame :: ICRSFrame
+  , referenceFrame :: ref
   }
-instance ToAsdf CelestialFrame where
+instance (ToAsdf ref) => ToAsdf (CelestialFrame ref) where
   schema _ = "tag:stsci.edu:gwcs/celestial_frame-1.0.0"
   toValue f =
     Object $
@@ -316,6 +317,89 @@
 instance ToAsdf ICRSFrame where
   schema _ = "tag:astropy.org:astropy/coordinates/frames/icrs-1.1.0"
   toValue _ = Object [("frame_attributes", toNode $ Object mempty)]
+
+
+data HelioprojectiveFrame = HelioprojectiveFrame
+  { coordinates :: Cartesian3D
+  , obstime :: UTCTime
+  , rsun :: Quantity
+  }
+instance ToAsdf HelioprojectiveFrame where
+  schema _ = "tag:sunpy.org:sunpy/coordinates/frames/helioprojective-1.0.0"
+  toValue frame =
+    Object
+      [("frame_attributes", fromValue attributes)]
+   where
+    observer = HelioObserver (CartesianRepresentation frame.coordinates) observation
+    observation = HelioObservation frame.obstime frame.rsun
+    attributes =
+      Object [("observer", toNode observer)] <> toValue observation
+
+
+-- reference_frame: !<tag:sunpy.org:sunpy/coordinates/frames/helioprojective-1.0.0>
+--   frame_attributes:
+--     observer: !<tag:sunpy.org:sunpy/coordinates/frames/heliographic_stonyhurst-1.1.0>
+--       data: !<tag:astropy.org:astropy/coordinates/representation-1.0.0>
+--         components:
+--           x: !unit/quantity-1.1.0 {unit: !unit/unit-1.0.0 m, value: 151741639088.53842}
+--           y: !unit/quantity-1.1.0 {unit: !unit/unit-1.0.0 m, value: 5050819.209579468}
+--           z: !unit/quantity-1.1.0 {unit: !unit/unit-1.0.0 m, value: -968701275.8656464}
+--         type: CartesianRepresentation
+--       frame_attributes:
+--         obstime: !time/time-1.1.0 2022-06-03T17:52:20.031
+--         rsun: !unit/quantity-1.1.0 {unit: !unit/unit-1.0.0 km, value: 695700.0}
+--     obstime: !time/time-1.1.0 2022-06-03T17:52:20.031
+--     rsun: !unit/quantity-1.1.0 {unit: !unit/unit-1.0.0 km, value: 695700.0}
+
+data HelioObserver = HelioObserver
+  { coordinates :: CartesianRepresentation Cartesian3D
+  , observation :: HelioObservation
+  }
+instance ToAsdf HelioObserver where
+  --     observer: !<tag:sunpy.org:sunpy/coordinates/frames/heliographic_stonyhurst-1.1.0>
+  --       data: !<tag:astropy.org:astropy/coordinates/representation-1.0.0>
+  --         components:
+  --           x: !unit/quantity-1.1.0 {unit: !unit/unit-1.0.0 m, value: 151741639088.53842}
+  --           y: !unit/quantity-1.1.0 {unit: !unit/unit-1.0.0 m, value: 5050819.209579468}
+  --           z: !unit/quantity-1.1.0 {unit: !unit/unit-1.0.0 m, value: -968701275.8656464}
+  --         type: CartesianRepresentation
+  --       frame_attributes:
+  --         obstime: !time/time-1.1.0 2022-06-03T17:52:20.031
+  --         rsun: !unit/quantity-1.1.0 {unit: !unit/unit-1.0.0 km, value: 695700.0}
+  schema _ = "tag:sunpy.org:sunpy/coordinates/frames/heliographic_stonyhurst-1.1.0"
+  toValue obs =
+    Object
+      [ ("data", toNode obs.coordinates)
+      , ("frame_attributes", toNode obs.observation)
+      ]
+
+
+data HelioObservation = HelioObservation
+  { obstime :: UTCTime
+  , rsun :: Quantity
+  }
+  deriving (Generic, ToAsdf)
+
+
+data Cartesian3D = Cartesian3D
+  { x :: Quantity
+  , y :: Quantity
+  , z :: Quantity
+  }
+  deriving (Generic, ToAsdf)
+
+
+data CartesianRepresentation dims = CartesianRepresentation dims
+instance (ToAsdf dims) => ToAsdf (CartesianRepresentation dims) where
+  -- data: !<tag:astropy.org:astropy/coordinates/representation-1.0.0>
+  --   components:
+  --     x: !unit/quantity-1.1.0 {unit: !unit/unit-1.0.0 m, value: 151741639088.53842}
+  --     y: !unit/quantity-1.1.0 {unit: !unit/unit-1.0.0 m, value: 5050819.209579468}
+  --     z: !unit/quantity-1.1.0 {unit: !unit/unit-1.0.0 m, value: -968701275.8656464}
+  --   type: CartesianRepresentation
+  schema _ = "tag:astropy.org:astropy/coordinates/representation-1.0.0"
+  toValue (CartesianRepresentation dims) =
+    Object [("components", toNode dims), ("type", "CartesianRepresentation")]
 
 
 data FrameAxis = FrameAxis
diff --git a/src/Telescope/Asdf/NDArray/Types.hs b/src/Telescope/Asdf/NDArray/Types.hs
--- a/src/Telescope/Asdf/NDArray/Types.hs
+++ b/src/Telescope/Asdf/NDArray/Types.hs
@@ -9,8 +9,6 @@
 import Telescope.Data.Binary
 
 
--- import Telescope.Asdf.Node
-
 {- | In-tree representation of an NDArray. You can parse a file as this and get it back. Not really what we want though
 but in haskell we can't easily just parse a multi-dimensional array
 we could do a simpler representation. Using an ADT
diff --git a/src/Telescope/Asdf/Reference.hs b/src/Telescope/Asdf/Reference.hs
--- a/src/Telescope/Asdf/Reference.hs
+++ b/src/Telescope/Asdf/Reference.hs
@@ -1,6 +1,5 @@
 module Telescope.Asdf.Reference where
 
-import Data.List ((!?))
 import Data.Text (Text)
 import Effectful
 import Telescope.Asdf.NDArray
@@ -37,10 +36,10 @@
       Just c -> pure c
 
   parseIndex :: Int -> [Node] -> Eff es Node
-  parseIndex n a =
-    case a !? n of
-      Nothing -> missingPointer (Index n) a
-      Just c -> pure c
+  parseIndex n nodes =
+    if n >= length nodes
+      then missingPointer (Index n) nodes
+      else pure $ nodes !! n
 
   missingPointer :: (Show ex, Show at) => ex -> at -> Eff es Node
   missingPointer expect at = parseFail $ "Could not locate pointer: " ++ show path ++ ". Expected " ++ show expect ++ " at " ++ show at
diff --git a/src/Telescope/Data/Array.hs b/src/Telescope/Data/Array.hs
--- a/src/Telescope/Data/Array.hs
+++ b/src/Telescope/Data/Array.hs
@@ -8,6 +8,8 @@
 import Data.ByteString.Lazy qualified as BL
 import Data.Massiv.Array (Array, Comp (..), D, Index, Ix1, Ix2 (..), Ix3, Ix4, Ix5, IxN (..), Lower, Prim, Source, Stream, Sz (..), Vector)
 import Data.Massiv.Array qualified as M
+import Data.Massiv.Array.IO (Linearity (..), Pixel (..), SRGB)
+import Data.Massiv.Array.IO qualified as M
 import Data.Word (Word8)
 import Telescope.Data.Axes
 import Telescope.Data.Binary
@@ -192,3 +194,16 @@
 -- | Get the Axes for a given array size
 sizeAxes :: (AxesIndex ix, Index ix) => Sz ix -> Axes Column
 sizeAxes (Sz ix) = toColumnMajor $ indexAxes ix
+
+
+-- approximates python's viridis
+heatmap :: forall n. (Ord n, RealFrac n) => Array D Ix2 n -> Array D Ix2 (Pixel (SRGB 'NonLinear) Word8)
+heatmap arr =
+  let maxn = maximum arr :: n
+   in fmap (color . (/ maxn)) arr
+ where
+  color x =
+    let r = floor (255 * abs (x - 0.3))
+        g = floor (255 * x)
+        b = floor (255 * (1 - x))
+     in M.PixelSRGB r g b
diff --git a/src/Telescope/Data/DataCube.hs b/src/Telescope/Data/DataCube.hs
new file mode 100644
--- /dev/null
+++ b/src/Telescope/Data/DataCube.hs
@@ -0,0 +1,150 @@
+{-# LANGUAGE UndecidableInstances #-}
+
+module Telescope.Data.DataCube where
+
+import Data.Kind
+import Data.Massiv.Array as M hiding (mapM)
+import Data.Proxy
+import GHC.TypeLits (natVal)
+import Telescope.Data.Array (AxesIndex (..))
+import Telescope.Data.Axes (Axes, Major (Row))
+
+
+-- Results ------------------------------------------------------------------------------
+
+newtype DataCube (as :: [Type]) = DataCube
+  { array :: Array D (IndexOf as) Double
+  }
+
+
+instance (Index (IndexOf as)) => Eq (DataCube as) where
+  DataCube arr == DataCube arr2 = arr == arr2
+
+
+instance (Ragged L (IndexOf as) Double) => Show (DataCube as) where
+  show (DataCube a) = show a
+
+
+class HasIndex (as :: [Type]) where
+  type IndexOf as :: Type
+
+
+instance HasIndex '[] where
+  type IndexOf '[] = Ix0
+instance HasIndex '[a] where
+  type IndexOf '[a] = Ix1
+instance HasIndex '[a, b] where
+  type IndexOf '[a, b] = Ix2
+instance HasIndex '[a, b, c] where
+  type IndexOf '[a, b, c] = Ix3
+instance HasIndex '[a, b, c, d] where
+  type IndexOf '[a, b, c, d] = Ix4
+instance HasIndex '[a, b, c, d, e] where
+  type IndexOf '[a, b, c, d, e] = Ix5
+
+
+outerList
+  :: forall a as
+   . (Lower (IndexOf (a : as)) ~ IndexOf as, Index (IndexOf as), Index (IndexOf (a : as)))
+  => DataCube (a : as)
+  -> [DataCube as]
+outerList (DataCube a) = foldOuterSlice row a
+ where
+  row :: Array D (IndexOf as) Double -> [DataCube as]
+  row r = [DataCube r]
+
+
+transposeMajor
+  :: (IndexOf (a : b : xs) ~ IndexOf (b : a : xs), Index (Lower (IndexOf (b : a : xs))), Index (IndexOf (b : a : xs)))
+  => DataCube (a : b : xs)
+  -> DataCube (b : a : xs)
+transposeMajor (DataCube arr) = DataCube $ transposeInner arr
+
+
+transposeMinor4
+  :: DataCube [a, b, c, d]
+  -> DataCube [a, b, d, c]
+transposeMinor4 (DataCube arr) = DataCube $ transposeOuter arr
+
+
+transposeMinor3
+  :: DataCube [a, b, c]
+  -> DataCube [a, c, b]
+transposeMinor3 (DataCube arr) = DataCube $ transposeOuter arr
+
+
+-- Slice along the 1st major dimension
+sliceM0
+  :: ( Lower (IndexOf (a : xs)) ~ IndexOf xs
+     , Index (IndexOf xs)
+     , Index (IndexOf (a : xs))
+     )
+  => Int
+  -> DataCube (a : xs)
+  -> DataCube xs
+sliceM0 a (DataCube arr) = DataCube (arr !> a)
+
+
+-- Slice along the 2nd major dimension
+sliceM1
+  :: forall a b xs
+   . ( Lower (IndexOf (a : b : xs)) ~ IndexOf (a : xs)
+     , Index (IndexOf (a : xs))
+     , Index (IndexOf (a : b : xs))
+     )
+  => Int
+  -> DataCube (a : b : xs)
+  -> DataCube (a : xs)
+sliceM1 b (DataCube arr) =
+  let dims = fromIntegral $ natVal @(Dimensions (IndexOf (a : b : xs))) Proxy
+   in DataCube $ arr <!> (Dim (dims - 1), b)
+
+
+-- Slice along the 3rd major dimension
+sliceM2
+  :: forall a b c xs
+   . ( Lower (IndexOf (a : b : c : xs)) ~ IndexOf (a : b : xs)
+     , Index (IndexOf (a : b : xs))
+     , Index (IndexOf (a : b : c : xs))
+     )
+  => Int
+  -> DataCube (a : b : c : xs)
+  -> DataCube (a : b : xs)
+sliceM2 c (DataCube arr) =
+  let dims = fromIntegral $ natVal @(Dimensions (IndexOf (a : b : c : xs))) Proxy
+   in DataCube $ arr <!> (Dim (dims - 2), c)
+
+
+splitM0
+  :: forall a xs m
+   . ( Index (IndexOf (a : xs))
+     , MonadThrow m
+     )
+  => Int
+  -> DataCube (a : xs)
+  -> m (DataCube (a : xs), DataCube (a : xs))
+splitM0 a (DataCube arr) = do
+  let dims = fromIntegral $ natVal @(Dimensions (IndexOf (a : xs))) Proxy
+  (arr1, arr2) <- M.splitAtM (Dim dims) a arr
+  pure (DataCube arr1, DataCube arr2)
+
+
+splitM1
+  :: forall a b xs m
+   . ( Index (IndexOf (a : xs))
+     , Index (IndexOf (a : b : xs))
+     , MonadThrow m
+     )
+  => Int
+  -> DataCube (a : b : xs)
+  -> m (DataCube (a : b : xs), DataCube (a : b : xs))
+splitM1 b (DataCube arr) = do
+  let dims = fromIntegral $ natVal @(Dimensions (IndexOf (a : xs))) Proxy
+  (arr1, arr2) <- M.splitAtM (Dim dims) b arr
+  pure (DataCube arr1, DataCube arr2)
+
+
+dataCubeAxes :: (Index (IndexOf as), AxesIndex (IndexOf as)) => DataCube as -> Axes Row
+dataCubeAxes (DataCube arr) =
+  let Sz ix = M.size arr
+   in indexAxes ix
diff --git a/src/Telescope/Data/Parser.hs b/src/Telescope/Data/Parser.hs
--- a/src/Telescope/Data/Parser.hs
+++ b/src/Telescope/Data/Parser.hs
@@ -1,3 +1,5 @@
+{-# LANGUAGE UndecidableInstances #-}
+
 module Telescope.Data.Parser where
 
 import Control.Monad.Catch (Exception)
@@ -11,28 +13,25 @@
 
 data Parser :: Effect where
   ParseFail :: String -> Parser m a
-  PathAdd :: Ref -> m a -> Parser m a
+  PathMod :: (Path -> Path) -> m a -> Parser m a
 
 
 type instance DispatchOf Parser = 'Dynamic
 
 
 runParser
-  :: (Error ParseError :> es)
-  => Eff (Parser : es) a
-  -> Eff es a
-runParser = reinterpret (runReader @Path mempty) $ \env -> \case
+  :: Eff (Parser : es) a
+  -> Eff es (Either ParseError a)
+runParser = reinterpret (runErrorNoCallStack @ParseError . runReader @Path mempty) $ \env -> \case
   ParseFail e -> do
     path <- ask @Path
     throwError $ ParseFailure path e
-  PathAdd p m -> do
-    -- copied from Effectful.Reader.Dynamic
-    localSeqUnlift env $ \unlift -> local (<> Path [p]) (unlift m)
+  PathMod mp m -> do
+    localSeqUnlift env $ \unlift -> local mp (unlift m)
 
 
-runPureParser :: Eff '[Parser, Error ParseError] a -> Either ParseError a
-runPureParser eff = runPureEff . runErrorNoCallStack @ParseError $ runParser eff
-
+-- copied from Effectful.Reader.Dynamic
+-- localSeqUnlift env $ \unlift -> local (<> Path [p]) (unlift m)
 
 data ParseError
   = ParseFailure Path String
@@ -85,4 +84,4 @@
 -- | Add a child to the parsing 'Path'
 parseAt :: (Parser :> es) => Ref -> Eff es a -> Eff es a
 parseAt p parse = do
-  send $ PathAdd p parse
+  send $ PathMod (<> Path [p]) parse
diff --git a/src/Telescope/Data/WCS.hs b/src/Telescope/Data/WCS.hs
--- a/src/Telescope/Data/WCS.hs
+++ b/src/Telescope/Data/WCS.hs
@@ -2,7 +2,7 @@
 
 module Telescope.Data.WCS where
 
-import Data.Text
+import Data.Text (Text, pack)
 import GHC.Generics (Generic)
 import Telescope.Data.Axes (AxisOrder (..))
 import Telescope.Data.KnownText
diff --git a/src/Telescope/Fits.hs b/src/Telescope/Fits.hs
--- a/src/Telescope/Fits.hs
+++ b/src/Telescope/Fits.hs
@@ -6,7 +6,7 @@
 Stability:   experimental
 Portability: portable
 
-Read, Generate, and Write FITS (Flexible Image Transport System) files
+Read, and Write FITS (Flexible Image Transport System) files
 
 @
 import Data.ByteString qualified as BS
@@ -16,8 +16,7 @@
 test = do
   inp <- BS.readFile "samples/simple2x3.fits"
   f <- decode inp
-  print f.primaryHDU.dataArray.axes
-  print f.primaryHDU.dataArray.bitpix
+  print f.primaryHDU
   print $ lookupKeyword \"BTYPE\" f.primaryHDU.header
 
   a <- decodeArray @Ix2 @Int f.primaryHDU.dataArray
@@ -40,16 +39,17 @@
     -- * Parsing Headers
   , FromHeader (..)
   , FromKeyword (..)
+  , Parser
+  , runParser
 
-    -- * Creating Headers
+    -- * Writing Headers
   , ToKeyword (..)
   , ToHeader (..)
-  , Parser
 
     -- * Types
   , Fits (..)
-  , PrimaryHDU (..)
-  , ImageHDU (..)
+  , KeywordRecord (..)
+  , DataHDU (..)
   , BinTableHDU (..)
   , DataArray (..)
   , Extension (..)
@@ -60,12 +60,16 @@
 
     -- * Generate
   , addComment
-  , keyword
   , dataArray
   , emptyDataArray
 
+    -- * Visualize
+  , heatmap
+  , writeImage
+
     -- * Exports from Data.Massiv.Array
   , Array
+  , D
   , Ix1
   , Ix2
   , Ix3
@@ -81,10 +85,14 @@
   -- , test
   ) where
 
-import Data.Massiv.Array (Array, Dim (..), Ix1, Ix2, Ix3, Ix4, Ix5, size, (!>), (!?>), (<!), (<!>), (<!?))
-import Telescope.Data.Parser (Parser)
+import Data.Massiv.Array (Array, D, Dim (..), Ix1, Ix2, Ix3, Ix4, Ix5, size, (!>), (!?>), (<!), (<!>), (<!?))
+import Data.Massiv.Array.IO (writeImage)
+import Telescope.Data.Array (heatmap)
+import Telescope.Data.Axes
+import Telescope.Data.Parser (Parser, runParser)
+import Telescope.Fits.BitPix
 import Telescope.Fits.DataArray
 import Telescope.Fits.Encoding
-import Telescope.Fits.Header (FromHeader (..), FromKeyword (..), ToHeader (..), ToKeyword (..), addComment, keyword, lookupKeyword)
-import Telescope.Fits.Types
+import Telescope.Fits.HDU
+import Telescope.Fits.Header
 
diff --git a/src/Telescope/Fits/BitPix.hs b/src/Telescope/Fits/BitPix.hs
new file mode 100644
--- /dev/null
+++ b/src/Telescope/Fits/BitPix.hs
@@ -0,0 +1,49 @@
+{-# LANGUAGE AllowAmbiguousTypes #-}
+
+module Telescope.Fits.BitPix where
+
+import GHC.Int
+import Prelude hiding (lookup)
+
+
+data BitPix
+  = BPInt8
+  | BPInt16
+  | BPInt32
+  | BPInt64
+  | BPFloat
+  | BPDouble
+  deriving (Show, Eq)
+
+
+bitPixBits :: BitPix -> Int
+bitPixBits bp = bitPixBytes bp * 8
+
+
+bitPixBytes :: BitPix -> Int
+bitPixBytes bp =
+  case bp of
+    BPInt8 -> 1
+    BPInt16 -> 2
+    BPInt32 -> 4
+    BPInt64 -> 8
+    BPFloat -> 4
+    BPDouble -> 8
+
+
+class IsBitPix a where
+  bitPix :: BitPix
+
+
+instance IsBitPix Int8 where
+  bitPix = BPInt8
+instance IsBitPix Int16 where
+  bitPix = BPInt16
+instance IsBitPix Int32 where
+  bitPix = BPInt32
+instance IsBitPix Int64 where
+  bitPix = BPInt64
+instance IsBitPix Float where
+  bitPix = BPFloat
+instance IsBitPix Double where
+  bitPix = BPDouble
diff --git a/src/Telescope/Fits/Checksum.hs b/src/Telescope/Fits/Checksum.hs
--- a/src/Telescope/Fits/Checksum.hs
+++ b/src/Telescope/Fits/Checksum.hs
@@ -2,7 +2,6 @@
 
 import Data.Bits (complement, shiftR, (.&.))
 import Data.ByteString.Internal
-import Data.Fits (Value (..))
 import Data.Text (Text, pack)
 import Data.Word
 import Foreign.C.String
@@ -10,6 +9,7 @@
 import Foreign.ForeignPtr (withForeignPtr)
 import Foreign.Ptr
 import GHC.IO
+import Telescope.Fits.Header.Value
 
 
 -- | Generate the Checksum per the FITS spec
diff --git a/src/Telescope/Fits/DataArray.hs b/src/Telescope/Fits/DataArray.hs
--- a/src/Telescope/Fits/DataArray.hs
+++ b/src/Telescope/Fits/DataArray.hs
@@ -1,26 +1,53 @@
-module Telescope.Fits.DataArray
-  ( DataArray (..)
-  , dataArray
-  , decodeDataArray
-  , encodeDataArray
-  )
-where
+module Telescope.Fits.DataArray where
 
 import Control.Monad.Catch (MonadCatch)
 import Data.ByteString (ByteString)
-import Data.Fits qualified as Fits
-import Data.Massiv.Array as M hiding (isEmpty, product)
+import Data.ByteString qualified as BS
+import Data.List qualified as L
+import Data.Massiv.Array as M hiding (Dimensions, isEmpty, product)
 import System.ByteOrder
 import Telescope.Data.Array
 import Telescope.Data.Axes
 import Telescope.Data.Binary
-import Telescope.Fits.Types
+import Telescope.Fits.BitPix
 
 
+data Dimensions = Dimensions
+  { bitpix :: BitPix
+  , axes :: Axes Column
+  }
+  deriving (Show, Eq)
+
+
+dataSizeBytes :: Dimensions -> Int
+dataSizeBytes (Dimensions bitpix axes) =
+  bitPixBytes bitpix * count axes
+ where
+  count (Axes []) = 0
+  count (Axes ax) = fromIntegral $ product ax
+
+
+-- | Raw HDU Data. See 'Telescope.Fits.DataArray'
+data DataArray = DataArray
+  { bitpix :: BitPix
+  , axes :: Axes Column
+  , rawData :: BS.ByteString
+  }
+
+
+instance Show DataArray where
+  show d =
+    L.intercalate
+      "\n"
+      [ "  data: " <> show (BS.length d.rawData) <> " bytes"
+      , "  dimensions: "
+      , "    format: " <> L.drop 2 (show d.bitpix)
+      , "    axes: " <> show d.axes.axes
+      ]
+
+
 -- > {-# LANGUAGE TypeApplications #-}
 -- > import Data.Massiv.Array
--- > import Data.Fits.Image
--- > import Data.Fits
 -- >
 -- > decodeExample :: BL.ByteString -> Either String Int
 -- > decodeExample bs = do
@@ -64,22 +91,7 @@
    in DataArray{bitpix, axes, rawData}
 
 
--- | Create a DataArray from raw Fits info
-dataArray :: Fits.Dimensions -> ByteString -> DataArray
-dataArray dim dat =
-  DataArray
-    { bitpix = bitpix dim._bitpix
-    , axes = axes dim._axes
-    , rawData = dat
-    }
- where
-  bitpix :: Fits.BitPixFormat -> BitPix
-  bitpix Fits.EightBitInt = BPInt8
-  bitpix Fits.SixteenBitInt = BPInt16
-  bitpix Fits.ThirtyTwoBitInt = BPInt32
-  bitpix Fits.SixtyFourBitInt = BPInt64
-  bitpix Fits.ThirtyTwoBitFloat = BPFloat
-  bitpix Fits.SixtyFourBitFloat = BPDouble
-
-  axes :: Fits.Axes -> Axes Column
-  axes = Axes
+-- -- | Create a DataArray from raw Fits info
+dataArray :: Dimensions -> ByteString -> DataArray
+dataArray Dimensions{bitpix, axes} rawData =
+  DataArray{bitpix, axes, rawData}
diff --git a/src/Telescope/Fits/Encoding.hs b/src/Telescope/Fits/Encoding.hs
--- a/src/Telescope/Fits/Encoding.hs
+++ b/src/Telescope/Fits/Encoding.hs
@@ -1,72 +1,137 @@
-module Telescope.Fits.Encoding
-  ( -- * Decoding
-    decode
-
-    -- * Encoding
-  , encode
-  , encodePrimaryHDU
-  , encodeImageHDU
-  , encodeExtension
-  , encodeHDU
-  , replaceKeywordLine
-
-    -- * Parser
-  , nextParser
-  , parseFits
-  , parseMainData
-  , HDUError (..)
-  )
-where
+module Telescope.Fits.Encoding where
 
-import Control.Exception (Exception)
 import Control.Monad.Catch (MonadThrow, throwM)
 import Data.ByteString (ByteString)
 import Data.ByteString qualified as BS
-import Data.Fits qualified as Fits
-import Data.Fits.MegaParser qualified as Fits
-import Data.Fits.Read (FitsError (..))
 import Data.Text (Text)
 import Data.Text.Encoding qualified as TE
 import Effectful
-import Effectful.Error.Static
+import Effectful.Dispatch.Dynamic
 import Effectful.State.Static.Local
+import Telescope.Data.Parser
 import Telescope.Fits.Checksum
-import Telescope.Fits.DataArray (dataArray)
+import Telescope.Fits.DataArray (DataArray (..), Dimensions (..), dataArray, dataSizeBytes)
+import Telescope.Fits.Encoding.MegaHeader qualified as MH
 import Telescope.Fits.Encoding.Render
-import Telescope.Fits.Types
-import Text.Megaparsec qualified as M
-import Text.Megaparsec.State qualified as M
+import Telescope.Fits.HDU
+import Telescope.Fits.Header
+import Text.Megaparsec (lookAhead)
+import Text.Megaparsec qualified as MP
 
 
-{- | Decode a FITS file read as a strict 'ByteString'
+{- | Decode a FITS file read from a strict 'ByteString'
 
 >  decode =<< BS.readFile "samples/simple2x3.fits"
 -}
 decode :: forall m. (MonadThrow m) => ByteString -> m Fits
-decode inp = do
-  let res = runPureEff $ runErrorNoCallStack @HDUError $ evalState inp parseFits
-  case res of
-    Left e -> throwM e
-    Right a -> pure a
+decode = fitsParseThrow parseFits
 
 
 {- | Encode a FITS file to a strict 'ByteString'
 
-> BS.writeFile $ encdoe fits
+> BS.writeFile $ encode fits
 -}
 encode :: Fits -> ByteString
 encode f =
-  let primary = encodePrimaryHDU f.primaryHDU
+  let prim = encodePrimaryHDU f.primaryHDU
       exts = fmap encodeExtension f.extensions
-   in mconcat $ primary : exts
+   in mconcat $ prim : exts
 
 
-encodePrimaryHDU :: PrimaryHDU -> ByteString
+fitsParseThrow :: (MonadThrow m) => Eff '[State ByteString, Parser] a -> ByteString -> m a
+fitsParseThrow parse inp =
+  case runFitsParse parse inp of
+    Left e -> throwM e
+    Right a -> pure a
+
+
+runFitsParse :: Eff '[State ByteString, Parser] a -> BS.ByteString -> Either ParseError a
+runFitsParse parse inp = do
+  runPureEff $ runParser $ evalState inp parse
+
+
+parseFits :: forall es. (State ByteString :> es, Parser :> es) => Eff es Fits
+parseFits = do
+  p <- primary
+  es <- parseAt (Child "extensions") $ extensions 1
+  pure $ Fits p es
+
+
+primary :: (State ByteString :> es, Parser :> es) => Eff es DataHDU
+primary = do
+  parseAt (Child "primary") $ do
+    (dm, hd) <- runMega "Primary Header" $ do
+      dm <- MP.lookAhead MH.parsePrimaryKeywords
+      hd <- MH.parseHeader
+      pure (dm, hd)
+    darr <- mainData dm
+    pure $ DataHDU hd darr
+
+
+extensions :: (State ByteString :> es, Parser :> es) => Int -> Eff es [Extension]
+extensions n = do
+  inp <- get @ByteString
+  -- recurse until you run out of input or error
+  case inp of
+    "" -> pure []
+    _ -> do
+      e <- parseAt (Index n) extension
+      es <- extensions (n + 1)
+      pure (e : es)
+
+
+extension :: (State ByteString :> es, Parser :> es) => Eff es Extension
+extension = do
+  resImg <- runParser image
+  resTbl <- runParser binTable
+  case (resImg, resTbl) of
+    (Right i, _) -> pure $ Image i
+    (_, Right b) -> pure $ BinTable b
+    -- (Left (_, FormatError ie), Left (_, FormatError be)) -> throwM $ InvalidHDU [ie, be]
+    (Left _, Left (ParseFailure p e)) -> do
+      send $ PathMod (<> p) $ parseFail e
+
+
+image :: (State ByteString :> es, Parser :> es) => Eff es DataHDU
+image = do
+  (dm, hd) <- runMega "Image Header" $ do
+    dm <- MP.lookAhead MH.parseImageKeywords
+    hd <- MH.parseHeader
+    pure (dm, hd)
+  darr <- mainData dm
+  pure $ DataHDU hd darr
+
+
+binTable :: (State ByteString :> es, Parser :> es) => Eff es BinTableHDU
+binTable = do
+  (dm, pcount, hd) <- do
+    runMega "BinTable Header" $ do
+      (dm, pcount) <- lookAhead MH.parseBinTableKeywords
+      hd <- MH.parseHeader
+      pure (dm, pcount, hd)
+
+  darr <- mainData dm
+  rest <- get
+  let heap = BS.take pcount rest
+  put $ BS.dropWhile (== 0) $ BS.drop pcount rest
+  pure $ BinTableHDU hd pcount heap darr
+
+
+mainData :: (State ByteString :> es) => Dimensions -> Eff es DataArray
+mainData dm = do
+  rest <- get
+  let len = dataSizeBytes dm
+  let dat = dataArray dm (BS.take len rest)
+  put $ BS.dropWhile (== 0) $ BS.drop len rest
+  pure dat
+
+
+encodePrimaryHDU :: DataHDU -> ByteString
 encodePrimaryHDU p =
   encodeHDU (renderPrimaryHeader p.header p.dataArray) p.dataArray.rawData
 
 
-encodeImageHDU :: ImageHDU -> ByteString
+encodeImageHDU :: DataHDU -> ByteString
 encodeImageHDU p =
   encodeHDU (renderImageHeader p.header p.dataArray) p.dataArray.rawData
 
@@ -101,98 +166,13 @@
    in start <> newKeyLine <> BS.drop 80 rest
 
 
-parseFits :: forall es. (State ByteString :> es, Error HDUError :> es) => Eff es Fits
-parseFits = do
-  p <- primary
-  es <- extensions
-  pure $ Fits p es
- where
-  primary :: Eff es PrimaryHDU
-  primary = do
-    (dm, hd) <- nextParser "Primary Header" $ do
-      dm <- Fits.parsePrimaryKeywords
-      hd <- Fits.parseHeader
-      pure (dm, hd)
-
-    darr <- parseMainData dm
-    pure $ PrimaryHDU hd darr
-
-  extensions :: Eff es [Extension]
-  extensions = do
-    inp <- get @ByteString
-    case inp of
-      "" -> pure []
-      _ -> do
-        e <- extension
-        es <- extensions
-        pure (e : es)
-
-  extension :: Eff es Extension
-  extension = do
-    -- this consumes input!
-    resImg <- tryError @HDUError image
-    resTbl <- tryError @HDUError binTable
-    case (resImg, resTbl) of
-      (Right i, _) -> pure $ Image i
-      (_, Right b) -> pure $ BinTable b
-      (Left (_, FormatError ie), Left (_, FormatError be)) -> throwM $ InvalidHDU [ie, be]
-      (Left _, Left (_, be)) -> throwM be
-
-  image :: Eff es ImageHDU
-  image = do
-    (dm, hd) <- imageHeader
-    darr <- parseMainData dm
-    pure $ ImageHDU hd darr
-
-  imageHeader :: Eff es (Fits.Dimensions, Header)
-  imageHeader = do
-    nextParser "Image Header" $ do
-      dm <- Fits.parseImageKeywords
-      hd <- Fits.parseHeader
-      pure (dm, hd)
-
-  binTable :: Eff es BinTableHDU
-  binTable = do
-    (dm, pcount, hd) <- binTableHeader
-    darr <- parseMainData dm
-    rest <- get
-    let heap = BS.take pcount rest
-    put $ BS.dropWhile (== 0) $ BS.drop pcount rest
-    pure $ BinTableHDU hd pcount heap darr
-
-  binTableHeader :: Eff es (Fits.Dimensions, Int, Header)
-  binTableHeader = do
-    nextParser "BinTable Header" $ do
-      (dm, pcount) <- Fits.parseBinTableKeywords
-      hd <- Fits.parseHeader
-      pure (dm, pcount, hd)
-
-
-parseMainData :: (State ByteString :> es) => Fits.Dimensions -> Eff es DataArray
-parseMainData dm = do
-  rest <- get
-  let len = Fits.dataSize dm
-  let dat = dataArray dm (BS.take len rest)
-  put $ BS.dropWhile (== 0) $ BS.drop len rest
-  pure dat
-
-
--- | Parse HDUs by running MegaParsec parsers one at a time, tracking how much of the ByteString we've consumed
-nextParser :: (Error HDUError :> es, State ByteString :> es) => String -> Fits.Parser a -> Eff es a
-nextParser src parse = do
-  bs <- get
-  let st1 = M.initialState src bs
-  case M.runParser' parse st1 of
-    (st2, Right a) -> do
-      -- only consumes input if it succeeds
-      put $ BS.drop st2.stateOffset bs
+-- | Parse HDUs by running MegaParsec parsers one at a time, and tracking how much of the ByteString we've consumed
+runMega :: (Parser :> es, State ByteString :> es) => String -> MH.Parser a -> Eff es a
+runMega src parse = do
+  inp <- get
+  case MH.runNextParser src inp parse of
+    Right (a, rest) -> do
+      put rest
       pure a
-    (_, Left err) -> throwError $ FormatError $ ParseError err
-
-
-data HDUError
-  = InvalidExtension String
-  | MissingPrimary
-  | FormatError FitsError
-  | InvalidHDU [FitsError]
-  deriving (Show, Exception)
+    Left err ->
+      parseFail $ MH.showParseError err
diff --git a/src/Telescope/Fits/Encoding/MegaHeader.hs b/src/Telescope/Fits/Encoding/MegaHeader.hs
new file mode 100644
--- /dev/null
+++ b/src/Telescope/Fits/Encoding/MegaHeader.hs
@@ -0,0 +1,358 @@
+{-# LANGUAGE OverloadedStrings #-}
+
+-- {-# OPTIONS_HADDOCK hide #-}
+
+module Telescope.Fits.Encoding.MegaHeader where
+
+import Control.Monad (void)
+import Data.ByteString (ByteString)
+import Data.ByteString qualified as BS
+import Data.ByteString.Char8 qualified as C8
+import Data.Char (ord)
+import Data.List qualified as L
+import Data.List.NonEmpty qualified as NE
+import Data.Text (Text)
+import Data.Text qualified as T
+import Data.Text.Encoding qualified as TE
+import Data.Void (Void)
+import Data.Word (Word8)
+import Telescope.Data.Axes
+import Telescope.Fits.BitPix
+import Telescope.Fits.DataArray
+import Telescope.Fits.HDU.Block (hduRecordLength)
+import Telescope.Fits.Header hiding (FromKeyword (..), parseHeader, parseKeyword)
+import Text.Megaparsec (ParseErrorBundle, Parsec, (<|>))
+import Text.Megaparsec qualified as M
+import Text.Megaparsec.Byte qualified as M
+import Text.Megaparsec.Byte.Lexer qualified as MBL
+import Text.Megaparsec.Error (ParseError (..), ParseErrorBundle (..))
+import Text.Megaparsec.Error qualified as MP
+import Text.Megaparsec.Pos qualified as MP
+import Text.Megaparsec.State qualified as M
+
+
+type Parser = Parsec Void ByteString
+type ParseErr = ParseErrorBundle ByteString Void
+
+
+-- | Runs a single parser and returns the remainder of the input
+runNextParser :: String -> BS.ByteString -> Parser a -> Either ParseErr (a, BS.ByteString)
+runNextParser src inp parse = do
+  let st1 = M.initialState src inp
+  case M.runParser' parse st1 of
+    (st2, Right a) -> do
+      -- only consumes input if it succeeds
+      let rest = BS.drop st2.stateOffset inp
+      pure (a, rest)
+    (_, Left err) -> Left err
+
+
+showParseError :: ParseErr -> String
+showParseError bundle =
+  let pos = bundle.bundlePosState
+      inp = pos.pstateInput
+   in L.intercalate "\n HI " (NE.toList (fmap (showError inp) bundle.bundleErrors))
+ where
+  showError :: ByteString -> ParseError ByteString Void -> String
+  showError inp err =
+    showCurrent inp (MP.errorOffset err) <> "\n  " <> L.intercalate "\n  " (lines (MP.parseErrorPretty err))
+
+  showCurrent inp off =
+    let line = floor @Float @Int (fromIntegral off / fromIntegral hduRecordLength)
+        col = off - (line * hduRecordLength)
+     in "HDU Header "
+          <> show line
+          <> " column "
+          <> show col
+          <> ": \n"
+          <> show (BS.take hduRecordLength $ BS.drop (line * hduRecordLength) inp)
+
+
+toWord :: Char -> Word8
+toWord = fromIntegral . ord
+
+
+wordsText :: [Word8] -> Text
+wordsText = TE.decodeUtf8 . BS.pack
+
+
+-- | Consumes ALL header blocks until end, then all remaining space
+parseHeader :: Parser Header
+parseHeader = do
+  pairs <- M.manyTill parseRecordLine (M.string' "end")
+  M.space -- consume space padding all the way to the end of the next 2880 bytes header block
+  return $ Header pairs
+
+
+parseRecordLine :: Parser HeaderRecord
+parseRecordLine = do
+  M.try (Keyword <$> parseKeywordRecord)
+    <|> M.try (BlankLine <$ parseLineBlank)
+    <|> M.try (Comment <$> parseLineComment)
+    <|> M.try (Comment <$> parseLineCommentSpaces)
+    <|> (History <$> parseLineHistory)
+
+
+parseLineHistory :: Parser Text
+parseLineHistory = do
+  lineStart <- parsePos
+  _ <- M.string' "HISTORY "
+  M.space
+  untilLineEnd lineStart M.anySingle
+
+
+parseKeywordRecord :: Parser KeywordRecord
+parseKeywordRecord = do
+  ((k, v), mc) <- withComments parseKeywordValue
+  pure $ KeywordRecord k v mc
+
+
+-- | Parses the specified keyword
+parseKeywordRecord' :: ByteString -> Parser a -> Parser a
+parseKeywordRecord' k pval = ignoreComments $ do
+  _ <- M.string' k
+  parseEquals
+  pval
+
+
+-- | Combinator to allow for parsing a record with inline comments
+withComments :: Parser a -> Parser (a, Maybe Text)
+withComments parse = do
+  -- assumes we are at the beginning of the line
+  lineStart <- parsePos
+  a <- parse
+  mc <- parseLineEnd lineStart
+  return (a, mc)
+
+
+ignoreComments :: Parser a -> Parser a
+ignoreComments parse = do
+  (a, _) <- withComments parse
+  pure a
+
+
+parseKeywordValue :: Parser (Text, Value)
+parseKeywordValue = do
+  key <- parseKeyword
+  parseEquals
+  val <- parseValue
+  return (key, val)
+
+
+parseLineEnd :: Int -> Parser (Maybe Text)
+parseLineEnd lineStart = do
+  M.try (Nothing <$ spacesToLineEnd lineStart) <|> (Just <$> parseInlineComment lineStart)
+
+
+untilLineEnd :: Int -> Parser Word8 -> Parser Text
+untilLineEnd lineStart parseChar = do
+  curr <- parsePos
+  let used = curr - lineStart
+  bs <- M.count (hduRecordLength - used) parseChar
+  return $ wordsText bs
+
+
+spacesToLineEnd :: Int -> Parser ()
+spacesToLineEnd lineStart = do
+  _ <- untilLineEnd lineStart (M.char $ toWord ' ')
+  pure ()
+
+
+parseInlineComment :: Int -> Parser Text
+parseInlineComment lineStart = do
+  M.space
+  _ <- M.char $ toWord '/'
+  _ <- M.optional charSpace
+  T.strip <$> untilLineEnd lineStart M.anySingle
+ where
+  charSpace = M.char $ toWord ' '
+
+
+parseLineComment :: Parser Text
+parseLineComment = do
+  let kw = "COMMENT " :: ByteString
+  _ <- M.string' kw
+  c <- M.count (hduRecordLength - BS.length kw) M.anySingle
+  return $ wordsText c
+
+
+parseLineCommentSpaces :: Parser Text
+parseLineCommentSpaces = do
+  -- Invalid comments used by JWST
+  lineStart <- parsePos
+  -- exactly 8 spaces
+  _ <- M.count 8 (M.satisfy (== toWord ' '))
+  T.strip <$> untilLineEnd lineStart M.anySingle
+
+
+parseLineBlank :: Parser ()
+parseLineBlank = do
+  _ <- withComments $ M.count' 0 80 (M.satisfy (== toWord ' '))
+  pure ()
+
+
+-- | Anything but a space or equals
+parseKeyword :: Parser Text
+parseKeyword = wordsText <$> M.some (M.noneOf $ fmap toWord [' ', '='])
+
+
+parseValue :: Parser Value
+parseValue =
+  -- try is required here because Megaparsec doesn't automatically backtrack if the parser consumes anything
+  M.try (Float <$> parseFloat)
+    <|> M.try (Integer <$> parseInt)
+    <|> (Logic <$> parseLogic)
+    <|> (String <$> parseStringContinue)
+
+
+parseInt :: (Num a) => Parser a
+parseInt = MBL.signed M.space MBL.decimal
+
+
+parseFloat :: Parser Double
+parseFloat = MBL.signed M.space MBL.float
+
+
+parseLogic :: Parser LogicalConstant
+parseLogic = do
+  T <$ M.string' "T" <|> F <$ M.string' "F"
+
+
+parseStringContinue :: Parser Text
+parseStringContinue = do
+  t <- parseStringValue
+
+  mc <- M.optional $ M.try parseContinue
+
+  case mc of
+    Nothing -> return t
+    Just tc -> return $ T.dropWhileEnd (== '&') t <> tc
+
+
+parseContinue :: Parser Text
+parseContinue = do
+  M.space
+  lineStart <- parsePos
+  _ <- M.string' "CONTINUE"
+  M.space
+  more <- parseStringContinue
+  _ <- parseLineEnd lineStart
+  pure more
+
+
+parseStringValue :: Parser Text
+parseStringValue = do
+  -- The rules are weird, NULL means a NULL string, '' is an empty
+  -- string, a ' followed by a bunch of spaces and a close ' is
+  -- considered an empty string, and trailing whitespace is ignored
+  -- within the quotes, but not leading spaces.
+  ls <- M.between (M.char quote) (M.char quote) $ M.many $ M.anySingleBut quote
+  return (T.stripEnd $ wordsText ls)
+ where
+  -- consumeDead :: Parser ()
+  -- consumeDead = M.space >> skipEmpty
+  quote = toWord '\''
+
+
+skipEmpty :: Parser ()
+skipEmpty = void (M.many $ M.satisfy (toWord '\0' ==))
+
+
+-- parseEnd :: Parser ()
+-- parseEnd = M.string' "end" >> M.space <* M.eof
+
+parseEquals :: Parser ()
+parseEquals = M.space >> M.char (toWord '=') >> M.space
+
+
+parsePos :: Parser Int
+parsePos = MP.unPos . MP.sourceColumn <$> M.getSourcePos
+
+
+parseBitPix :: Parser BitPix
+parseBitPix = do
+  v <- parseKeywordRecord' "BITPIX" parseValue
+  toBitpix v
+ where
+  toBitpix (Integer 8) = return BPInt8
+  toBitpix (Integer 16) = return BPInt16
+  toBitpix (Integer 32) = return BPInt32
+  toBitpix (Integer 64) = return BPInt64
+  toBitpix (Integer (-32)) = return BPFloat
+  toBitpix (Integer (-64)) = return BPDouble
+  toBitpix _ = fail "Invalid BITPIX header"
+
+
+parseNaxes :: Parser (Axes Column)
+parseNaxes = do
+  n <- parseKeywordRecord' "NAXIS" parseInt
+  Axes <$> mapM parseN [1 .. n]
+ where
+  parseN :: Int -> Parser Int
+  parseN n = parseKeywordRecord' (C8.pack $ "NAXIS" <> show n) parseInt
+
+
+-- | We don't parse simple here, because it isn't required on all HDUs
+parseDimensions :: Parser Dimensions
+parseDimensions = do
+  bp <- parseBitPix
+  Dimensions bp <$> parseNaxes
+
+
+parsePrimaryKeywords :: Parser Dimensions
+parsePrimaryKeywords = do
+  _ <- parseKeywordRecord' "SIMPLE" parseLogic
+  parseDimensions
+
+
+parseImageKeywords :: Parser Dimensions
+parseImageKeywords = do
+  _ <- ignoreComments $ M.string' "XTENSION= 'IMAGE   '"
+  parseDimensions
+
+
+parseBinTableKeywords :: Parser (Dimensions, Int)
+parseBinTableKeywords = do
+  _ <- ignoreComments $ M.string' "XTENSION= 'BINTABLE'"
+  sz <- parseDimensions
+  pc <- parseKeywordRecord' "PCOUNT" parseInt
+  return (sz, pc)
+
+-- parsePrimary :: Parser DataHDU
+-- parsePrimary = do
+--   -- do not consume the headers used for the dimensions
+--   dm <- M.lookAhead parsePrimaryKeywords
+--   hd <- parseHeader
+--   dt <- parseMainData dm
+--   return $ DataHDU hd (dataArray dm dt)
+
+-- parseImage :: Parser DataHDU
+-- parseImage = do
+--   -- do not consume the headers used for the dimensions
+--   dm <- M.lookAhead parseImageKeywords
+--   hd <- parseHeader
+--   dt <- parseMainData dm
+--   return $ DataHDU hd (dataArray dm dt)
+-- parseBinTable :: Parser BinTableHDU
+-- parseBinTable = do
+--   (dm, pc) <- M.lookAhead parseBinTableKeywords
+--   hd <- parseHeader
+--   dt <- parseMainData dm
+--   hp <- parseBinTableHeap
+--   return $ BinTableHDU hd pc hp (dataArray dm dt)
+--  where
+--   parseBinTableHeap = return ""
+
+-- parseMainData :: Dimensions -> Parser ByteString
+-- parseMainData size = do
+--   let len = dataSizeBytes size
+--   M.takeP (Just ("Data Array of " <> show len <> " Bytes")) (fromIntegral len)
+--
+--
+-- parseExtensions :: Parser [Extension]
+-- parseExtensions = do
+--   M.many parseExtension
+--  where
+--   parseExtension :: Parser Extension
+--   parseExtension =
+--     Image <$> parseImage <|> BinTable <$> parseBinTable
diff --git a/src/Telescope/Fits/Encoding/Render.hs b/src/Telescope/Fits/Encoding/Render.hs
--- a/src/Telescope/Fits/Encoding/Render.hs
+++ b/src/Telescope/Fits/Encoding/Render.hs
@@ -1,15 +1,23 @@
+{-# OPTIONS_GHC -fno-warn-orphans #-}
+
 module Telescope.Fits.Encoding.Render where
 
 import Data.ByteString (ByteString)
 import Data.ByteString qualified as BS
 import Data.ByteString.Builder
 import Data.ByteString.Lazy qualified as BL
+import Data.ByteString.Lazy.Char8 qualified as BL8
 import Data.Char (toUpper)
+import Data.List qualified as L
 import Data.String (IsString (..))
 import Data.Text (Text, isPrefixOf, pack, unpack)
 import Data.Text qualified as T
+import Telescope.Data.Axes
+import Telescope.Fits.BitPix
 import Telescope.Fits.Checksum
-import Telescope.Fits.Types
+import Telescope.Fits.DataArray
+import Telescope.Fits.HDU.Block (hduBlockSize)
+import Telescope.Fits.Header
 
 
 renderDataArray :: ByteString -> ByteString
@@ -29,7 +37,7 @@
       , renderKeywordLine "PCOUNT" (Integer 0) Nothing
       , renderKeywordLine "GCOUNT" (Integer 1) Nothing
       , renderDatasum dsum
-      , renderOtherKeywords h
+      , renderRecords $ nonSystemRecords h
       , renderEnd
       ]
 
@@ -42,7 +50,7 @@
       , renderDataKeywords d.bitpix d.axes
       , renderKeywordLine "EXTEND" (Logic T) Nothing
       , renderDatasum dsum
-      , renderOtherKeywords h
+      , renderRecords $ nonSystemRecords h
       , renderEnd
       ]
 
@@ -86,24 +94,37 @@
     BPDouble -> -64
 
 
+renderRecords :: [HeaderRecord] -> BuilderBlock
+renderRecords hrs =
+  mconcat $ fmap renderHeaderRecord hrs
+
+
 -- | 'Header' contains all other keywords. Filter out any that match system keywords so they aren't rendered twice
-renderOtherKeywords :: Header -> BuilderBlock
-renderOtherKeywords (Header ks) =
-  mconcat $ map toLine $ filter (not . isSystemKeyword) ks
+nonSystemRecords :: Header -> [HeaderRecord]
+nonSystemRecords (Header ks) =
+  filter (not . isSystemKeyword) ks
  where
-  toLine (Keyword kr) = renderKeywordLine kr._keyword kr._value kr._comment
-  toLine (Comment c) = pad 80 $ string $ "COMMENT " <> unpack c
-  toLine BlankLine = pad 80 ""
-  isSystemKeyword (Keyword kr) =
-    let k = kr._keyword
-     in k == "BITPIX"
-          || k == "EXTEND"
-          || k == "DATASUM"
-          || k == "CHECKSUM"
-          || "NAXIS" `isPrefixOf` k
-  isSystemKeyword _ = False
+  isSystemKeyword hr =
+    case hr of
+      Keyword kr ->
+        let k = kr.keyword
+         in k == "BITPIX"
+              || k == "EXTEND"
+              || k == "DATASUM"
+              || k == "CHECKSUM"
+              || k == "SIMPLE"
+              || "NAXIS" `isPrefixOf` k
+      _ -> False
 
 
+renderHeaderRecord :: HeaderRecord -> BuilderBlock
+renderHeaderRecord = \case
+  Keyword kr -> renderKeywordLine kr.keyword kr.value kr.comment
+  Comment c -> pad 80 $ string $ "COMMENT " <> unpack c
+  History h -> pad 80 $ string $ "HISTORY " <> unpack h
+  BlankLine -> pad 80 ""
+
+
 -- | Fill out the header or data block to the nearest 2880 bytes
 fillBlock :: (Int -> BuilderBlock) -> BuilderBlock -> BuilderBlock
 fillBlock fill b =
@@ -120,10 +141,10 @@
 renderKeywordLine :: Text -> Value -> Maybe Text -> BuilderBlock
 renderKeywordLine k v mc =
   let kv = renderKeywordValue k v
-   in pad 80 $ addComment kv mc
+   in pad 80 $ insertComment kv mc
  where
-  addComment kv Nothing = kv
-  addComment kv (Just c) =
+  insertComment kv Nothing = kv
+  insertComment kv (Just c) =
     let mx = 80 - kv.length
      in kv <> renderComment mx c
 
@@ -203,3 +224,13 @@
 
 string :: String -> BuilderBlock
 string s = builderBlock (length s) (stringUtf8 s)
+
+
+instance Show Header where
+  show h =
+    unlines $ L.unfoldr chunk $ BL8.unpack $ runRender $ renderRecords h.records
+   where
+    chunk :: String -> Maybe (String, String)
+    chunk "" = Nothing
+    chunk inp =
+      Just $ splitAt 80 inp
diff --git a/src/Telescope/Fits/HDU.hs b/src/Telescope/Fits/HDU.hs
new file mode 100644
--- /dev/null
+++ b/src/Telescope/Fits/HDU.hs
@@ -0,0 +1,72 @@
+module Telescope.Fits.HDU where
+
+import Data.ByteString qualified as BS
+import Data.List qualified as L
+import Telescope.Data.Axes
+import Telescope.Fits.BitPix
+import Telescope.Fits.DataArray
+import Telescope.Fits.Header.Header
+
+
+data Fits = Fits
+  { primaryHDU :: DataHDU
+  , extensions :: [Extension]
+  }
+
+
+instance Show Fits where
+  show f =
+    show f.primaryHDU
+      <> "\n"
+      <> L.intercalate "\n" (fmap show f.extensions)
+
+
+data DataHDU = DataHDU
+  { header :: Header
+  , dataArray :: DataArray
+  }
+
+
+instance Show DataHDU where
+  show p = showHDU "DataHDU" p.header p.dataArray
+
+
+data BinTableHDU = BinTableHDU
+  { header :: Header
+  , pCount :: Int
+  , heap :: BS.ByteString
+  , dataArray :: DataArray
+  }
+
+
+instance Show BinTableHDU where
+  show p = showHDU "BinTableHDU" p.header p.dataArray
+
+
+showHDU :: String -> Header -> DataArray -> String
+showHDU name h d =
+  L.intercalate
+    "\n"
+    [ name
+    , showHeader h
+    , show d
+    ]
+
+
+showHeader :: Header -> String
+showHeader h =
+  "  Header: " <> show (length $ keywords h) <> " records"
+
+
+emptyDataArray :: DataArray
+emptyDataArray = DataArray BPInt8 (Axes []) ""
+
+
+data Extension
+  = Image DataHDU
+  | BinTable BinTableHDU
+
+
+instance Show Extension where
+  show (Image i) = "\nImage: " <> show i
+  show (BinTable b) = "\nBinTable: " <> show b
diff --git a/src/Telescope/Fits/HDU/Block.hs b/src/Telescope/Fits/HDU/Block.hs
new file mode 100644
--- /dev/null
+++ b/src/Telescope/Fits/HDU/Block.hs
@@ -0,0 +1,25 @@
+module Telescope.Fits.HDU.Block where
+
+
+{- | The size of an HDU block is fixed at thirty-six eighty byte words. In
+    other words 2,880 bytes. These blocks are padded with zeros to this
+    boundary.
+-}
+hduBlockSize :: Int
+hduBlockSize = 2880
+
+
+-- \| A single record in the HDU is an eighty byte word.
+{-@ type HDURecordLength = {v:Int | v = 80} @-}
+{-@ hduRecordLength :: HDURecordLength @-}
+hduRecordLength :: Int
+hduRecordLength = 80
+
+
+-- \| The maximum amount of eighty byte records is thirty-six per the
+--        standard.
+--
+{-@ type HDUMaxRecords = {v:Int | v = 36} @-}
+{-@ hduMaxRecords :: HDUMaxRecords @-}
+hduMaxRecords :: Int
+hduMaxRecords = 36
diff --git a/src/Telescope/Fits/Header.hs b/src/Telescope/Fits/Header.hs
--- a/src/Telescope/Fits/Header.hs
+++ b/src/Telescope/Fits/Header.hs
@@ -17,27 +17,18 @@
   , lookupKeyword
   , findKeyword
   , isKeyword
+  , keywords
   , HeaderFor (..)
 
     -- * Re-exports
   , LogicalConstant (..)
-  , getKeywords
   , HeaderRecord (..)
   , KeywordRecord (..)
   ) where
 
-import Data.Fits hiding (isKeyword, lookup)
-import Data.Text (Text)
 import Telescope.Fits.Header.Class
-import Telescope.Fits.Header.Keyword (findKeyword, isKeyword, lookupKeyword)
+import Telescope.Fits.Header.Header
+import Telescope.Fits.Header.Keyword
+import Telescope.Fits.Header.Value
 import Prelude hiding (lookup)
 
-
--- | Construct a keyword HeaderRecord
-keyword :: Text -> Value -> Maybe Text -> HeaderRecord
-keyword k v mc = Keyword $ KeywordRecord k v mc
-
-
--- | Set the comment of a KeywordRecrod
-addComment :: Text -> KeywordRecord -> KeywordRecord
-addComment c kr = kr{_comment = Just c}
diff --git a/src/Telescope/Fits/Header/Class.hs b/src/Telescope/Fits/Header/Class.hs
--- a/src/Telescope/Fits/Header/Class.hs
+++ b/src/Telescope/Fits/Header/Class.hs
@@ -1,15 +1,18 @@
 module Telescope.Fits.Header.Class where
 
-import Data.Fits as Fits hiding (isKeyword)
-import Data.Text (Text, pack)
+import Data.Text (Text, pack, unpack)
 import Data.Text qualified as T
+import Data.Time.Clock (UTCTime)
+import Data.Time.Format.ISO8601 (iso8601ParseM, iso8601Show)
 import Effectful
 import GHC.Generics
 import Telescope.Data.Axes (AxisOrder (..))
 import Telescope.Data.KnownText
 import Telescope.Data.Parser
 import Telescope.Data.WCS (CType (..), CUnit (..), WCSAxis (..), toWCSAxisKey)
-import Telescope.Fits.Header.Keyword (lookupKeyword)
+import Telescope.Fits.Header.Header (Header (..), HeaderRecord (..), lookupKeyword)
+import Telescope.Fits.Header.Keyword
+import Telescope.Fits.Header.Value
 import Text.Casing (fromHumps, toSnake)
 
 
@@ -17,13 +20,16 @@
   toKeywordValue :: a -> Value
 
 
-  -- Can ignore the selector name, modify it, etc
   toKeywordRecord :: Text -> a -> KeywordRecord
   default toKeywordRecord :: Text -> a -> KeywordRecord
   toKeywordRecord key a =
     KeywordRecord key (toKeywordValue a) Nothing
 
 
+class FromKeyword a where
+  parseKeywordValue :: (Parser :> es) => Value -> Eff es a
+
+
 instance ToKeyword Int where
   toKeywordValue = Integer
 instance FromKeyword Int where
@@ -33,10 +39,10 @@
 
 
 instance ToKeyword Float where
-  toKeywordValue = Float
+  toKeywordValue = Float . realToFrac
 instance FromKeyword Float where
   parseKeywordValue = \case
-    Float n -> pure n
+    Float n -> pure $ realToFrac n
     v -> expected "Float" v
 
 
@@ -57,6 +63,17 @@
     v -> expected "Logic" v
 
 
+instance ToKeyword UTCTime where
+  toKeywordValue utc = String $ pack $ iso8601Show utc
+instance FromKeyword UTCTime where
+  parseKeywordValue = \case
+    String t -> do
+      case iso8601ParseM $ unpack t of
+        Nothing -> expected "UTCTime" t
+        Just utc -> pure utc
+    v -> expected "UTCTime" v
+
+
 instance ToKeyword CUnit where
   toKeywordValue (CUnit t) = toKeywordValue t
 instance FromKeyword CUnit where
@@ -73,10 +90,6 @@
     v -> expected "CType" v
 
 
-class FromKeyword a where
-  parseKeywordValue :: (Parser :> es) => Value -> Eff es a
-
-
 class ToHeader a where
   toHeader :: a -> Header
   default toHeader :: (Generic a, GToHeader (Rep a)) => a -> Header
@@ -131,7 +144,7 @@
 
 parseKeyword :: (FromKeyword a, Parser :> es) => Text -> Header -> Eff es a
 parseKeyword k h =
-  case Fits.lookup k h of
+  case lookupKeyword k h of
     Nothing -> parseFail $ "Missing key: " ++ show k
     Just v -> parseAt (Child k) $ parseKeywordValue v
 
diff --git a/src/Telescope/Fits/Header/Header.hs b/src/Telescope/Fits/Header/Header.hs
new file mode 100644
--- /dev/null
+++ b/src/Telescope/Fits/Header/Header.hs
@@ -0,0 +1,55 @@
+module Telescope.Fits.Header.Header where
+
+import Data.List qualified as L
+import Data.Maybe (mapMaybe)
+import Data.Text (Text)
+import Telescope.Fits.Header.Keyword
+import Telescope.Fits.Header.Value
+
+
+{- | The header part of the HDU is vital carrying not only authorship
+    metadata, but also specifying how to make sense of the binary payload
+    that starts 2,880 bytes after the start of the 'HeaderData'.
+-}
+newtype Header = Header {records :: [HeaderRecord]}
+  deriving newtype (Eq, Semigroup, Monoid)
+
+
+
+
+{- | Headers contain lines that are any of the following
+
+ > KEYWORD = VALUE / inline comment
+ > COMMENT full line comment
+ > (blank)
+-}
+data HeaderRecord
+  = Keyword KeywordRecord
+  | Comment Text
+  | History Text
+  | BlankLine
+  deriving (Show, Eq)
+
+
+-- | Manually look up a keyword from the header
+lookupKeyword :: Text -> Header -> Maybe Value
+lookupKeyword k = findKeyword (isKeyword k)
+
+
+findKeyword :: (KeywordRecord -> Bool) -> Header -> Maybe Value
+findKeyword p h = do
+  kr <- L.find p (keywords h)
+  pure kr.value
+
+
+-- | Return all 'KeywordRecord's from the header, filtering out full-line comments and blanks
+keywords :: Header -> [KeywordRecord]
+keywords h = mapMaybe toKeyword h.records
+ where
+  toKeyword (Keyword k) = Just k
+  toKeyword _ = Nothing
+
+
+-- | Construct a keyword HeaderRecord
+keyword :: Text -> Value -> Maybe Text -> HeaderRecord
+keyword k v mc = Keyword $ KeywordRecord k v mc
diff --git a/src/Telescope/Fits/Header/Keyword.hs b/src/Telescope/Fits/Header/Keyword.hs
--- a/src/Telescope/Fits/Header/Keyword.hs
+++ b/src/Telescope/Fits/Header/Keyword.hs
@@ -1,21 +1,25 @@
 module Telescope.Fits.Header.Keyword where
 
-import Data.Fits as Fits hiding (isKeyword)
-import Data.List qualified as L
 import Data.Text (Text)
 import Data.Text qualified as T
-
-
--- | Manually look up a keyword from the header
-lookupKeyword :: Text -> Header -> Maybe Value
-lookupKeyword k = findKeyword (isKeyword k)
+import Telescope.Fits.Header.Value
 
 
-findKeyword :: (KeywordRecord -> Bool) -> Header -> Maybe Value
-findKeyword p h = do
-  kr <- L.find p (getKeywords h)
-  pure kr._value
+{- | A single 80 character header keyword line of the form: KEYWORD = VALUE / comment
+    KEYWORD=VALUE
+-}
+data KeywordRecord = KeywordRecord
+  { keyword :: Text
+  , value :: Value
+  , comment :: Maybe Text
+  }
+  deriving (Show, Eq)
 
 
 isKeyword :: Text -> KeywordRecord -> Bool
 isKeyword k (KeywordRecord k2 _ _) = T.toLower k == T.toLower k2
+
+
+-- | Set the comment of a KeywordRecrod
+addComment :: Text -> KeywordRecord -> KeywordRecord
+addComment c kr = kr{comment = Just c}
diff --git a/src/Telescope/Fits/Header/Value.hs b/src/Telescope/Fits/Header/Value.hs
new file mode 100644
--- /dev/null
+++ b/src/Telescope/Fits/Header/Value.hs
@@ -0,0 +1,17 @@
+module Telescope.Fits.Header.Value where
+
+import Data.Text (Text)
+
+
+-- | `Value` datatype for discriminating valid FITS KEYWORD=VALUE types in an HDU.
+data Value
+  = Integer Int
+  | Float Double
+  | String Text
+  | Logic LogicalConstant
+  deriving (Show, Eq)
+
+
+-- | Direct encoding of a `Bool` for parsing `Value`
+data LogicalConstant = T | F
+  deriving (Show, Eq)
diff --git a/src/Telescope/Fits/Types.hs b/src/Telescope/Fits/Types.hs
deleted file mode 100644
--- a/src/Telescope/Fits/Types.hs
+++ /dev/null
@@ -1,166 +0,0 @@
-{-# LANGUAGE AllowAmbiguousTypes #-}
-
-module Telescope.Fits.Types
-  ( Fits (..)
-  , PrimaryHDU (..)
-  , ImageHDU (..)
-  , BinTableHDU (..)
-  , DataArray (..)
-  , Extension (..)
-  , Axis
-  , Axes (..)
-  , Major (..)
-  , BitPix (..)
-  , bitPixBits
-  , Header (..)
-  , getKeywords
-  , HeaderRecord (..)
-  , KeywordRecord (..)
-  , Value (..)
-  , LogicalConstant (..)
-  , hduBlockSize
-  , emptyDataArray
-  , IsBitPix (..)
-  ) where
-
-import Data.ByteString (ByteString)
-import Data.ByteString qualified as BS
-import Data.Fits (Header (..), HeaderRecord (..), KeywordRecord (..), LogicalConstant (..), Value (..), getKeywords, hduBlockSize)
-import Data.List qualified as L
-import GHC.Int
-import Telescope.Data.Axes
-
-
-data PrimaryHDU = PrimaryHDU
-  { header :: Header
-  , dataArray :: DataArray
-  }
-
-
-instance Show PrimaryHDU where
-  show p = showHDU "PrimaryHDU" p.header p.dataArray
-
-
-data ImageHDU = ImageHDU
-  { header :: Header
-  , dataArray :: DataArray
-  }
-
-
-instance Show ImageHDU where
-  show p = showHDU "ImageHDU" p.header p.dataArray
-
-
-data BinTableHDU = BinTableHDU
-  { header :: Header
-  , pCount :: Int
-  , heap :: ByteString
-  , dataArray :: DataArray
-  }
-
-
-instance Show BinTableHDU where
-  show p = showHDU "BinTableHDU" p.header p.dataArray
-
-
--- | Raw HDU Data. See 'Telescope.Fits.DataArray'
-data DataArray = DataArray
-  { bitpix :: BitPix
-  , axes :: Axes Column
-  , rawData :: BS.ByteString
-  }
-
-
-instance Show DataArray where
-  show d =
-    L.intercalate
-      "\n"
-      [ "  data: " <> show (BS.length d.rawData) <> " bytes"
-      , "  dimensions: "
-      , "    format: " <> L.drop 2 (show d.bitpix)
-      , "    axes: " <> show d.axes.axes
-      ]
-
-
-showHDU :: String -> Header -> DataArray -> String
-showHDU name h d =
-  L.intercalate
-    "\n"
-    [ name
-    , showHeader h
-    , show d
-    ]
-
-
-showHeader :: Header -> String
-showHeader h =
-  "  Header: " <> show (length $ getKeywords h)
-
-
-emptyDataArray :: DataArray
-emptyDataArray = DataArray BPInt8 (Axes []) ""
-
-
--- data BinaryTable = BinaryTable
---   { pCount :: Int
---   , heap :: ByteString
---   }
-
-data Extension
-  = Image ImageHDU
-  | BinTable BinTableHDU
-
-
-instance Show Extension where
-  show (Image i) = show i
-  show (BinTable b) = show b
-
-
-data Fits = Fits
-  { primaryHDU :: PrimaryHDU
-  , extensions :: [Extension]
-  }
-
-
-instance Show Fits where
-  show f =
-    show f.primaryHDU
-      <> "\n"
-      <> L.intercalate "\n" (fmap show f.extensions)
-
-
-data BitPix
-  = BPInt8
-  | BPInt16
-  | BPInt32
-  | BPInt64
-  | BPFloat
-  | BPDouble
-  deriving (Show, Eq)
-
-
-bitPixBits :: BitPix -> Int
-bitPixBits BPInt8 = 8
-bitPixBits BPInt16 = 16
-bitPixBits BPInt32 = 32
-bitPixBits BPInt64 = 64
-bitPixBits BPFloat = 32
-bitPixBits BPDouble = 64
-
-
-class IsBitPix a where
-  bitPix :: BitPix
-
-
-instance IsBitPix Int8 where
-  bitPix = BPInt8
-instance IsBitPix Int16 where
-  bitPix = BPInt16
-instance IsBitPix Int32 where
-  bitPix = BPInt32
-instance IsBitPix Int64 where
-  bitPix = BPInt64
-instance IsBitPix Float where
-  bitPix = BPFloat
-instance IsBitPix Double where
-  bitPix = BPDouble
diff --git a/telescope.cabal b/telescope.cabal
--- a/telescope.cabal
+++ b/telescope.cabal
@@ -5,7 +5,7 @@
 -- see: https://github.com/sol/hpack
 
 name:               telescope
-version:            0.2.0
+version:            0.3.0
 synopsis:           Astronomical Observations (FITS, ASDF, WCS, etc)
 description:        Work with astronomical observations from modern telescopes
 license:            BSD3
@@ -20,6 +20,7 @@
     README.md
     src/checksum.c
     src/checksum.h
+    example/output/hubble1.png
 
 library
     exposed-modules:
@@ -38,22 +39,28 @@
         Telescope.Data.Array
         Telescope.Data.Axes
         Telescope.Data.Binary
+        Telescope.Data.DataCube
         Telescope.Data.KnownText
         Telescope.Data.Parser
         Telescope.Data.WCS
         Telescope.Fits
+        Telescope.Fits.BitPix
         Telescope.Fits.Checksum
         Telescope.Fits.DataArray
         Telescope.Fits.Encoding
+        Telescope.Fits.Encoding.MegaHeader
         Telescope.Fits.Encoding.Render
+        Telescope.Fits.HDU
+        Telescope.Fits.HDU.Block
         Telescope.Fits.Header
         Telescope.Fits.Header.Class
+        Telescope.Fits.Header.Header
         Telescope.Fits.Header.Keyword
-        Telescope.Fits.Types
+        Telescope.Fits.Header.Value
     other-modules:
         Paths_telescope
     build-depends:
-        base >=4.7 && <5,
+        base >=4.18 && <5,
         binary >=0.8.9 && <0.9,
         byte-order >=0.1.3 && <0.2,
         bytestring >=0.11 && <0.13,
@@ -61,9 +68,9 @@
         conduit ==1.3.*,
         effectful >=2.3 && <3,
         exceptions ==0.10.*,
-        fits-parse >=0.4.2 && <0.5,
         libyaml >=0.1.4 && <0.2,
         massiv ==1.0.*,
+        massiv-io ==1.0.*,
         megaparsec ==9.6.*,
         resourcet-effectful >=1.0.1 && <1.1,
         scientific ==0.3.*,
@@ -90,7 +97,7 @@
         src/checksum.c
     default-language: GHC2021
 
-test-suite spec
+test-suite test
     type: exitcode-stdio-1.0
     main-is: Spec.hs
     other-modules:
@@ -105,6 +112,7 @@
         Test.Fits.ChecksumSpec
         Test.Fits.ClassSpec
         Test.Fits.EncodingSpec
+        Test.Fits.MegaHeaderSpec
         Paths_telescope
     hs-source-dirs:
         test/
@@ -126,7 +134,7 @@
     build-tool-depends:
         skeletest:skeletest-preprocessor
     build-depends:
-        base >=4.7 && <5,
+        base >=4.18 && <5,
         binary,
         byte-order >=0.1.3 && <0.2,
         bytestring >=0.11 && <0.13,
@@ -135,9 +143,9 @@
         containers,
         effectful >=2.3 && <3,
         exceptions ==0.10.*,
-        fits-parse >=0.4.2 && <0.5,
         libyaml >=0.1.4 && <0.2,
         massiv ==1.0.*,
+        massiv-io ==1.0.*,
         megaparsec ==9.6.*,
         resourcet-effectful >=1.0.1 && <1.1,
         scientific ==0.3.*,
diff --git a/test/Test/ArraySpec.hs b/test/Test/ArraySpec.hs
--- a/test/Test/ArraySpec.hs
+++ b/test/Test/ArraySpec.hs
@@ -10,7 +10,6 @@
 import System.ByteOrder
 import Telescope.Data.Array
 import Telescope.Data.Axes
-import Telescope.Fits.DataArray
 
 
 spec :: Spec
diff --git a/test/Test/Asdf/DecodeSpec.hs b/test/Test/Asdf/DecodeSpec.hs
--- a/test/Test/Asdf/DecodeSpec.hs
+++ b/test/Test/Asdf/DecodeSpec.hs
@@ -7,7 +7,6 @@
 import Data.Massiv.Array qualified as M
 import Data.Text (Text, unpack)
 import Effectful
-import Effectful.Error.Static
 import GHC.Generics (Generic)
 import GHC.Int (Int64)
 import Skeletest
@@ -266,5 +265,9 @@
     pure $ noCleanup $ RefTreeFix tree
 
 
-parseIO :: Eff '[Parser, Error ParseError, IOE] a -> IO a
-parseIO p = runEff $ runErrorNoCallStackWith @ParseError throwM $ runParser p
+parseIO :: Eff '[Parser, IOE] a -> IO a
+parseIO p = do
+  res <- runEff $ runParser p
+  case res of
+    Left e -> throwM e
+    Right a -> pure a
diff --git a/test/Test/Asdf/EncodeSpec.hs b/test/Test/Asdf/EncodeSpec.hs
--- a/test/Test/Asdf/EncodeSpec.hs
+++ b/test/Test/Asdf/EncodeSpec.hs
@@ -283,3 +283,4 @@
 --     String s -> pure $ PointyName s
 --     InternalRef p -> parsePointer p
 --     other -> expected "PointyName Ref" other
+--
diff --git a/test/Test/Asdf/GWCSSpec.hs b/test/Test/Asdf/GWCSSpec.hs
--- a/test/Test/Asdf/GWCSSpec.hs
+++ b/test/Test/Asdf/GWCSSpec.hs
@@ -25,7 +25,7 @@
 
 toAsdfSpec :: Spec
 toAsdfSpec = do
-  describe "Coorindate Frames" $ do
+  describe "Coordinate Frames" $ do
     it "should auto number axes order" $ do
       o <- expectObject $ toValue frame
       lookup "axes_order" o `shouldBe` Just (toNode @[Int] [0, 1])
@@ -49,7 +49,7 @@
 
 transformSpec :: Spec
 transformSpec = do
-  describe "linear" $ withMarkers ["focus"] $ do
+  describe "linear" $ do
     it "should calculate intercept as 0th value for crval = 0" $ do
       let wcs = WCSAxis{cunit = CUnit "cunit", ctype = CType "Ctype", crpix = 12.0, crval = 0, cdelt = 0.1}
       let Intercept int = wcsIntercept wcs
@@ -60,7 +60,7 @@
       let Intercept int = wcsIntercept wcs
       int `shouldSatisfy` P.approx P.tol (-0.1)
 
-  describe "(<&>) concatenate" $ withMarkers ["focus"] $ do
+  describe "(<&>) concatenate" $ do
     it "should combine two inputs" $ do
       let tx = shift 10 :: Transform (Pix X) (Shift X)
       let ty = shift 20 :: Transform (Pix Y) (Shift Y)
@@ -81,7 +81,7 @@
       tyt `shouldBe` ty.transformation
       tzt `shouldBe` tz.transformation
 
-  describe "concat and compose" $ withMarkers ["focus"] $ do
+  describe "concat and compose" $ do
     it "should compose with higher priority first" $ do
       let tx = shift 10 :: Transform (Pix X) (Shift X)
       let tx2 = scale 10 :: Transform (Shift X) (Scale X)
diff --git a/test/Test/Fits/ChecksumSpec.hs b/test/Test/Fits/ChecksumSpec.hs
--- a/test/Test/Fits/ChecksumSpec.hs
+++ b/test/Test/Fits/ChecksumSpec.hs
@@ -3,10 +3,11 @@
 import Data.ByteString qualified as BS
 import Data.Word
 import Skeletest
+import Telescope.Data.Axes
 import Telescope.Fits
 import Telescope.Fits.Checksum
 import Telescope.Fits.Encoding
-import Telescope.Fits.Types
+import Telescope.Fits.Header
 
 
 spec :: Spec
@@ -44,7 +45,7 @@
 
   describe "encoding" $ do
     it "checksum of empty HDU should be (-0)" $ do
-      let empty = PrimaryHDU (Header []) (DataArray BPInt8 (Axes []) "")
+      let empty = DataHDU (Header []) (DataArray BPInt8 (Axes []) "")
       let out = encodePrimaryHDU empty
       checksum out `shouldBe` Checksum maxBound
 
diff --git a/test/Test/Fits/ClassSpec.hs b/test/Test/Fits/ClassSpec.hs
--- a/test/Test/Fits/ClassSpec.hs
+++ b/test/Test/Fits/ClassSpec.hs
@@ -1,10 +1,8 @@
 module Test.Fits.ClassSpec where
 
 import Control.Monad.Catch (throwM)
-import Data.Fits as Fits hiding (isKeyword)
 import Data.Text (Text)
 import Effectful
-import Effectful.Error.Static
 import GHC.Generics
 import Skeletest
 import Skeletest.Predicate qualified as P
@@ -53,15 +51,15 @@
       toKeywordValue True `shouldBe` Logic T
 
     it "should parseKeywordValue" $ do
-      runPureParser (parseKeywordValue @Int $ Integer 23) `shouldBe` Right 23
-      runPureParser (parseKeywordValue @Text $ String "woot") `shouldBe` Right "woot"
-      runPureParser (parseKeywordValue @Bool $ Logic F) `shouldBe` Right False
+      runPureEff (runParser (parseKeywordValue @Int $ Integer 23)) `shouldBe` Right 23
+      runPureEff (runParser (parseKeywordValue @Text $ String "woot")) `shouldBe` Right "woot"
+      runPureEff (runParser (parseKeywordValue @Bool $ Logic F)) `shouldBe` Right False
 
   describe "ToHeader" $ do
     it "should uppercase and snake keywords" $ do
       let h = toHeader (Test 40 "Alice")
-      length h._records `shouldBe` 2
-      [Keyword (KeywordRecord k1 _ _), Keyword (KeywordRecord k2 _ _)] <- pure h._records
+      length h.records `shouldBe` 2
+      [Keyword (KeywordRecord k1 _ _), Keyword (KeywordRecord k2 _ _)] <- pure h.records
       k1 `shouldBe` "AGE"
       k2 `shouldBe` "FIRST_NAME"
 
@@ -74,12 +72,12 @@
     let h = toHeader $ Parent (HeaderFor $ Test 40 "Alice") (HeaderFor $ Woot 123.456)
     lookupKeyword "age" h `shouldBe` Just (Integer 40)
     lookupKeyword "first_name" h `shouldBe` Just (String "Alice")
-    lookupKeyword "woot" h `shouldSatisfy` P.just (P.con (Float (P.eq 123.456)))
+    lookupKeyword "woot" h `shouldSatisfy` P.just (P.con (Float (P.approx P.tol 123.456)))
 
   describe "FromHeader" $ do
     it "should convert datatype" $ do
       let h = toHeader (Test 40 "Alice")
-      let et = runPureParser $ parseHeader h
+      let et = runPureEff $ runParser $ parseHeader h
       et `shouldSatisfy` P.right (P.con Test{age = P.eq 40})
       et `shouldSatisfy` P.right (P.con Test{firstName = P.eq "Alice"})
 
@@ -96,14 +94,14 @@
 
 wcsSpec :: Spec
 wcsSpec = do
-  describe "axis To/From Header" $ withMarkers ["focus"] $ do
+  describe "axis To/From Header" $ do
     it "should incorporate axis order into keywords" $ do
       let hx = toHeader wcsX
-      Fits.lookup "CRPIX1" hx `shouldBe` Just (Float 1.0)
+      lookupKeyword "CRPIX1" hx `shouldBe` Just (Float 1.0)
 
     it "should incorporate wcsalt into keywords" $ do
       let h = toHeader wcsAY
-      Fits.lookup "CRVAL2A" h `shouldBe` Just (Float 5.0)
+      lookupKeyword "CRVAL2A" h `shouldBe` Just (Float 5.0)
 
     it "should roundtrip" $ do
       let h = toHeader wcsAY
@@ -114,5 +112,9 @@
   wcsAY = WCSAxis (CType "Y") (CUnit "M") 4 5 6 :: WCSAxis 'A Y
 
 
-parseIO :: Eff '[Parser, Error ParseError, IOE] a -> IO a
-parseIO p = runEff $ runErrorNoCallStackWith @ParseError throwM $ runParser p
+parseIO :: Eff '[Parser, IOE] a -> IO a
+parseIO p = do
+  res <- runEff $ runParser p
+  case res of
+    Left e -> throwM e
+    Right a -> pure a
diff --git a/test/Test/Fits/EncodingSpec.hs b/test/Test/Fits/EncodingSpec.hs
--- a/test/Test/Fits/EncodingSpec.hs
+++ b/test/Test/Fits/EncodingSpec.hs
@@ -3,17 +3,23 @@
 module Test.Fits.EncodingSpec where
 
 import Data.ByteString qualified as BS
+import Data.ByteString.Char8 qualified as C8
 import Data.ByteString.Lazy qualified as BL
-import Data.ByteString.Lazy.Char8 qualified as C8
-import Data.Massiv.Array (Ix2)
 import Data.Massiv.Array qualified as M
 import Data.Text (pack)
+import Effectful.State.Static.Local
 import Skeletest
-import Telescope.Fits.DataArray
-import Telescope.Fits.Encoding
+import Telescope.Data.Axes
+import Telescope.Fits
+import Telescope.Fits.BitPix (bitPixBytes)
+import Telescope.Fits.DataArray (Dimensions (..))
+import Telescope.Fits.Encoding (fitsParseThrow, mainData, runMega)
+import Telescope.Fits.Encoding qualified as Enc
+import Telescope.Fits.Encoding.MegaHeader qualified as MH
 import Telescope.Fits.Encoding.Render hiding (justify, pad, spaces)
+import Telescope.Fits.HDU.Block (hduBlockSize)
 import Telescope.Fits.Header
-import Telescope.Fits.Types
+import Test.Fits.MegaHeaderSpec (flattenKeywords)
 
 
 spec :: Spec
@@ -23,13 +29,23 @@
   describe "render data" testRenderData
   describe "encode primary" testEncodePrimary
   describe "round trip" testRoundTrip
+  dataArraySpec
+  fullHDUs
+  simple2x3
+  sampleNSO
 
 
 testDecodeFits :: Spec
 testDecodeFits = do
   describe "simple2x3.fits" $ do
+    it "should parse primary" $ do
+      Simple2x3Raw bs <- getFixture
+      dm <- flip fitsParseThrow bs $ runMega "Primary Header" MH.parsePrimaryKeywords
+      dm.axes `shouldBe` Axes [3, 2]
+
     it "should load metadata" $ do
-      f <- decode =<< BS.readFile "samples/simple2x3.fits"
+      Simple2x3Raw bs <- getFixture
+      f <- decode bs
       let dat = f.primaryHDU.dataArray
           hds = f.primaryHDU.header
       dat.axes `shouldBe` Axes [3, 2]
@@ -37,7 +53,8 @@
       lookupKeyword "CUSTOM" hds `shouldBe` Just (Integer 123456)
 
     it "should load data array" $ do
-      f <- decode =<< BS.readFile "samples/simple2x3.fits"
+      Simple2x3Raw bs <- getFixture
+      f <- decode bs
       arr <- decodeDataArray @Ix2 @Int f.primaryHDU.dataArray
       M.toLists arr `shouldBe` [[0, 1, 2], [3, 4, 5]]
 
@@ -152,7 +169,7 @@
 
 
 run :: BuilderBlock -> String
-run = C8.unpack . runRender
+run = C8.unpack . BL.toStrict . runRender
 
 
 headers :: [String] -> String
@@ -245,12 +262,12 @@
           h2 = f2.primaryHDU.header
       lookupKeyword "NAXIS" h2 `shouldBe` Just (Integer 2)
 
-      let ks = getKeywords hs :: [KeywordRecord]
-          k2 = getKeywords h2 :: [KeywordRecord]
+      let ks = keywords hs :: [KeywordRecord]
+          k2 = keywords h2 :: [KeywordRecord]
 
       length (filter (matchKeyword "BITPIX") k2) `shouldBe` length (filter (matchKeyword "BITPIX") ks)
 
-      length h2._records `shouldBe` length hs._records
+      length h2.records `shouldBe` length hs.records
 
     it "should keep naxes order preserved" $ do
       Simple2x3Fix fs <- getFixture
@@ -259,8 +276,8 @@
 
   describe "image fits" $ do
     it "should roundtrip image extensions" $ do
-      let prim = PrimaryHDU mempty emptyDataArray
-      let img = Image $ ImageHDU mempty emptyDataArray
+      let prim = DataHDU mempty emptyDataArray
+      let img = Image $ DataHDU mempty emptyDataArray
       let out = encode $ Fits prim [img]
       fits2 <- decode out
       length fits2.extensions `shouldBe` 1
@@ -268,6 +285,100 @@
   matchKeyword k (KeywordRecord k2 _ _) = k == k2
 
 
+dataArraySpec :: Spec
+dataArraySpec = describe "data array" $ do
+  it "should parse fake data" $ do
+    let fakeData = "1234" -- Related to NAXIS!
+    d <- fitsParseThrow (mainData (Dimensions BPInt8 (Axes [1, 4]))) fakeData
+    d.rawData `shouldBe` fakeData
+
+  it "should grab correct data array" $ do
+    let fakeData = "1234" -- Related to NAXIS!
+    h <- fitsParseThrow Enc.primary $ flattenKeywords ["SIMPLE = T", "BITPIX = 8", "NAXIS=2", "NAXIS1=2", "NAXIS2=2", "TEST='hi'"] <> "       " <> fakeData
+    h.dataArray.rawData `shouldBe` fakeData
+
+
+fullHDUs :: Spec
+fullHDUs = describe "Full HDUs" $ do
+  it "should include required headers in the keywords" $ do
+    let fakeData = "1234" -- Related to NAXIS!
+    h <- fitsParseThrow Enc.primary $ flattenKeywords ["SIMPLE = T", "BITPIX = 8", "NAXIS=2", "NAXIS1=2", "NAXIS2=2", "TEST='hi'"] <> fakeData
+    length (keywords h.header) `shouldBe` 6
+    lookupKeyword "NAXIS" h.header `shouldBe` Just (Integer 2)
+
+  it "should parse full extension" $ do
+    bt <- fitsParseThrow Enc.binTable $ flattenKeywords ["XTENSION= 'BINTABLE'", "BITPIX = -32", "NAXIS=0", "PCOUNT=0", "GCOUNT=1"]
+    bt.pCount `shouldBe` 0
+    bt.heap `shouldBe` ""
+
+
+simple2x3 :: Spec
+simple2x3 = do
+  describe "simple2x3.fits" $ do
+    it "should parse primary" $ do
+      Simple2x3Raw bs <- getFixture
+      p <- fitsParseThrow Enc.primary bs
+      p.dataArray.axes `shouldBe` Axes [3, 2]
+
+
+sampleNSO :: Spec
+sampleNSO = do
+  describe "NSO Sample FITS Parse" $ do
+    it "should parse primary HDU" $ do
+      DKISTFitsRaw bs <- getFixture
+      p <- fitsParseThrow Enc.primary bs
+      BS.length p.dataArray.rawData `shouldBe` 0
+
+    it "should parse BINTABLE" $ do
+      DKISTFitsRaw bs <- getFixture
+      (bt :: BinTableHDU, rest) <- flip fitsParseThrow bs $ do
+        bt <- Enc.primary >> Enc.binTable
+        rest <- get @BS.ByteString
+        pure (bt, rest)
+      lookupKeyword "INSTRUME" bt.header `shouldBe` Just (String "VISP")
+      lookupKeyword "NAXIS" bt.header `shouldBe` Just (Integer 2)
+
+      let countedHeaderBlocks = 11 -- this was manually counted... until end of all headers
+          payloadLength :: Int = BS.length bt.dataArray.rawData
+          headerLength :: Int = countedHeaderBlocks * hduBlockSize
+          -- sizeOnDisk = 161280
+          heapLength :: Int = bt.pCount
+
+      -- Payload size is as expected
+      payloadLength `shouldBe` 32 * 998 * fromIntegral (bitPixBytes BPInt8)
+      bt.pCount `shouldBe` 95968
+
+      if C8.all (/= '\0') (C8.take 100 (C8.drop (headerLength + payloadLength + heapLength - 100) rest))
+        then pure ()
+        else failTest "The end of the heap has some null data"
+
+      if C8.all (== '\0') (C8.drop (headerLength + payloadLength + heapLength) rest)
+        then pure ()
+        else failTest "The remainder of the file contains real data"
+
+
+newtype FitsDecodedFix = FitsDecodedFix Fits
+instance Fixture FitsDecodedFix where
+  fixtureAction = do
+    FitsEncodedFix enc <- getFixture
+    f <- decode enc
+    pure $ noCleanup $ FitsDecodedFix f
+
+
+newtype Simple2x3Raw = Simple2x3Raw BS.ByteString
+instance Fixture Simple2x3Raw where
+  fixtureAction = do
+    bs <- BS.readFile "samples/simple2x3.fits"
+    pure $ noCleanup $ Simple2x3Raw bs
+
+
+newtype DKISTFitsRaw = DKISTFitsRaw BS.ByteString
+instance Fixture DKISTFitsRaw where
+  fixtureAction = do
+    bs <- BS.readFile "./samples/nso_dkist.fits"
+    pure $ noCleanup $ DKISTFitsRaw bs
+
+
 newtype Simple2x3Fix = Simple2x3Fix Fits
 instance Fixture Simple2x3Fix where
   fixtureAction = do
@@ -285,12 +396,4 @@
     primary =
       let heads = Header [Keyword $ KeywordRecord "WOOT" (Integer 123) Nothing]
           dat = DataArray BPInt8 (Axes [3, 2]) $ BS.pack [0 .. 5]
-       in PrimaryHDU heads dat
-
-
-newtype FitsDecodedFix = FitsDecodedFix Fits
-instance Fixture FitsDecodedFix where
-  fixtureAction = do
-    FitsEncodedFix enc <- getFixture
-    f <- decode enc
-    pure $ noCleanup $ FitsDecodedFix f
+       in DataHDU heads dat
diff --git a/test/Test/Fits/MegaHeaderSpec.hs b/test/Test/Fits/MegaHeaderSpec.hs
new file mode 100644
--- /dev/null
+++ b/test/Test/Fits/MegaHeaderSpec.hs
@@ -0,0 +1,444 @@
+{-# LANGUAGE OverloadedStrings #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+
+module Test.Fits.MegaHeaderSpec where
+
+import Control.Monad (forM_)
+import Data.ByteString (ByteString)
+import Data.ByteString qualified as BS
+import Data.ByteString.Char8 qualified as C8
+import Data.Text (Text)
+import Data.Text qualified as T
+import Data.Text.Encoding qualified as TE
+import Skeletest
+import System.IO
+import Telescope.Data.Axes
+import Telescope.Fits.BitPix
+import Telescope.Fits.DataArray (Dimensions (..))
+import Telescope.Fits.Encoding.MegaHeader
+import Telescope.Fits.Header (Header (..), HeaderRecord (..), KeywordRecord (..), LogicalConstant (..), Value (..), keywords, lookupKeyword)
+import Text.Megaparsec qualified as M
+import Text.Megaparsec.Char qualified as M
+import Prelude hiding (lookup)
+
+
+spec :: Spec
+spec = do
+  basicParsing
+  keywordValueLines
+  comments
+  continue
+  fullRecord
+  fullRecordLine
+  headerMap
+  requiredHeaders
+  sampleNSOHeaders
+  sampleHubbleHeaders
+  sampleJWSTHeaders
+
+
+parse :: Parser a -> ByteString -> IO a
+parse p inp =
+  case M.parse p "Test" inp of
+    Left e -> fail $ showParseError e
+    Right v -> pure v
+
+
+flattenKeywords :: [ByteString] -> ByteString
+flattenKeywords ts = mconcat (map pad ts) <> "END"
+
+
+pad :: ByteString -> ByteString
+pad m =
+  let n = 80 - BS.length m
+   in m <> C8.replicate n ' '
+
+
+basicParsing :: Spec
+basicParsing = describe "Basic Parsing" $ do
+  it "should parse a string" $ do
+    res <- parse parseStringValue "'hello there' "
+    res `shouldBe` "hello there"
+
+  it "should parse a number value" $ do
+    res <- parse parseValue "42 "
+    res `shouldBe` Integer 42
+
+  it "should parse a keyword" $ do
+    res <- parse parseKeyword "WS_TEMP ="
+    res `shouldBe` "WS_TEMP"
+
+  it "should handle keyword symbols" $ do
+    res <- parse parseKeyword "OBSGEO-X=   -5466045.256954942 / [m]"
+    res `shouldBe` "OBSGEO-X"
+
+
+keywordValueLines :: Spec
+keywordValueLines = describe "parse keyword=value" $ do
+  it "should parse an integer" $ do
+    res <- parse parseKeywordValue "KEY=42 "
+    res `shouldBe` ("KEY", Integer 42)
+
+  it "should parse a string" $ do
+    res <- parse parseKeywordValue "KEY='value'"
+    res `shouldBe` ("KEY", String "value")
+
+  it "should absorb spaces" $ do
+    res <- parse parseKeywordValue "KEY   = 'value'   "
+    res `shouldBe` ("KEY", String "value")
+
+  it "should parse a float" $ do
+    res <- parse parseKeywordValue "KEY =   44.88 "
+    res `shouldBe` ("KEY", Float 44.88)
+
+  it "should parse a negative number" $ do
+    res <- parse parseKeywordValue "KEY = -44.88"
+    res `shouldBe` ("KEY", Float (-44.88))
+
+  it "should parse a logical constant" $ do
+    res <- parse parseKeywordValue "KEYT=     T "
+    res `shouldBe` ("KEYT", Logic T)
+
+    res2 <- parse parseKeywordValue "KEYF=     F "
+    res2 `shouldBe` ("KEYF", Logic F)
+
+  it "should ignore comments" $ do
+    res <- parse parseKeywordValue "SIMPLE  =                    T / conforms to FITS standard"
+    res `shouldBe` ("SIMPLE", Logic T)
+
+  it "should strip trailing spaces from strings" $ do
+    res <- parse parseKeywordValue "INSTRUME= 'VISP    '"
+    res `shouldBe` ("INSTRUME", String "VISP")
+
+
+fullRecord :: Spec
+fullRecord = describe "parseKeywordRecord" $ do
+  it "should parse an 80 character record" $ do
+    res <- parse parseKeywordRecord (flattenKeywords ["KEYWORD = 12345"])
+    res `shouldBe` KeywordRecord "KEYWORD" (Integer 12345) Nothing
+
+  it "should parse an a record and comment" $ do
+    res <- parse parseKeywordRecord (flattenKeywords ["KEYWORD = 12345 / this is a comment"])
+    res `shouldBe` KeywordRecord "KEYWORD" (Integer 12345) (Just "this is a comment")
+
+  it "should parse a record, comment, followed by next keyword" $ do
+    res <- parse parseKeywordRecord $ flattenKeywords ["SIMPLE  =                    T / conforms to FITS standard"]
+    res `shouldBe` KeywordRecord "SIMPLE" (Logic T) (Just "conforms to FITS standard")
+
+  it "should handle keyword symbols" $ do
+    res <- parse parseKeywordRecord $ flattenKeywords ["OBSGEO-X=   -5466045.256954942 / [m]"]
+    res `shouldBe` KeywordRecord "OBSGEO-X" (Float (-5466045.256954942)) (Just "[m]")
+
+  it "should handle extension" $ do
+    res <- parse parseKeywordRecord $ flattenKeywords ["XTENSION= 'IMAGE   '"]
+    res `shouldBe` KeywordRecord "XTENSION" (String "IMAGE") Nothing
+
+  it "should not consume blank lines" $ do
+    let inp = flattenKeywords ["SIMPLE  =                    T", " "]
+    res <- flip parse inp $ do
+      kv <- parseKeywordRecord
+      _ <- parseLineBlank
+      pure kv
+    res `shouldBe` KeywordRecord "SIMPLE" (Logic T) Nothing
+
+
+fullRecordLine :: Spec
+fullRecordLine = describe "parseRecordLine" $ do
+  it "should parse a normal line" $ do
+    res <- parse parseRecordLine "NAXIS1  =                  100 / [pix]                                          END"
+    res `shouldBe` Keyword (KeywordRecord "NAXIS1" (Integer 100) (Just "[pix]"))
+
+  it "should parse a comment line" $ do
+    res <- parse parseRecordLine "COMMENT ------------------------------ Telescope -------------------------------END"
+    res `shouldBe` Comment "------------------------------ Telescope -------------------------------"
+
+  it "should parse a blank line" $ do
+    res <- parse parseRecordLine $ flattenKeywords [" "]
+    res `shouldBe` BlankLine
+
+
+comments :: Spec
+comments = do
+  describe "Full-line comments" $ withMarkers ["focus"] $ do
+    it "should parse full-line comments" $ do
+      res <- parse parseLineComment $ flattenKeywords ["COMMENT --------------------------- VISP Instrument ----------------------------"]
+      res `shouldBe` "--------------------------- VISP Instrument ----------------------------"
+
+    it "should parse comments with text" $ do
+      res <- parse parseLineComment $ flattenKeywords ["COMMENT  Keys describing the pointing and operation of the telescope. Including "]
+      res `shouldBe` " Keys describing the pointing and operation of the telescope. Including "
+
+    it "should parse blank comments" $ do
+      res <- parse parseLineComment $ flattenKeywords ["COMMENT                                                                         "]
+      res `shouldBe` "                                                                        "
+
+    it "should parse silly JWST blank-line comments" $ do
+      res <- parse parseLineCommentSpaces $ flattenKeywords ["        Association information                                                "]
+      res `shouldBe` "Association information"
+
+    it "should parse silly JWST blank-line comments as a header" $ do
+      let hs =
+            flattenKeywords
+              [ "                                                                                "
+              , "        Association information                                                 "
+              , "                                                                                "
+              ]
+
+      h <- parse parseHeader hs
+      length h.records `shouldBe` 3
+
+    it "fails space-line comments on blank lines" $ do
+      res <- parse parseRecordLine (flattenKeywords ["                                                                                "])
+      res `shouldBe` BlankLine
+
+  describe "inline comments" $ do
+    it "should parse comment" $ do
+      res <- parse (parseInlineComment 0) $ " / Telescope" <> C8.replicate 68 ' '
+      res `shouldBe` "Telescope"
+
+  describe "parse to line end" $ do
+    it "should parse comment line end" $ do
+      res <- parse (parseLineEnd 0) $ " / Telescope" <> C8.replicate 68 ' '
+      res `shouldBe` Just "Telescope"
+
+    it "should ignore blanks" $ do
+      res <- parse (parseLineEnd 0) $ C8.replicate 80 ' '
+      res `shouldBe` Nothing
+
+    it "should parse comment after blanks" $ do
+      res <- parse (parseLineEnd 0) $ "          / comment " <> C8.replicate 60 ' '
+      res `shouldBe` Just "comment"
+
+  describe "withComments" $ do
+    it "should parse a number, ignoring" $ do
+      res <- parse (withComments parseValue) $ "12345 / Woot" <> C8.replicate 68 ' '
+      res `shouldBe` (Integer 12345, Just "Woot")
+
+    it "should parse no comment" $ do
+      res <- parse (withComments parseValue) $ "12345       " <> C8.replicate 68 ' '
+      res `shouldBe` (Integer 12345, Nothing)
+
+    it "should ignore comments on bintable" $ do
+      -- parse (withComments parseValue) $ "12345       " <> C8.replicate 68 ' '
+      let inp = flattenKeywords ["XTENSION= 'BINTABLE'           / binary table extension"]
+      (_, mc) <- flip parse inp $ do
+        withComments $ M.string' "XTENSION= 'BINTABLE'"
+      mc `shouldBe` Just "binary table extension"
+
+
+continue :: Spec
+continue = describe "Continue Keyword" $ withMarkers ["continue"] $ do
+  it "parses a normal string" $ do
+    let input = "'normal string  '     END"
+    t <- parse parseStringContinue input
+    t `shouldBe` "normal string"
+
+  it "parses a CONTINUE with space after" $ do
+    let input = "'has some sp&'CONTINUE  'ace'                                                                  "
+    t <- parse parseStringContinue input
+    t `shouldBe` "has some space"
+
+  it "parses a CONTINUE with comment" $ do
+    let input = "'has a comm'CONTINUE  'ent'    / hello world                                                 "
+    t <- parse parseStringContinue input
+    t `shouldBe` "has a comment"
+
+  it "parses a short CONTINUE" $ do
+    let input = "'stops before the e'  CONTINUE  'nd'                                                                  "
+    t <- parse parseStringContinue input
+    t `shouldBe` "stops before the end"
+
+  it "combines continue into previous keyword" $ do
+    let hs =
+          [ "CAL_URL = 'https://docs.dkist.nso.edu/projects/visp/en/v2.0.1/l0_to_l1_visp.ht&'"
+          , "CONTINUE  'ml'                                                                  "
+          ]
+
+    h <- parse parseHeader $ flattenKeywords hs
+    lookupKeyword "CAL_URL" h `shouldBe` Just (String "https://docs.dkist.nso.edu/projects/visp/en/v2.0.1/l0_to_l1_visp.html")
+
+  it "parses funny hubble full-line comments" $ do
+    -- Hubble has a comment KEYWORD at the end, which doesn't go fully to the end of the line
+    let input =
+          flattenKeywords
+            [ "COMMENT = 'If you have used HLA files for your research, please include the &'  "
+            , "CONTINUE  'following acknowledgment:&'                                          "
+            ]
+
+    h <- parse parseHeader input
+    h.records `shouldBe` [Keyword $ KeywordRecord "COMMENT" (String "If you have used HLA files for your research, please include the following acknowledgment:&") Nothing]
+
+
+headerMap :: Spec
+headerMap = describe "full header" $ do
+  it "should parse single header" $ do
+    h <- parse parseHeader $ flattenKeywords ["KEY1='value'"]
+    length (keywords h) `shouldBe` 1
+    lookupKeyword "KEY1" h `shouldBe` Just (String "value")
+
+  it "should parse multiple headers " $ do
+    res <- parse parseHeader $ flattenKeywords ["KEY1='value'", "KEY2=  23"]
+    length (keywords res) `shouldBe` 2
+    lookupKeyword "KEY2" res `shouldBe` Just (Integer 23)
+
+  it "should ignore comments" $ do
+    res <- parse parseHeader $ flattenKeywords ["KEY1='value' / this is a comment"]
+    length (keywords res) `shouldBe` 1
+    lookupKeyword "KEY1" res `shouldBe` Just (String "value")
+
+  it "should handle xtension" $ do
+    res <- parse parseHeader $ flattenKeywords ["XTENSION= 'IMAGE   '"]
+    length (keywords res) `shouldBe` 1
+    lookupKeyword "XTENSION" res `shouldBe` Just (String "IMAGE")
+
+  it "should parse blank line before keyword" $ do
+    res <- parse parseHeader $ flattenKeywords [" ", "KEY2    = 22"]
+    res.records `shouldBe` [BlankLine, Keyword (KeywordRecord "KEY2" (Integer 22) Nothing)]
+
+  it "should parse line after keyword" $ do
+    res <- parse parseHeader $ flattenKeywords ["KEY1    = 11", " ", "KEY2    = 22"]
+    res.records `shouldBe` [Keyword (KeywordRecord "KEY1" (Integer 11) Nothing), BlankLine, Keyword (KeywordRecord "KEY2" (Integer 22) Nothing)]
+
+
+requiredHeaders :: Spec
+requiredHeaders =
+  describe "required headers" $ do
+    it "should parse bitpix" $ do
+      res <- parse parseBitPix $ flattenKeywords ["BITPIX = 16"]
+      res `shouldBe` BPInt16
+
+    it "should parse NAXES" $ do
+      res <- parse parseNaxes $ flattenKeywords ["NAXIS   =3", "NAXIS1  =1", "NAXIS2  =2", "NAXIS3  =3"]
+      res `shouldBe` Axes [1, 2, 3]
+
+    it "should parse size" $ do
+      res <- parse parseDimensions $ flattenKeywords ["BITPIX = -32", "NAXIS=2", "NAXIS1=10", "NAXIS2=20"]
+      res.bitpix `shouldBe` BPFloat
+      res.axes `shouldBe` Axes [10, 20]
+
+    it "should parse data headers with data" $ do
+      let fakeData = "1234" -- Related to NAXIS!
+      h <- parse parseHeader $ flattenKeywords ["SIMPLE = T", "BITPIX = 8", "NAXIS=2", "NAXIS1=2", "NAXIS2=2", "TEST='hi'"] <> fakeData
+      length (keywords h) `shouldBe` 6
+      lookupKeyword "NAXIS" h `shouldBe` Just (Integer 2)
+
+
+sampleNSOHeaders :: Spec
+sampleNSOHeaders = do
+  describe "NSO Stripped Headers" $ do
+    it "should parse comment block" $ do
+      let hs =
+            [ "DATASUM = '550335088'          / data unit checksum updated 2023-04-22T04:10:59 "
+            , "   "
+            , "COMMENT ------------------------------ Telescope -------------------------------"
+            , "COMMENT  Keys describing the pointing and operation of the telescope. Including "
+            , "COMMENT     the FITS WCS keys describing the world coordinates of the array.    "
+            ]
+      h <- parse parseHeader $ flattenKeywords hs
+      length (keywords h) `shouldBe` 1
+
+    describe "sample header file" $ do
+      it "should parse all keywords individually" $ do
+        SampleNSO ts <- getFixture
+        forM_ (zip [1 :: Int ..] ts) $ \(_, t) -> do
+          _ <- parse parseRecordLine $ flattenKeywords [TE.encodeUtf8 t]
+          pure ()
+
+      it "should parse bin table keywords" $ do
+        DKISTHeaders bs <- getFixture
+        (sz, pc) <- parse parseBinTableKeywords (mconcat $ C8.lines bs)
+        sz.axes `shouldBe` Axes [32, 998]
+        pc `shouldBe` 95968
+
+
+sampleHubbleHeaders :: Spec
+sampleHubbleHeaders = do
+  describe "Hubble Headers" $ do
+    it "should parse weird comment lines" $ do
+      let hs =
+            [ "SIMPLE  =                    T / conforms to FITS standard                      "
+            , "BITPIX  =                    8 / array data type                                "
+            , "NAXIS   =                    0 / number of array dimensions                     "
+            , "EXTEND  =                    T                                                  "
+            , "DATE    = '2009-11-10'         / date this file was written (yyyy-mm-dd)        "
+            , "FILETYPE= 'SCI      '          / type of data found in data file                "
+            , "                                                                                "
+            , "TELESCOP= 'HST'                / telescope used to acquire data                 "
+            , "INSTRUME= 'WFPC2 '             / identifier for instrument used to acquire data "
+            , "EQUINOX =               2000.0 / equinox of celestial coord. system             "
+            , "                                                                                "
+            , "              / WFPC-II DATA DESCRIPTOR KEYWORDS                                "
+            , "                                                                                "
+            , "ROOTNAME= 'hst_8599_53_wfpc2_f814w_wf' / rootname of the observation set        "
+            ]
+
+      h <- parse parseHeader $ flattenKeywords hs
+      length (keywords h) `shouldBe` 10
+
+    it "unknown error reproduce" $ do
+      let hs =
+            [ "EQUINOX =              '2000.0'                                                 "
+            , "                                                                                "
+            , "              / WFPC-II DATA DESCRIPTOR KEYWORDS                                "
+            , "ROOTNAME= 'hst_8599_53_wfpc2_f814w_wf' / rootname of the observation set        "
+            ]
+
+      h <- parse parseHeader $ flattenKeywords hs
+      length (keywords h) `shouldBe` 2
+
+    it "should parse hubble headers" $ do
+      let hs1 =
+            flattenKeywords
+              [ "SIMPLE  =                    T / conforms to FITS standard                      "
+              , "BITPIX  =                    8 / array data type                                "
+              , "NAXIS   =                    0 / number of array dimensions                     "
+              , "EXTEND  =                    T                                                  "
+              , "DATE    = '2009-11-10'         / date this file was written (yyyy-mm-dd)        "
+              ]
+      HubbleHeaders bs <- getFixture
+      let input = mconcat $ C8.lines bs
+      BS.take 400 input <> "END" `shouldBe` hs1
+
+      h <- parse parseHeader $ BS.take (601 * 80) input <> "END"
+      lookupKeyword "DATE-OBS" h `shouldBe` Just (String "2001-04-07")
+      lookupKeyword "FILTROT" h `shouldBe` Just (Float 0)
+
+
+sampleJWSTHeaders :: Spec
+sampleJWSTHeaders = do
+  describe "JWST Headers" $ withMarkers ["focus"] $ do
+    it "should parse jwst headers" $ do
+      JWSTHeaders bs <- getFixture
+      let input = mconcat $ C8.lines bs
+      h <- parse parseHeader input
+      length (keywords h) `shouldBe` 214
+
+
+newtype DKISTHeaders = DKISTHeaders BS.ByteString
+instance Fixture DKISTHeaders where
+  fixtureAction = do
+    bs <- BS.readFile "./samples/nso_dkist_headers.txt"
+    pure $ noCleanup $ DKISTHeaders bs
+
+
+newtype HubbleHeaders = HubbleHeaders BS.ByteString
+instance Fixture HubbleHeaders where
+  fixtureAction = do
+    bs <- BS.readFile "./samples/hubble_headers.txt"
+    pure $ noCleanup $ HubbleHeaders bs
+
+
+newtype JWSTHeaders = JWSTHeaders BS.ByteString
+instance Fixture JWSTHeaders where
+  fixtureAction = do
+    bs <- BS.readFile "./samples/jwst_headers.txt"
+    pure $ noCleanup $ JWSTHeaders bs
+
+
+newtype SampleNSO = SampleNSO [Text]
+instance Fixture SampleNSO where
+  fixtureAction = do
+    DKISTHeaders bs <- getFixture
+    pure $ noCleanup $ SampleNSO $ filter (not . ignore) $ T.lines $ TE.decodeUtf8 bs
+   where
+    ignore t = T.isPrefixOf "CONTINUE" t || T.isPrefixOf "END" t