Impor semua paket yang dibutuhkan

from skyfield import api
from skyfield import almanac
from skyfield.nutationlib import iau2000b
from skyfield.units import Angle
from datetime import timedelta

import astropy.units as u

from datetime import datetime
from pytz import timezone

Buat timescale dan unduh ephemeris. Kita akan menggunakan ephemeris dari JPL NASA versi de421.bsp.

timescale = api.load.timescale()
ephemeris = api.load('de421.bsp')

Tentukan koodinat di permukaan Bumi di mana prediksi akan dilakukan, sebagai contoh kita gunakan Observatorium Bosscha. Koodinat yang digunakan dalam derajat.

# Observatorium Bosscha 107:37:00 BT; -6:49:50 LS
longitude = 107 + 37/60
latitude = -(6 + 49/60 + 50/3600)

Perhitungan posisi Bulan dan Matahari bergantung pada posisi pengamat di Bumi. Oleh karena itu kita perlu menentukan lokasi topografi Skyfield untuk koordinat di atas.

bosschaTopo = api.Topos(longitude=longitude, latitude=latitude)

Waktu Matahari dan Bulan Terbenam

Misalkan ingin menghitung waktu terbenam antara tanggal 5 April dan 6 April 2019.

t1 = timescale.utc(2019, 4, 5, 0)
t2 = timescale.utc(2019, 4, 6, 0)

Paket Skyfield telah menyediakan fungsi untuk menghitung waktu terbenam Matahari. Definisi matahari terbenam yang digunakan paket Skyfield sesuai dengan definisi Matahari terbenam dari USNO, yaitu ketika pusat piringan Matahari berada pada jarak zenit $90^o 50′$.

sunriset, sunBol = almanac.find_discrete(t1, t2, almanac.sunrise_sunset(ephemeris, bosschaTopo))

Hasilnya adalah sebuah array waktu, dan array boolean yang akan bernilai True jika matahari terbit dan bernilai False jika matahari terbenam.

Untuk menghitung waktu terbenam Bulan dengan definisi di atas kita dapat membuat fungsi sebagai berikut.

def moonrise_moonset(ephemeris, topos):
    moon = ephemeris['moon']
    topos_at = (ephemeris['earth'] + topos).at

    def is_moon_up_at(t):
        t._nutation_angles = iau2000b(t.tt)
        return topos_at(t).observe(moon).apparent().altaz()[0].degrees > -50/60

    is_moon_up_at.rough_period = 0.5
    return is_moon_up_at
moonriset, moonBol = almanac.find_discrete(t1, t2, moonrise_moonset(ephemeris, bosschaTopo))

Keluarannya akan mirip dengan format waktu terbenam Matahari.

Waktu terbenam dalam Julian Date

sunriset[~sunBol]
<Time tt=[2458578.95413447]>
moonriset[~moonBol]
<Time tt=[2458578.96095373]>

Waktu terbenam dalam standar UTC

sunriset[~sunBol].utc_iso()
['2019-04-05T10:52:48Z']
moonriset[~moonBol].utc_iso()
['2019-04-05T11:02:37Z']

Posisi Matahari dan Bulan saat Matahari Terbenam

Definisikan Bumi, Bulan, dan Matahari dari ephemeris di atas.

earth = ephemeris['earth']
sun = ephemeris['sun']
moon = ephemeris['moon']

Lokasi pengamat adalah posisi pusat Bumi ditambahkan posisi lokal topografi.

bosscha = earth + bosschaTopo

Koordinat RA dan Dec saat Matahari terbenam

sunradec = bosscha.at(sunriset[~sunBol]).observe(sun).apparent()
moonradec = bosscha.at(sunriset[~sunBol]).observe(moon).apparent()

RA dan Dec epoch saat itu.

