Well I've been dropping a few images into the megathread over at shatters.net, but I thought I'd drop a bit more info here.

The source dataset is the Catalog of field W UMa systems, and I am using a Perl script to generate models of the Roche equipotentials. The model I'm using is a first approximation - two point masses in orbits with semimajor axis worked out from Kepler's third law, masses and orbital period. Obviously since the stars fill their Roche lobes and are highly distorted this cannot be accurate, but I don't really have the resources to solve for the actual structure, taking into account the two fusion cores, the magnetohydrodynamics, etc. etc.

This is going to be the largest add-on, bytewise, that I have done. At present, I have the CMOD files (in ASCII form, I'm trying to get a working copy of cmodfix to convert the models to binary format before I release the add-on) and the .stc file which would enable it to work, however the code I'm using needs to be made more robust before I am satisfied.

Here's a screenshot showing six of the 43 systems that will be included...

Yes, Andrew, I followed that already in shatters with interest ! Very nice.

So it might seem that soon we will have a third category of binary systems in Celestia!

It would be perhaps good if you could write a few introductory sentences for people who wouldn't know since childhood what Roche lobes are and how safely they might be calculated...

How do you fix the parameters in these models?

Fridger

Ok, some more info on Roche lobes and W UMa systems for the uninitiated...

Start off with a binary star system with two stars in circular orbits around their centre of mass, treating each star as a point mass because it makes things a lot easier. Then switch to a reference frame that rotates with the binary stars. There are now three forces acting on a particle: the gravitational attraction of each star, and a centrifugal force. (Shock horror I said the "c" word... nevertheless since we're dealing with a rotating reference frame, it's perfectly ok, and less of the inertial frame chauvinism please.)

A star's Roche lobe is the teardrop shaped region of space where the net force results in the matter being gravitationally bound to that star. The Roche lobes of the two stars meet at the Lagrange-1 point (at L1 the forces balance so that a point mass would not move with respect to the two stars, however this is an unstable equilibrium) In the image above, 44i Boötis is depicted with both stars just filling the Roche lobe and touching at Lagrange-1*. If the star itself expands so much that it overflows its Roche lobe, it transfers matter through the Lagrange-1 point to the other star, resulting in lots of interesting phenomena.

In a W UMa system, both stars overflow their Roche lobes, resulting in them filling the figure-of-8 shape volumes depicted in the illustration... each star has its own fusion core but they have a shared photosphere. If the stars expand further, matter eventually gets transferred out of the system through the Lagrange-2 point.

Regarding the parameters in the models, I calculate the semimajor axis from the masses and the orbital period. (As I said, this is a first approximation since the system fairly obviously the stars are not point masses!) The value of the photosphere equipotential V is calculated from the fill-out factor, which is given by f = (V - V_1) / (V_2 - V_1), where V_1 and V_2 are the potentials at L1 and L2 respectively. The fill-out factors are listed in the catalogue.

More details on how to calculate the coordinates of the surface after I improve my algorithms...

* The properties of 44i Boötis are somewhat suspicious: the W UMa binary is itself the dimmer component of a wide binary with a period of around 200 years, which makes the measurements difficult.

Tx Andew for this introduction!

Well I promised more details on how to calculate the surface. For the following, I am going to take the convention where the rotation axis is the z-axis, and the two stars are located along the x-axis. This differs from the convention used in the CMOD files where the rotation axis is the y-axis, but I find it easier to work with.

The first stage is to calculate the Lagrange points L1 and L2. This problem ends up requiring the solution to quintic equations (insert discussion about trying to get analytic solutions to polynomials of degree >4 here), so I use the Newton-Raphson method to calculate the solution numerically. Once I have the positions of these points I can calculate the inner and outer critical potentials, and combine them with the fillout factor to get the surface potential. I then calculate the extents of the surface along the x-axis. Newton-Raphson doesn't work so well here (you have to be very careful about the starting point), so I use interval bisection which is slower but more reliable.

To find out the coordinates of points on the surface I work in cylindrical coordinates (x, rho, theta) where rho and theta are modified from their usual definitions to take account of the fact that the cylinder's axis is the x axis. The problem then is finding rho given x and theta. Since there are possibly two solutions and we want the inner one, I use a modified form of interval bisection to find the surface, taking advantage of the form of the potential to figure out which side of the potential barrier the tested point is on. To save time I exploit the mirror symmetry of the system about the y and z axes.

And in other news, the add-on is basically finished in terms of functionality, still need to do a bit of work on the README, etc. etc.

Andrew,

great that you're getting into more detail now!

ajtribick wrote:

Well I promised more details on how to calculate the surface. For the following, I am going to take the convention where the rotation axis is the z-axis, and the two stars are located along the x-axis. This differs from the convention used in the CMOD files where the rotation axis is the y-axis, but I find it easier to work with. The first stage is to calculate the Lagrange points L1 and L2. This problem ends up requiring the solution to quintic equations (insert discussion about trying to get analytic solutions to polynomials of degree >4 here), so I use the Newton-Raphson method to calculate the solution numerically.

Technically, this would be most easily done and investigated by means of Maple or Mathematica. Do you use one of those?

The main question would be, how you make sure that the particular numerical solution corresponds to the right one among several (5) existing ones!? Do you do this via prescribing admissible ranges for the variables in N-R?

Next, are the equations non-linear? If yes, then the situation with solution ambiguities is even worse.

Quote:

Once I have the positions of these points I can calculate the inner and outer critical potentials, and combine them with the fillout factor to get the surface potential. I then calculate the extents of the surface along the x-axis. Newton-Raphson doesn't work so well here (you have to be very careful about the starting point), so I use interval bisection which is slower but more reliable.

Again, in its general fsolve(...) package, Maple automatically calls various improvement alternatives in case of worsening convergence...

Quote:

To find out the coordinates of points on the surface I work in cylindrical coordinates (x, rho, theta) where rho and theta are modified from their usual definitions to take account of the fact that the cylinder's axis is the x axis. The problem then is finding rho given x and theta. Since there are possibly two solutions and we want the inner one, I use a modified form of interval bisection to find the surface, taking advantage of the form of the potential to figure out which side of the potential barrier the tested point is on. To save time I exploit the mirror symmetry of the system about the y and z axes.

Another physicist's question at this point would be about symmetries of the problem, typically related to conservation laws (angular momentum, energy, ...). Symmetries can help tremendously to break complications down...

Quote:

And in other news, the add-on is basically finished in terms of functionality, still need to do a bit of work on the README, etc. etc.

Looking forward to it!

Fridger

Would like to add my thanks for the explanations here too.
Very interesting to the Brain-Dead, and sure to be
fascinating in Celestia.

No I don't have access to Maple or Mathematica. However for the N-R, a quick bit of curve sketching and gnuplot showed that N-R is pretty well behaved for finding the Lagrange points, and it is fairly easy to find a good starting point for the iterations that won't lead to things ending up in a singularity. Similarly for the setup for the other iterations. Maybe not the most rigorous way to do it but as far as I can tell it works.

I estimate I'll have the whole thing finished next weekend - need to modify the .stc file parser a bit (as far as I can tell no .stc files need to be parsed, but it's better to future-proof the code)

Apologies about the delays... unfortunately I haven't managed to get as much done as I thought I would. But rest assured this is going to be released soon.

ajtribick wrote:

Apologies about the delays...

Hey don't worry! if someone is not happy, just send it to me!

I'm pleased to announce this has now been released. Go here for information, some pretty pictures and the download.