Shop OBEX P1 Docs P2 Docs Learn Events
How to deal with large numbers in the PROP? Calculating Julian Day — Parallax Forums

How to deal with large numbers in the PROP? Calculating Julian Day

tosjduenfstosjduenfs Posts: 37
edited 2013-10-13 20:45 in Propeller 1
I'm working on a program that converts Right Ascension and Declination coordinates to Altitude and Azimuth coordinates for use in an auto tracking telescope. As my calculations go on the results get necessarily larger and larger. I'm stuck now calculating days since Julian Date J2000. When I figure the number on my calculator I get 2456577.51006944 days but when I run the calcs on the prop I get 2456578. The decimal portion is very import to further calculations. I realize I'm probably running in to the limit of 32bit precision, is there a way to overcome these issues?

If someone could tell me how to add a code block in my post I will post my code.

Thanks

Mike

Comments

  • Phil Pilgrim (PhiPi)Phil Pilgrim (PhiPi) Posts: 23,514
    edited 2013-10-12 11:44
    Why not do your calculations in Julian seconds instead of Julian days? That way you won't have to worry about the stuff after the decimal point.

    -Phil
  • Duane DegnDuane Degn Posts: 10,588
    edited 2013-10-12 11:45
    I don't think there's going to a simple way around this.

    If you search the forum for Julian Date, I think you'll find some recent discussion on this issue.

    There's a link to Rich's forum search engine at the top of post #1 of my index (see my signature for link).

    You may need to use floating point numbers in your calculations. I'm sure I have a link to the F32 thread (a better faster floating point object) in my index IIRC, it's in post #5 (but I'm not sure).

    There's a link to a tutorial on how to post code in the forum in post #3 of my index.
  • tosjduenfstosjduenfs Posts: 37
    edited 2013-10-12 12:37
    I can look into Julian seconds but I'm pretty sure that would put me way over the limit of a signed 32 bit number. Also I am using floating point numbers and the Float32Full object.
  • Tracy AllenTracy Allen Posts: 6,664
    edited 2013-10-12 13:10
    The solar calculator by Prophead100 has to deal with similar issues...
    http://obex.parallax.com/object/551

    Also see its thread...
    http://forums.parallax.com/showthread.php/136178-Solar-Tracking-Object
  • Tracy AllenTracy Allen Posts: 6,664
    edited 2013-10-12 13:18
    Julian seconds is good for over 100 years in 32 bits. A lifetime!

    I've written ordinal time methods for real time data logging, all integer math.
    [SIZE=1][FONT=courier new]CON
    {The following are ordinal time methods, from RTC_ISL12020   (c) emesystems release under MIT license}
    PUB daySecondNow
      ' returns the current second of day from 0 to 86399 read from RTC current time
        return jdsCalc(timeNowBCD)
    
    PUB jdsCalc(propPnt) | idx
    ' point to a buffer that contains BCD seconds, minutes, hours
    ' return the number of seconds since midnight represented by that time value,
       'repeat idx from 2 to 0
          result := BCD2INT(byte[propPnt][2]) * 3600 + BCD2INT(byte[propPnt][1]) * 60 + BCD2INT(byte[propPnt])
    
    PUB yearDayNow
      ' returns the current day of the year, read from RCT current time
        return jydCalc(timeNowBCD+3)
    
    
    PUB jydCalc(proppnt) | yy,mm,dd
    ' point to a buffer that contains BCD day, month, year  in that order, 3 BCD bytes
    ' returns with the sequential number of that day of year, zero to 365 or 366
      yy := BCD2INT(byte[propPnt][2])
      mm := BCD2INT(byte[propPnt][1])
      dd := BCD2INT(byte[propPnt])
      return ((mm-1)*30) + ((mm/9+mm)/2) - ((mm <# 3) / 3)*((yy//4 <# 1)+ 1) + dd
    
    PUB epochSecondNow
      ' returns the current number of seconds since midnight on Jan. 1, 2001
      jysCalc(timeNowBCD)
    
    PUB jysCalc(proppnt) | jd, js, yy
    ' point to a buffer that contains BCD second, minute, hour, day, month, year  in that order, 6 BCD bytes
    ' returns with the sequential number of seconds with number one at one second after midnight on Jan. 1, 2001
    ' does not account for leap seconds, good through 2099.
      js := jdsCalc(proppnt)
      jd := jydCalc(proppnt+3)
      yy := BCD2INT(byte[propPnt][5])
      return (((yy-1)*365) +((yy-1)/4) + jd) * 86400 + js[/FONT][/SIZE]
    
  • tosjduenfstosjduenfs Posts: 37
    edited 2013-10-12 18:46
    I tried julian seconds but the value is far too large for a 32 bit integer. Julian seconds for 10/11/2013 19:14:30 (when I started writing my code) is 212,248,296,870 which is too large of a number to calculate with 32 bit math.

    In the end what I really need is Greenwich Mean Sidereal Time or my Local Hour Angle. The only way I know how to get GMST or LHA is a calculation starting with Julian days or seconds.

    Tracy, I'm sorry if this is an ignorant question but I've only really starting learning about Julian dates, sidereal time etc very recently. I'm not exactly sure how your code works. I see that you calculate seconds from midnight Jan. 1, 2001 but shouldn't a Julian date be calculated from J2000 or Noon Jan. 1, 2000? I realize it would be trivial to find the difference between the two dates but even if you do calculate from midnight 2001 rather than noon 2000 I don't see how you can get less than 212 billion seconds.
  • Phil Pilgrim (PhiPi)Phil Pilgrim (PhiPi) Posts: 23,514
    edited 2013-10-12 19:00
    tosjduents wrote:
    I tried julian seconds but the value is far too large for a 32 bit integer. ...
    231 - 1 = 2_147_483_647. That many seconds is 2_147_483_647 / (365.25 * 24 * 60 * 60) = 68 years. From J2000, that wil get you to 2068.

    -Phil
  • tosjduenfstosjduenfs Posts: 37
    edited 2013-10-12 21:26
    Ok. I'm beginning to realize that my problems stem more from my lack of knowledge of what the Julian date actually is. I think I have come up with a solution, while more tedious, will allow me to keep the numbers reasonable. I've been using a formula I found to find the Julian date, in the course of calculations it uses some really enormous numbers. I'll use a more direct approach and simply figure it out myself.
  • Tracy AllenTracy Allen Posts: 6,664
    edited 2013-10-12 22:02
    I was dividing into 232 to come up with 136 years, but true, 231 and 68 years as Phil pointed out is more manageable for the Prop.

    I use ordinal time in my programs as a means to start data logging at a certain date/time, for a certain interval, and to stop at another date/time. Ordinal time avoids the complications that come from the odd number system of transitioning from day to day, month to month, year to year. The starting point does not matter. It could be like the PC starting in 1900, or the Mac in 1904, or Unix epoch 1970. I think I started mine in 2001 because I wrote my first version of it for the BASIC Stamp in 2001 and needed only to keep track of future time. I really shouldn't call it julian day number. It is just a scheme of ordinal time that can be converted to other schemes by addition of a constant.

    Correct me if I'm wrong here, but for the astronomer, julian day #1 started at noon on January 1, 4713 BC. In that system, the new year's party on Jan 1 2001 occurred at julian day number 2451910.5. Multiplying that by the number of seconds in a day does in fact exceed 32 bits by a long shot. However, it does fit easily into 63 bits, 19- decimal digits. Phil made an object, "Umath" that contains various 64 bit integer methods. AFAIK, there is no 64 bit floating point for the Prop. A few threads have popped up here about sidereal time calculations, but I'm out of my element, not much help. Maybe you can deal with the integer and fractional parts separately?
  • tosjduenfstosjduenfs Posts: 37
    edited 2013-10-12 22:35
    Tracy, yes the formulas I was using to find the number of days from J2000 started by calculating the days from the first julian day Jan 1, 4713BC and then subtracted the number of days from 4173BC to J2000. Those formulas are intended to be used when calculating days from any epoch past or present. That's where I was running into problems with the 212 billions seconds. Now I am calculating how many days it has been since Jan 1, 2000 myself with no concern for future or past epochs. Why the hell not since it will be accurate enough for the rest of my lifetime.

    Thank you guys for your help in pushing me in the right direction my calculations are coming out much better now.
  • RaymanRayman Posts: 14,662
    edited 2013-10-13 08:18
    For something math intensive, you could try using PropGCC instead of Spin...

    One advantage is that you could test the math on a PC and then use the exact same code on the Prop...
  • YanomaniYanomani Posts: 1,524
    edited 2013-10-13 13:42
    tosjduenfs wrote: »
    Ok. I'm beginning to realize that my problems stem more from my lack of knowledge of what the Julian date actually is. I think I have come up with a solution, while more tedious, will allow me to keep the numbers reasonable. I've been using a formula I found to find the Julian date, in the course of calculations it uses some really enormous numbers. I'll use a more direct approach and simply figure it out myself.

    tosjduenfs

    Your problem and these posts acted as some kind of wayback machine has catched me up, in a time travel till early 70's.
    It was one of the first programs I'd coded, to calculate days between dates; kind of bond trading machine, for a banking institution.
    Good times, good remembrances. Thank you for posting this.

    Perhaps the following link can add a bit of understanding to the whole problem, including concerns about correct mathematical representation of values and scales.
    There you'll also find a Julian Date Converter, to check your values and assumptions.

    http://aa.usno.navy.mil/data/docs/JulianDate.php

    I hope it can help you a bit

    Yanomani
  • tosjduenfstosjduenfs Posts: 37
    edited 2013-10-13 20:45
    Thanks Yanomani. the USNO is an excellent resource. At this point I have a properly functioning spreadsheet that can calculate the values I need. Now I am just working on doing the same thing in spin.
Sign In or Register to comment.