Gravity Simulator http://www.orbitsimulator.com/cgi-bin/yabb/YaBB.pl General >> Discussion >> Eccentricity in binaries - simulation 2 http://www.orbitsimulator.com/cgi-bin/yabb/YaBB.pl?num=1162145224 Message started by Mal on 10/29/06 at 10:07:04

 Title: Eccentricity in binaries - simulation 2 Post by Mal on 10/29/06 at 10:07:04 New simulation here:-----------------------Simulation (2) Primary Star: 1 solar mass, 1 solar radius (yellow) Planet: SMA = 5 AU, mass = 1 Jupiter mass, radius = 69,911 km , e=0.000, i=0.000°, long of asc. node = 0°, arg of perifocus = 0°, mean anomaly = 16°. (purple) Companion: SMA = 250 AU, mass = 0.9 solar masses, radius = 626400 km (0.9 solar), e = 0.800, i = 75°, long of asc. node = 0°, arg of perifocus = 0°, mean anomaly = 140°. (red) Timestep = 512 secs.Start time = year 1 1 0 , time 0 0 0 -----------------------Should be quicker to run, hopefully the timestep is low enough to rule out any problems... going to run it for 500,000 years, the kozai cycle should be 335,795 years.

Title: Re: Eccentricity in binaries - simulation 2
Post by frankuitaalst on 10/29/06 at 11:45:44

Mal wrote:
 New simulation here:-----------------------Simulation (2) Primary Star: 1 solar mass, 1 solar radius (yellow) Planet: SMA = 5 AU, mass = 1 Jupiter mass, radius = 69,911 km , e=0.000, i=0.000°, long of asc. node = 0°, arg of perifocus = 0°, mean anomaly = 16°. (purple) Companion: SMA = 250 AU, mass = 0.9 solar masses, radius = 626400 km (0.9 solar), e = 0.800, i = 75°, long of asc. node = 0°, arg of perifocus = 0°, mean anomaly = 140°. (red) Timestep = 512 secs.Start time = year 1 1 0 , time 0 0 0 -----------------------Should be quicker to run, hopefully the timestep is low enough to rule out any problems... going to run it for 500,000 years, the kozai cycle should be 335,795 years.

I observed in my sim that the time to reach the first max. ecc  is appr. 2.5 times the cycle of the Kozai . This could mean that you'll need about 800.000 years to reach max ecc and then the cycle begins .

Frank

 Title: Re: Eccentricity in binaries - simulation 2 Post by Mal on 10/29/06 at 15:41:12 Well this is a bit odd. I'm taking samples every 100 years, and this is what I get for planet inclination and eccentricity:http://www.evildrganymede.net/temp/sim2_512_a.gifThose spikes are not always equidistant - they're 2900 years +/- 100 years apart, with no apparent pattern to the variation (sometimes it's 2800, sometimes it's 2900, sometimes it's 3000). I thought they might be the orbital period of the companion (ie due to close approaches when it's at perihelion?) but if I've done the calculation right the companion's orbital period is 3952 years. That said, I'm not sure what happens if the masses of the primary and compaion are similar, does that shorten the orbital period?EDIT: I got GravSim to calculate the period and it agrees that it's about the same as I calculated above - so I dunno what that 3000 year period is!Note that this is after only 30,000 years or so, so it's still early days.

 Title: Re: Eccentricity in binaries - simulation 2 Post by Mal on 10/29/06 at 18:30:56 More oddness. http://www.evildrganymede.net/temp/sim2_512_b.gifTop show the planet's orbit, bottom shows the companion's. Note that the spikes in eccentricity in both bodies align. There's also sharp dips in the Arg of Peri of the planet too at the same time (and smaller upward kinks in the LAN of the planet).Wonder what's going on??

 Title: Re: Eccentricity in binaries - simulation 2 Post by Mal on 10/29/06 at 19:47:35 Got it!I was watching the simulation around the time of the next encounter (33300 years), and it IS the closest approach of the companion. At 33300 the companion was at closest approach to the star, and the planet was directly between it and the primary. So it's obviously getting a sharp gravitational tug from this, which bumps up its eccentricity, but once the companion has passed the eccentricity isn't as low as it was before the encounter. So the planet eccentricity gets ratcheted higher with each pass. I've only seen this one encounter - I'm guessing the different height of the spikes is because the planet and companion aren't exactly lined up in each encounter?Though I'm still not sure why it's every 2800-3000 years if the orbital period of the companion is longer than that...

Title: Re: Eccentricity in binaries - simulation 2
Post by Tony on 10/29/06 at 20:25:23

