; $Id: spher_harm.pro,v 1.10 2004/01/21 15:55:02 scottm Exp $
;
; Copyright (c) 2000-2004, Research Systems, Inc.  All rights reserved.
;	Unauthorized reproduction prohibited.
;+
; NAME:
;	SPHER_HARM
;
; PURPOSE:
;	This function returns the value of the spherical harmonic
;   Y(L,M)[theta,phi], |M| <= L, which is a function of two coordinates
;   (theta, phi) on a spherical surface.
;
; CATEGORY:
;	Special Math Functions
;
; CALLING SEQUENCE:
;
;	Result = SPHER_HARM(Theta, Phi, L, M [, /DOUBLE])
;
; INPUTS:
;	Theta: The value of the polar (colatitudinal) coordinate at which Y(L,M)
;       is evaluated. Theta can be either a scalar or an array.
;
;   Phi: The value of the azimuthal (longitudinal) coordinate at which Y(L,M)
;       is evaluated. Phi can be either a scalar or an array.
;
;       If Theta and Phi are both arrays then they must have the same number
;       of elements.
;
;   L: An integer specifying the order L of Y(L,M). If L is of type float
;      then it will be truncated.
;
;   M: An integer, -L <= M <= L, specifying the azimuthal order M of Y(L,M).
;      If M is of type float then it will be truncated.
;
; KEYWORD PARAMETERS:
;	DOUBLE:	Set this keyword to force the computation to be done in
;       double-precision arithmetic.
;
; OUTPUT:
;	The result returned by SPHER_HARM is a complex array that has the
;   same dimensions as the input arrays.
;
; PROCEDURE:
;	Uses LEGENDRE().
;
; EXAMPLE:
;	To visualize the electron probability density for the hydrogen atom
;   in state 3d0. (Feynman, Leighton, and Sands, 1965:
;   The Feynman Lectures on Physics, Calif. Inst. Tech, Ch. 19).
;
; ;define a data cube (N x N x N)
;   n = 41L
;   a = 60*FINDGEN(n)/(n-1) - 29.999  ; [-1,+1]
;   x = REBIN(a, n, n, n)              ; X-coordinates of cube
;   y = REBIN(REFORM(a,1,n), n, n, n)  ; Y-coordinates
;   z = REBIN(REFORM(a,1,1,n), n, n, n); Z-coordinates
;
; ;convert from rectangular (x,y,z) to spherical (phi, theta, r)
;   spherCoord = CV_COORD(FROM_RECT=TRANSPOSE([[x[*]],[y[*]],[z[*]]]), /TO_SPHERE)
;   phi = REFORM(spherCoord[0,*], n, n, n)
;   theta = REFORM(!PI/2 - spherCoord[1,*], n, n, n)
;   r = REFORM(spherCoord[2,*], n, n, n)
;
; ;find electron probability density for hydrogen atom in state 3d0
;
; ;Angular component
;   L = 2   ; state "d" is electron spin L=2
;   M = 0   ; Z-component of spin is zero
;   angularState = SPHER_HARM(theta, phi, L, M)
;
; ;Radial component for state n=3, L=2
;   radialFunction = EXP(-r/2)*(r^2)
;
;   waveFunction = angularState*radialFunction
;   probabilityDensity = ABS(waveFunction)^2
;
;   SHADE_VOLUME, probabilityDensity, 0.1*MEAN(probabilityDensity), vertex, poly
;   oPolygon = OBJ_NEW('IDLgrPolygon', vertex, POLYGON=poly, COLOR=[180,180,180])
;   XOBJVIEW, oPolygon
;
; MODIFICATION HISTORY:
; 	Written by:	CT, RSI, March 2000.
;
;-

FUNCTION spher_harm, theta, phi, Linput, Minput, $
	DOUBLE=double

	COMPILE_OPT strictarr

; error checking
	ON_ERROR, 2
	IF (N_PARAMS() NE 4) THEN MESSAGE, 'Incorrect number of arguments.'
	IF ((N_ELEMENTS(Linput) NE 1) OR (N_ELEMENTS(Minput) NE 1)) THEN $
		MESSAGE,'L and M must be scalars.'
	L = LONG(Linput[0]) ; convert 1-element array to scalar
	M1 = LONG(Minput[0]) ; convert 1-element array to scalar
	IF (L LT 0) THEN MESSAGE, $
		'Argument L must be greater than or equal to zero.'

	M = ABS(M1)
	IF (M GT L) THEN MESSAGE, 'Argument M must be in the range [-L, L].'

	IF (SIZE(theta,/N_DIM) GT 0) AND (SIZE(phi,/N_DIM) GT 0) THEN $
		IF (N_ELEMENTS(theta) NE N_ELEMENTS(phi)) THEN MESSAGE, $
		'Theta or Phi must be scalar, or have the same number of values.'

; keyword checking
	thetaDouble = (SIZE(theta,/TNAME) EQ 'DOUBLE')
	phiDouble = (SIZE(phi,/TNAME) EQ 'DOUBLE')
	double = (N_ELEMENTS(double) GT 0) ? KEYWORD_SET(double) : $
		(thetaDouble OR phiDouble)

; normalizing factor (always double precision because of FACTORIAL)
	spherHarm = (M1 LT 0) ? (-1)^M : 1   ; handle negative M
	spherHarm = spherHarm*SQRT((2L*L + 1)/(4*!DPI)* $
		FACTORIAL(L - M)/FACTORIAL(L + M))

; define constants 1, I for single or double precision
	one = 1D
	I = DCOMPLEX(0,1)

	IF (NOT double) THEN BEGIN
		one = 1.0
		I = COMPLEX(0,1)
		spherHarm = FLOAT(spherHarm)  ; convert back
	ENDIF

	IF (M NE M1) THEN I = -I  ; complex conjugate for -M

; Polar coordinate
	cosTheta = (double) ? COS(DOUBLE(theta)) : COS(FLOAT(theta))
	Plm = LEGENDRE(TEMPORARY(cosTheta), L, M, DOUBLE=double)  ; P(L,M)
	spherHarm = spherHarm*TEMPORARY(Plm)

; Azimuthal coordinate
	mPhi = (double) ? M*DOUBLE(phi) : M*FLOAT(phi)

; Spherical harmonic
	realPart = spherHarm*COS(Mphi)
	imagPart = TEMPORARY(spherHarm)*SIN(TEMPORARY(Mphi))
	IF (M1 LT 0) THEN imagPart = -TEMPORARY(imagPart)  ; complex conjugate

	RETURN, COMPLEX(TEMPORARY(realPart), TEMPORARY(imagPart), DOUBLE=double)

END