## Plotting effective potential (Lagrangian/centrifugal)

**Moderators:** gmalivuk, Moderators General, Prelates

- sevenperforce
**Posts:**658**Joined:**Wed Feb 04, 2015 8:01 am UTC

### Plotting effective potential (Lagrangian/centrifugal)

Hey all...

Trying to figure out the proper equations to create a contour plot capable of showing the inflection points of effective potential in a rotating reference frame. For example, when a planet is orbiting a star, the contour plot looks like this:

From what I understand, the function U(x,y) is the sum of the gravitational potential of the star, the gravitational potential of the planet, and the term for the centrifugal force arising from putting the whole thing in the rotating reference frame of the planet. Unfortunately, I can't figure out what the last term is supposed to be. I've tried a few different possibilities, but nothing seems to work.

I want to be able to look at how this contour surface changes when there is more than one planet orbiting, as in the case of two co-trojans. But to do that, I'll need to know how to recreate the original plot in the first place.

Any help would be appreciated.

Trying to figure out the proper equations to create a contour plot capable of showing the inflection points of effective potential in a rotating reference frame. For example, when a planet is orbiting a star, the contour plot looks like this:

From what I understand, the function U(x,y) is the sum of the gravitational potential of the star, the gravitational potential of the planet, and the term for the centrifugal force arising from putting the whole thing in the rotating reference frame of the planet. Unfortunately, I can't figure out what the last term is supposed to be. I've tried a few different possibilities, but nothing seems to work.

I want to be able to look at how this contour surface changes when there is more than one planet orbiting, as in the case of two co-trojans. But to do that, I'll need to know how to recreate the original plot in the first place.

Any help would be appreciated.

- gmalivuk
- GNU Terry Pratchett
**Posts:**26091**Joined:**Wed Feb 28, 2007 6:02 pm UTC**Location:**Here and There-
**Contact:**

### Re: Plotting effective potential (Lagrangian/centrifugal)

What last term? The potential of the rotating frame?sevenperforce wrote:Unfortunately, I can't figure out what the last term is supposed to be. I've tried a few different possibilities, but nothing seems to work.

You're treating it as uniform circular motion, so L

^{2}= (mr

^{2}w)

^{2}(w = omega), and the rotational contribution to potential should just be 1/2 mr

^{2}w

^{2}based on the equation on Wikipedia.

- sevenperforce
**Posts:**658**Joined:**Wed Feb 04, 2015 8:01 am UTC

### Re: Plotting effective potential (Lagrangian/centrifugal)

gmalivuk wrote:What last term? The potential of the rotating frame?sevenperforce wrote:Unfortunately, I can't figure out what the last term is supposed to be. I've tried a few different possibilities, but nothing seems to work.

You're treating it as uniform circular motion, so L^{2}= (mr^{2}w)^{2}(w = omega), and the rotational contribution to potential should just be 1/2 mr^{2}w^{2}based on the equation on Wikipedia.

It's the w = omega part of things that is throwing me. This should be time-invariant, so I'm not sure why I would need to put in values for time. I'm setting G=m=r=w=1 because I'm only interested in the shape of the surface.

What I currently have (with M = 500m and r=m=G=w=1)is:

1/(x

^{2}+ y

^{2}) - 500/(x

^{2}+ y

^{2})

^{0.5}-1/((x-1)

^{2}+ y

^{2})

^{0.5}

But that's not giving me the right contour plot.

- gmalivuk
- GNU Terry Pratchett
**Posts:**26091**Joined:**Wed Feb 28, 2007 6:02 pm UTC**Location:**Here and There-
**Contact:**

### Re: Plotting effective potential (Lagrangian/centrifugal)

It is time-invariant if omega is constant.

- sevenperforce
**Posts:**658**Joined:**Wed Feb 04, 2015 8:01 am UTC

### Re: Plotting effective potential (Lagrangian/centrifugal)

gmalivuk wrote:It is time-invariant if omega is constant.

If omega is not a constant, then what is it a function of? And, more to the point, how do I add this function as a term in a cartesian coordinate system?

- gmalivuk
- GNU Terry Pratchett
**Posts:**26091**Joined:**Wed Feb 28, 2007 6:02 pm UTC**Location:**Here and There-
**Contact:**

### Re: Plotting effective potential (Lagrangian/centrifugal)

