Cannon shot!

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

#define g 9.81
#define EARTH 0
#define TOWER 100
#define POSITIVE 9
#define NEVER_MIND 10
#define YES 1
#define NO 2

struct pos_vel{
  double x;
  double y;
  double v;
} xyAndv0={0,0,20}, xyAndv;

struct velocities{
  double x;
  double y;
} velocity;

void energy(FILE*, struct pos_vel xyAndv, struct pos_vel xyAndv0, struct velocities velocity, int);

void heavy_motion( struct pos_vel, struct pos_vel, struct velocities, FILE*, FILE*, int, int, double, int);

void relative_error( double, double, struct pos_vel, struct pos_vel, struct velocities, double, double, int);

struct pos_vel only_motion_rk4( double, double, struct pos_vel, struct pos_vel, struct velocities, int, double, FILE*, int, int );

void gnuplot_command( FILE* );

void energy(FILE *fp, struct pos_vel xyAndv, struct pos_vel xyAndv0, struct velocities velocity, int choice ){
  double dt=1000;
  int i;
  for(i=0; i<7; i++){
    xyAndv=only_motion_rk4( dt, 45., xyAndv0, xyAndv, velocity, 0, 0, NULL, NO, choice );
    fprintf(fp,"%lf          %lf\n", dt, ((1./2.*pow(xyAndv.v,2.)+g*xyAndv.y)-(1./2.*pow(xyAndv0.v,2.)+g*xyAndv0.y))/(1./2.*pow(xyAndv0.v,2.)+g*xyAndv0.y) );
    dt*=0.1;
  }
}

struct pos_vel only_motion_rk4( double dt, double theta, struct pos_vel xyAndv0, struct pos_vel xyAndv, struct velocities velocity, int h, double A, FILE*fp, int answer, int answer2){
  int i;
  double x, y, x1, x2, y1, y2, vx1, vx2, vy1, vy2, rad, yn;
  struct pos_vel xyvSUM, s1, s2, s3, s4;
  struct velocities vxvySUM, v1, v2, v3, v4;
  rad=(theta*M_PI)/(180.);
  xyAndv0.y=h;
  yn=xyAndv0.y;
  xyAndv=xyAndv0;
  velocity.x=xyAndv.v*cos(rad);
  velocity.y=xyAndv.v*sin(rad);
  for(i=0; xyAndv.y>=0 ; i++){
    if(answer==YES){
      fprintf(fp, "%lf         %lf         %lf        %lf\n", i*dt, xyAndv.x, yn, 1./2.*pow(xyAndv.v,2.)+g*yn );
    }
    x=xyAndv.x;
    y=xyAndv.y;
    s1.x=velocity.x*dt;
    s1.y=velocity.y*dt;
    v1.x=(-A*xyAndv.v*velocity.x)*dt;
    v1.y=(-g-A*xyAndv.v*velocity.y)*dt;
    xyvSUM.x=xyAndv.x+0.5*s1.x;
    xyvSUM.y=xyAndv.y+0.5*s1.y;
    vxvySUM.x=velocity.x+0.5*v1.x;
    vxvySUM.y=velocity.y+0.5*v1.y;
    xyvSUM.v=sqrt(pow(vxvySUM.x,2.)+pow(vxvySUM.y,2.));
    s2.x=vxvySUM.x*dt;
    s2.y=vxvySUM.y*dt;
    v2.x=(-A*xyvSUM.v*vxvySUM.x)*dt;
    v2.y=(-g-A*xyvSUM.v*vxvySUM.y)*dt;
    xyvSUM.x=xyAndv.x+0.5*s2.x;
    xyvSUM.y=xyAndv.y+0.5*s2.y;
    vxvySUM.x=velocity.x+0.5*v2.x;
    vxvySUM.y=velocity.y+0.5*v2.y;
    s3.x=vxvySUM.x*dt;
    s3.y=vxvySUM.y*dt;
    v3.x=(-A*xyvSUM.v*vxvySUM.x)*dt;
    v3.y=(-g-A*xyvSUM.v*vxvySUM.y)*dt;
    xyvSUM.x=xyAndv.x+s3.x;
    xyvSUM.y=xyAndv.y+s3.y;
    vxvySUM.x=velocity.x+v3.x;
    vxvySUM.y=velocity.y+v3.y;
    s4.x=vxvySUM.x*dt;
    s4.y=vxvySUM.y*dt;
    v4.x=(-A*xyvSUM.v*vxvySUM.x)*dt;
    v4.y=(-g-A*xyvSUM.v*vxvySUM.y)*dt;
    xyAndv.x+=1./6.*(s1.x+2.*s2.x+2.*s3.x+s4.x);
    xyAndv.y+=1./6.*(s1.y+2.*s2.y+2.*s3.y+s4.y);
    velocity.x+=1./6.*(v1.x+2.*v2.x+2*v3.x+v4.x);
    velocity.y+=1./6.*(v1.y+2.*v2.y+2*v3.y+v4.y);
    xyAndv.v=sqrt(pow(velocity.x,2.)+pow(velocity.y,2.));
    yn=xyAndv.y;
    if(xyAndv.y<0 && answer2==YES){
      xyAndv.x=x+((xyAndv.x-x)*y)/(y-xyAndv.y);
      yn=0;
    }
  }
  if(answer==YES){
    fprintf(fp, "%lf       %lf        %lf         %lf\n", i*dt, xyAndv.x, yn, 1./2.*pow(xyAndv.v,2.)+g*yn );
  };
  xyAndv.y=yn;
  return xyAndv;
}

