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

#define UE 1
#define rou 0.002377
#define miu 3.737e-7
#define dx 0.025
#define dy 0.001
#define dt 0.001
#define sor 0.1

int i, j, k, steps=0, steps_sum;
double p[22][12] = { 0 }, u[23][12] = { 0 }, v[24][13] = { 0 };
double mass_u[23][12] = { 0 }, mass_v[24][13] = { 0 };

void show()
{
    printf("\n************************************%d************************************\n",steps);
    printf("\nu\n");
    for (j = 11; j >=1; j--)
    {
        for (i = 1; i <= 22; i++)
            printf("%.3lf ", u[i][j]);
        printf("\n");
    }        
    printf("\nv\n");
    for (j = 12; j >= 1; j--)
    {
        for (i = 1; i <= 23; i++)
            printf("%.3lf ", v[i][j]);
        printf("\n");
    }
    printf("\np\n");
    for (j = 11; j >= 1; j--)
    {
        for (i = 1; i <= 21; i++)
            printf("%.3lf ", v[i][j]);
        printf("\n");
    }
/*    printf("\nmass_u\n");
    for (j = 11; j >= 1; j--)
    {
        for (i = 1; i <= 22; i++)
            printf("%.3lf ", mass_u[i][j]);
        printf("\n");
    }
    printf("\nmass_v\n");
    for (j = 11; j >= 1; j--)
    {
        for (i = 1; i <= 22; i++)
            printf("%.3lf ", mass_v[i][j]);
        printf("\n");
    }*/
}

void initialize()
{
    for (i = 1; i <= 22; i++)
        u[i][11] = UE;
    v[15][5] = 0; //disturbancr
    for (i = 1; i <= 22; i++)
        mass_u[i][11] = rou*UE;
    mass_v[15][5] = rou*v[15][5];
}

void compute_mass()
{
    double A, B, v_1, v_2, u_1, u_2;
        for (i = 2; i <= 21; i++)
            for (j = 2; j <= 10; j++)
            {
                v_1 = 0.5*(v[i][j + 1] + v[i + 1][j + 1]);
                v_2 = 0.5*(v[i][j] + v[i + 1][j]);
                A = -((rou*u[i + 1][j] * u[i + 1][j] - rou*u[i - 1][j] * u[i - 1][j]) / (2 * dx) + (rou*u[i][j + 1] * v_1 - rou*u[i][j - 1] * v_2) / (2 * dy))
                    + miu*((u[i + 1][j] - 2 * u[i][j] + u[i - 1][j]) / (dx*dx) + (u[i][j + 1] - 2 * u[i][j] + u[i][j - 1]) / (dy*dy));
                mass_u[i][j] = rou*u[i][j] + A*dt - (dt / dx)*(p[i][j] - p[i - 1][j]);
            }
        for (i = 2; i <= 22; i++)
            for (j = 2; j <= 11; j++)
            {
                u_1 = 0.5*(u[i][j - 1] + u[i][j]);
                u_2 = 0.5*(u[i - 1][j - 1] + u[i - 1][j]);
                B = -((rou*v[i + 1][j] * u_1 - rou*v[i - 1][j] * u_2) / (2 * dx) + (rou*v[i][j + 1] * v[i][j + 1] - rou*v[i][j - 1] * v[i][j - 1]) / (2 * dy))
                    + miu*((v[i + 1][j] - 2 * v[i][j] + v[i - 1][j]) / (dx*dx) + (v[i][j + 1] - 2 * v[i][j] + v[i][j - 1]) / (dy*dy));
                mass_v[i][j] = rou*v[i][j] + B*dt - (dt / dy)*(p[i - 1][j] - p[i - 1][j - 1]);
            }
}

void pressure_correct()
{
    double a, b, c;
    double d[22][12] = { 0 };
    double _p[22][12] = { 0 };//pressure correct
    double _p_;//mid
    a = 2 * (dt / (dx*dx) + dt / (dy*dy));
    b = -dt / (dx*dx);
    c = -dt / (dy*dy);
    for (i = 2; i <= 20; i++)
        for (j = 2; j <= 10; j++)
            d[i][j] = (mass_u[i + 1][j] - mass_u[i][j]) / dx + (mass_v[i + 1][i + 1] - mass_v[i + 1][j]) / dy;    
    for (k = 1; k <= 100; k++)
    {
        for (i = 2; i <= 20; i++)
            for (j = 2; j <= 10; j++)
            {
                _p_ = -(b*_p[i + 1][j] + b*_p[i - 1][j] + c*_p[i][j + 1] + c*_p[i][j - 1] + d[i][j]) / a;
                _p[i][j] = _p[i][j] + 0.1*(_p_ - _p[i][j]);
            }

        
    }
    for (i = 2; i <= 20; i++)
        for (j = 2; j <= 20; j++)
            p[i][j] = p[i][j]+(sor*_p[i][j]);
}

void compute_velocity()
{
    for (i = 2; i <= 21; i++)
        for (j = 2; j <= 10; j++)
            u[i][j] = mass_u[i][j] / rou;
    for (i = 2; i <= 22; i++)
        for (j = 2; j <= 11; j++)
            v[i][j] = mass_v[i][j] / rou;
    for (j = 1; j <= 11; j++)
    {
        u[1][j] = u[2][j]; u[22][j] = u[21][j];
    }
    for (j = 1; j <= 12; j++)
        v[23][j] = v[22][j];
}

void main()
{
    double store[6][12] = { 0 };
    int m = 1;
    initialize();
    show();
    steps_sum = 20;
    for (steps = 1; steps <= steps_sum; steps++)
    {
        compute_mass();
        compute_velocity();        
        pressure_correct(); show();
        if (steps == 4 || steps == 20 || steps == 50 || steps == 150 || steps == 300)
        {
            show();
            for (j = 1; j <= 11; j++)
                store[m][j] = u[15][j];
            m += 1;
        }
    }
    printf("\n");
    for (j = 11; j >= 1; j--)
    {
        for (m = 1; m <= 5; m++)
            printf("%.9lf  ", store[m][j]);
        printf("\n");
    }
}