add meson, correct ecef2geodetic negative altitude

Author: scivision

CMakeLists.txt

@@ -10,10 +10,6 @@ if(NOT realbits)
   set(realbits 64)
 endif()
 
-if(NOT WIN32)
-  option(BUILD_SHARED_LIBS "Shared libraries .so. Feel free to turn off." ON)
-endif()
-
 # --- Maptran library
 add_library(maptran src/maptran.F90 src/vallado.F90)
 target_compile_options(maptran PRIVATE ${FFLAGS})

README.md

@@ -8,7 +8,7 @@ Similar to Python
 
 ## Install
 
-Requires Fortran 2003+ compiler, such as `gfortran`, `ifort`, PGI, `nagfor`, `flang`, etc.
+Requires Fortran 2008 compiler, such as `gfortran`, `ifort`, PGI, `nagfor`, `flang`, Cray, IBM XL, etc.
 
 ```sh
 cd bin

cmake/compilers.cmake

@@ -6,7 +6,7 @@ endif()
 
 if(CMAKE_Fortran_COMPILER_ID STREQUAL Intel)
 
-  list(APPEND FFLAGS -traceback -warn)
+  list(APPEND FFLAGS -traceback -warn -implicitnone)
   if(CMAKE_BUILD_TYPE STREQUAL Debug)
     list(APPEND FFLAGS -check all -debug extended)
   endif()
@@ -16,15 +16,14 @@ elseif(CMAKE_Fortran_COMPILER_ID STREQUAL GNU)
     list(APPEND FFLAGS -std=f2018)
   endif()
 
-  list(APPEND FFLAGS -march=native -Wall -Wextra -Wpedantic -Warray-bounds)# -Wfatal-errors)
+  list(APPEND FFLAGS -fimplicit-none -march=native -Wall -Wextra -Wpedantic -Warray-bounds)# -Wfatal-errors)
   if(CMAKE_BUILD_TYPE STREQUAL Debug)
-    list(APPEND FFLAGS -ffpe-trap=zero,overflow)
+    list(APPEND FFLAGS -ffpe-trap=overflow)
   endif()
 elseif(CMAKE_Fortran_COMPILER_ID STREQUAL PGI)
 
 elseif(CMAKE_Fortran_COMPILER_ID STREQUAL Flang)
   list(APPEND FFLAGS -Mallocatable=03)
-  list(APPEND FLIBS -static-flang-libs)
 elseif(CMAKE_Fortran_COMPILER_ID STREQUAL NAG)
   list(APPEND FFLAGS -f2008 -C -colour -gline -nan -info -u)
 endif()

meson.build

@@ -0,0 +1,19 @@
+project('MapTran', 'fortran', default_options : ['default_library=static'])
+
+REALBITS = '-DREALBITS='+get_option('realbits')
+
+fc = meson.get_compiler('fortran')
+if fc.get_id() == 'gcc'
+  add_global_arguments('-fimplicit-none', '-Wall', '-Wextra', '-Wpedantic', language : 'fortran')
+endif
+
+# --- Maptran library
+maptran = library('maptran', 'src/maptran.F90', 'src/vallado.F90',
+  fortran_args : REALBITS)
+  
+# --- testing
+mtexe = executable('testmaptran', 'tests/test_mod.f90', 'src/assert.F90',
+  link_with : maptran,
+  fortran_args : [REALBITS, '-DF08=1'])
+
+test('Maptran', mtexe)

meson_options.txt

@@ -0,0 +1 @@
+option('realbits', type : 'string', value : '64')

src/assert.F90

@@ -40,7 +40,10 @@ elemental logical function isclose(actual, desired, rtol, atol, equal_nan)
 
   real(wp) :: r,a
   logical :: n
-  ! this is appropriate INSTEAD OF merge(), since non present values aren't defined.
+  
+  isclose = .false. !< ensure it's defined
+  
+  !> INSTEAD OF merge(), since non present values aren't defined.
   r = 1e-5_wp
   a = 0._wp
   n = .false.
@@ -56,8 +59,10 @@ elemental logical function isclose(actual, desired, rtol, atol, equal_nan)
   !isclose = (actual == desired)
   !if (isclose) return
 !--- equal nan
-  isclose = n.and.(ieee_is_nan(actual).and.ieee_is_nan(desired))
-  if (isclose) return
+  if (n) then ! fortran is NOT short circuit logic in general
+    isclose = (ieee_is_nan(actual) .and. ieee_is_nan(desired)) 
+    if (isclose) return
+  endif
 !--- Inf /= -Inf, unequal NaN
   if (.not.ieee_is_finite(actual) .or. .not.ieee_is_finite(desired)) return
 !--- floating point closeness check

src/maptran.F90

@@ -58,11 +58,8 @@ impure elemental subroutine lookAtSpheroid(lat0, lon0, h0, az, tilt, lat, lon, s
 
 nan = ieee_value(0._wp, ieee_quiet_nan)
 
-if (present(spheroid)) then
-  ell = spheroid
-else
-  ell = wgs84Ellipsoid
-endif
+ell = wgs84Ellipsoid
+if (present(spheroid)) ell = spheroid
 
 d=.true.
 if (present(deg)) d = deg
@@ -119,13 +116,10 @@ elemental subroutine ecef2geodetic(x, y, z, lat, lon, alt, spheroid, deg)
 
   real(wp) :: ea, eb, r, E, u, Q, huE, Beta, eps
   type(Ellipsoid) :: ell
-  logical :: d
+  logical :: d, inside
 
-  if (present(spheroid)) then
-    ell = spheroid
-  else
-    ell = wgs84Ellipsoid
-  endif
+  ell = wgs84Ellipsoid
+  if (present(spheroid)) ell = spheroid
 
   ea = ell%SemimajorAxis
   eb = ell%SemiminorAxis
@@ -142,7 +136,7 @@ elemental subroutine ecef2geodetic(x, y, z, lat, lon, alt, spheroid, deg)
   huE = hypot(u, E)
 
   ! eqn. 4b
-  Beta = atan(huE / u * z / hypot(x, y))
+  Beta = atan(huE / u * z, hypot(x, y))
 
   ! eqn. 13
   eps = ((eb * u - ea * huE + E**2) * sin(Beta)) / (ea * huE * 1 / cos(Beta) - E**2 * cos(Beta))
@@ -151,11 +145,18 @@ elemental subroutine ecef2geodetic(x, y, z, lat, lon, alt, spheroid, deg)
 ! final output
   lat = atan(ea / eb * tan(Beta))
 
-  lon = atan(y/x)
+  lon = atan(y, x)
 
 ! eqn. 7
-  if (present(alt)) alt = sqrt((z - eb * sin(Beta))**2 + (Q - ea * cos(Beta))**2)
-
+  if (present(alt)) then
+    alt = hypot(z - eb * sin(Beta), Q - ea * cos(Beta))
+    
+    !> inside ellipsoid?
+    inside = x**2 / ea**2 + y**2 / ea**2 + z**2 / eb**2 < 1._wp
+    if (inside) alt = -alt
+  endif
+  
+  
   d=.true.
   if (present(deg)) d = deg
 
@@ -193,11 +194,8 @@ elemental subroutine geodetic2ecef(lat,lon,alt, x,y,z, spheroid, deg)
   d = .true.  ! not merge, Flang gives segfault
   if (present(deg)) d = deg
 
-  if (present(spheroid)) then
-    ell = spheroid
-  else
-    ell = wgs84Ellipsoid
-  endif 
+  ell = wgs84Ellipsoid
+  if (present(spheroid)) ell = spheroid
 
   if (d) then
     lat = radians(lat)
@@ -468,8 +466,8 @@ elemental subroutine enu2aer(east, north, up, az, elev, slantRange, deg)
   r = hypot(east, north)
   slantRange = hypot(r, up)
   ! radians
-  elev = atan2(up, r)
-  az = modulo(atan2(east, north), 2._wp * atan2(0._wp, -1._wp))
+  elev = atan(up, r)
+  az = modulo(atan(east, north), 2._wp * atan(0._wp, -1._wp))
 
   d=.true.
   if (present(deg)) d = deg

tests/test_mod.f90

@@ -7,7 +7,7 @@ program Test
 
 implicit none
 
-!type(Ellipsoid), parameter :: spheroid = wgs84Ellipsoid
+type(Ellipsoid), parameter :: spheroid = wgs84Ellipsoid
 
 real(wp), parameter :: lat = 42, lon= -82, alt = 200, &
                        az = 33, el=70, rng= 1e3_wp, &
@@ -45,8 +45,6 @@ program Test
 
 nan = ieee_value(0._wp, ieee_quiet_nan)
 
-!print*,'Default WGS84 Ellipsoid:',spheroid
-
 ! --------- scalar degrees      
 az5 = [0._wp, 10._wp, 125._wp];
 tilt = [30, 45, 90];
@@ -64,7 +62,7 @@ program Test
 call assert_allclose(rng5, [230.9413173_wp, 282.84715651_wp, nan], rtol=0.01_wp, equal_nan=.true.)
 
 
-!-------
+!! Scalar degrees
 
 call geodetic2ecef(lat,lon,alt, x1,y1,z1)
 call assert_allclose([x1,y1,z1],[x0,y0,z0], &
@@ -76,10 +74,23 @@ program Test
 call aer2ecef(az,el,rng,lat,lon,alt,x2,y2,z2)
 call assert_allclose([x2,y2,z2],[xl,yl,zl])
 
+!> ECEF2GEODETIC tests
 call ecef2geodetic(x1,y1,z1,lat2,lon2,alt2)
 call assert_allclose([lat2,lon2,alt2],[lat,lon,alt], &
                     rtol=0.01_wp, err_msg='ecef2geodetic-degrees')
 
+call ecef2geodetic(spheroid%SemimajorAxis-1._wp, 0._wp, 0._wp, lat2, lon2, alt2)
+call assert_allclose([lat2,lon2,alt2],[0._wp, 0._wp, -1._wp], &
+                    rtol=0.01_wp, err_msg='ecef2geodetic-degrees')
+
+call ecef2geodetic(0._wp, spheroid%SemimajorAxis-1._wp, 0._wp, lat2, lon2, alt2)
+call assert_allclose([lat2,lon2,alt2],[0._wp, 90._wp, -1._wp], &
+                    rtol=0.01_wp, err_msg='ecef2geodetic-degrees')
+
+call ecef2geodetic(0._wp, 0._wp, spheroid%SemiminorAxis-1._wp, lat2, lon2, alt2)
+call assert_allclose([lat2,lon2,alt2],[90._wp, 0._wp, -1._wp], &
+                    rtol=0.01_wp, err_msg='ecef2geodetic-degrees')
+
 call enu2aer(e1,n1,u1, az2, el2, rng2)
 call assert_allclose([az2,el2,rng2],[az,el,rng],err_msg='enu2aer-degrees')