void relative_error( double dt, double distance_max, struct pos_vel xyAndv0, struct pos_vel xyAndv, struct velocities velocity, double theta_max, double A, int h){
  int i;
  double max_distance, d1, d2, d3, precisione;
  for( dt=1; fabs((max_distance-d1)/((max_distance+d1)/2))>pow(10,-4); dt*=0.5 ){
    xyAndv=only_motion_rk4(2.*dt,theta_max, xyAndv0, xyAndv, velocity, h, A, NULL, NO, YES);
    d1=xyAndv.x;
    xyAndv=only_motion_rk4(dt,theta_max, xyAndv0, xyAndv, velocity, h, A, NULL, NO, YES);
    max_distance=xyAndv.x;
    while(max_distance+d1==0){
      dt*=0.5;
      xyAndv=only_motion_rk4(2.*dt,theta_max, xyAndv0, xyAndv, velocity, h, A, NULL, NO, YES);
      d1=xyAndv.x;
      xyAndv=only_motion_rk4(dt,theta_max, xyAndv0, xyAndv, velocity, h, A, NULL, NO, YES);
      max_distance=xyAndv.x;
    }
  }
  xyAndv=only_motion_rk4(2.*0.001,theta_max, xyAndv0, xyAndv, velocity, h, A, NULL, NO, YES);
  d2=xyAndv.x;
  xyAndv=only_motion_rk4(0.001,theta_max, xyAndv0, xyAndv, velocity, h, A, NULL, NO, YES);
  d3=xyAndv.x;
  precisione=fabs(d2-d3)*100.;
  printf("Fissati h = %d m e A = %lf 1/m e facendo variare il dt, l'errore relativo commesso sulla gittata e' : %lf, per un dt di %lf s.\n La precisione che si ha se ha un dt = 0.001 e' del %lf percento. \n\n", h, A, fabs((max_distance-d1)/((max_distance+d1)/2)), dt, precisione);
}

void heavy_motion( struct pos_vel xyAndv0, struct pos_vel xyAndv, struct velocities velocity, FILE *fp, FILE *mp, int h, int ang, double A, int answer2){
  int i;
  double T, x, y, xm, vxp, vyp, rad, yn, distance_max=0, theta_max, theta, dt=0.001;
  for(theta=ang; theta<91; theta++){
    xyAndv=only_motion_rk4(dt,theta, xyAndv0, xyAndv, velocity, h, A, NULL, NO, answer2);
    if(mp!=NULL){
      fprintf(mp, "%lf       %lf\n", theta, xyAndv.x);
    }
    if(distance_max<xyAndv.x){
      distance_max=xyAndv.x;
      theta_max=theta;
    }
  }
  theta=theta_max;
  for(i=0; i<800; i++){
    xyAndv=only_motion_rk4(dt,(theta-1+i*0.0025), xyAndv0, xyAndv, velocity, h, A, NULL, NO, answer2);
    if(distance_max<xyAndv.x){
      distance_max=xyAndv.x;
      theta_max=theta-1+i*0.0025;
    }
  }
  if(answer2==YES){
    printf("Dati iniziali: h = %d - A = %lf.\nDati finali: Rmax = %lf - Theta = %lf°. \n", h, A, distance_max, theta_max);
    relative_error(dt,distance_max,xyAndv0,xyAndv,velocity,theta_max,A,h);
  }
  xyAndv=only_motion_rk4(dt,theta_max, xyAndv0, xyAndv, velocity, h, A, fp, YES, answer2);
}

