pianeti

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#define dt 0.1
#define M 332000
#define teu 39,4784176
///////////////////////////////////////////////////////////struct
struct vett{
    double x,y;
    };

struct punto{
    struct vett p,v;
    double m;
};

struct posizioni{
    struct punto A,B;
}pos;
//////////////////////////////////////////////////////////funzioni
struct vett a(double vx,double vy,double ventx,double venty, double fA);
struct posizioni verlet(struct posizioni);
//////////////////////////////////////////////////////////main
main(){
int i;
FILE *fp;
fp=fopen("pianeti.c","w+");
pos.A.m=1;
pos.B.m=0.5;
pos.A.p.x=60;
pos.A.p.y=60;
pos.B.p.x=30;
pos.B.p.y=40;
pos.A.v.x=-30;
pos.A.v.y=-30;
pos.B.v.x=-50;
pos.B.v.y=-60;


for(i=0;i<1000;i++){
    fprintf(fp,"%lf %lf %lf %lf %lf\n",i*dt,pos.A.p.x,pos.A.p.y,pos.B.p.x,pos.B.p.y);
    pos=verlet(pos);
}

}
//////////////////////////////////////////////////////////funzioni
struct vett a(double xa,double ya,double xb,double yb, double m){
    double acc;
    acc=teu*xa/pow(sqrt(xa*xa+ya*ya),3)+teu*(-m/M)*(xa-xb)/pow((xa-xb)*(xa-xb)+(ya-yb)*(ya-yb)),3);
    return acc;
}

struct posizioni verlet(struct posizioni){
    //ottengo x(n+1) per ottenere l'accelerazione di v(n+1): per farlo creo delle nuove struct "temporanee" A1 e B1
    struct posizioni A1, B1;
    pos.A1.p.x=pos.A.p.x+pos.A.v.x*dt+0.5*a(pos.A.p.x,pos.A.p.y,pos.B.p.x,pos.A.p.y,pos.B.m)*dt*dt;
    pos.A1.p.y=pos.A.p.y+pos.A.v.y*dt+0.5*a(pos.A.p.y,pos.A.p.x,pos.B.p.y,pos.A.p.x,pos.B.m)*dt*dt;
    pos.B1.p.x=pos.B.p.x+pos.B.v.x*dt+0.5*a(pos.A.p.x,pos.A.p.y,pos.B.p.x,pos.A.p.y,pos.A.m)*dt*dt;
    pos.B1.p.y=pos.B.p.y+pos.B.v.y*dt+0.5*a(pos.A.p.y,pos.A.p.x,pos.B.p.y,pos.A.p.x,pos.A.m)*dt*dt;
    //ottengo le nuove velocità-  sfruttando le A1 e B1 precedentemente ricavate - sovrascrivendole direttamente alle vecchie
    pos.A.v.x=pos.A.v.x+0.5+(a(pos.A.p.x,pos.A.p.y,pos.B.p.x,pos.A.p.y,pos.A.m)+a(pos.A1.p.x,pos.A1.p.y,pos.B1.p.x,pos.A1.p.y,pos.B.m))*dt;
    pos.A.v.y=pos.A.v.y+0.5+(a(pos.A.p.y,pos.A.p.x,pos.B.p.y,pos.A.p.x,pos.A.m)+a(pos.A1.p.y,pos.A1.p.x,pos.B1.p.y,pos.A1.p.x,pos.B.m))*dt;
    pos.B.v.x=pos.B.v.x+0.5+(a(pos.A.p.x,pos.A.p.y,pos.B.p.x,pos.A.p.y,pos.A.m)+a(pos.A1.p.x,pos.A1.p.y,pos.B1.p.x,pos.A1.p.y,pos.A.m))*dt;
    pos.B.v.y=pos.B.v.y+0.5+(a(pos.A.p.y,pos.A.p.x,pos.B.p.y,pos.A.p.x,pos.A.m)+a(pos.A1.p.y,pos.A1.p.x,pos.B1.p.y,pos.A1.p.x,pos.A.m))*dt;
    //ottengo le nuove posizioni sovrascrivendole direttamente alle vecchie
    pos.A.p.x=pos.A.p.x+pos.A.v.x*dt+0.5*a(pos.A.p.x,pos.A.p.y,pos.B.p.x,pos.A.p.y,pos.B.m)*dt*dt;
    pos.A.p.y=pos.A.p.y+pos.A.v.y*dt+0.5*a(pos.A.p.y,pos.A.p.x,pos.B.p.y,pos.A.p.x,pos.B.m)*dt*dt;
    pos.B.p.x=pos.B.p.x+pos.B.v.x*dt+0.5*a(pos.A.p.x,pos.A.p.y,pos.B.p.x,pos.A.p.y,pos.A.m)*dt*dt;
    pos.B.p.y=pos.B.p.y+pos.B.v.y*dt+0.5*a(pos.A.p.y,pos.A.p.x,pos.B.p.y,pos.A.p.x,pos.A.m)*dt*dt;
    return pos;
}