Mal wrote:
 Though I'm still not sure why it's every 2800-3000 years if the orbital period of the companion is longer than that...

It's because you've helped uncover a bug in the beta's period formula  :-[

The forumla it uses is the generalized formula for when M1 >> M2:

http://orbitsimulator.com/gravity/dd1.GIF

In this case M1 is not >> M2, and it needs to use M1+M2 for M in the formula.

The real answer for period in of the companion is .... drumroll please.....

2868 years

You can use this formula in Excel to get the correct period.  1.9 represents 1.9 solar masses which is M1 + M2:

=2*3.14159 *SQRT(( K2*1000)^3/(0.000000000066725985*(1.98891691172467E+30*1.9)))

And here's a utility you might find helpful.  It's a javascript period calculator:
http://orbitsimulator.com/gravity/articles/PeriodCalculator.html

Also, you can have Gravity Simulator compute an object's period using the "Rotating Frame Adjustment" window.  Choose the object and it will display its period in seconds.  It gets this one right as it uses M1+M2.  Then just hit cancel so you don't go into Rotating Frame mode.

 Title: Re: Eccentricity in binaries - simulation 2 Post by Mal on 10/29/06 at 22:16:30 Aha, I thought the formula was a bit different when the masses of the bodies were similar. Huzzah for experimentation! :)This is cool, I feel like I'm doing real Science! here ;).

Title: Re: Eccentricity in binaries - simulation 2
Post by Tony on 10/29/06 at 22:23:36

Mal wrote:
 This is cool, I feel like I'm doing real Science! here ;).

You are!  How many people in the world, besides the authors of the quoted papers, do you think know more about the Kozai Mechanism than you?

The real science will be when we can figure out why the formulas break down at large distances.  Does it have to do with the Kozai Mechanism, or with our methods?

Your 750 AU sim was qualitatively correct, but the timing was seriously off.  My 400 AU sim was qualitatively correct, and the timing was roughly correct.  All sims less than 400 AU had good to perfect timings.

 Title: Re: Eccentricity in binaries - simulation 2 Post by Mal on 10/29/06 at 22:31:02 Well I looked at the close approach at 62000 years and the companion and planet weren't lined up then. So it's definitely not repeating the exact configuration each time. I wonder what the inverted spikes are though - sometimes the eccentricity suddenly drops and then shoots up and back down a bit to get to the next step. I'm guessing the planet is on the opposite side of the sun to the companion's pericentre there?So was this happening in the 750 AU sim too, but we couldn't see it because of the timestep? And is this what the Kozai mechanism really is? Because this "tugging at closest approach" thing seems more like a normal direct orbital mechanics situation, I thought the Kozai mechanism was not such a direct thing. I'm wondering if the timings were off in the more distant ones because of the timesteps. For the closer ones we can use shorter timesteps and not die of old age waiting for th results. For the further ones we have to use larger timesteps if we want to see the results in a few days, and that might be messing things up?

Title: Re: Eccentricity in binaries - simulation 2
Post by Tony on 10/29/06 at 22:37:10

Mal wrote:
 So was this happening in the 750 AU sim too, but we couldn't see it because of the timestep?

Not the time step of the simulation, but the time step of the data output.  You're now asking for output to be plotted every 100 years instead of every several thousand years.  So the spikes are just completely skipped over, and even the stair-stepping is just reduced to fuzz on the curve plotted ever few thousand years.  By asking for data every 100 years, you're zooming in on the fuzz.

Title: Re: Eccentricity in binaries - simulation 2
Post by Tony on 10/29/06 at 22:45:10

Mal wrote:
 I'm wondering if the timings were off in the more distant ones because of the timesteps. For the closer ones we can use shorter timesteps and not die of old age waiting for the results. For the further ones we have to use larger timesteps if we want to see the results in a few days, and that might be messing things up?

Patience is a virtue.  But that said, I think the correct approach is to start with short doable sims to try to spot a trend.  Maybe make a graph of how accurate the formula is vs. AU of companion.  Then derive a formula based on the graph and use it to predict the Kozai period for a high AU secondary.  Then perform a sim that might just take a few weeks to run.

Just some random thoughts...
Gravity Simulator takes up 50% of your computer's CPU cycles.  I've asked on coding forums how to increase this, but have never got a solution.  But because of this, you can easily check e-mail, surf the web, etc with gravity simulator running in the background.

