preserving as much precision as possible

tosjduenfs · 2013-10-16 06:03

I've been working on code that converts a stellar objects RA/DEC coordinates to ALT/AZ coordinates. Basically I first calculate the number of days it has been since Jan 1, 2000 (J2000) and then through a series of calculations I can determine and track where an object should be in the sky given my geographical coordinates and my local time.

In my first attempt I calculated the time since J2000 in decimal days using the float32 object. Unfortunately due to the size of the numbers involved and the lack of precision of a 32 bit floating point number my program would only generate new ALT/AZ coordinates every ~50 seconds or so (when enough time had passed that the last digit of JD changed) rather than every second as the program is written.

So the plan now is to convert the code to use integer seconds rather than decimal days. I've successfully changed the time calculation to accurately give seconds as an integer rather than decimal days but now I am not sure on how to continue using integer math. The next step is to use the following formula to find Greenwich Mean Sidereal Time:

GMST = (280.46061837 + 360.9856473662862*D) mod 360 ' where D is days since J2000 and GMST is in degrees (as of writing this post D = 5036.795139 days)

This formula is accurate to 0.1 seconds over a century but all I really need is for it to be accurate to a few minutes of actual GMST and to have a precision of 1 second so that I can generate new ALT/AZ coordinates every second for tracking.

If I were still using floating point math I would have no problem converting the formula but I'm not exactly sure how to do it with integer math while preserving as much precision as possible.

Any tips would be appreciated. I can post my code later tonight when I get home if that would be more helpful.

Mike

Bean · 2013-10-16 07:30

I depends on what resolution you need to have.

If GMST is in degrees, I assume you would want better resolution than that. If you want seconds(nav) then multiply everything by 3600.

GMST(nav sec) = (1009658.226132 + 1299476.33051863032 * D) MOD 1296000

Since you are using Julian seconds instead of days you would use:

GMST(nav sec) = (1009658.226132 + 15.040235306928591667 * S) MOD 1296000

That is off the top of my head. Does that final formula work ? If it does I can work out some integer math way to compute it.

Bean

Mark_T · 2013-10-16 08:09

tosjduenfs wrote: »

GMST = (280.46061837 + 360.9856473662862*D) mod 360 ' where D is days since J2000 and GMST is in degrees (as of writing this post D = 5036.795139 days)

Well I'd first convert the formula to use seconds rather than days - ie divide the 360.98... constant by 86400, replace D by S...

Then I'd select, as has been mentioned, a granularity for the GMST, and do that scaling for both constants (x 3600 for arc-seconds).

You then have a formula to calculate mod 1296000 (arc-seconds per revolution). The constant used to multiply S needs to
be approximated by a integer ratio (with at least one of numerator or denominator being full 32 bits for max precision).

So you'll need a 32x32->64 bit multiply and a 64bit modular reduction to compute it.

Alternatively you could choose your granularity to be 2^-32 revolutions, then the mod is for free, but you'd have to
scale the result back to arc-seconds or whatever.

Bean · 2013-10-16 08:26

Mark_T wrote: »

Alternatively you could choose your granularity to be 2^-32 revolutions, then the mod is for free, but you'd have to
scale the result back to arc-seconds or whatever.

That is a great idea. That is the path I would follow.

Bean

tosjduenfs · 2013-10-16 09:17

Bean wrote: »

I depends on what resolution you need to have.

If GMST is in degrees, I assume you would want better resolution than that. If you want seconds(nav) then multiply everything by 3600.

GMST(nav sec) = (1009658.226132 + 1299476.33051863032 * D) MOD 1296000

Since you are using Julian seconds instead of days you would use:

GMST(nav sec) = (1009658.226132 + 15.040235306928591667 * S) MOD 1296000

That is off the top of my head. Does that final formula work ? If it does I can work out some integer math way to compute it.

Bean

This makes sense to me, it hadn't occurred to me to use arc seconds for gmst. Arc seconds would be plenty accurate for my purpose.

tosjduenfs · 2013-10-16 09:47

Mark_T wrote: »

