/*  impdiff.c
	implicit diffusion simulation

	J. Watlington, 10/16/95
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define DEBUG
#define DEBUG_PERIOD       100.0

#define NUM_POINTS         50
#define NUM_BUFFERS        2
#define MIN_TOTAL_CHANGE   0.0000001
#define MAX_TIME           2000.0
#define DIFF_COEFF         1.0
#define DELTA_X            1.0
#define DEFAULT_DELTA_T    1.0
#define START_IMPULSE_AMPL 1.0
#define START_IMPULSE_LOC  ( NUM_POINTS / 2 )

typedef long double sim_type;
void display_lattice( sim_type *lattice );

void main( int argc, char **argv )
{
	int          x, buf;
	int          cur_buf = 1;
	int          prev_buf = 0;

	long double  t = 0.0;
	long double  delta_t = DEFAULT_DELTA_T;
	long double  total_change = 10.0 * MIN_TOTAL_CHANGE;

	sim_type     prev_contrib;
	sim_type     self_component;
	sim_type     alpha, neg_alpha;

	sim_type  	 state[ NUM_BUFFERS ][ NUM_POINTS ];
	sim_type     temp_c[ NUM_POINTS ];
	sim_type     temp_d[ NUM_POINTS ];  /* not strictly necessary  */

	/* Parse arguments, if provided  */

	if( argc > 2 )
	{
		printf( "usage: %s [<time_step_size>]\n", argv[0] );
		exit( -1 );
	}
	else if( argc > 1 )
		sscanf( argv[1], "%Lf", &delta_t );

	printf( "running till total change less than %f\n", MIN_TOTAL_CHANGE );
	printf( "using a %f second time step\n", delta_t );

	/*  Initialize the simulation state to initial conditions   */

	for( buf = 0; buf < NUM_BUFFERS; buf++ )
		for( x = 0; x < NUM_POINTS; x++ )
			state[ buf ][ x ] = 0.0;

	state[ prev_buf ][ START_IMPULSE_LOC ] = START_IMPULSE_AMPL;

	alpha = (DIFF_COEFF * (sim_type)delta_t) / (DELTA_X * DELTA_X);
	neg_alpha = - alpha;
	self_component = 1.0 + 2.0 * alpha;

#ifdef DEBUG
	printf( "\nAt time t = %f, here is the state of the system :\n", t );
	display_lattice( state[ prev_buf ] ); 
#endif

	/*  And loop until the simulation stops changing...  */

	while( (total_change > MIN_TOTAL_CHANGE) && (t < MAX_TIME) )
	{
		/*  Init the array storing the interesting parts of the tri-diag mat */

		for( x = 1; x < NUM_POINTS - 1; x++ )
			temp_c[ x ] = neg_alpha;

		/*  Perform forward pass through tri-diagonal matrix.    */

		prev_contrib = 1.0 /( 1.0 + alpha );
		temp_c[ 0 ] = neg_alpha * prev_contrib;
		temp_d[ 0 ] = state[ prev_buf ][ 0 ] * prev_contrib;
		
		for( x = 1; x < NUM_POINTS - 1; x++ )
		{
			prev_contrib = 1.0 / (self_component + alpha * temp_c[ x - 1 ]);
			temp_c[ x ] = temp_c[ x ] * prev_contrib;
			temp_d[ x ] = (state[ prev_buf ][ x ] + alpha * temp_d[ x - 1 ])
								* prev_contrib;
		}

		/* never store temp_d[ NUM_POINTS - 1 ] (calc'ed below) to memory */
		prev_contrib = (state[ prev_buf ][ x ] + alpha * temp_d[ x - 1 ])
								  / (1.0 + alpha + alpha * temp_c[ x - 1 ]);
		
		/*  Perform reverse pass through tri-diagonal matrix.     */
		
		state[ cur_buf ][ NUM_POINTS - 1 ] = prev_contrib;
		total_change = 0.0;
		
		for( x = NUM_POINTS - 2; x >= 0; x-- )
		{
			prev_contrib = temp_d[ x ] - temp_c[ x ] * prev_contrib;
			state[ cur_buf ][ x ] = prev_contrib;

			total_change += fabs( prev_contrib - state[ prev_buf ][ x ] );
		}

		/*  Flip buffers   */
		if( ++cur_buf >= NUM_BUFFERS ) cur_buf = 0;
		if( ++prev_buf >= NUM_BUFFERS ) prev_buf = 0;
		
		t += delta_t;

#ifdef DEBUG
	    if( fmod( t, DEBUG_PERIOD ) == 0.0 )
		{
			printf( "\nAt time t = %f, here are the results :\n", t );
			display_lattice( state[ prev_buf ] ); 
		}
#endif
	}

	if( t >= MAX_TIME )
		printf( "\nThe solution never stabilized !!!\n" );
	else
		printf( "\nThe solution stabilized at time t = %f :\n", t );
	display_lattice( state[ prev_buf ] );
}


#define NUM_COLS   5
#define NUM_ROWS   ( NUM_POINTS / NUM_COLS )

void display_lattice( sim_type *lattice )
{
	int   x, y, index;
	
	for( index = 0, y = 0; y < NUM_ROWS; y++ )
	{
		printf( "%d:", (NUM_COLS * y) );
		for( x = 0; x < NUM_COLS; x++, index++ )
			printf( "  %1.5f", lattice[ index ] );
		printf( "\n" );
	}
}