Output: The Jacobian of the cost functional
n = len(u0b)
#determine the assimilation window
t = t[:ind_m[-1]+1] #cut the time till the last observatino point
nt = len(t)-1
dt = t[1] - t[0]
ub = np.zeros([n,nt+1]) #base trajectory
lam = np.zeros([n,nt+1]) #lambda sequence
fk = np.zeros([n,len(ind_m)])
Ri = np.linalg.inv(R)
ub[:,0] = u0b
if opt == 0: #Euler
#forward model
for k in range(nt):
ub[:,k+1] = euler(rhs,ub[:,k],dt,*args)
#backward adjoint
k = ind_m[-1]
fk[:,-1] = (JObsOp(ub[:,k])).T @ Ri @ (w[:,-1]-ObsOp(ub[:,k]))
lam[:,k] = fk[:,-1] #lambda_N=f_N
km = len(ind_m)-2
for k in range(ind_m[-1],0,-1):
DM = Jeuler(rhs,Jrhs,ub[:,k-1],dt,*args)
lam[:,k-1] = (DM).T @ lam[:,k]
if k-1 == ind_m[km]:
fk[:,km] =(JObsOp(ub[:,k-1])).T @ Ri @ (w[:,km]-ObsOp(ub[:,k-1]))
lam[:,k-1] = lam[:,k-1] + fk[:,km]
km = km - 1
elif opt == 1: #RK4