Well I'd first convert the formula to use seconds rather than days - ie divide the 360.98... constant by 86400, replace D by S...

Then I'd select, as has been mentioned, a granularity for the GMST, and do that scaling for both constants (x 3600 for arc-seconds).

You then have a formula to calculate mod 1296000 (arc-seconds per revolution). The constant used to multiply S needs to
be approximated by a integer ratio (with at least one of numerator or denominator being full 32 bits for max precision).

So you'll need a 32x32->64 bit multiply and a 64bit modular reduction to compute it.

Alternatively you could choose your granularity to be 2^-32 revolutions, then the mod is for free, but you'd have to
scale the result back to arc-seconds or whatever.

I'm trying to wrap my head around this one. By granularity I'm assuming you mean that 1/(2^32) is the smallest increment that GMST could vary and then the mod would be"free" because 2^32 is the maximum size of a 32 bit integer?

Also I suppose it might be helpful to explain where it goes after this.

After getting GMST I calculate the Local Hour Angle:

LHA = GMST + Longitude - RA ' where LHA is degrees, Longitude is your geographical longitude also in degrees and RA is the right ascension of the object you are trying to locate also in degrees

Then you can calculate ALT/AZ

sin(ALT) = sin(DEC)*sin(LAT)+cos(DEC)*cos(LAT)*cos(LHA)

cos(Az) = [sin(DEC) - sin(ALT)*sin(LAT)] / [cos(ALT)*cos(LAT)] ' where Dec is the declination of the object you are locating, ALT is the altitude calculated above, and LAT is your geographical latitude

All of the inputs to the above equations had to be converted to radians first. Obviously Arcsin and ArcCos are applied to get the final results which are then converted back to degrees Which will be fed into a PI motor controller. Ultimately they will control two gear motors with encoders which will generate 67200 pulses per revolution of the telescope in ALT and AZ. Would it be best to work backwards from 67200?

Tracy Allen · 2013-10-16 11:42

I'm curious, what does GMST = 0 represent? Is that the transit of some particular star (which?) above Greenwich latitude 0?

The precision of the multiplier is 13 decimal digits, around 44 bits binary. If you use lower precision, error will accumulate over time and the granularity will increase. But at 31 bits, it should still be near 0.01 second per year, and at the 24 bits (precision of the floating point significand) it should be no more than 2 seconds per year. In light of that, I don't understand your statement, "Unfortunately due to the size of the numbers involved and the lack of precision of a 32 bit floating point number my program would only generate new ALT/AZ coordinates every ~50 seconds or so (when enough time had passed that the last digit of JD changed) rather than every second as the program is written. "

For better precision, the object F32.spin uses CORDIC for the trig functions out to full 24 bits. Float32 on the other hand uses linear interpolation to 16 bits in the 11 bit HUB table.

tosjduenfs · 2013-10-16 13:19

Tracy Allen wrote: »

I'm curious, what does GMST = 0 represent? Is that the transit of some particular star (which?) above Greenwich latitude 0?

The precision of the multiplier is 13 decimal digits, around 44 bits binary. If you use lower precision, error will accumulate over time and the granularity will increase. But at 31 bits, it should still be near 0.01 second per year, and at the 24 bits (precision of the floating point significand) it should be no more than 2 seconds per year. In light of that, I don't understand your statement, "Unfortunately due to the size of the numbers involved and the lack of precision of a 32 bit floating point number my program would only generate new ALT/AZ coordinates every ~50 seconds or so (when enough time had passed that the last digit of JD changed) rather than every second as the program is written. "

For better precision, the object F32.spin uses CORDIC for the trig functions out to full 24 bits. Float32 on the other hand uses linear interpolation to 16 bits in the 11 bit HUB table.

I will post my original code when I get home tonight so you can see the behavior I am talking about. I am using F32 now but I'm not sure If i was using it when I first had this problem so I will test that. I think GMST = 0 would mark the vernal equinox, 280.46... is the initial angle at j2000 and 360.98... is the angle the earth rotates in 24 hours.

