-
Notifications
You must be signed in to change notification settings - Fork 103
/
SHFindLWin.f95
91 lines (74 loc) · 2.57 KB
/
SHFindLWin.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
function SHFindLWin(theta0, m, alpha, taper_number)
!------------------------------------------------------------------------------
!
! This function will output the spherical harmonic bandwidth for
! the space concentrated taper that yields a concetration factor
! equal to or greater than the input value of alpha. By default, the
! concentration factor of the first taper is computed, but this can be
! modified by inputing the optional arguement taper_number.
!
! Calling Parameters
!
! IN
! theta0 Angular radius of the concentration spherical
! cap in RADIANS.
! alpha Concentration factor of the window.
! m Angular order of the space concetrated tapers.
!
! IN (OPTIONAL)
! taper_number Taper number used to calculate concentration
! factors.
!
! Copyright (c) 2005-2019, SHTOOLS
! All rights reserved.
!
!------------------------------------------------------------------------------
use SHTOOLS, only: ComputeDm, EigvalSym
use ftypes
implicit none
integer(int32) :: SHFindLWin
real(dp), intent(in) :: theta0, alpha
integer(int32), intent(in) :: m
integer(int32), intent(in), optional :: taper_number
real(dp), allocatable :: dllm(:,:), eval(:)
real(dp) :: pi, alpha1
integer(int32) :: l, astat(2), tn
if (alpha < 0.0_dp .or. alpha > 1.0_dp) then
print*, "Error --- SHFindLWin"
print*, "The concentration factor alpha must be between 0 and 1."
print*, "Input value is ", alpha
stop
end if
pi = acos(-1.0_dp)
if (present(taper_number)) then
if (taper_number < 1) then
print*, "Error --- SHFindLWin"
print*, "TAPER_NUMBER must be greater than 0."
print*, "Input value is ", taper_number
stop
end if
tn = taper_number
else
tn = 1
end if
l = tn
do
l = l + 1
allocate (dllm(l+1,l+1), stat = astat(1)) ; dllm = 0.0_dp
allocate (eval(l+1), stat = astat(2)) ; eval = 0.0_dp
if (astat(1) /=0 .or. astat(2) /=0) then
print*, "Error --- SHFindLWin"
print*, "Probelm allocating arrays."
stop
end if
call ComputeDm(dllm, l, abs(m), theta0)
call EigValSym(dllm, l+1, eval(1:l+1))
alpha1 = eval(tn)
deallocate (dllm)
deallocate (eval)
if (alpha1 >= alpha) then
SHFindLWin = l
return
end if
end do
end function SHFindLWin