NTP with GPS and PPS working without network clock

/in /by

Is it possible to have a PPS disciplined NTP server without a network reference clock or network connection?

While doing some camera synchronisation and timestamp testing, I setup a Raspberry Pi NTP server with PPS using an Adafruit GPS hat so that I could sync my cameras’ time to NTP, trigger them using the PPS signal, generate UTC timestamps and check the timestamps against the expected UTC second boundaries.

As part of this the question once again arose – Can I configure an NTP server with GPS time and PPS without having to configure a network reference clock along with the GPS and PPS sources?  This was important as my test setup didn’t have access to an external network, and yet I wanted typical NTP offsets of < 100us.

Well with quite a bit of messing I finally got it working, at this stage there are lots of pages that detail how to setup a Pi as an NTP server with PPS, some mention that you must have a network reference and some mention that you don’t – I had to do a lot of digging and a lot of messing to get to the bottom of things.

Amusingly, NTP configuration is something that ChatGPT is simply terrible at – don’t use it for this, it will just waste your time!

I had lots of problems, but the 2 most persistent was:

  1. The PPS source being marked as a false ticker (with an x to its left in the ntpq -p output)
  2. The PPS source not being referenced at all, “when” would remain at “-“

For ages, I lurched between problems 1 and 2…

Eventually I got it working for my setup.  The trick to get it working for me is that:

  • I needed 3 reference clocks in my setup for the PPS source to be used: PPS time, GPS time and a 3rd reference, normally this is a network NTP server but I found that configuring the local clock (127.127.1.0) worked instead.
  • I needed to be careful with my stratum settings, what worked for me was to put the local clock at stratum 15, GPS at stratum 10 and PPS at stratum 0.
  • Marked the PPS source as prefer (and not the GPS source as is often recommended?)

My working configuration:

tos mindist 0.500
# Local Clock
server 127.127.1.0
fudge 127.127.1.0 stratum 15
# GPS/NMEA Source from Hat
server 127.127.28.0 minpoll 4 maxpoll 4 prefer iburst
fudge 127.127.28.0 time1 0.000 refid GPS stratum 10
# PPS Source, stratum 0 and prefer
server 127.127.22.0 minpoll 4 maxpoll 4
fudge 127.127.22.0 refid PPS stratum 0 flag3 1 flag4 1 prefer

I deliberately left the time offset for the GPS source at 0.0 rather than trying to chase down the time packet offset so that I could see an “honest” offset value in the ntpq output:

 remote refid st t when poll reach delay offset jitter
==============================================================================
LOCAL(0) .LOCL. 15 l 58m 64 0 0.000 0.000 0.000
*SHM(0) .GPS. 10 l 16 16 377 0.000 -106.45 1.145
oPPS(0) .PPS. 0 l 15 16 377 0.000 -0.001 0.001

Flutter App no longer building after upgrading to Android Studio Electric Eel (Java errors)

/in /by

I decided to upgrade to Android Studio Electric Eel today in order to get WiFi pairing and debugging with my android phone – happily WiFi pairing now works, but unhappily my App refused to build anymore with Java version errors.

If when you run flutter doctor and get an error like this:

 

 

 

 

 

 

Take a look at this thread, linking from jbr to jre directories fixed it for me (windows 10), although I had to delete the empty C:\Program Files\Android\Android Studio\jre folder first!

From a CMD session with Admin privileges run:

mklink /D "jre" "jbr"

 

Flutter on Web does not support HTTP streaming (e.g. no support for MJPEG streaming)

/in /by

I have been looking into the possibility of using Flutter and Dart to implement a control GUI for a camera, the main target is a web app, I am looking at flutter as an alternative to  web frameworks like React.  The flutter platform is very impressive and Dart itself is a nice little language.  Flutter promises that many platforms can be targeted with one code-base – iOS, Android, Windows, Mac and Web.  Impressive claims indeed, although I am normally quite suspicious of such single-codebase claims as I have been very disappointed in the past!

As part of my feasibility project I wanted to get an MJPEG video stream displayed  in the UI to show a live camera image, this was quite easy to do using the  nice flutter_mjpeg package, it was fairly easy getting it to work targeting windows, but it wouldn’t work at all on chrome, after some research I discovered that flutter does not support HTTP/s streaming and keep-alives on web at all the HttpClient just times-out waiting for the response to complete, which of course it wont when streaming – so if you are struggling with it you may as well give up for now!

This is quite disappointing as this functionality is easy to achieve with a “traditional” javascript based Web app, hopefully it will be implemented in the future!  In the meatime I will continue working with flutter and dart as it looks looks like a very capable platform.

Some more info:

