Tiga Benda Bermasalah

Salah satu penyederhanaan dalam studi masalah tiga benda yaitu konsep masalah tiga benda terbatas.^[1] Dalam konsep ini, terdapat dua benda bermassa m_2 dan m_1 yang saling mengitari dalam orbit lingkaran, dan satu benda lain yang massanya diabaikan, yang bergerak dalam medan gravitasi m_2 dan m_1.

Merupakan penyederhanaan lain yaitu konsep sumbu korotasi. Dianggap bahwa sumbu koordinat sistem bergerak bersama m_2 dan m_1. Jadi, dua benda tersebut terletak pada sumbu-x.

Gambar 1. Deskripsi grafis sumbu korotasi dalam masalah tiga benda terbatas.

Sebagai unit satuan massa yaitu total massa sehingga massa m_2 dan m_1 yaitu \mu dan \left(1-\mu\right), dengan \mu\le1/2. Sumbu-x akan berotasi dengan kecepatan sudut \omega. Posisi m_2 dan m_1 yaitu \left(x_2,0,0\right) dan \left(x_1,0,0\right), dengan x_1 negatif. Sebagai unit satuan jarak yaitu jarak pisah m_1 dan m_2, atau \left(-x_1+x_2\right).

Benda m_3 berada pada posisi \left(x,y,z\right). Jadi, berlaku hubungan berikut.

\\ \left(x-x_1\right)^2+y^2+z^2=r_1^2\cdots\cdots\left(1a\right) \\ \left(x-x_2\right)^2+y^2+z^2=r_2^2\cdots\cdots\left(1b\right)

Jika v yaitu kelajuan m_3, maka berlaku hubungan berikut.

v^2=\left(\frac{dx}{dt}\right)^2+\left(\frac{dy}{dt}\right)^2+\left(\frac{dz}{dt}\right)^2\cdots\cdots\left(2\right)

Dengan demikian, berikut bentuk integral energi.

v^2=x^2+y^2+\frac{2\left(1-\mu\right)}{r_1}+\frac{2\mu}{r_2}-C\cdots\cdots\left(3\right)

Besaran C merupakan tetapan. Persamaan (3) dikenal dengan integral Jacobi.

Zilch Speed

Jika besaran v dalam persamaan (3) dianggap bernilai nol, didapat bentuk berikut.

x^2+y^2+\frac{2\left(1-\mu\right)}{r_1}+\frac{2\mu}{r_2}=C\cdots\cdots\left(4\right)

Nilai C tertentu mendeskripsikan batas tertentu di bidang orbit m_1 dan m_2. Dari persamaan (3), dapat disimpulkan bahwa v^2 berubah tanda jika melewati batas ini. Dengan demikian, gerakan m_3 hanya bisa berlangsung di daerah tertentu di salah satu sisi batas tersebut.

Untuk nilai C yang berbeda, bentuk daerah batas ini berbeda pula. Penggabungan daerah batas untuk bermacam-macam nilai C akan menghasilkan permukaan kecepatan relatif nol.

Gambar 2. Contoh kurva permukaan kecepatan relatif nol untuk sistem Saturnus-Dione. Plot dibuat dengan Microsoft Excel.

Fortran Frenzy

Saya telah membuat sebuah program untuk menghitung nilai C dalam persamaan (4). Nilai-nilai C yang didapat dikumpulkan dalam sebuah berkas teks. Berkas teks tersebut dapat digunakan sebagai masukan perangkat plot. Berikut kodenya. Program ini saya buat dalam bahasa Fortran, jadi maaf kalau nampak kuno sekali.

c here is the formula of distance of the third body
c from a massive body
      double precision function r(x,xn,y)
      double precision x, xn, y
      r = sqrt((x - xn)**2 + y**2)
      end

c main operations go from here
c variable r is a formula, defined in module above
c matrix c is used to store the results of
c calculation, since there are so many values
c generated, depend on x and y
c variable i and j is used to control the iteration
c variable a and b is used as indices of matrix c
c variable k is used to trace the iteration
      program zero_velocity
      double precision m1, m2, mu, d, com, x1, x2, r
      double precision x, y, step
      double precision x_ini, y_ini, c(0:30, 0:30)
      integer i, j, a, b, k