void gnuplot_command( FILE *fp ){
  fprintf(fp,"set term pos col \n");
  fprintf(fp,"set out 'grafici.ps' \n");
  fprintf(fp,"plot 'moto_grave0.dat' u 2:3 \n");
  fprintf(fp,"plot 'moto_grave0.dat' u 1:4 \n");
  fprintf(fp,"plot 'moto_grave100.dat' u 2:3 \n");
  fprintf(fp,"plot 'moto_grave100.dat' u 1:4 \n");
  fprintf(fp,"plot 'moto_grave_A01.dat' u 2:3 \n");
  fprintf(fp,"plot 'moto_grave_A01.dat' u 1:4 \n");
  fprintf(fp,"plot 'moto_grave_A1.dat' u 2:3 \n");
  fprintf(fp,"plot 'moto_grave_A1.dat' u 1:4 \n");
  fprintf(fp,"plot 'errore_energia.dat' \n");
  fprintf(fp,"plot 'errore_energia1.dat' \n");
  fprintf(fp,"set terminal x11");
 }

main(){
  int i;
  FILE *hm, *hmh, *fm, *fmA, *inter, *command, *ene, *ene1, *dRdT;
  if((hm=fopen("moto_grave0.dat", "w"))==NULL){
    printf("Errore nell'apertura file moto_grave0.dat.\n");
    exit(EXIT_FAILURE);
  }
  if((hmh=fopen("moto_grave100.dat", "w"))==NULL){
    printf("Errore nell'apertura file moto_grave100.dat.\n");
    exit(EXIT_FAILURE);
  }
  if((fm=fopen("moto_grave_A01.dat", "w"))==NULL){
    printf("Errore nell'apertura file moto_grave_A01.dat.\n");
    exit(EXIT_FAILURE);
  }
  if((fmA=fopen("moto_grave_A1.dat", "w"))==NULL){
    printf("Errore nell'apertura file moto_grave_A1.dat.\n");
    exit(EXIT_FAILURE);
  }
  if((inter=fopen("interpolazione.dat", "w"))==NULL){
    printf("Errore nell'apertura file interpolazione.dat.\n");
    exit(EXIT_FAILURE);
  }
  if((command=fopen("myCommand.dat", "w"))==NULL){
    printf("Errore nell'apertura file myCommnd.dat.\n");
    exit(EXIT_FAILURE);
  }
 if((ene=fopen("errore_energia.dat", "w"))==NULL){
    printf("Errore nell'apertura file errore_energia.dat.\n");
    exit(EXIT_FAILURE);
  }
 if((ene1=fopen("errore_energia1.dat", "w"))==NULL){
    printf("Errore nell'apertura file errore_energia1.dat.\n");
    exit(EXIT_FAILURE);
  }
 if((dRdT=fopen("variazione_R_in_theta.dat", "w"))==NULL){
    printf("Errore nell'apertura file variazione_R_in_theta.dat.\n");
    exit(EXIT_FAILURE);
  }
  fprintf(hm, "#Tempo         Ascissa         Ordinata         Energia\n");
  fprintf(hmh, "#Tempo         Ascissa         Ordinata          Energia\n");
  fprintf(fm, "#Tempo         Ascissa         Ordinata         Energia\n");
  fprintf(fmA, "#Tempo         Ascissa         Ordinata          Energia\n");
  fprintf(inter, "#Tempo         Ascissa         Ordinata          Energia\n");
  fprintf(ene, "#dt         DE/E0 \n");
  fprintf(ene1, "#dt         DE/E0 \n");
  fprintf(dRdT, "#theta         R(theta)\n");
  printf("Il programma simula il moto bidimensionale di un grave, sparato da un cannone con velocita' iniziale v0=20 m/s.\n");
  printf("\n");
  heavy_motion(xyAndv0, xyAndv, velocity, hm, NULL, EARTH, 0, 0.,YES);
  heavy_motion(xyAndv0, xyAndv, velocity, inter, NULL, EARTH, 0, 0.,NO);
  heavy_motion(xyAndv0, xyAndv, velocity, hmh, NULL, TOWER, -90, 0.,YES);
  heavy_motion(xyAndv0, xyAndv, velocity, fm, NULL, EARTH, 0, 0.1,YES);
  heavy_motion(xyAndv0, xyAndv, velocity, fmA, dRdT, EARTH, 0, 1.,YES);
  energy(ene, xyAndv, xyAndv0, velocity, NO);
  energy(ene1, xyAndv, xyAndv0, velocity, YES);
  gnuplot_command(command);
  fclose(hm);
  fclose(hmh);
  fclose(fm);
  fclose(fmA);
  fclose(inter);
  fclose(command);
  fclose(ene);
  fclose(ene1);
  i=system("gnuplot <myCommand.dat");
  if(i!=0){
    printf("Errore di utilizzo di gnuplot.\n");
  }
}