#if you want to figure out how many sigma tensions are between two Gaussian-error measurements, this is all what you need.
#by Zhiqi Huang (huangzhq25@mail.sysu.edu.cn)

import numpy as np
from scipy.special import gamma, erfc


def gamint(n, t):
    return t**(n/2-1)*np.exp(-t)

    
def tension_sigma(x1, x2, cov1, cov2):
    dim = len(x1)
    x = np.matrix(x1 - x2)
    cov = cov1+cov2
    invcov = cov **(-1)
    s = x * invcov * x.T
    t1d = np.sqrt(s[0, 0])
    if(dim == 1 or t1d > 20.):
        return t1d
    t = s[0,0]/2.
    dt = 1.e-3
    p = gamint(dim, t)/2.
    for i in range(6000):
        t += dt
        p += gamint(dim, t)
    p = p*dt
    t += dt/2.
    p += gamint(dim, t)
    p = p/gamma(dim/2.)
    t = t1d/np.sqrt(2.)
    p -= erfc(t)
    dt = 2.e-3
    p = p * np.sqrt(np.pi)/2./dt
    if( p == 0. ):
        return t*np.sqrt(2.)
    while(p < 0.):
        t += dt
        p += np.exp(-(t-dt/2.)**2)
    while p > 0.:
        t -= dt
        p -= np.exp(-(t+dt/2.)**2)
    pp = p + np.exp(-(t+dt/2.)**2)
    tt = t+dt
    return (tt*(-p)+t*pp)/(pp-p)*np.sqrt(2.)
    

x1 = np.array([180., 180.]) #measurement 1
x2 = np.array([140., 240.]) #measurement 2
cov1 = np.matrix([[100., 10.], [10., 100.]])  #covariance matrix for measurement 1
cov2 = np.matrix([[2., 0.], [0., 2.]]) #covariance matrix for measurement 2

print("there is a ", tension_sigma(x1, x2, cov1, cov2),"-sigma tension between measurement 1 and 2") #this is the number of sigmas for the tension between measurement 1 and 2