//Hanna Szreter, Mateusz Wasilewski, Karol Żala sem.VII IŚ
//Metody Sanitarne
//Zadanie Komin
#include <fstream>
#include <iostream>
#include <cmath>
using namespace std;

float g=9.80665,Cd=0.44,Fc,ro=1.2,rom=2500.0,ms=0.5,V,vx,vx2,vy,vy2,vx0,vy0=4.0,dt=0.001,u=3.0,d=0.0001,e=0.00001,vt,z1,z2,x,y=50.0;
int a,j;

float vsx(float vx0, float vx)//ruch poziomy
{
    vt=vx0+dt/2.0*z1*((u-vx0)*(u-vx0)+(u-vx)*(u-vx));
    return vt;
}

float vsyw(float vy0, float vy)//wznoszenie
{
    vt=vy0+dt/2.0*(-z2-z1*(vy0*vy0+vy*vy));
    return vt;
}

float vsyo(float vy0, float vy)//opadanie
{
    vt=vy0+dt/2.0*(z2-z1*(vy0*vy0+vy*vy));
    return vt;
}

int main()
{
    ofstream dane("dane.txt");//otworzenie pliku
    dane<<"Dane wejsciowe:\nVwiatru=3.0\nVwylotowe=4.0\nCd=0.44\nd=0.0001"<<endl;
    dane<<"Czas X Y Vsx Vsy"<<endl;
    dane<<0<<" "<<x<<" "<<y<<" "<<0<<" "<<vy0<<endl;
    V=4.0/3.0*M_PI*d*d*d/8.0;//objętosc czastki
    Fc=d*d*M_PI/4.0;//pole przekroju czastki
    z1=Cd*Fc*ro/(2.0*(rom+ms*ro)*V);//stała pomocnicza
    z2=2.0*(rom-ro)*g/(rom+ms*ro);//stała pomocnicza
        do
        {
            j++;//zliczanie kolejnych kroków
            do
            {
                vx2=vx;//przypisanie obecnego przybliżenia Vsx do zmiennej pomocniczej (do porównania kolejnych wartosci)
                vx=vsx(vx0,vx);//wyliczenie kolejnego przybliżenia metodą trapezową
            }
            while (abs(vx-vx2)>=e);//warunek zakończenia kolejnych przybliżeń Vsx
            if (a==0)//warunek dla wznoszenia
            {
                do
                {
                    vy2=vy;//przypisanie obecnego przybliżenia Vsy do zmiennej pomocniczej (do porównania kolejnych wartosci)
                    vy=vsyw(vy0,vy);//wyliczenie kolejnego przybliżenia metodą trapezową
                }
                while (abs(vy-vy2)>=e);//warunek zakończenia kolejnych przybliżeń Vsy
            }
            if (vy<=0)//waunek rozpoczęcia opadania
            {
                a++;//zmiana wartosci a odpowiedzialnej za wznoszenie/opadanie
                vy=vy2;//anulowanie poprzedniego wyliczenia Vsy-cząstka zaczyna spadać, a liczylismy wznoszenie
            }
            if (a==1)//warunek dla opadania
            {
                do
                {
                    vy2=vy;//przypisanie obecnego przybliżenia Vsy do zmiennej pomocniczej (do porównania kolejnych wartosci)
                    vy=vsyo(vy0,vy);//wyliczenie kolejnego przybliżenia metodą trapezową
                }
                while (abs(vy-vy2)>=e);//warunek zakończenia kolejnych przybliżeń Vsy
            }
            x+=dt*vx;//położenie x w nowym punkcie
            if(a==0) y+=dt*vy;//położenie y w nowym punkcie (wznoszenie)
            else y-=dt*vy;//położenie y w nowym punkcie (opadanie)
            vx0=vx;//nowa prędkosć początkowa w kierunku x dla kolejnego punktu
            vy0=vy;//nowa prędkosć początkowa w kierunku y dla kolejnego punktu
            if(y>=0) dane<<dt*j<<" "<<x<<" "<<y<<" "<<vx<<" "<<vy<<endl;//zapis na bieżąco do pliku .txt
        }
        while (y>=0);//warunek zakończenia obliczeń (y osiąga 0 czyli cząstka spada na grunt)
    dane.close();//zamykanie pliku
    return 0;
}