;+
; NAME:
;  vt3d_cn_tridc
; PURPOSE: (one line)
;  Upper-Lower decomposition of the substrate matrix, S''
; DESCRIPTION:
;  Upper-Lower decomposition of the substrate matrix, S''
; CATEGORY:
;  Volatile Transport
; CALLING SEQUENCE:
;  spp_tridc = vt3d_cn_tridc(alpha_i, beta_i, spp_indx)
; INPUTS:
;  alpha_i: dependance of temp[1..J] on temp[0..J-1], unitless, float[n_j]
;  beta_i:  dependance of temp[1..J] on temp[2..J+1], unitless, float[n_j]
; OUTPUTS:
;   spp_tridc - the tridiagonal LU decomposition of the substrate matrix, S.
;    spp_tridc[0,0:n_j-2] = spp_al 
;      lower bidiagonal array.  float[n_j]
;    spp_tridc[1,0:n_j-1] = spp_a  
;      diagonal elements of the upper array. float[n_j+1]
;    spp_tridc[2,0:n_j-2] = spp_au 
;      superdiagonal elements of the upper array. float[n_j]
;    spp_tridc[3,0:n_j-3] = u2 
;      Second superdiagonal of the upper array. float(n_j-1)
;   spp_index An output vector that records the row permutations which occurred as a result of partial pivoting. 
; PROCEDURE:
;  Solve Crank Nickleson implicit  equation
;  Young, L. A. 2016, Volatile transport on inhomogeneous surfaces: II. Numerical calculations (VT3D)
;   Resubmitted to Icarus.
;   This is used as part of the Crank-Nicholson time step.  See:
;   Table 4 & 5, Figure 3-6, Figure 3-8
;  The LU decomposition is done using LA_TRIDC
; MODIFICATION HISTORY:
;  Written 2011 Feb 6, by Leslie Young, SwRI
;  2016 Mar 25 LAY. Modified for inclusion in vt3d library.
;-

function vt3d_cn_tridc, alpha_i, beta_i, spp_index

    n_j = n_elements(alpha_i)

    lower = -[alpha_i[1:n_j-1]]/2
    diag = 1 + (alpha_i + beta_i)/2
    upper = -beta_i[0:n_j-2]/2
    LA_TRIDC, lower, diag, upper, u2, spp_index
    spp_tridc = dblarr(4,n_j)
    spp_tridc[0,0:n_j-2] = lower 
    spp_tridc[1,0:n_j-1] = diag  
    spp_tridc[2,0:n_j-2] = upper
    spp_tridc[3,0:n_j-3] = u2

    return, spp_tridc
    

end