pianeti(lab1)

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


#define M 1.989*pow(10,33)
#define UA 1.496*pow(10,13)
#define anno 3.156*pow(10,7)
#define gamma (4*pow(M_PI,2)*pow(UA,3))/pow(anno,2)


//prototipi

struct coordinate {

  double x;
  double y;

};

struct rvf {

  double ra;
  double rb;
  double va;
  double vb;
  double fia;
  double fib;

};


double stampa(char *messaggio);

void verletauto(struct coordinate a0,struct coordinate b0, struct coordinate v0a, struct coordinate v0b,double ma,double mb,double dt,double e,long int numeropassi);


main(){

  double ma,mb,tempototale,dt,e;
  struct coordinate a0,b0,v0a,v0b;
  long int numeropassi;

  // chiediamo i dati iniziali

  tempototale=stampa("inserire il tempo totale di integrazione");
  dt=stampa("inserire dt");
  numeropassi=(long int)(tempototale/dt);

  ma=stampa("inserire la massa del primo pianeta");
  mb=stampa("inserire la massa del secondo pianeta");

  a0.x=stampa("inserire la coordinata x del primo pianeta");
  a0.y=stampa("inserire la coordinata y del primo pianeta");
  b0.x=stampa("inserire la coordinata x del secondo pianeta");
  b0.y=stampa("inserire la coordinata y del secondo pianeta");

  v0a.x=stampa("inserire la coordinata x della velocita' del primo pianeta");
  v0a.y=stampa("inserire la coordinata y della velocita' del primo pianeta");
  v0b.x=stampa("inserire la coordinata x della velocita' del secondo pianeta");
  v0b.y=stampa("inserire la coordinata y della velocita' del secondo pianeta");

  e=stampa("inserire la costante di accoppiamento");

  // chiamiamo la funzione che integra la traiettoria dei pianeti

  verletauto(a0,b0,v0a,v0b,ma,mb,dt,e,numeropassi);

}


double stampa(char *messaggio){

  double input;

  printf("%s\n",messaggio);
  fflush(stdout);

  scanf("%lf",&input);

  printf("il valore inserito e'%lf\n\n\n",input);

  return input;
}


void verletauto(struct coordinate a0,struct coordinate b0, struct coordinate v0a, struct coordinate v0b,double ma,double mb,double dt,double e,long int numeropassi){

   long int i;
   double mua,mub,ra0,rb0,rab0,rabn,anx,any,bnx,bny,abnx,abny,vnax,vnay,vnbx,vnby;
  struct coordinate ab0,fia0,fib0,fian,fibn;
  struct rvf fine;
  FILE *fp;

  fp=fopen("pianeti.dat","w");

  if((fp=fopen("pianeti.dat","w"))==NULL){
    printf("errore nell'apertura del file pianeti.dat\n");
    exit(EXIT_FAILURE);
  }

  fprintf(fp,"# tempo     ax     ay     bx      by       ra      rb       va     vb        fia      fib\n");

  mua=ma/M;
  mub=mb/M;


  for(i=0;i<numeropassi;i++){

    ab0.x=b0.x-a0.x;
    ab0.y=b0.y-a0.y;

    ra0=sqrt(pow(a0.x,2)+pow(a0.y,2));
    rb0=sqrt(pow(b0.x,2)+pow(b0.y,2));

    rab0=rb0-ra0;

    fia0.x=((-gamma*a0.x)/pow(ra0,3))+(e*((gamma*mub*ab0.x)/pow(rab0,3)));
    fia0.y=((-gamma*a0.y)/pow(ra0,3))+(e*((gamma*mub*ab0.y)/pow(rab0,3)));

    fib0.x=((-gamma*b0.x)/pow(rb0,3))-(e*((gamma*mua*ab0.x)/pow(rab0,3)));
    fib0.y=((-gamma*b0.y)/pow(rb0,3))-(e*((gamma*mua*ab0.y)/pow(rab0,3)));


    anx=a0.x+(v0a.x*dt)+(0.5*fia0.x*dt*dt);
    any=a0.y+(v0a.y*dt)+(0.5*fia0.y*dt*dt);

    bnx=b0.x+(v0b.x*dt)+(0.5*fib0.x*dt*dt);
    bny=b0.y+(v0b.y*dt)+(0.5*fib0.y*dt*dt);


    abnx=bnx-anx;
    abny=bny-any;

    fine.ra=sqrt(pow(anx,2)+pow(any,2));
    fine.rb=sqrt(pow(bnx,2)+pow(bny,2));

    rabn=fine.rb-fine.ra;


    fian.x=((-gamma*anx)/pow(fine.ra,3))+(e*((gamma*mub*(abnx))/pow(rabn,3)));
    fian.y=((-gamma*any)/pow(fine.ra,3))+(e*((gamma*mub*(abny))/pow(rabn,3)));

    fibn.x=((-gamma*bnx)/pow(fine.rb,3))-(e*((gamma*mua*(abnx))/pow(rabn,3)));
    fibn.y=((-gamma*bny)/pow(fine.rb,3))-(e*((gamma*mua*(abny))/pow(rabn,3)));

    vnax=v0a.x+(0.5*(fia0.x+fian.x))*dt;
    vnay=v0a.y+(0.5*(fia0.y+fian.y))*dt;

    vnbx=v0b.x+(0.5*(fib0.x+fibn.x))*dt;
    vnby=v0b.y+(0.5*(fib0.y+fibn.y))*dt;

    fine.va=sqrt(pow(vnax,2)+pow(vnay,2));
    fine.vb=sqrt(pow(vnbx,2)+pow(vnby,2));

    fine.fia=sqrt(pow(fian.x,2)+pow(fian.y,2));
    fine.fib=sqrt(pow(fibn.x,2)+pow(fibn.y,2));

    fprintf(fp,"%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf\n",(double)i*dt,a0.x,a0.y,b0.x,b0.y,fine.ra,fine.rb,fine.va,fine.vb,fine.fia,fine.fib);

    a0.x=anx;
    a0.y=any;

    b0.x=bnx;
    b0.y=bny;

    ab0.x=abnx;
    ab0.y=abny;

    v0a.x=vnax;
    v0a.y=vnay;

    v0b.x=vnbx;
    v0b.y=vnby;

  }

  fclose(fp);

}