Guest User

Гравит. сис.

a guest
Jun 30th, 2015
255
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 6.88 KB | None | 0 0
  1. #include <iostream>
  2. #include <cmath>
  3. #include <fstream>
  4. #include <math.h>
  5.  
  6. using namespace std;
  7.  
  8.  
  9.  
  10.  
  11. int main ()
  12. {
  13.     double x1, y1, z1, x2, y2, z2, x3, y3, z3, R1, R2, R3, G=6.67191*pow(10,-11), eps, eps0;
  14.     double Vx1, Vy1, Vz1, Vx2, Vy2, Vz2, Vx3, Vy3, Vz3, ax1, ay1, az1, ax2, ay2, az2, ax3, ay3, az3, m1, m2, m3, a1, a2, a3, V1, V2, V3;
  15.     double delta_t, t_nach, t_konech, t=0, R11, R21, R31, R12, R22, R32;
  16.     int N, i, j, shag=0;
  17.     ifstream vvod ("vvod.txt");
  18.     vvod>>m1>>x1>>y1>>z1>>R1>>Vx1>>Vy1>>Vz1>>m2>>x2>>y2>>z2>>R2>>Vx2>>Vy2>>Vz2>>m3>>x3>>y3>>z3>>R3>>Vx3>>Vy3>>Vz3>>delta_t>>t_nach>>t_konech>>eps0;
  19.     ofstream koordinatii ("koorfinatii.txt");
  20.     ofstream skorosti ("skorosti.txt");
  21.     ofstream yskoreniia ("yskoreniia.txt");
  22.     do
  23.     {
  24.         N=(t_konech-t_nach)/delta_t+1;
  25.         for(i=0;i<N;i++)
  26.         {
  27.             ax1=G*m2*(x2-x1)/(fabs(pow(x2-x1,3)))+G*m3*(x3-x1)/(fabs(pow(x3-x1,3)));
  28.             ay1=G*m2*(y2-y1)/(fabs(pow(y2-y1,3)))+G*m3*(y3-y1)/(fabs(pow(y3-y1,3)));
  29.             az1=G*m2*(z2-z1)/(fabs(pow(z2-z1,3)))+G*m3*(z3-z1)/(fabs(pow(z3-z1,3)));
  30.             a1=sqrt(ax1*ax1+ay1*(ay1)+az1*(az1));
  31.  
  32.             ax2=G*m1*(x1-x2)/(fabs(pow(x1-x2,3)))+G*m3*(x3-x2)/(fabs(pow(x3-x2,3)));
  33.             ay2=G*m1*(y1-y2)/(fabs(pow(y1-y2,3)))+G*m3*(y3-y2)/(fabs(pow(y3-y2,3)));
  34.             az2=G*m1*(z1-z2)/(fabs(pow(z1-z2,3)))+G*m3*(z3-z2)/(fabs(pow(z3-z2,3)));
  35.             a2=sqrt(ax2*(ax2)+ay2*(ay2)+az2*(az2));
  36.  
  37.             ax3=G*m1*(x1-x3)/(fabs(pow(x1-x3,3)))+G*m2*(x2-x3)/(fabs(pow(x2-x3,3)));
  38.             ay3=G*m1*(y1-y3)/(fabs(pow(y1-y3,3)))+G*m2*(y2-y3)/(fabs(pow(y2-y3,3)));
  39.             az3=G*m1*(z1-z3)/(fabs(pow(z1-z3,3)))+G*m2*(z2-z3)/(fabs(pow(z2-z3,3)));
  40.             a3=sqrt(ax3*(ax3)+ay3*(ay3)+az3*(az3));
  41.  
  42.             yskoreniia<<a1<<'\t'<<a2<<'\t'<<a3<<endl;
  43.  
  44.             x1+=Vx1*delta_t+ax1*delta_t*(delta_t)/2;
  45.             y1+=Vy1*delta_t+ay1*delta_t*(delta_t)/2;
  46.             z1+=Vz1*delta_t+az1*delta_t*(delta_t)/2;
  47.  
  48.             x2+=Vx2*delta_t+ax2*delta_t*(delta_t)/2;
  49.             y2+=Vy2*delta_t+ay2*delta_t*(delta_t)/2;
  50.             z2+=Vz2*delta_t+az2*delta_t*(delta_t)/2;
  51.  
  52.             x3+=Vx3*delta_t+ax3*delta_t*(delta_t)/2;
  53.             y3+=Vy3*delta_t+ay3*delta_t*(delta_t)/2;
  54.             z3+=Vz3*delta_t+az3*delta_t*(delta_t)/2;
  55.  
  56.             koordinatii<<x1<<'/t'<<y1<<'/t'<<z1<<'/t'<<x2<<'/t'<<y2<<'/t'<<z2<<'/t'<<x3<<'/t'<<y3<<'/t'<<z3<<'/t'<<endl;
  57.  
  58.             Vx1+=ax1*delta_t;
  59.             Vy1+=ay1*delta_t;
  60.             Vz1+=az1*delta_t;
  61.             V1=sqrt(Vx1*(Vx1)+Vy1*(Vy1)+Vz1*(Vz1));
  62.  
  63.             Vx2+=ax2*delta_t;
  64.             Vy2+=ay2*delta_t;
  65.             Vz2+=az2*delta_t;
  66.             V2=sqrt(Vx2*(Vx2)+Vy2*(Vy2)+Vz2*(Vz2));
  67.  
  68.             Vx3+=ax3*delta_t;
  69.             Vy3+=ay3*delta_t;
  70.             Vz3+=az3*delta_t;
  71.             V3=sqrt(Vx3*(Vx3)+Vy3*(Vy3)+Vz3*(Vz3));
  72.  
  73.             skorosti<<V1<<'/t'<<V2<<'/t'<<V3<<endl;
  74.  
  75.             if(sqrt(pow(x1-x2,2)+pow(y1-y2,2)+pow(z1-z2,2))<=(R1+R2))
  76.             {
  77.                 cout<<"Tela 1 i 2 stolknylis'";
  78.                 i=N;
  79.             }
  80.             else
  81.             {
  82.                 if(sqrt(pow(x1-x3,2)+pow(y1-y3,2)+pow(z1-z3,2))<=(R1+R3))
  83.                 {
  84.                     cout<<"Tela 1 i 3 stolknylis'";
  85.                     i=N;
  86.                 }
  87.                 else
  88.                 {
  89.                     if(sqrt(pow(x2-x3,2)+pow(y2-y3,2)+pow(z2-z3,2))<=(R2+R3))
  90.                     {
  91.                         cout<<"Tela 2 i 3 stolknylis'";
  92.                         i=N;
  93.                     }
  94.                 }
  95.             }
  96.  
  97.         shag++;
  98.         cout<<shag<<endl;
  99.         }
  100.         R11=sqrt(pow(x1,2)+pow(y1,2)+pow(z1,2));
  101.         R21=sqrt(pow(x2,2)+pow(y2,2)+pow(z2,2));
  102.         R31=sqrt(pow(x3,2)+pow(y3,2)+pow(z3,2));
  103.  
  104.         delta_t=delta_t/2;
  105.         N=(t_konech-t_nach)/delta_t+1;
  106.         for(i=0;i<N;i++)
  107.         {
  108.             ax1=G*m2*(x2-x1)/(fabs(pow(x2-x1,3)))+G*m3*(x3-x1)/(fabs(pow(x3-x1,3)));
  109.             ay1=G*m2*(y2-y1)/(fabs(pow(y2-y1,3)))+G*m3*(y3-y1)/(fabs(pow(y3-y1,3)));
  110.             az1=G*m2*(z2-z1)/(fabs(pow(z2-z1,3)))+G*m3*(z3-z1)/(fabs(pow(z3-z1,3)));
  111.             a1=sqrt(ax1*ax1+ay1*(ay1)+az1*(az1));
  112.  
  113.             ax2=G*m1*(x1-x2)/(fabs(pow(x1-x2,3)))+G*m3*(x3-x2)/(fabs(pow(x3-x2,3)));
  114.             ay2=G*m1*(y1-y2)/(fabs(pow(y1-y2,3)))+G*m3*(y3-y2)/(fabs(pow(y3-y2,3)));
  115.             az2=G*m1*(z1-z2)/(fabs(pow(z1-z2,3)))+G*m3*(z3-z2)/(fabs(pow(z3-z2,3)));
  116.             a2=sqrt(ax2*(ax2)+ay2*(ay2)+az2*(az2));
  117.  
  118.             ax3=G*m1*(x1-x3)/(fabs(pow(x1-x3,3)))+G*m2*(x2-x3)/(fabs(pow(x2-x3,3)));
  119.             ay3=G*m1*(y1-y3)/(fabs(pow(y1-y3,3)))+G*m2*(y2-y3)/(fabs(pow(y2-y3,3)));
  120.             az3=G*m1*(z1-z3)/(fabs(pow(z1-z3,3)))+G*m2*(z2-z3)/(fabs(pow(z2-z3,3)));
  121.             a3=sqrt(ax3*(ax3)+ay3*(ay3)+az3*(az3));
  122.  
  123.             yskoreniia<<a1<<'\t'<<a2<<'\t'<<a3<<endl;
  124.  
  125.             x1+=Vx1*delta_t+ax1*delta_t*(delta_t)/2;
  126.             y1+=Vy1*delta_t+ay1*delta_t*(delta_t)/2;
  127.             z1+=Vz1*delta_t+az1*delta_t*(delta_t)/2;
  128.  
  129.             x2+=Vx2*delta_t+ax2*delta_t*(delta_t)/2;
  130.             y2+=Vy2*delta_t+ay2*delta_t*(delta_t)/2;
  131.             z2+=Vz2*delta_t+az2*delta_t*(delta_t)/2;
  132.  
  133.             x3+=Vx3*delta_t+ax3*delta_t*(delta_t)/2;
  134.             y3+=Vy3*delta_t+ay3*delta_t*(delta_t)/2;
  135.             z3+=Vz3*delta_t+az3*delta_t*(delta_t)/2;
  136.  
  137.             koordinatii<<x1<<'/t'<<y1<<'/t'<<z1<<'/t'<<x2<<'/t'<<y2<<'/t'<<z2<<'/t'<<x3<<'/t'<<y3<<'/t'<<z3<<'/t'<<endl;
  138.  
  139.             Vx1+=ax1*delta_t;
  140.             Vy1+=ay1*delta_t;
  141.             Vz1+=az1*delta_t;
  142.             V1=sqrt(Vx1*(Vx1)+Vy1*(Vy1)+Vz1*(Vz1));
  143.  
  144.             Vx2+=ax2*delta_t;
  145.             Vy2+=ay2*delta_t;
  146.             Vz2+=az2*delta_t;
  147.             V2=sqrt(Vx2*(Vx2)+Vy2*(Vy2)+Vz2*(Vz2));
  148.  
  149.             Vx3+=ax3*delta_t;
  150.             Vy3+=ay3*delta_t;
  151.             Vz3+=az3*delta_t;
  152.             V3=sqrt(Vx3*(Vx3)+Vy3*(Vy3)+Vz3*(Vz3));
  153.  
  154.             skorosti<<V1<<'/t'<<V2<<'/t'<<V3<<endl;
  155.  
  156.             if(sqrt(pow(x1-x2,2)+pow(y1-y2,2)+pow(z1-z2,2))<=(R1+R2))
  157.             {
  158.                 cout<<"Tela 1 i 2 stolknylis'";
  159.                 i=N;
  160.             }
  161.             else
  162.             {
  163.                 if(sqrt(pow(x1-x3,2)+pow(y1-y3,2)+pow(z1-z3,2))<=(R1+R3))
  164.                 {
  165.                     cout<<"Tela 1 i 3 stolknylis'";
  166.                     i=N;
  167.                 }
  168.                 else
  169.                 {
  170.                     if(sqrt(pow(x2-x3,2)+pow(y2-y3,2)+pow(z2-z3,2))<=(R2+R3))
  171.                     {
  172.                         cout<<"Tela 2 i 3 stolknylis'";
  173.                         i=N;
  174.                     }
  175.                 }
  176.             }
  177.         }
  178.         R12=sqrt(pow(x1,2)+pow(y1,2)+pow(z1,2));
  179.         R22=sqrt(pow(x2,2)+pow(y2,2)+pow(z2,2));
  180.         R32=sqrt(pow(x3,2)+pow(y3,2)+pow(z3,2));
  181.         eps=fabs((R11-R12)/R12)+fabs((R21-R22)/R22)+fabs((R31-R32)/R32);
  182.  
  183.     }while (eps>eps0);
  184.  
  185.  
  186.     return 0;
  187. }
Advertisement
Add Comment
Please, Sign In to add comment