pianeti1

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

#define dt 0.1
#define M 1989000
#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
double a(double vx,double vy,double ventx,double venty, double fA);
struct posizioni verlet(struct posizioni pos);
//////////////////////////////////////////////////////////main
main(){
int i;
FILE *fp, *fd;
fd=fopen("dati_png.gnu","w");
fp=fopen("pianeti.dat","w+");
fprintf(fd, "plot 'pianeti.dat' u 2:3 \n set term png \n set output 'eull.png' \n replot \n  \n ");
fprintf(fd, "plot 'pianeti.dat' u 4:5 \n set term png \n set output 'eull.png' \n replot \n  \n ");
fclose(fd);
pos.A.m=5; //terra
pos.B.m=0.6; //marte
pos.A.p.x=1.4;
pos.A.p.y=0.2;
pos.B.p.x=1.9;
pos.B.p.y=1.8;
pos.A.v.x=0;//143040
pos.A.v.y=0;
pos.B.v.x=0;//-100000
pos.B.v.y=0;


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);
}
fclose(fp);
}
//////////////////////////////////////////////////////////funzioni
double 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(sqrt((xa-xb)*(xa-xb)+(ya-yb)*(ya-yb)),3));  return acc;
}

struct posizioni verlet(struct posizioni pos){
    //ottengo x(n+1) per ottenere l'accelerazione di v(n+1): per farlo creo una nuova struct "temporanee" pos1
    struct posizioni pos1;
    pos1.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.B.p.y,pos.B.m)*dt*dt;
    pos1.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.B.p.x,pos.B.m)*dt*dt;
    pos1.B.p.x=pos.B.p.x+pos.B.v.x*dt+0.5*a(pos.B.p.x,pos.B.p.y,pos.A.p.x,pos.A.p.y,pos.A.m)*dt*dt;
    pos1.B.p.y=pos.B.p.y+pos.B.v.y*dt+0.5*a(pos.B.p.y,pos.B.p.x,pos.A.p.y,pos.A.p.x,pos.A.m)*dt*dt;
    //ottengo le nuove velocita'-  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.B.p.y,pos.B.m)+a(pos1.A.p.x,pos1.A.p.y,pos1.B.p.x,pos1.B.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.B.p.x,pos.B.m)+a(pos1.A.p.y,pos1.A.p.x,pos1.B.p.y,pos1.B.p.x,pos.B.m))*dt;
    pos.B.v.x=pos.B.v.x+0.5*(a(pos.B.p.x,pos.B.p.y,pos.A.p.x,pos.A.p.y,pos.A.m)+a(pos1.B.p.x,pos1.B.p.y,pos1.A.p.x,pos1.A.p.y,pos.A.m))*dt;
    pos.B.v.y=pos.B.v.y+0.5*(a(pos.B.p.y,pos.B.p.x,pos.A.p.y,pos.A.p.x,pos.A.m)+a(pos1.B.p.y,pos1.B.p.x,pos1.A.p.y,pos1.A.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.B.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.B.p.x,pos.B.m)*dt*dt;
    pos.B.p.x=pos.B.p.x+pos.B.v.x*dt+0.5*a(pos.B.p.x,pos.B.p.y,pos.A.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.B.p.y,pos.B.p.x,pos.A.p.y,pos.A.p.x,pos.A.m)*dt*dt;
    return pos;
}