Catatan – 2 benda bermasalah

Maksud tulisan ini adalah ingin sedikit membahas masalah dua benda (two-body problem) atau lebih tepatnya gravitational two-body problem. Kita tahu bahwa solusi dari gerak 2 benda dalam pengaruh gravitasi adalah gerak Keplerian, dimana dua benda akan saling mengorbit satu sama lain, dan bentuk orbitnya hanya ditentukan oleh energi (E) dan momentum sudut (L).

Mari kita tinjau kasus benda massless yang dipengaruhi oleh satu benda masif. Misalkan benda yang massless memiliki kecepatan awal v yang kurang dari kecepatan lepas, atau dengan kata lain E negatif, maka orbit yang kita peroleh adalah orbit elips atau lingkaran (e < 1). Sebagai contoh, saya punya benda masif di (0,0), lalu benda kedua yang massless saya lempar dari posisi (5, 0) dengan kecepatan awal yang sama, namun arahnya berbeda (momentum sudut — L berbeda), maka orbit yang dihasilkan pasti elips, dengan eksentrisitas berbeda tergantung arah tadi. Energi yang sama ditunjukkan dengan semimajor (a) yang sama dalam gambar berikut.
ellips

—– ** —–

Lalu pertanyaan muncul, bagaimana dengan bentuk orbit parabola dan hiperbola ditinjau dengan cara yang sama? Jujur susah membayangkan bentuknya, sebagai contoh jika kita melempar batu dari permukaan Bumi dengan kecepatan lepas di permukaan Bumi, namun saya lempar ke arah berbeda, bagaimana cara menggambar orbit parabola dari sudut awal lemparan saya tadi, harusnya berbeda namun tetap parabola bukan?

maka dari itu saya lakukan kerjaan ini.., ternyata hasilnya sebagai berikut:
parabola

untuk hiperbola:
hyperbola

Catatan Penting! kita harus berhati-hati dalam memplot posisi, software plot biasa memplot dengan membuat sumbu-x dan y dengan skala berbeda.. namun untuk plot “posisi” sebaiknya selalu gunakan skala yang sama.. karena akan berpengaruh pada interpretasi dan gambaran di kepala kita. Sebagai contoh, parabola di atas kalau di plot dengan tanpa mengatur skala yang sama, hasilnya akan aneh seperti dibawah ini (dalam kasus ini kita akan melihat massa pusat tidak berada di fokus orbit!).
parabola_aneh

Uniform random points in disk, annulus (ring), cylinder, and sphere

[Postingan ini dibuat karena ada yang nanya... :D haha]

Suppose we have u, v, and w as uniform random from 0 to 1. Here’s how we make uniform random points in several cases in corresponds coordinate system.

Disk
in Polar coordinate:
\theta = 2 \pi u
r = \sqrt{R^2 v} = R \sqrt{v}

x = r \cos \theta
y = r \sin \theta

// Copyleft (c) 2014. Ridlo W. Wibowo
#include<iostream>
#include<stdlib.h>
#include<math.h>
#include<stdio.h>
#include<time.h>
#include<fstream>
#define _USE_MATH_DEFINES
using namespace std;
 
 
double unirand(){return (double) rand()/(double) RAND_MAX;}
 
int main(){
    double rho, theta;
    double r = 1.0;
    double x, y;
    int n=10000;
 
    srand(time(NULL));
    ofstream out("circle.txt");
    for (int i=0;i<n;i++){
        rho = r*sqrt(unirand()); // uniform in r
        // rho = r*unirand(); // non-uniform
        theta = 2.*M_PI*unirand();
 
        x = rho*cos(theta);
        y = rho*sin(theta);
         
        out << x << " " << y << endl;
    }
    out.close();
     
    return 0;
}
Wrong

Wrong

Right

Right


Annulus
in Polar coordinate:
\theta = 2 \pi u
r = \sqrt{(R_{end}^{2} - R_{start}^{2}) v + R_{start}^{2}}

x = r \cos \theta
y = r \sin \theta

annulus

annulus2


