```{-# OPTIONS_GHC -Wall #-}

module Atmosphere.Atmosphere( Atmos(..)
, siAtmosphere
, usAtmosphere
, siAtmosphere'
, usAtmosphere'
, atmosphere
) where

import Atmosphere.Constants

data Atmos a = Atmos { atmosTemperature :: a
, atmosPressure :: a
, atmosDensity :: a
, atmosSpeedOfSound :: a
, atmosViscosity :: a
, atmosKinematicViscosity :: a
}

{- |
atmosphere in SI units

Input: altitude in meters

Output: (pressure, density, speed of sound, viscosity, kinematic viscosity)

> pressure            - N/m^2
> density             - kg/m^3
> speed of sound      - m/s
> viscosity           - N-s/m^2
> kinematic viscosity - m^2/s
-}
siAtmosphere :: (Floating a, Ord a) => a -> (a,a,a,a,a,a)
siAtmosphere alt_m = (temp, pressure, density, asound, viscosity, kinematicViscosity)
where
alt_km = 0.001*alt_m
(sigma, delta, theta) = atmosphere alt_km
temp = _TZERO * theta
pressure = _PZERO * delta
density = _RHOZERO * sigma
asound = _AZERO * sqrt theta
viscosity = metricViscosity theta
kinematicViscosity = viscosity/density

{- |
atmosphere in imperial units

Input: altitude in ft

Output: (pressure, density, speed of sound, viscosity, kinematic viscosity)

> pressure            - lb/ft^2
> density             - slugs/ft^3
> speed of sound      - ft/s
> viscosity           - slugs/(ft-s)
> kinematic viscosity - ft^2/s
-}
usAtmosphere :: (Floating a, Ord a) => a -> (a,a,a,a,a,a)
usAtmosphere alt_ft = (temp, pressure, density, asound, viscosity, kinematicViscosity)
where
alt_km = 0.001*_FT2METERS*alt_ft
(sigma, delta, theta) = atmosphere alt_km
temp = _KELVIN2RANKINE*_TZERO*theta
pressure = _PZERO*delta/47.88
density = _RHOZERO*sigma/515.379
asound = (_AZERO/_FT2METERS)*sqrt theta
viscosity=(1.0/_PSF2NSM)*metricViscosity theta
kinematicViscosity = viscosity/density

metricViscosity :: (Floating a, Ord a) => a -> a
metricViscosity theta = _BETAVISC*sqrt(t*t*t)/(t+_SUTH)
where
t = theta * _TZERO

-- | atmosphere in SI units with ADT output
siAtmosphere' :: (Floating a, Ord a) => a -> Atmos a
siAtmosphere' alt = Atmos temp pressure density asound viscosity kinematicViscosity
where
(temp, pressure, density, asound, viscosity, kinematicViscosity) = siAtmosphere alt

-- | atmosphere in imperial units with ADT output
usAtmosphere' :: (Floating a, Ord a) => a -> Atmos a
usAtmosphere' alt = Atmos temp pressure density asound viscosity kinematicViscosity
where
(temp, pressure, density, asound, viscosity, kinematicViscosity) = usAtmosphere alt

{- |
Compute temperature, density, and pressure in standard atmosphere.

Correct to 86 km.  Only approximate thereafter.

Input: alt geometric altitude, km.

Output: (sigma, delta, theta)

> sigma - density/sea-level standard density
> delta - pressure/sea-level standard pressure
> theta - temperature/sea-level std. temperature
-}
atmosphere :: (Floating a, Ord a) => a -> (a,a,a)
atmosphere alt = (sigma, delta, theta)
where
_REARTH = 6369.0             -- radius of the Earth (km)
_GMR = 34.163195

h = alt*_REARTH/(alt+_REARTH) -- geometric to geopotential altitude

(htabI, tbase, ptabI, tgradI) = getI htpgTable
where
getI [htab'] = htab'
getI (htab0:htab1@(h',_,_,_):htabs)
| h > h'    = getI (htab1:htabs)
| otherwise = htab0
getI [] = error "something went wrong"

deltah = h - htabI              -- height above local base
tlocal = tbase + tgradI*deltah  -- local temperature

theta  =  tlocal/_TZERO    -- temperature ratio

delta
| 0.0 == tgradI = ptabI*exp(-_GMR*deltah/tbase)