-
Notifications
You must be signed in to change notification settings - Fork 0
/
tlengt_classic.c
103 lines (89 loc) · 2.24 KB
/
tlengt_classic.c
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
/**
# Length of a loop of particles
`plength` is more accurate than simple distance.
The length of two curves is estimated from the particles placed on those
curves.
~~~gnuplot
set key left outside
set size ratio -1
plot 'log' w l lw 2 t 'curve 1', 'circ' w l lw 2 t 'curve 2'
~~~
The estimates are compared against Wolfram-Alpha's reference value
~~~gnuplot Error Curve 1
reset
set size square
set logscale xy
set xr [50: 25600]
set grid
set xlabel 'Particles'
set ylabel 'Error in length'
plot 'out' , 'out' u 1:3, 50000000*x**(-4), 1000*x**(-2)
~~~
The higher-order method is less robust for unresolved features.
~~~gnuplot Error Curve 2
reset
set size square
set logscale xy
set xr [1.5: 256]
set grid
set xlabel 'Particles'
set ylabel 'Error in length'
plot 'out' u 4:5 , 'out' u 4:6, 50*x**(-4), 10*x**(-2)
~~~
But only needs a few point to resolve a circle: The 6-point estimate is less than 1% less off.
*/
#define X(T) (sin(T) + 0.5*erf(2*cos(3*T)))
#define Y(T) (cos(T) + 0.4*exp(2*sin(4*T)))
#define ANS (27.6374980235184900345111)
#include "particle.h"
#include "run.h"
// Simple length:
double length (Particles p) {
double lenp = 0;
coord p1 = {0,0}, pp = {0};
foreach_particle_in(p) {
if (p1.x == 0 && p1.y == 0)
p1 = (coord){x, y};
else {
double d = 0;
foreach_dimension()
d += sq(p().x - pp.x);
lenp += sqrt(d);
}
pp = (coord){x, y};
}
double d = 0;
foreach_dimension()
d += sq(p1.x - pp.x);
lenp += sqrt(d);
return lenp;
}
int np;
int main() {
for (np = 100; np <= 12800; np*= 2)
run();
}
event init (t = 0) {
Particles part1 = new_particles(np);
Particles part2 = new_particles(np/33);
foreach_particle_in(part1) {
double T = j*2.*pi/(np);
p().x = X(T);
p().y = Y(T);
}
foreach_particle_in(part2) {
double Theta = j*2.*pi/(np/33);
p().x = cos(Theta);
p().y = sin(Theta);
}
printf ("%d %g %g %d %g %g\n", np, fabs(ANS-plength(part1)),
fabs(ANS-length(part1)), np/33, fabs(2*pi-plength(part2)),
fabs(2*pi-length(part2)));
if (np == 3200) {
foreach_particle_in(part1)
fprintf (stderr, "%g %g\n", x, y);
FILE * fp = fopen("circ", "w");
foreach_particle_in(part2)
fprintf (fp, "%g %g\n", x, y);
}
}