Does anyone know how to get around rounding problems?

[Egil]

I have started to make routines that hopefully will end in a small satellite tracking program to keep track of satellites with radiobeacons and open transponders in the radioamateur bands. Planning to do the needed routines for a ground track first.
When searching for formulas, I managed to get around all the confusion on internet, mixing the different "types" of miles. I even found online calulators out there stating they operated with kilometers and nautical miles, but inspecting the webpage source revealed that their scripts in fact used Statue Miles only!

Since the only length variable used in most of navigation formulas is the Earth Radius, it is easy to get set this straight. Earth Radius given in kilometers, produce results in kilometers. Earth Radius given in nautical miles to get results in nautical miles, and so on.
But I ran into problems when doing the routines for calulating directions (bearings). I think rounding problems are the main reason for this. And maybe the way IWB calculates the PI constant.

Most textbooks agree on the following constants:
Ratio of a circle's circumference to its diameter: PI = 3.1415926535897932384626433832795 .....
Factor to convert decimal degrees to radians:  0.01745329252
Factor to convert decimal degrees to radians: 57.29577951308

I made some test code to see what I can achieve with a precision of 18 decimals. (precision greater than 17 only produce zeros)
A screendump of the output is posted below. Also posted is the code I used for the test.

Some may ask if such insane precision really is needed. The answer is yes. Since the distances to these satellites are so great, I need great precision to get correct calulations for bearing and elevation angles.

So if anyone have ideas or suggestions on how to get around these problems, I'll be very grateful.

Code Select Expand


'
' trig_test.iwb
'
SETPRECISION 18

OPENCONSOLE
print
print " Degrees to Radians constant:   ",Deg2Rad(1)
print
print " Radians to Degrees constant:  ",Rad2deg(1)
print
print " PI:                            ",4*atan(1)
print:print
print "press any key to end"

WAITCON
CLOSECONSOLE
END

'
SUB Rad2deg(rad:FLOAT),float
'----------------------------------------------------------------------------------------
' Convert Degrees to Radians
' PI = 4 * atan(1)
' Returns Degrees in decimal format
'----------------------------------------------------------------------------------------
RETURN rad*180/(4*atan(1))
ENDSUB

'
SUB Deg2Rad(deg:FLOAT),float
'----------------------------------------------------------------------------------------
' Convert Degrees to Radians
' PI = 4 * atan(1)
' deg = Degrees in decimal format
' Returns Radians
'----------------------------------------------------------------------------------------
RETURN (deg*(4*atan(1)/180))
ENDSUB


[LarryMc]

Egil
I looked at the source code for the atan functions
They are all written in assembly  and look like they are using the math coprocessor in the computer.

All I can suggest is using an approximation routine to calculate what you want.

Hope someone else can give you a better answer.

[ckoehn]

** Ah yes.. I see what you are saying.. it is in the trig function not yielding a proper result.  I'm leaving this here anyway for those who would want to print large precision numbers  **


The print statement is what is rounding it off.  Using SETPRECISION makes no difference above 18.

The maximum value for type float is:  3.40282e+038
The maximum value for type double is:  1.79769e+308
The maximum value for type int is:  2147483647
The maximum value for type short int is:  32767

Try this...

Code Select Expand


'
' trig_test.iwb
'
SETPRECISION 18

OPENCONSOLE
print
print " Degrees to Radians constant:  ",:mPrint(Deg2Rad(1),30)
print
print " Radians to Degrees constant:  ",:mPrint(Rad2deg(1),30)
print
print " PI:                           ",:mPrint(4*atan(1),30)
print:print
print "press any key to end"

WAITCON
CLOSECONSOLE
END


SUB mPrint(num:DOUBLE,prec:BYTE)
Dim i:INT
Dim s:string

s=LTRIM$(str$(int(num)))+"."
num=num-int(num)

