def orbit_angle(x,y):
    nt  = len(x)
    ang = np.zeros(nt)
    for it in range(1,nt):
        l0   = np.sqrt(x[it-1]**2+y[it-1]**2)
        l1   = np.sqrt(x[it]**2+y[it]**2)
        dang = np.arcsin(np.abs(x[it]*y[it-1]-y[it]*x[it-1])/(l0*l1))
        assert dang < 1.0, 'Error in orbit angle determination: angular step too large!'
        ang[it] = ang[it-1] + dang
    return ang