################################################
#
# Solve the Second Order Differential Equation
#
# array
#         z: Vertical Coordinate, equi-distant, at 0 to K level
#                     0   : lower boundary, included
#                     K+1 : upper boundary, included
#         Q: Heating, length K+1
#         J: Forcing, length K+1
#
# boundary condition
#         wb : lower, fixed
#         ab : lower, dw/dz + ab w = bb
#         bb : lower, dw/dz + ab w = bb
#         
#         wt : upper, fixed
#         ab : upper, dw/dz + ab w = bb
#         bb : upper, dw/dz + ab w = bb
#
################################################
SolveSecondOrder <- function( z, Q, J,
                             wb=NULL, ab=NULL, bb=NULL,
                             wt=NULL, at=NULL, bt=NULL){
  delta <- z[2]-z[1]
  K     <- length(z)

  Qv    <- Q(z)
  Jv    <- J(z)
  a     <- delta^2 * Qv - 2
  b     <- delta^2 * Jv
  
  alph  <- a * 0
  beta  <- a * 0
  wres  <- a * 0
                                        # Lower Boundary
  if( is.null( wb )){
    alph[1] <- a[1]/2 + delta * ab
    beta[1] <- b[1]/2 + delta * bb
    alph[2] <- a[2] - 1       / alph[1]
    beta[2] <- b[2] - beta[1] / alph[1]
  }else if( is.null( ab )){
    alph[2] <- a[2]
    beta[2] <- b[2] - wb
  }else{
    print("Wrong Boundary Condition")
  }
                                        # alpha, beta, upward sweep
  for( k in 3:K ){
    alph[k] <- a[k] - (1        /alph[k-1])
    beta[k] <- b[k] - (beta[k-1]/alph[k-1])
  }
                                        # Upper Boundary
  if( is.null( wt )){
    wres[K] <- ( beta[K-1] - alph[K-1] * ( b[K]/2 - delta * bt ) ) / ( 1 - alph[K-1] * ( a[K]/2 - delta * at ) )
  }else if( is.null( at )){
    wres[K] <- ( beta[K] - wt ) / alph[K]
  }else{
    print("Wrong Boundary Condition")
  }
                                        # Downward Sweep
  for( k in seq(K-1,2,-1) ){
    wres[k] = ( beta[k] - wres[k+1] ) / alph[k]
  }
  if( is.null( wb )){
    wres[1] = ( beta[1] - wres[2] ) / alph[1]
  }else if( is.null( ab )){
    wres[1] <- wb
  }

  return(wres)
}

#################################
#
#  Example:
#
#    Value at the bottom :
#        sin(0.5) = 0.479425538604
#    Gradient at the bottom :
#        cos(0.5) = 0.87758256189
#
#    Note that "ab" is -1 times the negative value of the gradient
#################################
Q <- function(z){ return(rep(1, length(z)))}
J <- function(z){ return(rep(0, length(z)))}
z <- seq(0, 10, length=10001)
data <- SolveSecondOrder2(z, Q, J, wb= 0.479425538604, ab=0, bb=0.87758256189)
plot(z, data, pch='.')
lines(z, sin(z+.5), col='red')