ssht_sampling_mod.f90

Go to the documentation of this file.
00001 ! SSHT package to perform spin spherical harmonic transforms
00002 ! Copyright (C) 2011  Jason McEwen
00003 ! See LICENSE.txt for license details
00004 
00005 
00006 !------------------------------------------------------------------------------
00007 ! ssht_sampling_mod  -- SSHT library sampling class
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   ! Subroutine and function scope
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   ! Global variables
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   ! Data types
00069   !---------------------------------------
00070 
00071   ! None.
00072 
00073 
00074   !----------------------------------------------------------------------------
00075 
00076 contains
00077 
00078 
00079   !============================================================================
00080   ! Sampling weights
00081   !============================================================================
00082 
00083 
00084   !--------------------------------------------------------------------------
00085   ! ssht_sampling_weight_dh
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   ! ssht_sampling_weight_mw
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        ! Even case
00143        w = 2d0 / (1d0 - p**2)
00144     else
00145        ! Odd case (|p| /= 1)
00146        w = 0d0
00147     end if
00148 
00149   end function ssht_sampling_weight_mw
00150 
00151 
00152   !--------------------------------------------------------------------------
00153   ! ssht_sampling_gl_thetas_weights
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   ! gauleg
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   ! Sampling relations
00247   !============================================================================
00248 
00249 
00250   !----------------------------------------------------------------------------
00251   ! ssht_sampling_dh_t2theta
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   ! ssht_sampling_dh_p2phi
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   ! ssht_sampling_gl_p2phi
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   ! ssht_sampling_mw_t2theta
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   ! ssht_sampling_mw_p2phi
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   ! ssht_sampling_mweo_t2theta
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   ! ssht_sampling_mweo_p2phi
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   ! Harmonic index relations
00462   !============================================================================
00463 
00464   
00465   !----------------------------------------------------------------------------
00466   ! ssht_sampling_elm2ind
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   ! ssht_sampling_ind2elm
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
Generated on Mon Oct 31 01:20:05 2011 by  doxygen 1.6.3