API tools faq
paste
pietrod21

pianeti(lab)

pietrod21
Dec 30th, 2011
89
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C 9.35 KB | None | 0 0
raw download clone embed print report
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <math.h>
  4. #include <float.h>
  5.  
  6. #define G 6.67e-11 //m^3/(kg*s^2)
  7. #define UA 15e10 //m unita astronomica
  8. #define MS 2e30 //kg massa sole
  9. #define ANNO 31556736 //secondi anno solare della terra
  10.  
  11.  
  12. #define NOMEFILE "moto.dat"
  13.  
  14. /****************************************/
  15. //struttura usata per i corpi
  16. struct corpo{
  17.   double x, y, m, vx, vy, x1, y1;
  18. };
  19.  
  20.   typedef struct corpo pianeta;
  21.  
  22. struct status{
  23.   pianeta p1,p2;
  24.   double gamma, tn, dt, t1, epsilon;
  25.   char integrazione;
  26. };
  27.  
  28.   typedef struct status stato;
  29. /****************************************/
  30.  
  31.  
  32. /****************************************/
  33. //funzioni generali
  34. double modulo(double x, double y);
  35. double acc(stato* s, int n, char asse);
  36. double dinput(double min, double max);
  37. void cordinput(pianeta* p);
  38. char cinput(int n, char* caratteri);
  39. /****************************************/
  40.  
  41. /****************************************/
  42. //funzioni principali
  43. void input(stato* s, char norm);
  44. void integra(stato* s);
  45. void inizializza(stato* s);
  46. void passo_integra(stato* s);
  47. /****************************************/
  48.  
  49. /****************************************/
  50. //metodi d'integrazione
  51. void verlet(stato* s);
  52. void rungekutta(stato* s);
  53. /****************************************/
  54.  
  55.  
  56. /****************************************/
  57. /****************************************/
  58. /****************************************/
  59.  
  60. //                MAIN
  61.  
  62. int main(void){
  63.   stato s;
  64.  
  65.   printf("calcolare in coordinate normalizzate al sistema terra sole\n(y/n) = ");
  66.   char norm = cinput(4, (char[]){'y', 'Y', 'n', 'N'});
  67.  
  68.   input(&s, norm);
  69.  
  70.   integra(&s);
  71.  
  72.   return 0;
  73. }
  74.  
  75. /****************************************/
  76. /****************************************/
  77. /****************************************/
  78.  
  79.  
  80. /****************************************/
  81. //funzioni generali
  82. double modulo(double x, double y){
  83.   return sqrt(pow(x, 2) + pow(y, 2));
  84. }
  85.  
  86. double acc(stato* s, int n, char asse){
  87.   if(n==1){
  88.     if(asse=='x'){
  89.       return 0 - s->gamma * (s->p1).x / pow(modulo((s->p1).x, (s->p1).y), 3) + s->epsilon * (s->p2).m * s->gamma * ((s->p2).x - (s->p1).x) / pow(modulo((s->p2).x - (s->p1).x, (s->p2).y - (s->p1).y), 3);
  90.     }else if(asse=='y'){
  91.       return 0 - s->gamma * (s->p1).y / pow(modulo((s->p1).x, (s->p1).y), 3) + s->epsilon * (s->p2).m * s->gamma * ((s->p2).y - (s->p1).y) / pow(modulo((s->p2).x - (s->p1).x, (s->p2).y - (s->p1).y), 3);
  92.     }else{
  93.       printf("\n\n\nERROR ACC asse ne x ne y\n\n\n");
  94.       exit(EXIT_FAILURE);
  95.     }
  96.   }else if(n==2){
  97.     if(asse=='x'){
  98.       return 0 - s->gamma * (s->p2).x / pow(modulo((s->p2).x, (s->p2).y), 3) + s->epsilon * (s->p1).m * s->gamma * ((s->p1).x - (s->p2).x) / pow(modulo((s->p2).x - (s->p1).x, (s->p2).y - (s->p1).y), 3);
  99.     }else if(asse=='y'){
  100.       return 0 - s->gamma * (s->p2).y / pow(modulo((s->p2).x, (s->p2).y), 3) + s->epsilon * (s->p1).m * s->gamma * ((s->p1).y - (s->p2).y) / pow(modulo((s->p2).x - (s->p1).x, (s->p2).y - (s->p1).y), 3);
  101.     }else{
  102.       printf("\n\n\nERROR ACC asse ne x ne y\n\n\n");
  103.       exit(EXIT_FAILURE);
  104.     }
  105.   }else{
  106.     printf("\n\n\nERROR ACC n ne 1 ne 2\n\n\n");
  107.     exit(EXIT_FAILURE);
  108.   }
  109. }
  110.  
  111. double dinputs(double min, double max){
  112.   double num;
  113.   do{
  114.     scanf("%lf", &num);
  115.     if(num <= min || num >= max){
  116.       printf("il valore deve essere compreso tra %lf e %lf\nreinserisci il valore = ", min, max);
  117.     }
  118.   }while(num <= min || num >= max);
  119.   return num;
  120. }
  121.  
  122. double dinput(double min, double max){
  123.   double num;
  124.   do{
  125.     scanf("%lf", &num);
  126.     if(num < min || num > max){
  127.       printf("il valore deve essere compreso tra %lf e %lf\nreinserisci il valore = ", min, max);
  128.     }
  129.   }while(num < min || num > max);
  130.   return num;
  131. }
  132.  
  133. void cordinput(pianeta* p){
  134.   do{
  135.     printf("x = ");
  136.     p->x = dinput(-DBL_MAX, DBL_MAX);
  137.     printf("y = ");
  138.     p->y = dinput(-DBL_MAX, DBL_MAX);
  139.     if(p->x==0 && p->y==0){
  140.       printf("il pianeta non puo stare nel sole, dammi altre coordinate\n");
  141.     }
  142.   }while(p->x==0 && p->y==0);
  143. }
  144.  
  145. char cinput(int n, char* caratteri){
  146.   char input;
  147.   for(;;){
  148.     scanf("%c", &input);
  149.     int i;
  150.     for(i=0;i<n;i++){
  151.       if(input==caratteri[i]){
  152.     return input;
  153.       }
  154.     }
  155.     printf("gli unici caratteri accettati sono");
  156.     for(i=0;i<n;i++){
  157.       printf(" \'%c\'", caratteri[i]);
  158.     }
  159.     printf("\ninput = ");
  160.   }
  161. }
  162. /****************************************/
  163.  
  164. /****************************************/
  165. //funzioni principali
  166. void input(stato* s, char norm){
  167.   printf("scegli il metodo di integrazione\nv = verlet\nr = runge kutta\nscegli = ");
  168.   s->integrazione = cinput(2, (char[]){'v', 'r'});
  169.  
  170.   printf("dammi la massa del sole (kg)\nM = ");
  171.   double M = dinputs(0.0, DBL_MAX);
  172.   s->gamma = M * G;
  173.  
  174.  
  175.   printf("dammi la massa del primo pianeta (kg)\nm1 = ");
  176.   (s->p1).m = dinputs(0.0, DBL_MAX) / M;
  177.  
  178.   printf("dammi la massa del secondo pianeta (kg)\nm2 = ");
  179.   (s->p2).m = dinputs(0.0, DBL_MAX) / M;
  180.  
  181.   do{
  182.     printf("dammi la posizione del primo pianeta in coordinate cartesiane (m)\n");
  183.     cordinput(&((s->p1)));
  184.    
  185.     printf("dammi la posizione del secondo pianeta in coordinate cartesiane (m)\n");
  186.     cordinput(&((s->p2)));
  187.  
  188.     if((s->p1).x==(s->p2).x && (s->p1).y==(s->p2).y){
  189.       printf("i pianeti non possono stare alle stesse coordinate\n");
  190.     }
  191.   }while((s->p1).x==(s->p2).x && (s->p1).y==(s->p2).y);
  192.  
  193.  
  194.   printf("dammi le velocita iniziali del primo pianeta (m/s)\nvx = ");
  195.   (s->p1).vx = dinput(-DBL_MAX, DBL_MAX);
  196.   printf("vy = ");
  197.   (s->p1).vy = dinput(-DBL_MAX, DBL_MAX);
  198.   printf("dammi le velocita iniziali del secondo pianeta (m/s)\nvx = ");
  199.   (s->p2).vx = dinput(-DBL_MAX, DBL_MAX);
  200.   printf("vy = ");
  201.   (s->p2).vy = dinput(-DBL_MAX, DBL_MAX);
  202.  
  203.  
  204.  
  205.   printf("dammi epsilon\neps = ");
  206.   s->epsilon = dinput(0.0, 1.0);
  207.  
  208.   printf("dammi la durata (s)\nt1 = ");
  209.   s->t1 = dinputs(0.0, DBL_MAX);
  210.  
  211.   printf("dt (s) = ");
  212.   s->dt = dinputs(0.0, DBL_MAX);
  213.  
  214.   if(norm=='y' || norm=='Y'){
  215.     s->gamma *= (pow(ANNO, 2) / pow(UA, 3));
  216.     (s->p1).x /= UA;
  217.     (s->p1).y /= UA;
  218.     (s->p2).x /= UA;
  219.     (s->p2).y /= UA;
  220.     (s->p1).vx *= ANNO / UA;
  221.     (s->p1).vy *= ANNO / UA;
  222.     (s->p2).vx *= ANNO / UA;
  223.     (s->p2).vy *= ANNO / UA;
  224.     s->t1 /= ANNO;
  225.     s->dt /= ANNO;
  226.   }
  227.  
  228.  
  229.   FILE *file;
  230.   file = fopen(NOMEFILE, "w");
  231.   fprintf(file, "%c M = %lf, m1%cM = %lf, m2%cM = %lf\n", '#', M, '/', (s->p1).m, '/', (s->p2).m);
  232.   if(fclose(file)==EOF){
  233.     printf("\n\nERRORE CHIUSURA FILE\n\n");
  234.     exit(EXIT_FAILURE);
  235.   }
  236.  
  237. }
  238.  
  239. void integra(stato* s){
  240.   FILE *file;
  241.   file = fopen(NOMEFILE, "a");
  242.  
  243.   if(file==NULL){
  244.     printf("\n\nERRORE APERTURA FILE\n\n");
  245.     exit(EXIT_FAILURE);
  246.   }
  247.  
  248.   fprintf(file, "%ctn         x1          y1          x2           y2\n", '#');
  249.  
  250.   inizializza(s);
  251.  
  252.   for(s->tn = 0; s->tn <= s->t1; s->tn += s->dt){
  253.     passo_integra(s);
  254.     fprintf(file, "%lf    %lf    %lf    %lf     %lf\n", s->tn, (s->p1).x, (s->p1).y, (s->p2).x, (s->p2).y);
  255.   }
  256.  
  257.   if(fclose(file)==EOF){
  258.     printf("\n\nERRORE CHIUSURA FILE\n\n");
  259.     exit(EXIT_FAILURE);
  260.   }
  261. }
  262.  
  263. void inizializza(stato* s){
  264.   if(s->integrazione=='v'){
  265.     (s->p1).x1 = (s->p1).x;
  266.     (s->p1).y1 = (s->p1).y;
  267.     (s->p2).x1 = (s->p2).x;
  268.     (s->p2).y1 = (s->p2).y;
  269.    
  270.     (s->p1).x = (s->p1).x1 + (s->p1).vx * s->dt;
  271.     (s->p1).y = (s->p1).y1 + (s->p1).vy * s->dt;
  272.     (s->p2).x = (s->p2).x1 + (s->p2).vx * s->dt;
  273.     (s->p2).y = (s->p2).y1 + (s->p2).vy * s->dt;
  274.   }else if(s->integrazione!='r'){
  275.     printf("\n\nERRORE S INTEGRAZIONE\n\n");
  276.     exit(EXIT_FAILURE);
  277.   }
  278. }
  279.  
  280. void passo_integra(stato* s){
  281.   if(s->integrazione=='v'){
  282.     verlet(s);
  283.   }else if(s->integrazione=='r'){
  284.     rungekutta(s);
  285.   }else{
  286.     printf("\n\nERRORE S INTEGRAZIONE\n\n");
  287.     exit(EXIT_FAILURE);
  288.   }
  289. }
  290. /****************************************/
  291.  
  292. /****************************************/
  293. //metodi d'integrazione
  294. void verlet(stato* s){
  295.   double temp1, temp2, temp3, temp4;
  296.  
  297.   temp1 = 2 * (s->p1).x - (s->p1).x1 + acc(s, 1, 'x') * pow(s->dt, 2);
  298.   temp2 = 2 * (s->p1).y - (s->p1).y1 + acc(s, 1, 'y') * pow(s->dt, 2);
  299.   temp3 = 2 * (s->p2).x - (s->p2).x1 + acc(s, 2, 'x') * pow(s->dt, 2);
  300.   temp4 = 2 * (s->p2).y - (s->p2).y1 + acc(s, 2, 'y') * pow(s->dt, 2);
  301.  
  302.  
  303.   (s->p1).x1 = (s->p1).x;
  304.   (s->p1).x = temp1;
  305.  
  306.  
  307.   (s->p1).y1 = (s->p1).y;
  308.   (s->p1).y = temp2;
  309.  
  310.  
  311.   (s->p2).x1 = (s->p2).x;
  312.   (s->p2).x = temp3;
  313.  
  314.  
  315.   (s->p2).y1 = (s->p2).y;
  316.   (s->p2).y = temp4;
  317. }
  318.  
  319. void rungekutta(stato* s){
  320.   double v1x = acc(s, 1, 'x') * s->dt / 2 + (s->p1).vx;
  321.   double v1y = acc(s, 1, 'y') * s->dt / 2 + (s->p1).vy;
  322.   double v2x = acc(s, 2, 'x') * s->dt / 2 + (s->p2).vx;
  323.   double v2y = acc(s, 2, 'y') * s->dt / 2 + (s->p2).vy;
  324.  
  325.   (s->p1).x1 = (s->p1).x;
  326.   (s->p1).y1 = (s->p1).y;
  327.   (s->p2).x1 = (s->p2).x;
  328.   (s->p2).y1 = (s->p2).y;
  329.  
  330.   (s->p1).x += (s->p1).vx * s->dt / 2;
  331.   (s->p1).y += (s->p1).vy * s->dt / 2;
  332.   (s->p2).x += (s->p2).vx * s->dt / 2;
  333.   (s->p2).y += (s->p2).vy * s->dt / 2;
  334.  
  335.   (s->p1).vx += acc(s, 1, 'x') * s->dt;
  336.   (s->p1).vy += acc(s, 1, 'y') * s->dt;
  337.   (s->p2).vx += acc(s, 2, 'x') * s->dt;
  338.   (s->p2).vy += acc(s, 2, 'y') * s->dt;
  339.  
  340.   (s->p1).x = (s->p1).x1 + v1x * s->dt;
  341.   (s->p1).y = (s->p1).y1 + v1y * s->dt;
  342.   (s->p2).x = (s->p2).x1 + v2x * s->dt;
  343.   (s->p2).y = (s->p2).y1 + v2y * s->dt;
  344. }  
  345. /****************************************/
Advertisement
Add Comment
Please, Sign In to add comment
Public Pastes
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!