anchor.diff:  
   
   
  -  function Eigfunc(Ak,alpha,lamc,xpos)
  -  ! Used in KLrxi_point to call reconstructed realization at given point
  +  function Eigfunc(Ak,alpha,lamc,xpos,orderin)
  +  ! Evaluates Eigenfunction term for any order (0+) derivative of the eigenfunction
     real(8) :: Ak,alpha,lamc,xpos,Eigfunc
  -
  -  Eigfunc = Ak * ( sin(alpha*xpos) + lamc*alpha*cos(alpha*xpos) )
  +  integer, optional :: orderin
  +  integer :: order
  +
  +  order = merge(orderin,0,present(orderin))
  +  select case (mod(order,4))
  +    case (0)
  +      Eigfunc = Ak * ( sin(alpha*xpos) + lamc*alpha*cos(alpha*xpos) )
  +    case (1)
  +      Eigfunc = Ak * (-cos(alpha*xpos) + lamc*alpha*sin(alpha*xpos) )
  +    case (2)
  +      Eigfunc = Ak * (-sin(alpha*xpos) - lamc*alpha*cos(alpha*xpos) )
  +    case (3)
  +      Eigfunc = Ak * ( cos(alpha*xpos) - lamc*alpha*sin(alpha*xpos) )
  +  end select
  +  Eigfunc = Eigfunc/(alpha**order)
   
     end function Eigfunc
  +
  +
  +
  +  subroutine KLr_setmeans(order,tmeanadjust,tsigsmeanadjust,tsigameanadjust,tsigave,tsigscatave,tsigabsave)
  +  !This subroutine sets values for non-x-dependent terms based on derivative order
  +  use genRealzvars, only: sigave, sigscatave, sigabsave
  +  use KLvars, only: meanadjust, sigsmeanadjust, sigameanadjust
  +
  +  integer :: order
  +  real(8) :: tmeanadjust,tsigsmeanadjust,tsigameanadjust,tsigave,tsigscatave,tsigabsave
  +  if(order==0) then
  +    tmeanadjust     = meanadjust
  +    tsigsmeanadjust = sigsmeanadjust
  +    tsigameanadjust = sigameanadjust
  +    tsigave         = sigave
  +    tsigscatave     = sigscatave
  +    tsigabsave      = sigabsave
  +  else
  +    tmeanadjust     = 0
  +    tsigsmeanadjust = 0
  +    tsigameanadjust = 0
  +    tsigave         = 0
  +    tsigscatave     = 0
  +    tsigabsave      = 0
  +  endif
  +
  +  end subroutine KLr_setmeans
  +
   
   
   end module KLreconstruct
  
positive.code:   elseif(flmatbasedxs) then
      !cross section values
      if(chxstype .ne. 'scatter') &
        siga = sigabsave  + sigameanadjust + sqrt(Coabs)  * KL_sum
      if(chxstype .ne. 'absorb') &
        sigs = sigscatave + sigsmeanadjust + sqrt(Coscat) * KL_sum
  
      !determine point value
      select case (chxstype)
        case ("total")
  
negative.code:   real(8) :: KL_point
    character(*) :: chxstype
    integer, optional :: tnumEigsin
    integer, optional :: orderin
  
    real(8) :: sigt, siga, sigs, KL_sum, Eigfterm
    integer :: curEig,tnumEigs,order
    real(8) :: tmeanadjust,tsigsmeanadjust,tsigameanadjust,tsigave,tsigscatave,tsigabsave
  
    !load any optional values or their default
    tnumEigs = merge(tnumEigsin,numEigs,present(tnumEigsin))
    order    = merge(orderin   ,0      ,present(orderin)   )
  
    !solve summation of KL terms to use below
    KL_sum = 0d0
    do curEig=1,tnumEigs
      Eigfterm = Eigfunc(Ak(curEig),alpha(curEig),lamc,xpos,order)
      KL_sum   = KL_sum + sqrt(Eig(curEig)) * Eigfterm * KLrxivals(j,curEig)
    enddo
  
    !set non-x-dependent values based order
    call KLr_setmeans(order,tmeanadjust,tsigsmeanadjust,tsigameanadjust,tsigave,tsigscatave,tsigabsave)
  
    !solve value at point
    if(.not.flmatbasedxs) then
      !determine underlying value
      sigt = tsigave + tmeanadjust + KL_sum
      !determine point value
      select case (chxstype)
        case ("total")
  