Or... you can run 2 sims at the same time with little or no comprimise in the speed vs. running only 1 sim.  And e-mail, etc, isn't impossible, only sluggish, when running 2 instances of Gravity Simulator at the same time.

 Title: Re: Eccentricity in binaries - simulation 2 Post by Mal on 10/30/06 at 07:05:58 156,000 years:http://www.evildrganymede.net/temp/sim2_512_c.gifI hope the eccentricity isn't flattening out yet - we're nearly at halfway through the cycle and it's still nowhere near 0.94... and the inclination is at 32 degrees and still rising.I think what SHOULD be happening is that the eccentricity goes up to 0.94 around 170,000 years, and then the inclination reaches about 35 degrees at the same time and starts going down when eccentricity reaches its peak. I'm not sure if that's going to happen here...

 Title: Re: Eccentricity in binaries - simulation 2 Post by Mal on 10/30/06 at 20:51:28 306,200 years (planet):http://www.evildrganymede.net/temp/sim2_512_d.gifWell this looks a bit promising. The eccentricity peak's at the about the right height (0.934-ish). I presume the next peak is going to be 335,000 years from that point. The LAN's obviously dropped through 0 degrees and looped round to 360 degrees and still dropping, which explains the vertical line there. Though I have no clue whatsoever what is going on with the Inclination (yellow) or Arg of Peri (cyan). The Inc peaks at about 136 degrees and then starts dropping down again.The companion inclination and eccentricity isn't doing much, still just varying slightly around its starting values. I got nothin'. Anyone got any ideas what's going on here?Is there any way we can calculate the relative inclination here?

 Title: Re: Eccentricity in binaries - simulation 2 Post by Tony on 10/30/06 at 23:05:41 I doubt the companion is going to be affected much by the planet.  Though it would be interesting to see what would happen if two stars orbited at 2.5 AU, and a Jupiter-mass companion orbited the pair, out of plane from a distance.  I don't think the formulas account for such a scenerio.

Title: Re: Eccentricity in binaries - simulation 2
Post by Mal on 10/30/06 at 23:10:14

Tony wrote:
 I doubt the companion is going to be affected much by the planet.  Though it would be interesting to see what would happen if two stars orbited at 2.5 AU, and a Jupiter-mass companion orbited the pair, out of plane from a distance.  I don't think the formulas account for such a scenerio.

I wouldn't expect it to be affected by the planet, but I don't think we're seeing the cyclic variation in inclination that we're supposed to be seeing here. But what's weird is that the eccentricity seems to be doing what it's supposed to now.

Could there be something missing in gravity simulator's code here that prevents it from doing the inclination cycle properly? Or is the inclination actually doing the right thing here, but we're just not seeing it because we're not looking at the relative inclination (which presumably is calculated in a more complex way than just "companion inclination minus the planet inclination").

 Title: Re: Eccentricity in binaries - simulation 2 Post by Tony on 10/30/06 at 23:44:41 There's nothing in the code that tells it to do inclination cycles.  The code only tells the objects to follow Newton's laws.  Inclination cycles, just like ecc cycles should be a natural consequence of obeying Newton's laws.Gravity Simulator is giving you all your orbital elements based on the ecliptic, or more specifically, the x-y plane.The inclination you're reading for the planet is from the ecliptic, not from the plane of the secondary's orbit.  And the cycle you're expecting to see is probably based on the plane of the secondary.  But you can't just subtract one from the other.  Converting coordinates from one plane to the other requires some tricky geometry.  (It's probably not that tricky, but off the top of my head I don't know how to do it).What if you tried putting the planet in a 75 degree inclination orbit and the secondary in a 0 degree orbit.  Then you wouldn't have to do any conversions since the secondary would be in the ecliptic, and the inclination numbers you plot would now be the planet's plane with respect to the secondary's plane.  I'd bet you'd get the cycles you were looking for.

Title: Re: Eccentricity in binaries - simulation 2
Post by Mal on 10/30/06 at 23:52:37

Tony wrote:
 There's nothing in the code that tells it to do inclination cycles.  The code only tells the objects to follow Newton's laws.  Inclination cycles, just like ecc cycles should be a natural consequence of obeying Newton's laws.

Ok then.

Quote:
 Gravity Simulator is giving you all your orbital elements based on the ecliptic, or more specifically, the x-y plane.The inclination you're reading for the planet is from the ecliptic, not from the plane of the secondary's orbit.  And the cycle you're expecting to see is probably based on the plane of the secondary.  But you can't just subtract one from the other.  Converting coordinates from one plane to the other requires some tricky geometry.  (It's probably not that tricky, but off the top of my head I don't know how to do it).

