#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#define g 9.81
#define epsilon 0.1
#define dt 0.01

///////////////////////////////////////////////////////////struct
struct vett{
    double x,y;
    }vento;    //creo la struttura vettore e contestualmente creo un vettore vento (si può fare aggiungendo typedef)


struct punto{
    struct vett p,v,a;  //creo struttura che segue il punto in tutti i suoi valori (//)
}punt;

//////////////////////////////////////////////////////////funzioni
struct vett aa(double vx,double vy,double ventx,double venty, double fA);
struct punto kutta(struct punto fpunt, struct vett fvento, double fA);

//---------------------
main(){
int i=0;
double A=0.1, modv, theta=M_PI/5, flag=0;
vento.x=-8;
vento.y=0;
modv=0;
flag=0;
        do{
            punt.v.x=modv*cos(theta);
            punt.v.y=modv*sin(theta);
            punt.p.x=0;
            punt.p.y=0;
                do{
                    punt=kutta(punt, vento, A);
if((fabs(punt.p.y-20)<epsilon)&&(fabs(punt.p.x-100)<epsilon)){flag=1;}
                }while(punt.p.y>0);
            modv+=epsilon;
        }while(flag=0);

printf("la velocità con cui a 45 gradi raggiunge il bersaglio in (100,20) con approssimazione di %lf equivale a %lf \n", epsilon, (modv-epsilon));







//-----------------------------test-----------------------

            punt.p.x=0;
            punt.p.y=0;
            punt.v.x=modv*cos(theta);
            punt.v.y=modv*sin(theta*M_PI/180);
    FILE* fp=fopen("bers.dat","w");
                    punt=kutta(punt, vento, A);
        fprintf(fp, " %lf %lf %lf %lf %lf %lf \n", punt.p.x,punt.p.y,punt.v.x,punt.v.y,punt.a.x,punt.a.y);
}
//---------------------


struct vett aa(double vx,double vy,double ventx,double venty, double fA){ //prende le velocità del proiettile e del vento e scrive - tramite puntatore - l'accelerazione risultante nella struct punt.a

    struct vett acc;
    acc.x=-fA*fabs(vx-ventx)*(vx-ventx);
    acc.y=-fA*fabs(vy-venty)*(vy-venty)-g;
        //printf(" vettore acc %lf %lf \n", acc.x, acc.y);
    return acc;
}

struct punto kutta(struct punto fpunt, struct vett fvento, double fA){

    double x[5],y[5],vx[5],vy[5],ax[5],ay[5],ventx,venty;
    struct vett a[5];

    x[0]=fpunt.p.x;  //deposito il contenuto della struct fpunt in più array
    y[0]=fpunt.p.y;
    vx[0]=fpunt.v.x;
    vy[0]=fpunt.v.y;
    ventx=fvento.x;
    venty=fvento.y;

//inizio algoritmo

    a[0]=aa(vx[0],vy[0],ventx,venty,fA);  //eseguo aa(...) per ottenere l'accelerazione, che verrà scritta nell'array struct a[]
            //printf("A0=%lf %lf \n", a[0].x,a[0].y);
    vx[1]=a[0].x*dt;
    x[1]=vx[0]*dt;
    vy[1]=a[0].y*dt;
    y[1]=vy[0]*dt;

    a[1]=aa(vx[0]+0.5*vx[1],vy[0]+0.5*vy[1],ventx,venty,fA);
    vx[2]=a[1].x*dt;
    x[2]=(vx[0]+0.5*vx[1])*dt;
    vy[2]=a[1].y*dt;
    y[2]=(vy[0]+0.5*vy[1])*dt;

    a[2]=aa(vx[0]+0.5*vx[2],vy[0]+0.5*vy[2],ventx,venty,fA);
    vx[3]=a[2].x*dt;
    x[3]=(vx[0]+0.5*vx[2])*dt;
    vy[3]=a[2].y*dt;
    y[3]=(vy[0]+0.5*vy[2])*dt;

    a[3]=aa(vx[0]+vx[3],vy[0]+vy[3],ventx,venty,fA);
    vx[4]=a[3].x*dt;
    x[4]=(vx[0]+vx[3])*dt;
    vy[4]=a[3].y*dt;
    y[4]=(vy[0]+vy[3])*dt;

    fpunt.p.x=x[0]+(x[1]+2*x[2]+2*x[3]+x[4])/6;
    fpunt.p.y=y[0]+(y[1]+2*y[2]+2*y[3]+y[4])/6;
    fpunt.v.x=vx[0]+(vx[1]+2*vx[2]+2*vx[3]+vx[4])/6;
    fpunt.v.y=vy[0]+(vy[1]+2*vy[2]+2*vy[3]+vy[4])/6;
    fpunt.a.x=ax[0];
    fpunt.a.y=ay[0];
//printf("px=%lf,py=%lf,vx=%lf,vy=%lf,ax=%lf,ay=%lf \n", punt.p.x,punt.p.y,punt.v.x,punt.v.y,punt.a.x,punt.a.y);
    return fpunt;

}