Flecks of Rust #6: Math fun in Rust with geodesy

Math in Rust is a little bit different in Rust than other languages you may be familiar with. For example in Python if you want to compute the sine for an angle of 55 degrees, you would first import the math module, convert the angle to radians and then use the value of the math.sin() function:

>>> # NOTE: Python example! Not Rust!
>>> import math
>>> angle = 55
>>> angle_rad = math.radians(angle)
>>> math.sin(angle_rad)
0.8191520442889917

In Rust you use the methods of the numeric type you are working with. For example, the f32 type is good for math unless you absolutely need the added precision of the f64 type. So, you could declare your angle as f32. Then you would use the to_radians function of that type to convert degrees to radians, and finally you would take the sine of that:

fn main() {
    let angle = 55f32;
    let angle_rad = angle.to_radians();
    println!("{}", angle_rad.sin());
}

0.81915206

You will soon get used to starting with the value.

The Rust Cookbook has an example of computing the distance between two geographic locations using the Haversine formula. The result is accurate for the purpose, as it is rarely possible to follow the shortest path when travelling from one location to another. However, the example uses a fixed value for Earth's radius, but that can vary depending on the latitude.

Here is a function to compute the value of Earth's radius in kilometers at a given latitude. You could pass in the latitude component of either of the geographic coordinates, or perhaps their average. The function also depends on two constants defined as values in kilometers.

const EARTH_EQUATORIAL_RADIUS: f32 = 6378.1370;
const EARTH_POLAR_RADIUS: f32 = 6356.7523;

// Computes the geocentric radius of the Earth for the given latitude.
// See http://en.wikipedia.org/wiki/Earth_radius#Geocentric_radius
// See https://rust-lang-nursery.github.io/rust-cookbook/science/mathematics/trigonometry.html
fn earth_geocentric_radius(latitude_degrees: f32) -> f32 {
    let latitude = latitude_degrees.to_radians();

    let a = EARTH_EQUATORIAL_RADIUS;
    let b = EARTH_POLAR_RADIUS;

    let t1a = a.powi(2) * latitude.cos();
    let t1b = b.powi(2) * latitude.sin();
    let t1 = t1a.powi(2) + t1b.powi(2);

    let t2a = a * latitude.cos();
    let t2b = b * latitude.sin();
    let t2 = t2a.powi(2) + t2b.powi(2);

    (t1 / t2).sqrt()
}

If you replace the value of earth_radius_kilometer in the Rust Cookbook example with the value obtained from the earth_geocentric_radius function above (with the latitude given as the average of the latitudes of Paris and London), you will find that the result differs by 0.3 kilometers, which is hardly significant in this case:

Earth's geocentric radius for latitude 50.18097 degrees is 6365.565 kilometers
Distance between Paris and London on the surface of Earth is 334.7 kilometers

What is more interesting here is the way that the various mathematical functions are used in Rust. The standard library defines the usual trigonometic functions for the floating-point types, and many more.

Since the formula for computing Earth's geocentric radius requires some values to be squared, the function uses the powi function of the f32 type, which uses an integer exponent, and should be faster than the powf function which allows you to use a floating-point exponent.

Another special feature of Rust worth mentioning is the return value of the function. While Rust does have a return statement, the idiomatic way is to let the last expression in the function act as the return value. Here the expression (t1 / t2).sqrt() is that return value.

Finally, here is a complete example that implements the Haversine formula as a function. It gives you another idea of how the Rust mathematical functions are used.

// A coordinate point, with latitude and longitude in degrees.
struct Coordinates {
    latitude: f32,
    longitude: f32,
}

// Computes the distance between two geographical points in kilometers.
// The coordinate values are expressed in degrees.
// Uses the haversine formula (see http://en.wikipedia.org/wiki/Haversine_formula).
// This implementation is derived from the JavaScript code
// by Andrew Hedges at http://andrew.hedges.name/experiments/haversine/.
// The Earth radius defaults to the Earth's equatorial radius.
// For added accuracy, you may want pass an Earth radius value more suitable for the
// latitudes of the points you are processing.
fn haversine(point1: Coordinates, point2: Coordinates, earth_radius: f32) -> f32 {
    let lat1 = point1.latitude.to_radians();
    let lon1 = point1.longitude.to_radians();
    let lat2 = point2.latitude.to_radians();
    let lon2 = point2.longitude.to_radians();

    let dlon = lon2 - lon1;
    let dlat = lat2 - lat1;
    let a = ((dlat / 2.0).sin()).powi(2) + lat1.cos() * lat2.cos() * ((dlon / 2.0).sin()).powi(2);
    let c = 2.0 * a.sqrt().atan2((1.0 - a).sqrt());
    earth_radius * c
}

fn main() {
    let paris = Coordinates { latitude: 48.85341, longitude: -2.34880 };
    let london = Coordinates { latitude: 51.50853, longitude: -0.12574 };

    let average_latitude = (paris.latitude + london.latitude) / 2.0;
    let radius = earth_geocentric_radius(average_latitude);
    println!("Earth's geocentric radius for latitude {} degrees is {} kilometers", average_latitude, earth_geocentric_radius(average_latitude));

    let distance = haversine(paris, london, radius);
    println!("Distance between Paris and London on the surface of Earth is {:.1} kilometers", distance);
}