#include <iostream>
#include <cstdlib>
#include <fstream>

using namespace std;

int N; // size of the array


void initBoundaryCondition(double** a, double north, double south, double east, double west ){
    for (int i = 0; i < N; ++i){
        a[i][0]   = west;
        a[N-1][0] = east;
        a[0][i]   = north;
        a[N-1][i] = south;
    }
}

void collide(double** main, double** tmp){
    for (int iX = 1; iX < N-1; ++iX){
        for (int iY = 1; iY < N-1; ++iY){
            main[iX][iY] = (tmp[iX+1][iY]+tmp[iX-1][iY]+tmp[iX][iY+1]+tmp[iX][iY-1])/4.0;
        }
    }
}

void propagate(double** main, double** tmp){
    for (int iX = 1; iX < N-1; ++iX){
        for (int iY = 1; iY < N-1; ++iY){
            tmp[iX][iY] = main[iX][iY];
        }
    }
}

void writeResult(double** mat){
    ofstream out("result.dat");
    for (int iX = 1; iX < N-1; ++iX){
        for (int iY = 1; iY < N-1; ++iY){
            out << mat[iX][iY] << " ";
        }
        out << endl;
    }
    out << endl;
}

int main(int argc, char** argv){
    // the size of the arrays
    N = 400;
    
    // iteration number
    if (argc != 2){
        cout << "Usage: " << argv[0] << " iterations" << endl;
        return -1;
    }
    
    int maxIt = atoi(argv[1]);
    
    // the two arrays
    double** mainArray;
    double** tmpArray;
    
    // initialize
    mainArray = new double*[N];
    tmpArray = new double*[N];
    for (int i = 0; i < N; ++i){
        mainArray[i] = new double[N];
        tmpArray[i] = new double[N];
    }
    
    initBoundaryCondition(mainArray,1,0,0,1);
    initBoundaryCondition(tmpArray,1,0,0,1);

    // iteration number
    for (int iT = 0; iT < maxIt; ++iT){
        collide(mainArray,tmpArray);
        propagate(mainArray,tmpArray);
    }
    
    writeResult(mainArray);
    
    delete [] mainArray;
    delete [] tmpArray;
}