What? Omega *is* constant, if you're assuming uniform circular motion. It's just the angular speed, which you assume to be the same at all points in space because you're interested in things that co-rotate with the planet.

However, it's a constant that's based on the mass of the star and the planet and the distance between them. You can't just assume that it's 1.

Edit: w

With your other assumptions, you can treat it as 1 if M=1 instead of m=1, and M is sufficiently large compared to m.

Edit2: The r in this equation is the distance to the planet, which you set as 1. The r in the one above for potential is sqrt(x^2+y^2).

However, it's a constant that's based on the mass of the star and the planet and the distance between them. You can't just assume that it's 1.

Edit: w

^{2}= G(M+m) / r^{3}With your other assumptions, you can treat it as 1 if M=1 instead of m=1, and M is sufficiently large compared to m.

Edit2: The r in this equation is the distance to the planet, which you set as 1. The r in the one above for potential is sqrt(x^2+y^2).

- sevenperforce
**Posts:**658**Joined:**Wed Feb 04, 2015 8:01 am UTC

### Re: Plotting effective potential (Lagrangian/centrifugal)

gmalivuk wrote:What? Omega *is* constant, if you're assuming uniform circular motion. It's just the angular speed, which you assume to be the same at all points in space because you're interested in things that co-rotate with the planet.

However, it's a constant that's based on the mass of the star and the planet and the distance between them. You can't just assume that it's 1.

Oh.

So I just need to use Newton's equation for gravitational force along with the equation for centripetal acceleration to solve for w with G = m = a = 1 and M = 500m. Which gives me w = 22.4r

^{1/2}, so the rotational contribution to potential should be 1/(2r*(r*22.4)

^{2}) = 1/2000r

^{3}.

Unfortunately when I try to plot this:

U(x,y) = 1/(2000*(x

^{2}+ y

^{2})

^{3/2}) - 500/(x

^{2}+ y

^{2})

^{0.5}- 1/((x-1)

^{2}+ y

^{2})

^{0.5}

I get this:

which isn't at all what I was looking for. Kinda stumped here.

EDIT (sniped):

gmalivuk wrote:Edit: w^{2}= G(M+m) / r^{3}

With your other assumptions, you can treat it as 1 if M=1 instead of m=1, and M is sufficiently large compared to m.

Edit2: The r in this equation is the distance to the planet, which you set as 1. The r in the one above for potential is sqrt(x^2+y^2).

But if r in the above equation is the distance to the planet, which was set at 1, then w

^{2}= G(M+m) / r

^{3}= 1*(1)/1

^{3}= 1. Or am I missing something important?

- sevenperforce
**Posts:**658**Joined:**Wed Feb 04, 2015 8:01 am UTC

### Re: Plotting effective potential (Lagrangian/centrifugal)

And if I set M = 1, then m = 1/500, right? Otherwise that second potential well goes away entirely and I don't want that.

- gmalivuk
- GNU Terry Pratchett
**Posts:**26091**Joined:**Wed Feb 28, 2007 6:02 pm UTC**Location:**Here and There-
**Contact:**

### Re: Plotting effective potential (Lagrangian/centrifugal)

No, that's exactly what I said. If you set G=M=R=1, and assume m is negligible compared to M, then w=1 (approximately) as well. (This implies an orbital speed equal to the speed of light, but you're working with pure Newton so that isn't a problem and the shape should still look qualitatively the same.)sevenperforce wrote:gmalivuk wrote:Edit: w^{2}= G(M+m) / r^{3}

With your other assumptions, you can treat it as 1 if M=1 instead of m=1, and M is sufficiently large compared to m.

Edit2: The r in this equation is the distance to the planet, which you set as 1. The r in the one above for potential is sqrt(x^2+y^2).

But if r in the above equation is the distance to the planet, which was set at 1, then w^{2}= G(M+m) / r^{3}= 1*(1)/1^{3}= 1. Or am I missing something important?

Also, when I said the rotational part of the potential is 1/2 mr

^{2}w

^{2}, I wasn't lazily leaving out any parentheses. I meant 0.5*mr

^{2}w

^{2}.

- sevenperforce
**Posts:**658**Joined:**Wed Feb 04, 2015 8:01 am UTC

### Re: Plotting effective potential (Lagrangian/centrifugal)

That makes a difference.

Apparently not enough of one, though. Plotting this...