https://github.com/flutter/flutter/issues/84252

https://github.com/dart-lang/http/issues/595

C++ toolkit for analyzing GPS path data – distance, speed, heading etc.

/in /by

I have been working on and off on a simple C++ toolkit for working with GPS position path data, the type of data that might have been recorded on a mountain hike or run etc. the project is called gps_path_tools and can be found on GitHub.

The toolkit is a small collection of types and functions that allow you to do things like:

  • Calculate distance between two GPS locations
  • Calculate the length of a path
  • Calculate path “heading” at a point on the path
  • Find stationary points within a path
  • Find the closest point on the path to a given location
  • Find the speed of travel at a point on the path
  • Some basic coordinate transformation functions, e.g. to convert decimal-degrees-and-minutes (DDM) to decimal-degrees (DD).

It is a work in progress but has proven very useful for me so far when I have been analyzing previous mountain hikes.

The toolkit is contained in a single header (.h) file and has no dependencies on any libraries other than the standard C++ library.

Pylon EnumerateDevices() not finding Basler GigE cameras

/in /by

I hit a problem recently whereby CTlFactory::EnumerateDevices() would not find my connected Balser cameras while Pylon Viewer would find and connect to them OK.  On investigation I discovered that if WiFi was enabled on the host, then the GigE Vision Device discovery broadcast messages would not be sent over GigE to the cameras.  With WiFi disabled the cameras would be found but when it is enabled they wouldn’t be found!

After a good bit of googling I discovered that I should call IGigETransportLayer::EnumerateAllDevices() rather than the usual CTlFactory::EnumerateDevices(), here is some code that demonstrates the switch:

 

#include <pylon/gige/GigETransportLayer.h>
...
auto& factory = CTlFactory::GetInstance();
// This does not always work so we use EnumerateAllDevices() below instead..
// factory.EnumerateDevices();
auto* transport_layer = (IGigETransportLayer*)factory.CreateTl(Pylon::BaslerGigEDeviceClass);
    
transport_layer->EnumerateAllDevices(devices);

 

 

Code to calculate heading / bearing from two GPS latitude and longitude Points

/in /by

Example code to Estimate heading or course  of travel or bearing given two GPS locations taken from a recorded path.

A lot of GPS devices will report heading along with latitude and longitude, but what if you lust have a series of GPS locations recorded from a path, how can you estimate heading or course of travel at points along the path  given two consecutive GPS locations from the path?

Why might you want to know heading along a path?  Well imagine you have a mobile computer vision application, the effect of the Sun’s illumination may have a large effect on your imaging, if you can calculate the Sun’s position in the sky during the recorded path, and also calculate the camera’s “heading” then you can estimate the effect that the Sun will have on your images, for example are they back-lit, front-lit, will long shadows be cast etc.

The background to the calculation can be found here under the section called ‘bearing’.

Note that the terms ‘heading’, ‘bearing’ and ‘course’ have subtly different meanings especially for vehicles that might drift or yaw like boats or aircraft, but here we use them interchangeably to refer to the movement of the GPS measurement center.  As the points we choose are close together we don’t really have to worry about the difference between initial bearing and final bearing.

Here is some code that will calculate heading from 2 locations.  Pass in the latitude and longitude of both points, the order in which the points are passed will affect the calculated heading.  The two points that you choose should be close together to get a good estimate of ‘instantaneous’ heading, however if they are too close together the estimate may be noisy due to the positional accuracy limits of GPS, a distance of 100m to 200m could work well, experiment for your GPS device and path type.

 

#define _USE_MATH_DEFINES
#include <math.h>
#include "calculate_gps_heading.h"
//   Copyright 2022 Kevin Godden
//
//   Licensed under the Apache License, Version 2.0 (the "License");
//   you may not use this file except in compliance with the License.
//   You may obtain a copy of the License at
//
//       http://www.apache.org/licenses/LICENSE-2.0
//
//   Unless required by applicable law or agreed to in writing, software
//   distributed under the License is distributed on an "AS IS" BASIS,
//   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
//   See the License for the specific language governing permissions and
//   limitations under the License.
static inline double to_rad(double theta) {
    return (theta * M_PI) / 180.0;
}
static inline double to_degrees(double theta) {
    return (theta * 180.0) / M_PI;
}
//
// Calculate the heading in decimal degrees between 2 (preferably quite close) locations 
// (lat1, lon1) --> First location
// (lat2, lon2) --> Second location
//
double calculate_gps_heading(double lat1, double lon1, double lat2, double lon2) {
    // Convert degrees to radians
    lat1 = to_rad(lat1);
    lon1 = to_rad(lon1);
    lat2 = to_rad(lat2);
    lon2 = to_rad(lon2);    
        
    double dlon = lon2 - lon1;
    double X = cos(lat2) * sin(dlon);
    double Y = cos(lat1) * sin(lat2) - sin(lat1) * cos(lat2) * cos(dlon);
    
    double heading = atan2(X,Y);
    
    // We want heading in degrees, not radians.
    heading = to_degrees(heading);
  
    // We want a uniform heading of >=0 and <360
    if (heading < 0)
        heading = 360.0 + heading;
    
    return heading;
}

