/*
 *  =================================================================
 *  chimney.c : Finite difference computation for temperature 
 *              distribution in cross section of chimney.
 *
 *  Copyright (C) 1998 by Mark Austin and David Chancogne.
 *                                                                     
 *  This software is provided "as is" without express or implied warranty.
 *  Permission is granted to use this software on any computer system,
 *  and to redistribute it freely, subject to the following restrictions:
 * 
 *  1. The authors are not responsible for the consequences of use of
 *     this software, even if they arise from defects in the software.
 *  2. The origin of this software must not be misrepresented, either
 *     by explicit claim or by omission.
 *  3. Altered versions must be plainly marked as such, and must not
 *     be misrepresented as being the original software.
 *  4. This notice is to remain intact.
 *
 *  Written By : Mark Austin  
 *  =================================================================
 */

#include <stdio.h>

#define max(a, b) ((a) > (b) ? (a) : (b))

enum { Rows = 9, Cols = 9 };    /*  size of finite difference grid */
enum { InsideTemp = 200 };      /*  temperature inside chimney     */

int main( void ) {
float faaTemp[ Rows ][ Cols ];
float fTempSum, fMaxChange;
int iIndexi, iIndexj;

   /* [a] : Initialize temperatures. */

   for( iIndexi = 0; iIndexi < 5; iIndexi = iIndexi + 1 ) {
      faaTemp[ 4  ][ iIndexi ] = InsideTemp;
      faaTemp[ iIndexi ][  4 ] = InsideTemp;
   }

   for( iIndexi = 5; iIndexi < Rows; iIndexi = iIndexi + 1 ) {
      for( iIndexj = 0; iIndexj < iIndexi; iIndexj = iIndexj + 1 ) {
         faaTemp[ iIndexj ][ iIndexi ] = 0.0;
         faaTemp[ iIndexi ][ iIndexj ] = 0.0;
      }
   }

   /* [b] : Compute temperatures over interior points. */

   fMaxChange = InsideTemp;
   while( fMaxChange > 0.5 ) {
      fMaxChange = 0.0;

      for( iIndexi = 5; iIndexi <= 7; iIndexi = iIndexi + 1 ) {

         /* [b.1] : Temperature along edges 3 and 6 */

         fTempSum = faaTemp[ iIndexi+1 ][0] + faaTemp[ iIndexi-1 ][0] +
                    2*faaTemp[ iIndexi ][1];
         fMaxChange = max( fMaxChange,
                           fabs(faaTemp[ iIndexi ][ 0 ] - fTempSum/4.0) );

         faaTemp[ iIndexi ][ 0 ] = faaTemp[0][ iIndexi ] = fTempSum/4.0;

         /* [b.2] : Temperature inside region */

         for( iIndexj = 1; iIndexj <= iIndexi; iIndexj = iIndexj + 1 ) {
            fTempSum = faaTemp[ iIndexi+1 ][ iIndexj ] +
                       faaTemp[ iIndexi-1 ][ iIndexj ] +
                       faaTemp[ iIndexi ][ iIndexj+1 ] +
                       faaTemp[ iIndexi ][ iIndexj- 1];

            fMaxChange = max( fMaxChange,
                              fabs( faaTemp[iIndexi][iIndexj] -
                                    fTempSum/4.0) );
 
            faaTemp[ iIndexj ][ iIndexi ] = fTempSum/4.0;
            faaTemp[ iIndexi ][ iIndexj ] = fTempSum/4.0;
         }
      }
   }

   /* [c] : Print matrix of temperatures */

   printf(" temperature = [ \n");
   for( iIndexi = 0; iIndexi < Cols; iIndexi = iIndexi + 1 ) {
      for( iIndexj = 0; iIndexj < Rows; iIndexj = iIndexj + 1 ) {
          if( iIndexi < 4 && iIndexj < 4 )
              printf("   NaN ");
          else
              printf(" %5.1f ", faaTemp[ iIndexi ][ iIndexj ]);
      }
      printf(";\n");
   }
   printf(" ]; \n");
}