geodesy_module – fortran-astrodynamics-toolkit

Geodesy routines.

Uses

- numbers_module
- kind_module
Help

Graph Key

Nodes of different colours represent the following:

Solid arrows point from a submodule to the (sub)module which it is descended from. Dashed arrows point from a module or program unit to modules which it uses.

Used by

Help

Graph Key

Nodes of different colours represent the following:

Solid arrows point from a submodule to the (sub)module which it is descended from. Dashed arrows point from a module or program unit to modules which it uses.

Functions

public pure function great_circle_distance(r, long1, lat1, long2, lat2) result(d)

Author: Jacob Williams
Date: 7/13/2014

Great circle distance on a spherical body, using the Vincenty algorithm.

Arguments

Type	Intent		Name
real(kind=wp),	intent(in)	::	r	radius of the body [km]
real(kind=wp),	intent(in)	::	long1	longitude of first site [rad]
real(kind=wp),	intent(in)	::	lat1	latitude of the first site [rad]
real(kind=wp),	intent(in)	::	long2	longitude of the second site [rad]
real(kind=wp),	intent(in)	::	lat2	latitude of the second site [rad]

Return Value real(kind=wp)

great circle distance from 1 to 2 [km]

public pure function geocentric_radius(a, b, lat) result(r)

The distance from the center of a celestial body (e.g., the Earth) to a point on the spheroid surface at a specified geodetic latitude.

Arguments

Type	Intent		Name
real(kind=wp),	intent(in)	::	a	equatorial radius (km)
real(kind=wp),	intent(in)	::	b	polar radius of point (km)
real(kind=wp),	intent(in)	::	lat	geodetic latitude of point (rad)

Return Value real(kind=wp)

distance from center of body to point (km)

private pure function solve_polynomial(B, x0, error) result(x)

Numerical solution to polynomial equation using Newton-Raphson method

Arguments

Type	Intent	Attributes		Name
real(kind=wp),	intent(in),	dimension(0:6)	::	B	Polynomial `B = B(6) x^6 + B(5) x^5 + ... + B(0)`
real(kind=wp),	intent(in)		::	x0	Initial point
real(kind=wp),	intent(in)		::	error	Maximum error

Return Value real(kind=wp)

root found after applying Newton-Raphson method to B The function returns the value when the correction is smaller than error.

Subroutines

public pure subroutine heikkinen(rvec, a, b, h, lon, lat)

Author: Jacob Williams

Heikkinen routine for cartesian to geodetic transformation

Arguments

Type	Intent	Attributes		Name
real(kind=wp),	intent(in),	dimension(3)	::	rvec	position vector [km]
real(kind=wp),	intent(in)		::	a	geoid semimajor axis [km]
real(kind=wp),	intent(in)		::	b	geoid semiminor axis [km]
real(kind=wp),	intent(out)		::	h	geodetic altitude [km]
real(kind=wp),	intent(out)		::	lon	longitude [rad]
real(kind=wp),	intent(out)		::	lat	geodetic latitude [rad]

public pure subroutine olson(rvec, a, b, h, long, lat)

Author: Jacob Williams

Olson routine for cartesian to geodetic transformation.

Arguments

Type	Intent	Attributes		Name
real(kind=wp),	intent(in),	dimension(3)	::	rvec	position vector [km]
real(kind=wp),	intent(in)		::	a	geoid semimajor axis [km]
real(kind=wp),	intent(in)		::	b	geoid semiminor axis [km]
real(kind=wp),	intent(out)		::	h	geodetic altitude [km]
real(kind=wp),	intent(out)		::	long	longitude [rad]
real(kind=wp),	intent(out)		::	lat	geodetic latitude [rad]

public subroutine direct(a, f, glat1, glon1, faz, s, glat2, glon2, baz)

Author: Jacob Williams

