00001
00002
00003
00004
00005
00006
00007
00008
00017
00018 module ssht_sampling_mod
00019
00020 use ssht_types_mod
00021 use ssht_error_mod
00022
00023 implicit none
00024
00025 private
00026
00027
00028
00029
00030
00031
00032 public :: &
00033 ssht_sampling_elm2ind, &
00034 ssht_sampling_ind2elm, &
00035 ssht_sampling_dh_t2theta, &
00036 ssht_sampling_dh_p2phi, &
00037 ssht_sampling_gl_p2phi, &
00038 ssht_sampling_mw_t2theta, &
00039 ssht_sampling_mw_p2phi, &
00040 ssht_sampling_mweo_t2theta, &
00041 ssht_sampling_mweo_p2phi, &
00042 ssht_sampling_weight_dh, &
00043 ssht_sampling_weight_mw, &
00044 ssht_sampling_gl_thetas_weights
00045
00046
00047
00048
00049
00050
00052 integer, public, parameter :: SSHT_SAMPLING_DH = 1
00053
00055 integer, public, parameter :: SSHT_SAMPLING_GL = 2
00056
00058 integer, public, parameter :: SSHT_SAMPLING_MW = 3
00059
00061 integer, public, parameter :: SSHT_SAMPLING_MWEO = 4
00062
00064 integer, public, parameter :: SSHT_SAMPLING_DEFAULT = SSHT_SAMPLING_MW
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076 contains
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00099
00100 function ssht_sampling_weight_dh(theta_t, L) result(w)
00101
00102 real(dp), intent(in) :: theta_t
00103 integer, intent(in) :: L
00104 real(dp) :: w
00105
00106 integer :: k
00107
00108 w = 0d0
00109 do k = 0,L-1
00110 w = w + sin((2d0*k+1d0)*theta_t) / real(2d0*k+1d0,dp)
00111 end do
00112 w = (2d0/real(L,dp)) * sin(theta_t) * w
00113
00114 end function ssht_sampling_weight_dh
00115
00116
00117
00118
00119
00131
00132 function ssht_sampling_weight_mw(p) result(w)
00133
00134 integer, intent(in) :: p
00135 complex(dpc) :: w
00136
00137 if(p == 1) then
00138 w = I * PION2
00139 elseif(p == -1) then
00140 w = - I * PION2
00141 elseif(mod(p,2) == 0 ) then
00142
00143 w = 2d0 / (1d0 - p**2)
00144 else
00145
00146 w = 0d0
00147 end if
00148
00149 end function ssht_sampling_weight_mw
00150
00151
00152
00153
00154
00168
00169 subroutine ssht_sampling_gl_thetas_weights(thetas, weights, L)
00170
00171 integer, intent(in) :: L
00172 real(dp), intent(out) :: thetas(0:L-1)
00173 real(dp), intent(out) :: weights(0:L-1)
00174
00175 real(dp) :: znodes(0:L-1)
00176 integer :: p
00177
00178 call gauleg(-1d0, 1d0, znodes(0:L-1), weights(0:L-1), L)
00179
00180 do p = 0, L-1
00181 thetas(p) = acos(znodes(p))
00182 end do
00183
00184 end subroutine ssht_sampling_gl_thetas_weights
00185
00186
00187
00188
00189
00207
00208 subroutine gauleg(x1, x2, x, w, n)
00209
00210 integer, intent(in) :: n
00211 real(dp), intent(in) :: x1,x2
00212 real(dp), intent(out) :: x(1:n),w(1:n)
00213
00214 real(dp), parameter :: EPS = 1d-14
00215 integer :: i,j,m
00216 real(dp) :: p1,p2,p3,pp,xl,xm,z,z1
00217
00218 m=(n+1)/2
00219 xm=0.5d0*(x2+x1)
00220 xl=0.5d0*(x2-x1)
00221 do i=1,m
00222 z=cos(3.141592654d0*(i-.25d0)/(n+.5d0))
00223 do
00224 p1=1.d0
00225 p2=0.d0
00226 do j=1,n
00227 p3=p2
00228 p2=p1
00229 p1=((2.d0*j-1.d0)*z*p2-(j-1.d0)*p3)/j
00230 end do
00231 pp=n*(z*p1-p2)/(z*z-1.d0)
00232 z1=z
00233 z=z1-p1/pp
00234 if(abs(z-z1) <= EPS) exit
00235 end do
00236 x(i)=xm-xl*z
00237 x(n+1-i)=xm+xl*z
00238 w(i)=2.d0*xl/((1.d0-z*z)*pp*pp)
00239 w(n+1-i)=w(i)
00240 end do
00241
00242 end subroutine gauleg
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00268
00269 function ssht_sampling_dh_t2theta(t, L) result (theta)
00270
00271 integer, intent(in) :: t
00272 integer, intent(in) :: L
00273 real(dp) :: theta
00274
00275 theta = (2d0*t+1d0)*PI / (4d0*L)
00276
00277 end function ssht_sampling_dh_t2theta
00278
00279
00280
00281
00282
00298
00299 function ssht_sampling_dh_p2phi(p, L) result (phi)
00300
00301 integer, intent(in) :: p
00302 integer, intent(in) :: L
00303 real(dp) :: phi
00304
00305 phi = 2d0*p*PI / (2d0*L - 1d0)
00306
00307 end function ssht_sampling_dh_p2phi
00308
00309
00310
00311
00312
00328
00329 function ssht_sampling_gl_p2phi(p, L) result (phi)
00330
00331 integer, intent(in) :: p
00332 integer, intent(in) :: L
00333 real(dp) :: phi
00334
00335 phi = 2d0*p*PI / (2d0*L - 1d0)
00336
00337 end function ssht_sampling_gl_p2phi
00338
00339
00340
00341
00342
00358
00359 function ssht_sampling_mw_t2theta(t, L) result (theta)
00360
00361 integer, intent(in) :: t
00362 integer, intent(in) :: L
00363 real(dp) :: theta
00364
00365 theta = (2d0*t+1d0)*PI / (2d0*L - 1d0)
00366
00367 end function ssht_sampling_mw_t2theta
00368
00369
00370
00371
00372
00388
00389 function ssht_sampling_mw_p2phi(p, L) result (phi)
00390
00391 integer, intent(in) :: p
00392 integer, intent(in) :: L
00393 real(dp) :: phi
00394
00395 phi = 2d0*p*PI / (2d0*L - 1d0)
00396
00397 end function ssht_sampling_mw_p2phi
00398
00399
00400
00401
00402
00418
00419 function ssht_sampling_mweo_t2theta(t, L) result (theta)
00420
00421 integer, intent(in) :: t
00422 integer, intent(in) :: L
00423 real(dp) :: theta
00424
00425 theta = (2d0*t+1d0)*PI / (2d0*L - 1d0)
00426
00427 end function ssht_sampling_mweo_t2theta
00428
00429
00430
00431
00432
00448
00449 function ssht_sampling_mweo_p2phi(p, L) result (phi)
00450
00451 integer, intent(in) :: p
00452 integer, intent(in) :: L
00453 real(dp) :: phi
00454
00455 phi = (2d0*p+1d0)*PI / (2d0*L - 1d0)
00456
00457 end function ssht_sampling_mweo_p2phi
00458
00459
00460
00461
00462
00463
00464
00465
00466
00467
00485
00486 subroutine ssht_sampling_elm2ind(ind, el, m)
00487
00488 integer, intent(out) :: ind
00489 integer, intent(in) :: el, m
00490
00491 ind = el**2 + el + m
00492
00493 end subroutine ssht_sampling_elm2ind
00494
00495
00496
00497
00498
00516
00517 subroutine ssht_sampling_ind2elm(el, m, ind)
00518
00519 integer, intent(out) :: el, m
00520 integer, intent(in) :: ind
00521
00522 el = floor(sqrt(real(ind,dp)))
00523 m = ind - el**2 - el
00524
00525 end subroutine ssht_sampling_ind2elm
00526
00527
00528 end module ssht_sampling_mod