diff --git a/Project.toml b/Project.toml
index bf58126..4a16ab5 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "Geophysics"
 uuid = "74efdf00-b554-44e1-bd21-abb9281d951c"
 authors = ["Michael Reed"]
-version = "0.3.2"
+version = "0.3.3"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
diff --git a/docs/src/units.md b/docs/src/units.md
index c2cfb3e..fc2fd24 100644
--- a/docs/src/units.md
+++ b/docs/src/units.md
@@ -38,6 +38,10 @@ UnitSystem
 
 Other similar packages include [PhysicalConstants.jl](https://github.com/JuliaPhysics/PhysicalConstants.jl), [MathPhysicalConstants.jl](https://github.com/LaGuer/MathPhysicalConstants.jl), [Unitful.jl](https://github.com/PainterQubits/Unitful.jl.git), [UnitfulSystems.jl](https://github.com/chakravala/UnitfulSystems.jl), [UnitfulUS.jl](https://github.com/PainterQubits/UnitfulUS.jl), [UnitfulAstro.jl](https://github.com/JuliaAstro/UnitfulAstro.jl), [UnitfulAtomic.jl](https://github.com/sostock/UnitfulAtomic.jl), [NaturallyUnitful.jl](https://github.com/MasonProtter/NaturallyUnitful.jl), and [UnitfulMoles.jl](https://github.com/rafaqz/UnitfulMoles.jl).
 
+```@index
+Pages = ["units.md"]
+```
+
 ## Metric SI Unit Systems
 
 In the Systeme International d'Unites (the SI units) the `UnitSystem` constants are derived from the most accurate possible physical measurements and a few exactly defined constants.
diff --git a/src/Geophysics.jl b/src/Geophysics.jl
index 156c20f..c4fbd9c 100644
--- a/src/Geophysics.jl
+++ b/src/Geophysics.jl
@@ -6,15 +6,13 @@ module Geophysics
 using AbstractTensors, LinearAlgebra
 import Base: @pure, show, display
 
-export FluidState, Atmosphere, Weather
+export FluidState, Atmosphere, Weather, Planet
 
 export fluid, radius, gravity, layer
 
 decibel(I₁,I₂=intensity()) = 10log10(I₁/I₂)
 gage(P::Real,P0::Real=pressure()) = P-P0
 
-# thermodynamics
-
 using UnitSystems
 export Metric, English
 import UnitSystems: Constants, molarmass, units
@@ -25,190 +23,383 @@ end
 const Properties = (:molecularmass,:gasconstant,UnitSystems.Constants...)
 const Intrinsic = (:viscosity,:thermalconductivity,:heatvolume,:heatpressure,:heatratio,:prandtl,:sonicspeed,:freedom,:specificenergy,:specificenthalpy)
 
+export flattening, semimajor, semiminor, period, gravitation, frequency, angularfrequency
+export eccentricity, eccentricity2, lineareccentricity, aspectratio, radius, speed, mass
+export latitudegeodetic, latitudegeocentric, latitudeparametric
+export deflectiongeodetic, deflectiongeocentric, deflection
+export centripetal, oblateness, q0, q01, q, q1, dynamicformfactor, secondzonalharmonic
+export gravitycomponents, gravitygeodetic, radiusgeodetic
+
 include("chemistry.jl")
 
+# global geophysics
+
 """
-    FluidState{f}
+    Planet{f,a,t,Gm}
 
-Thermodynamic state of fluid `f` at temperature `T` and pressure `P`.
-Induces derived values `fluid`, `temperature`, `pressure`, `density`, `specificvolume`, `kinematic`, `heatcapacity`, `thermaldiffusivity`, `elasticity`, `specificimpedance`, `intensity`, and values associated with `f::AbstractMole` derivations.
+Celestial `Planet` ellipsoidal reference body with `flattening` of `f`, `semimajor` radius of `a`, sidereal `period` of `t`, and standard `gravitation` parameter `Gm`.
+Further derived parameters include `semiminor`, `frequency`, `angularfrequency`, `eccentricity`, `eccentricity2`, `lineareccentricity`, `aspectratio`, `oblateness`, `dynamicformfactor`, `secondzonalharmonic`, `latitudegeodetic`, `latitudegeocentric`, `latitudeparametric`, `radius`, `radiusgeodetic`, `gravity`, `gravitygeodetic`, and `gravitycomponents`.
 """
-struct FluidState{f,u}
-    T::Float64
-    P::Float64
-end
+struct Planet{f,a,t,Gm} end # WGS 84 spheroid
+const Earth = Planet{1/298.257223563,6378137.0,86164.098903691,3.986004418e14}()
 
-for Gas ∈ (:SutherlandGas,:Mixture,:AtomicGas,:DiatomicGas,:TriatomicGas)
-    @eval (G::$Gas)(T=288.15,P=atm,U=Metric) = FluidState{G,U}(T,P)
-end
+"""
+    flattening(P::Planet) = 1-semiminor(P)/semimajor(P)
 
-@pure units(::FluidState{f,u}) where {f,u} = u
-(U::UnitSystem)(F::FluidState{G,S}) where {G,S} = FluidState{G,U}(temperature(F,U),pressure(F,U))
+Oblate `flattening` parameter of a `Planet` reference ellipsoid (dimensionless).
 
+```Julia
+julia> 1/flattening(Earth)
+$(1/flattening(Earth))
+```
 """
-    fluid(x)
+@pure flattening(::Planet{f}=Earth) where {f} = f
 
-Specification of an `AbstractMole` fluid chemical instance.
 """
-@pure fluid(::FluidState{f}) where f = f
+    semimajor(P::Planet) = semiminor(P)/(1-flattening(P))
+
+Equatorial radius of oblate `Planet` reference ellipsoid (m).
 
+```Julia
+julia> semimajor(Earth) # m
+$(semimajor(Earth))
+```
 """
-    temperature(F::FluidState)
+@pure semimajor(P::Planet{f,a}=Earth,U::US=Metric) where {f,a} = length(a,U)
 
-Absolute thermodynamic temperature scale `T` at `F` (K or °R).
 """
-@pure temperature(F::FluidState) = F.T
-@pure temperature(F::FluidState,U::UnitSystem) = temperature(temperature(F),U,units(F))
+    period(P::Planet) = 2π/angularfrequency(P)
 
+Sidereal rotational `period` of `Planet` (s).
+
+```Julia
+julia> period(Earth) # s
+$(period(Earth))
+```
 """
-    pressure(F::FluidState)
+@pure period(::Planet{f,a,t}=Earth,U::US=Metric) where {f,a,t} = time(t,U)
 
-Absolute force per unit area `P` exerted by `F` (Pa or lb⋅ft⁻²).
 """
-@pure pressure(F::FluidState) = F.P
-@pure pressure(F::FluidState,U::UnitSystem) = pressure(pressure(F),U,units(F))
+    gravitation(P::Planet,U::UnitSystem) = mass(P)*newton(U)
 
-for op ∈ Properties
-    @eval @pure $op(F::FluidState,U::US=units(F)) = $op(fluid(F),U)
-end
-for op ∈ Intrinsic
-    @eval @pure $op(F::FluidState,U::US=units(F)) = $op(temperature(F,U),fluid(F),U)
-end
+Standard `gravitation` parameter for a celestial `Planet` body (m³⋅s⁻²).
 
-@doc """
-    viscosity(F::FluidState) = viscosity(temperature(F),fluid(F))
+```Julia
+julia> gravitation(Earth) # m³⋅s⁻²
+$(gravitation(Earth))
+```
+"""
+@pure gravitation(::Planet{f,a,t,Gm}=Earth) where {f,a,t,Gm} = Gm
+@pure gravitation(P::Planet,U::US) = gravitation(P)/(length(U)*specificenergy(U))
 
-Laminar dynamic vicsosity `μ` at the temperature of `F` (Pa⋅s or lb⋅s¹⋅ft⁻²).
-""" viscosity(F::FluidState)
+"""
+    mass(P::Planet,U::UnitSystem) = gravitation(P,U)/newton(U)
 
-@doc """
-    thermalconductivity(F::FluidState) = conductivity(temperature(F),fluid(F))
+Inertial `mass` of a celestial `Planet` body (kg).
 
-Laminar thermal conductivity `k` at the temperature `F` (W⋅m⁻¹⋅K⁻¹ or lb⋅s⁻¹⋅°R⁻¹).
-""" thermalconductivity(F::FluidState)
+```Julia
+julia> mass(Earth) # kg
+$(mass(Earth))
+```
+"""
+@pure mass(P::Planet,U::US=Metric) = gravitation(P,U)/newton(U)
 
-@doc """
-    heatvolume(F::FluidState) = heatvolume(temperature(F),fluid(F))
+"""
+    frequency(P::Planet) = 1/period(P)
 
-Specific heat `cᵥ` at the temperature of `F` (J⋅kg⁻¹⋅K⁻¹ or ft⋅lb⋅slug⁻¹⋅°R⁻¹).
-""" heatvolume(F::FluidState)
+Cyclic `frequency` of sidereal `Planet` rotation (Hz).
 
-@doc """
-    heatpressure(F::FluidState) = heatpressure(temperature(F),fluid(F))
+```Julia
+julia> frequency(Earth) # Hz
+$(frequency(Earth))
+```
+"""
+@pure frequency(P::Planet,U::US=Metric) = 1/period(P,U)
 
-Specific heat `cₚ` at the temperature of `F` (J⋅kg⁻¹⋅K⁻¹ or ft⋅lb⋅slug⁻¹⋅°R⁻¹).
-""" heatpressure(F::FluidState)
+"""
+    angularfrequency(P::Planet) = 2π/period(P)
 
-@doc """
-    heatratio(F::FluidState) = heatratio(temperature(F),fluid(F))
+Rate of `angularfrequency` of sidereal `Planet` rotation (rad⋅s⁻¹).
 
-Specific heat ratio `γ` at the temperature of `F` (dimensionless).
-""" heatratio(F::FluidState)
+```Julia
+julia> angularfrequency(Earth) # rad⋅s⁻¹
+$(angularfrequency(Earth))
+```
+"""
+@pure angularfrequency(P::Planet,U::US=Metric) = 2π/period(P,U)
 
-@doc """
-    specificenergy(F::FluidState) = specificenergy(temperature(F),fluid(F))
+@pure meanradius(P::Planet,U::US=Metric) = radius_fast(asin(sqrt(1/3)),P,U)
+@pure radius_fast(θ,P::Planet,U::US=Metric) = semimajor(P,U)*(1-flattening(P)*sin(θ)^2)
 
-Specific energy `e` at the temperature of `F` (J⋅kg⁻¹ or ft⋅lb⋅slug⁻¹).
-""" specificenergy(F::FluidState)
+"""
+    semiminor(P::Planet) = semimajor(P)*(1-flattening(P))
 
-@doc """
-    specificenthalpy(F::FluidState) = specificenthalpy(temperature(F),fluid(F))
+Polar radius of oblate `Planet` reference ellipsoid (m).
 
-Specific enthalpy `h` at the temperature of `F` (J⋅kg⁻¹ or ft⋅lb⋅slug⁻¹).
-""" specificenthalpy(F::FluidState)
+```Julia
+julia> semiminor(Earth) # m
+$(semiminor(Earth))
+```
+"""
+@pure semiminor(P::Planet,U::US=Metric) = radius_fast(π/2,P,U)
 
-@doc """
-    freedom(F::FluidState) = freedom(temperature(F),fluid(F))
+"""
+    eccentricity(P::Planet) = sqrt(flattening(P)*(2-flattening(P)))
 
-Degrees of freedom `f` at the temperature of `F` (dimensionless).
-""" freedom(F::FluidState)
+Elliptic `eccentricity` of oblate `Planet` reference ellipsoid (dimensionless).
 
-@doc """
-    prandtl(F::FluidState) = prandtl(temperature(F),fluid(F))
+```Julia
+julia> eccentricity(Earth)
+$(eccentricity(Earth))
+```
+"""
+@pure eccentricity(::Planet{f}) where f = sqrt(f*(2-f))
 
-Prandtl number ratio at the temperature of `F` (dimensionless).
-""" prandtl(F::FluidState)
+"""
+    eccentricity2(P::Planet) = eccentricity(P)/sqrt(1-eccentricity(P))
 
-@doc """
-    sonicspeed(F::FluidState) = sonicspeed(temperature(F),fluid(F))
+Second `eccentricity2` of oblate `Planet reference ellipsoid (dimensionless).
 
-Speed of sound wave disturbance at the temperature of `F` (m⋅s⁻¹ or ft⋅s⁻¹).
-""" sonicspeed(F::FluidState)
+```Julia
+julia> eccentricity2(Earth)
+$(eccentricity2(Earth))
+```
+"""
+@pure eccentricity2(P::Planet{f}) where f = eccentricity(P)/(1-f)
 
 """
-    density(F::FluidState) = pressure(F)/temperature(F)/gasconstant(F)
+    lineareccentricity(P::Planet) = semimajor(P)*eccentricity(P)
 
-Inertial mass per volume `ρ` of at a pressure and temperature (kg⋅m⁻³ or slugs⋅ft⁻³).
+Distance from center to foci of oblate `Planet` reference ellipsoid (m).
 """
-@pure density(F::FluidState,U::US=units(F)) = (pressure(F,U)/temperature(F,U))/gasconstant(F,U)
+@pure lineareccentricity(P::Planet,U::US=Metric) = semimajor(P,U)*eccentricity(P)
 
 """
-    specificvolume(F::FluidState) = 1/density(F)
+    aspectratio(P::Planet) = semiminor(P)/semimajor(P)
+
+Value of `aspectratio` or `semiminor` per `semimajor` of `Planet` (dimensionless).
 
-Specific volume per mass `v` at a pressure and temperature (m³⋅kg⁻¹ or ft³⋅slug⁻¹).
+```Julia
+julia> aspectratio(Earth)
+$(aspectratio(Earth))
+```
 """
-@pure specificvolume(F::FluidState,U::US=units(F)) = inv(density(F,U))
+@pure aspectratio(P::Planet) = semiminor(P)/semimajor(P)
 
 """
-    kinematic(F::FluidState) = viscosity(F)/density(F)
+    latitudegeodetic(θ,P::Planet) = atan(tan(θ)/(1-flattening(P))^2)
 
-Kinematic viscosity ratio `ν` at a pressure and temperature (m²⋅s⁻¹ or ft²⋅s⁻¹).
+Convert geocentric latitude `θ` to geodetic latitude on `Planet` surface (rad).
 """
-@pure kinematic(F::FluidState,U::US=units(F)) = viscosity(F,U)/density(F,U)
+@pure latitudegeodetic(θ,P::Planet=Earth) = atan(tan(θ)/(1-flattening(P))^2)
+@pure deflectiongeodetic(θ,P=Earth) = latitudegeodetic(θ,P)-θ
 
 """
-    heatcapacity(F::FluidState) = heatpressure(F)*density(F)
+    latitudegeocentric(ϕ,P::Planet) = atan(tan(ϕ)*(1-flattening(P))^2)
 
-Specific heat per mass at a pressure and temperature (J⋅m⁻³⋅K⁻¹ or lb⋅ft⁻²⋅°R⁻¹).
+Convert geodetic latitude `ϕ` to geocentric latitude on `Planet` surface (rad).
 """
-@pure heatcapacity(F::FluidState,U::US=units(F)) = heatpressure(F,U)*density(F,U)
+@pure latitudegeocentric(ϕ,P::Planet=Earth) = atan(tan(ϕ)*(1-flattening(P))^2)
+@pure deflectiongeocentric(ϕ,P=Earth) = latitudegeocentric(ϕ,P)-ϕ
 
 """
-    thermaldiffusivity(F::FluidState) = thermalconductivity(F)/heatcapacity(F)
+    latitudeparametric(ϕ,P::Planet) = atan(tan(ϕ)*(1-flattening(P)))
 
-Thermal diffusivity `α` at a pressure and temperature (m²⋅s⁻¹ or ft²⋅s⁻¹).
+Convert geodetic latitude `ϕ` to parametric latitude on surrounding sphere (rad).
 """
-@pure thermaldiffusivity(F::FluidState,U::US=units(F)) = thermalconductivity(F,U)/heatcapacity(F,U)
+@pure latitudeparametric(ϕ,P::Planet=Earth) = atan(tan(ϕ)*(1-flattening(P)))
 
 """
-    elasticity(F::FluidState) = heatratio(F)*pressure(F)
+    deflection(h,ϕ,P::Planet)
 
-Bulk modulus of elasticity `B` at a pressure and temperature (Pa or slug⋅ft⁻¹⋅s⁻²).
+Approximation for `deflection` angle between geodetic and geocentric latitudes (rad).
 """
-@pure elasticity(F::FluidState,U::US=units(F)) = heatratio(F,U)*pressure(F,U)
+@pure function deflection(h,ϕ,P::Planet=Earth,U::US=Metric)
+    f = flattening(P)
+    f*sin(2ϕ)*(1-f/2-h/radiusgeodetic(ϕ,P,U))
+end # geodetic: h, ϕ
 
 """
-    specificimpedance(F::FluidState) = density(F)*sonicspeed(F)
+    latitudegeocentric(h,ϕ,P::Planet) = ϕ-deflection(h,ϕ,P)
 
-Specific acoustic resistance at a pressure and temperature (kg⋅m⁻³⋅s⁻¹ or slug⋅ft⁻³⋅s⁻¹).
+Convert geodetic latitude `ϕ` to geocentric latitude at altitude `h` (rad).
 """
-@pure specificimpedance(F::FluidState,U::US=units(F)) = density(F,U)*sonicspeed(F,U)
+@pure latitudegeocentric(h,ϕ,P::Planet=Earth,U::US=Metric) = ϕ-deflection(h,ϕ,P,U)
 
 """
-    intensity(F::FluidState) = pressure(F)^2/impedance(F)
+    radius(θ,P::Planet) = 1/sqrt((cos(θ)/semimajor(P,U))^2+(sin(θ)/semiminor(P,U))^2)
 
-Instantaneous acoustic intensity `I` at a pressure and temperature (W⋅m⁻² or slug⋅s⁻³).
-""" # irradiance
-@pure intensity(F::FluidState,U::US=units(F)) = pressure(F,U)^2/impedance(F,U)
+Parametrized `radius` of `Planet` in terms of geocentric latitude `θ` coordinate (m).
+"""
+@pure radius(θ,P::Planet,U::US=Metric) = 1/sqrt((cos(θ)/semimajor(P,U))^2+(sin(θ)/semiminor(P,U))^2)
 
-# global geophysics
+"""
+    radiusgeodetic(ϕ,P::Planet)
+
+Approximated `radius` of `Planet` in terms of geodetic latitude `ϕ` coordinate (m).
+"""
+@pure function radiusgeodetic(ϕ,P::Planet=Earth,U::US=Metric)
+    f,a = flattening(P),semimajor(P,U)
+    a*(1-f/2*(1-cos(2ϕ))+5f^2/16*(1-cos(4ϕ)))
+end # approximation, geodetic: λd
+
+@pure _speed(θ,P::Planet,U::US=Metric) = radius(θ,P,U)*angularfrequency(P,U)
+"""
+    speed(θ,P::Planet) = cos(θ)*radius(θ,P)*angularfrequency(P)
+
+Velocity component due to sidereal rotation of `Planet` at latitude `θ` (m⋅s⁻¹).
+
+```Julia
+julia> speed(0,Earth)
+$(speed(0,Earth))
+```
+"""
+@pure speed(θ,P::Planet,U::US=Metric) = _speed(θ,P,U)*cos(θ)
+
+@pure _centripetal(θ,P::Planet=Earth,U::US=Metric) = _speed(θ,P,U)*angularfrequency(P,U)
+"""
+    centripetal(θ,P::Planet) = speed(θ,P)*angularfrequency(P)
+
+Acceleration due to `centripetal` force on `Planet` at geocentric latitude `θ` (m⋅s⁻²).
+"""
+@pure centripetal(θ,P::Planet=Earth,U::US=Metric) = _centripetal(θ,P,U)*cos(θ)
+
+"""
+    gravity(P::Planet) = gravitation(P,U)/semimajor(P,U)^2
 
+Spherical estimate of gravitational acceleration at equator of `Planet` (m⋅s⁻²)
 """
-    Atmosphere{n,U}
+@pure gravity(P::Planet,U::US=Metric) = gravitation(P,U)/semimajor(P,U)^2
 
-Temperature column of planet in unit system `U` having `n` thermal layers.
 """
-struct Atmosphere{n,U}
+    oblateness(θ,P::Planet) = radius(θ,P)*angularfrequency(P)^2/gravity(P)
+
+Constant of `oblateness` at latitude `θ` for gravitating `Planet` ellipse.
+"""
+@pure oblateness(θ,P::Planet=Earth,U::US=Metric) = _centripetal(θ,P,U)/gravity(P,U)
+
+"""
+    oblateness(P::Planet) = radius(π/2,P)*angularfrequency(P)^2/gravity(P)
+
+Constant of `oblateness` at latitude `π/2` for gravitating `Planet` ellipse.
+
+```Julia
+julia> oblateness(Earth) # θ ≡ π/2
+$(oblateness(Earth))
+```
+""" # normal gravity constant
+@pure oblateness(P::Planet=Earth,U::US=Metric) = oblateness(π/2,P,U)
+
+@pure q0(P::Planet) = (e2=eccentricity2(P); ((1+3/e2^2)*atan(e2)-3/e2)/2)
+@pure q01(P::Planet) = (e2=eccentricity2(P); 3((1+1/e2^2)*(1-atan(e2)/e2))-1)
+@pure q(u,P,U) = (E=lineareccentricity(P,U)/u; ((1+3/E^2)*atan(E)-3/E)/2)
+@pure q1(u,P,U) = (E=lineareccentricity(P,U)/u; 3((1+E^-2)*(1-atan(E)/E))-1)
+
+"""
+    dynamicormfactor(::Planet)
+
+Celestial body's second `dynamicformfactor` in nodal precession (dimensionless).
+
+```Julia
+julia> dynamicformfactor(Earth)
+$(dynamicformfactor(Earth))
+```
+"""
+@pure dynamicformfactor(P::Planet{f}=Earth) where f = (1-2oblateness(P)*eccentricity2(P)/15q0(P))*f*(2-f)/3
+
+"""
+    secondzonhalharmonic(P::Planet) = -dynamicformfactor(P)/sqrt(5)
+
+Celestial body's `secondzonalharmonic` in spherical approximation (dimensionless).
+
+```Julia
+julia> secondzonalharmonic(Earth)
+$(secondzonalharmonic(Earth))
+```
+"""
+@pure secondzonalharmonic(P::Planet) = -dynamicformfactor(P)/sqrt(5)
+
+@pure function _gravity(ϕ,P::Planet=Earth,U::US=Metric)
+    β = latitudeparametric(ϕ,P)
+    m,E = oblateness(P),lineareccentricity(P,U)
+    a,b = semimajor(P,U),semiminor(P,U)
+    q = m*eccentricity2(P)*q01(P)/3q0(P)
+    sβ,cβ = sin(β),cos(β)
+    g = gravitation(P,U)/(a*sqrt((a*sβ)^2+(b*cβ)^2))
+    return g*((1+q)*sβ^2 + (1-m-q/2)*cβ^2)
+end # geodetic gravity normal
+
+"""
+    gravity(ϕ,P::Planet)
+
+Calculate total `gravity` at geodetic latitude `ϕ` on `Planet` surface (m⋅s⁻²).
+
+```Julia
+julia> gravity(0,Earth)
+$(gravity(0,Earth))
+
+julia> gravity(π/2,Earth)
+$(gravity(π/2,Earth))
+```
+"""
+@pure function gravity(ϕ,P::Planet{f}=Earth,U::US=Metric) where f
+    sϕ2,ge,gp = sin(ϕ)^2,_gravity(0,P,U),_gravity(π/2,P,U)
+    ge*((1+(aspectratio(P)*(gp/ge)-1)*sϕ2)/sqrt(1-(f*(2-f))*sϕ2))
+end # somigliana(1.0111032235724*π/4)
+
+"""
+    gravitygeodetic(h,ϕ,P::Planet)
+
+Calculate total `gravity` at geodetic latitude `ϕ` and altitude `h` (m⋅s⁻²).
+"""
+@pure function gravitygeodetic(h,ϕ,P::Planet=Earth,U::US=Metric)
+    a,f,m = semimajor(P,U),flattening(P),oblateness(P)
+    gravity(ϕ,P,U)*(1-2*(1+f+m-2f*sin(ϕ)^2)*h/a+3*(h/a)^2)
+end
+
+"""
+    gravitycomponents(h,θ,P::Planet)
+
+Calculate components of `gravity` at geocentric latitude `θ` and altitude `h` (m⋅s⁻²).
+"""
+@pure function gravitycomponents(h,θ,P::Planet=Earth,U::US=Metric)
+    r = radius(θ,P,U)+h
+    J2ar = 3*dynamicformfactor(P)*(semimajor(P,U)/r)^2
+    sθ,cθ = sin(θ),cos(θ)
+    g,ω = gravitation(P,U)/r^2,angularfrequency(P,U)
+    Values(g*J2ar*sθ*cθ+r*ω^2*sθ*cθ,g*(1-J2ar/2*(3sθ^2-1))-r*(ω*cθ)^2)
+end
+
+@pure _gravity(h,θ,P::Planet=Earth,U::US=Metric) = norm(gravitycomponents(h,θ,P,U))
+
+"""
+    gravity(h,θ,P::Planet)
+
+Calculate total `gravity` at geocentric latitude `θ` and altitude `h` (m⋅s⁻²).
+"""
+@pure function gravity(h,θ,P::Planet=Earth,U::US=Metric)
+    gp,gp0 = _gravity(π/2,P,U),_gravity(0,π/2)
+    _gravity(h,θ,P,U)*(1+((gp-gp0)/3gp)*sin(θ)^2)
+end
+
+"""
+    Atmosphere{n,P,U}
+
+Temperature column of `P::Planet` with `U::UnitSystem` having `n` thermal layers.
+"""
+struct Atmosphere{n,P,U}
     a::Values{n,Float64} # lapse rate
     h::Values{n,Float64} # altitude
     m::Values{n,Float64} # molar rate
-    Atmosphere{U}(a::Values{n},h::Values{n},m::Values{n}) where {n,U} = new{n,U}(a,h,m)
+    Atmosphere{P,U}(a::Values{n},h::Values{n},m::Values{n}) where {n,P,U} = new{n,P,U}(a,h,m)
 end
 
-@pure units(::Atmosphere{n,U}) where {n,U} = U
-Atmosphere(a::Values{n},h::Values{n}) where n = Atmosphere{Metric}(a,h)
-Atmosphere{U}(a::Values{n},h::Values{n}) where {n,U} = Atmosphere{U}(a,h,zeros(Values{n}))
-(U::UnitSystem)(A::Atmosphere{n,S}) where {n,S} = Atmosphere{U}(lapserate.(A.a,Ref(U),Ref(S)),length.(A.h,Ref(U),Ref(S)))
+@pure Planet(::Atmosphere{n,P}) where {n,P} = P
+@pure units(::Atmosphere{n,P,U}) where {n,P,U} = U
+Atmosphere(a::Values{n},h::Values{n}) where n = Atmosphere{Earth}(a,h)
+Atmosphere{P}(a::Values{n},h::Values{n}) where {n,P} = Atmosphere{P,Metric}(a,h)
+Atmosphere{P,U}(a::Values{n},h::Values{n}) where {n,P,U} = Atmosphere{P,U}(a,h,zeros(Values{n}))
+(U::UnitSystem)(A::Atmosphere{n,P,S}) where {n,P,S} = Atmosphere{P,U}(lapserate.(A.a,Ref(U),Ref(S)),length.(A.h,Ref(U),Ref(S)))
 
 function display(A::Atmosphere)
     println(typeof(A))
@@ -220,29 +411,46 @@ end
 # local weather models
 
 """
-    Weather{r,g,f,n}
+    Weather{ϕ,f,n,P,U}
 
-Thermodynamic column state of fluid `f` at sea level radius `r` and gravitational acceleration `g` having `n` thermal `Atmosphere` layers.
+Thermodynamic column state of fluid `f` at geodetic latitude `ϕ`, having `n` thermal `Atmosphere` layers (of `P::Planet` with `U::UnitSystem`).
 Induces derived values `fluid`, `temperature`, `pressure`, `density`, `specificvolume`, `kinematic`, `heatcapacity`, `thermaldiffusivity`, `elasticity`, `specificimpedance`, `intensity`, `specificweight`, `geopotential`, and inherited values associated with `f::AbstractMole` derivations.
 """
-struct Weather{r,g,f,n,U} # radius, acceleration
-    A::Atmosphere{n,U} # altitude lapse rate
+struct Weather{ϕ,f,n,P,U} # ϕ ↦ radius, acceleration
+    A::Atmosphere{n,P,U} # altitude lapse rate
     T::Values{n,Float64} # temperature
     p::Values{n,Float64} # pressure
     ρ::Values{n,Float64} # density
-    Weather{r,g,f}(A::Atmosphere{n,U},T::Values{n},p::Values{n},ρ::Values{n}) where {r,g,f,n,U} = new{r,g,f,n,U}(A,T,p,ρ)
+    Tc::Float64
+    ha::Float64
+    Weather{ϕ,f}(A::Atmosphere{n,P,U},T::Values{n},p::Values{n},ρ::Values{n},Tc,ha) where {ϕ,f,n,P,U} = new{ϕ,f,n,P,U}(A,T,p,ρ,Tc,ha)
 end
 
-function Weather{r,g}(A::Atmosphere{n,U},F::FluidState{f}) where {r,g,f,n,U}
+function Weather{ϕ}(A::Atmosphere{n,P,U},F::FluidState{f}) where {ϕ,f,n,P,U}
     T = zeros(Variables{n,Float64})
     p = zeros(Variables{n,Float64})
     ρ = zeros(Variables{n,Float64})
+    μ = zeros(Variables{n,Float64})
+    k = zeros(Variables{n,Float64})
+    c = zeros(Variables{n,Float64})
+    Δμ = zeros(Variables{n,Float64})
+    Δk = zeros(Variables{n,Float64})
+    Δc = zeros(Variables{n,Float64})
     T0,p0 = temperature(F,U),pressure(F,U)
     R,γ = gasconstant(F,U),heatratio(F,U)
+    r,g = radiusgeodetic(ϕ,P,U),gravity(ϕ,P,U)
     T[1] = T0; p[1] = p0; ρ[1] = p0/(R*T0)
+    μ[1],k[1],c[1] = viscosity(F,U),thermalconductivity(F,U),heatvolume(F,U)
+    Tc,ha = 0.0,0.0
     for i ∈ 2:n
         Δh = A.h[i]-A.h[i-1]
-        T[i] = T[i-1]+A.a[i-1]*Δh
+        T[i] = if isinf(A.a[i-1])
+            Tc = (A.a[i]*Δh*T[i]+T[i-1]^2-T[i]^2)/(A.h[i]*Δh+2T[i-1]-2T[i])
+            ha = Δh*(T[i-1]-Tc)/sqrt((T[i-1]-Tc)^2-(T[i]-Tc)^2)
+            Tc+(T[i-1]-Tc)*sqrt(1-(Δh/ha)^2)
+        else
+            T[i-1]+A.a[i-1]*Δh
+        end
         vp,vρ = if iszero(A.a[i-1])
             val = exp((-g/R)*Δh/T[i])
             val,val
@@ -252,18 +460,29 @@ function Weather{r,g}(A::Atmosphere{n,U},F::FluidState{f}) where {r,g,f,n,U}
         end
         p[i] = p[i-1]*vp
         ρ[i] = ρ[i-1]*vρ
+        μ[i] = viscosity(T[i],f,U)
+        k[i] = thermalconductivity(T[i],f,U)
+        c[i] = heatvolume(T[i],f,U)
+        Δμ[i-1] = (μ[i]-μ[i-1])/Δh
+        Δk[i-1] = (k[i]-k[i-1])/Δh
+        Δc[i-1] = (c[i]-c[i-1])/Δh
+        if i == n
+            TT = T[end]+A.a[end]*Δh
+            Δμ[end] = (viscosity(TT,f,U)-μ[i])/Δh
+            Δk[end] = (thermalconductivity(TT,f,U)-k[i])/Δh
+            Δc[end] = (heatvolume(TT,f,U)-c[i])/Δh
+        end
     end
-    Weather{r,g,f}(A,Values(T),Values(p),Values(ρ))
+    Weather{ϕ,f}(A,Values(T),Values(p),Values(ρ),Tc,ha)
 end
 
-(U::UnitSystem)(W::Weather{r,g,f,n,S}) where {r,g,f,n,S} = Weather{r,g,f}(U(W.A),temperature.(W.T,Ref(U),Ref(S)),pressure.(W.p,Ref(U),Ref(S)),density.(W.ρ,Ref(U),Ref(S)))
-(A::Atmosphere)(F::FluidState=Air(288.16,atm,US(A)),r=6.356766e6,g=g₀) = Weather{r,g}(A,F)
-(A::Atmosphere)(T,p=atm,r=6.356766e6,g=g₀) = Weather{r,g}(A,Air(T,p,units(A)))
-(W::Weather)(hG::Real=0) = W(hG,layer(hG,W))
+(U::UnitSystem)(W::Weather{ϕ,f,n,P,S}) where {ϕ,f,n,P,S} = Weather{ϕ,f}(U(W.A),temperature.(W.T,Ref(U),Ref(S)),pressure.(W.p,Ref(U),Ref(S)),density.(W.ρ,Ref(U),Ref(S)))
+(A::Atmosphere)(F::FluidState=Air(288.16,atm,US(A)),r=6.356766e6,g=g₀) = Weather{ϕ}(A,F)
+(A::Atmosphere)(T,p=atm,ϕ=1.0111032235724*π/4) = Weather{ϕ}(A,Air(T,p,units(A)))
+(W::Weather)(h::Real=0) = (hG=altgeopotent(h,W); W(hG,layer(hG,W)))
 function (W::Weather)(hG::Real,i)
-    h = altgeopotent(hG,W)
-    T = _temperature(h,i,W)
-    FluidState{fluid(W),units(W)}(T,pressure(h,T,i,W))
+    T = temperature(hG,i,W)
+    FluidState{fluid(W),units(W)}(T,pressure(hG,T,i,W))
 end
 
 function display(W::Weather)
@@ -284,22 +503,30 @@ end
 @pure @inline layer(h::Real,W::Weather=Standard) = h≤W.A.h[1] ? 1 : (i=findfirst(x->x≥h,W.A.h); isnothing(i) ? length(W.A.h) : i-1)
 @pure @inline layer(h::Real,W::Weather,U::US) = layer(length(h,units(W),U),W)
 
-@pure units(::Weather{r,g,f,n,U}) where {r,g,f,n,U} = U
-@pure fluid(::Weather{r,g,f}=Standard) where {r,g,f} = f
+@pure Planet(::Weather{ϕ,f,n,P}) where {ϕ,f,n,P} = P
+@pure units(::Weather{ϕ,f,n,P,U}) where {ϕ,f,n,P,U} = U
+@pure fluid(::Weather{ϕ,f}=Standard) where {ϕ,f} = f
+
+"""
+    latitude(::Weather)
+
+Geodetic latitude `ϕ` at `Weather` column location (rad).
+"""
+@pure latitude(W::Weather{ϕ}=Standard) where ϕ = ϕ
 
 """
     radius(::Weather)
 
-Sea level radius `r` to planet's center of gravity at `Weather` column location (m or ft).
+Sea level radius `r` to planet's geodetic focus at `Weather` column location (m or ft).
 """
-@pure radius(W::Weather{r}=Standard,U::US=units(W)) where r = length(r,U,units(W))
+@pure radius(W::Weather=Standard,U::US=US(W)) = radiusgeodetic(latitude(W),Planet(W),U)
 
 """
     gravity(::Weather)
 
 Sea level gravitational acceleration `g` at `Weather` column location (m⋅s⁻² or ft⋅s⁻²).
 """
-@pure gravity(W::Weather{r,g}=Standard,U::US=units(W)) where {r,g} = acceleration(g,U,units(W))
+@pure gravity(W::Weather=Standard,U::US=units(W)) = gravity(latitude(W),Planet(W),U)
 
 @pure molecularmass(W::Weather=Standard,U::US=units(W)) = molecularmass(fluid(W),U)
 @pure gasconstant(W::Weather=Standard,U::US=units(W)) = gasconstant(fluid(W),U)
@@ -335,7 +562,13 @@ Geometric altitude `h` conversion from geopotential altitude (m or ft).
 
 Gravitational acceleration `g` at altitude `h` of `Weather` column (m⋅s⁻² or ft⋅s⁻²).
 """
-@pure gravity(h::Real,W::Weather=Standard,U::US=US(W)) = (gravity(W,U)*radius(W,U)^2)/altabs(h,W,U)^2
+@pure function gravity(h::Real,W::Weather=Standard,U::US=US(W))
+    if h ≤ 0.007radius(W)
+        (gravity(W,U)*radius(W,U)^2)/altabs(h,W,U)^2
+    else
+        gravitygeodetic(h,latitude(W),Planet(W),U)
+    end
+end
 @pure gravity(h::Real,W::Weather,U::US,S::US) = gravity(length(h,U,S),W,U)
 
 """
@@ -343,17 +576,27 @@ Gravitational acceleration `g` at altitude `h` of `Weather` column (m⋅s⁻² o
 
 Absolute temperature `T` at geometric altitude `h` of `Weather` location (K or °R).
 """
-@pure temperature(h::Real,i,W::Weather=Standard,U::US=US(W)) = _temperature(altgeopotent(h,W,U),i,W,U)
-@pure function _temperature(hG::Real,i,W::Weather=Standard,U::US=units(W))
+@pure function temperature(hG::Real,i,W::Weather=Standard,U::US=units(W))
     T0,a0,h0 = W[i,U]
-    return T0+a0*(hG-h0)
+    if isinf(a0) # 1976 upper atmosphere
+        Δh = hG-h0
+        if a0 < 0 # 1976 elliptic layer
+            W.Tc+(T0-W.Tc)*sqrt(1-(Δh/W.ha)^2)
+        else # 1976 exponential layer
+            r = radius(Earth1976)
+            ξ = Δh*((r+h0)/(r+hG))
+            1000-(1000-T0)*exp((-0.012/(1000-T0))*ξ)
+        end
+    else # standard
+        return iszero(a0) ? T0 : T0+a0*(hG-h0)
+    end
 end
 
 # k, μ, cᵥ, cₚ, γ, Pr, a, e, h
 
 for op ∈ Intrinsic
-    @eval @pure $op(h::Real,i,W::Weather=Standard,U::US=US(W)) = $op(temperature(h,i,W,U),fluid(W),U)
-    @eval @pure $op(h::Real,i,W::Weather,U::US,S::US) = $op(length(h,U,S),i,W,U)
+    @eval @pure $op(hG::Real,i,W::Weather=Standard,U::US=US(W)) = $op(temperature(hG,i,W,U),fluid(W),U)
+    #@eval @pure $op(h::Real,i,W::Weather,U::US,S::US) = $op(length(h,U,S),i,W,U)
 end
 
 @doc """
@@ -421,10 +664,7 @@ Speed of sound wave disturbance at altitude `h` of `Weather` location (m⋅s⁻
 
 Absolute force per unit area `P`  at altitude `h` of `Weather` column (Pa or slug⋅ft⁻¹⋅s⁻²).
 """
-@pure function pressure(h::Real,i,W::Weather=Standard,U::US=units(W))
-    hG = altgeopotent(h,W,U)
-    pressure(h,_temperature(hG,i,W,U),i,W,U)
-end
+@pure pressure(hG::Real,i,W::Weather=Standard,U::US=units(W)) = pressure(hG,temperature(hG,i,W,U),i,W,U)
 @pure function pressure(hG::Real,T,i,W::Weather=Standard,U::US=units(W))
     g,R = gravity(W,U),gasconstant(W,U)
     T0,a,h0,p = W[i,U]
@@ -440,11 +680,8 @@ end
 
 Inertial mass per volume `ρ` at altitude `h` of `Weather` location (kg⋅m⁻³ or slugs⋅ft⁻³).
 """
-@pure function density(h::Real,i,W::Weather=Standard,U::US=units(W))
-    hG = altgeopotent(h,W,U)
-    density(h,_temperature(hG,i,W,U),i,W,U)
-end
-@pure density(hG::Real,T,i,W::Weather,U::US,S::US) = density(length(h,U,S),temperature(T,U,S),i,W,U)
+@pure density(hG::Real,i,W::Weather=Standard,U::US=units(W)) = density(hG,temperature(hG,i,W,U),i,W,U)
+#@pure density(hG::Real,T,i,W::Weather,U::US,S::US) = density(length(h,U,S),temperature(T,U,S),i,W,U)
 @pure function density(hG::Real,T,i,W::Weather=Standard,U::US=units(W))
     g,R = gravity(W,U),gasconstant(W,U)
     T0,a,h0,_,ρ = W[i,U]
@@ -460,11 +697,8 @@ end
 
 Kinematic viscosity ratio `ν` at altitude `h` of `Weather` location (m²⋅s⁻¹ or ft²⋅s⁻¹).
 """
-@pure function kinematic(h::Real,i,W::Weather=Standard,U::US=units(W))
-    hG = altgeopotent(h,W,U)
-    kinematic(h,_temperature(hG,i,W,U),i,W,U)
-end
-@pure kinematic(hG::Real,T,i,W::Weather,U::US,S::US) = kinematic(length(h,U,S),temperature(T,U,S),i,W,U)
+@pure kinematic(hG::Real,i,W::Weather=Standard,U::US=units(W)) = kinematic(hG,temperature(hG,i,W,U),i,W,U)
+#@pure kinematic(hG::Real,T,i,W::Weather,U::US,S::US) = kinematic(length(h,U,S),temperature(T,U,S),i,W,U)
 @pure kinematic(hG::Real,T,i,W::Weather=Standard,U::US=units(W)) = viscosity(T,fluid(W),U)/density(hG,T,i,W,U)
 
 """
@@ -472,9 +706,8 @@ end
 
 Specific heat per mass at altitude `h` of `Weather` location (J⋅m⁻³⋅K⁻¹ or lb⋅ft⁻²⋅°R⁻¹).
 """
-@pure function heatcapacity(h::Real,i,W::Weather=Standard,U::US=units(W))
-    hG = altgeopotent(h,W,U)
-    T = _temperature(hG,i,W,U)
+@pure function heatcapacity(hG::Real,i,W::Weather=Standard,U::US=units(W))
+    T = temperature(hG,i,W,U)
     heatpressure(T,fluid(W),U)*density(hG,T,i,W,U)
 end
 
@@ -483,9 +716,9 @@ end
 
 Thermal diffusivity `α` at altitude `h` of `Weather` location (m²⋅s⁻¹ or ft²⋅s⁻¹).
 """
-@pure function thermaldiffusivity(h::Real,i,W::Weather=Standard,U::US=units(W))
-    F,hG = fluid(W),altgeopotent(h,W,U)
-    T = _temperature(hG,i,W,U)
+@pure function thermaldiffusivity(hG::Real,i,W::Weather=Standard,U::US=units(W))
+    F = fluid(W)
+    T = temperature(hG,i,W,U)
     thermalconductivity(T,F,U)/heatpressure(T,F,U)/density(hG,T,i,W,U)
 end
 
@@ -494,9 +727,8 @@ end
 
 Bulk modulus of elasticity `B` at altitude `h` of `Weather` location (Pa or slug⋅ft⁻¹⋅s⁻²).
 """
-@pure function elasticity(h,i,W::Weather=Standard,U::US=units(W))
-    hG = altgeopotent(h,W,U)
-    T = _temperature(hG,i,W,U)
+@pure function elasticity(hG,i,W::Weather=Standard,U::US=units(W))
+    T = temperature(hG,i,W,U)
     heatratio(T,fluid(W),U)*pressure(hG,T,i,W,U)
 end
 
@@ -505,9 +737,8 @@ end
 
 Specific acoustic resistance at altitude `h` of `Weather` (kg⋅m⁻³⋅s⁻¹ or slug⋅ft⁻³⋅s⁻¹).
 """
-@pure function specificimpedance(h::Real,i,W::Weather=Standard,U::US=units(W))
-    hG = altgeopotent(h,W,U)
-    T = _temperature(hG,i,W,U)
+@pure function specificimpedance(hG::Real,i,W::Weather=Standard,U::US=units(W))
+    T = temperature(hG,i,W,U)
     density(hG,T,i,W,U)*sonicspeed(T,fluid(W),U)
 end
 
@@ -516,9 +747,8 @@ end
 
 Instantaneous intensity `I` at altitude `h` of `Weather` at location (W⋅m⁻² or slug⋅s⁻³).
 """
-@pure function intensity(h::Real,i,W::Weather=Standard,U::US=units(W))
-    hG = altgeopotent(h,W,U)
-    T = _temperature(hG,i,W,U)
+@pure function intensity(hG::Real,i,W::Weather=Standard,U::US=units(W))
+    T = temperature(hG,i,W,U)
     g,R = gravity(W,U),gasconstant(W,U)
     T0,a,h0,p,ρ = W[i,U]
     (p^2/ρ)*(if iszero(a)
@@ -531,50 +761,57 @@ end
 
 # Grashof number
 
-@pure function grashof(h::Real,i,W::Weather=Standard,U::US=units(W))
-    hG = altgeopotent(h,W,U)
+#=@pure function grashof(hG::Real,i,W::Weather=Standard,U::US=units(W))
     T = _temperature(hG,i,W,U)
     gravity(h,W,U)*(temperature(W,U)-T)*(h^3)/(T*kinematic(hG,T,i,W,U)^2)
-end
+end=#
 
 """
     specificweight(h::Real=0,W::Weather=Standard) = density(h,W)*gravity(h,W)
 
 Specific weight at altitude `h` of `Weather` location (kg⋅m⁻²⋅s⁻² or slugs⋅ft⁻²⋅s⁻²).
 """
-@pure specificweight(h::Real,i,W::Weather=Standard,U::US=US(W)) = density(h,i,W,U)*gravity(h,W,U)
+@pure specificweight(hG::Real,i,W::Weather=Standard,U::US=US(W)) = density(hG,i,W,U)*gravity(hG,W,U)
 
 """
     specificvolume(h::Real=0,W::Weather=Standard) = specificvolume(W(h))
 
 Specific volume per mass `v` at altitude `h` of `Weather` location (m³⋅kg⁻¹, ft³⋅slug⁻¹).
 """
-@pure specificvolume(h::Real,i,W::Weather=Standard,U::US=US(W)) = inv(density(h,i,W,U))
+@pure specificvolume(hG::Real,i,W::Weather=Standard,U::US=US(W)) = inv(density(hG,i,W,U))
 
 """
     geopotential(h::Real=0,W::Weather=Standard) = gravity(h,W)*h
 
 Specifc gravitational potential energy `g` at altitude `h` of `Weather` (m²⋅s⁻², ft²⋅s⁻²).
 """
-@pure geopotential(h::Real,i,W::Weather=Standard,U::US=US(W)) = gravity(h,W,U)*h
+@pure geopotential(h::Real,W::Weather=Standard,U::US=US(W)) = gravity(h,W,U)*h
+@pure geopotential(h::Real,W::Weather,U::US,S::US) = geopotential(length(h,U,S),W,U)
+@pure geopotential(W::Weather=Standard,U::US=Metric) = geopotential(0,W,U)
 
 # common interface
 
-for op ∈ (:temperature,:pressure,:density,:specificweight,:specificvolume,:geopotential,:specificimpedance,:grashof,:thermaldiffusivity,:intensity,:heatcapacity,:kinematic,:elasticity,Intrinsic...)
+for op ∈ (:temperature,:pressure,:density,:specificweight,:specificvolume,:specificimpedance,:thermaldiffusivity,:intensity,:heatcapacity,:kinematic,:elasticity,Intrinsic...) # grashof
     opratio = Symbol(op,:ratio)
     @eval begin
         export $op
+        @pure function $op(h::Real,W::Weather=Standard,U::US=US(W))
+            hG = altgeopotent(h,W,U)
+            $op(hG,layer(hG,W,U),W,U)
+        end
         @pure $op(W::Weather=Standard,U::US=Metric) = $op(0,W,U)
-        @pure $op(h::Real,W::Weather=Standard,U::US=US(W)) = $op(h,layer(h,W,U),W,U)
         @pure $op(h::Real,W::Weather,U::US,S::US) = $op(length(h,U,S),W,U)
-        @pure $op(h,i,W::Weather,U::US,S::US) = $op(length(h,U,S),i,W,U)
+        #@pure $op(h,i,W::Weather,U::US,S::US) = $op(length(h,U,S),i,W,U)
     end
     op ∉ (:heatcapacity,:grashof,:geopotential) && @eval begin
         export $opratio
-        @pure $opratio(h::Real,W::Weather=Standard,U::US=US(W)) = $opratio(h,layer(h,W,U),W,U)
-        @pure $opratio(h::Real,i,W::Weather=Standard,U::US=US(W)) = $op(h,i,W,U)/$op(W,U)
+        @pure function $opratio(h::Real,W::Weather=Standard,U::US=US(W))
+            hG = altgeopotent(h,W,U)
+            $opratio(hG,layer(hG,W,U),W,U)
+        end
+        @pure $opratio(hG::Real,i,W::Weather=Standard,U::US=US(W)) = $op(hG,i,W,U)/$op(W,U)
         @pure $opratio(h::Real,W::Weather,U::US,S::US) = $opratio(length(h,U,S),W,U)
-        @pure $opratio(h::Real,i,W::Weather,U::US,S::US) = $opratio(length(h,U,S),i,W,U)
+        #@pure $opratio(h::Real,i,W::Weather,U::US,S::US) = $opratio(length(h,U,S),i,W,U)
     end
 end
 
diff --git a/src/chemistry.jl b/src/chemistry.jl
index 4b44fff..4663aad 100644
--- a/src/chemistry.jl
+++ b/src/chemistry.jl
@@ -264,3 +264,168 @@ function Base.:+(m::Mixture...)
     Mixture{N,Tuple(vcat(molecules.(m)...))}(Values{N}(vcat(fractions.(m)...)))
 end
 
+# thermodynamic states
+
+"""
+    FluidState{f}
+
+Thermodynamic state of fluid `f` at temperature `T` and pressure `P`.
+Induces derived values `fluid`, `temperature`, `pressure`, `density`, `specificvolume`, `kinematic`, `heatcapacity`, `thermaldiffusivity`, `elasticity`, `specificimpedance`, `intensity`, and values associated with `f::AbstractMole` derivations.
+"""
+struct FluidState{f,u}
+    T::Float64
+    P::Float64
+end
+
+for Gas ∈ (:SutherlandGas,:Mixture,:AtomicGas,:DiatomicGas,:TriatomicGas)
+    @eval (G::$Gas)(T=288.15,P=atm,U=Metric) = FluidState{G,U}(T,P)
+end
+
+@pure units(::FluidState{f,u}) where {f,u} = u
+(U::UnitSystem)(F::FluidState{G,S}) where {G,S} = FluidState{G,U}(temperature(F,U),pressure(F,U))
+
+"""
+    fluid(x)
+
+Specification of an `AbstractMole` fluid chemical instance.
+"""
+@pure fluid(::FluidState{f}) where f = f
+
+"""
+    temperature(F::FluidState)
+
+Absolute thermodynamic temperature scale `T` at `F` (K or °R).
+"""
+@pure temperature(F::FluidState) = F.T
+@pure temperature(F::FluidState,U::UnitSystem) = temperature(temperature(F),U,units(F))
+
+"""
+    pressure(F::FluidState)
+
+Absolute force per unit area `P` exerted by `F` (Pa or lb⋅ft⁻²).
+"""
+@pure pressure(F::FluidState) = F.P
+@pure pressure(F::FluidState,U::UnitSystem) = pressure(pressure(F),U,units(F))
+
+for op ∈ Properties
+    @eval @pure $op(F::FluidState,U::US=units(F)) = $op(fluid(F),U)
+end
+for op ∈ Intrinsic
+    @eval @pure $op(F::FluidState,U::US=units(F)) = $op(temperature(F,U),fluid(F),U)
+end
+
+@doc """
+    viscosity(F::FluidState) = viscosity(temperature(F),fluid(F))
+
+Laminar dynamic vicsosity `μ` at the temperature of `F` (Pa⋅s or lb⋅s¹⋅ft⁻²).
+""" viscosity(F::FluidState)
+
+@doc """
+    thermalconductivity(F::FluidState) = conductivity(temperature(F),fluid(F))
+
+Laminar thermal conductivity `k` at the temperature `F` (W⋅m⁻¹⋅K⁻¹ or lb⋅s⁻¹⋅°R⁻¹).
+""" thermalconductivity(F::FluidState)
+
+@doc """
+    heatvolume(F::FluidState) = heatvolume(temperature(F),fluid(F))
+
+Specific heat `cᵥ` at the temperature of `F` (J⋅kg⁻¹⋅K⁻¹ or ft⋅lb⋅slug⁻¹⋅°R⁻¹).
+""" heatvolume(F::FluidState)
+
+@doc """
+    heatpressure(F::FluidState) = heatpressure(temperature(F),fluid(F))
+
+Specific heat `cₚ` at the temperature of `F` (J⋅kg⁻¹⋅K⁻¹ or ft⋅lb⋅slug⁻¹⋅°R⁻¹).
+""" heatpressure(F::FluidState)
+
+@doc """
+    heatratio(F::FluidState) = heatratio(temperature(F),fluid(F))
+
+Specific heat ratio `γ` at the temperature of `F` (dimensionless).
+""" heatratio(F::FluidState)
+
+@doc """
+    specificenergy(F::FluidState) = specificenergy(temperature(F),fluid(F))
+
+Specific energy `e` at the temperature of `F` (J⋅kg⁻¹ or ft⋅lb⋅slug⁻¹).
+""" specificenergy(F::FluidState)
+
+@doc """
+    specificenthalpy(F::FluidState) = specificenthalpy(temperature(F),fluid(F))
+
+Specific enthalpy `h` at the temperature of `F` (J⋅kg⁻¹ or ft⋅lb⋅slug⁻¹).
+""" specificenthalpy(F::FluidState)
+
+@doc """
+    freedom(F::FluidState) = freedom(temperature(F),fluid(F))
+
+Degrees of freedom `f` at the temperature of `F` (dimensionless).
+""" freedom(F::FluidState)
+
+@doc """
+    prandtl(F::FluidState) = prandtl(temperature(F),fluid(F))
+
+Prandtl number ratio at the temperature of `F` (dimensionless).
+""" prandtl(F::FluidState)
+
+@doc """
+    sonicspeed(F::FluidState) = sonicspeed(temperature(F),fluid(F))
+
+Speed of sound wave disturbance at the temperature of `F` (m⋅s⁻¹ or ft⋅s⁻¹).
+""" sonicspeed(F::FluidState)
+
+"""
+    density(F::FluidState) = pressure(F)/temperature(F)/gasconstant(F)
+
+Inertial mass per volume `ρ` of at a pressure and temperature (kg⋅m⁻³ or slugs⋅ft⁻³).
+"""
+@pure density(F::FluidState,U::US=units(F)) = (pressure(F,U)/temperature(F,U))/gasconstant(F,U)
+
+"""
+    specificvolume(F::FluidState) = 1/density(F)
+
+Specific volume per mass `v` at a pressure and temperature (m³⋅kg⁻¹ or ft³⋅slug⁻¹).
+"""
+@pure specificvolume(F::FluidState,U::US=units(F)) = inv(density(F,U))
+
+"""
+    kinematic(F::FluidState) = viscosity(F)/density(F)
+
+Kinematic viscosity ratio `ν` at a pressure and temperature (m²⋅s⁻¹ or ft²⋅s⁻¹).
+"""
+@pure kinematic(F::FluidState,U::US=units(F)) = viscosity(F,U)/density(F,U)
+
+"""
+    heatcapacity(F::FluidState) = heatpressure(F)*density(F)
+
+Specific heat per mass at a pressure and temperature (J⋅m⁻³⋅K⁻¹ or lb⋅ft⁻²⋅°R⁻¹).
+"""
+@pure heatcapacity(F::FluidState,U::US=units(F)) = heatpressure(F,U)*density(F,U)
+
+"""
+    thermaldiffusivity(F::FluidState) = thermalconductivity(F)/heatcapacity(F)
+
+Thermal diffusivity `α` at a pressure and temperature (m²⋅s⁻¹ or ft²⋅s⁻¹).
+"""
+@pure thermaldiffusivity(F::FluidState,U::US=units(F)) = thermalconductivity(F,U)/heatcapacity(F,U)
+
+"""
+    elasticity(F::FluidState) = heatratio(F)*pressure(F)
+
+Bulk modulus of elasticity `B` at a pressure and temperature (Pa or slug⋅ft⁻¹⋅s⁻²).
+"""
+@pure elasticity(F::FluidState,U::US=units(F)) = heatratio(F,U)*pressure(F,U)
+
+"""
+    specificimpedance(F::FluidState) = density(F)*sonicspeed(F)
+
+Specific acoustic resistance at a pressure and temperature (kg⋅m⁻³⋅s⁻¹ or slug⋅ft⁻³⋅s⁻¹).
+"""
+@pure specificimpedance(F::FluidState,U::US=units(F)) = density(F,U)*sonicspeed(F,U)
+
+"""
+    intensity(F::FluidState) = pressure(F)^2/impedance(F)
+
+Instantaneous acoustic intensity `I` at a pressure and temperature (W⋅m⁻² or slug⋅s⁻³).
+""" # irradiance
+@pure intensity(F::FluidState,U::US=units(F)) = pressure(F,U)^2/impedance(F,U)
diff --git a/src/planets.jl b/src/planets.jl
index ca98efb..ccbdd4e 100644
--- a/src/planets.jl
+++ b/src/planets.jl
@@ -9,6 +9,24 @@ export Earth1959English, Earth1962English, Earth1966English, Earth1976English
 #   This file is part of Geophysics.jl. It is licensed under the AGPL license
 #   Geophysics Copyright (C) 2020 Michael Reed
 
+export Mercury, Venus, Earth, Mars
+export Jupiter, Saturn, Uranus, Neptune
+export Sun, Moon, Pluto, Ceres, Eris
+
+const Sun = Planet{0.00005,696342e3,25.38*24*60^2,1.32712440018e20}()
+const Mercury = Planet{0,2439.7e3,1407.5*60^2,2.2032e13}()
+const Venus = Planet{0,6051.8e3,-243.025*24*60^2,3.24859e14}()
+const Earth = Planet{1/298.257223563,6378137.0,86164.098903691,3.986004418e14}()
+const Moon = Planet{0.0012,1738.1e3,27.321661*24*60^2,4.9048695e12}()
+const Mars = Planet{0.00589,3396.2e3,1.025957*24*60^2,4.282837e13}()
+const Jupiter = Planet{0.06487,71492e3,9.925*60^2,1.26686534e17}()
+const Saturn = Planet{0.09796,60268e3,38018,3.7931187e16}()
+const Uranus = Planet{0.02293,25559e3,-0.71833*24*60^2,5.793939e15}()
+const Neptune = Planet{0.01708,24764e3,16.11*60^2,6.836529e15}()
+const Pluto = Planet{0,1188.3e3,6.38723*24*60^2,8.71e11}()
+const Ceres = Planet{0,469.73e3,9.074170*60^2,6.26325e10}()
+const Eris = Planet{0,1163e3,349.44*60^2,1.108e12}()
+
 export Nitrogen, Oxygen, CarbonDioxide, Methane, Hydrogen
 export Argon, Neon, Helium, Krypton, Xenon
 export N2, O2, CO2, CH4, H2, N₂, O₂, CO₂, CH₄, H₂, Ar, Ne, He, Kr, Xe
@@ -56,44 +74,44 @@ Nitrox else air end
 
 # Atmospheric temperature models
 
-const US22 = Atmosphere{Metric}(Values(-6.5e-3,0e-3),Values(-00e3,11e3))
-const US25 = Atmosphere{Metric}(Values(-6.5e-3,0e-3),Values(-00e3,10.76923e3))
-const US56 = Atmosphere{Metric}(
+const US22 = Atmosphere{Earth,Metric}(Values(-6.5e-3,0e-3),Values(-00e3,11e3))
+const US25 = Atmosphere{Earth,Metric}(Values(-6.5e-3,0e-3),Values(-00e3,10.76923e3))
+const US56 = Atmosphere{Earth,Metric}(
     Values(-6.5e-3,0e-3,3e-3,0e-3,-3.9e-3,0e-3,3.5e-3,10e-3,5.8e-3),
     Values(-00e3,11e3,25e3,47e3,53e3,75e3,90e3,126e3,175e3))
-const US59 = Atmosphere{Metric}(
+const US59 = Atmosphere{Earth,Metric}(
     Values(-6.5e-3,0e-3,3e-3,0e-3,-4.5e-3,0e-3,4e-3,20e-3,10e-3,5e-3,3.5e-3),
     Values(-00e3,11e3,25e3,47e3,53e3,79e3,90e3,105e3,160e3,170e3,200e3))
-#=const US62 = Atmosphere{Metric}(
+#=const US62 = Atmosphere{Earth,Metric}(
     Values(-6.5e-3,0.,1e-3,2.8e-3,0e3,-2e-3,-4e-3,0.,3e-3),
     Values(-00e3,11e3,20e3,32e3,47e3,52e3,61e3,79e3,90e3))=#
-const US62 = Atmosphere{Metric}(
+const US62 = Atmosphere{Earth,Metric}(
     Values(-6.5e-3,0.,1e-3,2.8e-3,0e3,-2e-3,-4e-3,0.,3e-3, 5e-3, 10e-3, 20e-3, 15e-3, 10e-3, 7e-3, 5e-3, 4e-3, 3.3e-3, 2.6e-3, 1.7e-3, 1.1e-3),
     Values(-00e3,11e3,20e3,32e3,47e3,52e3,61e3,79e3,90e3, 100e3, 110e3, 120e3, 150e3, 160e3, 170e3, 190e3, 230e3, 300e3, 400e3, 600e3, 700e3))
-const US66 = Atmosphere{Metric}(
+const US66 = Atmosphere{Earth,Metric}(
     Values(-6.5e-3,0.,1e-3,2.8e-3,0e3,-2e-3,-3.9e-3,0.,3e-3),
     Values(-00e3,11e3,20.1e3,32.2e3,47.3e3,52.4e3,61.6e3,80e3,90e3))
-const US76 = Atmosphere{Metric}(
+const US76 = Atmosphere{Earth,Metric}(
     Values(-6.5e-3,0.,1e-3,2.8e-3,0e3,-2.8e-3,-2e-3,0.,-Inf,12e-3,Inf),
     Values(-00e3,11e3,20e3,32e3,47e3,51e3,71e3,86e3,91e3,110e3,120e3))
-const US22E = Atmosphere{English}(Values(-3.5658e-3,0e-3),Values(-0e3,36.089e3))
-const US25E = Atmosphere{English}(Values(-3.5658e-3,0e-3),Values(-0e3,35.332e3))
-const US56E = Atmosphere{English}(
+const US22E = Atmosphere{Earth,English}(Values(-3.5658e-3,0e-3),Values(-0e3,36.089e3))
+const US25E = Atmosphere{Earth,English}(Values(-3.5658e-3,0e-3),Values(-0e3,35.332e3))
+const US56E = Atmosphere{Earth,English}(
     Values(-3.5658e-3,0.,1.64584e-3,0.,-2.1397e-3,0.,1.92024e-3,5.4864e-3,3.1821e-3),
     Values(-0.,36089.,82021.,154199.,173885.,246063.,295276.,413386.,574147.))
-const US59E = Atmosphere{English}(
+const US59E = Atmosphere{Earth,English}(
     Values(-3.5658e-3,0.,1.64584e-3,0.,-2.46876e-3,0.,2.19456e-3,10.9728e-3,5.4864e-3,2.7432e-3,1.92024e-3),
     Values(-0.,36089.,82021.,154199.,173885.,259176.,295276.,344488.,524934.,557743.,656168.))
-const US62E = Atmosphere{English}(
+const US62E = Atmosphere{Earth,English}(
     Values(-3.5658e-3,0.,0.54864e-3,1.53612e-3,0e3,-1.09728e-3,-2.1946e-3,0.,1.6459e-3),
     Values(-0.,36089.,65617.,104987.,154199.,170604.,200131.,259186.,295276.))
-const US66E = Atmosphere{English}(
+const US66E = Atmosphere{Earth,English}(
     Values(-3.5658e-3,0.,0.54864e-3,1.53612e-3,0e3,-1.09728e-3,-2.1397e-3,0.,1.6459e-3),
     Values(-0.,36089.,65945.,105643.,155184.,171916.,202.1e3,262467.,295276.))
-const US76E = Atmosphere{English}(
+const US76E = Atmosphere{Earth,English}(
     Values(-3.5658e-3,0e3,0.54864e-3,1.53612e-3,0e3,-1.53612e-3,-1.09728e-3),
     Values(-0e3,36.089e3,65.617e3,104.987e3,154.199e3,167.323e3,232.940e3))
-#=const MarsAtmosphere = Atmosphere{MarsAir}(
+#=const MarsAtmosphere = Atmosphere{Mars,MarsAir}(
     Values(-9.88e-4,-2.22e-3),
     Values(-0e3,7e3))=#
 
@@ -102,14 +120,15 @@ const ARDC,ARDCE = US59,US59E # deprecate?
 const layers = Values("Troposphere","Tropopause","Stratosphere","Stratosphere","Stratopause","Mesosphere","Mesosphere","Mesopause")
 
 # Standard atmosphere weather conditions (sea level 45° lat)
+# 1.0111032235724*π/4 # 6.356766e6, g₀ # 2.085553e7, lbm
 
-const Earth1922,Earth1922English = US22(288.16),US22E(518.69,2116.2,2.085553e7,lbm)
-const Earth1925,Earth1925English = US25(288.16),US25E(518.69,2116.2,2.085553e7,lbm)
-const Earth1956,Earth1956English = US56(288.16),US56E(518.69,2116.2,2.085553e7,lbm)
-const Earth1959,Earth1959English = US59(288.16),US59E(518.69,2116.2,2.085553e7,lbm)
-const Earth1962,Earth1962English = US62(288.15),US62E(518.67,2116.2,2.085553e7,lbm)
-const Earth1966,Earth1966English = US66(288.15),US66E(518.67,2116.2,2.085553e7,lbm)
-const Earth1976,Earth1976English = US76(288.15),US76E(518.67,2116.2,2.085553e7,lbm)
+const Earth1922,Earth1922English = US22(288.16),US22E(518.69,2116.2)
+const Earth1925,Earth1925English = US25(288.16),US25E(518.69,2116.2)
+const Earth1956,Earth1956English = US56(288.16),US56E(518.69,2116.2)
+const Earth1959,Earth1959English = US59(288.16),US59E(518.69,2116.2)
+const Earth1962,Earth1962English = US62(288.15),US62E(518.67,2116.2)
+const Earth1966,Earth1966English = US66(288.15),US66E(518.67,2116.2)
+const Earth1976,Earth1976English = US76(288.15),US76E(518.67,2116.2)
 #const MarsWeather = MarsAtmosphere(242.15,699.) # to-do
 const english = haskey(ENV,"GEOUNITS") && ENV["GEOUNITS"] == "english"
 const Standard = if haskey(ENV,"STDATM")