Wishful Coding

Didn't you ever wish your
computer understood you?

Simulating the solar system

I was reading an old book about astrology, which contained a section about how to calculate your horoscope. It used crude lookup tables for the planets, and it was only valid until 2000.

I had also wondered if two planets were still considered conjunct if they where above/below each other.

So I set out to compute my horoscope in 3D using real orbits instead of lookup tables. I read some things on Wikipedia, but nothing I could translate to running code.

Then I finally found a NASA page that contains all the orbital parameters and a PDF explaining the math.

It took me a long time to get it working for two reasons. The first was solving Kepler’s equation to obtain the eccentric anomaly.

They mix and match degrees and radians in confusing ways. While it might make sense to take the cosine of a value in degrees, it certainly does not make sense to math.cos.

I suggest any math code to convert degrees to radians at the edges, and deal with radians only. Degrees should only be used for input and display. They make sense to humans.

Once you have an approximation of Kepler’s formula solved for eccentric anomaly, it is good to verify that the original formula returns your mean anomaly.

def solve_kepler(eccentricity, mean_anomaly):
    # for the approximate formulae in the present context, tol = 10e-6 degrees is sufficient
    tolerance = 10e-6
    # E0 = M + e sin M
    eccentric_anomaly = mean_anomaly + (eccentricity * math.sin(mean_anomaly))
    # and itterate the following equations with n = 0,1,2,... unitl |delta E| <= tol
    while True:
        # delta M = M - (En - e sin En)
        delta_mean_anomaly = mean_anomaly - (eccentric_anomaly - (eccentricity * math.sin(eccentric_anomaly)))
        # delta E = delta M / (1 - e cos En)
        delta_eccentric_anomaly = delta_mean_anomaly / (1 - (eccentricity * math.cos(eccentric_anomaly)))
        # En+1 = En + delta E
        eccentric_anomaly += delta_eccentric_anomaly

        if abs(delta_eccentric_anomaly) <= tolerance:
            return eccentric_anomaly

The second problem was rotating the planet’s ellipse into the Sun’s coordinate system. Written out in Python it’s a screen full of trigonometric functions. Missing anything gives completely bogus results. It’s not hard, just tedious and error prone.

def ecliptic_coordinates(orbit_x, orbit_y, perihelion, longitude_ascending, inclination):
    term1 = math.cos(perihelion) * math.cos(longitude_ascending)
    term2 = math.sin(perihelion) * math.sin(longitude_ascending) * math.cos(inclination)
    term3 = math.sin(perihelion) * math.cos(longitude_ascending)
    term4 = math.cos(perihelion) * math.sin(longitude_ascending) * math.cos(inclination)
    x = (term1 - term2) * orbit_x + (-term3 - term4) * orbit_y

    term1 = math.cos(perihelion) * math.sin(longitude_ascending)
    term2 = math.sin(perihelion) * math.cos(longitude_ascending) * math.cos(inclination)
    term3 = math.sin(perihelion) * math.sin(longitude_ascending)
    term4 = math.cos(perihelion) * math.cos(longitude_ascending) * math.cos(inclination)
    y = (term1 + term2) * orbit_x + (-term3 + term4) * orbit_y

    term1 = math.sin(perihelion) * math.sin(inclination)
    term2 = math.cos(perihelion) * math.sin(inclination)
    z = term1 * orbit_x + term2 * orbit_y

    return x, y, z

And in the end I did not even use this code. Once I started using Processing, I also wanted to draw the orbits themselves. This should not be hard, given the semi-major axis and eccentricity.

I again hit two problems. The first is that Processing’s ellipse function takes a center point, height and width. But I solved this using Wikipedia and math I can understand.

def planet_ellipse(planet, jd):
    eccentricity = planet.eccentricity(jd)
    smajor = planet.semi_major_axis(jd)
    sminor = smajor * math.sqrt(1 - (eccentricity ** 2))
    center = smajor * eccentricity
    return smajor * 2, sminor * 2, center

The tougher problem is rotating the ellipse to match the planet. While I successfully translated the math to code, I had no idea what is going on. I still have no clue how this pile of trigonometry works.

Instead I substituted the function arguments to ecliptic_coordinates for mouseX and mouseY to experiment what is going on. After some failed experiments I found some rotation commands that would align the ellipse with the planet. Then I removed the trigonometric mess and simply drew the planet in the same plane as the ellipse.

def planet_rotation(planet, jd):
    perihelion = planet.longitude_perhelion(jd) - planet.longitude_ascending(jd)
    rotateZ(-planet.longitude_ascending(jd))
    rotateX(planet.inclination(jd))
    rotateZ(-perihelion)

That’s it. The rest is mostly just drawing code. You can run the code for yourself.

I really like how it turned out. If you pan around you can clearly see that Pluto might not be as conjunct as it looks in 2D. If you let things rotate, you can visually see Mercury speed up near the sun. It’s also intuitively clear how retrograde motion works.

