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

#define DEBUG 0
#define THRESHOLD (5e-3)

void init_phi( double * phi, int n );
void update( double *cur, double *next, int n );
int converged( double *cur, double *next, int n );

int main( int argc, char *argv[] ) 
{
    int n, niters, conv;
    double *phi_cur, *phi_next, *tmp;

    if (argc != 2)
    {
        fprintf( stderr, "Usage: %s <n>\n", argv[0] );
        exit( -1 );
    }

    n = atoi( argv[1] );

    phi_cur  = (double *) malloc ( n * n * sizeof(double) );
    phi_next = (double *) malloc ( n * n * sizeof(double) );
    
    init_phi(phi_cur, n);
    init_phi(phi_next, n);

    niters = 0;
    while( 1 )
    {
        niters++;
#if DEBUG
        printf("Iteration %d\n", niters );
        int i,j;
        for ( i = 0; i < n; i++ )
        {
            printf("[" );
            for ( j = 0; j < n; j++ )
                printf(" %10.6f ", phi_cur[j*n + i] );
            printf("]\n" );
        }
        sleep(1);
#endif
        // Compute next (new) phi from current (old) phi
        update( phi_cur, phi_next, n );
        // If converged, we are done
        conv = converged( phi_cur, phi_next, n );
        if ( conv )
            break;
        // Otherwise, swap pointers and continue
        tmp = phi_cur;
        phi_cur = phi_next;
        phi_next = tmp;
    }

    free( phi_cur );
    free( phi_next );

    printf( "Converged after %d iterations\n", niters );
}


void init_phi( double * phi, int n )
{
    int i, j;

    // Interior points initialized to 50 degrees
    for ( i = 1; i < n-1; i++ )
        for ( j = 1; j < n-1; j++ )
            phi[j*n+i] = 50.0;

    // Top, left, and right boundaries fixed at 100 degrees
    for ( i = 0; i < n; i++ )
    {
        phi[    0*n +    i] = 100.0;
        phi[(n-1)*n +    i] = 100.0;
        phi[    i*n +    0] = 100.0;
    }
    // Bottom boundary fixed at 0 degrees
    for ( i = 0; i < n; i++ )
        phi[    i*n +(n-1)] = 0.0;
}

void update( double *cur, double *next, int n ) 
{
    int i, j;

    for( j = 1; j < n-1; j++ )
        for( i = 1; i < n-1; i++ )
            next[ j*n+i ] = 
                ( cur[ (j-1)*n+i ] +
                  cur[ j*n+(i-1) ] +
                  cur[ j*n+(i+1) ] +
                  cur[ (j+1)*n+i ] ) / 4;
}

int converged( double *cur, double *next, int n ) 
{
    int i, j;

    for( j = 1; j < n-1; j++ )
        for( i = 1; i < n-1; i++ )
            if( fabs( next[ j*n + i ] - cur[ j*n + i ] ) > THRESHOLD )
                return 0;
    return 1;
}