ring/eight.html at main · bobjen/ring · GitHub

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
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
<!DOCTYPE html>
<html>
<head>
  <title>32-Ring Dyson Swarm</title>
  <script src="ring.js"></script>
</head>
<body bgcolor="#000000" text="#ffffff" style="font-family:monospace;">

<p><center><canvas id="cv" width="600" height="400" style="border:1px solid #ffffff;"></canvas></center></p>
<pre id="log" style="color:#00ff00;font-size:11px;max-height:300px;overflow:auto;"></pre>

<script>
    document.getElementById("cv").width = window.innerWidth*0.6;
    document.getElementById("cv").height = window.innerWidth*0.4;
(function() {
    var N = 32;
    var rm = 0.001;
    var colors = [];
    for (var ci = 0; ci < N; ci++) {
        var h = (ci * 360 / N) | 0;
        colors.push('hsl(' + h + ',100%,60%)');
    }

    // r and rm are user choices: r sets the orbit amplitude (both the Y and Z
    // extents of the orbit scale with r), rm sets the ring mass.
    // The solver only varies sunR to close the orbit for the given r and rm.
    var r = 0.4;   // orbit radius in YZ — sets the physical size; do not vary in solver

    var BASE = { increment: 0.01, work: 2, points: 8, trail: false };

    // Run a single ring around the oblate sun for one full period, recording its
    // trajectory.  Return 8 equally-time-spaced samples (k*T/8, k=0..7) as the
    // initial conditions for the 8 rings.  This naturally accounts for the ring
    // speeding up near the sun and slowing down at the extremes, unlike equally-
    // spaced angles which assume constant angular speed.
    //
    // States are stored as [y, z, vy_phys, vz_phys] (positions unscaled,
    // velocities unscaled by dividing out inc).  Period detection uses the same
    // vz zero-crossing interpolation as halfOrbitDy.
    function sampleOrbit(sunR) {
        var cosmos = new Cosmos2D(Object.assign({}, BASE, { moons: [
            { y: sunR, z: 0, vy: 0, vz: 0, mass: 1 - N*rm, l: 0, fixed: true },
            { y: Math.sqrt(1 - r*r), z: r, l: 1 - r*r, vy: -r, vz: 0, mass: rm }
        ]}));
        var ring0 = cosmos.rings[1];
        var inc   = cosmos.inc;

        // traj[i] = state after i multistep calls: [y, z, vy, vz] (physical units)
        var traj = [[ring0.p.y, ring0.p.z, ring0.v.y / inc, ring0.v.z / inc]];
        var prevVz = 0, wentNeg = false, wentPos = false, T = 0;

        for (var i = 0; i < 100000; i++) {
            cosmos.multistep();
            var vz = ring0.v.z;
            traj.push([ring0.p.y, ring0.p.z, ring0.v.y / inc, ring0.v.z / inc]);

            if (!wentNeg && vz < 0)             wentNeg = true;
            if ( wentNeg && !wentPos && vz >= 0) wentPos = true;
            if ( wentPos && vz < 0) {
                // Interpolate: prevVz > 0, vz < 0 → crossing at fraction alpha into this step
                var alpha = prevVz / (prevVz - vz);
                T = i + alpha;   // period in step units
                break;
            }
            prevVz = vz;
        }

        // Extend traj to cover a second full orbit for the seam blend below.
        var need = Math.ceil(2 * T) + 2;
        while (traj.length <= need) {
            cosmos.multistep();
            traj.push([ring0.p.y, ring0.p.z, ring0.v.y / inc, ring0.v.z / inc]);
        }

        // Sample N points blended between orbit 1 and orbit 2: ring k gets
        // (k/N)*orbit1[k] + (1-k/N)*orbit2[k], smoothing the seam at k=0/k=N-1.
        var samples = [];
        for (var k = 0; k < N; k++) {
            var w  = k / N;
            var t1 = k * T / N,     i1 = Math.floor(t1), f1 = t1 - i1;
            var t2 = T + k * T / N, i2 = Math.floor(t2), f2 = t2 - i2;
            var p1 = traj[i1], q1 = traj[i1 + 1] || traj[i1];
            var p2 = traj[i2], q2 = traj[i2 + 1] || traj[i2];
            samples.push({
                y:  w*(p1[0]+f1*(q1[0]-p1[0])) + (1-w)*(p2[0]+f2*(q2[0]-p2[0])),
                z:  w*(p1[1]+f1*(q1[1]-p1[1])) + (1-w)*(p2[1]+f2*(q2[1]-p2[1])),
                vy: w*(p1[2]+f1*(q1[2]-p1[2])) + (1-w)*(p2[2]+f2*(q2[2]-p2[2])),
                vz: w*(p1[3]+f1*(q1[3]-p1[3])) + (1-w)*(p2[3]+f2*(q2[3]-p2[3])),
            });
        }
        return samples;
    }

    // Build the full moon list from 8 sampled states + the oblate sun.
    function makeMoons(samples, sunR) {
        var moons = [];
        moons.push({ y: sunR, z: 0, vy: 0, vz: 0,
                     mass: 1 - N*rm, l: 0, fixed: true, color: '#ffff00', radius: 10 });
        for (var i = 0; i < N; i++) {
            moons.push({ y: samples[i].y, z: samples[i].z,
                         vy: samples[i].vy, vz: samples[i].vz,
                         l: 1 - r*r, mass: rm, color: colors[i], radius: 3 });
        }
        return moons;
    }

    // Like sampleOrbit but runs the full N-ring system so ring-ring interactions
    // are included.  Used to iterate samples to self-consistency.
    function sampleOrbitFull(sunR, samples) {
        var cosmos = new Cosmos2D(Object.assign({}, BASE,
                                  { moons: makeMoons(samples, sunR) }));
        var ring0 = cosmos.rings[1];
        var inc   = cosmos.inc;

        var traj = [[ring0.p.y, ring0.p.z, ring0.v.y / inc, ring0.v.z / inc]];
        var prevVz = 0, wentNeg = false, wentPos = false, T = 0;

        for (var i = 0; i < 100000; i++) {
            cosmos.multistep();
            var vz = ring0.v.z;
            traj.push([ring0.p.y, ring0.p.z, ring0.v.y / inc, ring0.v.z / inc]);

            if (!wentNeg && vz < 0)             wentNeg = true;
            if ( wentNeg && !wentPos && vz >= 0) wentPos = true;
            if ( wentPos && vz < 0) {
                var alpha = prevVz / (prevVz - vz);
                T = i + alpha;
                break;
            }
            prevVz = vz;
        }

        var need = Math.ceil(2 * T) + 2;
        while (traj.length <= need) {
            cosmos.multistep();
            traj.push([ring0.p.y, ring0.p.z, ring0.v.y / inc, ring0.v.z / inc]);
        }

        var out = [];
        for (var k = 0; k < N; k++) {
            var w  = k / N;
            var t1 = k * T / N,     i1 = Math.floor(t1), f1 = t1 - i1;
            var t2 = T + k * T / N, i2 = Math.floor(t2), f2 = t2 - i2;
            var p1 = traj[i1], q1 = traj[i1 + 1] || traj[i1];
            var p2 = traj[i2], q2 = traj[i2 + 1] || traj[i2];
            out.push({
                y:  w*(p1[0]+f1*(q1[0]-p1[0])) + (1-w)*(p2[0]+f2*(q2[0]-p2[0])),
                z:  w*(p1[1]+f1*(q1[1]-p1[1])) + (1-w)*(p2[1]+f2*(q2[1]-p2[1])),
                vy: w*(p1[2]+f1*(q1[2]-p1[2])) + (1-w)*(p2[2]+f2*(q2[2]-p2[2])),
                vz: w*(p1[3]+f1*(q1[3]-p1[3])) + (1-w)*(p2[3]+f2*(q2[3]-p2[3])),
            });
        }
        return out;
    }

    // Run the 8-ring system forward until ring 0 reaches the z-minimum (half orbit).
    // Returns Δy = y(T/2) - y(0); zero for a closed symmetric orbit.
    // Uses the sampled orbit for initial conditions and interpolates the vz crossing.
    function halfOrbitDy(sunR) {
        var samples = sampleOrbit(sunR);
        var cosmos  = new Cosmos2D(Object.assign({}, BASE,
                                   { moons: makeMoons(samples, sunR) }));
        var ring0   = cosmos.rings[1];
        var y0      = ring0.p.y;
        var wentNeg = false, prevVz = 0;

        for (var i = 0; i < 50000; i++) {
            cosmos.multistep();
            var vz = ring0.v.z;
            if (!wentNeg && vz < -1e-12) wentNeg = true;
            if ( wentNeg && vz >= 0) {
                var alpha  = -prevVz / (vz - prevVz);
                return ring0.p.y - (1 - alpha) * ring0.v.y - y0;
            }
            prevVz = vz;
        }
        return ring0.p.y - y0;
    }

    // Scan sunR in [0, maxSun] at nScan points to find a sign change,
    // then bisect to precision eps.  sunR is the only free parameter.
    // Returns the bracket midpoint; the bracket [lo,hi] is kept tight enough
    // that further Newton polish is unnecessary (the interpolation above means
    // dy noise is now O(dt²) ≈ 1e-5, so eps=1e-6 is safely below signal).
    function findSunR(maxSun, nScan, eps) {
        var log  = document.getElementById('log');
        var step = maxSun / nScan;
        var prev = { s: 0, dy: halfOrbitDy(0) };
        // log.textContent += 'Scanning sunR  (r=' + r.toFixed(4) + '  rm=' + rm.toFixed(4) + '):\n';
        // log.textContent += '  sunR=0.000000  dy=' + prev.dy.toExponential(4) + '\n';

        for (var i = 1; i <= nScan; i++) {
            var s  = i * step;
            var dy = halfOrbitDy(s);
            // log.textContent += '  sunR=' + s.toFixed(6) + '  dy=' + dy.toExponential(4) + '\n';
            if (prev.dy * dy <= 0) {
                var lo = prev.s, dlo = prev.dy;
                var hi = s;
                // log.textContent += 'Bisecting [' + lo.toFixed(6) + ', ' + hi.toFixed(6) + '] ...\n';
                for (var j = 0; j < 80; j++) {
                    var mid = (lo + hi) * 0.5;
                    var dm  = halfOrbitDy(mid);
                    if (dlo * dm <= 0) { hi = mid; }
                    else               { lo = mid; dlo = dm; }
                    if (hi - lo < eps) break;
                }
                var sunR = (lo + hi) * 0.5;
                // log.textContent += 'Converged: sunR=' + sunR.toFixed(10) +
                //                    '  dy=' + halfOrbitDy(sunR).toExponential(4) + '\n';
                return sunR;
            }
            prev = { s: s, dy: dy };
        }
        // log.textContent += 'No root in [0, ' + maxSun + '] — try different r or rm.\n';
        return 0;
    }

    var sunR    = findSunR(0.5, 20, 1e-6);
    var samples = sampleOrbit(sunR);

    // Refine samples to self-consistency under the full N-ring system.
    // sunR is kept fixed — ring-ring forces are symmetric so the required
    // sunR barely changes, but the sample positions do need correction.
    var lg = document.getElementById('log');
    // lg.textContent += '\nRefining samples with full ' + N + '-ring system:\n';
    for (var iter = 0; iter < 6; iter++) {
        var s2 = sampleOrbitFull(sunR, samples);
        var maxD = 0;
        for (var k = 0; k < N; k++)
            maxD = Math.max(maxD, Math.abs(s2[k].y - samples[k].y),
                                  Math.abs(s2[k].z - samples[k].z));
        samples = s2;
        // lg.textContent += '  iter ' + iter + ': maxDelta=' + maxD.toExponential(3) + '\n';
        if (maxD < 1e-9) break;
    }

    // Animate the result
    ring('cv', Object.assign({}, BASE, {
        background: '#000000',
        scale:    window.innerWidth*0.4,
        ymargin:  20,
        framerate: 20,
        trail:    false,
        traillen: 1,
        fade:     1,
        stop:     true,
        moons:    makeMoons(samples, sunR),
    }));

    document.getElementById('log').textContent +=
        '\nFinal: r=' + r.toFixed(8) + '  rm=' + rm.toFixed(6) +
        '  sunR=' + sunR.toFixed(8) + '\n';
})();
</script>

</body>
</html>