Aviator on the BBC Micro

Maths: Sine16Bit

```       Name: Sine16Bit                                               [Show more]
Type: Subroutine
Category: Maths
Summary: Calculate the sine of a 16-bit number
Deep dive: Trigonometry
Context: See this subroutine in context in the source code
References: This subroutine is called as follows:
* ProjectAxisAngle calls Sine16Bit

Sets (A Y) = sin(X W)

where (X W) is a 16-bit angle, with 0 to 65535 representing a quarter circle.

The sine lookup table contains 256 entries that represent a quarter circle, so
we work out the result by first looking up the sine for the high byte X, and
then interpolating W/256 of the way between the results for X and X + 1. It
might help to think of X being an integer (0 to 255) and W being a fraction
(0 to 255) and we're trying to map X.W onto the sine table by finding the
entry for X and doing linear interpolation between X and X + 1 for the
fractional amount.

Arguments:

(X W)                The 16-bit angle, representing a quarter-circle in the
range 0 to 255

Returns:

(A Y)                sin(X W)

(X W)                Unchanged

.Sine16Bit

STX H                  \ Store X in H for later

SEC                    \ Set (A I) = sin(X + 1) - sin(X)
LDA sinLo+1,X          \
SBC sinLo,X            \ starting with the low bytes
STA I

LDA sinHi+1,X          \ And then the high bytes
SBC sinHi,X

\ Let's call this dX, the difference between sin(X + 1)
\ sin(X) and, so (A I) = dX

LSR A                  \ The maximum differential value we can have is 402, so
ROR I                  \ we halve (A I) to fit the result into one byte, like
LDX I                  \ this:
\
\   (A I) = dX / 2
\
\ and then setting X to the low byte in I, so:
\
\   X = dX / 2

LDY W                  \ Set Y = W

JSR Multiply8x8        \ Set (A V) = X * Y
\           = (dX / 2) * W

LDX H                  \ Set X to the original argument we stored above

ASL V                  \ Set (A V) = (A V) * 2
ROL A                  \           = (dX / 2) * W * 2
\           = dX * W
\
\ and put bit 7 of (A V) into the C flag, so in effect
\ we have:
\
\   (C A V) = dX * W

PHP                    \ Store the C flag on the stack

CLC                    \ Set (A Y) = (C A) + sin(X)
ADC sinLo,X            \           = (dX * W / 256) + sin(X)
TAY                    \
\ starting with the low bytes

LDA #0                 \ And then adding the carry to the high byte, so now we
ADC sinHi,X            \ have the interim result for (0 A) + sin(X)

PLP                    \ And finally adding C (which we retrieve from the
ADC #0                 \ stack) to give the final result for (C A) + sin(X)

\ So at this point we have:
\
\   (A Y) = sin(X) + (dX * W / 256)
\
\ which is the result we want, as it takes sin(X) and
\ adds on W/256 of the difference between sin(X) and
\ sin(X + 1)

RTS                    \ Return from the subroutine
```