Solve the "direct" geodetic problem: given the latitude and longitude of one point and the azimuth and distance to a second point, determine the latitude and longitude of that second point. The solution is obtained using the algorithm by Vincenty.

Arguments

Type	Intent		Name
real(kind=wp),	intent(in)	::	a	semimajor axis of ellipsoid [m]
real(kind=wp),	intent(in)	::	f	flattening of ellipsoid [-]
real(kind=wp),	intent(in)	::	glat1	latitude of 1 [rad]
real(kind=wp),	intent(in)	::	glon1	longitude of 1 [rad]
real(kind=wp),	intent(in)	::	faz	forward azimuth 1->2 [rad]
real(kind=wp),	intent(in)	::	s	distance from 1->2 [m]
real(kind=wp),	intent(out)	::	glat2	latitude of 2 [rad]
real(kind=wp),	intent(out)	::	glon2	longitude of 2 [rad]
real(kind=wp),	intent(out)	::	baz	back azimuth 2->1 [rad]

public subroutine geodetic_to_cartesian(a, b, glat, lon, h, r)

Author: Jacob Williams

Geodetic latitude, longitude, and height to Cartesian position vector.

Arguments

Type	Intent	Attributes		Name
real(kind=wp),	intent(in)		::	a	geoid semimajor axis [km]
real(kind=wp),	intent(in)		::	b	geoid semiminor axis [km]
real(kind=wp),	intent(in)		::	glat	geodetic latitude [rad]
real(kind=wp),	intent(in)		::	lon	longitude [rad]
real(kind=wp),	intent(in)		::	h	geodetic altitude [km]
real(kind=wp),	intent(out),	dimension(3)	::	r	Cartesian position vector [x,y,z]

public subroutine inverse(a, rf, b1, l1, b2, l2, faz, baz, s, it, sig, lam, kind)

INVERSE computes the geodetic azimuth and distance between two points, given their geographic positions.

Arguments

Type	Intent		Name
real(kind=wp),	intent(in)	::	a	Equatorial semimajor axis
real(kind=wp),	intent(in)	::	rf	reciprocal flattening (1/f)
real(kind=wp),	intent(in)	::	b1	latitude of point 1 (rad, positive north)
real(kind=wp),	intent(in)	::	l1	longitude of point 1 (rad, positive east)
real(kind=wp),	intent(in)	::	b2	latitude of point 2 (rad, positive north)
real(kind=wp),	intent(in)	::	l2	longitude of point 2 (rad, positive east)
real(kind=wp),	intent(out)	::	faz	Forward azimuth (rad, clockwise from north)
real(kind=wp),	intent(out)	::	baz	Back azimuth (rad, clockwise from north)
real(kind=wp),	intent(out)	::	s	Ellipsoidal distance
integer,	intent(out)	::	it	iteration count
real(kind=wp),	intent(out)	::	sig	spherical distance on auxiliary sphere
real(kind=wp),	intent(out)	::	lam	longitude difference on auxiliary sphere
integer,	intent(out)	::	kind	solution flag: kind=1, long-line; kind=2, antipodal

public subroutine geodetic_to_cartesian_triaxial(ax, ay, b, phi, lambda, h, r)

Function computes the Cartesian coordinates given the geodetic latitude (phi), longitude (lambda) and height (h) of a point related to an ellipsoid defined by its three semiaxes ax, ay and b

Arguments

Type	Intent	Attributes		Name
real(kind=wp),	intent(in)		::	ax	semiaxes (0 < b <= ay <= ax)
real(kind=wp),	intent(in)		::	ay	semiaxes (0 < b <= ay <= ax)
real(kind=wp),	intent(in)		::	b	semiaxes (0 < b <= ay <= ax)
real(kind=wp),	intent(in)		::	phi	geodetic latitude (radians)
real(kind=wp),	intent(in)		::	lambda	geodetic longitude (radians)
real(kind=wp),	intent(in)		::	h	geodetic height
real(kind=wp),	intent(out),	dimension(3)	::	r	Cartesian position vector [x,y,z]