for i=1 to prec
num=num*10
s+=LTRIM$(STR$(int(num)))
num=num-int(num)
next i
print s
ENDSUB
'
SUB Rad2deg(rad:DOUBLE),DOUBLE
'----------------------------------------------------------------------------------------
' Convert Degrees to Radians
' PI = 4 * atan(1)
' Returns Degrees in decimal format
'----------------------------------------------------------------------------------------
RETURN rad*180/(4*atan(1))
ENDSUB

'
SUB Deg2Rad(deg:DOUBLE),DOUBLE
'----------------------------------------------------------------------------------------
' Convert Degrees to Radians
' PI = 4 * atan(1)
' deg = Degrees in decimal format
' Returns Radians
'----------------------------------------------------------------------------------------
RETURN (deg*(4*atan(1)/180))
ENDSUB


Later,
Clint

[ckoehn]

LarryMc,

How do you import a library like this one?

http://www.crbond.com/extended_precision.htm

Maybe Egil could use something like this.

Later,
Clint

[Egil]

Thanks to both of you!

This really gave me something to work on. I really appreciate it.
The most important value in the angular calculations is the PI constant, it is used for setting every variable in the formulas. So maybe I shall define it as a constant with as many desimals I can find in a textbook and drop the atan function.
When the calculated values are ready to be printed, I only need one decimal. So I was relieved to know that I can work internally with very high precision.

And Clint, thanks for the URL. Looks promising, so will study that webpage more thoroughly tomorrow. And the mPrint function really opened my eyes.

[LarryMc]

That math package is based on a dll.
To make the necessary lib file all you have to do is create it via Tools > Create Import Library from the IDE's Maiin Menu.
The difficult part will be to create the inc file with the proper declarations.

[Bill-Bo]

LarryMc,

How do you select all the text in a code window like
the one above? You posted it once, but I cannot find
it. I ever tried different searches.

Bill

[LarryMc]

LarryMc,

How do you select all the text in a code window like
the one above? You posted it once, but I cannot find
it. I ever tried different searches.

Bill

put the cursor at the beginning of the 1st line in the code window
hold down the CTRL key
use the mouse to scroll down the  code window
while still holding the CTRL key down click at the end of the last line in the code window
now everything should be selected
now hit CTRL-C to copy

[ckoehn]

It looks like there is too much rounding error in double precision numbers.  Even using a constant doesn't help because I don't think everything is converted to double.

Double-precision examples
3ff0 0000 0000 0000   = 1
3ff0 0000 0000 0001   â‰ˆ 1.0000000000000002, the smallest number > 1
3ff0 0000 0000 0002   â‰ˆ 1.0000000000000004


Later,
Clint

[Bill-Bo]

LarryMc,

I tried your way on reply #2. The code
is not highlighted and the mouse does
not highlight. All that happens is a blue
rectangle outlines the whole reply. I can
only copy the code by using the mouse
only and clicking at the top of the code
window and drag down to the bottom of
the code, then copy.

Bill

[Egil]

It looks like there is too much rounding error in double precision numbers.  Even using a constant doesn't help because I don't think everything is converted to double.

That's what I suspect too. But I have been thinking of the way I used to do satellite calulations way before I had a computer. It was done using a slide rule, paper and pencil. Had  to break the formulas into several parts to ease work though, but the final results were very accurate.
So  tomorrow I'll try to rewrite the code, using that approach. 

[ckoehn]

I'm trying to get a library ported over but I'm having some trouble.  Wish someone would do a tutorial (where is Sapero  ) on the conversion.

LarryMc, can you tell me why I'm getting these errors when it tries to link?

I made a library with IWBasic and it installed it in the lib folder.

Then I modified the inc files that came with it (I have them below with the lib) and put $USE "NUMS.lib" at the top of the "NUMS.inc" file.

What am I doing wrong?  Does it also need the DLL somewhere?  When I tried to regsvr32 the dll it said it loaded but had no entry point.


Compiling...
TrigTest.iwb
No Errors