print(sunradec.radec(epoch='date'))
(<Angle 00h 56m 36.75s>, <Angle +06deg 03' 05.6">, <Distance [1.00036989] au>)
print(moonradec.radec(epoch='date'))
(<Angle 01h 04m 12.04s>, <Angle +01deg 56' 37.8">, <Distance [0.00265776] au>)

Koordinat dalam Altitud dan Azimut

print(sunradec.altaz())
(<Angle -00deg 49' 59.9">, <Angle 275deg 59' 42.8">, <Distance [1.00036989] au>)
print(moonradec.altaz())
(<Angle 01deg 32' 08.7">, <Angle 272deg 08' 33.2">, <Distance [0.00265776] au>)

Elongasi

sunradec.separation_from(moonradec)
<Angle 04deg 31' 21.0">

Fraksi Iluminasi Bulan saat Matahari terbenam

almanac.fraction_illuminated(ephemeris, 'moon', sunriset[~sunBol])[0]*100
0.19602608000200417

Waktu konjungsi

Waktu konjungsi dapat ditentukan dengan mencari waktu saat fase bulan bernilai 0.

t0 = timescale.utc(2019, 4, 4, 0)
t1 = timescale.utc(2019, 4, 6, 0)
t, y = almanac.find_discrete(t0, t1, almanac.moon_phases(ephemeris))
t[y==0].utc_iso()
['2019-04-05T08:50:29Z']

Fungsi

Untuk memudahkan menjalankan seluruh prosedur di atas, kita dapat membuat fungsi yang mencakup seluruh prosedur sekaligus membuat plot posisi Bulan dan Matahari.

%matplotlib inline
from matplotlib import pyplot as plt
plt.style.use('ggplot')
longitude = 107+37/60
latitude = -(6+49/60+50/3600)
bosschaTopo = api.Topos(longitude=longitude, latitude=latitude)

def ucupkiyah(date, plot=False, day=0):
    
    wib = timezone('Asia/Jakarta')

    date = wib.localize(date)
    t0 = timescale.utc(date + timedelta(days=day))
    
    t1 = timescale.utc(t0.astimezone(wib) - timedelta(days=28))
    t2 = timescale.utc(t0.astimezone(wib) + timedelta(days=28))
    t, y = almanac.find_discrete(t1, t2, almanac.moon_phases(ephemeris))
    
    # t0 = t[y==0]
    t1 = timescale.utc(t0.astimezone(wib) - timedelta(days=0))
    t2 = timescale.utc(t0.astimezone(wib) + timedelta(days=1))
    
    sunriset, sunBol = almanac.find_discrete(t1, t2, almanac.sunrise_sunset(ephemeris, bosschaTopo))
    moonriset, moonBol = almanac.find_discrete(t1, t2, moonrise_moonset(ephemeris, bosschaTopo))

    sunset = sunriset[~sunBol].astimezone(wib)[0]
    moonset = moonriset[~moonBol].astimezone(wib)[0]
    
    moonAge = sunriset[~sunBol] - t[y==0]
    moonAge = moonAge[0]*24*3600
    
    lagTime = (moonriset[~moonBol][0]-sunriset[~sunBol][0])*24*3600
    
    earth = ephemeris['earth']
    sun = ephemeris['sun']
    moon = ephemeris['moon']

    bosscha = earth + bosschaTopo
    
    sunradec = bosscha.at(sunriset[~sunBol]).observe(sun).apparent()
    moonradec = bosscha.at(sunriset[~sunBol]).observe(moon).apparent()
    
    sunra, sundec, _ = sunradec.radec(epoch='date')
    moonra, moondec, _ = moonradec.radec(epoch='date')
    sunal, sunaz, _ = sunradec.altaz()
    moonal, moonaz, _ = moonradec.altaz()
    
    relal = moonal.degrees[0] - sunal.degrees[0]
    relaz = moonaz.degrees[0] - sunaz.degrees[0]
        
    elongation = sunradec.separation_from(moonradec)
    illumination = almanac.fraction_illuminated(ephemeris, 'moon', sunriset[~sunBol])[0]*100

    print('By the Name of Allah')
    print('Ucupkiyah Calculator by Muhammad Yusuf')
    print('')
    
    print('G. Conjunction Time\t : {}'.format(t[y==0].astimezone(wib)[0]))
    print('Sunset Time \t \t : {}'.format(sunset))
    print('Moonset Time \t \t : {}'.format(moonset))
    print('')
    
    print('G. Moon Age\t \t : {:02.0f}h {:02.0f}m {:02.2f}s'.format(moonAge//3600, moonAge%3600//60, moonAge%60))
    print('Lag Time \t \t : {:02.0f}h {:02.0f}m {:02.2f}s'.format(lagTime//3600, lagTime%3600//60, lagTime%60)) 
    print('')
    
    print('Sun RA \t \t \t : {}'.format(sunra))
    print('Moon RA \t \t : {}'.format(moonra))
    print('')
    
    print('Sun Dec \t \t : {}'.format(sundec))
    print('Moon Dec \t \t : {}'.format(moondec))
    print('')
    
    print('Sun Azimuth \t \t : {}'.format(sunaz))
    print('Moon Azimuth \t \t : {}'.format(moonaz))
    print('')
    
    print('Sun Altitude \t \t : {}'.format(sunal))
    print('Moon Altitude \t \t : {}'.format(moonal))
    print('')
    
    print('Relative Altitude \t : {}'.format(Angle(degrees=relal).dstr()))
    print('Relative Azimuth \t : {}'.format(Angle(degrees=relaz).dstr()))
    print('')
    
    print('Moon Elongation \t : {}'.format(elongation))
    print('Moon Illumination \t : {:02.2f}%'.format(illumination))
    
    if plot:
        lingkaranMatahari = plt.Circle((sunaz.to(u.deg)[0].value, sunal.to(u.deg)[0].value), 
                        (0.5*u.deg).value,
                        fc='yellow')
        lingkaranBulan = plt.Circle((moonaz.to(u.deg)[0].value, moonal.to(u.deg)[0].value), 
                                 (0.5*u.deg).value,
                                 fc='black', alpha=.5)
        plt.figure(figsize=(15, 15))
        ax = plt.subplot(aspect=1)
        ax.add_patch(lingkaranMatahari)
        ax.add_patch(lingkaranBulan)
        if sunaz.to(u.deg)[0] <= moonaz.to(u.deg)[0]:
            xlimlow = sunaz.to(u.deg)[0].value
            xlimhigh = moonaz.to(u.deg)[0].value
        else:
            xlimlow = moonaz.to(u.deg)[0].value
            xlimhigh = sunaz.to(u.deg)[0].value

        if sunal.to(u.deg)[0] <= moonal.to(u.deg)[0]:
            ylimlow = sunal.to(u.deg)[0].value
            ylimhigh = moonal.to(u.deg)[0].value
        else:
            ylimlow = moonal.to(u.deg)[0].value
            ylimhigh = sunal.to(u.deg)[0].value

        plt.xlim(xlimlow-2, xlimhigh+2)
        plt.ylim(ylimlow-2, ylimhigh+2)
        plt.text(sunaz.to(u.deg)[0].value-0.28, sunal.to(u.deg)[0].value-0.05, 'Matahari', fontsize=15)
        plt.text(moonaz.to(u.deg)[0].value-0.2, moonal.to(u.deg)[0].value-0.05, 'Bulan', fontsize=15, color='white')
        plt.axhline(y=0)

        plt.plot([moonaz.to(u.deg)[0].value, sunaz.to(u.deg)[0].value], 
                 [moonal.to(u.deg)[0].value, moonal.to(u.deg)[0].value],
                 label='$\Delta$ Azimut = {}'.format(Angle(degrees=relaz).dstr()))
        plt.plot([sunaz.to(u.deg)[0].value, sunaz.to(u.deg)[0].value],
                 [moonal.to(u.deg)[0].value, sunal.to(u.deg)[0].value],
                 label='$\Delta$ Altitud = {}'.format(Angle(degrees=relal).dstr()))
        plt.plot([moonaz.to(u.deg)[0].value, sunaz.to(u.deg)[0].value], 
                 [moonal.to(u.deg)[0].value, sunal.to(u.deg)[0].value],
                 label='Elongasi = {}'.format(elongation))


        plt.xlabel('Azimut (derajat)')
        plt.ylabel('Altitud (derajat)')
        plt.legend()
        plt.title(sunset)

Prediksi Hilal

Syaban 1440H

ucupkiyah(datetime(2019, 4, 5), True)

By the Name of Allah
Ucupkiyah Calculator by Muhammad Yusuf

G. Conjunction Time  : 2019-04-05 15:50:28.624000+07:00
Sunset Time          : 2019-04-05 17:52:48.034000+07:00
Moonset Time         : 2019-04-05 18:02:37.218000+07:00

G. Moon Age      : 02h 02m 19.41s
Lag Time         : 00h 09m 49.18s

Sun RA           : 00h 56m 36.75s
Moon RA          : 01h 04m 12.04s

Sun Dec          : +06deg 03' 05.6"
Moon Dec         : +01deg 56' 37.8"

Sun Azimuth          : 275deg 59' 42.8"
Moon Azimuth         : 272deg 08' 33.2"

Sun Altitude         : -00deg 49' 59.9"
Moon Altitude        : 01deg 32' 08.7"

Relative Altitude    : 02deg 22' 08.6"
Relative Azimuth     : -03deg 51' 09.6"

Moon Elongation      : 04deg 31' 21.0"
Moon Illumination    : 0.20%

Simulasi Hilal Sya'ban 1440H
Simulasi Hilal Sya'ban 1440H

Ramadhan 1440H

ucupkiyah(datetime(2019, 5, 5), True, 0)

By the Name of Allah
Ucupkiyah Calculator by Muhammad Yusuf

G. Conjunction Time  : 2019-05-05 05:45:29.812000+07:00
Sunset Time          : 2019-05-05 17:41:46.344000+07:00
Moonset Time         : 2019-05-05 18:08:49.400000+07:00

G. Moon Age      : 11h 56m 16.53s
Lag Time         : 00h 27m 3.06s

Sun RA           : 02h 48m 47.59s
Moon RA          : 03h 13m 50.25s

Sun Dec          : +16deg 14' 05.4"
Moon Dec         : +13deg 52' 39.4"

Sun Azimuth          : 286deg 15' 06.9"
Moon Azimuth         : 284deg 43' 23.0"

Sun Altitude         : -00deg 49' 59.9"
Moon Altitude        : 05deg 28' 23.6"

Relative Altitude    : 06deg 18' 23.5"
Relative Azimuth     : -01deg 31' 43.9"

Moon Elongation      : 06deg 29' 19.4"
Moon Illumination    : 0.42%

Simulasi Kondisi Hilal Ramadhan 1440H
Simulasi Kondisi Hilal Ramadhan 1440H

Syawal 1440H

ucupkiyah(datetime(2019, 6, 3), True, 0)

By the Name of Allah
Ucupkiyah Calculator by Muhammad Yusuf

G. Conjunction Time  : 2019-06-03 17:01:56.654000+07:00
Sunset Time          : 2019-06-03 17:40:01.059000+07:00
Moonset Time         : 2019-06-03 17:40:28.707000+07:00

G. Moon Age      : 00h 38m 4.40s
Lag Time         : 00h 00m 27.65s

Sun RA           : 04h 44m 30.26s
Moon RA          : 04h 43m 28.73s

Sun Dec          : +22deg 18' 09.4"
Moon Dec         : +19deg 30' 21.5"

Sun Azimuth          : 292deg 21' 54.8"
Moon Azimuth         : 289deg 33' 36.1"

Sun Altitude         : -00deg 49' 59.9"
Moon Altitude        : -00deg 43' 46.0"

Relative Altitude    : 00deg 06' 13.9"
Relative Azimuth     : -02deg 48' 18.8"

Moon Elongation      : 02deg 48' 24.7"
Moon Illumination    : 0.07%

Simulasi Kondisi Hilal Syawal 1440H
Simulasi Kondisi Hilal Syawal 1440H

Jupyter Notebook