//
//  main.cpp
//  Kolejna_proba
//
//  Created by Syn Kury on 13/10/2019.
//  Copyright © 2019 Syn Kury. All rights reserved.
//

#include <iostream>
#include<cmath>

   //************************************
bool czyZgodne(double tab1[], double tab2[], const int n)
{
    double delta = 1e-3;
    if (std::abs(tab1[0] - tab2[2]) <= abs(delta * tab1[0])
        and
        std::abs(tab1[n-1]-tab2[n-3])  <=abs(delta * tab1[n-1]) )
    {
        return true;
    }
    else {
        return false;
    }

}

//************************************

int main(int argc, const char* argv[]) {

    double ro = 1.225;
    double ni = 23.13 * 1e-6;
    double g = 10;
    double dd = 0.1;
    const int n = 41; //ilość węzłów
    double h = dd / n;
    double dt = 0.0001;

    double a = -(dt / (2 * h));
    double b = dt / (h);

    double u[n]{};  //wektor prędkości w kroku f
    double u1[n]{}; //wektor prędkości w kroku f+1

    double p[n]{};  //wektor ciśnienia w kroku f
    double p1[n]{}; //wektor ciśnienia w kroku f+1
    double deltp[n]{};

    double D[n]{};
    double S[n]{};

    long iteracja = 0;
    long iterp = 0;

    double maxPressureIteration = 10000;
    double maxVelocityIteration = 10000;

    //************************************
    // zadanie warunków brzegowych

    u[0] = 10;
    u[n - 1] = 10;

    u1[0] = 10;
    u1[n - 1] = 10;

    //************************************
    // zadanie prędkości początkowych


    for (int i = 2; i < n - 2; i += 2)
    {
        u[i] = 0;
    }


    //************************************
    // rozpoczęcie algorytmu iteracyjnego

    do
    {
        //************************************
        // obliczenie prędkości w węzłach k=2,4...,n-2

        for (int i = 2; i < n - 2; i += 2)
        {
            u1[i] = u[i] + ni * dt * (u[i - 2] - 2 * u[i] + u[i + 2]) / (4 * h * h) - dt * ((u[i] + u[i + 2]) * (u[i] + u[i + 2]) - (u[i - 2] + u[i]) * (u[i - 2] + u[i])) / (8 * h) - dt / ro * (p[i + 1] - p[i - 1]) / (2 * h) + dt * g;
        }

        //************************************
        // obliczenie wektora D oraz S w węzłach k=3,5,...n-3
        for (int i = 3; i < n - 3; i += 2)
        {
            D[i] = u[i + 1] - u[i - 1];
            S[i] = (-((u[i + 1] + u[i + 3]) * (u[i + 1] + u[i + 3])) + 2 * (u[i - 1] + u[i + 1]) * (u[i - 1] + u[i + 1]) - ((u[i - 3] + u[i - 1]) * (u[i - 3] + u[i - 1]))) / (8 * h);

        }

        D[1] = D[3];
        S[1] = S[3];
        D[n - 2] = D[n - 4];
        S[n - 2] = D[n - 4];

        //************************************
        // obliczenie wektora p w węzłach k=3,5,...,n-3
        iterp = 0;
        do
        {

            p[1] = p[3];
            p[n - 2] = p[n - 4];


            /* for(int i=3;i<n-3;i+=2)
             {
                 p1[i]=0.5*(p[i+2]+p[i-2])-h*ro/dt*D[i]-ni*ro/(4*h)*(D[i-2]-2*D[i]+D[i+2]);
             }*/
            for (int i = 3; i < n - 3; i += 2)
            {
                deltp[i] = (-a * (deltp[i - 2] + deltp[i + 2]) - (u1[i + 1] - u1[i - 1])) / b;
            }


            for (int i = 3; i < n - 3; i += 2)
            {
                p1[i] = p[i] + deltp[i];
            }

            for (int i = 3; i < n - 3; i += 2) {

                if (u1[i + 1] - u1[i - 1] < 0.00000001)
                {
                    break;
                }
                else {

                    for (int i = 3; i < n - 3; i += 2) {
                        p[i] = p1[i];
                    }
                }
            }


            ++iterp;
        } while (iterp != maxPressureIteration);//max->2
        //p[1]=p[3];
        //p[n-2]=p[n-4];


              //p[i]=p1[i];

        //************************************
        // korekcja prędkości

        for (int i = 2; i < n - 2; i += 2)
        {
            u1[i] = u[i] + ni * dt * (u[i - 2] - 2 * u[i] + u[i + 2]) / (4 * h * h) - dt * ((u[i] + u[i + 2]) * (u[i] + u[i + 2]) - (u[i - 2] + u[i]) * (u[i - 2] + u[i])) / (8 * h) - dt / ro * (p1[i + 1] - p1[i - 1]) / (2 * h) + dt * g;
        }

        //************************************
        // sprawdzenie warunku zaczymania algorytmu

        if (czyZgodne(u1, u, n))
        {
            break;      //wyjście z pętli
        }
        else
        {
            for (int i = 0; i < n; i++)
            {
                u[i] = u1[i];     //przypisanie wartości z kolejnego kroku (f+1) do tablicy (f)

            }
        }

        iteracja++;
        std::cout << "Licze dalej " << iteracja << "/" << maxVelocityIteration << std::endl;
    } while (iteracja != maxVelocityIteration);


    //************************************
    // wyświetlenie wyników

    std::cout << "Po osiągnięciu " << iteracja << " iteracji uzyskano następujące wyniki prędkości: " << std::endl;

    for (int i = 0; i < n; i += 2)
    {
        //std::cout << round(u1[i]) << std::endl;
        std::cout << u1[i] << std::endl;
    }

    return 0;
}