Linking...
IWBasic Linker v1.11 Copyright © 2011 Ionic Wind Software
Unresolved external __imp_nu_AllocVar
Error: C:\Users\Public\Documents\IWBasic\projects\My Projects\TrigTest\TrigTest.o - Unresolved extern __imp_nu_AllocVar
Error: C:\Users\Public\Documents\IWBasic\projects\My Projects\TrigTest\TrigTest.o - Unresolved extern __imp_nu_PutString_DS
Error: C:\Users\Public\Documents\IWBasic\projects\My Projects\TrigTest\TrigTest.o - Unresolved extern __imp_nu_ATan
Error: C:\Users\Public\Documents\IWBasic\projects\My Projects\TrigTest\TrigTest.o - Unresolved extern __imp_nu_Mul
Error: C:\Users\Public\Documents\IWBasic\projects\My Projects\TrigTest\TrigTest.o - Unresolved extern __imp_nu_GetAny_DS
Error: C:\Users\Public\Documents\IWBasic\projects\My Projects\TrigTest\TrigTest.o - Unresolved extern __imp_nu_UnAllocVar
Error(s) in linking C:\Users\Public\Documents\IWBasic\projects\My Projects\TrigTest\TrigTest.exe

Later,
Clint

[ckoehn]

OK..  found one problem.  All the declares needed to be uppercase.. moving on... moving on...

Now I need to fix some of the way the functions are called...

Later,
Clint

[ckoehn]

OK,  try this.  It does have a nag screen but I wanted to try to see if I could convert a dll to a lib.

The bottom PI: is the one that is computed by the precision math library.  Everything above it is going to be inaccurate.  Look at the code.. you'll know what I mean. 

LarryMc,  how do you pass a string byval instead of byref on a declare import,?  I have a few of those in there so some of the strings functions don't work.

Later,
Clint

[LarryMc]

... how do you pass a string byval instead of byref on a declare import,?  

When I've ran into that in the past I've had to declare the parameter as a UINT and then do the string like this

Code Select Expand

declare import mysub(mystring:uint,blah...)
string mystring ="aaaaaa"
mysub(&mystring, blah......)
maybe that will work this time

[ckoehn]

Thanks LarryMc,  will try that tomorrow.

Later,
Clint

[Egil]

What I see when starting the exe file looks very good, Clint.

Talked to my grandson last night. He meant that if I could use a atan2 function (taking two arguments), and using doubles, most of these problems will dissapear.
When inspecting your source code I find this statement: nu_ATan(a, b). Looks like that function works similar to the atan2 functions used in C/C++.
Have to study the NUMS documentation a little more...

[ckoehn]

Any number stored in double precision will be inaccurate, no matter how it gets there.

I have found an issue with this library with some of the functions not working, like nu_GetAny_Double(a).  It will not place the return value in a double like it is suppose to for some reason.  I am positive the declare is correct.

I haven't made any modifications to the string function calls yet.  I need to pass a literal string like so.. nu_PutString_DS("1",a) .. so I'm not sure using a pointer will work like you pointed out LarryMc.  Using "&" is basically sending a pointer to a location.

I apologize LarryMc.  Your way and Pauls way (below) will work I think. It is a round about way of doing which is necessary because of the byval.  

LarryMc, I tried it your and Pauls way and it didn't work.. 

Egil,  check the code below..

Code Select Expand

'
' trig_test.iwb
'
$INCLUDE "windowssdk.inc"

$INCLUDE "NUMS.inc"

LONG npi
LONG a
LONG b
LONG c
LONG i

double pi = 3.1415926535897932384626433832795
double ang, mult
ang=1
mult=4

i = nu_SetPrecision(nu_SetPrec_Digits, 35) '<-- set precision

npi = nu_AllocVar(nu_DT_DecFloat)
a = nu_AllocVar(nu_DT_DecFloat)
b = nu_AllocVar(nu_DT_DecFloat)
c = nu_AllocVar(nu_DT_DecFloat)

