Новичок
Профиль
Группа: Участник
Сообщений: 2
Регистрация: 17.4.2012
Репутация: нет
Всего: нет
|
Знаю, с такими просьбами уже задолбали( Помогите, пожалуйста!  У меня есть программа, в которой я задаю начальные координаты объекта(x, y, z, v_x, v_y, v_z), задаю интервал времени (T), и программа выдает координаты, которые принимает тело спустя заданное время. (Речь идет о движении звезд в Галактике) Программа работает правильно. Всё хорошо. | Код |
#include <iostream>
#include <cmath>
using namespace std;
void Runge_Kutta (int, double, double *, void (*f));
double dphi_dx(double, double, double);
double dphi_dy (double, double, double);
double dphi_dz (double,double, double);
double dphi_dr (double, double, double);
void diff_equi (int, double*);
void star(double ,double* , double , double , double, double, double, double );
int main () {
int n=6, i;
double input[6], x, y, z, v_x, v_y, v_z;
double T=1e7, a;
x=8.3;
y=0;
z=0;
v_x=0;
v_y=2.3e-7;
v_z=0;
input[0] = x;
input[1] = y;
input[2] = z;
input[3] = v_x;
input[4] = v_y;
input[5] = v_z;
star(T, &input[0], x, y, z, v_x, v_y, v_z );
for (i=0;i<6;i++)
cout<<input[i]<<endl;
cin>>a;
return 0; }
void Runge_Kutta (int n, double T, double * input, void (*f)(int n, double * in_fun)) {
double result [6], h;
double quant_it;
double frac_quant, step = 1000;
frac_quant = modf (T/step, &quant_it); //
if (quant_it < 0) {
quant_it *= -1;
h = -step; }
else
h = step;
double k_1[n], k_2[n], k_3[n], k_4[n];
double vect_for_sent[n];
double t = 0;
for (int i = 0; i < quant_it; i++) {
for (int k = 0; k < n; k++)
vect_for_sent[k] = input[k];
f(n, &vect_for_sent[0]);
for (int k = 0; k < n; k++) {
k_1[k] = vect_for_sent[k];
vect_for_sent[k] = input[k] + 0.5 * h * k_1[k];
}
f(n, &vect_for_sent[0]);
for (int k = 0; k < n; k++) {
k_2[k] = vect_for_sent[k];
vect_for_sent[k] = input[k] + 0.5 * h * k_2[k];
}
f(n, &vect_for_sent[0]);
for (int k = 0; k < n; k++) {
k_3[k] = vect_for_sent[k];
vect_for_sent[k] = input[k] + h * k_3[k];
}
f(n, &vect_for_sent[0]);
for (int k = 0; k < n; k++) {
k_4[k] = vect_for_sent[k];
input[k] += 1./6. * h * (k_1[k] + 2*k_2[k] + 2*k_3[k] + k_4[k]);
}
t+=h;
}
quant_it = 1;
h = step * frac_quant;
for (int i = 0; i < quant_it; i++) {
for (int k = 0; k < n; k++)
vect_for_sent[k] = input[k];
f(n, &vect_for_sent[0]);
for (int k = 0; k < n; k++) {
k_1[k] = vect_for_sent[k];
vect_for_sent[k] = input[k] + 0.5 * h * k_1[k];
}
f(n, &vect_for_sent[0]);
for (int k = 0; k < n; k++) {
k_2[k] = vect_for_sent[k];
vect_for_sent[k] = input[k] + 0.5 * h * k_2[k];
}
f(n, &vect_for_sent[0]);
for (int k = 0; k < n; k++) {
k_3[k] = vect_for_sent[k];
vect_for_sent[k] = input[k] + h * k_3[k];
}
f(n, &vect_for_sent[0]);
for (int k = 0; k < n; k++) {
k_4[k] = vect_for_sent[k];
input[k] += 1./6. * h * (k_1[k] + 2*k_2[k] + 2*k_3[k] + k_4[k]);
}
t+=h;
}
}
double G = 2.2608e-57, M_sol = 2e33;
double M_dh = 1.45e+11*M_sol, M_b = 9.3e+9*M_sol, M_n = 1.e+10*M_sol;
double beta_1 = 0.4, beta_2 = 0.5, beta_3 = 0.1;
double h_1 = 0.325, h_2 = 0.090, h_3 = 0.125;
double a_G = 2.4;
double b_dh = 5.5, b_b = 0.25, b_n = 1.5;
double dphi_dx (double x, double y, double z) {
double res;
res = (M_dh*x*G)/pow(pow(a_G+beta_3*sqrt(z*z+h_3*h_3)+beta_2*sqrt(z*z+h_2*h_2)+beta_1*sqrt(z*z+h_1*h_1),2)+y*y+x*x+b_dh*b_dh, 3./2.) + (M_b*x*G)/pow(y*y+x*x+b_b*b_b, 3./2.)+ (M_n*x*G)/pow(y*y+x*x+b_n*b_n, 3./2.);
return res;
}
double dphi_dy (double x, double y, double z) {
double res;
res = (M_dh*y*G)/pow(pow(a_G+beta_3*sqrt(z*z+h_3*h_3)+beta_2*sqrt(z*z+h_2*h_2)+beta_1*sqrt(z*z+h_1*h_1),2)+y*y+x*x+b_dh*b_dh, 3./2.) + (M_b*y*G)/pow(y*y+x*x+b_b*b_b, 3./2.)+ (M_n*y*G)/pow(y*y+x*x+b_n*b_n, 3./2.);
return res;
}
double dphi_dz (double x, double y, double z) {
double res;
res = (M_dh*((beta_3*z)/sqrt(z*z+h_3*h_3)+(beta_2*z)/sqrt(z*z+h_2*h_2)+(beta_1*z)/sqrt(z*z+h_1*h_1))*(a_G+beta_3*sqrt(z*z+h_3*h_3)+beta_2*sqrt(z*z+h_2*h_2)+beta_1*sqrt(z*z+h_1*h_1))*G)/pow(pow(a_G+beta_3*sqrt(z*z+h_3*h_3)+beta_2*sqrt(z*z+h_2*h_2)+beta_1*sqrt(z*z+h_1*h_1),2)+y*y+x*x+b_dh*b_dh, 3./2.);
return res;
}
double dphi_dr (double x, double y, double z) {
double res;
double r = sqrt(x*x+y*y);
res = (M_dh*r*G)/pow(pow(a_G+beta_3*sqrt(z*z+h_3*h_3)+beta_2*sqrt(z*z+h_2*h_2)+beta_1*sqrt(z*z+h_1*h_1),2)+r*r+b_dh*b_dh, 3./2.) + (M_b*r*G)/pow(r*r+b_b*b_b, 3./2.) + (M_n*r*G)/pow(r*r+b_n*b_n, 3./2.);
return res;
}
void diff_equi (int n, double * input) {
double result [6];
result[3] = - dphi_dx(input[0], input[1], input[2]);
result[4] = - dphi_dy(input[0], input[1], input[2]);
result[5] = - dphi_dz(input[0], input[1], input[2]);
result[0] = input [3];
result[1] = input [4];
result[2] = input [5];
for (int i = 0; i < n; i++)
input[i] = result[i];
}
void star(double T,double* input, double x, double y, double z, double v_x, double v_y, double v_z ) {
double result [6];
result [0] = x;
result [1] = y;
result [2] = z;
result [3] = v_x;
result [4] = v_y;
result [5] = v_z;
//cout<<"A neutron star is moving to "<<T<<endl;
Runge_Kutta (6, T, &input[0], &diff_equi);
x = result [0];
y = result [1];
z = result [2];
v_x = result [3];
v_y = result [4];
v_z = result [5];
}
|
Но, увы, задача поменялась. Для моделирования такого движения мне нужно множество точек(где-то 1000) Мне надо, чтобы программа выдавала координаты через равноотстоящие промежутки времени(каждую 1000 лет), причем выводила их все столбцами (столбец T, столбец x, ..) Была бы очень признательна за помощь. Астрономию обожаю, а программирование не дается( Как я понимаю, надо внести изменения в main , но что-то у меня не получается:( |
Я астроном, а не программист(, застряла с прогой 
Загрузка ...