The code along with some tests can be found in this Git repo.

Code Example: Convert Decimal Degrees and Minutes (DDM) to Decimal Degrees (DD)

/in /by

How to convert between Decimal Degrees & Minutes (DDM) to Decimal Degrees (DD).

I always prefer to represent latitude and longitude as full decimal degrees (DD), I find it much simpler and you can do calculations with the values.  However, other people and devices love to use the slightly (to me) esoteric Decimal Degrees and Minutes (DDM) form and I have to remember how to convert to and from decimal degrees – so here’s a short description of how it’s done along with some C++ sample code.

Decimal Degrees: Latitude and longitude are just represented as a decimal degree value.

Decimal Degrees & Minutes: Encodes a whole number of degrees in the first few digits, the rest of the digits represent a decimal minutes value with 2 digits to the left of the decimal place representing a whole minute value,  e.g. 5254.9098940 : ddmm.mmmmmm

For example, a DDM value of: 5254.9098940 represents 52 degrees and 54.9098940 minutes.

To convert this to decimal degrees, we do something like this:

dd = 52 + (54.9098940 / 60) = 52.9151649 degrees.

And in code:

//
double ddm_to_dd(double ddm) {
    double degrees = floor(ddm / 100.0);
    double minutes = ddm - degrees * 100.0;
    double decimal_degrees = degrees + minutes / 60.0;
    return decimal_degrees;
}

 

 

 

 

Simple Wrapper of boost::asio for receiving UDP datagrams

/in /by

More easily receive UDP Data via boost::asio providing just IP address and port umber.

The boost::asio classes provide a cross platform way to communicate via UDP without having to use the raw socket libraries, however they can be confusing and hard to use do to their designed-in flexibility.

Sometimes I just want to provide an IP address and a port number and just receive messages so I have written a wrapper class called boost_udp_receive_rar which can be found on GitHub here.

It has 4 functions that allow you to receive string or binary messages in a synchronous or asynchronous manner, for example to receive ASCII string data:

#include "boost_udp_send_faf.h"
#include <thread>
// Setup a receiver, specifying IP address and port
// Note the IP address is that of the receiving network
// interface.  As always with UDP and non-standard ports
// check that a firewall isn't blocking data transfer.
boost_udp_receive_rar rar("127.0.0.1", 8861);
// Receive a datagram synchronously as a string
//
std::string datagram = rar.receive_sync();
std::count << "received: " << datagram << std::endl;
// Asynchronously receive a datagram as a string
//
// We loop calling receive_async() until it returns
// a non-empty string. receive_async() returns quickly
// and we can do other work (or other async reads on different
// ports) while we wait for somethign to arrive.

do {
    datagram = rar.receive_async();
    // Do other stuff
    //
    // Better do a little sleep to be nice to the
    // CPU
    std::this_thread::sleep_for(std::chrono::milliseconds(100));
} while (datagram.empty());
std::count << "received: " << datagram << std::endl;

For an easier way to send UDP data see: boost_udp_send_faf

mktime() always uses local time zone and converts from UTC

/in /by

If you’re working with UTC times, you may be upset to find that mktime() always converts using your local time zone thus wrecking your UTC times.  Instead of mktime() use:

  • Linux/Posix: timegm()
  • Windows: _mkgmtime()

For example, let’s use timegm() to make a UTC time which we will pass to SolarAzEl() to calculate the Solar Elevation at that time, it must be passed a UTC time!  See this post if you are interested in SolarAzEl().

//
// Make UTC time for 10:16:00 on 5th. May 2022
// and use it to calculate the Solar Elevation
// at the time via SolarAzEl().  SolarAzEl() 
// requires a UTC time!
//
tm utc;
// tm_year is time since 1900
utc.tm_year = 2022 - 1900;
// Month is zero based, i.e. Jan is month 0, May is 4
utc.tm_mon = 5 - 1;
utc.tm_mday = 5;
utc.tm_hour = 10;
utc.tm_min = 16;
utc.tm_sec = 00;
utc.tm_isdst = 0;
// Get UTC time_t val, do not change for local time offset!
tim = timegm(&utc);	// or _mkgmtime() on windows
//
// Now calculate Solar Elevation at
// GPS coordinates 52.975, -6.0494
// at sea level using the UTC time.
//
double altitude = 0;
double Az = 0.0;
double El = 0.0;
double lat = 52.975;
double lon = -6.0494;
SolarAzEl(tim, lat, lon, 0, &Az, &El);
printf("Solar Azimuth: %f\n", Az);
printf("Solar Elevation: %f\n", El);
//

 

 

