Julien Date

This is another example to show the importance of having enough precision for your use case. It calculates the Julien Date, a continuous count of days and time from the beginning of the Julian period. This is used in astronomy, e.g. for calculating position of sun, moon, planets and stars for a given time and date.

The following snippet calculates the Julien Date both with the builtin precision of double variables and with fp64lib.

#include <fp64lib.h>
	
/* 
 * Another example for comparing precision between double and fp64lib,
 * calculating julien date for a given day and time. This is used as
 * base for many astronomical calculations, like determining position of the moon.
 *
 */

#define TIMEZONE	2		// delta in hours to UTC/ Greenwich time
#define LOCALDATE	20250613	// date in YYYYMMDD
#define LOCALTIME	180000		// time in HHMMSS (24 hours clock)

bool verbose = true;

// convert date to julien date
float64_t fp64_julienDate( int32_t date, int32_t timezone, uint32_t time ) {
	int day = date % 100; date /= 100;
	int month = date % 100;
	int year = date / 100;
	if( month <= 2 ) {
		month += 12;
		year--;
	}
	float64_t fday = fp64_uint32_to_float64( day );
	float64_t fmonth = fp64_add( fp64_uint32_to_float64( month ), float64_NUMBER_ONE	);
	float64_t fyear = fp64_uint32_to_float64( year );
	float64_t res1 = fp64_floor( fp64_mul( fyear, fp64_atof("365.25" ) ) );
	float64_t res2 = fp64_floor( fp64_mul( fmonth, fp64_atof("30.6001") ) );
	float64_t res3 = fp64_add( fday, fp64_atof("1720981.5") );
	float64_t res4 = fp64_div( fp64_uint32_to_float64( time ), fp64_uint32_to_float64(10000));
	res4 = fp64_div( fp64_sub( res4, fp64_uint32_to_float64( TIMEZONE )), fp64_uint32_to_float64(24));
	if( verbose ) {
		Serial.print("  fp64 res1:"); Serial.println( fp64_to_string( res1, 17, 15 ) );
		Serial.print("  fp64 res2:"); Serial.println( fp64_to_string( res2, 17, 15 ) );
		Serial.print("  fp64 res3:"); Serial.println( fp64_to_string( res3, 17, 15 ) );
		Serial.print("  fp64 res4:"); Serial.println( fp64_to_string( res4, 17, 15 ) );
	}
	float64_t res  = fp64_add(fp64_add(fp64_add(res1, res2),res3),res4);
	return res;
	}

double dbl_julienDate( int32_t date, int32_t timezone, uint32_t time ) {
	int day = date % 100; date /= 100;
	int month = date % 100;
	int year = date / 100;
	if( month <= 2 ) {
		month += 12;
		year--;
	}
	long res1 = year * 365.25;
	long res2 = (month+1) * 30.6001;
	double res3 = day + 1720981.5;
	double res4 = (time / 10000.0 - TIMEZONE) / 24.0;
	if( verbose ) {
		Serial.print("double res1:"); Serial.println(res1);
		Serial.print("double res2:"); Serial.println(res2);
		Serial.print("double res3:"); Serial.println(res3,15);
		Serial.print("double res4:"); Serial.println(res4,15);
	}
	double res  = res1 + res2 + res3 + res4;
	return res;
	}
	
void setup() {
	Serial.begin(57600);
	
	float64_t fp64_jul = fp64_julienDate( LOCALDATE, TIMEZONE, LOCALTIME );
	double 	  dbl_jul  = dbl_julienDate( LOCALDATE, TIMEZONE, LOCALTIME );
	
	Serial.println();
	Serial.print( "Julien date of " ); 
	Serial.print( LOCALDATE ); Serial.print(" "); 
	Serial.print( LOCALTIME ); Serial.print(" "); Serial.println( TIMEZONE );
	Serial.print("  fp64: "); Serial.println( fp64_to_string( fp64_jul, 17, 15 ) );
	Serial.print("double: "); Serial.println( dbl_jul, 15 );

	float64_t diff = fp64_sub( fp64_sd(dbl_jul), fp64_jul );
	Serial.print( "Difference " ); 
	Serial.print( fp64_to_string( diff, 17, 15 ) );
	float64_t diffMin = fp64_mul( diff, fp64_uint32_to_float64(24*60) );
	Serial.print( " = " ); Serial.print( fp64_to_string( diffMin, 17, 4 ) );
	Serial.println( " minutes!" );
	Serial.println();
}

void loop() {
}

Example output:

  fp64 res1:739631
fp64 res2:214
fp64 res3:1720994.5
fp64 res4:0.666666666666667
double res1:739631
double res2:214
double res3:1720994.500000000000000
double res4:0.666666698455810

Julien date of 20250613 180000 2
fp64: 2460840.166666667
double: 2460840.250000000000000
Difference 0.083333333488554 = 120.0000002235174 minutes!

Interestingly, the partial results res1 to res3 are identical, only res4 shows a difference after the 7th digit. However, the sum of the 4 partial results is wildly off. Precision might be improved slightly by changing the sequence of adding up the intermediate results.

However, precision of double (on the MegaAVR architecture) will never be sufficient to calculate Julien Dates down to a second or even a minute, thus ranking it out for usage in tracking stars (or even the moon) with a telescope – unless you heavily modify the algorithms.

The source code example provided with the library shows you, how the Julien Date changes on a per second level. Here is the start of it’s output:


20250613 180001 2: 2460840.166670833
20250613 180002 2: 2460840.166675
20250613 180003 2: 2460840.166679167
20250613 180004 2: 2460840.166683333
20250613 180005 2: 2460840.1666875
20250613 180006 2: 2460840.166691666
20250613 180007 2: 2460840.166695833
20250613 180008 2: 2460840.1667
20250613 180009 2: 2460840.166704167
20250613 180010 2: 2460840.166708333
20250613 180011 2: 2460840.1667125

As you can see, the change happens on the 6th decimal digit or 13th total digit. The best double/float can give you is 9 total digits.