What I was trying to say before is that my original code takes its input from the realtimeclock object and recalculates the Julian Day Number every second and then it calculates GMST then LHA then finds the ALT/AZ coordinates of the object I'm trying to track. For example, right now it has been about 5000 days since J2000 so with float32 I can get two or three decimal places (so it seams) after 5000 so 5000.000. The best precision I can get is a thousandth of a day, and with 86400 seconds in a day you can see that when I add 1 second to 5000.000 days it will get rounded to 5000.000. Only when enough seconds have passed will it get up to 5000.001 and thus change the rest of the calculations down the line. Combine that with other rounding errors throughout the code and I only get updated alt/az coordinates about every~50seconds even though they are recalculated, with updated time, every second.

If I could get 13 decimal digits of precision I think my code would run fine but that is not what I see when I send my variables to PST. I'll post the code when I get home tonight so you can see what I mean, maybe it will be a simple fix.

tosjduenfs · 2013-10-16 17:22

Here is my code, the ALT and AZ values only vary from the actual values by about 1/10 of a degree which is fine for my application since I will do an initial alignment of the telescope but as you can see they only change about every 50 seconds rather than every second like I need them to. In order for this to work you'll need the float to string object and the F32 object. Just load to ram and open the serial terminal.

CON
   _clkmode = xtal1 + pll16x 
  _xinfreq = 5_000_000


VAR
  long Year          'current year
  long Month         'current month
  long Day           'current day
  long Hour          'current hour
  long Minute        'current minute
  long Sec           'current second
  long DeltaT        'time diff from GMT in hours + or -
  long JD            'Julian date in decimal days
  long RaHour        'Right Ascension Hour of target
  long RaMin         'Right Ascension Minute
  long RaSec         'Right Ascension Second
  long RaDecimal     'Right Ascension decimal degrees
  long DecDeg        'Declination Degrees
  long DecMin        'Declination Minutes
  long DecSec        'Declination Seconds
  long DecDecimal    'Declination in decimal degrees
  long Alt
  long Az
  long longdeg       'longitude deg minutes second
  long longmin
  long longsec
  long longdecimal   'longitude in decimal degrees
  long latdeg        'latitude deg minutes seconds
  long latmin
  long latsec
  long latdecimal    'latitude in decimal degrees
  long GMST          'Greenwich Mean Sidereal Time
  long LHA           'Local Hour Angle
  long A             'A B C E F used in calculating Julian Date
  long B
  long C
  long E
  long F
  long D              'used in calc for GMST




OBJ


  RTC : "RealTimeClock"
  FLT : "F32"
  FTS : "FloatString"
  PST : "Parallax Serial Terminal"


  
PUB Converter


  RTC.start
  FLT.start
  PST.start(115200)


  RTC.WriteTimeReg(0,30)
  RTC.WriteTimeReg(1,52)
  RTC.WriteTimeReg(2,2)
  RTC.WriteTimeReg(3,13)
  RTC.WriteTimeReg(4,10)
  RTC.WriteTimeReg(5,13)


repeat


  sec := FLT.Ffloat(RTC.ReadTimeReg(0))
  minute := FLT.Ffloat(RTC.ReadTimeReg(1))
  hour := FLT.Ffloat(RTC.ReadTimeReg(2))
  day := FLT.Ffloat(RTC.ReadTimeReg(3))
  month := FLT.Ffloat(RTC.ReadTimeReg(4))
  year := FLT.Fadd(FLT.Ffloat(RTC.ReadTimeReg(5)),2000.0)
  
  'year := 2013.0
  'month := 10.0
  'day := 13.0
  'hour := 2.0
  'minute := 52.0
  'sec := 30.0
  DeltaT := -5.0


  RaHour := 5.0
  RaMin := 34.0
  RaSec := 30.0


  DecDeg := 22.0
  DecMin := 1.0
  DecSec := 5.0


  Latdeg := 41.0
  Latmin := 15.0
  Latsec := 0.46


  Longdeg := -87.0
  Longmin := 50.0
  Longsec := 17.81
  
