Моделирование молекул тела с помощью потенциалов Леннарда-Джонса, Ми или Морзе. Неконсервативная сила - простейшая диссипация(вязкое трение с коэффициентом)

  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
#include <iostream>
#include <time.h>
#include <cmath>
const int type = 3; //(вид потенциала) 1 - Леннард-Джонс, 2 - Ми, 3 - Морзе
const int N = 5; //количество точек
const double T = 1.; //шаг интегрирования
const double D = .1;
const double massa = 1.;
const double a = .2;
const double B = .5;
const double m = 2.;
const double n = 1.;
const double alpha = 1.;
typedef struct Point { //описание точки
double x, y; //положение точки
double v_x, v_y; //скорость точки
} Point;
typedef struct Speed { //описание скорости точки
double v_x, v_y;
} Speed;
typedef struct Place { //описание положения точки
double x, y;
} Place;
double Random_double(double min, double max) {
return min + ((double)rand() / RAND_MAX)*(max-min); //возвращает рандомное число от min до max
}
double r(Place a, Place b) {
return sqrt((b.x-a.x)*(b.x-a.x) + (b.y-a.y)*(b.y-a.y));
}
double r(Speed a, Speed b) {
return sqrt((b.v_x-a.v_x)*(b.v_x-a.v_x) + (b.v_y-a.v_y)*(b.v_y-a.v_y));
}
Speed Lennard_Djons(Place buf, Place tmp) {
Speed res;
double res_d = 12*D/(massa*a) * ( pow(a/r(tmp, buf), 13.) - pow(a/r(tmp, buf), 7.) );
res.v_x = res_d * (tmp.x - buf.x);
res.v_y = res_d * (tmp.y - buf.y);
return res;
}
Speed Mi(Place buf, Place tmp) {
Speed res;
double res_d = (n*m*D)/((n-m)*a) * ( pow(a/r(tmp, buf), n+1) - pow(a/r(tmp, buf), m+1) );
res.v_x = res_d * (tmp.x - buf.x);
res.v_y = res_d * (tmp.y - buf.y);
return res;
}
Speed Morze(Place buf, Place tmp) {
Speed res;
double res_d = 2*alpha*D * (exp(-2*alpha*(r(buf,tmp)-a)) - exp(-alpha*(r(buf,tmp)-a)));
res.v_x = res_d * (tmp.x - buf.x);
res.v_y = res_d * (tmp.y - buf.y);
return res;
}
Speed Fau(Speed w_k, Speed w_n, Place buf, Place tmp) {
Speed res;
double res_d = B * r(w_k, w_n);
res.v_x = res_d * (tmp.x - buf.x);
res.v_y = res_d * (tmp.y - buf.y);
return res;
}
void next_Point(Point *r_t, int i, Point *r_t_T) {
Speed w;
w.v_x = 0.;
w.v_y = 0.;
Place tmp;
tmp.x = r_t[i].x;
tmp.y = r_t[i].y;
for (int j = 0; j < N; j++) {
Place buf;
buf.x = r_t[j].x;
buf.y = r_t[j].y;
Speed tmp_w;
tmp_w.v_x = r_t[i].v_x;
tmp_w.v_y = r_t[i].v_y;
if ((tmp.x != buf.x) || (tmp.y != buf.y)) {
Speed w_tmp;
switch (type) {
case 1: w_tmp = Lennard_Djons(buf, tmp); break;
case 2: w_tmp = Mi(buf, tmp);
case 3: w_tmp = Morze(buf, tmp);
default: w_tmp = Lennard_Djons(buf, tmp);
}
w.v_x += w_tmp.v_x;
w.v_y += w_tmp.v_y;
Speed buf_w;
buf_w.v_x = r_t[j].v_x;
buf_w.v_y = r_t[j].v_y;
w_tmp = Fau(buf_w, tmp_w, buf, tmp);
w.v_x += w_tmp.v_x;
w.v_y += w_tmp.v_y;
}
}
Point tmp_p;
tmp_p.x = 2*r_t[i].x - r_t_T[i].x + w.v_x*T*T;
tmp_p.y = 2*r_t[i].y - r_t_T[i].y + w.v_y*T*T;
r_t_T[i].x = r_t[i].x;
r_t_T[i].y = r_t[i].y;
r_t[i].v_x = w.v_x;
r_t[i].v_y = w.v_y;
r_t[i].x = tmp_p.x;
r_t[i].y = tmp_p.y;
}
void next(Point *r_t, Point *r_t_T) {
for (int i = 0; i < N; i++)
next_Point(r_t, i, r_t_T);
}
int main(int argc, const char * argv[]) {
Point *r_t_T = new Point[N]; //положение точек в момент r(t-T)
Point *r_t = new Point[N]; //положение точек в момент r(t)
srand((unsigned)time(NULL));
for (int i = 0; i < N; i++) {
std::cout << "Введите координаты(2) для точки " << i+1 << ": ";
std::cin >> r_t_T[i].x >> r_t_T[i].y;
r_t_T[i].v_x = Random_double(-0.1, 1.);
r_t_T[i].v_y = Random_double(-0.1, 1.);
r_t[i] = r_t_T[i];
}
int c;
do {
next(r_t, r_t_T);
for (int i = 0; i < N; i++)
std::cout << r_t[i].x << " " << r_t[i].y << std::endl;
std::cout << "Для продолжения нажмите любую клавишу кроме точки(.)" << std::endl;
c = getchar();
} while (c != '.');
return 0;
}