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

module Atmosphere
( Atmos(..)
, siAtmosphere
, usAtmosphere
, atmosphere
, siAltitudeFromPressure
) 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 -> Atmos a
siAtmosphere alt_m =
Atmos { atmosTemperature = temp
, atmosPressure = pressure
, atmosDensity = density
, atmosSpeedOfSound = asound
, atmosViscosity = viscosity
, atmosKinematicViscosity = 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 -> Atmos a
usAtmosphere alt_ft =
Atmos { atmosTemperature = temp
, atmosPressure = pressure
, atmosDensity = density
, atmosSpeedOfSound = asound
, atmosViscosity = viscosity
, atmosKinematicViscosity = 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

{- |
Compute altitude at which the standard atmosphere has a certain pressure.

Input: Pressure, N/m^2

Output: Altitude in meters
-}
siAltitudeFromPressure :: (Floating a, Ord a) => a -> a
siAltitudeFromPressure pressureIn = 1000*alt
where
alt = _REARTH / (_REARTH/h - 1)
deltaIn = pressureIn / _PZERO

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

h
| 0.0 == tgradI = htabI - tbase / _GMR * (log (deltaIn / ptabI))

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

{- |
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
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)