'calculate decimal form of Right Ascension  
  RaDecimal := FLT.Fmul(FLT.Fadd(FLT.Fadd(RaHour,FLT.Fdiv(RaMin,60.0)),FLT.Fdiv(RaSec,3600.0)),15.0)


'calculate decimal form of Declination
  DecDecimal := FLT.Fadd(DecDeg,FLT.Fadd(FLT.Fdiv(DecMin,60.0),FLT.Fdiv(DecSec,3600.0)))


'calculate decimal form of latitude
  LatDecimal := FLT.Fadd(LatDeg,FLT.Fadd(FLT.Fdiv(LatMin,60.0),FLT.Fdiv(LatSec,3600.0)))


'calculate decimal form of longitude
  IF LongDeg < 0
    LongDeg := FLT.Fneg(LongDeg)
    LongDecimal := FLT.Fadd(LongDeg,FLT.Fadd(FLT.Fdiv(LongMin,60.0),FLT.Fdiv(LongSec,3600.0)))
    LongDecimal := FLT.Fneg(LongDecimal)
  Else
    LongDecimal := FLT.Fadd(LongDeg,FLT.Fadd(FLT.Fdiv(LongMin,60.0),FLT.Fdiv(LongSec,3600.0)))


'calculate Julian Date in decimal days since Noon Jan 1,2000, fractional parts of all multiplications and divisions are dropped in the following calculations


  
  IF (Month == 1.0) OR (Month == 2.0)
    Year := FLT.FSub(Year,1.0)
  ELSE 


  A := FLT.Floor(FLT.Fdiv(Year,100.0))


  B := FLT.Floor(FLT.Fdiv(A,4.0))


  C := FLT.Floor(FLT.Fadd(FLT.Fsub(2.0,A),B))
  
  E := FLT.Floor(FLT.Fmul(365.25,FLT.Fadd(Year,4716.0)))


  F := FLT.Floor(FLT.Fmul(30.6001,FLT.Fadd(Month,1.0)))


'Stop dropping fractions  (daylight savings time still needs to be accounted for)


  JD := FLT.Fsub(FLT.Fsub(FLT.Fadd(FLT.Fadd(FLT.Fadd(C,Day),E),F),1524.5),2451545.0)
   
'Calculate Greenwich Mean Sidereal Time (GMST) GMST = 280.46061837 + 360.98564736629 * JD, D is JD plus fraction of day from hour min sec and time diff


  D := FLT.Fsub(FLT.Fadd(FLT.Fadd(FLT.Fadd(JD,FLT.Fdiv(hour,24.0)),FLT.Fdiv(minute,1440.0)),FLT.Fdiv(sec,86400.0)),FLT.Fdiv(DeltaT,24.0))


  GMST := FLT.Fmod(FLT.Fadd(280.46061837,FLT.Fmul(360.98564736629,D)),360.0)


'Calculate Local Hour Angle LHA = GMST + Longitude - RA


  LHA := FLT.Fsub(FLT.Fadd(GMST,LongDecimal),RaDecimal)


   IF LHA < 0.0
     repeat until (LHA > 0.0) OR (LHA == 0.0)
       LHA := FLT.Fadd(LHA,360.0)
       
   IF LHA > 360.0
     repeat until (LHA < 360.0) OR (LHA == 360.0)
       LHA := FLT.Fsub(LHA,360.0)


