Erroneous mathematical processing?
I'm currently doing a project that involves 2 separate pic processors on 2 separate boards doing completely different things. One of the boards is more avionics based and one of its tasks is to process GPS data to establish a required heading and distance to a destination point. I'm using the Cordic tools posted here by Walter a few years ago to help with the trig processing.
Also, I like the idea of transitioning from previous 28 pin PICs to the 18F26K83 due to numerous virtues it has; however one thing I discovered attempting to implement it for this particular application is that the 18F26K83 and the Cordic routines don't see eye-to-eye and produce erroneous results. Well, for some reason the results (particularly for atan2) always appear to exceed the maximum.
If I change to an 18F26K22 with exactly the same code (except the osc settings which are different from the 28K83), the routines work perfectly. This is just straight math processing, no external inputs to factor for this debugging.
Troy
Re: Erroneous mathematical processing?
if thats not just a comment then start here
Re: Erroneous mathematical processing?
Richard,
No, it was just a comment. The solution for me was to use a different 28 pin PIC which I had no problem with.
Nevertheless, I think it's importand to put these issues on public record for anyone else who might experience similar issues. Saying that, I'm happy to post pertinent code and whatever else if anyone is interested in identifying the cause.
Cheers,
Troy
3 Attachment(s)
Re: Erroneous mathematical processing?
Okay here is the code for the PIC18F26K22:
Code:
include "trig18.inc"
Define OSC 16 'Clock will be configured to 16MHz
OSCCON.6 = 1 ' Set Internal Osc to 16Mhz
OSCCON.5 = 1 ' Set Internal Osc to 16Mhz
OSCCON.4 = 1 ' Set Internal Osc to 16Mhz
OSCCON.1 = 0 ' Use Internal Osc through Primary Clock (not Directly) to Access PLL Block
OSCCON.0 = 0 ' Use Internal Osc through Primary Clock (not Directly) to Access PLL Block
OSCTUNE.6 = 0 '16MHz osc - PLL Disabled
'Configure AN2 as Analog Inputs (PortA2)
ANSELA = %00000000 'Everything Digital
ANSELB = %00000000 'Everything Else Digital
ANSELC = %00000000 'Everything Else Digital
TRISB.7 = 0 'Outout
Serial_out var PORTB.7 'Serial out pin for debugging/developing comms
n var byte
main:
pause 2000
serout2 Serial_out,84, ["Test 123",13,10]
Pause 1000
ang=3000
serout2 Serial_out,84, ["Ang: ",dec ang,13,10]
'sincos example:
'input ang in degrees.dd example: ang = 3000 (= 30.00 degrees)
'call sincos
'result: x = 15004 , Y = 25981
'so, 15004/30000 = sin(ang) = 0.5001 and 25981/30000 = cos(ang) = 0.8660
call sincos
serout2 Serial_out,84, ["Sin ang: ",dec x,13,10]
serout2 Serial_out,84, ["Cos ang: ",dec y,13,10]
serout2 Serial_out,84, [13,10]
x = 12000
y = 8000
serout2 Serial_out,84, ["x input: ",sdec x,13,10]
serout2 Serial_out,84, ["y input: ",sdec y,13,10]
call atan2
'result: x = 14422 = hypotenuse or distance to x,y, y = angle to x,y in degrees.dd = 5631 (=56.31 degrees), ang = 6133 = radians 0 to 65535
serout2 Serial_out,84, ["hypotenuse: ",sdec x,13,10]
serout2 Serial_out,84, ["Ang in degrees.dd: ",sdec y/100,".",dec2 y,13,10]
serout2 Serial_out,84, ["Ang in Rad: ",dec ang,13,10]
serout2 Serial_out,84, [13,10]
x = -12000
y = 8000
serout2 Serial_out,84, ["x input: ",sdec x,13,10]
serout2 Serial_out,84, ["y input: ",sdec y,13,10]
call atan2
'result: x = 14422 = hypotenuse or distance to x,y, y = angle to x,y in degrees.dd = 3037 (=303.70 degrees), ang = 26634 = radians 0 to 65535
serout2 Serial_out,84, ["hypotenuse: ",sdec x,13,10]
serout2 Serial_out,84, ["Ang in degrees.dd: ",sdec y/100,".",dec2 y,13,10]
serout2 Serial_out,84, ["Ang in Rad: ",dec ang,13,10]
serout2 Serial_out,84, [13,10]
x = -12000
y = -8000
serout2 Serial_out,84, ["x input: ",sdec x,13,10]
serout2 Serial_out,84, ["y input: ",sdec y,13,10]
call atan2
'result: x = 14422 = hypotenuse or distance to x,y, y = angle to x,y in degrees.dd = 23632 (=236.32 degrees), ang = 38899 = radians 0 to 65535
serout2 Serial_out,84, ["hypotenuse: ",sdec x,13,10]
serout2 Serial_out,84, ["Ang in degrees.dd: ",sdec y/100,".",dec2 y,13,10]
serout2 Serial_out,84, ["Ang in Rad: ",dec ang,13,10]
serout2 Serial_out,84, [13,10]
x = 12000
y = -8000
serout2 Serial_out,84, ["x input: ",sdec x,13,10]
serout2 Serial_out,84, ["y input: ",sdec y,13,10]
call atan2
'result: x = 14422 = hypotenuse or distance to x,y, y = angle to x,y in degrees.dd = 12367 (=123.67 degrees), ang = 59405 = radians 0 to 65535
serout2 Serial_out,84, ["hypotenuse: ",sdec x,13,10]
serout2 Serial_out,84, ["Ang in degrees.dd: ",sdec y/100,".",dec2 y,13,10]
serout2 Serial_out,84, ["Ang in Rad: ",dec ang,13,10]
Goto main
Results of the 18F26K22 code running on a 18F26K22:
Attachment 8950
Here is the (COPY & PASTED) code for the PIC18F26K83 - note the only difference is the osc setting method:
Code:
include "trig18.inc"
'Configure OSC postscaler to provide 16mhz system clock
OSCFRQ = %0101
Define OSC 16
'Configure AN2 as Analog Inputs (PortA2)
ANSELA = %00000000 'Everything Digital
ANSELB = %00000000 'Everything Else Digital
ANSELC = %00000000 'Everything Else Digital
TRISB.7 = 0 'Outout
Serial_out var PORTB.7 'Serial out pin for debugging/developing comms
n var byte
main:
pause 2000
serout2 Serial_out,84, ["Test 123",13,10]
Pause 1000
ang=3000
serout2 Serial_out,84, ["Ang: ",dec ang,13,10]
'sincos example:
'input ang in degrees.dd example: ang = 3000 (= 30.00 degrees)
'call sincos
'result: x = 15004 , Y = 25981
'so, 15004/30000 = sin(ang) = 0.5001 and 25981/30000 = cos(ang) = 0.8660
call sincos
serout2 Serial_out,84, ["Sin ang: ",dec x,13,10]
serout2 Serial_out,84, ["Cos ang: ",dec y,13,10]
serout2 Serial_out,84, [13,10]
x = 12000
y = 8000
serout2 Serial_out,84, ["x input: ",sdec x,13,10]
serout2 Serial_out,84, ["y input: ",sdec y,13,10]
call atan2
'result: x = 14422 = hypotenuse or distance to x,y, y = angle to x,y in degrees.dd = 5631 (=56.31 degrees), ang = 6133 = radians 0 to 65535
serout2 Serial_out,84, ["hypotenuse: ",sdec x,13,10]
serout2 Serial_out,84, ["Ang in degrees.dd: ",sdec y/100,".",dec2 y,13,10]
serout2 Serial_out,84, ["Ang in Rad: ",dec ang,13,10]
serout2 Serial_out,84, [13,10]
x = -12000
y = 8000
serout2 Serial_out,84, ["x input: ",sdec x,13,10]
serout2 Serial_out,84, ["y input: ",sdec y,13,10]
call atan2
'result: x = 14422 = hypotenuse or distance to x,y, y = angle to x,y in degrees.dd = 3037 (=303.70 degrees), ang = 26634 = radians 0 to 65535
serout2 Serial_out,84, ["hypotenuse: ",sdec x,13,10]
serout2 Serial_out,84, ["Ang in degrees.dd: ",sdec y/100,".",dec2 y,13,10]
serout2 Serial_out,84, ["Ang in Rad: ",dec ang,13,10]
serout2 Serial_out,84, [13,10]
x = -12000
y = -8000
serout2 Serial_out,84, ["x input: ",sdec x,13,10]
serout2 Serial_out,84, ["y input: ",sdec y,13,10]
call atan2
'result: x = 14422 = hypotenuse or distance to x,y, y = angle to x,y in degrees.dd = 23632 (=236.32 degrees), ang = 38899 = radians 0 to 65535
serout2 Serial_out,84, ["hypotenuse: ",sdec x,13,10]
serout2 Serial_out,84, ["Ang in degrees.dd: ",sdec y/100,".",dec2 y,13,10]
serout2 Serial_out,84, ["Ang in Rad: ",dec ang,13,10]
serout2 Serial_out,84, [13,10]
x = 12000
y = -8000
serout2 Serial_out,84, ["x input: ",sdec x,13,10]
serout2 Serial_out,84, ["y input: ",sdec y,13,10]
call atan2
'result: x = 14422 = hypotenuse or distance to x,y, y = angle to x,y in degrees.dd = 12367 (=123.67 degrees), ang = 59405 = radians 0 to 65535
serout2 Serial_out,84, ["hypotenuse: ",sdec x,13,10]
serout2 Serial_out,84, ["Ang in degrees.dd: ",sdec y/100,".",dec2 y,13,10]
serout2 Serial_out,84, ["Ang in Rad: ",dec ang,13,10]
Goto main
(Erroneous) Results of the 18F26K83 code running on a 18F26K83:
Attachment 8951
I've also attached my slightly modifed Trig18.inc file - only modified to allow the use of longs for PIC18s. Same inc file is used for both examples.
Re: Erroneous mathematical processing?
i don't have one of those chips or the ability to sim it either but for a start
rather than trying to cram those extra vars into bank0 {unnecessarily} i would modify the include this way
for "LONG" compilation. if you did push the important vars out of bank 0 then it wont fly anymore
Code:
'/*******************************************************************************' TRIG.inc for PIC18 Long Version
'/*******************************************************************************
'* FUNCTION NAME: sincos
'*
'* ARGUMENTS: ang (angle deg.dd format example: 35999 = 359.99 deg)
'*
'* RETURNS: x = sin(ang) , y = cos(ang)
'*
'* DESCRIPTION: The angle is given in degrees.dd (see above)
'* The function simultaneously calculates the sine
'* and cosine of the angle as fractions of 30,000 (where 30,000
'* equates to 1 and -30,000 equates to -1) and returns them in
'* x, and y as x=sin(ang), y = cos(ang).
'*
'* EXAMPLE: ang = 3000 (30.00 degrees) x = 15004 , Y = 25981
'* so, 15004/30000 = 0.5001 and 25981/30000 = 0.8660
'*
'*******************************************************************************/
'/*******************************************************************************
'* FUNCTION NAME: atan2
'*
'* ARGUMENTS: int x (x-coordinate)
'* int y (y-coordinate)
'*
'* RETURNS: atan2 of x,y:(angle to x,y) and the hypotenuse (distance) of x,y
'*
'* DESCRIPTION: Given an ordered pair of coordinates, the function
'* simultaneously calculates the atan2 (the direction of the
'* position as y degrees, and ang radians) and the square root of
'* the sum of the squares of the coordinates (the magnitude of
'* the position vector) as x
'*
'*
'* NOTES: (1) The accuracy of the returned values increases as the
'* sizes of x and y increase. Consider multiplying both by a
'* scaling factor before calling the function.
'* (2) The function will fail for x and y values that result in
'* magnitues greater than 32,767 (the size of a signed int).
'*
'* EXAMPLE: atan2
'* x = 25980, y = 15000;
'* results:
'* ang = angle in radians 0 to 65535
'* x = distance = 30000 (same units as x and y)
'* y = 6000 (angle in degrees 60.00 degrees)
'*******************************************************************************/
;*******************************************************************************
; --- CORDIC TRIG LIBRARY ---
;http://www.chiefdelphi.com/media/papers/2016
; FILE NAME: trig.inc
; AUTHOR: Patrick Fairbank
; Last Modified Feb 15, 2012 - fixed bug in sincos which calculated wrong
; result for sin and cos in quadrants 2,3 and 4. Thanks Martin!
; FEB. 15, 2009 to make it PicBasic compatible and added degree conversion
; Modified by Walter Dunckel (Scale Robotics Inc.) with help from Darrel Taylor
; http://www.scalerobotics.com/79-cordic-for-picbasic.html
; DESCRIPTION: This file contains functions implementing the COORDIC
; algorithm, or how to get a 16 bit sin, cos, tan2 and hypotenuse result
;
; USAGE: Add this file to your PicBasic Pro project using INCLUDE "TRIG.inc"
; and add a line with main: directly below the include line
; Then fill x,y values for atan2, or fill ang value for sincos
; then either CALL sincos or CALL atan2
; LICENSE: Users are free to use, modify, and distribute this code
; as they see fit.
;
; sincos: input ang, output x = sin(ang) and y = cos(ang)
;
;
; atan2: input x and y coordinates
; Output the calculated angle and hypotenuse values
; as output: y = angle in degrees, ang = angle in radians, x = hypotenuse
;
;******************************************************************************/
' Example output using TRIG.inc using hserout
' HSEROUT ["Angle = ",dec a,", sin(ang) = ",sdec x,", cos(ang) = ",sdec y,13]
'Angle = 0, sin(ang) = 5, cos(ang) = 30000
'Angle = 10, sin(ang) = 5209, cos(ang) = 29544
'Angle = 20, sin(ang) = 10266, cos(ang) = 28188
'Angle = 30, sin(ang) = 15004, cos(ang) = 25981
'Angle = 40, sin(ang) = 19287, cos(ang) = 22981
'Angle = 50, sin(ang) = 22983, cos(ang) = 19285
'Angle = 60, sin(ang) = 25983, cos(ang) = 14998
'Angle = 70, sin(ang) = 28188, cos(ang) = 10264
'Angle = 80, sin(ang) = 29544, cos(ang) = 5207
'Angle = 90, sin(ang) = 30000, cos(ang) = 3
'Angle = 100, sin(ang) = 29543, cos(ang) = -5223
'Angle = 110, sin(ang) = 28189, cos(ang) = -10276
'Angle = 120, sin(ang) = 25974, cos(ang) = -15010
'Angle = 130, sin(ang) = 22973, cos(ang) = -19293
'Angle = 140, sin(ang) = 19274, cos(ang) = -22992
'Angle = 150, sin(ang) = 14991, cos(ang) = -25989
'Angle = 160, sin(ang) = 10249, cos(ang) = -28190
'Angle = 170, sin(ang) = 5192, cos(ang) = -29546
'Angle = 180, sin(ang) = -9, cos(ang) = -29999
'Angle = 190, sin(ang) = -5216, cos(ang) = -29544
'Angle = 200, sin(ang) = -10271, cos(ang) = -28184
'Angle = 210, sin(ang) = -15013, cos(ang) = -25977
'Angle = 220, sin(ang) = -19296, cos(ang) = -22976
'Angle = 230, sin(ang) = -22987, cos(ang) = -19275
'Angle = 240, sin(ang) = -25986, cos(ang) = -14992
'Angle = 250, sin(ang) = -28199, cos(ang) = -10258
'Angle = 260, sin(ang) = -29549, cos(ang) = -5195
'Angle = 270, sin(ang) = -30000, cos(ang) = 5
'Angle = 280, sin(ang) = -29544, cos(ang) = 5213
'Angle = 290, sin(ang) = -28188, cos(ang) = 10266
'Angle = 300, sin(ang) = -25981, cos(ang) = 15004
'Angle = 310, sin(ang) = -22981, cos(ang) = 19287
'Angle = 320, sin(ang) = -19283, cos(ang) = 22985
'Angle = 330, sin(ang) = -14998, cos(ang) = 25983
'Angle = 340, sin(ang) = -10260, cos(ang) = 28190
'Angle = 350, sin(ang) = -5207, cos(ang) = 29544
i var byte BANK0
j Var byte BANK0
quad var byte BANK0
x var word bank0
y var word bank0
ang var word bank0
dy var word bank0
dx var word bank0
atans var word[15] bank0
atans(0) = 16384
atans(1) = 9672
atans(2) = 5110
atans(3) = 2594
atans(4) = 1302
atans(5) = 652
atans(6) = 326
atans(7) = 163
atans(8) = 81
atans(9) = 41
atans(10) = 20
atans(11) = 10
atans(12) = 5
atans(13) = 3
atans(14) = 1
goto OverTrig
atan2:
asm
call atan2_sqrt
endasm
'convert to degrees.dd y is degrees
If ang < 16384 then y = 16383 - ang
if ang > 16383 then
y = 65535 - ang
y = y + 16383 'correct 90 degrees for radian
endif
'y = y * 256 'divides radians to get degrees within 57ppm
'y = div32 466 'degrees.dd is y, radians is ang
y = y ** 36002
return
sincos:
'use angle as deg.dd for example 35999 is 359.99 degrees
' ang_temp = ang ' Store ang.
if ang < 9001 then ang = 9000 - ang 'change degrees to radians
if ang > 9000 then ang = 45000 - ang 'change degrees to radians
'ang = ang * 466
'ang = div32 256
ang = ang */ 466
asm
call sin_cos
endasm
' if ang_temp > 9000 then y = ~y + 1 ' Perform 2's complement of y if ang > 90 for cos(ang).
return
asm
; Calculates the sine and cosine of the given angle
sin_cos:
; Initialize _x to 18218
movlw 0x2a
movwf _x
movlw 0x47
movwf _x+1
; Initialize _y to 0
clrf _y
clrf _y+1
; Initialize _quad to 0
clrf _quad
; Check if the angle is greater than 16383 (90°)
sc_check_greaterthan:
btfss _ang+1, 7
btfss _ang+1, 6
bra sc_check_lessthan
bra sc_adjust_quad2
; Check if the angle is less than -16384 (-90°)
sc_check_lessthan:
btfsc _ang+1, 7
btfsc _ang+1, 6
bra sc_setup_end
; If the angle is in quadrant 3, adjust it to quadrant 4
sc_adjust_quad3:
negf _ang
bc sc_negate_quad3
comf _ang+1
bra sc_adjust_end
; If the low byte negation causes a carry, negate the upper byte
sc_negate_quad3:
negf _ang+1
bra sc_adjust_end
; If the angle is in quadrant 2, adjust it to quadrant 1
sc_adjust_quad2:
comf _ang
comf _ang+1
; Toggle the sign bit and set the '_quad' flag
sc_adjust_end:
btg _ang+1, 7
setf _quad
; Multiply the angle by 2 to get better resolution
sc_setup_end:
bcf STATUS, 0
rlcf _ang
rlcf _ang+1
; Set up the main loop
sc_loop_start:
clrf _i
lfsr FSR0, _atans
; The main loop label
sc_loop:
movff _x, _dy
movff _x+1, _dy+1
movff _i, _j
movf _j
bz sc_bs_x_done
; Loop to shift _dy right
sc_bs_x_loop:
bcf STATUS, 0
rrcf _dy+1
rrcf _dy
btfsc _x+1, 7
bsf _dy+1, 7
decfsz _j
bra sc_bs_x_loop
; Calculate what needs to be added to _x
sc_bs_x_done:
movff _y, _dx
movff _y+1, _dx+1
movff _i, _j
movf _j
bz sc_do_rotation
; Loop to shift _dx right
sc_bs_y_loop:
bcf STATUS, 0
rrcf _dx+1
rrcf _dx
btfsc _y+1, 7
bsf _dx+1, 7
decfsz _j
bra sc_bs_y_loop
; Perform adding operations on _x, _y and _ang
sc_do_rotation:
btfss _ang+1, 7
bra sc_sub_angle
; If _ang is negative
movf POSTINC0, W
addwf _ang
movf POSTINC0, W
addwfc _ang+1
movf _dx, W
addwf _x
movf _dx+1, W
addwfc _x+1
movf _dy, W
subwf _y
movf _dy+1, W
subwfb _y+1
bra sc_loop_bottom
; If _ang is positive
sc_sub_angle:
movf POSTINC0, W
subwf _ang
movf POSTINC0, W
subwfb _ang+1
movf _dx, W
subwf _x
movf _dx+1, W
subwfb _x+1
movf _dy, W
addwf _y
movf _dy+1, W
addwfc _y+1
; Increment the counter and exit the loop if done
sc_loop_bottom:
incf _i
movlw 0x0f
cpfseq _i
bra sc_loop
; Negate _x if it was initially in quadrant 2 or 3
sc_finished:
btfss _quad, 7
bra sc_output
negf _x
bc sc_negate_x
comf _x+1
bra sc_output
; If the low byte negation causes a carry, negate the upper byte
sc_negate_x:
negf _x+1
; Output the calculated _x and _y values
sc_output:
return
; Calculates the magnitude and direction of the given ordered pair
atan2_sqrt:
; Initialize _ang to 0
clrf _ang
clrf _ang+1
; Initialize _quad to 0
clrf _quad
; If the point is in quadrant 2 or 3, make _x positive and set flag
as_check_negative:
btfss _x+1, 7
bra as_shift_x
setf _quad
negf _x
bc as_negate_x
comf _x+1
bra as_shift_x
; If the low byte negation causes a carry, negate the upper byte
as_negate_x:
negf _x+1
; Divide the _x coordinate by 2 to prevent overflowing
as_shift_x:
bcf STATUS, 0
rrcf _x+1
rrcf _x
; Divide the _y coordinate by 2 to prevent overflowing
as_shift_y:
bcf STATUS, 0
rrcf _y+1
rrcf _y
btfsc _y+1, 6
bsf _y+1, 7
; Set up the main loop
as_loop_start:
clrf _i
lfsr FSR0, _atans
; The main loop label
as_loop:
movff _x, _dy
movff _x+1, _dy+1
movff _i, _j
movf _j
bz as_bs_x_done
; Loop to shift _dy right
as_bs_x_loop:
bcf STATUS, 0
rrcf _dy+1
rrcf _dy
btfsc _x+1, 7
bsf _dy+1, 7
decfsz _j
bra as_bs_x_loop
; Calculate what needs to be added to _x
as_bs_x_done:
movff _y, _dx
movff _y+1, _dx+1
movff _i, _j
movf _j
bz as_do_rotation
; Loop to shift _dx right
as_bs_y_loop:
bcf STATUS, 0
rrcf _dx+1
rrcf _dx
btfsc _y+1, 7
bsf _dx+1, 7
decfsz _j
bra as_bs_y_loop
; Perform adding operations on _x, _y and _ang, shifting the _atans right one
as_do_rotation:
movff POSTINC0, PRODL
movff POSTINC0, PRODH
bcf STATUS, 0
rrcf PRODH
rrcf PRODL
btfsc _y+1, 7
bra as_sub_angle
; If _y is positive
movf PRODL, W
addwf _ang
movf PRODH, W
addwfc _ang+1
movf _dx, W
addwf _x
movf _dx+1, W
addwfc _x+1
movf _dy, W
subwf _y
movf _dy+1, W
subwfb _y+1
bra as_loop_bottom
; If _y is negative
as_sub_angle:
movf PRODL, W
subwf _ang
movf PRODH, W
subwfb _ang+1
movf _dx, W
subwf _x
movf _dx+1, W
subwfb _x+1
movf _dy, W
addwf _y
movf _dy+1, W
addwfc _y+1
; Increment the counter and exit the loop if done
as_loop_bottom:
incf _i
movlw 0x0e
cpfseq _i
bra as_loop
; Multiply the _x value by 19898 and divide by 2^14 to scale it
as_scale_x:
movff _x, _dx
movff _x+1, _dx+1
movlw 0xba
mulwf _dx
movff PRODH, _x
movlw 0x4d
mulwf _dx+1
movff PRODH, _dy
movff PRODL, _x+1
movlw 0xba
mulwf _dx+1
movf PRODL, W
addwf _x, F
movf PRODH, W
addwfc _x+1, F
clrf WREG
addwfc _dy, F
movlw 0x4d
mulwf _dx
movf PRODL, W
addwf _x, F
movf PRODH, W
addwfc _x+1, F
clrf WREG
addwfc _dy, F
movlw 0x06
movwf _j
as_scale_bs_loop:
bcf STATUS, 0
rrcf _dy
rrcf _x+1
rrcf _x
decfsz _j
bra as_scale_bs_loop
; Check if the quadrant was originally changed
as_check_quad:
btfss _quad, 7
bra as_output
btfss _ang+1,7
bra as_adjust_quad1
; If the angle is in quadrant 4, adjust it to quadrant 3
as_adjust_quad4:
negf _ang
bc as_negate_quad4
comf _ang+1
bra as_adjust_end
; If the low byte negation causes a carry, negate the upper byte
as_negate_quad4:
negf _ang+1
bra as_adjust_end
; If the angle is in quadrant 1, adjust it to quadrant 2
as_adjust_quad1:
comf _ang
comf _ang+1
; Toggle the sign bit
as_adjust_end:
btg _ang+1, 7
; Output the calculated angle and hypotenuse values
as_output:
return
endasm
OverTrig: 'jump over code
1 Attachment(s)
Re: Erroneous mathematical processing?
Richard,
Thanks for the tip. I've modified my trig file accordingly. Alas though, it didn't fix the issue with the ...K83 though.
Attachment 8953
Troy
Re: Erroneous mathematical processing?
maybe a config issue , there is not much i can see that would not work in the asm stuff unless the 8x8 mult works differently
you could try this , clutching at straws
the code just assumes bank0
Code:
atan2_sqrt:
banksel 0
; Initialize _ang to 0
clrf _ang
clrf _ang+1
.......
Code:
sin_cos:
banksel 0
; Initialize _x to 18218
movlw 0x2a
movwf _x
movlw 0x47
movwf _x+1
...............
Re: Erroneous mathematical processing?
The K83 has its SFR registers located at the top of ram, and the common ones are in bank 63 ($3F00).
The MOVFF instruction only uses 12-bits of addressing, so the highest address it can reach is $0FFF.
To reach the upper registers past that you have to use the MOVFFL instruction instead. It's a 3-word instruction
(vs 2 for MOVFF) that can reach the entire ram space. I mention that since that could require changing any
jumps or gotos in the asm code as well.
This would effect any code like the following:
Code:
movff POSTINC0, PRODL
movff POSTINC0, PRODH
For the K83 (and K42) that would have to be:
Code:
movffl POSTINC0, PRODL
movffl POSTINC0, PRODH
You also have to be using a version of MPASM that supports this instruction (MPASM v5.71 or later)
Re: Erroneous mathematical processing?
Thanks tumbleweed
Troy try this
Code:
'/*******************************************************************************
' TRIG.inc for ALL PIC18's with or without longs
'/*******************************************************************************
'* FUNCTION NAME: sincos
'*
'* ARGUMENTS: ang (angle deg.dd format example: 35999 = 359.99 deg)
'*
'* RETURNS: x = sin(ang) , y = cos(ang)
'*
'* DESCRIPTION: The angle is given in degrees.dd (see above)
'* The function simultaneously calculates the sine
'* and cosine of the angle as fractions of 30,000 (where 30,000
'* equates to 1 and -30,000 equates to -1) and returns them in
'* x, and y as x=sin(ang), y = cos(ang).
'*
'* EXAMPLE: ang = 3000 (30.00 degrees) x = 15004 , Y = 25981
'* so, 15004/30000 = 0.5001 and 25981/30000 = 0.8660
'*
'*******************************************************************************/
'/*******************************************************************************
'* FUNCTION NAME: atan2
'*
'* ARGUMENTS: int x (x-coordinate)
'* int y (y-coordinate)
'*
'* RETURNS: atan2 of x,y:(angle to x,y) and the hypotenuse (distance) of x,y
'*
'* DESCRIPTION: Given an ordered pair of coordinates, the function
'* simultaneously calculates the atan2 (the direction of the
'* position as y degrees, and ang radians) and the square root of
'* the sum of the squares of the coordinates (the magnitude of
'* the position vector) as x
'*
'*
'* NOTES: (1) The accuracy of the returned values increases as the
'* sizes of x and y increase. Consider multiplying both by a
'* scaling factor before calling the function.
'* (2) The function will fail for x and y values that result in
'* magnitues greater than 32,767 (the size of a signed int).
'*
'* EXAMPLE: atan2
'* x = 25980, y = 15000;
'* results:
'* ang = angle in radians 0 to 65535
'* x = distance = 30000 (same units as x and y)
'* y = 6000 (angle in degrees 60.00 degrees)
'*******************************************************************************/
;*******************************************************************************
; --- CORDIC TRIG LIBRARY ---
;http://www.chiefdelphi.com/media/papers/2016
; FILE NAME: trig.inc
; AUTHOR: Patrick Fairbank
; Last Modified Feb 15, 2012 - fixed bug in sincos which calculated wrong
; result for sin and cos in quadrants 2,3 and 4. Thanks Martin!
; FEB. 15, 2009 to make it PicBasic compatible and added degree conversion
; Modified by Walter Dunckel (Scale Robotics Inc.) with help from Darrel Taylor
; http://www.scalerobotics.com/79-cordic-for-picbasic.html
; DESCRIPTION: This file contains functions implementing the COORDIC
; algorithm, or how to get a 16 bit sin, cos, tan2 and hypotenuse result
;
; USAGE: Add this file to your PicBasic Pro project using INCLUDE "TRIG.inc"
; and add a line with main: directly below the include line
; Then fill x,y values for atan2, or fill ang value for sincos
; then either CALL sincos or CALL atan2
; LICENSE: Users are free to use, modify, and distribute this code
; as they see fit.
;
; sincos: input ang, output x = sin(ang) and y = cos(ang)
;
;
; atan2: input x and y coordinates
; Output the calculated angle and hypotenuse values
; as output: y = angle in degrees, ang = angle in radians, x = hypotenuse
;
;******************************************************************************/
' Example output using TRIG.inc using hserout
' HSEROUT ["Angle = ",dec a,", sin(ang) = ",sdec x,", cos(ang) = ",sdec y,13]
'Angle = 0, sin(ang) = 5, cos(ang) = 30000
'Angle = 10, sin(ang) = 5209, cos(ang) = 29544
'Angle = 20, sin(ang) = 10266, cos(ang) = 28188
'Angle = 30, sin(ang) = 15004, cos(ang) = 25981
'Angle = 40, sin(ang) = 19287, cos(ang) = 22981
'Angle = 50, sin(ang) = 22983, cos(ang) = 19285
'Angle = 60, sin(ang) = 25983, cos(ang) = 14998
'Angle = 70, sin(ang) = 28188, cos(ang) = 10264
'Angle = 80, sin(ang) = 29544, cos(ang) = 5207
'Angle = 90, sin(ang) = 30000, cos(ang) = 3
'Angle = 100, sin(ang) = 29543, cos(ang) = -5223
'Angle = 110, sin(ang) = 28189, cos(ang) = -10276
'Angle = 120, sin(ang) = 25974, cos(ang) = -15010
'Angle = 130, sin(ang) = 22973, cos(ang) = -19293
'Angle = 140, sin(ang) = 19274, cos(ang) = -22992
'Angle = 150, sin(ang) = 14991, cos(ang) = -25989
'Angle = 160, sin(ang) = 10249, cos(ang) = -28190
'Angle = 170, sin(ang) = 5192, cos(ang) = -29546
'Angle = 180, sin(ang) = -9, cos(ang) = -29999
'Angle = 190, sin(ang) = -5216, cos(ang) = -29544
'Angle = 200, sin(ang) = -10271, cos(ang) = -28184
'Angle = 210, sin(ang) = -15013, cos(ang) = -25977
'Angle = 220, sin(ang) = -19296, cos(ang) = -22976
'Angle = 230, sin(ang) = -22987, cos(ang) = -19275
'Angle = 240, sin(ang) = -25986, cos(ang) = -14992
'Angle = 250, sin(ang) = -28199, cos(ang) = -10258
'Angle = 260, sin(ang) = -29549, cos(ang) = -5195
'Angle = 270, sin(ang) = -30000, cos(ang) = 5
'Angle = 280, sin(ang) = -29544, cos(ang) = 5213
'Angle = 290, sin(ang) = -28188, cos(ang) = 10266
'Angle = 300, sin(ang) = -25981, cos(ang) = 15004
'Angle = 310, sin(ang) = -22981, cos(ang) = 19287
'Angle = 320, sin(ang) = -19283, cos(ang) = 22985
'Angle = 330, sin(ang) = -14998, cos(ang) = 25983
'Angle = 340, sin(ang) = -10260, cos(ang) = 28190
'Angle = 350, sin(ang) = -5207, cos(ang) = 29544
#IFDEF MOVFFL
DEFINE MOXVE MOVFFL
#ELSE
DEFINE MOXVE MOVFF
#ENDIF
i var byte BANK0
j Var byte BANK0
quad var byte BANK0
x var word bank0
y var word bank0
ang var word bank0
dy var word bank0
dx var word bank0
atans var word[15] bank0
atans(0) = 16384
atans(1) = 9672
atans(2) = 5110
atans(3) = 2594
atans(4) = 1302
atans(5) = 652
atans(6) = 326
atans(7) = 163
atans(8) = 81
atans(9) = 41
atans(10) = 20
atans(11) = 10
atans(12) = 5
atans(13) = 3
atans(14) = 1
goto OverTrig
atan2:
asm
call atan2_sqrt
endasm
'convert to degrees.dd y is degrees
If ang < 16384 then y = 16383 - ang
if ang > 16383 then
y = 65535 - ang
y = y + 16383 'correct 90 degrees for radian
endif
'y = y * 256 'divides radians to get degrees within 57ppm
'y = div32 466 'degrees.dd is y, radians is ang
y = y ** 36002
return
sincos:
'use angle as deg.dd for example 35999 is 359.99 degrees
' ang_temp = ang ' Store ang.
if ang < 9001 then ang = 9000 - ang 'change degrees to radians
if ang > 9000 then ang = 45000 - ang 'change degrees to radians
'ang = ang * 466
'ang = div32 256
ang = ang */ 466
asm
call sin_cos
endasm
' if ang_temp > 9000 then y = ~y + 1 ' Perform 2's complement of y if ang > 90 for cos(ang).
return
asm
; Calculates the sine and cosine of the given angle
sin_cos:
; Initialize _x to 18218
movlw 0x2a
movwf _x
movlw 0x47
movwf _x+1
; Initialize _y to 0
clrf _y
clrf _y+1
; Initialize _quad to 0
clrf _quad
; Check if the angle is greater than 16383 (90°)
sc_check_greaterthan:
btfss _ang+1, 7
btfss _ang+1, 6
bra sc_check_lessthan
bra sc_adjust_quad2
; Check if the angle is less than -16384 (-90°)
sc_check_lessthan:
btfsc _ang+1, 7
btfsc _ang+1, 6
bra sc_setup_end
; If the angle is in quadrant 3, adjust it to quadrant 4
sc_adjust_quad3:
negf _ang
bc sc_negate_quad3
comf _ang+1
bra sc_adjust_end
; If the low byte negation causes a carry, negate the upper byte
sc_negate_quad3:
negf _ang+1
bra sc_adjust_end
; If the angle is in quadrant 2, adjust it to quadrant 1
sc_adjust_quad2:
comf _ang
comf _ang+1
; Toggle the sign bit and set the '_quad' flag
sc_adjust_end:
btg _ang+1, 7
setf _quad
; Multiply the angle by 2 to get better resolution
sc_setup_end:
bcf STATUS, 0
rlcf _ang
rlcf _ang+1
; Set up the main loop
sc_loop_start:
clrf _i
lfsr FSR0, _atans
; The main loop label
sc_loop:
movff _x, _dy
movff _x+1, _dy+1
movff _i, _j
movf _j
bz sc_bs_x_done
; Loop to shift _dy right
sc_bs_x_loop:
bcf STATUS, 0
rrcf _dy+1
rrcf _dy
btfsc _x+1, 7
bsf _dy+1, 7
decfsz _j
bra sc_bs_x_loop
; Calculate what needs to be added to _x
sc_bs_x_done:
movff _y, _dx
movff _y+1, _dx+1
movff _i, _j
movf _j
bz sc_do_rotation
; Loop to shift _dx right
sc_bs_y_loop:
bcf STATUS, 0
rrcf _dx+1
rrcf _dx
btfsc _y+1, 7
bsf _dx+1, 7
decfsz _j
bra sc_bs_y_loop
; Perform adding operations on _x, _y and _ang
sc_do_rotation:
btfss _ang+1, 7
bra sc_sub_angle
; If _ang is negative
movf POSTINC0, W
addwf _ang
movf POSTINC0, W
addwfc _ang+1
movf _dx, W
addwf _x
movf _dx+1, W
addwfc _x+1
movf _dy, W
subwf _y
movf _dy+1, W
subwfb _y+1
bra sc_loop_bottom
; If _ang is positive
sc_sub_angle:
movf POSTINC0, W
subwf _ang
movf POSTINC0, W
subwfb _ang+1
movf _dx, W
subwf _x
movf _dx+1, W
subwfb _x+1
movf _dy, W
addwf _y
movf _dy+1, W
addwfc _y+1
; Increment the counter and exit the loop if done
sc_loop_bottom:
incf _i
movlw 0x0f
cpfseq _i
bra sc_loop
; Negate _x if it was initially in quadrant 2 or 3
sc_finished:
btfss _quad, 7
bra sc_output
negf _x
bc sc_negate_x
comf _x+1
bra sc_output
; If the low byte negation causes a carry, negate the upper byte
sc_negate_x:
negf _x+1
; Output the calculated _x and _y values
sc_output:
return
; Calculates the magnitude and direction of the given ordered pair
atan2_sqrt:
; Initialize _ang to 0
clrf _ang
clrf _ang+1
; Initialize _quad to 0
clrf _quad
; If the point is in quadrant 2 or 3, make _x positive and set flag
as_check_negative:
btfss _x+1, 7
bra as_shift_x
setf _quad
negf _x
bc as_negate_x
comf _x+1
bra as_shift_x
; If the low byte negation causes a carry, negate the upper byte
as_negate_x:
negf _x+1
; Divide the _x coordinate by 2 to prevent overflowing
as_shift_x:
bcf STATUS, 0
rrcf _x+1
rrcf _x
; Divide the _y coordinate by 2 to prevent overflowing
as_shift_y:
bcf STATUS, 0
rrcf _y+1
rrcf _y
btfsc _y+1, 6
bsf _y+1, 7
; Set up the main loop
as_loop_start:
clrf _i
lfsr FSR0, _atans
; The main loop label
as_loop:
movff _x, _dy
movff _x+1, _dy+1
movff _i, _j
movf _j
bz as_bs_x_done
; Loop to shift _dy right
as_bs_x_loop:
bcf STATUS, 0
rrcf _dy+1
rrcf _dy
btfsc _x+1, 7
bsf _dy+1, 7
decfsz _j
bra as_bs_x_loop
; Calculate what needs to be added to _x
as_bs_x_done:
movff _y, _dx
movff _y+1, _dx+1
movff _i, _j
movf _j
bz as_do_rotation
; Loop to shift _dx right
as_bs_y_loop:
bcf STATUS, 0
rrcf _dx+1
rrcf _dx
btfsc _y+1, 7
bsf _dx+1, 7
decfsz _j
bra as_bs_y_loop
; Perform adding operations on _x, _y and _ang, shifting the _atans right one
as_do_rotation:
MOXVE POSTINC0, PRODL
MOXVE POSTINC0, PRODH
bcf STATUS, 0
rrcf PRODH
rrcf PRODL
btfsc _y+1, 7
bra as_sub_angle
; If _y is positive
movf PRODL, W
addwf _ang
movf PRODH, W
addwfc _ang+1
movf _dx, W
addwf _x
movf _dx+1, W
addwfc _x+1
movf _dy, W
subwf _y
movf _dy+1, W
subwfb _y+1
bra as_loop_bottom
; If _y is negative
as_sub_angle:
movf PRODL, W
subwf _ang
movf PRODH, W
subwfb _ang+1
movf _dx, W
subwf _x
movf _dx+1, W
subwfb _x+1
movf _dy, W
addwf _y
movf _dy+1, W
addwfc _y+1
; Increment the counter and exit the loop if done
as_loop_bottom:
incf _i
movlw 0x0e
cpfseq _i
bra as_loop
; Multiply the _x value by 19898 and divide by 2^14 to scale it
as_scale_x:
movff _x, _dx
movff _x+1, _dx+1
movlw 0xba
mulwf _dx
MOXVE PRODH, _x
movlw 0x4d
mulwf _dx+1
MOXVE PRODH, _dy
MOXVE PRODL, _x+1
movlw 0xba
mulwf _dx+1
movf PRODL, W
addwf _x, F
movf PRODH, W
addwfc _x+1, F
clrf WREG
addwfc _dy, F
movlw 0x4d
mulwf _dx
movf PRODL, W
addwf _x, F
movf PRODH, W
addwfc _x+1, F
clrf WREG
addwfc _dy, F
movlw 0x06
movwf _j
as_scale_bs_loop:
bcf STATUS, 0
rrcf _dy
rrcf _x+1
rrcf _x
decfsz _j
bra as_scale_bs_loop
; Check if the quadrant was originally changed
as_check_quad:
btfss _quad, 7
bra as_output
btfss _ang+1,7
bra as_adjust_quad1
; If the angle is in quadrant 4, adjust it to quadrant 3
as_adjust_quad4:
negf _ang
bc as_negate_quad4
comf _ang+1
bra as_adjust_end
; If the low byte negation causes a carry, negate the upper byte
as_negate_quad4:
negf _ang+1
bra as_adjust_end
; If the angle is in quadrant 1, adjust it to quadrant 2
as_adjust_quad1:
comf _ang
comf _ang+1
; Toggle the sign bit
as_adjust_end:
btg _ang+1, 7
; Output the calculated angle and hypotenuse values
as_output:
return
endasm
OverTrig: 'jump over code
1 Attachment(s)
Re: Erroneous mathematical processing?
Copy and pasted into a new inc file but still no luck. The hyp is different but still wrong.
Attachment 8954
1 Attachment(s)
Re: Erroneous mathematical processing?
BTW: I tried compiling the same code but with the "Use Compiler Long Words" checkbox un-ticked this time in the Microcode Studios settings.
Still erroneous results, but curiously different:
Attachment 8955
Troy
