//
//  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)
{
   for(int i=2;i<n-2;i+=2)
    {
        if(tab1[i]-tab2[i]<0.000000000000001 and tab1[4]-tab2[4]<0.000000000000001) return true;
        else return false;
    }
    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=1;
    long iterp=1;

    //************************************
    // 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

           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!=10000);
    //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++;

    } while (iteracja!=1000000);


    //************************************
    // 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;
    }



    return 0;
}