1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
| #include <iostream>
#include <fstream>
#include <iomanip>
#include <math.h>
#include <ctime>
#include<time.h>
#include<stdlib.h>
using namespace std;
double const G = 6;
class TWOSTARS
{
public:
double F1x, F2x, F1y, F2y, F1z, F2z, ex1, ex2, ey1, ey2, ez1, ez2, exs, eys, ezs, FSy, FSz, FSx;
double rast1, rast2, rasts;
double star1_x,star1_y,star1_z; //координаты первой звезды
double vstar1_x,vstar1_y,vstar1_z;//координаты скорости первой звезды
double star2_x,star2_y,star2_z;//координаты второй звезды
double planet_x,planet_y,planet_z;//координаты планеты
double vstar2_x,vstar2_y,vstar2_z;//координаты скорости второй звезды
double vplanet_x,vplanet_y,vplanet_z;//координаты скороти планеты звезды
double m1,m2,M;//массы звезд и планеты
TWOSTARS ()//начальные координаты
{
star1_x=15 ;
star1_y=15 ;
star1_z= 15;
star2_x=20;
star2_y= 20 ;
star2_z= 20 ;
vstar1_x=1 ;
vstar1_y=1 ;
vstar1_z=-2 ;
vstar2_x= -2;
vstar2_y= -2;
vstar2_z= 4;
planet_x= 500;
planet_y= 500;
planet_z= 500;
vplanet_x= 1;
vplanet_y=1 ;
vplanet_z= 1;
m1= 20;
m2= 20;
M= 5;
}
void changeforce(double planet_x,double planet_y,double planet_z,double star1_x,double star1_y, double star1_z,double star2_x, double star2_y,double star2_z) //сила гравитационного взаимодействия между объектами
{
rast1 = sqrt( (planet_x-star1_x)*(planet_x-star1_x)+(planet_y-star1_y)*(planet_y-star1_y)+(planet_z-star1_z)*(planet_z-star1_z));// расстояние между планетой и первой звездой
rast2 = sqrt( (planet_x-star2_x)*(planet_x-star2_x)+(planet_y-star2_y)*(planet_y-star2_y)+(planet_z-star2_z)*(planet_z-star2_z));// расстояние между планетой и второй звездой
rasts = sqrt( (star1_x-star2_x)*(star1_x-star2_x)+(star1_y-star2_y)*(star1_y-star2_y)+(star1_z-star2_z)*(star1_z-star2_z));// расстояние между звездами
ex1=(planet_x-star1_x)/rast1;//проекция вектора соединяющего положение планеты и звезды 1 на ось x, нужна чтоб у проекций сил гравитации было направление
ey1=(planet_y-star1_y)/rast1;//проекция вектора соединяющего положение планеты и звезды 1 на ось y
ez1=(planet_z-star1_z)/rast1;//проекция вектора соединяющего положение планеты и звезды 1 на ось z
ex2=(planet_x-star2_x)/rast2;//проекция вектора соединяющего положение планеты и звезды 2 на ось x
ey2=(planet_y-star2_y)/rast2;//проекция вектора соединяющего положение планеты и звезды 2 на ось y
ez2=(planet_z-star2_z)/rast2;//проекция вектора соединяющего положение планеты и звезды 2 на ось z
exs=(star1_x-star2_x)/rasts;//проекция вектора соединяющего положение звезд на ось x
eys=(star1_y-star2_y)/rasts;//проекция вектора соединяющего положение звезд на ось y
ezs=(star1_z-star2_z)/rasts;//проекция вектора соединяющего положение звезд на ось z
F1x=(G*m1*M*ex1)/((planet_x-star1_x)*(planet_x-star1_x));//проекция силы гравитации между платнетой и звездой 1 на ось х
F1y=(G*m1*M*ey1)/((planet_y-star1_y)*(planet_y-star1_y));//проекция на ось y
F1z=(G*m1*M*ez1)/((planet_z-star1_z)*(planet_z-star1_z));//проекция на ось z
F2x=(G*m2*M*ex2)/((planet_x-star2_x)*(planet_x-star2_x));//проекция силы гравитации между платнетой и звездой 2 на ось х
F2y=(G*m2*M*ey2)/((planet_y-star2_y)*(planet_y-star2_y));//проекция на ось y
F2z=(G*m2*M*ez2)/((planet_z-star2_z)*(planet_z-star2_z));//проекция на ось z
FSx=(G*m1*m2*exs)/((star1_x-star2_x)*(star1_x-star2_x));//проекция силы гравитации между звездами на ось х
FSy=(G*m1*m2*eys)/((star1_y-star2_y)*(star1_y-star2_y));//проекция на ось y
FSz=(G*m1*m2*ezs)/((star1_z-star2_z)*(star1_z-star2_z));//проекция на ось z
}
void SaveFileStar1(char filename[])//сохранение в файл координат
{
ofstream fout(filename, ios::app);
fout<< setw(20) << star1_x << setw(20) << star1_y << setw(20) << star1_z << "\n";
}
void SaveFileStar2(char filename[])
{
ofstream fout(filename, ios::app);
fout<< setw(20) << star2_x << setw(20) << star2_y << setw(20) << star2_z << "\n";
}
void SaveFilePlanet(char filename[])
{
ofstream fout(filename, ios::app);
fout<< setw(20) << planet_x << setw(20) << planet_y << setw(20) << planet_z << "\n";
}
void Verle()//метод интегрирования Верле
{
double const step=1;//шаг по времени
int time = 1000;//все время в течение которого рассматриваем движение
int amount = (time/step);//сколько шагов нужно сделать в цикле
for (int t=0; t<amount; t++)
{
changeforce(planet_x, planet_y, planet_z, star1_x, star1_y, star1_z, star2_x, star2_y, star2_z);//изменяем скорость
vstar1_x=vstar1_x+(F1x-FSx)*step/m1;//здесь и есть сам метод Верле, FSx*step/m1 это формула для скорости после изменения силы
vstar1_y=vstar1_y+(F1y-FSy)*step/m1;// знак минус потому что вектор F направлен в противоположную сторону
vstar1_z=vstar1_z+(F1z-FSz)*step/m1;
vstar2_x=vstar2_x+(F2x+FSx)*step/m2;//то же самое для второй звезды
vstar2_y=vstar2_y+(F2x+FSy)*step/m2;
vstar2_z=vstar2_z+(F2x+FSz)*step/m2;
vplanet_x+=(-F1x-F2x)/M;//и для планеты
vplanet_y+=(-F1y-F2y)/M;
vplanet_z+=(-F1z-F2z)/M;
star1_x=star1_x+vstar1_x*step+0.5*FSx*step*step/m1;//тот же метод для координат, тут для звезды 1
star1_y=star1_y+vstar1_y*step+0.5*FSy*step*step/m1;
star1_z=star1_z+vstar1_z*step+0.5*FSz*step*step/m1;
star2_x=star2_x+vstar2_x*step+0.5*FSx*step*step/m2;//координаты звезды 2
star2_y=star2_y+vstar2_y*step+0.5*FSy*step*step/m2;
star2_z=star2_z+vstar2_z*step+0.5*FSz*step*step/m2;
planet_x+=vplanet_x*step;//координаты планеты
planet_y+=vplanet_y*step;
planet_z+=vplanet_z*step;
SaveFileStar1("Star1"); //сохранение в файл координат
SaveFileStar2("Star2");
SaveFilePlanet("Planet");
}
}
};
int main()
{
ofstream Star1("Star1", ios::trunc); //чистим файл если там что-то есть
Star1.close(); //закрываем его
ofstream Star2("Star2", ios::trunc);
Star2.close();
ofstream Planet("Planet", ios::trunc);
Planet.close();
TWOSTARS planet;//используем наш класс
planet.Verle();//вызываем в классе функцию которая просчитывает координаты
return 0;
} |
(
Сообщение от anpolol
