diff --git a/.gitignore b/.gitignore
index 6a1c302..73d7a20 100644
--- a/.gitignore
+++ b/.gitignore
@@ -3,4 +3,5 @@
 docs/build/*
 *.code-workspace
 .VSCodeCounter/*
-*png
\ No newline at end of file
+*png
+*json
\ No newline at end of file
diff --git a/Project.toml b/Project.toml
index 4e46694..48bac4a 100644
--- a/Project.toml
+++ b/Project.toml
@@ -4,6 +4,7 @@ authors = ["Prashanth Prakash <prash@mit.edu> and contributors"]
 version = "1.0.0"
 
 [deps]
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
diff --git a/README.md b/README.md
index 1f922f1..ef8c70c 100644
--- a/README.md
+++ b/README.md
@@ -10,4 +10,4 @@ This package is primarily intended for use in gas turbine cycle calculations wit
 
 This package also introduces a method for fast calculations of thermodynamic properties of mixtures with fixed composition by calculating an equivalent polynomial representation of the mixture. See docs for theory and performance details. 
 
- 
+This package uses parametric typing to make it compatible with automatic forward-mode differentiation.
\ No newline at end of file
diff --git a/docs/src/gas.md b/docs/src/gas.md
index f5a66bb..cf3e451 100644
--- a/docs/src/gas.md
+++ b/docs/src/gas.md
@@ -20,7 +20,6 @@ information about the gas mixture.
 ```@docs
 Gas
 Gas()
-Base.setproperty!(gas::Gas, sym::Symbol, val::Float64)
 ```
 
 ## Single component gases
diff --git a/docs/src/index.md b/docs/src/index.md
index 2b9c26f..5b6f32f 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -5,8 +5,7 @@ of the gas/mixture is only a function of its temperature, i.e., $c_p(T)$ ,
 $h(T)$, and $s(T,p)$ (note the entropy is a function of both pressure and temperature).
 
 One of the important features of `IdealGasThermo.jl` is the ability to represent
-a mixture of gases as a single [`composite_species`](@ref vitiated).
-
+a mixture of gases as a single [`composite_species`](@ref vitiated). This package uses parametric typing to make it compatible with forward-mode automatic differentiation.
 
 ## Getting started
 
diff --git a/src/Gas.jl b/src/Gas.jl
index 2d2258a..6f39d1e 100644
--- a/src/Gas.jl
+++ b/src/Gas.jl
@@ -1,22 +1,20 @@
 abstract type AbstractGas end
 """
-    Gas{N}
+    Gas{N, R}
 
 A type that represents an ideal gas that is calorically perfect 
 i.e. ``c_p(T)``, ``h(T)``, ``\\phi(T)`` and ``s(T,P)``.
 """
-mutable struct Gas{N} <: AbstractGas
-    P::Float64 # [Pa]
-    T::Float64 # [K]
-    Tarray::MVector{8,Float64} # Temperature array to make calcs allocation free
-
-    cp::Float64 #[J/mol/K]
-    cp_T::Float64 # derivative dcp/dT
-    h::Float64  #[J/mol]
-    ϕ::Float64  #[J/mol/K] Entropy complement fn ϕ(T) = ∫ cp(T)/T dT from Tref to T
-
-    Y::MVector{N,Float64} # Mass fraction of species
-    MW::Float64 # Molecular weight [g/mol]
+mutable struct Gas{N, R<:Real} <: AbstractGas
+    P::R                   # Pressure [Pa]
+    T::R                   # Temperature [K]
+    Tarray::MVector{8, R}  # Preallocated temp array
+    cp::R                  # Heat capacity [J/mol/K]
+    cp_T::R                # Derivative dcp/dT
+    h::R                   # Enthalpy [J/mol]
+    ϕ::R                   # Entropy complement function
+    Y::MVector{N, R}       # Mass fractions
+    MW::R                  # Molecular weight [g/mol]
 end
 
 """
@@ -25,10 +23,10 @@ end
 Constructs `Gas` with given composition `Y`
 
 """
-function Gas(Y::AbstractVector)
-    gas = Gas()
-    gas.Y = convert(Vector{Float64}, Y)
-    set_TP!(gas, Tstd, Pstd) #setting temperature and pressure to recalculate thermodynamic properties
+function Gas(Y::AbstractVector{R}) where {R<:Real}
+  
+    gas = Gas(R(Tstd), R(Pstd))
+    gas.Y = SVector{Nspecies, R}(Y)
     return gas
 end
 
@@ -64,22 +62,52 @@ function Gas()
     Air = spdict[i]
     Y = zeros(Nspecies)
     Y[i] = 1.0
-    Gas{Nspecies}(
+    Gas{Nspecies, Float64}(
         Pstd,
         Tstd,
         Tarray(Tstd),
         Cp(Tstd, Air),
         (Cp(Tstd + 1.0, Air) - Cp(Tstd - 1.0, Air)) / 2.0, #finite diff dCp/dT
         h(Tstd, Air),
-        s(Tstd, Pstd, Air),
+        𝜙(Tstd, Air),
         Y,
-        Air.MW,
+        Air.MW
     )
 
 end
 
+"""
+    Gas(T, P)
+
+Constructor that returns a `Gas` type representing 
+Air at a desired temperature and pressure. It inherits
+the field types from the inputs and is recommended for 
+use with ForwardDiff.
+
+See also [`Gas`](@ref).
+"""
+function Gas(T::R, P::R) where R<:Real
+    i = findfirst(x -> x == "Air", spdict.name)
+    Air = spdict[i]
+
+    Y = zeros(R, Nspecies)
+    Y[i] = one(R)
+
+    Gas{Nspecies, R}(
+        P,
+        T,
+        Tarray(T),  # must return Vector{R}
+        Cp(T, Air),  # must return R
+        (Cp(T + one(R), Air) - Cp(T - one(R), Air)) / (2one(R)),  # finite diff in R
+        h(T, Air),   # must return R
+        𝜙(T, Air), # must return R
+        SVector{Nspecies, R}(Y),
+        Air.MW
+    )
+end
+
 # Overload Base.getproperty for convenience
-function Base.getproperty(gas::Gas, sym::Symbol)
+function Base.getproperty(gas::Gas{N, R}, sym::Symbol) where {N, R<:Real}
     if sym === :h_T # dh/dT
         return getfield(gas, :cp)
     elseif sym === :ϕ_T # dϕ/dT
@@ -110,24 +138,29 @@ function Base.getproperty(gas::Gas, sym::Symbol)
         # H += getproperty(gas, :h) * getproperty(gas, :MW)/1000.0
         return H
     elseif sym === :R #specific gas constant
-        return Runiv / getproperty(gas, :MW) * 1000.0
+        return R(Runiv / gas.MW * 1000.0)::R #Took some effort to make this stable!
     elseif sym === :γ
-        R = getproperty(gas, :R)
+        Rgas = R(getproperty(gas, :R))::R
         cp = getproperty(gas, :cp)
-        return cp / (cp - R)
+        return R(cp / (cp - Rgas))::R
     elseif sym === :gamma
         return getproperty(gas, :γ)
     elseif sym === :ρ
-        R = getproperty(gas, :R)
+        Rgas = R(getproperty(gas, :R))::R
         T = getfield(gas, :T)
         P = getfield(gas, :P)
-        return P / (R * T)
+        return R(P / (Rgas * T))::R
     elseif sym === :rho
         return getproperty(gas, :ρ)
     elseif sym === :ν
-        return 1 / getproperty(gas, :ρ)
+        return R(1 / getproperty(gas, :ρ))::R
     elseif sym === :nu
         return getproperty(gas, :ν)
+    elseif sym === :a #Speed of sound
+        Rgas = R(getproperty(gas, :R))::R
+        T = getfield(gas, :T)
+        γ = getproperty(gas, :γ)
+        return R(sqrt(γ * Rgas * T))::R
     elseif sym === :X # Get mole fractions
         Y = getfield(gas, :Y)
         MW = spdict.MW
@@ -148,148 +181,102 @@ function Base.getproperty(gas::Gas, sym::Symbol)
     end
 end
 
-
-"""
-    Base.setproperty!(gas::Gas, s::Symbol, val)
-
-"""
-function Base.setproperty!(gas::Gas, sym::Symbol, val::Float64)
-    ## Setting Temperature
+function Base.setproperty!(gas::Gas{N, R}, sym::Symbol, val::Real) where {N, R<:Real}
     if sym === :T
-        setfield!(gas, :T, val) # first set T
-        setfield!(gas, :Tarray, Tarray!(val, getfield(gas, :Tarray))) # update Tarray
-        TT = view(getfield(gas, :Tarray), :) # Just convenience
-        # Next set the cp, h and s of the gas
-        ## Get the right coefficients 
-        ## (assumes Tmid is always 1000.0. Check performed in readThermo.jl.):
-        if val < 1000.0 #Value given is the desired temperature
-            A = view(spdict.alow, :)
-        else
-            A = view(spdict.ahigh, :)
-        end
-        ## Initialize temporary vars
-        cptemp = 0.0
-        htemp = 0.0
-        ϕtemp = 0.0
-        cp_Ttemp = 0.0
-
-        P = getfield(gas, :P)
-        Y = getfield(gas, :Y)
-        MW = view(spdict.MW, :) # g/mol
-        # Go through every species where mass fraction is not zero
-        @inbounds for (Yᵢ, a, m) in zip(Y, A, MW)
-            if Yᵢ != 0.0
-                cptemp = cptemp + Yᵢ * Cp(TT, a) / m
-                htemp = htemp + Yᵢ * h(TT, a) / m
-                ϕtemp = ϕtemp + Yᵢ * 𝜙(TT, a) / m
-                cp_Ttemp = cp_Ttemp + Yᵢ * dCpdT(TT, a) / m
+        setfield!(gas, :T, val)
+        setfield!(gas, :Tarray, Tarray!(val, gas.Tarray))
+        TT = view(gas.Tarray, :)
+
+        A = val < 1000 ? view(spdict.alow, :) : view(spdict.ahigh, :)
+
+        cptemp = zero(R)
+        htemp = zero(R)
+        ϕtemp = zero(R)
+        cp_Ttemp = zero(R)
+
+        for (Yᵢ, a, m) in zip(gas.Y, A, spdict.MW)
+            if Yᵢ != 0
+                cptemp   += Yᵢ * Cp(TT, a) / m
+                htemp    += Yᵢ * h(TT, a) / m
+                ϕtemp    += Yᵢ * 𝜙(TT, a) / m
+                cp_Ttemp += Yᵢ * dCpdT(TT, a) / m
             end
         end
 
-        setfield!(gas, :cp, cptemp * 1000.0)
-        setfield!(gas, :h, htemp * 1000.0)
-        setfield!(gas, :ϕ, ϕtemp * 1000.0)
-        setfield!(gas, :cp_T, cp_Ttemp * 1000.0)
+        setfield!(gas, :cp, cptemp * 1000)
+        setfield!(gas, :h, htemp * 1000)
+        setfield!(gas, :ϕ, ϕtemp * 1000)
+        setfield!(gas, :cp_T, cp_Ttemp * 1000)
 
-        ## Setting Pressure
     elseif sym === :P
         setfield!(gas, :P, val)
-        TT = view(getfield(gas, :Tarray), :) # Just convenience
-        # Next set s of the gas
-        ## Get the right coefficients 
-        ## (assumes Tmid is always 1000.0. Check performed in readThermo.jl.):
-        if TT[4] < 1000.0 #TT[4] is T
-            A = view(spdict.alow, :)
-        else
-            A = view(spdict.ahigh, :)
-        end
-        ## Initialize temporary vars
-        ϕtemp = 0.0
-
-        P = val
-        Y = view(getfield(gas, :Y), :)
-        MW = view(spdict.MW, :) # g/mol
-        # Go through every species where mass fraction is not zero
-        @inbounds for (Yᵢ, a, m) in zip(Y, A, MW)
-            if Yᵢ != 0.0
-                ϕtemp = ϕtemp + Yᵢ * 𝜙(TT, a) / m
+        TT = gas.Tarray
+        A = TT[4] < 1000 ? view(spdict.alow, :) : view(spdict.ahigh, :)
+        ϕtemp = zero(R)
+
+        for (Yᵢ, a, m) in zip(gas.Y, A, spdict.MW)
+            if Yᵢ != 0
+                ϕtemp += Yᵢ * 𝜙(TT, a) / m
             end
         end
 
-        setfield!(gas, :ϕ, ϕtemp * 1000.0)
+        setfield!(gas, :ϕ, ϕtemp * 1000)
 
     elseif sym === :h
         set_h!(gas, val)
     elseif sym === :TP
         set_TP!(gas, val[1], val[2])
     else
-        error("You tried setting gas.", sym, " to a ", typeof(val))
+        setfield!(gas, sym, val)
     end
-    # Note: intentionally not including other variables to prevent 
-    # users from trying to directly set h, s, cp, MW etc.
-    nothing
+    return nothing
 end
 
-function Base.setproperty!(gas::Gas, sym::Symbol, val::AbstractVector{Float64})
-    if sym === :Y # directly set mass fractions Y
-        n = length(getfield(gas, :Y))
-        setfield!(gas, :Y, MVector{n}(val))
-        setfield!(gas, :MW, MW(gas)) # Update the MW of the gas mixture
-        setfield!(gas, :T, gas.T) # Force update of h
-
-    elseif sym === :X # directly set mole fractions Y
-        n = length(getfield(gas, :Y))
+function Base.setproperty!(gas::Gas{N, R1}, sym::Symbol, val::AbstractVector{<:R2}) where {N, R1<:Real, R2<:Real}
+    if sym === :Y
+        setfield!(gas, :Y, MVector{N, R1}(val))
+        setfield!(gas, :MW, MW(gas))
+        gas.T = gas.T #Reset T
+    elseif sym === :X
         Y = X2Y(val)
-        setfield!(gas, :Y, MVector{n}(Y))
-        setfield!(gas, :MW, MW(gas)) # Update the MW of the gas mixture
-        setfield!(gas, :T, gas.T) # Force update of h
+        setfield!(gas, :Y, MVector{N, R1}(Y))
+        setfield!(gas, :MW, MW(gas))
+        gas.T = gas.T #Reset T
     else
-        error(
-            "Only mass factions Y can be set with an array.",
-            "\n       You tried to set gas.",
-            sym,
-            " to a ",
-            typeof(val),
-        )
+        error("Only mass fractions Y/X can be set with a vector. Tried: $sym")
     end
-    nothing
+    return nothing
 end
 
-function Base.setproperty!(gas::Gas, sym::Symbol, val::AbstractDict{String,Float64})
-    if sym === :Y # directly set mass fractions Y
-        # If dict provided set each species in the right order
-        names = spdict.name
-        Y = zeros(MVector{length(names),Float64})
+function Base.setproperty!(gas::Gas{N, R1}, sym::Symbol, val::AbstractDict{String, R2}) where {N, R1<:Real, R2<:Real}
+    names = spdict.name
+    Y = zeros(MVector{N, R1})
+
+    if sym === :Y
         for (key, value) in val
-            index = findfirst(x -> x == key, names)
-            Y[index] = value
+            idx = findfirst(==(key), names)
+            Y[idx] = value
         end
-        setfield!(gas, :Y, Y)
-        setfield!(gas, :MW, MW(gas)) # Update the MW of the gas mixture
+        gas.Y = Y
+        gas.T = gas.T #Reset T
+        gas.MW = MW(gas)
+
     elseif sym === :X
-        names = spdict.name
-        X = zeros(MVector{length(names),Float64})
-        Y = zeros(MVector{length(names),Float64})
-        S = 0.0
+        X = zeros(MVector{N, R1})
+        S = zero(R1)
         for (key, value) in val
-            index = findfirst(x -> x == key, names)
+            idx = findfirst(==(key), names)
+            X[idx] = value
             S += value
-            X[index] = value
         end
-        X = X ./ S
-        Y .= X2Y(X)
-        setfield!(gas, :Y, Y)
-        setfield!(gas, :MW, MW(gas))
+        X ./= S
+        gas.Y = X2Y(X)
+        gas.T = gas.T #Reset T
+        gas.MW = MW(gas)
     else
-        error(
-            "Only mass factions Y can be set with a dict.",
-            "\n       You tried to set gas.",
-            sym,
-            " to a ",
-            typeof(val),
-        )
+        error("Only mass fractions Y/X can be set with a Dict. Tried: $sym")
     end
-    nothing
+    return nothing
 end
 
 """
@@ -329,7 +316,7 @@ with composition:
      ΣYᵢ     1.000     28.965
 ```
 """
-function set_h!(gas::AbstractGas, hspec::Float64)
+function set_h!(gas::AbstractGas, hspec::Real)
     T = gas.T
     dT = T
 
@@ -371,7 +358,7 @@ Sets the gas state based on a specified change in enthalpy (Δh) [J/mol],
 and a given polytropic efficiency. This represents adding or removing some work
 from the gas.
 """
-function set_Δh!(gas::AbstractGas, Δhspec::Float64, ηp::Float64 = 1.0)
+function set_Δh!(gas::AbstractGas, Δhspec::Real, ηp::Real = 1.0)
     P0 = gas.P
     ϕ0 = gas.ϕ
     hf = gas.h + Δhspec
@@ -384,7 +371,7 @@ end
 
 Calculates state of the gas given enthalpy and pressure (h,P)
 """
-function set_hP!(gas::AbstractGas, hspec::Float64, P::Float64)
+function set_hP!(gas::AbstractGas, hspec::Real, P::Real)
     set_h!(gas, hspec)
     gas.P = P
     return gas
@@ -416,7 +403,7 @@ with composition:
 ```
 
 """
-function set_TP!(gas::AbstractGas, T::Float64, P::Float64)
+function set_TP!(gas::AbstractGas, T::Real, P::Real)
     gas.T = T
     gas.P = P
     return gas
diff --git a/src/Gas1D.jl b/src/Gas1D.jl
index 0abcc69..7463f55 100644
--- a/src/Gas1D.jl
+++ b/src/Gas1D.jl
@@ -1,18 +1,18 @@
 """
-    Gas1D
+    Gas1D{R}
 
 Type that represents single component gases.
 """
-mutable struct Gas1D <: AbstractGas
-    comp_sp::composite_species
-    P::Float64 # [Pa]
-    T::Float64 # [K]
-    Tarray::MVector{8,Float64} # Temperature array to make calcs allocation free
-
-    cp::Float64 #[J/kg/K]
-    cp_T::Float64 # derivative dcp/dT
-    h::Float64  #[J/kg]
-    ϕ::Float64  #[J/kg/K] Entropy complement fn ϕ(T) = ∫ cp(T)/T dT from Tref to T
+mutable struct Gas1D{R<:Real} <: AbstractGas
+    comp_sp::composite_species{R}
+    P::R # [Pa]
+    T::R # [K]
+    Tarray::MVector{8,R} # Temperature array to make calcs allocation free
+
+    cp::R #[J/kg/K]
+    cp_T::R # derivative dcp/dT
+    h::R  #[J/kg]
+    ϕ::R  #[J/kg/K] Entropy complement fn ϕ(T) = ∫ cp(T)/T dT from Tref to T
     #Intentionally not storing s since log is expensive so only calculated when requested
 
 end
@@ -26,7 +26,7 @@ Dry Air at standard conditions
 See also [`Gas1D`](@ref).
 """
 function Gas1D()
-    Gas1D(
+    Gas1D{Float64}(
         DryAir,
         Pstd,
         Tstd,
@@ -34,7 +34,7 @@ function Gas1D()
         Cp(Tstd, DryAir),
         (Cp(Tstd + 1.0, DryAir) - Cp(Tstd - 1.0, DryAir)) / 2.0,
         h(Tstd, DryAir),
-        s(Tstd, Pstd, DryAir),
+        𝜙(Tstd, DryAir),
     )
 end
 
@@ -47,7 +47,7 @@ Dry Air at standard conditions
 See also [`Gas1D`](@ref).
 """
 function Gas1D(sp::composite_species)
-    Gas1D(
+    Gas1D{Float64}(
         sp,
         Pstd,
         Tstd,
@@ -55,11 +55,37 @@ function Gas1D(sp::composite_species)
         Cp(Tstd, sp),
         (Cp(Tstd + 1.0, sp) - Cp(Tstd - 1.0, sp)) / 2.0,
         h(Tstd, sp),
-        s(Tstd, Pstd, sp),
+        𝜙(Tstd, sp),
+    )
+end
+
+
+"""
+    Gas1D(T, P)
+
+Constructor that returns a `Gas` type representing 
+Dry Air at a desired temperature and pressure. It inherits the 
+field types from the inputs and is recommended for 
+use with ForwardDiff.
+
+See also [`Gas1D`](@ref).
+"""
+function Gas1D(T::R, P::R) where R<:Real
+    Xvec = Xidict2Array(Xair)
+    DA = generate_composite_species(R.(Xvec), "Dry Air")
+    Gas1D{R}(
+        DA,
+        P,
+        T,
+        Tarray(T),
+        Cp(T, DA),
+        (Cp(T + one(R), DA) - Cp(T - one(R), DA)) / 2.0,
+        h(T, DA),
+        𝜙(T, DA),
     )
 end
 
-function Base.setproperty!(gas::Gas1D, sym::Symbol, val::Float64)
+function Base.setproperty!(gas::Gas1D{R}, sym::Symbol, val::Real) where {R<:Real}
     ## Setting Temperature
     if sym === :T
         setfield!(gas, :T, val) # first set T
@@ -76,7 +102,7 @@ function Base.setproperty!(gas::Gas1D, sym::Symbol, val::Float64)
         end
 
         P = getfield(gas, :P)
-        MW = sp.MW / 1000.0 # kg/mol
+        MW = R(sp.MW / 1000.0) # kg/mol
 
         setfield!(gas, :cp, Cp(TT, A) / MW)
         setfield!(gas, :h, h(TT, A) / MW)
@@ -98,8 +124,8 @@ function Base.setproperty!(gas::Gas1D, sym::Symbol, val::Float64)
     nothing
 end
 
-# Overload Base.getproperty for convinence
-function Base.getproperty(gas::Gas1D, sym::Symbol)
+# Overload Base.getproperty for convenience
+function Base.getproperty(gas::Gas1D{R}, sym::Symbol) where R<:Real
     if sym === :h_T # dh/dT
         return getfield(gas, :cp)
     elseif sym === :ϕ_T # dϕ/dT
@@ -111,30 +137,35 @@ function Base.getproperty(gas::Gas1D, sym::Symbol)
     elseif sym === :TP
         return [getfield(gas, :T), getfield(gas, :P)]
     elseif sym === :s
-        Rgas = Runiv / getfield(gas, :comp_sp).MW * 1000.0
+        Rgas = R(Runiv) / getfield(gas, :comp_sp).MW * 1000.0
         return getfield(gas, :ϕ) - Rgas * log(getfield(gas, :P) / Pstd)
     elseif sym === :MW
         sp = getfield(gas, :comp_sp)
         return sp.MW
     elseif sym === :R #specific gas constant
-        return Runiv / getproperty(gas, :MW) * 1000.0
+        return R(Runiv / getproperty(gas, :MW) * 1000.0)::R
     elseif sym === :γ
-        R = getproperty(gas, :R)
+        Rgas = getproperty(gas, :R)
         cp = getproperty(gas, :cp)
-        return cp / (cp - R)
+        return R(cp / (cp - Rgas))::R
     elseif sym === :gamma
         return getproperty(gas, :γ)
     elseif sym === :ρ
-        R = getproperty(gas, :R)
+        Rgas = getproperty(gas, :R)
         T = getfield(gas, :T)
         P = getfield(gas, :P)
-        return P / (R * T)
+        return R(P / (Rgas * T))::R
     elseif sym === :rho
         return getproperty(gas, :ρ)
     elseif sym === :ν
-        return 1 / getproperty(gas, :ρ)
+        return R(1 / getproperty(gas, :ρ))::R
     elseif sym === :nu
         return getproperty(gas, :ν)
+    elseif sym === :a #Speed of sound
+        Rgas = R(getproperty(gas, :R))::R
+        T = getfield(gas, :T)
+        γ = getproperty(gas, :γ)
+        return R(sqrt(γ * Rgas * T))::R
     else
         return getfield(gas, sym)
     end
diff --git a/src/IdealGasThermo.jl b/src/IdealGasThermo.jl
index 2949284..9db4274 100644
--- a/src/IdealGasThermo.jl
+++ b/src/IdealGasThermo.jl
@@ -8,6 +8,7 @@ const default_thermo_path = joinpath(__Gasroot__, "data/thermo.inp")
 using LinearAlgebra
 using StaticArrays
 using Printf
+using ForwardDiff
 
 export Gas, set_h!, set_hP!, set_TP!, set_Δh!
 
@@ -34,5 +35,4 @@ gas.X = Xair
 const DryAir = generate_composite_species(gas.X, "Dry Air")
 export DryAir
 include("humidity.jl")
-
 end
diff --git a/src/atmosphere.jl b/src/atmosphere.jl
index 7a422c8..e37335c 100644
--- a/src/atmosphere.jl
+++ b/src/atmosphere.jl
@@ -5,14 +5,16 @@ US Standard Atmosphere. Model valid between -5 and 86 km in geometric altitude.
 geometric altitude ("alt_type == "geometric") or geopotential altitude ("alt_type == "geopotential"). 
 Model from: `https://ntrs.nasa.gov/citations/19770009539`
 """
-function standard_atmosphere(z, alt_type = "geometric")
+function standard_atmosphere(z::R, alt_type = "geometric") where R<:Real
 
     g = 9.80665 #Acceleration of gravity on Earth surface (z = 0)
     R_e = 6.356766e6 #Radius of Earth
 
     #Initialize gas (composition does not change)
-    gas = Gas1D()
-    R = gas.R
+    T = R(Tstd)
+    P = R(Pstd)
+    gas = Gas(T, P)
+    Rg = gas.R
 
     #Calculate geopotential altitude
     if alt_type == "geometric"
@@ -44,12 +46,12 @@ function standard_atmosphere(z, alt_type = "geometric")
     #Find temperature and pressure using the lapse rate
     #Use closed-form solutions
     if λ != 0
-        P = P0 * ((T0 + λ * (H - H0)) / T0)^(-g / (R * λ))
+        P = P0 * ((T0 + λ * (H - H0)) / T0)^(-g / (Rg * λ))
         T = T0 + λ * (H - H0)
 
     elseif λ == 0
         T = T0
-        P = P0 * exp(-g * (H - H0) / (R * T))
+        P = P0 * exp(-g * (H - H0) / (Rg * T))
     end
 
     #set the gas to the correct T and P
diff --git a/src/combustion.jl b/src/combustion.jl
index e2169ef..3ff5583 100644
--- a/src/combustion.jl
+++ b/src/combustion.jl
@@ -184,7 +184,7 @@ julia> IdealGasThermo.stoich_molar_FOR(CH4)
 """
 function stoich_molar_FOR(fuel::AbstractSpecies, oxidizer::AbstractSpecies = DryAir)
     molFuelOxyRatio = stoich_molar_fuel_oxy_ratio(fuel.name)
-    if typeof(oxidizer) == species
+    if oxidizer isa species
         Xin = Dict(oxidizer.name => 1.0)
         if oxidizer.name == "Air"
             Xin = Xair
@@ -218,8 +218,8 @@ stoich_FOR(fuel::AbstractString, oxi::AbstractString) =
     stoich_FOR(species_in_spdict(fuel), species_in_spdict(oxi))
 
 """
-    vitiated_mixture(fuel::AbstractSpecies, oxidizer::AbstractSpecies, 
-    FAR::Float64, ηburn::Float64=1.0)
+    vitiated_mixture(fuel, oxidizer, 
+    FAR, ηburn=1.0)
 
 Calculates the composition of a burnt gas mixture. Defaults to stoichiometric 
 conditions if FAR is not specified. `vitiated_mixture` returns the number of 
@@ -255,11 +255,11 @@ See [here](@ref vitiated) for some explanation of the background.
 function vitiated_mixture(
     fuel::AbstractSpecies,
     oxidizer::AbstractSpecies,
-    FAR::Float64,
-    ηburn::Float64 = 1.0,
-)
+    FAR::R1,
+    ηburn::R2 = 1.0,
+) where {R1 <: Real, R2 <: Real}
 
-    if typeof(oxidizer) == species
+    if oxidizer isa species
         Xin = Dict(oxidizer.name => 1.0)
         if oxidizer.name == "Air"
             Xin = Xair
@@ -279,7 +279,7 @@ function vitiated_mixture(
 
     Xdict = Dict(zip(names, ΔX))
 
-    Xvitiated::Dict{String,Float64} = mergewith(+, Xin, Xdict)
+    Xvitiated::Dict{String, R1} = mergewith(+, Xin, Xdict)
 
     for key in keys(Xvitiated) #Normalize such that sum(Xi) = 1
         Xvitiated[key] = Xvitiated[key] / (1 + ηburn * molFAR)
@@ -292,16 +292,16 @@ vitiated_mixture(fuel, oxidizer) =
     vitiated_mixture(fuel, oxidizer, stoich_FOR(fuel, oxidizer))
 
 """
-    vitiated_mixture(fuel::AbstractString, oxidizer::AbstractString, 
-    FAR::Float64, ηburn::Float64=1.0)
+    vitiated_mixture(fuel, oxidizer, 
+    FAR, ηburn=1.0)
 
 Convenience function that finds fuel and oxidizer from thermo database
 """
 function vitiated_mixture(
     fuel::AbstractString,
     oxidizer::AbstractString,
-    FAR::Float64,
-    ηburn::Float64 = 1.0,
+    FAR::Real,
+    ηburn::Real = 1.0,
 )
 
     fuel = species_in_spdict(fuel)
@@ -373,6 +373,7 @@ function AFT(fuel::AbstractSpecies, oxidizer::AbstractSpecies = DryAir)
 
     gas.X = dict
     gas.T = 298.15
+    h0 = gas.h
     Hfr = gas.Hf
     # Hr298 = gas.h*gas.MW/1000.0
 
@@ -382,15 +383,15 @@ function AFT(fuel::AbstractSpecies, oxidizer::AbstractSpecies = DryAir)
 
     #Energy Balance: Hrf + H(298.15) - Hr298.15 = Hpf + H(T) - Hp298.15
     #                Hrf - Hpf = ΔH(T)
-    ΔHpsens = Hfr - Hfp
+    hf = Hfr - Hfp + h0
     # set_h!(gas, Hpsens)
-    set_Δh!(gas, ΔHpsens, 0.0)
+    set_h!(gas, hf)
 
     return gas.T
 end  # function AFT
 """
-    vitiated_species(fuel::AbstractSpecies, oxidizer::AbstractSpecies, 
-    FAR::Float64, ηburn::Float64=1.0, name::AbstractString="vitiated species")
+    vitiated_species(fuel, oxidizer, 
+    FAR, ηburn=1.0, name="vitiated species")
 
 Returns a [`composite_species`](@ref) that represents the burnt gas mixture
 at the specified FAR. If no FAR is provided stoichiometeric conditions are
@@ -415,14 +416,14 @@ with composition:
 function vitiated_species(
     fuel::AbstractSpecies,
     oxidizer::AbstractSpecies,
-    FAR::Float64;
-    ηburn::Float64 = 1.0,
+    FAR::R1;
+    ηburn::R2 = 1.0,
     name::AbstractString = "vitiated species",
-)
+) where {R1 <: Real, R2 <: Real}
 
     Xdict = vitiated_mixture(fuel, oxidizer, FAR, ηburn)
 
-    X = zeros(Float64, Nspecies)
+    X = zeros(R1, Nspecies)
     Xidict2Array!(Xdict, X)
 
     return generate_composite_species(X, name)
@@ -431,8 +432,8 @@ end  # function vitiated_species
 function vitiated_species(
     fuel::AbstractString,
     oxidizer::AbstractString,
-    FAR::Float64;
-    ηburn::Float64 = 1.0,
+    FAR::Real;
+    ηburn::Real = 1.0,
     name::AbstractString = "vitiated species",
 )
 
@@ -449,9 +450,9 @@ vitiated_species(f::AbstractString, o::AbstractString) =
     vitiated_species(f, o, stoich_FOR(f, o))
 
 """
-    fixed_fuel_vitiated_species(fuel, oxidizer, ηburn::Float64=1.0)
+    fixed_fuel_vitiated_species(fuel, oxidizer, ηburn=1.0)
 
-Returns a function `burntgas(FAR::Float64)` that is specific to the fuel and oxidizer
+Returns a function `burntgas(FAR::Real)` that is specific to the fuel and oxidizer
 combination provided. This gives a highly performant function that can 
 simply be called at any given FAR for that specific fuel+oxidizer combo.
 
@@ -481,13 +482,13 @@ Gas1D(burntgas(CH4 + Dry Air; 0.02); MW = 28.51473501878705 g/mol)
 at T = 298.15 K; P = 101.325 kPa
 ```
 """
-function fixed_fuel_vitiated_species(fuel, oxidizer, ηburn::Float64 = 1.0)
+function fixed_fuel_vitiated_species(fuel, oxidizer, ηburn::R = 1.0) where R <: Real
 
     nCO2, nN2, nH2O, nO2 = ηburn .* reaction_change_molar_fraction(fuel.name)
     nFuel = (1.0 - ηburn)
     massratio = oxidizer.MW / fuel.MW
-    if typeof(oxidizer) == species
-        Xin::Dict{String,Float64} = Dict(oxidizer.name => 1.0)
+    if oxidizer isa species
+        Xin::Dict{String,Real} = Dict(oxidizer.name => 1.0)
         if oxidizer.name == "Air"
             Xin = Xair
         end
@@ -495,11 +496,11 @@ function fixed_fuel_vitiated_species(fuel, oxidizer, ηburn::Float64 = 1.0)
         Xin = oxidizer.composition
     end
     names::Vector{String} = [fuel.name, "CO2", "H2O", "N2", "O2"]
-    ΔX::Vector{Float64} = [nFuel, nCO2, nH2O, nN2, nO2]
-    Xdict::Dict{String,Float64} = Dict(zip(names, ΔX))
+    ΔX::Vector{R} = [nFuel, nCO2, nH2O, nN2, nO2]
+    Xdict::Dict{String,Real} = Dict(zip(names, ΔX))
 
-    ΔX_array = zeros(Float64, Nspecies)
-    Xin_array = zeros(Float64, Nspecies)
+    ΔX_array = zeros(R, Nspecies)
+    Xin_array = zeros(R, Nspecies)
     allnames = view(spdict.name, :)
 
     for (key, value) in Xdict
@@ -516,7 +517,7 @@ function fixed_fuel_vitiated_species(fuel, oxidizer, ηburn::Float64 = 1.0)
        burntgas(FAR::Float64)
     Inner function that will be returned
     """
-    function burntgas(FAR::Float64)
+    function burntgas(FAR::Real)
         molFAR = FAR .* massratio
         X = ΔX_array .* molFAR
         @. X = (X + Xin_array) / (1 + molFAR)
@@ -528,43 +529,41 @@ function fixed_fuel_vitiated_species(fuel, oxidizer, ηburn::Float64 = 1.0)
 end  # function fixed_fuel_vitiated_species
 
 """
-    fuel_combustion(gas_ox::AbstractGas, fuel::String, Tf::Float64, FAR::Float64, 
-    ηburn::Float64 = 1.0, hvap::Float64 = 0.0)
+    fuel_combustion(gas_ox, fuel, Tf, FAR, 
+    ηburn = 1.0, hvap = 0.0)
 
 This function returns an `AbstractGas` with the combustion products at the combustor exit temperature
 for a given fuel type, oxidizer gas at a given enthalpy, and FAR. It includes the combustion 
 efficiency and the fuel enthalpy of vaporization as optional inputs. 
 """
 function fuel_combustion(
-    gas_ox::AbstractGas,
+    gas_ox::Gas{N, R},
     fuel::String,
-    Tf::Float64,
-    FAR::Float64,
-    ηburn::Float64 = 1.0,
-    hvap::Float64 = 0.0,
-)
-
+    Tf::Real,
+    FAR::Real,
+    ηburn::Real = 1.0,
+    hvap::Real = 0.0
+) where {N, R<:Real}
     #Create variables corresponding to the oxidizer and fuel species and mixtures
     fuel_sps = species_in_spdict(fuel)
 
     #Find oxidizer composition
-    if typeof(gas_ox) == Gas1D
-        gas_sps = gas_ox.comp_sp
+    if gas_ox isa Gas1D
+        gas_sps = gas_ox.comp_sp::composite_species{R}
     else
         if "Air" in keys(gas_ox.Xdict)
-            Xin = Xair
+            Xin = Xair_typed(R)::Dict{String, R}
         else
-            Xin = gas_ox.Xdict
+            Xin = gas_ox.Xdict::Dict{String, R}
         end
-        gas_sps = generate_composite_species(IdealGasThermo.Xidict2Array(Xin))
+        gas_sps = generate_composite_species(IdealGasThermo.Xidict2Array(Xin))::composite_species{R}
     end
 
     #Find the vectors with the fuel mole and mass fractions
-    Xfuel = Xidict2Array(Dict([(fuel, 1.0)])) #Mole fraction
+    Xfuel = Xidict2Array(Dict{String, R}([(fuel, R(1.0))])) #Mole fraction
     Yfuel = X2Y(Xfuel) #Mass fraction
-    gas_fuel = Gas(Yfuel) #Create a fuel gas to calculate enthalpy
-
-    set_TP!(gas_fuel, Tf, gas_ox.P) #Set the fuel gas at the correct conditions
+    gas_fuel = Gas(Tf, gas_ox.P) #Create a fuel gas to calculate enthalpy
+    gas_fuel.Y = Yfuel
 
     #Find product composition and mole and mass fractions
     Xprod_dict = vitiated_mixture(fuel_sps, gas_sps, FAR, ηburn)
@@ -572,55 +571,55 @@ function fuel_combustion(
     Yprod = X2Y(Xprod)
 
     #Initialize output 
-    gas_prod = Gas(Yprod)
+    gas_prod = Gas(gas_ox.T, gas_ox.P) #Initialize with oxidizer properties
+    gas_prod.Y = Yprod #Set correct composition
 
     # Combustion enthalpy calculations
     #enthalpy before reaction = enthalpy after reaction
     h0 = (gas_ox.h + FAR * (gas_fuel.h - abs(hvap))) / (1 + FAR)
 
-    #Set the products at the correct enthalpy and pressure
-    set_hP!(gas_prod, h0, gas_ox.P)
+    #Set the products at the correct enthalpy
+    set_h!(gas_prod, h0)
     return gas_prod
 end
 
 """
-    gas_burn(gas_ox::AbstractGas, fuel::String, Tf::Float64, Tburn::Float64, 
-    ηburn::Float64 = 1.0, hvap::Float64 = 0.0)
+    gas_burn(gas_ox, fuel, Tf, Tburn, 
+    ηburn = 1.0, hvap = 0.0)
 
 This function returns a tuple containing the FAR and an `AbstractGas` with the combustion 
 products for a given fuel type, oxidizer gas at a given enthalpy, and desired combustor exit temperature. 
 It includes the combustion efficiency and the fuel enthalpy of vaporization as optional inputs. 
 """
 function gas_burn(
-    gas_ox::AbstractGas,
+    gas_ox::Gas{N, R},
     fuel::String,
-    Tf::Float64,
-    Tburn::Float64,
-    ηburn::Float64 = 1.0,
-    hvap::Float64 = 0.0,
-)
+    Tf::Real,
+    Tburn::Real,
+    ηburn::Real = 1.0,
+    hvap::Real = 0.0
+) where {N, R<:Real}
 
     #Create variables corresponding to the oxidizer and fuel species and mixtures
     fuel_sps = species_in_spdict(fuel)
 
     #Extract composite species with oxidizer gas composition
-    if typeof(gas_ox) == Gas1D
+    if gas_ox isa Gas1D
         gas_sps = gas_ox.comp_sp
     else
         if "Air" in keys(gas_ox.Xdict)
-            Xin = Xair
+            Xin = Xair_typed(R)::Dict{String, R}
         else
-            Xin = gas_ox.Xdict
+            Xin = gas_ox.Xdict::Dict{String, R}
         end
-        gas_sps = generate_composite_species(IdealGasThermo.Xidict2Array(Xin))
+        gas_sps = generate_composite_species(IdealGasThermo.Xidict2Array(Xin))::composite_species{R}
     end
 
     #Find the vectors with the fuel mole and mass fractions
-    Xfuel = Xidict2Array(Dict([(fuel, 1.0)])) #Mole fraction
+    Xfuel = Xidict2Array(Dict{String, R}([(fuel, 1.0)])) #Mole fraction
     Yfuel = X2Y(Xfuel) #Mass fraction
-    gas_fuel = Gas(Yfuel) #Create a fuel gas to calculate enthalpy
-
-    set_TP!(gas_fuel, Tf, gas_ox.P) #Set the fuel gas at the correct conditions
+    gas_fuel = Gas(Tf, gas_ox.P)
+    gas_fuel.Y = Yfuel #Create a fuel gas to calculate enthalpy
 
     #Store enthalpies of oxidizer and fuel at original temperatures
     ho = gas_ox.h
@@ -631,14 +630,14 @@ function gas_burn(
 
     names = ["CO2", "H2O", "N2", "O2"]
     ΔX = [nCO2, nH2O, nN2, nO2]
-    Xdict = Dict(zip(names, ΔX))
+    Xdict = Dict{String, R}(zip(names, ΔX))
 
     Xc = Xidict2Array(Xdict)
     Yc = X2Y(Xc) #mass fraction change in combustion for FAR = 1
 
-    gas_c = Gas(Yc) #Create a "virtual" gas with the changes in combustion, for enthalpy
+    gas_c = Gas(Tburn, gas_ox.P) #Create a "virtual" gas with the changes in combustion, for enthalpy
     #calculations
-    set_TP!(gas_c, Tburn, gas_ox.P)
+    gas_c.Y = Yc
 
     hc = gas_c.h #Enthalpy change for FAR = 1
 
@@ -658,10 +657,10 @@ function gas_burn(
     Yprod = X2Y(Xprod)
 
     #Initialize output 
-    gas_prod = Gas(Yprod)
+    gas_prod = Gas(Tburn, gas_ox.P)
 
-    #Set the products at the correct enthalpy and pressure
-    set_TP!(gas_prod, Tburn, gas_ox.P)
+    #Set the correct composition
+    gas_prod.Y = Yprod
 
     return FAR, gas_prod
-end
+end
\ No newline at end of file
diff --git a/src/constants.jl b/src/constants.jl
index 0b69e5d..629220a 100644
--- a/src/constants.jl
+++ b/src/constants.jl
@@ -4,21 +4,22 @@ const Tstd = 298.15 # K
 const ϵ = sqrt(eps()) #standard tolerance 
 
 # Air composition
-const Xair = Dict(
-    "N2" => 0.78084,
-    "Ar" => 0.009365,
-    "H2O" => 0.0,
-    "CO2" => 0.000319,
-    "O2" => 0.209476,
-)
+function Xair_typed(R)
+    Xair = Dict{String, R}(
+        "N2" => 0.78084,
+        "Ar" => 0.009365,
+        "H2O" => 0.0,
+        "CO2" => 0.000319,
+        "O2" => 0.209476
+    )
+    return Xair
+end
+
+const Xair = Xair_typed(Float64)
+
+const XwetAir = Dict{String, Float64}(
 
-const XwetAir = Dict(
-    "N2" => 0.78084,
-    "Ar" => 0.009365,
-    "H2O" => 0.018722,
-    "CO2" => 0.000319,
-    "O2" => 0.209476,
 )
 
 const Ytasopt =
-    Dict("N2" => 0.7532, "O2" => 0.2315, "CO2" => 0.0006, "H2O" => 0.0020, "Ar" => 0.0127)
+    Dict{String, Float64}("N2" => 0.7532, "O2" => 0.2315, "CO2" => 0.0006, "H2O" => 0.0020, "Ar" => 0.0127)
diff --git a/src/humidity.jl b/src/humidity.jl
index f66a4b1..f1317c1 100644
--- a/src/humidity.jl
+++ b/src/humidity.jl
@@ -48,12 +48,12 @@ function generate_humid_air(
     RH::type,
     T::type = Tstd,
     P::type = Pstd,
-) where {type<:AbstractFloat}
+) where {type<:Real}
     q = specific_humidity(RH, T, P)
     Xwater = q / ε
     Xdict = mergewith(+, Xair, Dict("H2O" => Xwater))
 
-    X = zeros(Float64, Nspecies)
+    X = zeros(type, Nspecies)
     Xidict2Array!(Xdict, X)
 
     return generate_composite_species(X, "Wet air with RH = $RH at ($T K; $(P/1000.0) kPa)")
diff --git a/src/readThermo.jl b/src/readThermo.jl
index 5f6d378..7ba389f 100644
--- a/src/readThermo.jl
+++ b/src/readThermo.jl
@@ -70,7 +70,7 @@ function readThermo(filename)
                 [coeffs[(i-1)*16+1:i*16] for i in range(1, 10, step = 1) if i != 8],
             )
 
-        Species[i] = species(sp, Tmid, a[1, :], a[2, :], MW, Hf, formula)
+        Species[i] = species(String(sp), Tmid, a[1, :], a[2, :], MW, Hf, formula)
     end
     return StructArray(Species)
 end
diff --git a/src/species.jl b/src/species.jl
index 9407d91..2d740d8 100644
--- a/src/species.jl
+++ b/src/species.jl
@@ -29,14 +29,14 @@ defining ``c_p``, ``h``, and ``s``.
 
 See [here](@ref gas1dthermo) for a more detailed explanation.
 """
-struct composite_species <: AbstractSpecies
+struct composite_species{R<:Real} <: AbstractSpecies
     name::String
-    Tmid::Float64
-    alow::Array{Float64,1}
-    ahigh::Array{Float64,1}
-    MW::Float64
-    Hf::Float64
-    composition::Dict{String,Float64}
+    Tmid::R
+    alow::Array{R,1}
+    ahigh::Array{R,1}
+    MW::R
+    Hf::R
+    composition::Dict{String,R}
 end
 """
     generate_composite_species(Xi::AbstractVector, name::AbstractString="composite species")
@@ -45,15 +45,15 @@ Generates a composite psuedo-species to represent a gas mixture given the
 mole fraction `Xi` of its constitutents.
 """
 function generate_composite_species(
-    Xi::AbstractVector,
+    Xi::AbstractVector{R},
     name::AbstractString = "composite species",
-)
+) where R<:Real
     ALOW = reduce(hcat, spdict.alow)
     AHIGH = reduce(hcat, spdict.ahigh)
     alow = ALOW * Xi
     ahigh = AHIGH * Xi
-    MW = dot(spdict.MW, Xi)
-    Hf = dot(spdict.Hf, Xi)
+    MW = R(dot(spdict.MW, Xi))
+    Hf = R(dot(spdict.Hf, Xi))
     # Need to account for the entropy of mixing:
     Δs_mix = 0.0
     for i in eachindex(Xi)
@@ -66,7 +66,7 @@ function generate_composite_species(
     ahigh[end] = ahigh[end] - Δs_mix
     Tmid = 1000.0 #Assumed to always be the mid (this is checked for all input thermo when read)
     # comp
-    d = Dict()
+    d = Dict{String, R}()
     if !(sum(Xi) ≈ 1.0)
         error("Gas mixture composition is not well defined. Sum of Xi = $(sum(Xi)) != 1.0")
     end
@@ -78,6 +78,6 @@ function generate_composite_species(
             push!(d, spdict.name[i] => Xi[i])
         end
     end
-    return composite_species(name, Tmid, alow, ahigh, MW, Hf, d)
+    return composite_species{R}(name, Tmid, alow, ahigh, MW, Hf, d)
 
 end  # function generate_composite_species
diff --git a/src/thermoProps.jl b/src/thermoProps.jl
index 0717f04..87e31ba 100644
--- a/src/thermoProps.jl
+++ b/src/thermoProps.jl
@@ -8,7 +8,7 @@ Calculates cp of the given species in J/K/mol
 Cp0/R = a₁T⁻² + a₂T⁻¹ + a₃ + a₄T + a₅T² + a₆T³ + a₇T⁴
 ```
 """
-function Cp(Tarray::AbstractVector{T}, a::AbstractVector{T}) where {T}
+function Cp(Tarray::AbstractVector{T}, a::AbstractVector{T1}) where {T, T1}
     Cp_R = dot(view(a, 1:7), view(Tarray, 1:7))
     Cp = Cp_R * Runiv
     return Cp #J/K/mol
@@ -22,7 +22,7 @@ Returns the derivative dcp/dT [J/K²/mol]
 dCp0/dT = R(-2a1*T^-3 -a2*T^-2 + a4 + 2a5*T + 3a6*T^2 + 4a7*T^3)
 ```
 """
-function dCpdT(TT::AbstractVector{T}, a::AbstractVector{T}) where {T}
+function dCpdT(TT::AbstractVector{T}, a::AbstractVector{T1}) where {T, T1}
     dcp_RdT =
         -2 * a[1] * TT[1] * TT[2] - a[2] * TT[1] +
         a[4] +
@@ -44,7 +44,7 @@ H0/RT = -a1*T^-2 + a2*T^-1*ln(T) + a3 + a4*T/2 + a5*T^2/3 + a6*T^3/4 + a7*T^4/5
       = -a1*T₁   + a2*T₂*T₈      + a3 + a4*T₄/2 + a5*T₅/3  + a6*T₆/4  + a7*T₇/5  + a₈*T₂
 ```
 """
-function h(TT::AbstractVector{type}, a::AbstractVector{type}) where {type}
+function h(TT::AbstractVector{type}, a::AbstractVector{type1}) where {type, type1}
     h_RT = -a[1] * TT[1] +
             a[2] * TT[8] * TT[2] +
             a[3] +
@@ -69,7 +69,7 @@ S0/R = -a1*T^-2/2 - a2*T^-1 + a3*ln(T) + a4*T + a5*T^2/2 + a6*T^3/3.0 + a7*T^4/4
      = -a1*T₁/2   - a2*T₂   + a3*T₈    + a4*T₄+ a5*T₅/2  + a6*T₆/3.0  + a7*T₇/4  + a₉   
 ```
 """
-function 𝜙(TT::AbstractVector{type}, a::AbstractVector{type}) where {type}
+function 𝜙(TT::AbstractVector{type}, a::AbstractVector{type1}) where {type, type1}
     so_R =
         -0.5 * a[1] * TT[1] - a[2] * TT[2] +
         a[3] * TT[8] +
@@ -134,3 +134,21 @@ function s(T, P, sp::AbstractSpecies)
     sᵒ = 𝜙(TT, a) - Runiv * log(P / Pstd)
     return sᵒ * 1000.0 / sp.MW
 end
+
+"""
+    𝜙(T, P, sp::AbstractSpecies)
+
+Calculates the entropy complement function 𝜙 for a species in J/K/kg
+"""
+function 𝜙(T::Real, sp::AbstractSpecies)
+    TT = Tarray(T)
+    if T < 1000.0
+        s = :alow
+    else
+        s = :ahigh
+    end
+    a = getfield(sp, s)
+    phi = 𝜙(TT, a)
+    return phi * 1000.0 / sp.MW
+end
+
diff --git a/src/turbo.jl b/src/turbo.jl
index 4f0ed76..3a9316f 100644
--- a/src/turbo.jl
+++ b/src/turbo.jl
@@ -1,11 +1,11 @@
 # Turbomachinery related functions
 
 """
-    PressureRatio(gas::AbstractGas, PR::Float64, ηp::Float64=1.0,)
+    PressureRatio!(gas, PR, ηp=1.0,)
 
 Generic pressure ratio conversion. See [`compress`](@ref) and [`expand`](@ref)
 """
-function PressureRatio(gas::AbstractGas, PR::Float64, ηp::Float64 = 1.0)
+function PressureRatio!(gas::AbstractGas, PR::Real, ηp::Real)
 
     T0 = gas.T
     ϕ0 = gas.ϕ
@@ -45,44 +45,44 @@ function PressureRatio(gas::AbstractGas, PR::Float64, ηp::Float64 = 1.0)
         )
     end
 
-    return gas
-
 end
 
 """
-   compress(gas::AbstractGas, PR::Float64, ηp::Float64=1.0,)
+   compress!(gas, PR, ηp=1.0,)
 
 Compression with an optional polytropic efficiency.PR should be ≥ 1.0.
 See also [`expand`](@ref).
 """
-function compress(gas::AbstractGas, PR::Float64, ηp::Float64 = 1.0)
+function compress!(gas::AbstractGas, PR::Real, ηp::Real = 1.0)
+
     if PR < 1.0
         error("The specified pressure ratio (PR) to compress by needs to be ≥ 1.0.
         Provided PR = $PR. Did you mean to use `expand`?")
     end
-    return PressureRatio(gas, PR, ηp)
+    PressureRatio!(gas, PR, ηp)
 end
 
 """
-   expand(gas::AbstractGas, PR::Float64, ηp::Float64=1.0,)
+   expand!(gas, PR, ηp=1.0,)
 
 Expansion at a given polytropic efficiency. PR should be ≤ 1.0.
 See also [`compress`](@ref).
 """
-function expand(gas::AbstractGas, PR::Float64, ηp::Float64 = 1.0)
+function expand!(gas::AbstractGas, PR::Real, ηp::Real = 1.0)
     if PR > 1.0
         error("The specified pressure ratio (PR) to compress by needs to be ≤ 1.0.
         Provided PR = $PR. Did you mean to use `compress`?")
     end
-    return PressureRatio(gas, PR, 1 / ηp)
+    PressureRatio!(gas, PR, 1 / ηp)
 end
 
 """
-   gas_Mach!(gas::AbstractGas, M0::Float64, M::Float64, ηp::Float64 = 1.0)
+   gas_Mach!(gas, M0, M, ηp = 1.0)
 
 Calculates the gas state for a change in Mach number with an optional polytropic efficiency.
 """
-function gas_Mach!(gas::AbstractGas, M0::Float64, M::Float64, ηp::Float64 = 1.0)
+function gas_Mach!(gas::AbstractGas, M0::Real, M::Real, ηp::Real = 1.0)
+
     itmax = 10
     ttol = 0.000001
 
@@ -121,15 +121,15 @@ function gas_Mach!(gas::AbstractGas, M0::Float64, M::Float64, ηp::Float64 = 1.0
 end
 
 """
-   gas_mixing(gas1::AbstractGas, gas2::AbstractGas, mratio::Float64)
+   gas_mixing(gas1, gas2, mratio)
 
 Calculates the resulting gas after two gases (gas1 and gas2) are mixed at constant pressure, with a mass ratio
 mratio = mass of gas2 / mass gas1.
 """
-function gas_mixing(gas1::AbstractGas, gas2::AbstractGas, mratio::Float64)
+function gas_mixing(gas1::AbstractGas, gas2::AbstractGas, mratio::Real)
 
     #Extract dictionaries with gas molar fractions
-    if typeof(gas1) == Gas1D
+    if gas1 isa Gas1D
         X1 = gas1.comp_sp.composition
     else
         if "Air" in keys(gas1.Xdict)
@@ -139,7 +139,7 @@ function gas_mixing(gas1::AbstractGas, gas2::AbstractGas, mratio::Float64)
         end
     end
 
-    if typeof(gas2) == Gas1D
+    if gas2 isa Gas1D
         X2 = gas2.comp_sp.composition
     else
         if "Air" in keys(gas2.Xdict)
diff --git a/src/utils.jl b/src/utils.jl
index 67d2447..13fbe56 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -96,7 +96,7 @@ end  # function X2Y
 Convert a mole fraction dictonary into the given array X with the right order of
 compounds.
 """
-function Xidict2Array!(Xdict::Dict{String,Float64}, X::AbstractVector)
+function Xidict2Array!(Xdict::Dict{String,R}, X::AbstractVector) where R <: Real
     names = spdict.name
     for (key, value) in Xdict
         index = findfirst(x -> x == key, names)
@@ -111,9 +111,9 @@ end  # function Xidict2Array
     Xidict2Array(Xdict::Dict{String, Float64})
 Converts the dict into a new array with mole fractions in the right order
 """
-function Xidict2Array(Xdict::Dict{String,Float64})
+function Xidict2Array(Xdict::Dict{String, R}) where R <: Real
     names = spdict.name
-    X = zeros(Float64, Nspecies)
+    X = zeros(R, Nspecies)
     for (key, value) in Xdict
         index = findfirst(x -> x == key, names)
         X[index] = value
diff --git a/test/runtests.jl b/test/runtests.jl
index bc74f9c..d3c0201 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -13,4 +13,5 @@ const IdealGasThermo = IdealGasThermo
     include("unit_test_turbo.jl")
     include("unit_test_atmos.jl")
     include("unit_test_utils.jl")
+    include("unit_test_autodiff.jl")
 end
diff --git a/test/unit_test_autodiff.jl b/test/unit_test_autodiff.jl
new file mode 100644
index 0000000..aab0dcd
--- /dev/null
+++ b/test/unit_test_autodiff.jl
@@ -0,0 +1,291 @@
+using ForwardDiff
+
+function FDjacobian(fun, x0)
+    delta = 1e-6 #Finite diff perturbation
+    f0 = fun(x0)
+    jac_FD = zeros(length(f0), length(x0))
+    for i in 1:length(x0)
+        x = deepcopy(x0)
+        x[i] = x[i]*(1+delta)
+        fi = fun(x)
+        jac_FD[:, i] .= (fi.-f0)./(x[i]-x0[i])
+    end
+    return jac_FD
+end
+@testset "autodiff" begin 
+
+    @testset "thermo properties" begin 
+
+        #Simple equation that returns enthalpy
+        function h(x)
+            T = x[1]
+            P = x[2]
+            
+            g = Gas(T, P)
+            h = g.h
+            return h
+        end
+
+        x0 = [288.15, 101325.0]
+        grad = ForwardDiff.gradient(h, x0)
+
+        #Compare gradient to exact solution
+        g = Gas(x0[1], x0[2])
+        grad_exact = [g.cp, 0.0]
+        for (i,g) in enumerate(grad)
+            @test g ≈ grad_exact[i] rtol = 1e-8
+        end
+
+        #Elaborate function of thermodynamic parameters
+        function thermofun(x)
+            T = x[1]
+            P = x[2]
+            
+            g = Gas(T, P)
+            s = g.s
+            h = g.h
+            return s^2/h
+        end
+
+        x0 = [288.15, 101325.0]
+        grad = ForwardDiff.gradient(thermofun, x0)
+
+        #Compare gradient to finite differences
+        grad_FD = FDjacobian(thermofun, x0)
+        for (i,g) in enumerate(grad)
+            @test g ≈ grad_FD[i] rtol = 1e-4
+        end
+    end
+
+    @testset "combustion" begin 
+        #Function that returns FAR
+        function FARfun(x)
+            T = x[1]
+            P = x[2]
+            Tf = x[3]
+            Tb = x[4]
+            etab = x[5]
+            hvap = x[6]
+            
+            g = Gas(T, P)
+            FAR,_ = IdealGasThermo.gas_burn(g,
+                "CH4",
+                Tf,
+                Tb, etab, hvap)
+            return FAR
+        end
+
+        x0 = [298.15, 101325.0, 298.15, 1000.0, 0.99, 1e5]
+        grad = ForwardDiff.gradient(FARfun, x0)
+
+        #Compare gradient to finite differences
+        grad_FD = FDjacobian(FARfun, x0)
+        for (i,g) in enumerate(grad)
+            @test g ≈ grad_FD[i] rtol = 1e-4
+        end
+
+        #Function that returns the product temperature
+        function Tflame(x)
+            T = x[1]
+            P = x[2]
+            Tf = x[3]
+            FAR = x[4]
+            etab = x[5]
+            hvap = x[6]
+
+            g = Gas(T, P)
+            g_prod = IdealGasThermo.fuel_combustion(g,
+                "CH4", Tf, FAR,
+                etab, hvap) 
+            return g_prod.T
+        end
+
+        x0 = [298.15, 101325.0, 298.15, 0.01, 0.99, 1e5]
+        grad = ForwardDiff.gradient(Tflame, x0)
+
+        #Compare gradient to finite differences
+        grad_FD = FDjacobian(Tflame, x0)
+        for (i,g) in enumerate(grad)
+            @test g ≈ grad_FD[i] rtol = 1e-4
+        end
+    end
+
+    @testset "turbo" begin 
+        function PRfun(x)
+            T = x[1]
+            P = x[2]
+            PR = x[3]
+            eta = x[4]
+
+            gas = Gas(T, P)
+
+            IdealGasThermo.PressureRatio!(gas, PR, eta)
+            return [gas.T, gas.P]
+        end
+
+        x0 = [298.15, 101325.0, 10.0, 0.9]
+        J = ForwardDiff.jacobian(PRfun, x0)
+
+        J_FD = FDjacobian(PRfun, x0)
+        for i in 1:2
+            for j in 1:length(x0)
+                @test J[i, j] ≈ J_FD[i, j] rtol = 1e-4
+            end
+        end
+
+        function Mach(x)
+            T = x[1]
+            P = x[2]
+            M0 = x[3]
+            M = x[4]
+
+            gas = Gas(T, P)
+            IdealGasThermo.gas_Mach!(gas, M0, M)
+            return [gas.T, gas.P]
+        end
+        x0 = [298.15, 101325.0, 0.1, 0.8]
+        J = ForwardDiff.jacobian(Mach, x0)
+
+        J_FD = FDjacobian(Mach, x0)
+        for i in 1:2
+            for j in 1:length(x0)
+                @test J[i, j] ≈ J_FD[i, j] rtol = 1e-4
+            end
+        end
+    end
+
+    @testset "Gas1D" begin 
+        #Elaborate function of thermodynamic parameters
+        function thermofun1D(x)
+            T = x[1]
+            P = x[2]
+            
+            g = Gas1D(T, P)
+            s = g.s
+            h = g.h
+            return s^2/h
+        end
+        x0 = [288.15, 101325.0]
+        grad = ForwardDiff.gradient(thermofun1D, x0)
+
+        #Compare gradient to finite differences
+        grad_FD = FDjacobian(thermofun1D, x0)
+        for (i,g) in enumerate(grad)
+            @test g ≈ grad_FD[i] rtol = 1e-4
+        end
+    end
+
+    @testset "atmos" begin 
+        #Returns temperature and pressure at a given altitude
+        function atmosfun(x)
+            z = x[1]
+            
+            g = IdealGasThermo.standard_atmosphere(z)
+            return [g.T, g.P]
+        end
+        x0 = [1e3]
+        J = ForwardDiff.jacobian(atmosfun, x0)
+
+        #Compare gradient to finite differences
+        J_FD = FDjacobian(atmosfun, x0)
+        for i in 1:2
+            for j in 1:length(x0)
+                @test J[i, j] ≈ J_FD[i, j] rtol = 1e-4
+            end
+        end
+    end
+
+    @testset "Gas1D" begin 
+        #Elaborate function of thermodynamic parameters
+        function thermofun1D(x)
+            T = x[1]
+            P = x[2]
+            
+            g = Gas1D(T, P)
+            s = g.s
+            h = g.h
+            return s^2/h
+        end
+        x0 = [288.15, 101325.0]
+        grad = ForwardDiff.gradient(thermofun1D, x0)
+
+        #Compare gradient to finite differences
+        grad_FD = FDjacobian(thermofun1D, x0)
+        for (i,g) in enumerate(grad)
+            @test g ≈ grad_FD[i] rtol = 1e-4
+        end
+    end
+
+    @testset "engine" begin 
+        #Basic Brayton thermodynamic cycle to test derivatives
+        function JetEngineResiduals(x, p)
+            PR_t = x[1]
+            mdot_c = x[2]
+            R = typeof(PR_t)
+
+            z0 = p[1]
+            M0 = p[2]
+            PR = R(p[3])
+            Tt4 = R(p[4])
+            F_N = p[5]
+            
+            #Intake
+            g0 = IdealGasThermo.standard_atmosphere(z0)
+            g0 = Gas(R(g0.T), R(g0.P))
+            u0 = M0*sqrt(g0.T*g0.gamma*g0.R)
+            gt2 = deepcopy(g0) #Copy so as to not modify it
+            IdealGasThermo.gas_Mach!(gt2, M0, 0.0) #stagnation properties
+            #Compressor
+            gt3 = deepcopy(gt2) #Copy so as to not modify it
+            IdealGasThermo.compress!(gt3, PR)
+            #Burner
+            FAR,_ = IdealGasThermo.gas_burn(gt3,
+                "CH4",
+                R(298.15),
+                Tt4, 1.0, 0.0)
+            gt4 = IdealGasThermo.fuel_combustion(gt3,
+                "CH4", R(298.15), FAR,
+                1.0, 0.0)
+            #Turbine
+            gt5 = deepcopy(gt4) #Make a copy to avoid modifying it
+            IdealGasThermo.expand!(gt5, PR_t) #Expand in turbine
+            #Nozzle
+            g6 = deepcopy(gt5) #Make a copy to avoid modifying it
+            PR_n = gt2.P / gt5.P #Nozzle pressure ratio
+            IdealGasThermo.expand!(g6, PR_n)
+
+            #Residuals
+            u = sqrt(2 * (gt5.h - g6.h))
+            Fsp = (1 + FAR)*u - u0
+            
+            return [F_N - mdot_c*Fsp; gt3.h - gt2.h - (1+FAR)*(gt4.h - gt5.h)] #Thrust and shaft power residuals
+        end
+        x0 = [0.5, 1.0]
+        p = [1e4, 0.8, 10.0, 1700.0, 1e4]
+        f(x) = JetEngineResiduals(x, p)
+        J = ForwardDiff.jacobian(f, x0)
+
+        #Compare gradient to finite differences
+        J_FD = FDjacobian(f, x0)
+        for i in 1:2
+            for j in 1:length(x0)
+                @test J[i, j] ≈ J_FD[i, j] rtol = 1e-4
+            end
+        end
+
+        #Solve problem for PR_t and mdot_c
+        eps = 1e-6
+        f0 = [1.0, 1.0]
+        x = x0
+        while maximum(abs.(f0)) > eps #Use a basic Newton's method
+            f0 = f(x)
+            J = ForwardDiff.jacobian(f, x)
+            dx = -J\f0
+            x = x + dx
+        end
+        x_check = [0.6183259629412734, 10.459646582243932] #Validated against a basic turbojet model
+        for (i,xi) in enumerate(x)
+            @test xi ≈ x_check[i] rtol = 1e-6
+        end
+    end
+end
\ No newline at end of file
diff --git a/test/unit_test_combustion.jl b/test/unit_test_combustion.jl
index 3587d2d..0026a19 100644
--- a/test/unit_test_combustion.jl
+++ b/test/unit_test_combustion.jl
@@ -1,11 +1,11 @@
 @testset "combustion" begin
     CH4 = species_in_spdict("CH4")
-    @test IdealGasThermo.LHV(CH4) == 5.002736488044851e7
-    @test IdealGasThermo.LHV("CH4") == 5.002736488044851e7
+    @test IdealGasThermo.LHV(CH4) ≈ 5.002736488044851e7 rtol = 1e-6
+    @test IdealGasThermo.LHV("CH4") ≈ 5.002736488044851e7 rtol = 1e-6
 
-    @test IdealGasThermo.AFT(CH4) == 2376.6102617357988
+    @test IdealGasThermo.AFT(CH4) ≈ 2376.6102617357988 rtol = 1e-6
     O2 = species_in_spdict("O2")
-    @test IdealGasThermo.AFT(CH4, O2) == 5280.53933225877
+    @test IdealGasThermo.AFT(CH4, O2) ≈ 5280.53933225877 rtol = 1e-6
 end
 @testset "burnt gas funcs." begin
     CH4 = species_in_spdict("CH4")
diff --git a/test/unit_test_turbo.jl b/test/unit_test_turbo.jl
index 474216b..9f05977 100644
--- a/test/unit_test_turbo.jl
+++ b/test/unit_test_turbo.jl
@@ -3,25 +3,25 @@
     gas.T = Tstd = IdealGasThermo.Tstd
     gas.P = Pstd = IdealGasThermo.Pstd
 
-    IdealGasThermo.compress(gas, 2.0, 1.0)
+    IdealGasThermo.compress!(gas, 2.0, 1.0)
     @test gas.P == 2 * Pstd
 
     #Test that compress throws right errors
     PR = 0.5
     err_msg = "The specified pressure ratio (PR) to compress by needs to be ≥ 1.0.
         Provided PR = $PR. Did you mean to use `expand`?"
-    @test_throws ErrorException(err_msg) IdealGasThermo.compress(gas, PR)
+    @test_throws ErrorException(err_msg) IdealGasThermo.compress!(gas, PR)
     PR = 1.5
     err_msg = "The specified pressure ratio (PR) to compress by needs to be ≤ 1.0.
         Provided PR = $PR. Did you mean to use `compress`?"
-    @test_throws ErrorException(err_msg) IdealGasThermo.expand(gas, PR)
+    @test_throws ErrorException(err_msg) IdealGasThermo.expand!(gas, PR)
 
     #Test gas mixing
     set_TP!(gas, Tstd, Pstd)
     gas2 = deepcopy(gas)
     set_TP!(gas2, 3 * Tstd, Pstd)
     gas3 = IdealGasThermo.gas_mixing(gas, gas2, 2.0)
-    @test gas3.T == 703.6767764998808
+    @test gas3.T ≈ 703.6767764998808 rtol = 1e-6
 
     @testset "tasopt comparisons" begin
         gas = Gas()
@@ -38,11 +38,11 @@
         @test gas.T ≈ 329.99 rtol = 1e-3
 
         set_TP!(gas, Tstd, Pstd)
-        IdealGasThermo.compress(gas, 2.0, 1.0)
+        IdealGasThermo.compress!(gas, 2.0, 1.0)
         @test gas.T ≈ 363.29 atol = 1e-1
 
         set_TP!(gas, Tstd, Pstd)
-        IdealGasThermo.expand(gas, 0.5, 1.0)
+        IdealGasThermo.expand!(gas, 0.5, 1.0)
         @test gas.T ≈ 244.547 atol = 1e-1
 
         set_TP!(gas, Tstd, Pstd)
diff --git a/test/unit_test_vitiated.jl b/test/unit_test_vitiated.jl
index 7168515..287d9ef 100644
--- a/test/unit_test_vitiated.jl
+++ b/test/unit_test_vitiated.jl
@@ -139,8 +139,8 @@ end
     set_TP!(gas, Tstd, Pstd)
     set_TP!(gas1, Tstd, Pstd)
     PR = 10.0
-    IdealGasThermo.compress(gas, PR, ηp)
-    IdealGasThermo.compress(gas1, PR, ηp)
+    IdealGasThermo.compress!(gas, PR, ηp)
+    IdealGasThermo.compress!(gas1, PR, ηp)
     @test gas.P ≈ Pstd * PR
     @test gas1.P ≈ Pstd * PR
     @test gas.h ≈ gas1.h
@@ -152,8 +152,8 @@ end
     set_TP!(gas, Tstd, Pstd)
     set_TP!(gas1, Tstd, Pstd)
     PR = 0.5
-    IdealGasThermo.expand(gas, PR, ηp)
-    IdealGasThermo.expand(gas1, PR, ηp)
+    IdealGasThermo.expand!(gas, PR, ηp)
+    IdealGasThermo.expand!(gas1, PR, ηp)
     @test gas.P ≈ Pstd * PR
     @test gas1.P ≈ Pstd * PR
     @test gas.T < Tstd