massima gittata kutta

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

///////////////////////////////////////////////////////////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, temp;
//////////////////////////////////////////////////////////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=20, theta=0, thetaMAX, grad=(M_PI/2)/dth;

punt.p.x=0;
punt.p.y=0;
vento.x=0;
vento.y=0;

    do{
    punt.v.x=modv*cos(theta);
    punt.v.y=modv*sin(theta);
        do{
        temp=kutta(punt, vento, A);
        printf(" %lf %lf \n", theta*180/M_PI,punt.p.x);
        }while(punt.p.y>0);
            if(temp.p.x>punt.p.x){thetaMAX=theta;}
        punt=temp;
        theta+=dth;
    }while(theta<M_PI/2);
printf("massima gittata raggiunta a %lf gradi \n", thetaMAX*180/M_PI);
}
//---------------------
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[5],venty[5];
    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[0]=fvento.x;
    venty[0]=fvento.y;

//inizio algoritmo

    a[0]=aa(vx[0],vy[0],ventx[0],venty[0],fA);  //eseguo aa(...) per ottenere l'accelerazione, che verrà scritta nell'array struct a[]

    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[0]+0.5*ventx[1],venty[0]+0.5*venty[1],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[0]+0.5*ventx[2],venty[0]+0.5*venty[2],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[0]+ventx[3],venty[0]+venty[3],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]+(ax[1]+2*ax[2]+2*ax[3]+ax[4])/6;
    fpunt.a.y=ay[0]+(ay[1]+2*ay[2]+2*ay[3]+ay[4])/6;
    return fpunt;
}