Cylinder
in Cylindrical coordinate:
\theta = 2 \pi u
\rho = R \sqrt{v}

x = r \cos \theta
y = r \sin \theta
z = (z_{2} - z_{1}) w + z_{1}

cylinder01

cylinder02


Direction in 3D (spherical surface)
in Spherical coordinate:
\phi = 2 \pi u
\theta = \arccos(2 v - 1)
with a contant r or r = 1

x = r \sin \theta \cos \phi
y = r \sin \theta \sin \phi
z = r \cos \theta

not uniform distribution

not a uniform distribution

uniform distribution

uniform distribution


Sphere
in Spherical coordinate:
\phi = 2 \pi u
\theta = \arccos(2 v - 1)
r = \sqrt[3]{R^{3} w} = R \sqrt[3]{w}

x = r \sin \theta \cos \phi
y = r \sin \theta \sin \phi
z = r \cos \theta

// Copyleft (c) 2014. Ridlo W. Wibowo
#include<iostream>
#include<stdlib.h>
#include<math.h>
#include<stdio.h>
#include<time.h>
#include<fstream>
#define _USE_MATH_DEFINES
using namespace std;
 
 
double unirand(){return (double) rand()/(double) RAND_MAX;}
 
int main(){
    double r, phi, theta;
    double x, y, z;
    int n=10000;
 
    srand(time(NULL));
    ofstream out("bola.txt");
    for (int i=0;i<n;i++){
        phi = 2.*M_PI*unirand();
        theta = acos(1. - 2.*unirand());
        r = pow(unirand(), 0.33333333333);

        x = r*sin(theta)*cos(phi);
        y = r*sin(theta)*sin(phi);
        z = r*cos(theta);
         
        out << x << " " << y << " " << z << " " << endl;
    }
    out.close();
     
    return 0;
}

spheris

Gnuplot script

# gnuplot input file
set term wxt size 800,600
set parametric
set angle degree
set urange [0:360]
set vrange [-90:90]
set isosample 12,11
set ticslevel 0
#set xtics -0.06,0.02,0.06
#set ytics -0.06,0.02,0.06
set size 0.7,1.0
set nokey
a=1.0
fx(u,v)=cos(u)*cos(v)
fy(u,v)=sin(u)*cos(v)
fz(v)=sin(v)
splot a*fx(u,v),a*fy(u,v),a*fz(v), "bola.txt" w d lc 3

Gadget2 and HDF5

Kemarin saya iseng menginstal Gadget 2A code for cosmological simulations of structure formation.  Cara instalasi dapat dilihat pada website http://astrobites.org/2011/04/02/installing-and-running-gadget-2/.

Setelah berhasil menginstal, lalu dicoba examplenya.. (btw saya instal di server Perseus, Prodi Astronomi, karena dependensinya sudah terinstal semua, jadi mudah, haha). Example yang sudah saya coba adalah tumbukan galaksi. Boleh dibilang problem SPH (Smoothed Particle Hydrodynamics) klasik, dengan dua jenis partikel, perhitungan gravitasinya menggunakan metode Tree. Untuk tujuan pembelajaran saya sekaligus keluarkan outputnya dalam format hdf5, lalu cari tahu cara membaca format file tersebut. HDF5 adalah format file yang katanya paling efektif untuk data besar dan kompleks; banyak bidang sudah mengimplemetasikan format file ini antara lain di meteorologi dan astronomi. :D

Hasil output hdf5 coba saya baca dengan script python sederhana.. ingin coba multiprocess di python namun tidak sempat, ya sudahlah kapan2 saja. Setelah mengulik cara pakainya akhirnya berhasil, haha… berikut script sederhana untuk membaca dan plot tersebut:

#!/usr/bin/python -tt

import matplotlib
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import math
import h5py
import numpy as numpy
from pylab import *

DR = 180./math.pi
matplotlib.rcParams.update({'font.size': 32})
TWOPI = 2.*math.pi