nu_PutDouble(ang, a)

'nu_Put_String_DS("1",a)

nu_ATan(a, b)
nu_PutDouble(mult, a)
nu_Mul(a, b, npi)

SETPRECISION 30

OPENCONSOLE
CLS 'clears the nag screen
print
print " Deg to Rad         accurate:  ",:mPrint(Deg2Rad(1),40)
print
print " Rad to Deg         accurate:  ",:mPrint(Rad2deg(1),40)
print
print " PI: accurate                 ",Nu_GetAny_DS(npi)
print " PI: not accurate              ",:mPrint(pi,40)
print
print "press any key to end"

WAITCON

CLOSECONSOLE

nu_UnAllocVar(npi)
nu_UnAllocVar(a)
nu_UnAllocVar(b)
nu_UnAllocVar(c)

END

SUB mPrint(num:DOUBLE,prec:BYTE)
Dim i:INT
Dim s:string

s=LTRIM$(str$(int(num)))+"."
num=num-int(num)

for i=1 to prec
num=num*10
s+=LTRIM$(STR$(int(num)))
num=num-int(num)
next i
print s
ENDSUB
'
SUB Rad2deg(rad:DOUBLE),DOUBLE
'----------------------------------------------------------------------------------------
' Convert Degrees to Radians
' PI = 4 * atan(1)
' Returns Degrees in decimal format
'----------------------------------------------------------------------------------------
double nrad = 180

nu_PutDouble(nrad, a) 'store 180 in a
nu_PutDouble(rad, b) 'store radian in b
nu_Mul(a, b, c) 'store rad * 180  in c
nu_Div(c, npi, a) 'store c/pi in a

print nu_GetAny_DS(a) '<--- this returns a good number as a string
print "                 not accurate:  ",

nrad = VAL(nu_GetAny_DS(a)) '<--- this one is inaccurate because of double

RETURN nrad
ENDSUB

'
SUB Deg2Rad(deg:DOUBLE),DOUBLE
'----------------------------------------------------------------------------------------
' Convert Degrees to Radians
' PI = 4 * atan(1)
' deg = Degrees in decimal format
' Returns Radians
'----------------------------------------------------------------------------------------
double nrad = 180

nu_PutDouble(nrad, a) 'store 180 in a
nu_Div(npi, a, c) 'store pi/180 in c
nu_PutDouble(deg, b) 'store degree in b
nu_Mul(b, c, a) 'store deg * c in a

print nu_GetAny_DS(a) '<--- this returns a good number as a string
print "                not accurate:   ",

nrad = VAL(nu_GetAny_DS(a)) '<--- this one is inaccurate because of double

RETURN nrad
ENDSUB


Later,
Clint

[ckoehn]

LarryMc,

There was one thing I found on the forum but it was reguarding EBasic.  I think it was Paul that said to do it this way..

TYPE BVSTRING
   STRING s
ENDTYPE

And then in the call..   declare import, calltomake(sstr as BVSTRING BYVAL).

That did not work for IWBasic.  It gave a "cannot convert string to structure" error.

Still digging..

Later,
Clint

[LarryMc]

Clint
change this

Code Select Expand

declare import, calltomake(sstr as BVSTRING BYVAL)
to this

Code Select Expand

declare import, calltomake(sstr as BVSTRING)
which is what was done in the post that covered using BVSTRING

[ckoehn]

Nothing seems to work in trying to get a string to pass by value.

I did find another function which passes by ref... nu_PutString_PAZ("1",a) ' a is an initialized precision variable

You can also pass a double without declaring a variable first.. nu_PutDouble(1,a) ' a is an initialized precision variable

I guess I need to let this go for awhile and get back to work..

Later,
Clint

[ckoehn]

I'm working on importing the GMP.dll when I'm not busy and should have the necessary function done toward the end of the week.

