proiettile/vento (cartesiano))

#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, k, m;
printf("velocità iniziale su x e' %lf su y %lf \n",vento.x,vento.y);

printf("inserimento dati iniziali del proiettile \n");////////////////proiettile
printf("inserire posizione iniziale punto su x \n");
scanf("%lf", &punt.p.x);
printf("inserire posizione iniziale punto su y \n");
scanf("%lf", &punt.p.y);
printf("inserire velocità iniziale punto su x \n");
scanf("%lf", &punt.v.x);
printf("inserire velocità iniziale punto su y \n");
scanf("%lf", &punt.v.y);
printf("inserimento dati iniziali del vento \n");/////////////////////vento
printf("velocità iniziale vento su x \n");
scanf("%lf", &vento.x);
printf("velocità iniziale vento su y \n");
scanf("%lf", &vento.y);
printf("controllo dati inseriti");/////////////////////////////////controllo
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);
double hstart=punt.p.y;

    FILE* fp=fopen("proiettile.dat","w");
        fprintf(fp, "px=        py=        vx=       vy=     ax=         ay= \n");
    do{
        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);
        punt=kutta(punt, vento, A);
    }while(punt.p.y>hstart);


}
//---------------------


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
        //printf("scrittura aa aa aa %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 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[]
            //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[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];
    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;

}