Моделирование орбит Земли, Луны и Солнца вместе в Фортране

Question

Это ломает мой мозг уже несколько часов. Я постараюсь объяснить как можно больше. Я знаю, что этот код ужасно неблагодарен, но я новичок в Фортране, и сделать вещи более эффективными нелегко для меня.

Я пытаюсь смоделировать небесное движение, которое я указал в заголовке. Я собираюсь иметь (до перехода в совместно движущуюся систему) Солнце в начале координат, Землю на расстоянии x = 149597870d+3 метров (среднее расстояние от Солнца до Земли, которое я считаю), y = 0 Что касается Луны, она у меня на расстоянии от Луны до Земли (384400d+3 метра) плюс расстояние от Земли до Солнца, так что у нас, по сути, есть ситуация, когда все планеты расположены на y = 0 линия, где Земля является астрономической единицей на оси X, а Луна на расстоянии, равном расстоянию между Землей и Луной.

Я попытался проиллюстрировать ситуацию, а также диаграммы свободных тел на следующих двух изображениях.

Оттуда я определил 12-мерный массив, sol, который содержит положения x и y и скорости каждого тела, так что sol = (x1, y1, x2, y2, x3, y3, dx1 / dt, dy1 / dt, dx2 / dt, dy2 / dt, dx3 / dt , Dy3 / дт). Затем я инициализировал массив: у Солнца не будет начальных скоростей, у Земли будет начальная скорость y, равная его средней орбитальной скорости вокруг Солнца без начальной скорости x, а у Луны будет начальная скорость y, равная ее среднему значению. орбитальная скорость вокруг Земли и отсутствие начальной x-скорости.

 inits(1:2) = 0d0
 inits(3) = distE
 inits(4) = 0d0
 inits(5) = distE+distM
 inits(6:9) = 0d0
 inits(10) = vE
 inits(11) = 0d0
 inits(12) = vM
 sol = inits

Затем я пытаюсь привести все к движущейся системе отсчета на основе следующего уравнения:

 mf = 0
 Mtot = mSun + mEarth + mMoon
 mf(1:2) = mSun/(Mtot) * mf(1:2) + mEarth/(Mtot) * sol(3:4) + mMoon/(Mtot) * sol(5:6)
 mf(3:4) = mEarth/(Mtot) * mf(3:4) + mSun/(Mtot) * sol(1:2) + mMoon/(Mtot) * sol(5:6)
 mf(5:6) = mMoon/(Mtot) * mf(5:6) + mSun/(Mtot) * sol(1:2) + mEarth/(Mtot) * sol(3:4)
 mf(7:8) = mSun/(Mtot) * mf(7:8) + mEarth/(Mtot) * sol(9:10) +  mMoon/(Mtot) * sol(11:12)
 mf(9:10) = mEarth/(Mtot) * mf(9:10) + mSun/(Mtot) * sol(7:8) + mMoon/(Mtot) * sol(11:12)
 mf(11:12) = mMoon/(Mtot) * mf(11:12) + mSun/(Mtot) * sol(7:8) + mEarth/(Mtot) * sol(9:10)
 sol = inits - mf

Затем я использую переменный временной шаг для моих расчетов на основе уравнения:

real*8 function fun_tstepq8(arr)
use importants 
implicit none
real*8,dimension(12), intent(in) :: arr
real*8::alp,part1,part2,part3,distEtoS,distEtoM,distStoM
real*8 :: distEdottoMdot,distEdottoSdot,distSdottoMdot
alp = 1d-4
distEtoS = SQRT((sol(3)-sol(1))**2 + (sol(4)-sol(2))**2)**3
distEtoM = SQRT((sol(5)-sol(3))**2 + (sol(6)-sol(4))**2)**3
distStoM = SQRT((sol(5)-sol(1))**2 + (sol(6)-sol(2))**2)**3
distEdottoSdot = SQRT((sol(9)-sol(7))**2 + (sol(10)-sol(8))**2)
distEdottoMdot = SQRT((sol(11)-sol(9))**2 + (sol(12)-sol(10))**2)
distSdottoMdot = SQRT((sol(11)-sol(7))**2 + (sol(12)-sol(8))**2)
part1= distEtoS/distEdottoSdot
part2= distEtoM/distEdottoMdot
part3= distStoM/distSdottoMdot
fun_tstepq8 = alp * MIN(part1,part2,part3)
end function 

Собрав все это вместе, я использую метод Форвард-Эйлера и вычисляю функцию f из y0 = y0 + hf, используя подпрограмму:

subroutine q8RHS(sol)
use importants, ONLY: mSun, mEarth,mMoon, F, G
implicit none
real*8 :: distEtoS,distEtoM,distStoM
real*8,dimension(12) :: sol
integer :: i
distEtoS = SQRT((sol(3)-sol(1))**2 + (sol(4)-sol(2))**2)**3
distEtoM = SQRT((sol(5)-sol(3))**2 + (sol(6)-sol(4))**2)**3
distStoM = SQRT((sol(5)-sol(1))**2 + (sol(6)-sol(2))**2)**3
do i = 1,12    
    if (i < 7) then
        F(i) = sol(i+6)
    elseif (i < 9) then
        F(i) = G * (mEarth * (sol(i-4) - sol(i-6))/distEtoS + mMoon * (sol(i-2) - sol(i-6))/distStoM)

    elseif (i < 11) then
        F(i) =  G * mSun * (sol(i-6) - sol(i-8))/distEtoS - mMoon * G *(sol(i-4) - sol(i-6))/distEtoM
    else
        F(i) =  -G * mSun * (sol(i-6) - sol(i-10))/distStoM -G*mEarth * (sol(i-6) - sol(i-8))/distEtoM

    endif
enddo
end subroutine  

В терминах получения подпрограммы q8RHS я начал с: 1043 *

И вот моя работа по выводу механики для одной из планет. Заметьте, что здесь я использую метод Эйлера и извлекаю подходящие функции для формы f (x, t) = dx / dt, необходимой для метода Эйлера. Заметьте, что я не мог использовать знаменатель, который я получил, как в какой-то момент x2 = x1, в моем моделировании, если я использую это, поэтому проблема, о которой я говорил, с временным шагом, становящимся крошечным, когда моделирование останавливается.

В целом моя программа выглядит так:

module importants
implicit none
integer,parameter::Ndim = 3
integer:: N,N2
real*8, parameter :: G = 6.67408e-11
integer,parameter :: Nbodies = 2
integer::specific_Nbodies = 2*Nbodies
real*8,dimension(12) :: inits,sol,F
real*8 :: d_j = 778.547200d+9
real*8 :: v_j = 13.1d+3
integer :: rank = Nbodies * Ndim * 2
real*8 :: mSun = 1.989d+30    
real*8, dimension(12) :: mf
real*8 :: mJup = 1.89819d+27
real*8, parameter :: pi= DACOS(-1.d0)
real*8 :: a,P,FF
real*8 :: pJ = 374080032d0
real*8, dimension(2) :: QEcompare,QLcompare,QPcompare
real*8 :: mEarth = 5.9722d+24
real*8 :: mMoon = 7.34767309d+22
real*8 :: distE = 149597870d+3 ! Dist from sun to earth
real*8 :: distM = 384400d+3 ! Dist from earth to moon
real*8 :: distMS = 152d+9
real*8 :: vE = 29.78d+3
real*8 :: vM = 1.022d+3
end module

program exercise3
use importants
implicit none
print*,'two'
!call solve()
!call q3()
!call q4()
!call q5()
!call q5circ()
!call q6ecc()
!call q6pt2()
!call q7()
call q8()
end program

subroutine q8()
 use importants
 implicit none
 !real*8,external :: Epot,Ekin, L
 real*8,external :: fun_tstepq8
 real*8 :: tstep,time,tracker,Mtot,t
 !real*8,dimension(2) :: Etot,L1,L2
 !real*8, dimension(3,2) :: r2  
 integer :: j,i
 print*, 'Starting Q8. -------------------------------------------------------'
 inits(1:2) = 0d0
 inits(3) = distE
 inits(4) = 0d0
 inits(5) = distE+distM
 inits(6:9) = 0d0
 inits(10) = vE
 inits(11) = 0d0
 inits(12) = vM
 sol = inits
 mf = 0
 Mtot = mSun + mEarth + mMoon
 mf(1:2) = mSun/(Mtot) * mf(1:2) + mEarth/(Mtot) * sol(3:4) + mMoon/(Mtot) * sol(5:6)
 mf(3:4) = mEarth/(Mtot) * mf(3:4) + mSun/(Mtot) * sol(1:2) + mMoon/(Mtot) * sol(5:6)
 mf(5:6) = mMoon/(Mtot) * mf(5:6) + mSun/(Mtot) * sol(1:2) + mEarth/(Mtot) * sol(3:4)
 mf(7:8) = mSun/(Mtot) * mf(7:8) + mEarth/(Mtot) * sol(9:10) +  mMoon/(Mtot) * sol(11:12)
 mf(9:10) = mEarth/(Mtot) * mf(9:10) + mSun/(Mtot) * sol(7:8) + mMoon/(Mtot) * sol(11:12)
 mf(11:12) = mMoon/(Mtot) * mf(11:12) + mSun/(Mtot) * sol(7:8) + mEarth/(Mtot) * sol(9:10)
 sol = inits - mf
 print*, '---------------------------------------------------------------------'
 print*, 'Beginning fixed t timestep. '
 t = 0d0
 P = 3.154d+7
 time=0
 tracker=0
 print*, 'P understood as :', P
 open(70,file='q8.dat')
 time=0
 tracker=0
 write(70,*) sol(1),sol(2),sol(3),sol(4),sol(5),sol(6)
 do while(time.le.P)
     call q8RHS(sol)
     tstep=fun_tstepq8(sol)
     time=time+tstep
     tracker=tracker+1
     print*, tstep
     do i=1,12
       sol(i) = sol(i) + tstep * F(i)
     enddo
 write(70,*) sol(1),sol(2),sol(3),sol(4),sol(5),sol(6)

 enddo