U(x,y) = (x

...gives me this...

...which still isn't right.

If I make the centrifugal potential negative, then I finally get a better-looking contour, but there's only one maxima and it's in the wrong place:

I'm sure there's just one simple thing that I'm getting wrong but I have no idea what it might be.

Apparently not enough of one, though. Plotting this...

U(x,y) = (x

^{2}+ y^{2})/2 - (x^{2}+ y^{2})^{-0.5}- ((x-1)^{2}+ y^{2})^{-0.5}...gives me this...

...which still isn't right.

If I make the centrifugal potential negative, then I finally get a better-looking contour, but there's only one maxima and it's in the wrong place:

I'm sure there's just one simple thing that I'm getting wrong but I have no idea what it might be.

- sevenperforce
**Posts:**658**Joined:**Wed Feb 04, 2015 8:01 am UTC

### Re: Plotting effective potential (Lagrangian/centrifugal)

Changing to this...

U(x,y) = -(x

...still gives me this...

...which is the same shape.

U(x,y) = -(x

^{2}+ y^{2})/2 - (x^{2}+ y^{2})^{-0.5}- 0.002*((x-1)^{2}+ y^{2})^{-0.5}...still gives me this...

...which is the same shape.

- gmalivuk
- GNU Terry Pratchett
**Posts:**26091**Joined:**Wed Feb 28, 2007 6:02 pm UTC**Location:**Here and There-
**Contact:**

### Re: Plotting effective potential (Lagrangian/centrifugal)

Why did you switch the sign on the rotational term?

- sevenperforce
**Posts:**658**Joined:**Wed Feb 04, 2015 8:01 am UTC

### Re: Plotting effective potential (Lagrangian/centrifugal)

gmalivuk wrote:Why did you switch the sign on the rotational term?

Because otherwise I don't get any local extrema at all:

- gmalivuk
- GNU Terry Pratchett
**Posts:**26091**Joined:**Wed Feb 28, 2007 6:02 pm UTC**Location:**Here and There-
**Contact:**

### Re: Plotting effective potential (Lagrangian/centrifugal)

Ah, right, centrifugal potential should indeed decrease with increasing radius.

The main remaining problem is mostly that you're plotting with too low a resolution to see the points where the potential is level. There's also the fact that we treated M+m as equal to M, and we're pretending the center of mass of the system is exactly at the center of the central body.

The main remaining problem is mostly that you're plotting with too low a resolution to see the points where the potential is level. There's also the fact that we treated M+m as equal to M, and we're pretending the center of mass of the system is exactly at the center of the central body.

- sevenperforce
**Posts:**658**Joined:**Wed Feb 04, 2015 8:01 am UTC

### Re: Plotting effective potential (Lagrangian/centrifugal)

I'm not sure that plotting resolution is the issue; I have gone arbitrarily deep to no avail. But the center-of-mass issue and M+m issue might be the kicker. Recalculating...

For M = 50*m, the displacement from 0,0 for the primary will be d = a/(1 + 50) or at x= -0.0196, meaning that the secondary will be at x = 0.9804 opposite. Omega will now be 1.0396. This changes the equation a bit...

U(x,y) = -1.0396*(x

Success!

Looks like it was that tiny offset -- basically, I was wrapping the centrifugal potential around the primary rather than around the barycentre, so there was nowhere for the extrema to appear.

EDIT: I did all that to try and do this:

Two bodies of near-equal mass can co-orbit as each other's trojans as long as the primary is more than 26 times their combined mass. So I was wondering whether additional trojan points would arise relative to those. Apparently the answer is yes.

For M = 50*m, the displacement from 0,0 for the primary will be d = a/(1 + 50) or at x= -0.0196, meaning that the secondary will be at x = 0.9804 opposite. Omega will now be 1.0396. This changes the equation a bit...

U(x,y) = -1.0396*(x

^{2}+ y^{2})/2 - ((x+0.0196)^{2}+ y^{2})^{-0.5}- 0.02*((x-0.9804)^{2}+ y^{2})^{-0.5}Success!

Looks like it was that tiny offset -- basically, I was wrapping the centrifugal potential around the primary rather than around the barycentre, so there was nowhere for the extrema to appear.

EDIT: I did all that to try and do this:

Two bodies of near-equal mass can co-orbit as each other's trojans as long as the primary is more than 26 times their combined mass. So I was wondering whether additional trojan points would arise relative to those. Apparently the answer is yes.

### Who is online

Users browsing this forum: No registered users and 11 guests