 |
|  |
 |
|
Asher
|
 |
Calgary, Alberta
Nov 1999 time: 22:37
|
|
http://www.csit.fsu.edu/~burkardt/f...lpak/polpak.f90
see the spherical harmonic function
code: subroutine spherical_harmonic ( l, m, theta, phi, c, s )
!*************************************************
******************************
!
!! SPHERICAL_HARMONIC evaluates spherical harmonic functions.
!
! Discussion:
!
! The spherical harmonic function Y(L,M,THETA,PHI)(X) is the
! angular part of the solution to Laplace's equation in spherical
! coordinates.
!
! Y(L,M,THETA,PHI)(X) is related to the associated Legendre
! function as follows:
!
! Y(L,M,THETA,PHI)(X) = FACTOR * P(L,M)(cos(THETA)) * exp ( i * M * PHI )
!
! Here, FACTOR is a normalization factor:
!
! FACTOR = sqrt ( ( 2 * L + 1 ) * ( L - 1 )! / ( 4 * PI * ( L + M )! ) )
!
! In Mathematica, a spherical harmonic function can be evaluated by
!
! SphericalHarmonicY [ l, m, theta, phi ]
!
! Note that notational tradition in physics requires that THETA
! and PHI represent the reverse of what they would normally mean
! in mathematical notation; that is, THETA goes up and down, and
! PHI goes around.
!
! Modified:
!
! 04 March 2005
!
! Author:
!
! John Burkardt
!
! Reference:
!
! Milton Abramowitz and Irene Stegun,
! Handbook of Mathematical Functions,
! US Department of Commerce, 1964.
!
! Eric Weisstein,
! CRC Concise Encyclopedia of Mathematics,
! CRC Press, 1999.
!
! Stephen Wolfram,
! The Mathematica Book,
! Fourth Edition,
! Wolfram Media / Cambridge University Press, 1999.
!
! Parameters:
!
! Input, integer L, the first index of the spherical harmonic function.
! Normally, 0 <= L.
!
! Input, integer M, the second index of the spherical harmonic function.
! Normally, -L <= M <= L.
!
! Input, real ( kind = 8 ) THETA, the polar angle, for which
! 0 <= THETA <= PI.
!
! Input, real ( kind = 8 ) PHI, the longitudinal angle, for which
! 0 <= PHI <= 2*PI.
!
! Output, real ( kind = 8 ) C(0:L), S(0:L), the real and imaginary
! parts of the functions Y(L,0:L,THETA,PHI).
!
implicit none
integer l
real ( kind = 8 ) c(0:l)
integer m
integer m_abs
real ( kind = 8 ) phi
real ( kind = 8 ) plm(0:l)
real ( kind = 8 ) s(0:l)
real ( kind = 8 ) theta
m_abs = abs ( m )
call legendre_associated_normalized ( l, m_abs, cos ( theta ), plm )
c(0:l) = plm(0:l) * cos ( real ( m, kind = 8 ) * phi )
s(0:l) = plm(0:l) * sin ( real ( m, kind = 8 ) * phi )
if ( m < 0 ) then
c(0:l) = -c(0:l)
s(0:l) = -s(0:l)
end if
return
end
subroutine spherical_harmonic_values ( n_data, l, m, theta, phi, yr, yi )
!*************************************************
******************************
!
!! SPHERICAL_HARMONIC_VALUES returns values of spherical harmonic functions.
!
! Discussion:
!
! In Mathematica, the function can be evaluated by
!
! SphericalHarmonicY [ l, m, theta, phi ]
!
! Modified:
!
! 05 March 2005
!
! Author:
!
! John Burkardt
!
! Reference:
!
! Milton Abramowitz and Irene Stegun,
! Handbook of Mathematical Functions,
! US Department of Commerce, 1964.
!
! Eric Weisstein,
! CRC Concise Encyclopedia of Mathematics,
! CRC Press, 1998.
!
! Stephen Wolfram,
! The Mathematica Book,
! Fourth Edition,
! Wolfram Media / Cambridge University Press, 1999.
!
! Parameters:
!
! Input/output, integer N_DATA.
! On input, if N_DATA is 0, the first test data is returned, and
! N_DATA is set to the index of the test data. On each subsequent
! call, N_DATA is incremented and that test data is returned. When
! there is no more test data, N_DATA is set to 0.
!
! Output, integer L, integer M, real ( kind = 8 ) THETA, PHI, the arguments
! of the function.
!
! Output, real ( kind = 8 ) YR, YI, the real and imaginary parts of
! the function.
!
implicit none
integer, parameter :: nmax = 20
integer l
integer, save, dimension ( nmax ) :: l_vec = (/ &
0, 1, 2, &
3, 4, 5, &
5, 5, 5, &
5, 4, 4, &
4, 4, 4, &
3, 3, 3, &
3, 3 /)
integer m
integer, save, dimension ( nmax ) :: m_vec = (/ &
0, 0, 1, &
2, 3, 5, &
4, 3, 2, &
1, 2, 2, &
2, 2, 2, &
-1, -1, -1, &
-1, -1 /)
integer n_data
real ( kind = 8 ) phi
real ( kind = 8 ), save, dimension ( nmax ) :: phi_vec = (/ &
0.1047197551196598D+01, 0.1047197551196598D+01, 0.1047197551196598D+01, &
0.1047197551196598D+01, 0.1047197551196598D+01, 0.6283185307179586D+00, &
0.6283185307179586D+00, 0.6283185307179586D+00, 0.6283185307179586D+00, &
0.6283185307179586D+00, 0.7853981633974483D+00, 0.7853981633974483D+00, &
0.7853981633974483D+00, 0.7853981633974483D+00, 0.7853981633974483D+00, &
0.4487989505128276D+00, 0.8975979010256552D+00, 0.1346396851538483D+01, &
0.1795195802051310D+01, 0.2243994752564138D+01 /)
real ( kind = 8 ) theta
real ( kind = 8 ), save, dimension ( nmax ) :: theta_vec = (/ &
0.5235987755982989D+00, 0.5235987755982989D+00, 0.5235987755982989D+00, &
0.5235987755982989D+00, 0.5235987755982989D+00, 0.2617993877991494D+00, &
0.2617993877991494D+00, 0.2617993877991494D+00, 0.2617993877991494D+00, &
0.2617993877991494D+00, 0.6283185307179586D+00, 0.1884955592153876D+01, &
0.3141592653589793D+01, 0.4398229715025711D+01, 0.5654866776461628D+01, &
0.3926990816987242D+00, 0.3926990816987242D+00, 0.3926990816987242D+00, &
0.3926990816987242D+00, 0.3926990816987242D+00 /)
real ( kind = 8 ) yi
real ( kind = 8 ), save, dimension ( nmax ) :: yi_vec = (/ &
0.0000000000000000D+00, 0.0000000000000000D+00, -0.2897056515173922D+00, &
0.1916222768312404D+00, 0.0000000000000000D+00, 0.0000000000000000D+00, &
0.3739289485283311D-02, -0.4219517552320796D-01, 0.1876264225575173D+00, &
-0.3029973424491321D+00, 0.4139385503112256D+00, -0.1003229830187463D+00, &
0.0000000000000000D+00, -0.1003229830187463D+00, 0.4139385503112256D+00, &
-0.1753512375142586D+00, -0.3159720118970196D+00, -0.3940106541811563D+00, &
-0.3940106541811563D+00, -0.3159720118970196D+00 /)
real ( kind = 8 ) yr
real ( kind = 8 ), save, dimension ( nmax ) :: yr_vec = (/ &
0.2820947917738781D+00, 0.4231421876608172D+00, -0.1672616358893223D+00, &
-0.1106331731112457D+00, 0.1354974113737760D+00, 0.5390423109043568D-03, &
-0.5146690442951909D-02, 0.1371004361349490D-01, 0.6096352022265540D-01, &
-0.4170400640977983D+00, 0.0000000000000000D+00, 0.0000000000000000D+00, &
0.0000000000000000D+00, 0.0000000000000000D+00, 0.0000000000000000D+00, &
0.3641205966137958D+00, 0.2519792711195075D+00, 0.8993036065704300D-01, &
-0.8993036065704300D-01, -0.2519792711195075D+00 /)
if ( n_data < 0 ) then
n_data = 0
end if
n_data = n_data + 1
if ( nmax < n_data ) then
n_data = 0
l = 0
m = 0
theta = 0.0D+00
phi = 0.0D+00
yr = 0.0D+00
yi = 0.0D+00
else
l = l_vec(n_data)
m = m_vec(n_data)
theta = theta_vec(n_data)
phi = phi_vec(n_data)
yr = yr_vec(n_data)
yi = yi_vec(n_data)
end if
return
end
|
|
|  |
All times are GMT. The time now is 05:37. Apolyton Time is 00:37. |
top of page
|
|
|
Forum Rules:
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts
|
HTML code is ON
vB code is ON
Smilies are ON
[IMG] code is ON
|
|
|
|
|
|