#!/usr/bin/env python
# coding: utf-8

# # Gaussian Error Propagation with SymPy

# Klaus Reygers, 2020

# In[19]:


from sympy import *
from IPython.display import display, Latex


# In[20]:


def gaussian_error_propagation(f, vars):
    """
    f: formula (sympy expression)
    vars: list of independent variables and corresponding uncertainties 
    [(x1, sigma_x1), (x2, sigma_x2), ...]
    """
    sum = sympify("0") # empty sympy expression
    for (x, sigma) in vars:
        sum += diff(f, x)**2 * sigma**2 
    return sqrt(simplify(sum))


# Show usage for a simple example: Volume of a cylinder with radius $r$ and height $h$:

# In[21]:


r, h, sigma_r, sigma_h = symbols('r, h, sigma_r, sigma_h', positive=True)
V = pi * r**2 * h # volume of a cylinder


# In[22]:


sigma_V = gaussian_error_propagation(V, [(r, sigma_r), (h, sigma_h)])
display(Latex(f"$V = {latex(V)}, \, \sigma_V = {latex(sigma_V)}$"))


# Plug in some numbers and print the calculated volume with its uncertaity:

# In[23]:


r_meas = 3 # cm
sigma_r_meas = 0.1 # cm
h_meas = 5 # cm
sigma_h_meas = 0.1 # cm


# In[24]:


central_value = V.subs([(r,r_meas), (h, h_meas)]).evalf()
sigma = sigma_V.subs([(r, r_meas), (sigma_r, sigma_r_meas), (h, h_meas), (sigma_h, sigma_h_meas)]).evalf()
display(Latex(f"$$V = ({central_value:0.1f} \pm {sigma:.1f}) \, \mathrm{{cm}}^3$$"))