The only thing missing is the moon. I could not find orbital parameters for it.

Published on

31c3 Lightning Talk

I’ve been putting off writing this post for far to long, because I’m to busy hacking on other things. So here are the video, slides, and code. Enjoy.

Slides

import usb.core
import binascii
import re

dev = usb.core.find(idProduct=0x001e, idVendor=0x0b0c)
dev.set_configuration()

def hex_print(s):
    print(binascii.hexlify(s).decode(), re.sub("[^a-zA-Z0-9]", ".", str(s, 'ascii', 'replace')))

def write(data):
    s = binascii.unhexlify(data)
    hex_print(s)
    dev.write(0x02, s)

def read(t=100):
    try:
        while True:
            s = dev.read(0x81, 8, t)
            hex_print(s)
            t = 100
    except usb.core.USBError:
        pass
    print('---')

def message(data, display=0x01):
    #write = print
    data = chr(display).encode() + chr(len(data)).encode() + data
    write(b"010305" + binascii.hexlify(chr(len(data)).encode()) + b"00000000")
    for i in range(int(len(data)/6)):
        s = data[i*6:i*6+6]
        b = binascii.hexlify(s)
        write(b"0006" + b)

    i += 1
    s = data[i*6:i*6+6]
    b = binascii.hexlify(s.ljust(6))
    write(b'0004' + b)

# always the same lenth
def confirm_login(data, lang='nl'):
    data = b'\x03' + data
    write(b"0103081500000000")
    for i in range(3):
        s = data[i*6:i*6+6]
        b = binascii.hexlify(s.ljust(6))
        write(b"0006" + b)

    write(b'000300' + binascii.hexlify(lang) + b'000000')

if __name__ == '__main__':
    write(b"0209000000000000") # shield
    read()
    write(b"0103020000000000") # version
    read()
    write(b"0103010200000000") # insert card
    write(b"00026e6c00000000")
    read(10000)
    write(b"0103030000000000") # card info
    read()
    write(b"0103040000000000") # ask pin
    read(60000)
    message(b'abbalalalala', 0x00) # sign data
    read()

    message(b"Never gonna give you up Never gonna let you down                    ")
    read(10000)

    write(b"0103060000000000") # cryptogram
    read()

    confirm_login(b'You where drunk', b'en')
    read()
Published on

Branch-free FizzBuzz in Assembly

I came across this post that discusses ways to write FizzBuzz in Clojure without using conditionals. However, most if not all of the solutions still do a lot of branching behind the scenes. Think of hash lookups for example.

So I asked to myself, how can I write a FizzBuzz solution with no branches at all? Probably not in Clojure; you can’t easily tell where it is branching or not.

The only way to be absolutely sure is to write it in assembly. So I did. I never did assembly before, so it might be terrible code.

I used an array of 15 pointers to either “fizz”, “buzz”, “fizzbuzz”, or a number buffer. I then filled the number buffer with the current number in ascii and printed whatever I get from the array.

One thing I struggeled with is how to stop. At first I had one condition to see if I reached 100. Now I use a lookup table that calls sys_time 99 times and then sys_exit.

section .data
  f db "fizz    "
  b db "buzz    "
  fb db "fizzbuzz"
  n db 10 ; newline string
  cycle dq num, num, f, num, b, f, num, num, f, b, num, f, num, num, fb
  callid dq 13, 1 ; sys_time, sys_exit

section .bss
  num resb 8 ; number buffer

section .text
global _start

print:      ;write rcx
  mov rax,4 ;sys_write
  mov rbx,1 ;stdout
  mov rdx,8
  int 0x80
  ret

newline:    ;write newline
  mov rax,4 ;sys_write
  mov rbx,1 ;stdout
  mov rcx,n
  mov rdx,1
  int 0x80
  ret

itoa:      ;convert rax to str
  mov byte[num+1],0x30
  mov rdx,0
  mov rcx,10
  div rcx
  add [num+1],rdx
  mov byte[num],0x30
  mov rdx,0
  mov rcx,10
  div rcx
  add [num],rdx
  ret

_start:
  mov r12,0

hundredtimes:
  ; initialise number buffer
  mov rax,r12
  inc rax
  call itoa

  ; mod 15 the number
  mov rax,r12
  mov rdx,0
  mov rcx,15
  div rcx

  ; look up the number in cycle
  ; prints the num buffer or any of the strings
  mov rcx,[cycle+rdx*8]
  call print
  call newline
  
  ; next...
  inc r12

  ; devide the number by 100
  mov rax,r12
  mov rdx,0
  mov rcx,100
  div rcx

  ; get the time or exit
  mov rax,[callid+rax*8]
  int 0x80

  ; jump to the top of the loop
  jmp hundredtimes

To compile on a 64 bit machine:

nasm -f elf64 fizzbuzz.asm
ld -o fb fizzbuzz.o
./fb
Published on