Cooldt Code

/*
 * Cooldt cool-down time calculator for multi-layers pipe
 *
 * Copyright 2010, 2013, 2014, Benjamin DEGLO DE BESSES. All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without modification, are
 * permitted provided that the following conditions are met:
 *
 *    1. Redistributions of source code must retain the above copyright notice, this list of
 *       conditions and the following disclaimer.
 *
 *    2. Redistributions in binary form must reproduce the above copyright notice, this list
 *       of conditions and the following disclaimer in the documentation and/or other materials
 *       provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY Benjamin DEGLO DE BESSES "AS IS" AND ANY EXPRESS OR IMPLIED
 * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
 * FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL Benjamin DEGLO DE BESSES OR
 * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
 * ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
 * ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 *
 */

#include "cooldt.h"

/*
 * Solve a tridiagonal system with LU factorization (AKA Thomas algorithm)
 * !! this implementation will modify the arrays C and D
 * see http://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm for details
 * here it doesn't matter because A,B,C & D are recalculated by
 * ImplicitSystemCoefs before every TriDiagonalSolve calls
 */
void TriDiagonalSolve(double *A, double *B, double *C, double *D, int I_inf, int I_sup, double *X)
{
    int i ;

    /*
     * Step 1 : Triangulation
     */
    C[I_inf] /= B[I_inf] ;
    D[I_inf] /= B[I_inf] ;

    for( i = (I_inf + 1) ; i <= I_sup ; i++ )
    {
        double beta = 1.0 / (B[i] - C[i-1] * A[i]);
        C[i] *= beta;
        D[i] = (D[i] - D[i-1] * A[i]) * beta ;
    }

    /*
     * Step 2 : Backsubstitution
     */
    X[I_sup] = D[I_sup];
    for( i = (I_sup - 1) ; i >= I_inf ; i-- )
        X[i] = D[i] - C[i] * X[i + 1];

}

/* End of algebra.c */