#COMPILE EXE
#DIM ALL

FUNCTION PBMAIN() AS LONG
    LOCAL q1_w AS DOUBLE
    LOCAL q1_x AS DOUBLE
    LOCAL q1_y AS DOUBLE
    LOCAL q1_z AS DOUBLE
    LOCAL q2_w AS DOUBLE
    LOCAL q2_x AS DOUBLE
    LOCAL q2_y AS DOUBLE
    LOCAL q2_z AS DOUBLE
    LOCAL r_w AS DOUBLE
    LOCAL r_x AS DOUBLE
    LOCAL r_y AS DOUBLE
    LOCAL r_z AS DOUBLE
    LOCAL conj_w AS DOUBLE
    LOCAL conj_x AS DOUBLE
    LOCAL conj_y AS DOUBLE
    LOCAL conj_z AS DOUBLE
    LOCAL v_x AS DOUBLE
    LOCAL v_y AS DOUBLE
    LOCAL v_z AS DOUBLE
    LOCAL v_rot_x AS DOUBLE
    LOCAL v_rot_y AS DOUBLE
    LOCAL v_rot_z AS DOUBLE
    LOCAL len_q AS DOUBLE
    LOCAL temp_w AS DOUBLE
    LOCAL temp_str AS STRING

    PRINT "Quaternion Operations"
    PRINT "===================="
    PRINT ""

    ' Q1 = 1 + 2i + 3j + 4k
    q1_w = 1.0
    q1_x = 2.0
    q1_y = 3.0
    q1_z = 4.0

    ' Q2 = 5 + 6i + 7j + 8k
    q2_w = 5.0
    q2_x = 6.0
    q2_y = 7.0
    q2_z = 8.0

    ' Multiply q1 * q2
    r_w = q1_w * q2_w - q1_x * q2_x - q1_y * q2_y - q1_z * q2_z
    r_x = q1_w * q2_x + q1_x * q2_w + q1_y * q2_z - q1_z * q2_y
    r_y = q1_w * q2_y - q1_x * q2_z + q1_y * q2_w + q1_z * q2_x
    r_z = q1_w * q2_z + q1_x * q2_y - q1_y * q2_x + q1_z * q2_w
    PRINT "Q1 * Q2 = (" + STR$(r_w) + ") + (" + STR$(r_x) + ")i + (" + STR$(r_y) + ")j + (" + STR$(r_z) + ")k"

    ' Conjugate
    conj_w = q1_w
    conj_x = -q1_x
    conj_y = -q1_y
    conj_z = -q1_z
    PRINT "Conjugate(Q1) = (" + STR$(conj_w) + ") + (" + STR$(conj_x) + ")i + (" + STR$(conj_y) + ")j + (" + STR$(conj_z) + ")k"

    ' Norm
    len_q = SQR(q1_w * q1_w + q1_x * q1_x + q1_y * q1_y + q1_z * q1_z)
    PRINT "|Q1| = " + STR$(len_q)

    ' Rotate vector (1,0,0) by 90 degrees around Y axis using quaternion
    ' q = cos(45) + 0i + sin(45)j + 0k
    q1_w = 0.70710678
    q1_x = 0.0
    q1_y = 0.70710678
    q1_z = 0.0

    v_x = 1.0
    v_y = 0.0
    v_z = 0.0

    ' v' = q * v * q^(-1)
    ' First: q * v (v as pure quaternion)
    r_w = -q1_x * v_x - q1_y * v_y - q1_z * v_z
    r_x = q1_w * v_x + q1_y * v_z - q1_z * v_y
    r_y = q1_w * v_y - q1_x * v_z + q1_z * v_x
    r_z = q1_w * v_z + q1_x * v_y - q1_y * v_x

    ' conj(q) for inverse (unit quaternion)
    conj_w = q1_w
    conj_x = -q1_x
    conj_y = -q1_y
    conj_z = -q1_z

    ' Then (q*v) * conj(q)
    temp_w = r_w * conj_w - r_x * conj_x - r_y * conj_y - r_z * conj_z
    v_rot_x = r_w * conj_x + r_x * conj_w + r_y * conj_z - r_z * conj_y
    v_rot_y = r_w * conj_y - r_x * conj_z + r_y * conj_w + r_z * conj_x
    v_rot_z = r_w * conj_z + r_x * conj_y - r_y * conj_x + r_z * conj_w

    PRINT ""
    PRINT "Rotate (1,0,0) by 90 around Y axis:"
    PRINT "  Result: (" + STR$(v_rot_x) + ", " + STR$(v_rot_y) + ", " + STR$(v_rot_z) + ")"

    ' Identity quaternion
    q1_w = 1.0
    q1_x = 0.0
    q1_y = 0.0
    q1_z = 0.0
    PRINT ""
    PRINT "Identity quaternion: (1) + (0)i + (0)j + (0)k"

    ' Slerp demo (simplified)
    q1_w = 1.0
    q1_x = 0.0
    q1_y = 0.0
    q1_z = 0.0
    q2_w = 0.0
    q2_x = 1.0
    q2_y = 0.0
    q2_z = 0.0

    ' Dot product
    dot_val = q1_w * q2_w + q1_x * q2_x + q1_y * q2_y + q1_z * q2_z
    PRINT "Dot(Q1, Q2) = " + STR$(dot_val)

    PRINT ""
    PRINT "Quaternion operations complete."

    FUNCTION = 0
END FUNCTION
