"""
    intensity_sppa(pressures::Array{Unitful.Pressure, N}, medium::Medium, excitation::Excitation)

Spatial Peak, Pulse Averaged intensity [W/cm²]

`pressures` is expected to be an array with up to 4 dimensions. The first dimension is time,
and the remaining dimensions are space
"""
function intensity_sppa(pressures::Array{T, N}, medium, excitation) where {T <: Unitful.Pressure, N}
    # account for sampling frequency by multiplying pressure values with dt
    dt = 1/excitation.fs
    intensities = map(p -> intensity(p, medium), pressures)
    i_sppa = sum(intensities * dt; dims=1) / excitation.pulse_duration
    return format_intensity.(squeeze(i_sppa))
end


"""
    intensity(pressure::Unitful.Pressure, medium::Medium)

Instantaneous intensity [W/cm²]

`pressure` is expected to be a single number
"""
function intensity(pressure::Unitful.Pressure, medium::Medium)
    return format_intensity(pressure^2 / (medium.density * medium.c))
end


format_intensity(input) = uconvert(u"W/cm^2", input)


"""
    intensity_spta(pressures::Array{Unitful.Pressure, N}, medium::Medium, excitation::Excitation)

Spatial Peak, Time Averaged intensity [W/cm²]

`pressures` is expected to be an array with up to 4 dimensions. The first dimension is time,
and the remaining dimensions are space
"""
function intensity_spta(pressures::Array{T, N}, medium::Medium, excitation::Excitation) where {T <: Unitful.Pressure, N}
    return intensity_sppa(pressures, medium, excitation) * excitation.duty_cycle
end


"""
    energy_spta(pressures::Array{Unitful.Pressure, N}, medium::Medium, excitation::Excitation)

Spatial Peak, Time Averaged energy [J/cm²]
"""
function energy_spta(pressures::Array{T, N}, medium::Medium, excitation::Excitation) where {T <: Unitful.Pressure, N}
    return intensity_spta(pressures, medium, excitation) * excitation.total_duration
end


"""
    mechanical_index(pressures::Array{Unitful.Pressure, N}, excitation::Excitation)

Mechanical index [dimensionless]

`pressures` is expected to be an array with up to 4 dimensions. The first dimension is time,
and the remaining dimensions are space
"""
function mechanical_index(pressures::Array{T, N}, excitation::Excitation) where {T <: Unitful.Pressure, N}
    peak_rarefaction_pressure = pn_pressure(uconvert.(u"MPa", pressures))
    center_frequency = sqrt(uconvert(u"MHz", excitation.f0))
    return ustrip(peak_rarefaction_pressure / center_frequency)
end


"""
    pn_pressure(pressures::Array{Unitful.Pressure, N})

Peak Negative pressure [units unchanged]

`pressures` is expected to be an array with up to 4 dimensions. The first dimension is time,
and the remaining dimensions are space
"""
function pn_pressure(pressures::Array{T, N}) where {T <: Unitful.Pressure, N}
    peak_neg_pressure = abs.(minimum(pressures; dims=1))
    return squeeze(peak_neg_pressure)
end


"""
    pp_pressure(pressures::Array{Unitful.Pressure, N})

Peak Positive pressure [units unchanged]

`pressures` is expected to be an array with up to 4 dimensions. The first dimension is time,
and the remaining dimensions are space
"""
function pp_pressure(pressures::Array{T, N}) where {T <: Unitful.Pressure, N}
    peak_pos_pressure = abs.(maximum(pressures; dims=1))
    return squeeze(peak_pos_pressure)
end


"""
    ptp_pressure(pressures::Array{Unitful.Pressure, N})

Peak To Peak pressure [units unchanged]

`pressures` is expected to be an array with up to 4 dimensions. The first dimension is time,
and the remaining dimensions are space
"""
function ptp_pressure(pressures::Array{T, N}) where {T <: Unitful.Pressure, N}
    ptp_pressure = pp_pressure(pressures) .+ pn_pressure(pressures)
    return ptp_pressure
end