'Calculate Altitude and Azimuth values


  ALT := FLT.Degrees(FLT.Asin(FLT.Fadd(FLT.Fmul(FLT.Sin(FLT.Radians(LatDecimal)),FLT.Sin(FLT.Radians(DecDecimal))),FLT.Fmul(FLT.Fmul(FLT.Cos(FLT.Radians(LatDecimal)),FLT.Cos(FLT.Radians(DecDecimal))),FLT.Cos(FLT.Radians(LHA))))))
 
  AZ := FLT.Degrees(FLT.Acos(FLT.Fdiv((FLT.Fsub(FLT.Sin(FLT.Radians(DecDecimal)),FLT.Fmul(FLT.Sin(FLT.Radians(LatDecimal)),FLT.Sin(FLT.Radians(ALT))))),(FLT.Fmul(FLT.Cos(FLT.Radians(LatDecimal)),FLT.Cos(FLT.Radians(ALT)))))))  
   
    IF FLT.Sin(FLT.radians(LHA)) > 0.0
      AZ := FLT.Fsub(360.0,AZ)


  
  'Repeat
   pst.str(string("sec = "))
   pst.str(FTS.FloatToString(sec))
   pst.newline
   pst.str(string("minute = ")) 
   pst.str(FTS.FloatToString(minute))
   pst.newline
   pst.str(string("hour = ")) 
   pst.str(FTS.FloatToString(hour))
   pst.newline
   pst.str(string("day = ")) 
   pst.str(FTS.FloatToString(day))
   pst.newline
   pst.str(string("month = ")) 
   pst.str(FTS.FloatToString(month))
   pst.newline
   pst.str(string("year = ")) 
   pst.str(FTS.FloatToString(year))
   pst.newline
   pst.str(string("Julian Day = ")) 
   pst.str(FTS.FloatToString(D))
   pst.Newline
   pst.str(string("GMST = ")) 
   pst.str(FTS.FloatToString(GMST))
   pst.Newline
   pst.str(string("LHA = ")) 
   pst.str(FTS.FloatToString(LHA))
   pst.Newline
   pst.str(string("AZ = "))    
   pst.str(FTS.FloatToString(AZ))
   pst.newline
   pst.str(string("ALT = ")) 
   pst.str(FTS.FloatToString(ALT))
   waitcnt(clkfreq/1 + cnt)
   pst.clear

RA-DEC to ALT-AZ Converter 201310131437.spin

davekorpi · 2013-10-16 22:22

I want to do FP math too and wonder why you can not just use Floats? They are 32 bit in the shared code by Cam Thompson...
Float32.spin...

Here is the header...

''=============================================================================
'' Copyright (C) 2006-2007 Parallax, Inc. All rights reserved
''
'' @file Float32
'' @target Propeller
''
'' IEEE 754 compliant 32-bit floating point math routines.
'' This object supports the basic set of floating point routines using one cog.
''
'' ──Float32
''
'' @author Cam Thompson, Micromega Corporation
'' @version
'' V1.2 - March 26, 2007
'' - fixed Pow to handle negative base values
'' V1.1 - Oct 5, 2006
'' - corrected constant value for Radians and Degrees
'' V1.0 - May 17, 2006
'' - original version
''=============================================================================

I want super precision math too in a fifth order polynomial I am doing and just have to wrestle with the MAXIMUM number the float can return...

In the end can you use floats here??

Lawson · 2013-10-17 00:44

tosjduenfs wrote: »

I've been working on code that converts a stellar objects RA/DEC coordinates to ALT/AZ coordinates. Basically I first calculate the number of days it has been since Jan 1, 2000 (J2000) and then through a series of calculations I can determine and track where an object should be in the sky given my geographical coordinates and my local time.

In my first attempt I calculated the time since J2000 in decimal days using the float32 object. Unfortunately due to the size of the numbers involved and the lack of precision of a 32 bit floating point number my program would only generate new ALT/AZ coordinates every ~50 seconds or so (when enough time had passed that the last digit of JD changed) rather than every second as the program is written.

So the plan now is to convert the code to use integer seconds rather than decimal days. I've successfully changed the time calculation to accurately give seconds as an integer rather than decimal days but now I am not sure on how to continue using integer math. The next step is to use the following formula to find Greenwich Mean Sidereal Time:

GMST = (280.46061837 + 360.9856473662862*D) mod 360 ' where D is days since J2000 and GMST is in degrees (as of writing this post D = 5036.795139 days)

This formula is accurate to 0.1 seconds over a century but all I really need is for it to be accurate to a few minutes of actual GMST and to have a precision of 1 second so that I can generate new ALT/AZ coordinates every second for tracking.

If I were still using floating point math I would have no problem converting the formula but I'm not exactly sure how to do it with integer math while preserving as much precision as possible.

Any tips would be appreciated. I can post my code later tonight when I get home if that would be more helpful.

