Numerical integration seems odd

[won3d]

https://github.com/pablo-mayrgundter/celestiary/blob/5c8e5f2b32272cf07a95f450987bb16596acbba8/js/Galaxy.js#L116-L124

What are lines 118 and 119 for?

On the plus side, lines 120-122 seem to be doing semi-implicit Euler integration, which is a reasonable choice. Basically, the order you accumulate velocity and position matters. Doing:

position += velocity
velocity += acceleration

Isn't as stable as vanilla Euler:

velocity += acceleration
position += velocity

To be clear, that seems to be what you are already doing, but you should note that the order of accumulation matters.

A bonus thought:

A lot of physical simulations will do something slightly different called Verlet integration. I think it has some nicer properties than the Euler variants without being too much more complicated. You could even argue that it is conceptually simpler because it doesn't explicitly even use velocity. This is nice because a lot of more simulations will have constraints on positions, and when you go to enforce them you have to make sure your velocities are also corrected. Not having explicit velocities avoids that .Or maybe you want to draw streaks behind the particles, where you'll need the old positions. But most importantly, Verlet is just a more accurate numerical integrator.

Verlet works by reconstructing velocity from positions. You want to do something like:

pos_after = pos_now + vel + accel

You can recover vel with finite differences:

vel == pos_now - pos_before

When you substitute you get:

pos_after = 2*pos_now - pos_before + accel

In practice, pos_old and pos_next can be the same storage -- you're effectively double-buffering. In the context of ThreeJS, I don't know if that's a good thing or a bad thing. It is possible that ThreeJS tightly couples GPU and CPU vertex buffer storage layouts, which makes updating a buffer not as simple as just setting a flag. There are lots of ways around it but I'd have to know how ThreeJS worked better, and what you would want to accomplish.

[pablo-mayrgundter]

"What are lines 118 and 119 for?" Ah darn, I meant that to be inside an "if (this.first) {}" block. In other words, set the particles spinning normal and proportional to their inward force. I think I added that, saw it worked and forgot it should only be for the first step. Of course, it gives me exactly the effect I want, but for the wrong reasons. Now that I have it in the conditional, I'm having a much harder time balancing the forces. I suppose this is where things actually need to be correct! :) Thanks for the explanation of Verlet integration. I'll also try that.

…

On Tue, Jan 12, 2021 at 11:47 AM won3d ***@***.***> wrote: https://github.com/pablo-mayrgundter/celestiary/blob/5c8e5f2b32272cf07a95f450987bb16596acbba8/js/Galaxy.js#L116-L124 What are lines 118 and 119 for? On the plus side, lines 120-122 seem to be doing semi-implicit Euler integration, which is a reasonable choice. Basically, the order you accumulate velocity and position matters. Doing: position += velocity velocity += acceleration Isn't as stable as vanilla Euler: velocity += acceleration position += velocity To be clear, that seems to be what you are already doing, but you should note that the order of accumulation matters. A bonus thought: A lot of physical simulations will do something slightly different called *Verlet integration*. I think it has some nicer properties than the Euler variants without being too much more complicated. You could even argue that it is conceptually simpler because it doesn't explicitly even use velocity. This is nice because a lot of more simulations will have constraints on positions, and when you go to enforce them you have to make sure your velocities are also corrected. Not having explicit velocities avoids that .Or maybe you want to draw streaks behind the particles, where you'll need the old positions. But most importantly, Verlet is just a more accurate numerical integrator. Verlet works by reconstructing velocity from positions. You want to do something like: pos_after = pos_now + vel + accel You can recover vel with finite differences: vel == pos_now - pos_before When you substitute you get: pos_after = 2*pos_now - pos_before + accel In practice, pos_old and pos_next can be the same storage -- you're effectively double-buffering. In the context of ThreeJS, I don't know if that's a good thing or a bad thing. It is possible that ThreeJS tightly couples GPU and CPU vertex buffer storage layouts, which makes updating a buffer not as simple as just setting a flag. There are lots of ways around it but I'd have to know how ThreeJS worked better, and what you would want to accomplish. — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <https://github.com/pablo-mayrgundter/celestiary/issues/17>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAS5V36EX2QHNEBTQMSMCR3SZSDMVANCNFSM4V7SDWFQ> .

-- *interweb homepage <http://sites.google.com/site/pablomayrgundter>*

[won3d]

Instead of if (this.first) {} can't you set it up when the positions get initialized?

I'd save the Verlet stuff for later if it gives you any trouble. Semi-implicit Euler isn't too bad and you'll want to get things working before you start making structural changes, however minor.

[pablo-mayrgundter]

Yeah, but I think I need the accelerations to set the initial velocities, e.g. using this: https://en.wikipedia.org/wiki/Orbital_speed#Mean_orbital_speed So just have it here until I get them working right

…

On Wed, Jan 13, 2021 at 8:55 PM won3d ***@***.***> wrote: Instead of if (this.first) {} can't you set it up when the positions get initialized? I'd save the Verlet stuff for later if it gives you any trouble. Semi-implicit Euler isn't too bad and you'll want to get things working before you start making structural changes, however minor. — You are receiving this because you commented. Reply to this email directly, view it on GitHub <https://github.com/pablo-mayrgundter/celestiary/issues/17#issuecomment-759891877>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAS5V3YRDW6K7CIU7XLRNGTSZZMLRANCNFSM4V7SDWFQ> .

-- *interweb homepage <http://sites.google.com/site/pablomayrgundter>*

[won3d]

It looks like the "a" in that equation isn't acceleration but axis-length (aka radius)?

[pablo-mayrgundter]

I moved the init code fully into the ctor as suggested. I have a few test points setup for debug and it looks good by eye, but can't find a tangential velocity that works for both of the inner two points. If it works for them, it doesn't work for the outer, and vice-versa. I've been reading the wiki page linked in the comments and some solutions online. I think the idea is the tangential velocity needs to be set to balance the sum effect of the inward accelerations by the time the point gets to pi/2, i.e. if it starts at 1,0, it needs to have velocity 0,1 per OrbitPeriod/4. Like I said, I can tune this manually, but haven't solved it analytically. Will try that next unless you can think of something obv.

…

On Thu, Jan 14, 2021 at 8:59 PM won3d ***@***.***> wrote: It looks like the "a" in that equation isn't acceleration but axis-length (aka radius)? — You are receiving this because you commented. Reply to this email directly, view it on GitHub <https://github.com/pablo-mayrgundter/celestiary/issues/17#issuecomment-760613066>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAS5V3Y5OM5BPI5GWSBO3L3SZ6VSPANCNFSM4V7SDWFQ> .

-- *interweb homepage <http://sites.google.com/site/pablomayrgundter>*

[won3d]

Oh I see, yeah that makes sense! I don't think I ever thought to do it this way, even though in retrospect this is clearly the right way to do it. I always just set velocity as a function of radius, and kind of just approximated what the expected inward acceleration ought to be. You already have a full simulation function, no need to try to do anything analytically, I'd say.

As a tangent, there's something that almost applies: Gauss's Shell Theorem: https://en.wikipedia.org/wiki/Shell_theorem

If you have some force that is inverse-square, like gravity or charge, then spherically symmetric distributions have a nice simplification: a particle can ignore all of the mass that is further than the center of mass than itself. Unfortunately, I don't think this applies to mass distributions around a disk. So some analog exists, but I don't think it would be as nice. But I could be wrong! It would be interesting to plot the results of the force magnitude vs. radius.