void setInitialConditions(Array<float,2>& c, Array<float,2>& P1,
    Array<float,2>& P2, Array<float,2>& P3, int N)
{
    // Set the velocity field
    c = 0.2;

    // Solid block with which the pulse collides
    int blockLeft = 0;
    int blockRight = 2*N/5.0-1;
    int blockTop = N/3-1;
    int blockBottom = 2*N/3.0-1;
    c(Range(blockTop,blockBottom),Range(blockLeft,blockRight)) = 0.5;

    // Channel directing the pulse leftwards
    int channelLeft = 4*N/5.0-1;
    int channelRight = N-1;
    int channel1Height = 3*N/8.0-1;
    int channel2Height = 5*N/8.0-1;
    c(channel1Height,Range(channelLeft,channelRight)) = 0.0;
    c(channel2Height,Range(channelLeft,channelRight)) = 0.0;

    // Initial pressure distribution: gaussian pulse inside the channel
    using namespace blitz::tensor;
    int cr = N/2-1;
    int cc = 7.0*N/8.0-1;
    float s2 = 64.0 * 9.0 / pow2(N/2.0);
    P1 = 0.0;
    P2 = exp(-(pow2(i-cr)+pow2(j-cc)) * s2);
    P3 = 0.0;
}