Mike

When doing floating point addition it's easy to loose bits of precision if the inputs have different magnitudes. The given equation will suffer badly from this. (D*360.xxx >> 280.xxx) Also, Fmod( a, b) will have reduced precision when 'a' is much larger than 'b'. I'd try moving the modulus function inside the equation. Looks like this will require re-figuring the three constants of the equation though.

If you don't care much about speed, I've written a spin only floating point library with all the same functions as F32. It's a bit more precise than F32 with Log and Exp, but isn't standards compliant with all the edge cases. (see the links to FME.spin in my signature)

After making FME I've wondered weather it would be worth modifying it so most of the functions worked on unpacked floats pulled off of a stack? You could get ~30 bits of mantissa, expanded exponent range, plus a nice speed boost. The main cost is that floats would take 3-longs during equation evaluation. (and creating full precision constants would need some extra code)

Marty

Bean · 2013-10-17 04:48

I haven't looked at the code you posted, but here is what I came up with for the original post:

CON
  _clkmode = xtal1 + pll16x
  _xinfreq = 5_000_000

VAR

OBJ
  Terminal   : "Parallax Serial Terminal"  

PUB Main
  Terminal.Start(115200)
  
  waitcnt(cnt+clkfreq*5) ' Wait for serial terminal
  Calc(435_179_100) ' Julian Seconds 5036.795139 * 86400
  Calc(435_179_101)


PUB Calc(Given) | jSec, sum, temp
  jSec := Given
  Terminal.Char(13)
  Terminal.Dec(jSec)
  Terminal.Str(String(" J-Seconds", 13))

  ' sum = jSec * 49846.3718463343  ' ((360.9856473662 / 360) / 86400) * 2^32 
  sum :=  jSec * 49846
  sum +=  jSec ** 1597067845

  ' sum = sum + 3346025510 ' (280.4606184 / 360) * 2^32
  sum += 3346025510

  Terminal.Dec(sum)
  Terminal.Str(String(" MOD 2^32", 13))
                              
  ' Convert from MOD 2^32 to MOD 1296000  ' 360 * 60 * 60
  temp := sum ** 1296000
  if temp < 0
    temp += 1296000
  Terminal.Dec(temp) ' Display arc-seconds
  Terminal.Str(String(" Arc-seconds", 13))

  ' Convert from MOD 2^32 to MOD 67200
  temp := sum ** 67200
  if temp < 0
    temp += 67200
  Terminal.Dec(temp) ' Display encoder pulses
  Terminal.Str(String(" Encoder pulses", 13))

And here is the output:

435179100 J-Seconds
1565451767 MOD 2^32
472372 Arc-seconds
24493 Encoder pulses

435179101 J-Seconds
1565501613 MOD 2^32
472387 Arc-seconds
24494 Encoder pulses

So 1 Julian second change resulted in a 15 arc-second change.

Bean

tosjduenfs · 2013-10-17 06:24

Thank you for taking the time to write that out Bean, I'll have to do some studying of your code. I think I follow the conversions you make but there are some things I don't quite understand. For instance:

sum := jSec * 49846
sum += jSec ** 1597067845

so sum = (jsec * 49846) + (jsec ** 1597067845)

I'm not sure I understand the ** operation. It returns the upper 32 bits of a 32x32bit multiplication. So the upper 32 bits is the part that would normally be truncated? As in the upper 32 bits is the precision that is normally lost after the decimal? Also where does the 1597067845 come from and what exactly do these two lines accomplish?

I wish it were just as easy as going from GMST to encoder pulses but there are also trig functions these number must pass through first.

Again, thanks for your time.

Bean · 2013-10-17 06:41

The multiplication of jSec * 49846.3718463343 needs to be done in two parts.
First the whole value is multiplied.

sum := jSec * 49846

Then the fractional part. 0.3718463343
To multiply by a faction you take the fractional part and multiply it by 2^32 then use the ** operator.
So 0.3718463343 * 2^32 = 1597067845