Estimate Solar Azimuth and Elevation given GPS position and time.

/in /by

C++ code to estimate Solar Azimuth and Elevation given GPS position and time.

For reasons that I won’t go into here I found it necessary to estimate Solar Azimuth and Elevation given GPS position and time.

Now there is quite a bit of information on how to this on the old inter-web, none of it could be described as being ‘easy’, and most methods use the ‘equation of time’ and require you to know your current time zone in order to to calculate Local Solar Time (LST) so that you can then calculate the current hour angle.

However dynamically determining which time zone you’re in is very difficult, most approaches involve using a web API – but what can you do if you have no internet connection?

Lucky I found some Mathlab code written by Darin C. Koblick which calculates Solar Azimuth and Elevation using just latitude, longitude, UTC time and altitude. It models the Sun orbiting the Earth and does not require to know your time zone! It seems to work very well and I am very happy with the results so far (thanks Darin!)

I ported the code over to C++, trying to change as little from the original Mathlab code as possible, I have included the code here, please note the license text at the end.

The code can be found in this Git repo.

For example, to estimate Azimuth and Elevation for Lat/Lon 52.975/-6.0494 at sea level for the current time:

/*
Copyright(c) 2010, Darin Koblick
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met :
*Redistributions of source code must retain the above copyright notice, this
list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright notice,
this list of conditions and the following disclaimer in the documentation
and/or other materials provided with the distribution
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED.IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
#include <time.h>
#include <math.h>
#include <stdio.h>
#ifndef M_PI
#define M_PI (3.14159265358979323846264338327950288)
#endif /* M_PI */
// Programed by Darin C.Koblick 2 / 17 / 2009
//
//              Darin C.Koblick 4 / 16 / 2013 Vectorized for Speed
//                                         Allow for MATLAB Datevec input in
//                                         addition to a UTC string.
//                                         Cleaned up comments and code to
//                                         avoid warnings in MATLAB editor.
//
//				Kevin Godden 9/1/2020      Ported from Matlab to C++, tried to change as little as possible.
//                                         this is a non-vectorised port.
//
//--------------------------------------------------------------------------
//
// External Function Call Sequence :
//
// double lat = 52.975;
// double lon = -6.0494;
// double altitude = 0;
//
// double Az = 0.0;
// double El = 0.0;
// SolarAzEl(time(NULL), lat, lon, 0, &Az, &El);
// 
// printf("Azimuth: %f\n", Az);
// printf("Elevation: %f\n", El);
//
// Or to calculate Az & El for an arbitary UTC time:
//
//
// tm utc;
// tm_year is time since 1900
// utc.tm_year = y - 1900;
// Month is zero based, i.e. Jan is month 0
// utc.tm_mon = m - 1;
// utc.tm_mday = d;
// utc.tm_hour = 10;
// utc.tm_min = 16;
// utc.tm_sec = 00;
// utc.tm_isdst = 0;
// 
// Get UTC time_t val
// tim = timegm(&utc);	// or _mkgmtime() on windows
// 
// double altitude = 0;
// double Az = 0.0;
// double El = 0.0;
// 
// double lat = 52.975;
// double lon = -6.0494;
// 
// SolarAzEl(tim, lat, lon, 0, &Az, &El);
// 
// printf("Az: %f\n", Az);
// printf("El: %f\n", El);
// 
//
// Function Description :
//
// SolarAzEl will ingest a Universal Time, and specific site location on earth
// it will then output the solar Azimuth and Elevation angles relative to that
// site.
//
// Input Description :
//
// utc_time_point : time_t containing target time for sun position calculations.
//
// Lat : Site Latitude in degrees -90:90->S(-) N(+)
//
// Lon : Site Longitude in degrees -180:180 W(-) E(+)
//
// Alt : Altitude of the site above sea level(Km)
//
// Output Description :
//  Az    Azimuth location of the sun(deg)
//  El    Elevation location of the sun(deg)
//
//
// Source References :
// Solar Position obtained from :
// http ://stjarnhimlen.se/comp/tutorial.html#5
//
double julian_day(time_t utc_time_point);
void SolarAzEl(time_t utc_time_point, double Lat, double Lon, double Alt, double* Az, double* El) {
    double jd = julian_day(utc_time_point);
    double d = jd - 2451543.5;
    
    // Keplerian Elements for the Sun(geocentric)
    double w = 282.9404 + 4.70935e-5*d; // (longitude of perihelion degrees)
    // a = 1.000000; % (mean distance, a.u.)
    double e = 0.016709 - 1.151e-9*d; // (eccentricity)
    double M = fmod(356.0470 + 0.9856002585*d, 360.0); // (mean anomaly degrees)
        
    double L = w + M; // (Sun's mean longitude degrees)
    double oblecl = 23.4393 - 3.563e-7*d; // (Sun's obliquity of the ecliptic)
    // auxiliary angle
    double  E = M + (180 / M_PI)*e*sin(M*(M_PI / 180))*(1 + e*cos(M*(M_PI / 180)));
    // rectangular coordinates in the plane of the ecliptic(x axis toward perhilion)
    double x = cos(E*(M_PI / 180)) - e;
    double y = sin(E*(M_PI / 180))*sqrt(1 - pow(e, 2));
    // find the distance and true anomaly
    double r = sqrt(pow(x,2) + pow(y,2));
    double v = atan2(y, x)*(180 / M_PI);
    // find the longitude of the sun
    double lon = v + w;
    // compute the ecliptic rectangular coordinates
    double xeclip = r*cos(lon*(M_PI / 180));
    double yeclip = r*sin(lon*(M_PI / 180));
    double zeclip = 0.0;
    //rotate these coordinates to equitorial rectangular coordinates
    double xequat = xeclip;
    double yequat = yeclip*cos(oblecl*(M_PI / 180)) + zeclip * sin(oblecl*(M_PI / 180));
    double zequat = yeclip*sin(23.4406*(M_PI / 180)) + zeclip * cos(oblecl*(M_PI / 180));
    // convert equatorial rectangular coordinates to RA and Decl:
    r = sqrt(pow(xequat, 2) + pow(yequat, 2) + pow(zequat, 2)) - (Alt / 149598000); //roll up the altitude correction
    double RA = atan2(yequat, xequat)*(180 / M_PI);
    double delta = asin(zequat / r)*(180 / M_PI);
    
    // Following the RA DEC to Az Alt conversion sequence explained here :
    // http ://www.stargazing.net/kepler/altaz.html
    //	Find the J2000 value
    //	J2000 = jd - 2451545.0;
    //hourvec = datevec(UTC);
    //UTH = hourvec(:, 4) + hourvec(:, 5) / 60 + hourvec(:, 6) / 3600;
    // Get UTC representation of time / C++ Specific
    tm *ptm;
    ptm = gmtime(&utc_time_point);
    double UTH = (double)ptm->tm_hour + (double)ptm->tm_min / 60 + (double)ptm->tm_sec / 3600;
    // Calculate local siderial time
    double GMST0 = fmod(L + 180, 360.0) / 15;
    double SIDTIME = GMST0 + UTH + Lon / 15;
    
    // Replace RA with hour angle HA
    double HA = (SIDTIME*15 - RA);
    // convert to rectangular coordinate system
    x = cos(HA*(M_PI / 180))*cos(delta*(M_PI / 180));
    y = sin(HA*(M_PI / 180))*cos(delta*(M_PI / 180));
    double z = sin(delta*(M_PI / 180));
    // rotate this along an axis going east - west.
    double xhor = x*cos((90 - Lat)*(M_PI / 180)) - z*sin((90 - Lat)*(M_PI / 180));
    double yhor = y;
    double zhor = x*sin((90 - Lat)*(M_PI / 180)) + z*cos((90 - Lat)*(M_PI / 180));
    
    // Find the h and AZ
    *Az = atan2(yhor, xhor)*(180 / M_PI) + 180;
    *El = asin(zhor)*(180 / M_PI);
}
double julian_day(time_t utc_time_point) {
    // Extract UTC Time
    struct tm* tm = gmtime(&utc_time_point);
    double year = tm->tm_year + 1900;
    double month = tm->tm_mon + 1;
    double day = tm->tm_mday;
    double hour = tm->tm_hour;
    double min = tm->tm_min;
    double sec = tm->tm_sec;
    if (month <= 2) {
        year -= 1;
        month += 12;
    }
    double jd = floor(365.25*(year + 4716.0)) + floor(30.6001*(month + 1.0)) + 2.0 -
        floor(year / 100.0) + floor(floor(year / 100.0) / 4.0) + day - 1524.5 +
        (hour + min / 60 + sec / 3600) / 24;
    return jd;
}