made j2000_to_frame a function

also added frame_to_j2000

Author: jacobwilliams

src/moon_frame_module.F90

@@ -34,6 +34,7 @@ module moon_frame_module
         procedure,public :: initialize => initialize_moon_frame_interpolater
         procedure,public :: destroy    => destroy_moon_frame_interpolater
         procedure,public :: j2000_to_frame
+        procedure,public :: frame_to_j2000
     end type moon_frame_interpolater
 
 contains
@@ -96,10 +97,11 @@ subroutine initialize_moon_frame_interpolater(me, filename, k, extrapolate)
 
     end subroutine initialize_moon_frame_interpolater
 
-    subroutine j2000_to_frame(me, et, rot)
+    function j2000_to_frame(me, et) result(rot)
+        !! rotation matrix from J2000 to the frame
         class(moon_frame_interpolater), intent(inout) :: me
         real(wp), intent(in) :: et
-        real(wp), intent(out) :: rot(3, 3)
+        real(wp) :: rot(3, 3)
 
         real(wp) :: roll, pitch, yaw_x, yaw_y, yaw
         integer :: iflag
@@ -113,7 +115,15 @@ subroutine j2000_to_frame(me, et, rot)
         ! convert to rotation matrix:
         call rpy_to_rot(roll, pitch, yaw, rot)
 
-    end subroutine j2000_to_frame
+    end function j2000_to_frame
+
+    function frame_to_j2000(me, et) result(rot)
+        !! rotation matrix from the frame to j2000
+        class(moon_frame_interpolater), intent(inout) :: me
+        real(wp), intent(in) :: et
+        real(wp) :: rot(3, 3)
+        rot = transpose(me%j2000_to_frame(et))
+    end function frame_to_j2000
 
     pure subroutine rpy_to_rot(roll, pitch, yaw, r)
         !! roll, patch, yaw to rotation matrix

test/interp_test.f90

@@ -11,7 +11,7 @@ program interp_test
     call moon_pa%initialize('data/moon_pa_2000_2100.csv')
 
     et = 0.0_wp
-    call moon_pa%j2000_to_frame(et, rot)
+    rot = moon_pa%j2000_to_frame(et)
 
     write(*,*) ''
     write(*,*) 'rot at et=0.0:'

test/test.f90

@@ -32,14 +32,14 @@ program test
     r(:,1) = 1737.4_wp * unit([1000.0_wp, 1000.0_wp, 1000.0_wp])  ! point on surface of moon
 
     et = 0.0_wp
-    call moon_pa%j2000_to_frame(et, rot)
+    rot = moon_pa%j2000_to_frame(et)
     call from_j2000_to_moon_pa(real(et, real64), rot_true)
     write(*,*) ''
     write(*,*) 'rot error at et=0', rot_true - rot
     write(*,*) 'pos error: ', norm2(matmul(rot_true, r) - matmul(rot, r))
 
     et = 21600.0_wp ! halfway between first two points
-    call moon_pa%j2000_to_frame(et, rot)
+    rot = moon_pa%j2000_to_frame(et)
     call from_j2000_to_moon_pa(real(et, real64), rot_true)
     write(*,*) ''
     write(*,*) 'rot error at et=21600', rot_true - rot
@@ -49,7 +49,7 @@ program test
         ! quintic: -1.0719543411141785E-004   8.4190997773703202E-004  -1.9197702022211161E-004
 
     et = (129600.0_wp + 86400.0_wp) / 2.0_wp
-    call moon_pa%j2000_to_frame(et, rot)
+    rot = moon_pa%j2000_to_frame(et)
     call from_j2000_to_moon_pa(real(et, real64), rot_true)
     write(*,*) ''
     write(*,*) 'rot error at et', rot_true - rot
@@ -66,7 +66,7 @@ program test
     do
         et = et + 6.0_wp *3600.0_wp
         if (et > 3187296000.0_wp) exit
-        call moon_pa%j2000_to_frame(et, rot)          ! splined version
+        rot = moon_pa%j2000_to_frame(et)  ! splined version
         call from_j2000_to_moon_pa(real(et, real64), rot_true)      ! true version from spice : max error: 5.5664815594496319E-004 km
         ! call from_j2000_to_iau_moon(et, rot_true)   ! compare to iau_moon :     max error: 0.98148100150030848 km
         r_error = norm2(matmul(rot_true, r) - matmul(rot, r))
@@ -94,7 +94,7 @@ program test
     do
         et = et + 6.0_wp *3600.0_wp
         if (et > 3187296000.0_wp) exit
-        call moon_pa%j2000_to_frame(et, rot)
+        rot = moon_pa%j2000_to_frame(et)
     end do
     call cpu_time(tend)
     write(*,'(A,F6.3,A)') 'spline time: ', (tend-tstart), 'sec'