On Christmas morning, I was watching the launch of the JWST. As they talked about it, they said it would orbit at the L2 Lagrange point to stay at a position opposite the Sun from the Earth. I looked up the info on it, but I wasn’t sure why it would stay at that point. I got my 3 body orbit program and figured I’d see if it could do the same thing. So I set the initial values to position the Sun, Earth, and the JWST. It took quite awhile of altering the values of the JWST, but I finally got it to orbit for about 2.5 Earth orbits before it became unstable and went off on its own. After watching it, I could see what was really happening. The JWST was at an orbit position so that it took about one year to orbit the Earth. Since the Earth took one year to orbit the Sun, the JWST would remain in the same position relative to the Sun Earth line, keeping it opposite the Sun.

So here’s the program. It only works for about 2.5 Earth orbits before it goes unstable because of the accuracy of the calculations and scale. If you uncomment the code at the end, that shows the same thing but with just a simple sin cos calculation of points revolving around in a circle.

```
viewer.mode=FULLSCREEN
function setup()
xx=0
yy=0
-- object 1
gs1 = 1000 -- gravity of object 1
xs1 = WIDTH/2 -- start x position
ys1 = HEIGHT/2 -- start y position
dxs1 = 0 -- initial change in x position
dys1 = 0 -- initial change in y position
-- object 2
gs2 = 10 -- gravity of object 2
xs2 = WIDTH/2 -- start x position
ys2 = HEIGHT/2+200 -- start y position
dxs2 = -2.236 -- initial change in x position
dys2 = 0 -- initial change in y position
-- object 3
gs3 = 0 -- gravity of object 3
xs3 = WIDTH/2 -- start x position
ys3 = HEIGHT/2+230 -- start y position
dxs3 = -2.6475032398507 -- initial change in x position
dys3 = 0 -- initial change in y position
end
function draw()
background(0)
-- gravity object 2 to object 1
d = gs2 / math.pow(( (xs1 - xs2) * (xs1 - xs2) + (ys1 - ys2) * (ys1 - ys2) ), 1.5 )
dxs1 = dxs1 - (( xs1 - xs2 ) * d)
dys1 = dys1 - (( ys1 - ys2 ) * d)
-- gravity object 1 to object 2
d = gs1 / math.pow(( (xs2 - xs1) * (xs2 - xs1) + (ys2 - ys1) * (ys2 - ys1) ), 1.5 )
dxs2 = dxs2 - (( xs2 - xs1 ) * d)
dys2 = dys2 - (( ys2 - ys1 ) * d)
-- gravity object 2 to object 3
d = gs2 / math.pow(( (xs3 - xs2) * (xs3 - xs2) + (ys3 - ys2) * (ys3 - ys2) ), 1.5 )
dxs3 = dxs3 - (( xs3 - xs2 ) * d)
dys3 = dys3 - (( ys3 - ys2 ) * d)
-- gravity object 3 to object 2
d = gs3 / math.pow(( (xs2 - xs3) * (xs2 - xs3) + (ys2 - ys3) * (ys2 - ys3) ), 1.5 )
dxs2 = dxs2 - (( xs2 - xs3 ) * d)
dys2 = dys2 - (( ys2 - ys3 ) * d)
-- gravity object 1 to object 3
d = gs1 / math.pow(( (xs3 - xs1) * (xs3 - xs1) + (ys3 - ys1) * (ys3 - ys1) ), 1.5 )
dxs3 = dxs3 - (( xs3 - xs1 ) * d)
dys3 = dys3 - (( ys3 - ys1 ) * d)
-- gravity object 3 to object 1
d = gs3 / math.pow(( (xs1 - xs3) * (xs1 - xs3) + (ys1 - ys3) * (ys1 - ys3) ), 1.5 )
dxs1 = dxs1 - (( xs1 - xs3 ) * d)
dys1 = dys1 - (( ys1 - ys3 ) * d)
xs1 = xs1 + dxs1
ys1 = ys1 + dys1
xs2 = xs2 + dxs2
ys2 = ys2 + dys2
xs3 = xs3 + dxs3
ys3 = ys3 + dys3
fill(255, 255, 0)
ellipse(xs1-xx,ys1-yy,25)
fill(0,255,0)
ellipse(xs2-xx,ys2-yy,10)
fill(255)
ellipse(xs3-xx,ys3-yy,5)
end
--[[
function setup()
ang=0
end
function draw()
background(0)
fill(255, 255, 0)
ellipse(WIDTH/2,HEIGHT/2,25)
fill(0,255,0)
x=math.cos(math.rad(ang))*200
y=math.sin(math.rad(ang))*200
ellipse(WIDTH/2+x,HEIGHT/2+y,10)
fill(255)
x=math.cos(math.rad(ang))*230
y=math.sin(math.rad(ang))*230
ellipse(WIDTH/2+x,HEIGHT/2+y,5)
ang=ang+1
end
--]]
```