def plotter(ifile):
    f = h5py.File(ifile+".hdf5", 'r')
    halo = f['/PartType1']
    disk = f['/PartType2']

    figure(figsize=(45,15))
    ax1 = subplot(131, aspect='equal', xlim=[-300, 300], ylim=[-300, 300])
    ax1.scatter(halo['Coordinates'][:,0], halo['Coordinates'][:,1], color='red', s=0.01)
    ax1.scatter(disk['Coordinates'][:,0], disk['Coordinates'][:,1], color='blue', s=0.3)
    ax1.set_xlabel(r"$x$")
    ax1.set_ylabel(r"$y$")

    ax2 = subplot(132, aspect='equal', xlim=[-300, 300], ylim=[-300, 300])
    ax2.scatter(halo['Coordinates'][:,0], halo['Coordinates'][:,2], color='red', s=0.01)
    ax2.scatter(disk['Coordinates'][:,0], disk['Coordinates'][:,2], color='blue', s=0.3)
    ax2.set_xlabel(r"$x$")
    ax2.set_ylabel(r"$z$")

    ax3 = subplot(133, aspect='equal', xlim=[-300, 300], ylim=[-300, 300])
    ax3.scatter(halo['Coordinates'][:,1], halo['Coordinates'][:,2], color='red', s=0.01)
    ax3.scatter(disk['Coordinates'][:,1], disk['Coordinates'][:,2], color='blue', s=0.3)
    ax3.set_xlabel(r"$y$")
    ax3.set_ylabel(r"$z$")

    tight_layout()
    savefig(ifile+'.png')
    f.close() # don't forget to close it

if __name__ == "__main__":
    n_snap = 151
    for i in range(0, n_snap, 1):
        if i < 10:
            j = "00"+str(i)
        elif i>=10 and i<100:
            j = "0"+str(i)
        else:
            j = str(i)
        ifile = "snapshot_"+j # without .hdf5
        plotter(ifile)

Setelah digabungkan 150 “screenshot” berupa image menjadi movie, maka hasilnya sebagai berikut

ffmpeg -f image2 -r 15 -i 'snapshot_%03d.png' col.mp4

warna biru dan merah mewakili partikel disk dan halo. Plot 150 gambar tersebut membutuhkan waktu 18 menit dengan program di atas. Hasilnya cukup menggambarkan yang di dapat di alam seperti ini (galaksi antennae.. mirip bentuk cintah.. <3 :v ):
antena_collision
Credit & Copyright: Star Shadows Remote Observatory and PROMPT/CTIO
(Jack Harvey, Steve Mazlin, Rick Gilbert, and Daniel Verschatse)

berikutnya ingin saya coba yang simulasi kosmologi.. :D

The Precession of Mercury’s Perihelion (simulation)

One of the classical test of general relativity is the perihelion precession of Mercury (also gravitational lensing & gravitational redshift).

I1_1
Image source: http://www.relativity.li/en/epstein2/read/i0_en/i1_en/

Under Newtonian physics, a two-body system consisting of a lone object orbiting a spherical mass would trace out an ellipse with the spherical mass at a focus. The point of closest approach, called the periapsis (or, as the central body in our Solar System is the sun, perihelion), is fixed. A number of effects in our solar system cause the perihelia of planets to precess (rotate) around the sun. The principal cause is the presence of other planets which perturb each other’s orbit. Another (much less significant) effect is solar oblateness.

In 1859, the French mathematician and astronomer Urbain Le Verrier reported that the slow precession of Mercury’s orbit around the Sun could not be completely explained by “Newtonian mechanics” and “perturbations by the known planets“.

Observation: The perihelion precession of Mercury is 574.10±0.65 arcseconds per century relative to the inertial ICFR (International Celestial Reference Frame: barycenter of solar system).

Classic calculation: Newtonian mechanics, taking into account all the effects from the other planets, predicts a precession 531.63±0.69 arcseconds per century.

In the early 20th century, Albert Einstein’s General Theory of Relativity provided the explanation for the observed precession. The effect is very small: the Mercurian relativistic perihelion advance excess is just 42.98 arcseconds per century. Similar, but much smaller, effects operate for other planets: 8.62 arcseconds per century for Venus, 3.84 for Earth, 1.35 for Mars, and 10.05 for 1566 Icarus.