end subroutine 

Я знаю, что есть множество переменных, которые здесь не используются, но использовались для других подпрограмм, поскольку это немного неуклюжий зверь из программы, которую я создал. Я просто знаю, что предыдущие подпрограммы и все связанные с ними переменные работали нормально, и у меня нет с ними проблем, и они не мешают коду, который я имею в виду.

В любом случае, когда я запускаю программу, временной шаг сразу становится огромным, а затем программа завершается за одну или две итерации, так как знаменатели в моей функции временного шага очень быстро уменьшают свои знаменатели. Перед тем, как показать код, который у меня есть, я сначала заметил, что Земля и Луна не взаимодействуют, и Луна просто движется к Солнцу, в то время как Земля вращается вокруг Солнца нормально. Выглядело это так (обратите внимание, что это с другой подпрограммой RHS, не той, что я перечислил, а только для контекста):

Я пытался затем сначала посмотреть на Землю и саму Луну с фиксированным временным шагом. В конце концов я заставил их «увидеть друг друга» (взаимодействовать с их гравитационными полями так, чтобы это имело смысл), но как только Луна приблизилась к Земле, она отмахнулась от нее (думаю, как только она прошла мимо Земля так, что она была чуть ближе к Солнцу, чем Земля) и просто продолжала стрелять. Я пытался использовать различный временной шаг, но тогда у него был такой маленький временной шаг, все замедлялось до половины, когда Луна подходила достаточно близко к Земле. Для этого использовалось уравнение временного шага, в котором учитывались только координаты x Луны и Земли. Затем я попытался отрегулировать временной шаг так, чтобы он выглядел так, как я его перечислил, глядя на уравнение, которое я дольше изобразил, что вызвало огромные временные шаги, которые я обсуждал. Затем я сдался и возвратил присутствие Солнца, которое является кодом, который я перечислил здесь. Я абсолютно в тупике. Вот как выглядят траектории с имеющимся у меня кодом:

Что не так с моим кодом? Что я не симулирую орбиту должным образом? Что я делаю не так?

ОБНОВЛЕНИЕ

По совету одного из комментаторов я сделал alp меньше. Вот результат:

и более увеличенный в версии.

Земля прямо отталкивается, а луна летит к Солнцу.

ОБНОВЛЕНИЕ 2

Я считаю, что исправил ошибку, так как считаю, что бит q8RHS, относящийся к функции Эйлера Земли, перевернул свои знаки.

        elseif (i < 11) then
        F(i) =  -G * mSun * (sol(i-6) - sol(i-8))/distEtoS + mMoon * G *(sol(i-4) - sol(i-6))/distEtoM

Поскольку гравитационное притяжение от Солнца должно быть отрицательным, а с Луны - положительным, учитывая то, как я это нарисовал. Вот самое новое изображение

ОБНОВЛЕНИЕ 3

Мне сообщили, что моя начальная орбитальная скорость поскольку Луна была неправа - это должна быть сумма ее орбитальной скорости вокруг Луны плюс орбитальная скорость Земли вокруг Солнца. Мои движения теперь:

Это было гораздо более обнадеживающим, но я не уверен, что Земля и Луна ведут себя правильно с С уважением друг к другу. При увеличении, кажется, что луна колеблется влево и вправо от траектории Солнца, поэтому, если бы ось z была видимой, это могло бы иметь больше смысла, но также не могло бы.

Поведение выше пути Луны, движущегося влево и вправо от Земли, постоянно повторяется, но я не уверен, что это все еще обязательно физически разумным.