I didn't get this last time you mentioned it, but now you say that it's relative to the x-y plane it makes more sense. But still... in this sim the planet starts in the x-y plane and the companion is inclined at 75° to it. As the sim progresses the planet inclination starts increasing, but that's still measured from the x-y plane isn't it? So if at one point the planet inclination is 35° and the companion inclination is 75° doesn't that mean the relative inclination is 40°, since they're both measured relative to the same x-y plane?

Quote:
 What if you tried putting the planet in a 75 degree inclination orbit and the secondary in a 0 degree orbit.  Then you wouldn't have to do any conversions since the secondary would be in the ecliptic, and the inclination numbers you plot would now be the planet's plane with respect to the secondary's plane.  I'd bet you'd get the cycles you were looking for.

Well I can try that next, I'll leave the current sim running for tonight to see where it goes. It's going to stop recording at 500,000 years anyway.

 Title: Re: Eccentricity in binaries - simulation 2 Post by Tony on 10/31/06 at 00:06:39 Only if your longitude of the ascending nodes are equal in both orbits can you simply add or subtract inclinations and get a useful number.  Consider the following:http://www.orbitsimulator.com/gravity/images/inclinations.GIFAll three orbits have inclinations of 30 degrees.  The purple and green orbits have the same longitude of ascending node.  Therefore you can do 30-30=0 and conclude that the purple orbit is inclined 0 degrees to the green orbit.  But the 30-30 trick doesn't work for the orange orbit.   It is clearly 60 degrees inclined from the green and purple orbits. Instead, you need 30+30 = 60.The green and purple orbits have a longitude of ascending node of 90 degrees and the orange orbit has a longitude of 270, 180 degrees away from the purple and green.You may be able to express the inclinations relative to each other in excel with this formula:i1+cos(deltaLAN)*i2Excel likes radians instead of degrees, so...=E2+cos((G2-O2)/180 * 3.14159)*M2Assuming my conversion logic is correct.

 Title: Re: Eccentricity in binaries - simulation 2 Post by Mal on 10/31/06 at 07:30:01 420,000 years (planet):http://www.evildrganymede.net/temp/sim2_512_e.gifSo I guess Arg of Peri, LAN, and Inclination are doing weird looking things because they're not in the same plane as the companion? Maybe it'd look more sensible if they were calculated relative to the companion instead?Here's what the companion's been doing:http://www.evildrganymede.net/temp/sim2_512_comp.gifLAN and Arg of Peri aren't shown, because they're oscilating slightly around 0 degrees - so they're going from 1 to 359 degrees and messing up the graph.

Title: Re: Eccentricity in binaries - simulation 2
Post by Tony on 10/31/06 at 15:33:09

Tony wrote:
 You may be able to express the inclinations relative to each other in excel with this formula:i1+cos(deltaLAN)*i2Excel likes radians instead of degrees, so...=E2+cos((G2-O2)/180 * 3.14159)*M2

This is wrong.  I think the correct way is this:

i3 = acos(sin(i1)sin(i2)cos(deltaLAN)+cos(i1)cos(i2))

and for Excel:

=ACOS(SIN(E2/180*PI())*SIN(M2/180*PI())*COS((G2-O2)/180*PI())+COS(E2/180*PI())*COS(M2/180*PI()))/PI()*180

Assuming:
i1 is in the E column
i2 is in the M column
LAN1 is in the G column
LAN2 is in the O column.

I'll run a sim to verify.

 Title: Re: Eccentricity in binaries - simulation 2 Post by Mal on 10/31/06 at 18:03:48 536,700 years:http://www.evildrganymede.net/temp/sim2_512_f.gifTwo peaks now, but the odd thing is that they're about 218,000 years apart - I thought the calculation said that the cycle should be 335,795 years?!

 Title: Re: Eccentricity in binaries - simulation 2 Post by Tony on 10/31/06 at 18:22:13 If instead of measuring peak-to-peak, your measure start to first peak and multiply by 2, the accuracy improves a little.  I think in ideal conditions (ie, much slower timestep), the eccentricity should have gone all the way back down to 0.Also keep in mind that the formulas in the paper don't use the equal sign.  They use the =~ sign.

 Title: Re: Eccentricity in binaries - simulation 2 Post by Tony on 11/01/06 at 16:40:59 I verified the formula I gave above by running 2 simulations.:1:  planet 2.5 AU, 1 Mj     Secondary 200 AU 0.9 Ms, 75 degrees2:  planet 2.5 AU, 1 Mj, 75 degrees     secondary 200 AU 0.9 MsTime step 32K for 2 million years.With a small modification the formula works flawlessly.  The inclination of the planes of the planet and the secondary are identical when compared to each other instead of against the ecliptic by using the Excel Formula.I edited the formula in the above post by placing a few missing multiplication symbols (*) and converting the answer to degrees.  Again, the Excel formula is:=ACOS(SIN(E2/180*PI())*SIN(M2/180*PI())*COS((G2-O2)/180*PI())+COS(E2/180*PI())*COS(M2/180*PI()))/PI()*180Assuming: i1 is in the E column i2 is in the M column LAN1 is in the G column LAN2 is in the O column.

 Title: Re: Eccentricity in binaries - simulation 2 Post by Mal on 11/01/06 at 17:19:58 789,800 yearsI'm calling it quits here. Top graph is the planet, bottom is the companion. I got no idea what all the properties of the planets' orbit that aren't eccentricity are doing... http://www.evildrganymede.net/temp/sim2_512_planet.gifhttp://www.evildrganymede.net/temp/sim2_512_comp2.gifI'll try those formulas when I get a chance...