public pure subroutine geodetic_to_cartesian_triaxial_2(a, b, c, lat, long, h, r)

Geodetic to Cartesian for Triaxial Ellipsoid.

Arguments

Type	Intent	Attributes		Name
real(kind=wp),	intent(in)		::	a	ellipsoid radii `a >= b >= c`
real(kind=wp),	intent(in)		::	b	ellipsoid radii `a >= b >= c`
real(kind=wp),	intent(in)		::	c	ellipsoid radii `a >= b >= c`
real(kind=wp),	intent(in)		::	lat	latitude (rad)
real(kind=wp),	intent(in)		::	long	longitude (rad)
real(kind=wp),	intent(in)		::	h	altitude
real(kind=wp),	intent(out),	dimension(3)	::	r	Cartesian coordinates (x,y,z)

public subroutine cartesian_to_geodetic_triaxial(ax, ay, b, r, tol, phi, lambda, h)

Function computes the geodetic latitude (phi), longitude (lambda) and height (h) of a point related to an ellipsoid defined by its three semiaxes ax, ay and b (0 < b <= ay <= ax) given Cartesian coordinates Xi, Yi, Zi and tolerance (tol). Latitude and longitude are returned in radians.

Arguments

Type	Intent	Attributes		Name
real(kind=wp),	intent(in)		::	ax	semiaxes (0 < b <= ay <= ax)
real(kind=wp),	intent(in)		::	ay	semiaxes (0 < b <= ay <= ax)
real(kind=wp),	intent(in)		::	b	semiaxes (0 < b <= ay <= ax)
real(kind=wp),	intent(in),	dimension(3)	::	r	Cartesian coordinates (x,y,z)
real(kind=wp),	intent(in)		::	tol	tolerance (may be set to zero)
real(kind=wp),	intent(out)		::	phi	geodetic latitude (radians)
real(kind=wp),	intent(out)		::	lambda	geodetic longitude (radians)
real(kind=wp),	intent(out)		::	h	geodetic height

private subroutine bisection_special_2(cz, Xo, Yo, tol, n, m, Gm)

Arguments

Type	Intent		Name
real(kind=wp),	intent(in)	::	cz
real(kind=wp),	intent(in)	::	Xo
real(kind=wp),	intent(in)	::	Yo
real(kind=wp),	intent(in)	::	tol
integer,	intent(out)	::	n
real(kind=wp),	intent(out)	::	m
real(kind=wp),	intent(out)	::	Gm

private subroutine bisection_special_3(cx, cy, Xo, Yo, Zo, tol, n, m, Hm)

Arguments

Type	Intent		Name
real(kind=wp),	intent(in)	::	cx
real(kind=wp),	intent(in)	::	cy
real(kind=wp),	intent(in)	::	Xo
real(kind=wp),	intent(in)	::	Yo
real(kind=wp),	intent(in)	::	Zo
real(kind=wp),	intent(in)	::	tol
integer,	intent(out)	::	n
real(kind=wp),	intent(out)	::	m
real(kind=wp),	intent(out)	::	Hm

private subroutine philambda_quadrant(x, y, z, phi, lambda)

Arguments

Type	Intent		Name
real(kind=wp),	intent(in)	::	x
real(kind=wp),	intent(in)	::	y
real(kind=wp),	intent(in)	::	z
real(kind=wp),	intent(inout)	::	phi
real(kind=wp),	intent(inout)	::	lambda

private subroutine xyz2philambda(ax, ay, b, x, y, z, phi, lambda)

Determination of the geodetic latitude and longitude

Arguments

Type	Intent		Name
real(kind=wp),	intent(in)	::	ax
real(kind=wp),	intent(in)	::	ay
real(kind=wp),	intent(in)	::	b
real(kind=wp),	intent(in)	::	x
real(kind=wp),	intent(in)	::	y
real(kind=wp),	intent(in)	::	z
real(kind=wp),	intent(inout)	::	phi
real(kind=wp),	intent(inout)	::	lambda