c constant parameters go here, in cgs system
c m1 is the more massive object and m2 is the lesser
c d is the separation of m1 and m2
c the orbit is assumed to be circle
      parameter(m1 = 1.989D+33, m2 = 5.977D+27,
     +          d = 1.49597892D+13)

c the matrix c is written into this file
      open(unit = 1, file = 'matrix.txt',
     +     status = 'unknown')

c here is the calculation of the position of the
c center of mass aka com
c here, com is measured from the center of m1
      com = m2*d/(m1 + m2)

c the calculations of x1, x2, and mu also go here
c as simplification, the separation of m1 and m2 is
c taken as one
c the total mass of m1 and m2 is taken as one too
c x1 is the coordinates of m1 from center of mass, so
c it is negative
c x2 is the coordinates of m2 from center of mass
c mu is the mass of m2, with simplication above
      x1 = (0D0 - com)/d
      x2 = (d - com)/d
      mu = m2/(m1 + m2)

c here go some initializations
c k indicates the numbers of iteration
c a and b determine the address of the storage of
c calculation results in matrix c
c x_ini and y_ini determine the initial x and y
c of course x_ini must be smaller than x1
c step determine the difference between two
c consecutive x or y
c here in this task, maximum of x minus minimum of x
c is three, and it applies too for y
      k = 0
      a = 0
      b = 0
      x_ini = x1 - 1D0
      y_ini = -1.5D0
      step = 1D-1

c here go the calculations of orbital energy
c c is calculated for each point on xy-plane
c each point is addressed by x and y, with zero at the
c center of mass
c the iteration scans vertical direction first (or
c y-direction)
c in fortran, first indice of two dimensional matrix
c addresses the row so, just to make it more
c intuitive, the first indice moves first, just like
c variables j and y
      do i = 0,30
      x = x_ini + i*step
        do j = 0,30
        k = k + 1
        a = j
        b = i
        y = y_ini + j*step
        c(a,b) = x**2 + y**2 + 2D0*(1 - mu)/r(x,x1,y)
     +           + 2D0*mu/r(x,x2,y)
        enddo
      enddo

c theoretically, at the center of m1 and m2, c is
c infinite
c the calculations and iteration above have take this
c into account
c however, just to make sure, here go two corrections
c for m1 and m2 location
c remember that the first indice corresponds to
c y-direction
c remember too that although y starts from minus one
c half, the first indice in c starts from zero
      y = 0
      c(15,10) = x1**2 + 2D0*(1 - mu)/r(x1,x1,y)
     +           + 2D0*mu/r(x1,x2,y)
      c(15,20) = x2**2 + 2D0*(1 - mu)/r(x2,x1,y)
     +           + 2D0*mu/r(x2,x2,y)

c here go some confirmations about the process,
c displayed on terminal
      write(*,*)
      write(*,*) 'initial x = ', x_ini
      write(*,*) 'initial y = ', y_ini
      write(*,*)
      write(*,*) '-----------------------------'
      write(*,*)
      write(*,*) 'final x = ', x_ini + 3D1*step
      write(*,*) 'final y = ', y_ini + 3D1*step
      write(*,*)
      write(*,*) '-----------------------------'
      write(*,*)
      write(*,*) 'number of x = ', i
      write(*,*) 'number of y = ', j
      write(*,*)
      write(*,*) '-----------------------------'
      write(*,*)
      write(*,*) 'C at the more massive body = '
     +           , c(15,10)
      write(*,*) 'C at the less massive body = '
     +           , c(15,20)

c here goes the writing of matrix c into the file,
c formatted below
      write(1,10) ((c(a,b), b = 0,30), a = 0,30)
   10 format(31E10.2)

c well done, we can now plot the curve from matrix.txt
c with any plot tool we have
      close(unit = 1)

      end

Rujukan

^[1] J.M.A. Danby, Fundamentals of Celestial Mechanics (Virginia: Willman-Bell, Inc., 1989), hlm. 253-259