Having some difficulty with the GMPR functions included in the dll.  The trig functions seem to work, but the gmpr_mul doesn't.

Copied the declares from a PowerBASIC forum and modified them.

Later,
Clint

[ckoehn]

I never did get the mpfr_mul to work.  You don't really need it.  I just transfered the results to mpf_mul and it works fine.

Other functions need to be "declared" yet but I currently don't have time to do that.  If one is needed it is easy to add.

Copy the "gmp.dll" to either your system32 folder or SysWow64 folder.
Copy the "gmp.inc" to iwbdev include folder
Copy the "gmp.lib" to iwbdev libs folder

This requires "windowssdk.inc" which the "gmp.inc" has a "$INCLUDE "windowssdk.inc" in it.

Code Select Expand

'
' trig_test.iwb
'
$INCLUDE "gmp.inc"

MPF npi
MPFR a_tan
MPFR b
MPF x
MPF y
MPF z
IString outb[1000]   'buffer for output

MPFRinit2(a_tan, 150) : MPFRinit2(b, 150)
MPFinit2(npi, 150) : MPFinit2(x, 150) : MPFinit2(y, 150) : MPFinit2(z, 150)

MPFRsetstr(b, "1", 10, 0) 'place 1 in b
MPFRatan(a_tan,b,0) 'atan(1)
MPAsprintf(outb, "%0.60Fe", a_tan) 'store in outb to transfer to MPF since MPFR doesn't work

MPFsetstr(x, outb, 10) 'x now hold atan(1)
MPFsetstr(y, "4", 10) 'y hold 4
MPFmul(npi, x, y) 'mul atan(1)*4

OPENCONSOLE
print
print " Deg to Rad :  ",:Print Deg2Rad(1)
print
print " Rad to Deg :  ",:Print Rad2deg(1)
print
MPAsprintf(outb, "%0.50Fe", npi)
print "          PI:  ",:Print outb
print
print "press any key to end"

WAITCON

CLOSECONSOLE

MPFRclear(a_tan) : MPFRclear(b)
MPFclear(npi) : MPFclear(x) : MPFclear(y) : MPFclear(z)
END


SUB Rad2deg(rad:DOUBLE),STRING
'----------------------------------------------------------------------------------------
' Convert Degrees to Radians
' PI = 4 * atan(1)
' Returns Degrees in decimal format
'----------------------------------------------------------------------------------------
MPFsetstr(x, "180", 10) 'store 180 in x
MPFsetstr(y, LTRIM$(STR$(rad)), 10) 'store radian in y
MPFmul(z, x, y) 'store rad * 180  in z
MPFdiv(x, z, npi) 'store z/pi in x
MPAsprintf(outb, "%0.50Fe", x)
RETURN outb
ENDSUB

'
SUB Deg2Rad(deg:DOUBLE),STRING
'----------------------------------------------------------------------------------------
' Convert Degrees to Radians
' PI = 4 * atan(1)
' deg = Degrees in decimal format
' Returns Radians
'----------------------------------------------------------------------------------------
MPFsetstr(x, "180", 0) 'store 180 in x
MPFdiv(z, npi, x) 'store pi/180 in z
MPFsetstr(y, LTRIM$(STR$(deg)), 0) 'store degree in y
MPFmul(x, y, z) 'store z * y in x
MPAsprintf(outb, "%0.50Fe", x)
RETURN outb
ENDSUB

Later,
Clint

[Egil]

I really appreciate your work Clint!
Will be gone away for the rest of the week, but will bring my laptop, so hope to be able to try out your work. Expect to return monday evening.
Going to have some surgery on my spine, the last of three, and hopefully it will help me get rid of most of the painkillers I have to use to be able to function. (I call them "brainkillers"...)
Hopefully I'll be able to do some coding when staying there. Do not know yet if they'll put me in a room with Internet access. If not I'll be silent for almost a week...

[LarryMc]

...Going to have some surgery on my spine...

Wish you the best and a speedy recovery.