private subroutine xyz2fl(ax, ay, b, x, y, z, latitude, longitude)

Computes the transformation of Cartesian to geodetic coordinates on the surface of the ellipsoid assuming x,y,z are all non-negative Angular coordinates in radians

Arguments

Type	Intent		Name
real(kind=wp),	intent(in)	::	ax
real(kind=wp),	intent(in)	::	ay
real(kind=wp),	intent(in)	::	b
real(kind=wp),	intent(in)	::	x
real(kind=wp),	intent(in)	::	y
real(kind=wp),	intent(in)	::	z
real(kind=wp),	intent(out)	::	latitude
real(kind=wp),	intent(out)	::	longitude

private pure subroutine horner(B, c, BB)

Horner's method to compute B(x-c) in terms of B(x).

Arguments

Type	Intent	Attributes		Name
real(kind=wp),	intent(in),	dimension(0:6)	::	B	Polynomial `B = B(6) x^6 + B(5) x^5 + ... + B(0)`
real(kind=wp),	intent(in)		::	c
real(kind=wp),	intent(out),	dimension(0:6)	::	BB	Polynomial `BB` such that `B(x-c) = BB(x)`

public subroutine CartesianIntoGeodeticI(ax, ay, az, r, latitude, longitude, altitude, error)

Cartesian to Geodetic I

Arguments

Type	Intent	Attributes		Name
real(kind=wp),	intent(in)		::	ax	semiaxes of the celestial body: ax>ay>az
real(kind=wp),	intent(in)		::	ay	semiaxes of the celestial body: ax>ay>az
real(kind=wp),	intent(in)		::	az	semiaxes of the celestial body: ax>ay>az
real(kind=wp),	intent(in),	dimension(3)	::	r	cartesian coordinates of the considered point in the first octant: xG, yG, zG with (xG,yG,zG)<>(0,0,0)
real(kind=wp),	intent(out)		::	latitude	geodetic coordinates of the considered point
real(kind=wp),	intent(out)		::	longitude	geodetic coordinates of the considered point
real(kind=wp),	intent(out)		::	altitude	geodetic coordinates of the considered point
real(kind=wp),	intent(in)		::	error	Values smaller than error treated as 0.0

public subroutine CartesianIntoGeodeticII(ax, ay, az, r, latitude, longitude, altitude, error)

Cartesian into Geodetic II

Arguments

Type	Intent	Attributes		Name
real(kind=wp),	intent(in)		::	ax	semiaxes of the celestial body: ax>ay>az
real(kind=wp),	intent(in)		::	ay	semiaxes of the celestial body: ax>ay>az
real(kind=wp),	intent(in)		::	az	semiaxes of the celestial body: ax>ay>az
real(kind=wp),	intent(in),	dimension(3)	::	r	cartesian coordinates of the considered point in the first octant: xG, yG, zG with (xG,yG,zG)<>(0,0,0)
real(kind=wp),	intent(out)		::	latitude	geodetic coordinates of the considered point
real(kind=wp),	intent(out)		::	longitude	geodetic coordinates of the considered point
real(kind=wp),	intent(out)		::	altitude	geodetic coordinates of the considered point
real(kind=wp),	intent(in)		::	error	Values smaller than error treated as 0.0

public subroutine cartesian_to_geodetic_triaxial_2(a, b, c, r, eps, phi, lambda, h)

Cartesian to geodetic for Triaxial Ellipsoid.

Arguments

