/
TestExpandDH.f95
104 lines (81 loc) · 3.53 KB
/
TestExpandDH.f95
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
Program TestExpandDH
!------------------------------------------------------------------------------
!
! This program will expand a grid containing n samples in both longitude and
! latitude into spherical harmonics using the sampling theorems given in
! Driscoll and Heally (1994).
!
! Copyright (c) 2005, SHTOOLS
! All rights reserved.
!
!------------------------------------------------------------------------------
use SHTOOLS
use ftypes
implicit none
integer(int32), parameter :: maxdeg = 120, maxgrid = 800
character(240) :: infile
real(dp) :: cilm(2, maxdeg+1, maxdeg+1), grid(maxgrid, maxgrid), &
griddh(maxgrid, maxgrid), cilm2(2, maxdeg+1, maxdeg+1), &
interval, error, maxerror
integer(int32) :: lmax, n, nlat, nlong, lmax2, l, m, i
! Path to example data files may be passed as first argument, or use a default.
if (command_argument_count() > 0) then
call get_command_argument(1, infile)
else
infile = "../../ExampleDataFiles"
end if
infile = trim(infile) // "/MarsTopo719.shape"
call SHRead(infile, cilm, lmax)
print*, "Lmax of data file = ", lmax
!--------------------------------------------------------------------------
!
! Create an equally sampled grid
!
!--------------------------------------------------------------------------
lmax = 89
print*, "For this test of SHExpandDH, the initial grid will be created with lmax = ", lmax
interval = 1.0_dp
print*, "Making grid using MakeGrid2D..."
call MakeGrid2D(grid, cilm, lmax, interval, nlat, nlong)
print*, "nlat, nlong = ", nlat, nlong
print*, "Minimum and maximum values of the input grid = ", minval(grid(1:nlat, 1:nlong)), maxval(grid(1:nlat,1:nlong))
n = nlat -1
print*, "Number of samples in latitude used with SHExpandDH = ", n
griddh(1:n,1:2*n) = grid(1:n, 1:2*n)
!--------------------------------------------------------------------------
!
! Expand equally spaced grid using Driscoll and Heally routine, and
! compute relative error between input and output spherical harmonic
! coefficients.
!
!--------------------------------------------------------------------------
call SHExpandDH(griddh, n, cilm2, lmax2, sampling = 2)
print*, "Lmax of Driscoll Heally expansion = ", lmax2
maxerror = 0.0_dp
do l = 0, lmax
error = (cilm2(1,l+1,1)-cilm(1,l+1,1)) / cilm(1,l+1,1)
if (error > maxerror) maxerror = error
do m = 1, l, 1
error = (cilm2(1,l+1,m+1)-cilm(1,l+1,m+1)) / cilm(1,l+1,m+1)
if (error > maxerror) maxerror = error
error = (cilm2(2,l+1,m+1)-cilm(2,l+1,m+1)) / cilm(2,l+1,m+1)
if (error > maxerror) maxerror = error
end do
end do
print*, "Maximum relative error between initial and output spherical harmonic coefficients = ", maxerror
!--------------------------------------------------------------------------
!
! Re-expand coefficients in the space domain, and compare
! with input grid.
!
!--------------------------------------------------------------------------
griddh = 0.0_dp
call MakeGridDH(griddh, n, cilm2, lmax2)
print*, "Number of samples in latitude and longitude = ", n
maxerror = 0.0_dp
do i = 1, n
error = maxval( abs( griddh(1:n,i) - grid(1:n, 2*i-1) ) )
if (error > maxerror) maxerror = error
end do
print*, "Maximum error (meters) in the space domain = ", maxerror
end program TestExpandDH