source: wikipedia hahahaha

Sooo…. Lets try this out! :) with numerical simulation..

First, of course the theoretical background, you can find the derivation of the General relativity correction in literature. :D ahaha
Ok, we can calculate that in gravitational system consist of many particles, the Post-Newtonian equation of motion is [ujug-ujug.. mak bedunduk]:
PN
then, become [klothaa..ak]:
(the derivation is a lot more complicated, LOL)
Post-newtonian
Source: Benitez and Gallardo, Celest Mech Dyn Astr 2008
I really want to calculate the gravity using above equation, but if we just want to know Mercury’s precession, that equation is overkill.. [halah]. (maybe next time, e.g. for binary Neutron star, globular cluster, or S-star around supermassive black hole :) )

Assume that only the central body give effect of general relativity, the general relativity correction given by this simple equation:
PN-correction

WOKE, lets do our simulation.. using Rebound integrator. Oh wait, there is no GR-correction in this integrator yet, so we can add it,
simply give “additional_force” like this:

void relativity_correction_force(){
    for (int i=1; i<N; i++){ // for all planet 
        struct particle* com = &(particles[i]);
    
        // relative to the Sun
        struct particle sun = particles[0];
        double dx = com->x - sun.x;
        double dy = com->y - sun.y;
        double dz = com->z - sun.z;
    
        double dvx = com->vx - sun.vx;
        double dvy = com->vy - sun.vy;
        double dvz = com->vz - sun.vz;

        double r = sqrt(dx*dx + dy*dy + dz*dz);
        double _r = 1./r; // reciprocal
    
        // Benitez and Gallardo 2008
        // alpha = k^2/(r^3*c^2); beta  = 4k^2/r - v.v; gamma = 4*(v.r)
        // delta_a = alpha*(beta*r + gamma*v)
        double alpha = G*_r*_r*_r*_c2;
        double beta  = 4.*G*_r - (dvx*dvx + dvy*dvy + dvz*dvz);
        double gamma = 4.*(dvx*dx + dvy*dy + dvz*dz);
    
        // add acceleration to comet
        com->ax += alpha * (beta*dx + gamma*dvx);
        com->ay += alpha * (beta*dy + gamma*dvy);
        com->az += alpha * (beta*dz + gamma*dvz);
    }
}

or similar with this code. :)

using this program and some initial data of planet from HORIZON JPL NASA we can make simulation of “Solar System” for example in 10000 years. We want to know how much the different of the Mercury’s orbit with and without GR correction, hopefully I can get the value given above. ahahaha

Result
Calculation:
Package Rebound
Integrator: IAS15 (similar with RADAU15)
time step: 0.1 days
integration time: 10000 years
including all planets, moon, and 3 massive asteroid.. LOL (why these asteroids here? because I’m too lazy to delete it)

computing time (OpenMP 8 threads)
with GR-correction 1011.007421 s
w/o (Classic) 948.855602 s

First of all lets remind you (if u already know) about Keplerian coordinate, hehe, its easier to look how the orbit change using this coordinate. :D
400px-Orbit1.svg

Mercury’s orbit change after 10000 years:
mercury_a

mercury_e

mercury_i

mercury_omega

mercury_Omega

Different between integration w/o and with GR-correction is clearly visible in argument of perihelion parameter.
mercury_omega2
assume that the other parameters are exactly the same (not change due to GR-correction) we can simply calculate the precession of Mercury’s perihelion using argument of perihelion. After 10000 years the different is 1.22397536º or 44.063 arcseconds/century. Uh.. not exactly as calculated before (by someone).. hahaha. I think 1 arcseconds is tolerable for this simulation :P
This different perhaps due to:
– that last calculation only using argument of perihelion
– need a smaller “integrator_epsilon” (error)
– the precession is not linear! (look in paper above)
– didn’t include correction from Oblateness of the Sun (quadrupole moment of gravity)

:P