What is going on with the Lagrange multipliers in LsSensorModel

Link to the Jupyter notebook on gist.github.com (which is easier to read):

Lagrangian Interpolation

Given a set of points (t_i, y_i), Lagrangian interpolation computes a polynomial which attempts to approximate the function y(t):

y(t) \approx L(t) = \sum_{j=0}^k y_j*L_j(t)

The L_j(t) functions are polynomials in their own right (called the fundamental polynomials), and expressed by the product formula:

L_j(t) = \prod_{j\neq i}\frac{t - t_j}{t_i - t_j}

Consider a specific example where t=[1,2,3,4], y =[2,4,6,8] (y = 2t): In lagrangeInterp(....) this corresponds to the order=4 case. Just focus on the numbers in the denominator and based upon the vector t that are hard-coded within the function:

L_0(t) = \frac{numer(t)}{(1-2)(1-3)(1-4)} = \frac{num(t)}{-6}
L_1(t) = \frac{numer(t)}{(2-1)(2-3)(2-4)} = \frac{num(t)}{2}
L_2(t) = \frac{numer(t)}{(3-1)(3-2)(3-4)} = \frac{num(t)}{-2}
L_3(t) = \frac{numer(t)}{(4-1)(4-2)(4-3)} = \frac{num(t)}{6}

Below is Python code which implements the same method (except that it is hopefully easier to understand):

---

import numpy as np

def weight(x,j):
    n = len(x)
    xj=x[j]
    w=1;
    a = [k for k in range(n) if k!=j]
    for i in a:
        w=w*(xj-x[i])
    print('j=%d'% j)
    print(w)
    if (abs(w) < 0.00000000000000001):
        return 0;
    else:
        return (1.0/w)

def compute_weights(x):
    n = len(x)
    weights=[]
    for i in range(n):
        weights.append(weight(x,i))
    return weights

def lagrange(xi,x,y,w):
    n = len(w)
    lx=1
    wx=[]
    for i in range(n):
        s=xi-x[i]
        lx=lx*s
        wx.append((w[i]/s)*y[i])

    return lx*sum(wx)

y=[1,4,9,16,25,36,49,64]
x=[1,2,3,4,5,6,7,8]

w=compute_weights(x)
                                                
x0=9.5;
y0=lagrange(x0,x,y,w)
print y0

---

This scheme only works in the case where the distance between the t_i are uniform and unchanging. If that isn't true then the approximation will not be valid. It's being called in the UsgsAstroLsSenorModel::getAdjSensorPosVel(...) function to approximate position and velocity given a vector of time points, which in turn is called by getSensorVelocity(...), losToECF(...), getSensorPosition(...). I think it's probably ok, although I would have went with cubic splines myself. There is also no maximum error calculation happening in the code, but that can be easily rectified.

In the code, the t_i are time-points. t is an a arbitrary time-point, and the uniform interval is chopped into pieces of length: (time-startTime)/dt.

In UsgsAstroLsStateData:

• dt = m_DtEphem
• startTime = m_T0Ephem

I think m_DtEphem is the exposure time for the image.