When you multiply two 32 bit numbers you get a 64 bit result. The * operator gives the lower 32 bits, the ** gives the upper 32 bits (or the overflow).

The 49846.3718463343 comes from converting from 360 modulus to 2^32 modulus (and also going from days to seconds as the units).

So 360.9856473662 / 86400 because we are using seconds instead of days as the units.
Then (360.9856473662 / 86400) / 360 to use a modulus of 1 instead of 360
Then ((360.9856473662 / 86400) / 360) * 2^32 to use a modulus of 2^32 instead of 1
And that gives us the 49846.3718463343

Because the modulus is 2^32 we don't have to even think about it because all math is performed mod 2^32 by the propeller.

As far a ** math goes, you can think of 2^32 as 1.0000000000, so to get a different modulus you just high-multiply by the modulus you desire.

I hope that helps,

Bean

Tracy Allen · 2013-10-17 10:18

I'm all in favor of integer math where possible. To amplify Bean's explanation of ** as a fraction...

Think of it as multiplying two numbers together, followed by shifting right 32 bits. That is the same as dropping the low 32 bits from the result of a 32x32 mulitply.
Z := X ** F
is mathmatically equivalent to
' Z = X * Y = X * F / 2³²
where Y<1 and the factor F = Y * 2³²
The fraction Y has been well approximated by the fraction F/2³²

If Y has both an integer part and a fractional part, they have to be done separately,
Z := X * (integer part of Y) + X ** (fractional part of Y * 2³²)

There is caveat that involves negative numbers. To use ** in this way, both factors X and F have to be positive numbers, (less than or equal to posx:2147483647) or else the Prop's twos complement math returns the wrong answer when interpreted as fractional multiply. A trick for fractions greater than or equal to 0.5 is to separate it out like this:
Z = X * (integer part of Y) + 2*X ** (fractional part of Y/2 * 2³²)

tosjduenfs · 2013-10-17 10:25

Thank you Bean and Tracy I understand now, that was very helpful. I'll try putting this in my code when I get home tonight. So now I have GMST in arcseconds rather than degrees. Is there any way to use arc seconds in a trigonometric function rather than converting it back to decimal degrees or radians? For instance the next calculation is:

ALT = arcsin[ sin(Declination) * sin(Latitude) + cos(Declination) * cos(Latitude) * cos(GMST + Longitude - Right Ascension) ]

Where Declination, Latitude, Longitude and Right Ascension are all given values.

Bean · 2013-10-17 11:34

tosjduenfs wrote: »

ALT = arcsin[ sin(Declination) * sin(Latitude) + cos(Declination) * cos(Latitude) * cos(GMST + Longitude - Right Ascension) ]

Where Declination, Latitude, Longitude and Right Ascension are all given values.

When you say "given values" do you mean the are constants ? If they are then the "sin" and "cos" can be precalculated. Leaving just a cos and arcsin function which could probably be done with CORDIC (not sure about arcsin though).

Bean

tosjduenfs · 2013-10-17 12:16

By given I meant they are input from the user but do not change with time. They can't be precalculated since they are subject to change when the user wants to locate a different object (RA/DEC) or change their location (LAT/LONG).

Tracy Allen · 2013-10-17 14:13

The usual way to deal with arcsin(x) is via the identity,
arcsin(x) = arctan(x / (1² - x²))

Note that 1 is squared in the formula, a detail, but important in implementation through a process of normalization as done in floating point.

F32 uses 24 bit CORDIC for the atan function, but it uses HUB table lookup/interpolation to 16 bits at best for sin, cos and tan. Thinking about Lawson's comments, I'm unclear on how much precision that last formula for ALT needs once you have an accurate value of GMST.

tosjduenfs · 2013-10-17 21:48

IT WORKS!!!! Tracy you were correct, having a more accurate value of GMST was all that I needed, thank you for your insight.

Bean, thank you so much for taking the time to write out examples and for showing me how to get an accurate value from the GMST formula.

Not only does it work but my final values are accurate to better than 1/100th of a degree often times less than that. That is more than adequate for my application.

Thanks again.

Mike

preserving as much precision as possible

Comments