Title: Re: Eccentricity in binaries - simulation 2
Post by Tony on 11/01/06 at 20:07:43

Quote:
 I got no idea what all the properties of the planets' orbit that aren't eccentricity are doing...

Your inclination is rising and falling with a period that is twice the period of the eccentricity's.

Your longitude of ascending node is changing with a period approximately equal to the period of your eccentricity.
You might want to create an additional column in your spreadsheet and add 180 to your LAN to center it on the graph so you don't get those vertical lines.  There's nothing special about 0 degrees or 90 degrees, etc.  They're simply degrees from the Vernal Equinox, which is a manmade reference point based on Earth's orbit.  Adding 180 degrees would give you the angle from the Spring Equinox, while still giving you the same trend, but unbroken.

It's a cool looking graph nonetheless.

Your graph for the secondary is very tight.  Ecc only ranges from .79 to .81.  This probably isn't too different than if you ran a simulation with only the star and the secondary, eliminating the planet.  Theoretically, in such a 2-body configuration these values should remain perfectly constant.  But noise generated by the numerical method will give you some fluxuation.  I suspect you're zoomed in on the noise.

Another orbital element you may be interested in is Longitude of Periapsis.  This differs from Argument of Periapsis in that "Argument" means degrees from the ascending node, and "Longitude" means degrees from the Vernal Equinox which is a fixed point.  Gravity Simulator doesn't output it, but its easy to compute.  It's Longitude of Ascending Node + Argument of Periapsis.  So you could make another column in Excel for Longitude of Periapsis.

Title: Re: Eccentricity in binaries - simulation 2
Post by Mal on 11/02/06 at 07:33:00

Tony wrote:
 Your inclination is rising and falling with a period that is twice the period of the eccentricity's.

I don't think it is. It looks like the inclination's period is the same as the eccentricity's, each time there's an ecc peak the inclination starts going in the other direction on the graph. You can see a cycle more clearly on the companion graph (where the bigger spikes are).

Quote:
 Your longitude of ascending node is changing with a period approximately equal to the period of your eccentricity.

Don't think so (see below)

Quote:
 Your graph for the secondary is very tight.  Ecc only ranges from .79 to .81.  This probably isn't too different than if you ran a simulation with only the star and the secondary, eliminating the planet.  Theoretically, in such a 2-body configuration these values should remain perfectly constant.  But noise generated by the numerical method will give you some fluxuation.  I suspect you're zoomed in on the noise.

There's definitely something corresponding to the eccentricity peaks in the companion inc and ecc graphs though (where the biggest spikes are).

Quote:
 Another orbital element you may be interested in is Longitude of Periapsis.  This differs from Argument of Periapsis in that "Argument" means degrees from the ascending node, and "Longitude" means degrees from the Vernal Equinox which is a fixed point.  Gravity Simulator doesn't output it, but its easy to compute.  It's Longitude of Ascending Node + Argument of Periapsis.  So you could make another column in Excel for Longitude of Periapsis.

Is there any reason why GravSim couldn't output this? It's a fairly useful orbital parameter after all.

Here's a graph with the LAN+180° and the Long of Peri too.

http://www.evildrganymede.net/temp/sim2_512_planet2.gif

And the period of the eccentricity is the peak-to-peak or trough-to-trough distance, which looks like it's about 200,000 years. It looks to me on this graph that the period of oscillation for the ecc and inc are about 200,000 years, and the period of oscillation of the LAN and LoP is double that (about 400,000 years). I can't tell much from the Arg of Peri because it's just going all over the place.