Type	Intent	Attributes		Name
real(kind=wp),	intent(in)		::	a	ellipsoid radii `a >= b >= c`
real(kind=wp),	intent(in)		::	b	ellipsoid radii `a >= b >= c`
real(kind=wp),	intent(in)		::	c	ellipsoid radii `a >= b >= c`
real(kind=wp),	intent(in),	dimension(3)	::	r	Cartesian coordinates (x,y,z)
real(kind=wp),	intent(in)		::	eps	convergence tolerance
real(kind=wp),	intent(out)		::	phi	latitude (rad)
real(kind=wp),	intent(out)		::	lambda	longitude (rad)
real(kind=wp),	intent(out)		::	h	altitude

private subroutine special_cases(a, b, c, x, y, z, phi, lambda, h, done)

Special cases for lat/lon/altitude

Arguments

Type	Intent		Name
real(kind=wp),	intent(in)	::	a	ellipsoid radii `a >= b >= c`
real(kind=wp),	intent(in)	::	b	ellipsoid radii `a >= b >= c`
real(kind=wp),	intent(in)	::	c	ellipsoid radii `a >= b >= c`
real(kind=wp),	intent(in)	::	x	Cartesian x coordinate
real(kind=wp),	intent(in)	::	y	Cartesian y coordinate
real(kind=wp),	intent(in)	::	z	Cartesian z coordinate
real(kind=wp),	intent(out)	::	phi	latitude (rad)
real(kind=wp),	intent(out)	::	lambda	longitude (rad)
real(kind=wp),	intent(out)	::	h	altitude
logical,	intent(out)	::	done	true if one of the special cases was computed

public subroutine direct_inverse_test()

Unit test for the direct and inverse geodetic routines.

Arguments

None

geodesy_module Module

Contents

Functions

Subroutines

Uses

Used by

Functions

public pure function great_circle_distance(r, long1, lat1, long2, lat2) result(d)

Arguments

Return Value real(kind=wp)

public pure function geocentric_radius(a, b, lat) result(r)

Arguments

Return Value real(kind=wp)

private pure function solve_polynomial(B, x0, error) result(x)

Arguments

Return Value real(kind=wp)

Subroutines

public pure subroutine heikkinen(rvec, a, b, h, lon, lat)

Arguments

public pure subroutine olson(rvec, a, b, h, long, lat)

Arguments

public subroutine direct(a, f, glat1, glon1, faz, s, glat2, glon2, baz)

Arguments

public subroutine geodetic_to_cartesian(a, b, glat, lon, h, r)

Arguments

public subroutine inverse(a, rf, b1, l1, b2, l2, faz, baz, s, it, sig, lam, kind)

Arguments

public subroutine geodetic_to_cartesian_triaxial(ax, ay, b, phi, lambda, h, r)

Arguments

public pure subroutine geodetic_to_cartesian_triaxial_2(a, b, c, lat, long, h, r)

Arguments

public subroutine cartesian_to_geodetic_triaxial(ax, ay, b, r, tol, phi, lambda, h)

Arguments

private subroutine bisection_special_2(cz, Xo, Yo, tol, n, m, Gm)

Arguments

private subroutine bisection_special_3(cx, cy, Xo, Yo, Zo, tol, n, m, Hm)

Arguments

private subroutine philambda_quadrant(x, y, z, phi, lambda)

Arguments

private subroutine xyz2philambda(ax, ay, b, x, y, z, phi, lambda)

Arguments

private subroutine xyz2fl(ax, ay, b, x, y, z, latitude, longitude)

Arguments

private pure subroutine horner(B, c, BB)

Arguments

public subroutine CartesianIntoGeodeticI(ax, ay, az, r, latitude, longitude, altitude, error)

Arguments

public subroutine CartesianIntoGeodeticII(ax, ay, az, r, latitude, longitude, altitude, error)

Arguments

public subroutine cartesian_to_geodetic_triaxial_2(a, b, c, r, eps, phi, lambda, h)

Arguments

private subroutine special_cases(a, b, c, x, y, z, phi, lambda, h, done)

Arguments

public subroutine direct_inverse_test()

Arguments