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
|
% file: "barnes.erl"
-module(barnes2).
-export([?MODULE/0]).
?MODULE() ->
Stars = create_scenario(1000, 1.0),
R = hd(loop(10,1000.0,Stars,0)),
Str = lists:flatten(io:lib_format("~s", [R])),
{R,Str =:= {1.00000,-1.92269e+4,-1.92269e+4,2.86459e-2,2.86459e-2}}.
create_scenario(N, M) ->
create_scenario0(0, 0, trunc(math:sqrt(N)), M).
create_scenario0(_X, _SN, _SN, _M) ->
[];
create_scenario0(SN, Y, SN, M) ->
create_scenario0(0, Y+1, SN, M);
create_scenario0(X, Y, SN, M) ->
XPos0 = (((20000.0 * 2) / SN) * X) - 20000.0,
YPos0 = (((20000.0 * 2) / SN) * Y) - 20000.0,
Calibrate = ((20000.0 * 2) / SN) / 2,
XPos = XPos0 + Calibrate,
YPos = YPos0 + Calibrate,
[{M, XPos, YPos, 0.0, 0.0} | create_scenario0(X+1, Y, SN, M)].
relpos_to_quadrant(DX, DY) when DX >= 0 ->
if
DY >= 0 -> 0;
true -> 3
end;
relpos_to_quadrant(_, DY) ->
if
DY >= 0 -> 1;
true -> 2
end.
quadrant_to_dx(0, D) -> D;
quadrant_to_dx(1, D) -> -D;
quadrant_to_dx(2, D) -> -D;
quadrant_to_dx(3, D) -> D.
quadrant_to_dy(Q,D) ->
if
Q < 2 -> D;
true -> -D
end.
create_tree(Stars) ->
create_tree0(Stars, empty).
create_tree0([],Tree) ->
Tree;
create_tree0([{M,X,Y,_,_} | Stars], Tree) ->
create_tree0(Stars, insert_tree_element(Tree, M, X, Y, 0.0, 0.0, 20000.0)).
insert_tree_element(empty, M, X, Y, _OX, _OY, _D) ->
{body,M,X,Y};
insert_tree_element({branch,M0,SubTree}, M, X, Y, OX, OY, D) ->
Q = relpos_to_quadrant(X-OX,Y-OY),
D2 = D / 2,
DX = quadrant_to_dx(Q,D2),
DY = quadrant_to_dy(Q,D2),
{branch,M0+M,setelement(Q+1,SubTree,
insert_tree_element(element(Q+1,SubTree),
M, X, Y, OX+DX, OY+DY,D2))};
insert_tree_element({body,M0,X0,Y0},M,X,Y,OX,OY,D) ->
resolve_body_conflict(M,X,Y,M0,X0,Y0,OX,OY,D).
resolve_body_conflict(M0, X0, Y0, M1, X1, Y1, OX, OY, D) ->
T = {empty,empty,empty,empty},
Q0 = relpos_to_quadrant(X0-OX,Y0-OY),
Q1 = relpos_to_quadrant(X1-OX,Y1-OY),
D2 = D / 2,
if
Q0 == Q1 -> DX = quadrant_to_dx(Q0,D2),
DY = quadrant_to_dy(Q1,D2),
{branch,M0+M1,setelement(Q0+1,T,
resolve_body_conflict(M0,X0,Y0,
M1,X1,Y1,
OX+DX,OY+DY,
D2))} ;
true -> {branch,M0+M1, setelement(Q1+1,
setelement(Q0+1,T,{body,M0,X0,Y0}),
{body,M1,X1,Y1})}
end.
compute_acceleration(empty, _, _, _, _, _,L) ->
{{0.0, 0.0},L+1};
compute_acceleration({body,BM,BX,BY}, _D, _OX, _OY, X, Y,L) ->
DX = BX - X,
DY = BY - Y,
R2 = (DX * DX) + (DY * DY),
Divisor = R2 * math:sqrt(R2),
if
Divisor < 0.000001 -> % was: Divisor < ?EPSILON ->
{{0.0, 0.0},L+1};
true ->
Expr = BM / Divisor,
{{DX * Expr, DY * Expr},L+1}
end;
compute_acceleration({branch,M,SubTree}, D, OX, OY, X, Y,L) ->
DX = OX - X,
DY = OY - Y,
R2 = (DX * DX) + (DY * DY),
DD = D*D,
R2_THETA2 = 0.09 * R2, % TRY 2.0 *R2, !!! was: R2_THETA2 = ?THETA2*R2,
if
% Ok to approximate?
DD < R2_THETA2 ->
Divisor = R2 * math:sqrt(R2),
if
Divisor < 0.000001 ->
{{0.0,0.0},L};
true ->
Expr = M / Divisor,
{{DX*Expr, DY*Expr},L+1}
end;
% Not ok to approximate...
true ->
D2 = D / 2,
{{AX0, AY0},L1} = compute_acceleration(element(1,SubTree),
D2, OX + quadrant_to_dx(0,D2),
OY + quadrant_to_dy(0,D2),X,Y,L),
{{AX1, AY1},L2} = compute_acceleration(element(2,SubTree),
D2, OX + quadrant_to_dx(1,D2),
OY + quadrant_to_dy(1,D2),X,Y,L1),
{{AX2, AY2},L3} = compute_acceleration(element(3,SubTree),
D2,OX + quadrant_to_dx(2,D2),
OY + quadrant_to_dy(2,D2),X,Y,L2),
{{AX3, AY3},L4} = compute_acceleration(element(4,SubTree),
D2, OX + quadrant_to_dx(3,D2),
OY + quadrant_to_dy(3,D2),X,Y,L3),
{{AX0+AX1+AX2+AX3, AY0+AY1+AY2+AY3},L4+1}
end.
compute_star_accelerations(_Tree,[],L) ->
{[],L};
compute_star_accelerations(Tree,[{_,X, Y,_,_}|Stars],L) ->
{A,AL} = compute_acceleration(Tree, 20000.0, 0.0, 0.0, X, Y,L),
{B,BL} = compute_star_accelerations(Tree, Stars,AL),
{[A | B],BL}.
compute_next_state([],_,_) ->
[];
compute_next_state([{M,X,Y,VX,VY}|Stars],[{AX,AY}|Accs],Time) ->
VX0 = VX + (AX * Time),
VY0 = VY + (AY * Time),
[{M,X+(VX*Time),Y+(VY*Time),VX0,VY0} | compute_next_state(Stars,Accs,Time)].
advance_time(Time,Stars,L) ->
Tree = create_tree(Stars),
{Acc,NL} = compute_star_accelerations(Tree, Stars,L),
{compute_next_state(Stars, Acc, Time),NL}.
loop(0,_Time,Stars,_L) ->
Stars;
loop(N,Time,Stars,L) ->
{NS,NL} = advance_time(Time,Stars,L),
loop(N-1,Time,NS,NL).
|
