|
CP2K 2.4 (Revision 12889)
|
00001 !-----------------------------------------------------------------------------! 00002 ! CP2K: A general program to perform molecular dynamics simulations ! 00003 ! Copyright (C) 2000 - 2013 CP2K developers group ! 00004 !-----------------------------------------------------------------------------! 00005 00006 ! ***************************************************************************** 00023 MODULE spherical_harmonics 00024 00025 USE f77_blas 00026 USE kinds, ONLY: dp,& 00027 dp_size 00028 USE mathconstants, ONLY: dfac,& 00029 fac,& 00030 maxfac,& 00031 pi 00032 USE termination, ONLY: stop_memory,& 00033 stop_program 00034 USE timings, ONLY: timeset,& 00035 timestop 00036 #include "cp_common_uses.h" 00037 00038 IMPLICIT NONE 00039 00040 PRIVATE 00041 00042 PUBLIC :: y_lm 00043 PUBLIC :: clebsch_gordon_deallocate, clebsch_gordon_init,& 00044 clebsch_gordon, clebsch_gordon_getm 00045 PUBLIC :: dPof1 , drvy_lm 00046 PUBLIC :: legendre, dlegendre 00047 00048 INTERFACE y_lm 00049 MODULE PROCEDURE rvy_lm, rry_lm, cvy_lm, ccy_lm 00050 END INTERFACE 00051 00052 INTERFACE clebsch_gordon 00053 MODULE PROCEDURE clebsch_gordon_real, clebsch_gordon_complex 00054 END INTERFACE 00055 00056 INTERFACE clebsch_gordon_getm 00057 MODULE PROCEDURE getm 00058 END INTERFACE 00059 00060 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'spherical_harmonics' 00061 00062 REAL(KIND=dp), DIMENSION(:,:,:), ALLOCATABLE :: cg_table 00063 INTEGER :: lmax=-1 00064 REAL(KIND=dp) :: osq2, sq2 00065 00066 CONTAINS 00067 00068 ! Clebsch-Gordon Coefficients 00069 00070 ! ***************************************************************************** 00071 SUBROUTINE clebsch_gordon_complex (l1,m1,l2,m2,clm) 00072 INTEGER, INTENT(IN) :: l1, m1, l2, m2 00073 REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: clm 00074 00075 CHARACTER(len=*), PARAMETER :: routineN = 'clebsch_gordon_complex', 00076 routineP = moduleN//':'//routineN 00077 00078 INTEGER :: icase, ind, l, lm, lp, n 00079 00080 l = l1 + l2 00081 IF (l > lmax) CALL clebsch_gordon_init ( l ) 00082 n = l/2 + 1 00083 IF (n > SIZE(clm)) CALL stop_program(routineN,moduleN,__LINE__,& 00084 "Array too small") 00085 00086 IF ( (m1 >= 0 .AND. m2 >= 0) .OR. (m1 < 0 .AND. m2 < 0) ) THEN 00087 icase = 1 00088 ELSE 00089 icase = 2 00090 END IF 00091 ind = order(l1,m1,l2,m2) 00092 00093 DO lp = MOD(l,2),l,2 00094 lm = lp/2+1 00095 clm(lm) = cg_table(ind,lm,icase) 00096 END DO 00097 00098 END SUBROUTINE clebsch_gordon_complex 00099 00100 ! ***************************************************************************** 00101 SUBROUTINE clebsch_gordon_real (l1,m1,l2,m2,rlm) 00102 INTEGER, INTENT(IN) :: l1, m1, l2, m2 00103 REAL(KIND=dp), DIMENSION(:, :), 00104 INTENT(OUT) :: rlm 00105 00106 CHARACTER(len=*), PARAMETER :: routineN = 'clebsch_gordon_real', 00107 routineP = moduleN//':'//routineN 00108 00109 INTEGER :: icase1, icase2, ind, l, lm, 00110 lp, mm(2), n 00111 REAL(KIND=dp) :: xsi 00112 00113 l = l1 + l2 00114 IF (l > lmax) CALL clebsch_gordon_init ( l ) 00115 n = l/2 + 1 00116 IF (n > SIZE(rlm,1)) CALL stop_program(routineN,moduleN,__LINE__,& 00117 "Array too small") 00118 00119 ind = order(l1,m1,l2,m2) 00120 mm = getm(m1,m2) 00121 IF ( (m1 >= 0 .AND. m2 >= 0) .OR. (m1 < 0 .AND. m2 < 0) ) THEN 00122 icase1 = 1 00123 icase2 = 2 00124 ELSE 00125 icase1 = 2 00126 icase2 = 1 00127 END IF 00128 00129 DO lp = MOD(l,2),l,2 00130 lm = lp/2+1 00131 xsi = get_factor(m1,m2,mm(1)) 00132 rlm(lm,1) = xsi*cg_table(ind,lm,icase1) 00133 xsi = get_factor(m1,m2,mm(2)) 00134 rlm(lm,2) = xsi*cg_table(ind,lm,icase2) 00135 END DO 00136 00137 END SUBROUTINE clebsch_gordon_real 00138 00139 ! ***************************************************************************** 00140 FUNCTION getm(m1,m2) RESULT(m) 00141 INTEGER, INTENT(IN) :: m1, m2 00142 INTEGER, DIMENSION(2) :: m 00143 00144 INTEGER :: mm, mp 00145 00146 mp = m1 + m2 00147 mm = m1 - m2 00148 IF ( m1*m2 < 0 .OR. (m1*m2==0 .AND. (m1<0 .OR. m2<0))) THEN 00149 mp = -ABS(mp) 00150 mm = -ABS(mm) 00151 ELSE 00152 mp = ABS(mp) 00153 mm = ABS(mm) 00154 END IF 00155 m(1) = mp 00156 m(2) = mm 00157 END FUNCTION getm 00158 00159 ! ***************************************************************************** 00160 FUNCTION get_factor(m1,m2,m) RESULT(f) 00161 INTEGER, INTENT(IN) :: m1, m2, m 00162 REAL(KIND=dp) :: f 00163 00164 INTEGER :: mx, my 00165 00166 f = 1.0_dp 00167 IF ( ABS(m1)>=ABS(m2) ) THEN 00168 mx = m1 00169 my = m2 00170 ELSE 00171 mx = m2 00172 my = m1 00173 ENDIF 00174 IF ( mx*my == 0 ) THEN 00175 f = 1.0_dp 00176 ELSE IF ( m == 0 ) THEN 00177 IF(ABS(mx) /= ABS(my)) WRITE(6,'(A,3I6)') " 1) Illegal Case ",m1,m2,m 00178 IF(mx*my > 0) THEN 00179 f = 1.0_dp 00180 ELSE 00181 f = 0.0_dp 00182 END IF 00183 ELSE IF ( ABS(mx)+ABS(my) == m ) THEN 00184 f = osq2 00185 IF(mx<0) f = -osq2 00186 ELSE IF ( ABS(mx)+ABS(my) == -m ) THEN 00187 f = osq2 00188 ELSE IF ( ABS(mx)-ABS(my) == -m ) THEN 00189 IF(mx*my > 0) WRITE(6,'(A,3I6)') " 2) Illegal Case ",m1,m2,m 00190 IF( mx > 0 ) f = -osq2 00191 IF( mx < 0 ) f = osq2 00192 ELSE IF ( ABS(mx)-ABS(my) == m ) THEN 00193 IF(mx*my < 0) WRITE(6,'(A,3I6)') " 3) Illegal Case ",m1,m2,m 00194 f = osq2 00195 ELSE 00196 WRITE(6,'(A,3I6)') " 4) Illegal Case ",m1,m2,m 00197 END IF 00198 END FUNCTION get_factor 00199 00200 ! ***************************************************************************** 00201 SUBROUTINE clebsch_gordon_deallocate() 00202 CHARACTER(len=*), PARAMETER :: routineN = 'clebsch_gordon_deallocate', 00203 routineP = moduleN//':'//routineN 00204 00205 INTEGER :: handle, ierr 00206 00207 CALL timeset(routineN,handle) 00208 00209 00210 IF ( ALLOCATED ( cg_table ) ) THEN 00211 DEALLOCATE ( cg_table, STAT=ierr ) 00212 IF (ierr /= 0) CALL stop_memory(routineN,moduleN,__LINE__,"cg_table") 00213 END IF 00214 CALL timestop(handle) 00215 00216 END SUBROUTINE clebsch_gordon_deallocate 00217 00218 ! ***************************************************************************** 00219 SUBROUTINE clebsch_gordon_init ( l ) 00220 INTEGER, INTENT(IN) :: l 00221 00222 CHARACTER(len=*), PARAMETER :: routineN = 'clebsch_gordon_init', 00223 routineP = moduleN//':'//routineN 00224 00225 INTEGER :: handle, i1, i2, ierr, ix, iy, 00226 l1, l2, lp, m, m1, m2, ml, 00227 mp, n 00228 00229 CALL timeset(routineN,handle) 00230 00231 sq2 = SQRT(2.0_dp) 00232 osq2 = 1.0_dp/sq2 00233 00234 IF ( l < 0 ) CALL stop_program(routineN,moduleN,__LINE__,"l < 0") 00235 IF ( ALLOCATED ( cg_table ) ) THEN 00236 DEALLOCATE ( cg_table, STAT=ierr ) 00237 IF (ierr /= 0) CALL stop_memory(routineN,moduleN,__LINE__,"cg_table") 00238 END IF 00239 ! maximum size of table 00240 n = (l**4 + 6*l**3 + 15*l**2 + 18*l + 8)/8 00241 m = l+1 00242 ALLOCATE ( cg_table(n,m,2), STAT=ierr ) 00243 IF (ierr /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00244 "cg_table",dp_size*2*n*m) 00245 lmax = l 00246 00247 DO l1=0,lmax 00248 DO m1=0,l1 00249 iy = (l1*(l1+1))/2 + m1 + 1 00250 DO l2=l1,lmax 00251 ml=0 00252 IF ( l1==l2 ) ml = m1 00253 DO m2=ml,l2 00254 ix = (l2*(l2+1))/2 + m2 + 1 00255 i1 = (ix*(ix-1))/2 + iy 00256 DO lp=MOD(l1+l2,2),l1+l2,2 00257 i2 = lp/2 + 1 00258 mp = m2 + m1 00259 cg_table(i1,i2,1) = cgc(l1,m1,l2,m2,lp,mp) 00260 mp = ABS ( m2 - m1 ) 00261 IF ( m2 >= m1 ) THEN 00262 cg_table(i1,i2,2) = cgc(l1,m1,lp,mp,l2,m2) 00263 ELSE 00264 cg_table(i1,i2,2) = cgc(l2,m2,lp,mp,l1,m1) 00265 END IF 00266 END DO 00267 END DO 00268 END DO 00269 END DO 00270 END DO 00271 00272 CALL timestop(handle) 00273 00274 END SUBROUTINE clebsch_gordon_init 00275 00276 ! ***************************************************************************** 00277 FUNCTION cgc(l1,m1,l2,m2,lp,mp) 00278 INTEGER, INTENT(IN) :: l1, m1, l2, m2, lp, mp 00279 REAL(KIND=dp) :: cgc 00280 00281 CHARACTER(len=*), PARAMETER :: routineN = 'cgc', 00282 routineP = moduleN//':'//routineN 00283 00284 INTEGER :: la, lb, ll, ma, mb, s, t, 00285 tmax, tmin, z1, z2, z3, z4, z5 00286 REAL(KIND=dp) :: f1, f2, pref 00287 00288 ! m1 >= 0; m2 >= 0; mp >= 0 00289 00290 IF ( m1<0 .OR. m2<0 .OR. mp<0 ) THEN 00291 WRITE(6,*) l1,l2,lp 00292 WRITE(6,*) m1,m2,mp 00293 CALL stop_program(routineN,moduleN,__LINE__,"Illegal input values") 00294 END IF 00295 IF ( l2 < l1 ) THEN 00296 la = l2 00297 ma = m2 00298 lb = l1 00299 mb = m1 00300 ELSE 00301 la = l1 00302 ma = m1 00303 lb = l2 00304 mb = m2 00305 END IF 00306 00307 IF ( MOD(la+lb+lp,2) == 0 .AND. la+lb >= lp .AND. lp >= lb-la & 00308 .AND. lb-mb >= 0 ) THEN 00309 ll = (2*lp+1) * (2*la+1) * (2*lb+1) 00310 pref = 1.0_dp/SQRT(4.0_dp*pi) * 0.5_dp * SQRT ( REAL(ll,dp) * 00311 (sfac(lp-mp)/sfac(lp+mp)) * 00312 (sfac(la-ma)/sfac(la+ma)) * (sfac(lb-mb)/sfac(lb+mb)) ) 00313 s = (la+lb+lp)/2 00314 tmin = MAX ( 0, -lb+la-mp ) 00315 tmax = MIN ( lb+la-mp, lp-mp, la-ma ) 00316 f1 = REAL( 2 * (-1)**(s-lb-ma),KIND=dp) * (sfac(lb+mb)/sfac(lb-mb)) * 00317 sfac(la+ma)/(sfac(s-lp)*sfac(s-lb)) * sfac(2*s-2*la)/sfac(s-la) * 00318 (sfac(s)/sfac(2*s+1)) 00319 f2 = 0.0_dp 00320 DO t=tmin,tmax 00321 z1 = lp+mp+t 00322 z2 = la+lb-mp-t 00323 z3 = lp-mp-t 00324 z4 = lb-la+mp+t 00325 z5 = la-ma-t 00326 f2 = f2 + (-1)**t * (sfac(z1)/(sfac(t)*sfac(z3)))* (sfac(z2)/(sfac(z4)*sfac(z5))) 00327 END DO 00328 cgc = pref * f1 * f2 00329 ELSE 00330 cgc = 0.0_dp 00331 END IF 00332 00333 END FUNCTION cgc 00334 00335 ! ***************************************************************************** 00336 FUNCTION sfac ( n ) RESULT (fval) 00337 INTEGER :: n 00338 REAL(KIND=dp) :: fval 00339 00340 INTEGER :: i 00341 00342 IF ( n > maxfac ) THEN 00343 fval = fac(maxfac) 00344 DO i = maxfac+1, n 00345 fval = REAL(i,dp)*fval 00346 END DO 00347 ELSE IF ( n >= 0 ) THEN 00348 fval = fac(n) 00349 ELSE 00350 fval = 1.0_dp 00351 END IF 00352 END FUNCTION sfac 00353 00354 ! ***************************************************************************** 00355 FUNCTION order(l1,m1,l2,m2) RESULT(ind) 00356 INTEGER, INTENT(IN) :: l1, m1, l2, m2 00357 INTEGER :: ind 00358 00359 INTEGER :: i1, i2, ix, iy 00360 00361 i1 = (l1*(l1+1))/2 + ABS(m1) + 1 00362 i2 = (l2*(l2+1))/2 + ABS(m2) + 1 00363 ix = MAX(i1,i2) 00364 iy = MIN(i1,i2) 00365 ind = (ix*(ix-1))/2+iy 00366 END FUNCTION order 00367 00368 ! Calculation of Spherical Harmonics 00369 00370 ! ***************************************************************************** 00371 SUBROUTINE rvy_lm ( r, y, l, m ) 00372 ! 00373 ! Real Spherical Harmonics 00374 ! _ _ 00375 ! | [(2l+1)(l-|m|)!] |1/2 m cos(m p) m>=0 00376 ! Y_lm ( t, p ) = |---------------------| P_l(cos(t)) 00377 ! |[2Pi(1+d_m0)(l+|m|)!]| sin(|m| p) m<0 00378 ! 00379 ! Input: r == (x,y,z) : normalised x^2 + y^2 + z^2 = 1 00380 ! 00381 REAL(KIND=dp), DIMENSION(:, :), 00382 INTENT(IN) :: r 00383 REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: y 00384 INTEGER, INTENT(IN) :: l, m 00385 00386 CHARACTER(len=*), PARAMETER :: routineN = 'rvy_lm', 00387 routineP = moduleN//':'//routineN 00388 00389 INTEGER :: i 00390 REAL(KIND=dp) :: cp, lmm, lpm, pf, plm, rxy, 00391 sp, t, z 00392 00393 SELECT CASE ( l ) 00394 CASE ( : -1 ) 00395 CALL stop_program(routineN,moduleN,__LINE__,"Negative l value") 00396 CASE ( 0 ) 00397 pf = SQRT ( 1.0_dp / ( 4.0_dp * pi ) ) 00398 IF ( m /= 0 ) CALL stop_program(routineN,moduleN,__LINE__,& 00399 "l = 0 and m value out of bounds") 00400 y(:) = pf 00401 CASE ( 1 ) 00402 SELECT CASE ( m ) 00403 CASE DEFAULT 00404 CALL stop_program(routineN,moduleN,__LINE__,& 00405 "l = 1 and m value out of bounds") 00406 CASE ( 1 ) 00407 pf = SQRT ( 3.0_dp / ( 4.0_dp * pi ) ) 00408 y ( : ) = pf * r ( 1, : ) 00409 CASE ( 0 ) 00410 pf = SQRT ( 3.0_dp / ( 4.0_dp * pi ) ) 00411 y ( : ) = pf * r ( 3, : ) 00412 CASE ( -1 ) 00413 pf = SQRT ( 3.0_dp / ( 4.0_dp * pi ) ) 00414 y ( : ) = pf * r ( 2, : ) 00415 END SELECT 00416 CASE ( 2 ) 00417 SELECT CASE ( m ) 00418 CASE DEFAULT 00419 CALL stop_program(routineN,moduleN,__LINE__,& 00420 "l = 2 and m value out of bounds") 00421 CASE ( 2 ) 00422 pf = SQRT ( 15.0_dp / ( 16.0_dp * pi ) ) 00423 y ( : ) = pf * ( r ( 1, : ) * r ( 1, : ) - r ( 2, : ) * r ( 2, : ) ) 00424 CASE ( 1 ) 00425 pf = SQRT ( 15.0_dp / ( 4.0_dp * pi ) ) 00426 y ( : ) = pf * r ( 3, : ) * r ( 1, : ) 00427 CASE ( 0 ) 00428 pf = SQRT ( 5.0_dp / ( 16.0_dp * pi ) ) 00429 y ( : ) = pf * ( 3.0_dp * r ( 3, : ) * r ( 3, : ) - 1.0_dp ) 00430 CASE ( -1 ) 00431 pf = SQRT ( 15.0_dp / ( 4.0_dp * pi ) ) 00432 y ( : ) = pf * r ( 3, : ) * r ( 2, : ) 00433 CASE ( -2 ) 00434 pf = SQRT ( 15.0_dp / ( 16.0_dp * pi ) ) 00435 y ( : ) = pf * 2.0_dp * r ( 1, : ) * r ( 2, : ) 00436 END SELECT 00437 CASE ( 3 ) 00438 SELECT CASE ( m ) 00439 CASE DEFAULT 00440 CALL stop_program(routineN,moduleN,__LINE__,& 00441 "l = 3 and m value out of bounds") 00442 CASE ( 3 ) 00443 pf = SQRT ( 35.0_dp / ( 32.0_dp * pi ) ) 00444 y ( : ) = pf * r ( 1, : ) * ( r ( 1, : )**2 - 3.0_dp * r ( 2, : )**2 ) 00445 CASE ( 2 ) 00446 pf = SQRT ( 105.0_dp / ( 16.0_dp * pi ) ) 00447 y ( : ) = pf * r ( 3, : ) * ( r ( 1, : )**2 - r ( 2, : )**2 ) 00448 CASE ( 1 ) 00449 pf = SQRT ( 21.0_dp / ( 32.0_dp * pi ) ) 00450 y ( : ) = pf * r ( 1, : ) * ( 5.0_dp * r ( 3, : ) * r ( 3, : ) - 1.0_dp ) 00451 CASE ( 0 ) 00452 pf = SQRT ( 7.0_dp / ( 16.0_dp * pi ) ) 00453 y ( : ) = pf * r ( 3, : ) * ( 5.0_dp * r ( 3, : ) * r ( 3, : ) - 3.0_dp ) 00454 CASE ( -1 ) 00455 pf = SQRT ( 21.0_dp / ( 32.0_dp * pi ) ) 00456 y ( : ) = pf * r ( 2, : ) * ( 5.0_dp * r ( 3, : ) * r ( 3, : ) - 1.0_dp ) 00457 CASE ( -2 ) 00458 pf = SQRT ( 105.0_dp / ( 16.0_dp * pi ) ) 00459 y ( : ) = pf * 2.0_dp * r ( 1, : ) * r ( 2, : ) * r ( 3, : ) 00460 CASE ( -3 ) 00461 pf = SQRT ( 35.0_dp / ( 32.0_dp * pi ) ) 00462 y ( : ) = pf * r ( 2, : ) * ( 3.0_dp * r ( 1, : )**2 - r ( 2, : )**2 ) 00463 END SELECT 00464 CASE DEFAULT 00465 IF ( m < -l .OR. m > l ) CALL stop_program(routineN,moduleN,__LINE__,& 00466 "m value out of bounds") 00467 lpm = fac ( l + ABS ( m ) ) 00468 lmm = fac ( l - ABS ( m ) ) 00469 IF ( m == 0 ) THEN 00470 t = 4.0_dp * pi 00471 ELSE 00472 t = 2.0_dp * pi 00473 END IF 00474 IF ( ABS ( lpm ) < EPSILON ( 1.0_dp ) ) THEN 00475 pf = REAL ( 2*l+1,KIND=dp) / t 00476 ELSE 00477 pf = ( REAL ( 2*l+1,KIND=dp) * lmm ) / ( t * lpm ) 00478 ENDIF 00479 pf = SQRT ( pf ) 00480 DO i = 1, SIZE ( r ,2 ) 00481 z = r ( 3, i ) 00482 ! plm = legendre ( z, l, m ) 00483 plm = legendre ( z, l, ABS ( m ) ) 00484 IF ( m == 0 ) THEN 00485 y ( i ) = pf * plm 00486 ELSE 00487 rxy = SQRT ( r ( 1, i ) ** 2 + r ( 2, i ) ** 2 ) 00488 IF ( rxy < EPSILON ( 1.0_dp ) ) THEN 00489 y ( i ) = 0.0_dp 00490 ELSE 00491 cp = r ( 1, i ) / rxy 00492 sp = r ( 2, i ) / rxy 00493 IF ( m > 0 ) THEN 00494 y ( i ) = pf * plm * cosn ( m, cp, sp ) 00495 ELSE 00496 y ( i ) = pf * plm * sinn ( -m , cp, sp ) 00497 END IF 00498 END IF 00499 END IF 00500 END DO 00501 END SELECT 00502 00503 END SUBROUTINE rvy_lm 00504 00505 ! ***************************************************************************** 00506 SUBROUTINE rry_lm ( r, y, l, m ) 00507 ! 00508 ! Real Spherical Harmonics 00509 ! _ _ 00510 ! | [(2l+1)(l-|m|)!] |1/2 m cos(m p) m>=0 00511 ! Y_lm ( t, p ) = |---------------------| P_l(cos(t)) 00512 ! |[2Pi(1+d_m0)(l+|m|)!]| sin(|m| p) m<0 00513 ! 00514 ! Input: r == (x,y,z) : normalised x^2 + y^2 + z^2 = 1 00515 ! 00516 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: r 00517 REAL(KIND=dp), INTENT(OUT) :: y 00518 INTEGER, INTENT(IN) :: l, m 00519 00520 CHARACTER(len=*), PARAMETER :: routineN = 'rry_lm', 00521 routineP = moduleN//':'//routineN 00522 00523 REAL(KIND=dp) :: cp, lmm, lpm, pf, plm, rxy, 00524 sp, t, z 00525 00526 SELECT CASE ( l ) 00527 CASE ( : -1 ) 00528 CALL stop_program(routineN,moduleN,__LINE__,"Negative l value") 00529 CASE ( 0 ) 00530 pf = SQRT ( 1.0_dp / ( 4.0_dp * pi ) ) 00531 IF ( m /= 0 ) CALL stop_program(routineN,moduleN,__LINE__,& 00532 "l = 0 and m value out of bounds") 00533 y = pf 00534 CASE ( 1 ) 00535 SELECT CASE ( m ) 00536 CASE DEFAULT 00537 CALL stop_program(routineN,moduleN,__LINE__,& 00538 "l = 1 and m value out of bounds") 00539 CASE ( 1 ) 00540 pf = SQRT ( 3.0_dp / ( 4.0_dp * pi ) ) 00541 y = pf * r ( 1 ) 00542 CASE ( 0 ) 00543 pf = SQRT ( 3.0_dp / ( 4.0_dp * pi ) ) 00544 y = pf * r ( 3 ) 00545 CASE ( -1 ) 00546 pf = SQRT ( 3.0_dp / ( 4.0_dp * pi ) ) 00547 y = pf * r ( 2 ) 00548 END SELECT 00549 CASE ( 2 ) 00550 SELECT CASE ( m ) 00551 CASE DEFAULT 00552 CALL stop_program(routineN,moduleN,__LINE__,& 00553 "l = 2 and m value out of bounds") 00554 CASE ( 2 ) 00555 pf = SQRT ( 15.0_dp / ( 16.0_dp * pi ) ) 00556 y = pf * ( r ( 1 ) * r ( 1 ) - r ( 2 ) * r ( 2 ) ) 00557 CASE ( 1 ) 00558 pf = SQRT ( 15.0_dp / ( 4.0_dp * pi ) ) 00559 y = pf * r ( 3 ) * r ( 1 ) 00560 CASE ( 0 ) 00561 pf = SQRT ( 5.0_dp / ( 16.0_dp * pi ) ) 00562 y = pf * ( 3.0_dp * r ( 3 ) * r ( 3 ) - 1.0_dp ) 00563 CASE ( -1 ) 00564 pf = SQRT ( 15.0_dp / ( 4.0_dp * pi ) ) 00565 y = pf * r ( 3 ) * r ( 2 ) 00566 CASE ( -2 ) 00567 pf = SQRT ( 15.0_dp / ( 16.0_dp * pi ) ) 00568 y = pf * 2.0_dp * r ( 1 ) * r ( 2 ) 00569 END SELECT 00570 CASE ( 3 ) 00571 SELECT CASE ( m ) 00572 CASE DEFAULT 00573 CALL stop_program(routineN,moduleN,__LINE__,"l = 3 and m value out of bounds") 00574 CASE ( 3 ) 00575 pf = SQRT ( 35.0_dp / ( 32.0_dp * pi ) ) 00576 y = pf * r ( 1 ) * ( r ( 1 )**2 - 3.0_dp * r ( 2 )**2 ) 00577 CASE ( 2 ) 00578 pf = SQRT ( 105.0_dp / ( 16.0_dp * pi ) ) 00579 y = pf * r ( 3 ) * ( r ( 1 )**2 - r ( 2 )**2 ) 00580 CASE ( 1 ) 00581 pf = SQRT ( 21.0_dp / ( 32.0_dp * pi ) ) 00582 y = pf * r ( 1 ) * ( 5.0_dp * r ( 3 ) * r ( 3 ) - 1.0_dp ) 00583 CASE ( 0 ) 00584 pf = SQRT ( 7.0_dp / ( 16.0_dp * pi ) ) 00585 y = pf * r ( 3 ) * ( 5.0_dp * r ( 3 ) * r ( 3 ) - 3.0_dp ) 00586 CASE ( -1 ) 00587 pf = SQRT ( 21.0_dp / ( 32.0_dp * pi ) ) 00588 y = pf * r ( 2 ) * ( 5.0_dp * r ( 3 ) * r ( 3 ) - 1.0_dp ) 00589 CASE ( -2 ) 00590 pf = SQRT ( 105.0_dp / ( 16.0_dp * pi ) ) 00591 y = pf * 2.0_dp * r ( 1 ) * r ( 2 ) * r ( 3 ) 00592 CASE ( -3 ) 00593 pf = SQRT ( 35.0_dp / ( 32.0_dp * pi ) ) 00594 y = pf * r ( 2 ) * ( 3.0_dp * r ( 1 )**2 - r ( 2 )**2 ) 00595 END SELECT 00596 CASE DEFAULT 00597 IF ( m < -l .OR. m > l ) CALL stop_program(routineN,moduleN,__LINE__,& 00598 "m value out of bounds") 00599 lpm = fac ( l + ABS ( m ) ) 00600 lmm = fac ( l - ABS ( m ) ) 00601 IF ( m == 0 ) THEN 00602 t = 4.0_dp * pi 00603 ELSE 00604 t = 2.0_dp * pi 00605 END IF 00606 IF ( ABS ( lpm ) < EPSILON ( 1.0_dp ) ) THEN 00607 pf = REAL ( 2*l+1,KIND=dp) / t 00608 ELSE 00609 pf = ( REAL ( 2*l+1,KIND=dp) * lmm ) / ( t * lpm ) 00610 ENDIF 00611 pf = SQRT ( pf ) 00612 z = r ( 3 ) 00613 plm = legendre ( z, l, m ) 00614 IF ( m == 0 ) THEN 00615 y = pf * plm 00616 ELSE 00617 rxy = SQRT ( r ( 1 ) ** 2 + r ( 2 ) ** 2 ) 00618 IF ( rxy < EPSILON ( 1.0_dp ) ) THEN 00619 y = 0.0_dp 00620 ELSE 00621 cp = r ( 1 ) / rxy 00622 sp = r ( 2 ) / rxy 00623 IF ( m > 0 ) THEN 00624 y = pf * plm * cosn ( m, cp, sp ) 00625 ELSE 00626 y = pf * plm * sinn ( -m , cp, sp ) 00627 END IF 00628 END IF 00629 END IF 00630 END SELECT 00631 00632 END SUBROUTINE rry_lm 00633 00634 ! ***************************************************************************** 00635 SUBROUTINE cvy_lm ( r, y, l, m ) 00636 ! 00637 ! Complex Spherical Harmonics 00638 ! _ _ 00639 ! | [(2l+1)(l-|m|)!] |1/2 m 00640 ! Y_lm ( t, p ) = |------------------| P_l(cos(t)) Exp[ i m p ] 00641 ! | [4Pi(l+|m|)!]| | 00642 ! 00643 ! Input: r == (x,y,z) : normalised x^2 + y^2 + z^2 = 1 00644 ! 00645 REAL(KIND=dp), DIMENSION(:, :), 00646 INTENT(IN) :: r 00647 COMPLEX(KIND=dp), DIMENSION(:), 00648 INTENT(OUT) :: y 00649 INTEGER, INTENT(IN) :: l, m 00650 00651 CHARACTER(len=*), PARAMETER :: routineN = 'cvy_lm', 00652 routineP = moduleN//':'//routineN 00653 00654 INTEGER :: i, n 00655 REAL(KIND=dp) :: cp, lmm, lpm, pf, plm, rxy, 00656 sp, t, ym, yp, z 00657 00658 n = SIZE ( r ,2 ) 00659 SELECT CASE ( l ) 00660 CASE ( : -1 ) 00661 CALL stop_program(routineN,moduleN,__LINE__,"Negative l value") 00662 CASE ( 0 ) 00663 pf = SQRT ( 1.0_dp / ( 4.0_dp * pi ) ) 00664 IF ( m /= 0 ) CALL stop_program(routineN,moduleN,__LINE__,& 00665 "l = 0 and m value out of bounds") 00666 y(:) = pf 00667 CASE ( 1 ) 00668 SELECT CASE ( m ) 00669 CASE DEFAULT 00670 CALL stop_program(routineN,moduleN,__LINE__,"l = 1 and m value out of bounds") 00671 CASE ( 1 ) 00672 pf = SQRT ( 3.0_dp / ( 8.0_dp * pi ) ) 00673 DO i = 1, n 00674 yp = r ( 1, i ) 00675 ym = r ( 2, i ) 00676 y ( i ) = pf * CMPLX ( yp, ym,KIND=dp) 00677 END DO 00678 CASE ( 0 ) 00679 pf = SQRT ( 3.0_dp / ( 4.0_dp * pi ) ) 00680 y ( : ) = pf * r ( 3, : ) 00681 CASE ( -1 ) 00682 pf = SQRT ( 3.0_dp / ( 8.0_dp * pi ) ) 00683 DO i = 1, n 00684 yp = r ( 1, i ) 00685 ym = r ( 2, i ) 00686 y ( i ) = pf * CMPLX ( yp, -ym,KIND=dp) 00687 END DO 00688 END SELECT 00689 CASE ( 2 ) 00690 SELECT CASE ( m ) 00691 CASE DEFAULT 00692 CALL stop_program(routineN,moduleN,__LINE__,"l = 2 and m value out of bounds") 00693 CASE ( 2 ) 00694 pf = SQRT ( 15.0_dp / ( 32.0_dp * pi ) ) 00695 DO i = 1, n 00696 yp = ( r ( 1, i ) * r ( 1, i ) - r ( 2, i ) * r ( 2, i ) ) 00697 ym = 2.0_dp * r ( 1, i ) * r ( 2, i ) 00698 y ( i ) = pf * CMPLX ( yp, ym,KIND=dp) 00699 END DO 00700 CASE ( 1 ) 00701 pf = SQRT ( 15.0_dp / ( 8.0_dp * pi ) ) 00702 DO i = 1, n 00703 yp = r ( 3, i ) * r ( 1, i ) 00704 ym = r ( 3, i ) * r ( 2, i ) 00705 y ( i ) = pf * CMPLX ( yp, ym,KIND=dp) 00706 END DO 00707 CASE ( 0 ) 00708 pf = SQRT ( 5.0_dp / ( 16.0_dp * pi ) ) 00709 y ( : ) = pf * ( 3.0_dp * r ( 3, : ) * r ( 3, : ) - 1.0_dp ) 00710 CASE ( -1 ) 00711 pf = SQRT ( 15.0_dp / ( 8.0_dp * pi ) ) 00712 DO i = 1, n 00713 yp = r ( 3, i ) * r ( 1, i ) 00714 ym = r ( 3, i ) * r ( 2, i ) 00715 y ( i ) = pf * CMPLX ( yp, -ym,KIND=dp) 00716 END DO 00717 CASE ( -2 ) 00718 pf = SQRT ( 15.0_dp / ( 32.0_dp * pi ) ) 00719 DO i = 1, n 00720 yp = ( r ( 1, i ) * r ( 1, i ) - r ( 2, i ) * r ( 2, i ) ) 00721 ym = 2.0_dp * r ( 1, i ) * r ( 2, i ) 00722 y ( i ) = pf * CMPLX ( yp, -ym,KIND=dp) 00723 END DO 00724 END SELECT 00725 CASE ( 3 ) 00726 SELECT CASE ( m ) 00727 CASE DEFAULT 00728 CALL stop_program(routineN,moduleN,__LINE__,"l = 3 and m value out of bounds") 00729 CASE ( 3 ) 00730 pf = SQRT ( 35.0_dp / ( 64.0_dp * pi ) ) 00731 DO i = 1, n 00732 yp = r ( 1, i ) * ( r ( 1, i )**2 - 3.0_dp * r ( 2, i )**2 ) 00733 ym = r ( 2, i ) * ( 3.0_dp * r ( 1, i )**2 - r ( 2, i )**2 ) 00734 y ( i ) = pf * CMPLX ( yp, ym,KIND=dp) 00735 END DO 00736 CASE ( 2 ) 00737 pf = SQRT ( 105.0_dp / ( 32.0_dp * pi ) ) 00738 DO i = 1, n 00739 yp = r ( 3, i ) * ( r ( 1, i )**2 - r ( 2, i )**2 ) 00740 ym = 2.0_dp * r ( 1, i ) * r ( 2, i ) * r ( 3, i ) 00741 y ( i ) = pf * CMPLX ( yp, ym,KIND=dp) 00742 END DO 00743 CASE ( 1 ) 00744 pf = SQRT ( 21.0_dp / ( 64.0_dp * pi ) ) 00745 DO i = 1, n 00746 yp = r ( 1, i ) * ( 5.0_dp * r ( 3, i ) * r ( 3, i ) - 1.0_dp ) 00747 ym = r ( 2, i ) * ( 5.0_dp * r ( 3, i ) * r ( 3, i ) - 1.0_dp ) 00748 y ( i ) = pf * CMPLX ( yp, ym,KIND=dp) 00749 END DO 00750 CASE ( 0 ) 00751 pf = SQRT ( 7.0_dp / ( 16.0_dp * pi ) ) 00752 y ( : ) = pf * r ( 3, : ) * ( 5.0_dp * r ( 3, : ) * r ( 3, : ) - 3.0_dp ) 00753 CASE ( -1 ) 00754 pf = SQRT ( 21.0_dp / ( 64.0_dp * pi ) ) 00755 DO i = 1, n 00756 yp = r ( 1, i ) * ( 5.0_dp * r ( 3, i ) * r ( 3, i ) - 1.0_dp ) 00757 ym = r ( 2, i ) * ( 5.0_dp * r ( 3, i ) * r ( 3, i ) - 1.0_dp ) 00758 y ( i ) = pf * CMPLX ( yp, -ym,KIND=dp) 00759 END DO 00760 CASE ( -2 ) 00761 pf = SQRT ( 105.0_dp / ( 32.0_dp * pi ) ) 00762 DO i = 1, n 00763 yp = r ( 3, i ) * ( r ( 1, i )**2 - r ( 2, i )**2 ) 00764 ym = 2.0_dp * r ( 1, i ) * r ( 2, i ) * r ( 3, i ) 00765 y ( i ) = pf * CMPLX ( yp, -ym,KIND=dp) 00766 END DO 00767 CASE ( -3 ) 00768 pf = SQRT ( 35.0_dp / ( 64.0_dp * pi ) ) 00769 DO i = 1, n 00770 yp = r ( 1, i ) * ( r ( 1, i )**2 - 3.0_dp * r ( 2, i )**2 ) 00771 ym = r ( 2, i ) * ( 3.0_dp * r ( 1, i )**2 - r ( 2, i )**2 ) 00772 y ( i ) = pf * CMPLX ( yp, -ym,KIND=dp) 00773 END DO 00774 END SELECT 00775 CASE DEFAULT 00776 IF ( m < -l .OR. m > l ) CALL stop_program(routineN,moduleN,__LINE__,& 00777 "m value out of bounds") 00778 lpm = fac ( l + ABS ( m ) ) 00779 lmm = fac ( l - ABS ( m ) ) 00780 t = 4.0_dp * pi 00781 IF ( ABS ( lpm ) < EPSILON ( 1.0_dp ) ) THEN 00782 pf = REAL ( 2*l+1,KIND=dp) / t 00783 ELSE 00784 pf = ( REAL ( 2*l+1,KIND=dp) * lmm ) / ( t * lpm ) 00785 ENDIF 00786 pf = SQRT ( pf ) 00787 DO i = 1, n 00788 z = r ( 3, i ) 00789 plm = legendre ( z, l, m ) 00790 IF ( m == 0 ) THEN 00791 y ( i ) = pf * plm 00792 ELSE 00793 rxy = SQRT ( r ( 1, i ) ** 2 + r ( 2, i ) ** 2 ) 00794 IF ( rxy < EPSILON ( 1.0_dp ) ) THEN 00795 y ( i ) = 0.0_dp 00796 ELSE 00797 cp = r ( 1, i ) / rxy 00798 sp = r ( 2, i ) / rxy 00799 IF ( m > 0 ) THEN 00800 yp = cosn ( m, cp, sp ) 00801 ym = sinn ( m, cp, sp ) 00802 y ( i ) = pf * plm * CMPLX( yp, ym,KIND=dp) 00803 ELSE 00804 yp = cosn ( -m, cp, sp ) 00805 ym = sinn ( -m, cp, sp ) 00806 y ( i ) = pf * plm * CMPLX( yp, -ym,KIND=dp) 00807 END IF 00808 END IF 00809 END IF 00810 END DO 00811 END SELECT 00812 00813 END SUBROUTINE cvy_lm 00814 00815 ! ***************************************************************************** 00816 SUBROUTINE ccy_lm ( r, y, l, m ) 00817 ! 00818 ! Complex Spherical Harmonics 00819 ! _ _ 00820 ! | [(2l+1)(l-|m|)!] |1/2 m 00821 ! Y_lm ( t, p ) = |------------------| P_l(cos(t)) Exp[ i m p ] 00822 ! | [4Pi(l+|m|)!]| | 00823 ! 00824 ! Input: r == (x,y,z) : normalised x^2 + y^2 + z^2 = 1 00825 ! 00826 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: r 00827 COMPLEX(KIND=dp), INTENT(OUT) :: y 00828 INTEGER, INTENT(IN) :: l, m 00829 00830 CHARACTER(len=*), PARAMETER :: routineN = 'ccy_lm', 00831 routineP = moduleN//':'//routineN 00832 00833 REAL(KIND=dp) :: cp, lmm, lpm, pf, plm, rxy, 00834 sp, t, ym, yp, z 00835 00836 SELECT CASE ( l ) 00837 CASE ( : -1 ) 00838 CALL stop_program(routineN,moduleN,__LINE__,"Negative l value") 00839 CASE ( 0 ) 00840 pf = SQRT ( 1.0_dp / ( 4.0_dp * pi ) ) 00841 IF ( m /= 0 ) CALL stop_program(routineN,moduleN,__LINE__,& 00842 "l = 0 and m value out of bounds") 00843 y = pf 00844 CASE ( 1 ) 00845 SELECT CASE ( m ) 00846 CASE DEFAULT 00847 CALL stop_program(routineN,moduleN,__LINE__,"l = 1 and m value out of bounds") 00848 CASE ( 1 ) 00849 pf = SQRT ( 3.0_dp / ( 8.0_dp * pi ) ) 00850 yp = r ( 1 ) 00851 ym = r ( 2 ) 00852 y = pf * CMPLX ( yp, ym,KIND=dp) 00853 CASE ( 0 ) 00854 pf = SQRT ( 3.0_dp / ( 4.0_dp * pi ) ) 00855 y = pf * r ( 3 ) 00856 CASE ( -1 ) 00857 pf = SQRT ( 3.0_dp / ( 8.0_dp * pi ) ) 00858 yp = r ( 1 ) 00859 ym = r ( 2 ) 00860 y = pf * CMPLX ( yp, -ym,KIND=dp) 00861 END SELECT 00862 CASE ( 2 ) 00863 SELECT CASE ( m ) 00864 CASE DEFAULT 00865 CALL stop_program(routineN,moduleN,__LINE__,"l = 2 and m value out of bounds") 00866 CASE ( 2 ) 00867 pf = SQRT ( 15.0_dp / ( 32.0_dp * pi ) ) 00868 yp = ( r ( 1 ) * r ( 1 ) - r ( 2 ) * r ( 2 ) ) 00869 ym = 2.0_dp * r ( 1 ) * r ( 2 ) 00870 y = pf * CMPLX ( yp, ym,KIND=dp) 00871 CASE ( 1 ) 00872 pf = SQRT ( 15.0_dp / ( 8.0_dp * pi ) ) 00873 yp = r ( 3 ) * r ( 1 ) 00874 ym = r ( 3 ) * r ( 2 ) 00875 y = pf * CMPLX ( yp, ym,KIND=dp) 00876 CASE ( 0 ) 00877 pf = SQRT ( 5.0_dp / ( 16.0_dp * pi ) ) 00878 y = pf * ( 3.0_dp * r ( 3 ) * r ( 3 ) - 1.0_dp ) 00879 CASE ( -1 ) 00880 pf = SQRT ( 15.0_dp / ( 8.0_dp * pi ) ) 00881 yp = r ( 3 ) * r ( 1 ) 00882 ym = r ( 3 ) * r ( 2 ) 00883 y = pf * CMPLX ( yp, -ym,KIND=dp) 00884 CASE ( -2 ) 00885 pf = SQRT ( 15.0_dp / ( 32.0_dp * pi ) ) 00886 yp = ( r ( 1 ) * r ( 1 ) - r ( 2 ) * r ( 2 ) ) 00887 ym = 2.0_dp * r ( 1 ) * r ( 2 ) 00888 y = pf * CMPLX ( yp, -ym,KIND=dp) 00889 END SELECT 00890 CASE ( 3 ) 00891 SELECT CASE ( m ) 00892 CASE DEFAULT 00893 CALL stop_program(routineN,moduleN,__LINE__,"l = 3 and m value out of bounds") 00894 CASE ( 3 ) 00895 pf = SQRT ( 35.0_dp / ( 64.0_dp * pi ) ) 00896 yp = r ( 1 ) * ( r ( 1 )**2 - 3.0_dp * r ( 2 )**2 ) 00897 ym = r ( 2 ) * ( 3.0_dp * r ( 1 )**2 - r ( 2 )**2 ) 00898 y = pf * CMPLX ( yp, ym,KIND=dp) 00899 CASE ( 2 ) 00900 pf = SQRT ( 105.0_dp / ( 32.0_dp * pi ) ) 00901 yp = r ( 3 ) * ( r ( 1 )**2 - r ( 2 )**2 ) 00902 ym = 2.0_dp * r ( 1 ) * r ( 2 ) * r ( 3 ) 00903 y = pf * CMPLX ( yp, ym,KIND=dp) 00904 CASE ( 1 ) 00905 pf = SQRT ( 21.0_dp / ( 64.0_dp * pi ) ) 00906 yp = r ( 1 ) * ( 5.0_dp * r ( 3 ) * r ( 3 ) - 1.0_dp ) 00907 ym = r ( 2 ) * ( 5.0_dp * r ( 3 ) * r ( 3 ) - 1.0_dp ) 00908 y = pf * CMPLX ( yp, ym,KIND=dp) 00909 CASE ( 0 ) 00910 pf = SQRT ( 7.0_dp / ( 16.0_dp * pi ) ) 00911 y = pf * r ( 3 ) * ( 5.0_dp * r ( 3 ) * r ( 3 ) - 3.0_dp ) 00912 CASE ( -1 ) 00913 pf = SQRT ( 21.0_dp / ( 64.0_dp * pi ) ) 00914 yp = r ( 1 ) * ( 5.0_dp * r ( 3 ) * r ( 3 ) - 1.0_dp ) 00915 ym = r ( 2 ) * ( 5.0_dp * r ( 3 ) * r ( 3 ) - 1.0_dp ) 00916 y = pf * CMPLX ( yp, -ym,KIND=dp) 00917 CASE ( -2 ) 00918 pf = SQRT ( 105.0_dp / ( 32.0_dp * pi ) ) 00919 yp = r ( 3 ) * ( r ( 1 )**2 - r ( 2 )**2 ) 00920 ym = 2.0_dp * r ( 1 ) * r ( 2 ) * r ( 3 ) 00921 y = pf * CMPLX ( yp, -ym,KIND=dp) 00922 CASE ( -3 ) 00923 pf = SQRT ( 35.0_dp / ( 64.0_dp * pi ) ) 00924 yp = r ( 1 ) * ( r ( 1 )**2 - 3.0_dp * r ( 2 )**2 ) 00925 ym = r ( 2 ) * ( 3.0_dp * r ( 1 )**2 - r ( 2 )**2 ) 00926 y = pf * CMPLX ( yp, -ym,KIND=dp) 00927 END SELECT 00928 CASE DEFAULT 00929 IF ( m < -l .OR. m > l ) CALL stop_program(routineN,moduleN,__LINE__,& 00930 "m value out of bounds") 00931 lpm = fac ( l + ABS ( m ) ) 00932 lmm = fac ( l - ABS ( m ) ) 00933 t = 4.0_dp * pi 00934 IF ( ABS ( lpm ) < EPSILON ( 1.0_dp ) ) THEN 00935 pf = REAL ( 2*l+1,KIND=dp) / t 00936 ELSE 00937 pf = ( REAL ( 2*l+1,KIND=dp) * lmm ) / ( t * lpm ) 00938 ENDIF 00939 pf = SQRT ( pf ) 00940 z = r ( 3 ) 00941 plm = legendre ( z, l, m ) 00942 IF ( m == 0 ) THEN 00943 y = pf * plm 00944 ELSE 00945 rxy = SQRT ( r ( 1 ) ** 2 + r ( 2 ) ** 2 ) 00946 IF ( rxy < EPSILON ( 1.0_dp ) ) THEN 00947 y = 0.0_dp 00948 ELSE 00949 cp = r ( 1 ) / rxy 00950 sp = r ( 2 ) / rxy 00951 IF ( m > 0 ) THEN 00952 yp = cosn ( m, cp, sp ) 00953 ym = sinn ( m , cp, sp ) 00954 y = pf * plm * CMPLX( yp, ym,KIND=dp) 00955 ELSE 00956 yp = cosn ( -m, cp, sp ) 00957 ym = sinn ( -m , cp, sp ) 00958 y = pf * plm * CMPLX( yp, -ym,KIND=dp) 00959 END IF 00960 END IF 00961 END IF 00962 END SELECT 00963 00964 END SUBROUTINE ccy_lm 00965 00966 ! ***************************************************************************** 00967 FUNCTION legendre ( x, l, m ) RESULT ( plm ) 00968 00969 REAL(KIND=dp), INTENT(IN) :: x 00970 INTEGER, INTENT(IN) :: l, m 00971 REAL(KIND=dp) :: plm 00972 00973 CHARACTER(len=*), PARAMETER :: routineN = 'legendre', 00974 routineP = moduleN//':'//routineN 00975 00976 INTEGER :: il, im, mm 00977 REAL(KIND=dp) :: fact, pll, pmm, pmmp1, somx2 00978 00979 IF ( ABS ( x ) > 1.0_dp ) CALL stop_program(routineN,moduleN,__LINE__,"x value > 1") 00980 SELECT CASE ( l ) 00981 CASE ( : -1 ) 00982 CALL stop_program(routineN,moduleN,__LINE__,"Negative l value") 00983 CASE ( 0 ) 00984 plm = 1.0_dp 00985 CASE ( 1 ) 00986 SELECT CASE ( ABS ( m ) ) 00987 CASE DEFAULT 00988 CALL stop_program(routineN,moduleN,__LINE__,"l = 1 and m value out of bounds") 00989 CASE ( 1 ) 00990 plm = SQRT ( 1.0_dp - x * x ) 00991 CASE ( 0 ) 00992 plm = x 00993 END SELECT 00994 CASE ( 2 ) 00995 SELECT CASE ( ABS ( m ) ) 00996 CASE DEFAULT 00997 CALL stop_program(routineN,moduleN,__LINE__,"l = 2 and m value out of bounds") 00998 CASE ( 2 ) 00999 plm = 3.0_dp * ( 1.0_dp - x * x ) 01000 CASE ( 1 ) 01001 plm = 3.0_dp * x * SQRT ( 1.0_dp - x * x ) 01002 CASE ( 0 ) 01003 plm = 1.5_dp * x * x - 0.5_dp 01004 END SELECT 01005 CASE ( 3 ) 01006 SELECT CASE ( ABS ( m ) ) 01007 CASE DEFAULT 01008 CALL stop_program(routineN,moduleN,__LINE__,"l = 3 and m value out of bounds") 01009 CASE ( 3 ) 01010 plm = 15.0_dp * ( 1.0_dp - x * x ) ** 1.5_dp 01011 CASE ( 2 ) 01012 plm = 15.0_dp * x * ( 1.0_dp - x * x ) 01013 CASE ( 1 ) 01014 plm = ( 7.5_dp * x * x - 1.5_dp ) * SQRT ( 1.0_dp - x * x ) 01015 CASE ( 0 ) 01016 plm = 2.5_dp * x**3 - 1.5_dp * x 01017 END SELECT 01018 CASE ( 4 ) 01019 SELECT CASE ( ABS ( m ) ) 01020 CASE DEFAULT 01021 CALL stop_program(routineN,moduleN,__LINE__,"l = 4 and m value out of bounds") 01022 CASE ( 4 ) 01023 plm = 105.0_dp * ( 1.0_dp - x * x ) ** 2 01024 CASE ( 3 ) 01025 plm = 105.0_dp * x * ( 1.0_dp - x * x ) ** 1.5_dp 01026 CASE ( 2 ) 01027 plm = ( 52.5_dp * x * x - 7.5_dp ) * ( 1.0_dp - x * x ) 01028 CASE ( 1 ) 01029 plm = ( 17.5_dp * x ** 3 - 7.5_dp * x ) * SQRT ( 1.0_dp - x * x ) 01030 CASE ( 0 ) 01031 plm = 4.375_dp * x ** 4 - 3.75_dp * x ** 2 + 0.375_dp 01032 END SELECT 01033 CASE ( 5 ) 01034 SELECT CASE ( ABS ( m ) ) 01035 CASE DEFAULT 01036 CALL stop_program(routineN,moduleN,__LINE__,"l = 5 and m value out of bounds") 01037 CASE ( 5 ) 01038 plm = 945.0_dp * ( 1.0_dp - x * x ) ** 2.5_dp 01039 CASE ( 4 ) 01040 plm = 945.0_dp * x * ( 1.0_dp - x * x ) ** 2 01041 CASE ( 3 ) 01042 plm = -( -472.5_dp * x * x + 52.5_dp ) * ( 1.0_dp - x * x ) ** 1.5_dp 01043 CASE ( 2 ) 01044 plm = ( 157.5_dp * x ** 3 - 52.5_dp * x ) * ( 1.0_dp - x * x ) 01045 CASE ( 1 ) 01046 plm = -( -39.375_dp * x ** 4 + 26.25_dp * x * x - & 01047 1.875_dp ) * SQRT ( 1.0_dp - x * x ) 01048 CASE ( 0 ) 01049 plm = 7.875_dp * x ** 5 - 8.75_dp * x ** 3 + 1.875_dp * x 01050 END SELECT 01051 CASE ( 6 ) 01052 SELECT CASE ( ABS ( m ) ) 01053 CASE DEFAULT 01054 CALL stop_program(routineN,moduleN,__LINE__,"l = 6 and m value out of bounds") 01055 CASE ( 6 ) 01056 plm = 10395.0_dp * ( 1.0_dp - x * x ) ** 3 01057 CASE ( 5 ) 01058 plm = 10395.0_dp * x * ( 1.0_dp - x * x ) ** 2.5_dp 01059 CASE ( 4 ) 01060 plm = ( 5197.5_dp * x * x - 472.5_dp ) * ( 1.0_dp - x * x ) ** 2 01061 CASE ( 3 ) 01062 plm = -( - 1732.5_dp * x ** 3 + 472.5_dp * x ) * & 01063 ( 1.0_dp - x * x ) ** 1.5_dp 01064 CASE ( 2 ) 01065 plm = ( 433.125_dp * x ** 4 - 236.25_dp * x ** 2 + & 01066 13.125_dp ) * ( 1.0_dp - x * x ) 01067 CASE ( 1 ) 01068 plm = -( -86.625_dp * x ** 5 + 78.75_dp * x ** 3 - & 01069 13.125_dp * x ) * SQRT ( 1.0_dp - x * x ) 01070 CASE ( 0 ) 01071 plm = 14.4375_dp * x ** 6 - 19.6875_dp * x ** 4 + & 01072 6.5625_dp * x ** 2 - 0.3125_dp 01073 END SELECT 01074 CASE DEFAULT 01075 mm = ABS ( m ) 01076 IF ( mm > l ) CALL stop_program(routineN,moduleN,__LINE__,"m out of bounds") 01077 ! use recurence from numerical recipies 01078 pmm = 1.0_dp 01079 IF ( mm > 0 ) THEN 01080 somx2 = SQRT((1.0_dp-x)*(1.0_dp+x)) 01081 fact = 1.0_dp 01082 DO im = 1, mm 01083 pmm = pmm * fact * somx2 01084 fact = fact + 2.0_dp 01085 END DO 01086 END IF 01087 IF ( l == mm ) THEN 01088 plm = pmm 01089 ELSE 01090 pmmp1 = x * REAL ( 2 * mm + 1,KIND=dp) * pmm 01091 IF ( l == mm + 1 ) THEN 01092 plm = pmmp1 01093 ELSE 01094 DO il = mm + 2, l 01095 pll = ( x * REAL ( 2 * il - 1,KIND=dp) * pmmp1 - 01096 REAL ( il + mm - 1,KIND=dp) * pmm ) / REAL ( il - mm,KIND=dp) 01097 pmm = pmmp1 01098 pmmp1 = pll 01099 END DO 01100 plm = pll 01101 END IF 01102 END IF 01103 END SELECT 01104 01105 END FUNCTION legendre 01106 01107 ! ***************************************************************************** 01108 FUNCTION dlegendre (x, l, m) RESULT (dplm) 01109 REAL(KIND=dp), INTENT(IN) :: x 01110 INTEGER, INTENT(IN) :: l, m 01111 REAL(KIND=dp) :: dplm 01112 01113 CHARACTER(len=*), PARAMETER :: routineN = 'dlegendre', 01114 routineP = moduleN//':'//routineN 01115 01116 INTEGER :: il, im, mm 01117 REAL(KIND=dp) :: dpll, dpmm, dpmmp1, dsomx2, 01118 fact, pll, plm, pmm, pmmp1, 01119 somx2 01120 01121 IF ( ABS ( x ) > 1.0_dp ) CALL stop_program(routineN,moduleN,__LINE__,"x value > 1") 01122 SELECT CASE ( l ) 01123 CASE ( 0 ) 01124 dplm = 0.0_dp 01125 CASE ( 1 ) 01126 SELECT CASE ( ABS ( m ) ) 01127 CASE DEFAULT 01128 CALL stop_program(routineN,moduleN,__LINE__,"l = 1 and m value out of bounds") 01129 CASE ( 1 ) 01130 dplm = -x / SQRT ( 1.0_dp - x * x ) 01131 CASE ( 0 ) 01132 dplm = 1.0_dp 01133 END SELECT 01134 CASE ( 2 ) 01135 SELECT CASE ( ABS ( m ) ) 01136 CASE DEFAULT 01137 CALL stop_program(routineN,moduleN,__LINE__,"l = 2 and m value out of bounds") 01138 CASE ( 2 ) 01139 dplm = -6.0_dp * x 01140 CASE ( 1 ) 01141 dplm = 3.0_dp * SQRT ( 1.0_dp - x * x ) - 3.0_dp * x * x / SQRT ( 1.0_dp - x * x ) 01142 CASE ( 0 ) 01143 dplm = 3.0_dp * x 01144 END SELECT 01145 CASE ( 3 ) 01146 SELECT CASE ( ABS ( m ) ) 01147 CASE DEFAULT 01148 CALL stop_program(routineN,moduleN,__LINE__,"l = 3 and m value out of bounds") 01149 CASE ( 3 ) 01150 dplm = -45.0_dp * SQRT( 1.0_dp - x * x ) * x 01151 CASE ( 2 ) 01152 dplm = 15.0_dp * ( 1.0_dp - x * x ) - 30.0_dp * x * x 01153 CASE ( 1 ) 01154 dplm = 15.0_dp * x * SQRT ( 1.0_dp - x * x ) - (x * ( 7.5_dp * x * x - 1.5_dp ))/ SQRT ( 1.0_dp - x * x ) 01155 CASE ( 0 ) 01156 dplm = 7.5_dp * x * x - 1.5_dp 01157 END SELECT 01158 CASE ( 4 ) 01159 SELECT CASE ( ABS ( m ) ) 01160 CASE DEFAULT 01161 CALL stop_program(routineN,moduleN,__LINE__,"l = 4 and m value out of bounds") 01162 CASE ( 4 ) 01163 dplm = - 420 * x * (1 - x*x) 01164 CASE ( 3 ) 01165 dplm = 105.0_dp * ( ( 1.0_dp - x * x ) ** 1.5_dp - 3.0_dp * x * x * ( 1.0_dp - x * x ) ** 0.5_dp) 01166 CASE ( 2 ) 01167 dplm = 105.0_dp * x * ( 1.0_dp - x * x ) -2.0_dp * x * ( 52.5_dp * x * x - 7.5_dp ) 01168 CASE ( 1 ) 01169 IF( ABS(x) - 1.0_dp < EPSILON(1.0_dp)) THEN 01170 dplm = 0.0_dp 01171 ELSE 01172 dplm = ( 17.5_dp * 3.0_dp * x ** 2 - 7.5_dp ) * SQRT ( 1.0_dp - x * x ) - & 01173 x * ( 17.5_dp * x ** 3 - 7.5_dp * x ) / SQRT ( 1.0_dp - x * x ) 01174 END IF 01175 CASE ( 0 ) 01176 dplm = 4.375_dp * 4.0_dp * x ** 3 - 3.75_dp * 2.0_dp * x 01177 END SELECT 01178 CASE ( 5 ) 01179 SELECT CASE ( ABS ( m ) ) 01180 CASE DEFAULT 01181 CALL stop_program(routineN,moduleN,__LINE__,"l = 5 and m value out of bounds") 01182 CASE ( 5 ) 01183 dplm = -945.0_dp * 5.0_dp * x * ( 1.0_dp - x * x ) ** 1.5_dp 01184 CASE ( 4 ) 01185 dplm = 945.0_dp * (( 1.0_dp - x * x ) ** 2 - 2.0_dp * x * x * ( 1.0_dp - x * x )) 01186 CASE ( 3 ) 01187 dplm = 945.0_dp * x * ( 1.0_dp - x * x ) ** 1.5_dp - & 01188 3.0_dp * x * ( 472.5_dp * x * x - 52.5_dp ) * ( 1.0_dp - x * x ) ** 0.5_dp 01189 CASE ( 2 ) 01190 dplm = (3.0_dp *157.5_dp * x ** 2 - 52.5_dp ) * ( 1.0_dp - x * x ) - & 01191 ( 157.5_dp * x ** 3 - 52.5_dp * x ) * ( - 2.0_dp * x ) 01192 CASE ( 1 ) 01193 IF( ABS(x) - 1.0_dp < EPSILON(1.0_dp)) THEN 01194 dplm = 0.0_dp 01195 ELSE 01196 dplm = - (-39.375_dp * 4.0_dp * x * x * x + 2.0_dp * 26.25_dp*x ) * SQRT ( 1.0_dp - x * x ) +& 01197 x * ( -39.375_dp * x **4 + 26.25_dp * x * x - 1.875_dp ) / SQRT ( 1.0_dp - x * x ) 01198 END IF 01199 CASE ( 0 ) 01200 dplm = 5.0_dp * 7.875_dp * x ** 4 - 3.0_dp * 8.75_dp * x ** 2 + 1.875_dp 01201 END SELECT 01202 CASE ( 6 ) 01203 SELECT CASE ( ABS ( m ) ) 01204 CASE DEFAULT 01205 CALL stop_program(routineN,moduleN,__LINE__,"l = 6 and m value out of bounds") 01206 CASE ( 6 ) 01207 dplm = -10395.0_dp * 6.0_dp * x * ( 1.0_dp - x * x ) ** 2 01208 CASE ( 5 ) 01209 dplm = 10395.0_dp * ( ( 1.0_dp - x * x ) ** 2.5_dp - 5.0_dp * x * x * ( 1.0_dp - x * x ) ** 1.5_dp ) 01210 CASE ( 4 ) 01211 dplm = 2.0_dp * 5197.5_dp * x * ( 1.0_dp - x * x ) ** 2 - & 01212 4.0_dp * x * ( 5197.5_dp * x * x - 472.5_dp ) * ( 1.0_dp - x * x ) 01213 CASE ( 3 ) 01214 dplm = - (-3.0_dp * 1732.5_dp * x * x + 472.5_dp ) * ( 1.0_dp - x * x ) ** 1.5_dp + & 01215 ( - 1732.5_dp * x ** 3 + 472.5_dp * x ) * 3.0_dp * x * ( 1.0_dp - x * x ) ** 0.5_dp 01216 CASE ( 2 ) 01217 dplm = (433.125_dp * 4.0_dp * x ** 3 -2.0_dp * 236.25_dp *x) * ( 1.0_dp - x * x ) - & 01218 2.0_dp * x * ( 433.125_dp * x ** 4 - 236.25_dp * x ** 2 + 13.125_dp ) 01219 CASE ( 1 ) 01220 IF( ABS(x) - 1.0_dp < EPSILON(1.0_dp)) THEN 01221 dplm = 0.0_dp 01222 ELSE 01223 dplm = -(-5.0_dp * 86.625_dp * x ** 4 + 3.0_dp *78.75_dp ** 2 - 13.125_dp ) * SQRT(1.0_dp-x*x) + & 01224 x * ( -86.625_dp * x ** 5 + 78.75_dp * x ** 3 - 13.125_dp * x ) / SQRT ( 1.0_dp - x * x ) 01225 END IF 01226 CASE ( 0 ) 01227 dplm = 14.4375_dp * 6.0_dp * x ** 5 - 19.6875_dp * 4.0_dp * x ** 3 + & 01228 6.5625_dp * 2.0_dp * x 01229 END SELECT 01230 CASE DEFAULT 01231 ! CALL stop_program(routineN,moduleN,__LINE__,"l > 6 dplm not implemented") 01232 mm = ABS ( m ) 01233 IF ( mm > l ) CALL stop_program(routineN,moduleN,__LINE__,"m out of bounds") 01234 ! use recurence from numerical recipies 01235 pmm = 1.0_dp 01236 dpmm = 0.0_dp 01237 IF ( mm > 0 ) THEN 01238 somx2 = SQRT((1.0_dp-x)*(1.0_dp+x)) 01239 dsomx2 = ((1e0_dp - x)*(1e0_dp + x))**(0.5e0_dp*m-1) 01240 dpmm = -m*x*dsomx2*dfac(2*m-1) 01241 fact = 1.0_dp 01242 DO im = 1, mm 01243 pmm = -pmm * fact * somx2 01244 dpmm = -dpmm 01245 fact = fact + 2.0_dp 01246 END DO 01247 END IF 01248 IF ( l == mm ) THEN 01249 plm = pmm 01250 dplm = dpmm 01251 ELSE 01252 pmmp1 = x * REAL ( 2 * mm + 1,KIND=dp) * pmm 01253 dpmmp1 = (2*mm+1)*pmm + x*(2*mm+1)*dpmm 01254 IF ( l == mm + 1 ) THEN 01255 plm = pmmp1 01256 dplm = dpmmp1 01257 ELSE 01258 DO il = mm + 2, l 01259 pll = ( x * REAL ( 2 * il - 1,KIND=dp) * pmmp1 - 01260 REAL ( il + mm - 1,KIND=dp) * pmm ) / REAL ( il - mm,KIND=dp) 01261 dpll = ((2*il-1)*pmmp1 + x*(2*il-1)*dpmmp1 - (il+mm-1)*dpmm)/(il-mm) 01262 pmm = pmmp1 01263 pmmp1 = pll 01264 dpmm = dpmmp1 01265 dpmmp1 = dpll 01266 END DO 01267 plm = pll 01268 dplm = dpll 01269 END IF 01270 END IF 01271 01272 END SELECT 01273 01274 END FUNCTION dlegendre 01275 01276 ! ***************************************************************************** 01277 FUNCTION dPof1(x, l) 01278 01279 REAL(KIND=dp), INTENT(IN) :: x 01280 INTEGER, INTENT(IN) :: l 01281 REAL(KIND=dp) :: dPof1 01282 01283 CHARACTER(len=*), PARAMETER :: routineN = 'dPof1', 01284 routineP = moduleN//':'//routineN 01285 01286 IF ( ABS ( x ) -1.0_dp > EPSILON(1.0_dp)) THEN 01287 CALL stop_program(routineN,moduleN,__LINE__,"|x| is not 1") 01288 END IF 01289 IF (x > 0.0_dp) THEN 01290 SELECT CASE ( l ) 01291 CASE ( : -1 ) 01292 CALL stop_program(routineN,moduleN,__LINE__,"Negative l value") 01293 CASE ( 0 ) 01294 dPof1 = 0.0_dp 01295 CASE ( 1 ) 01296 dPof1 = 1.0_dp 01297 CASE ( 2 ) 01298 dPof1 = 3.0_dp 01299 CASE ( 3 ) 01300 dPof1 = 6.0_dp 01301 CASE ( 4 ) 01302 dPof1 = 10.0_dp 01303 CASE ( 5 ) 01304 dPof1 = 15.0_dp 01305 CASE ( 6 ) 01306 dPof1 = 21.0_dp 01307 CASE ( 7 : ) 01308 CALL stop_program(routineN,moduleN,__LINE__,"Not implemented") 01309 END SELECT 01310 ELSE 01311 SELECT CASE ( l ) 01312 CASE ( : -1 ) 01313 CALL stop_program(routineN,moduleN,__LINE__,"Negative l value") 01314 CASE ( 0 ) 01315 dPof1 = 0.0_dp 01316 CASE ( 1 ) 01317 dPof1 = 1.0_dp 01318 CASE ( 2 ) 01319 dPof1 = -3.0_dp 01320 CASE ( 3 ) 01321 dPof1 = 6.0_dp 01322 CASE ( 4 ) 01323 dPof1 = -10.0_dp 01324 CASE ( 5 ) 01325 dPof1 = 15.0_dp 01326 CASE ( 6 ) 01327 dPof1 = -21.0_dp 01328 CASE ( 7 : ) 01329 CALL stop_program(routineN,moduleN,__LINE__,"Not implemented") 01330 END SELECT 01331 END IF 01332 01333 END FUNCTION dPof1 01334 01335 ! ***************************************************************************** 01336 FUNCTION choose ( n, k ) 01337 01338 INTEGER, INTENT(IN) :: n, k 01339 REAL(KIND=dp) :: choose 01340 01341 IF ( n >= k ) THEN 01342 choose = REAL ( NINT ( fac ( n ) / ( fac ( k ) * fac ( n - k ) ) ),KIND=dp) 01343 ELSE 01344 choose = 0.0_dp 01345 ENDIF 01346 01347 END FUNCTION choose 01348 01349 ! ***************************************************************************** 01350 FUNCTION cosn ( n, c, s ) 01351 01352 INTEGER, INTENT(IN) :: n 01353 REAL(KIND=dp), INTENT(IN) :: c, s 01354 REAL(KIND=dp) :: cosn 01355 01356 INTEGER :: i, j 01357 01358 cosn = 0.0_dp 01359 IF( ABS(c) < EPSILON ( 1.0_dp ) .OR. n==0 ) THEN 01360 IF ( MOD(n,2) == 0 ) THEN 01361 cosn = (-1.0_dp)**INT(n/2) 01362 ELSE 01363 cosn = 0.0_dp 01364 END IF 01365 ELSE 01366 DO i = n, 0, -2 01367 IF (i==0) THEN 01368 j = n 01369 IF (j==0) THEN 01370 cosn = cosn + choose ( n, j ) 01371 ELSE IF ( MOD(j/2,2) == 0 ) THEN 01372 cosn = cosn + choose ( n, j ) * s**j 01373 ELSE 01374 cosn = cosn - choose ( n, j ) * s**j 01375 END IF 01376 ELSE 01377 j = n - i 01378 IF (j==0) THEN 01379 cosn = cosn + choose ( n, j ) * c**i 01380 ELSE IF ( MOD(j/2,2) == 0 ) THEN 01381 cosn = cosn + choose ( n, j ) * c**i * s**j 01382 ELSE 01383 cosn = cosn - choose ( n, j ) * c**i * s**j 01384 END IF 01385 END IF 01386 END DO 01387 END IF 01388 01389 END FUNCTION cosn 01390 01391 ! ***************************************************************************** 01392 FUNCTION dcosn_dcp ( n, c, s ) 01393 01394 INTEGER, INTENT(IN) :: n 01395 REAL(KIND=dp), INTENT(IN) :: c, s 01396 REAL(KIND=dp) :: dcosn_dcp 01397 01398 INTEGER :: i, j 01399 01400 dcosn_dcp = 0.0_dp 01401 01402 IF( s < EPSILON ( 1.0_dp ) ) THEN 01403 dcosn_dcp = 0.0_dp 01404 ELSE 01405 DO i = n, 0, -2 01406 IF (i == 0) THEN 01407 dcosn_dcp = dcosn_dcp 01408 ELSE 01409 j = n - i 01410 IF ( MOD(j/2,2) == 0 ) THEN 01411 dcosn_dcp = dcosn_dcp + choose ( n, j ) * REAL(i,dp) * c**(i-1) * s**j 01412 ELSE 01413 dcosn_dcp = dcosn_dcp - choose ( n, j ) * REAL(i,dp) * c**(i-1) * s**j 01414 END IF 01415 END IF 01416 END DO 01417 END IF 01418 01419 END FUNCTION dcosn_dcp 01420 01421 ! ***************************************************************************** 01422 FUNCTION dcosn_dsp ( n, c, s ) 01423 01424 INTEGER, INTENT(IN) :: n 01425 REAL(KIND=dp), INTENT(IN) :: c, s 01426 REAL(KIND=dp) :: dcosn_dsp 01427 01428 INTEGER :: i, j 01429 01430 dcosn_dsp = 0.0_dp 01431 IF( c < EPSILON ( 1.0_dp ) .OR. s < EPSILON ( 1.0_dp ) ) THEN 01432 dcosn_dsp = 0.0_dp 01433 ELSE 01434 DO i = n, 0, -2 01435 j = n - i 01436 IF (j == 0) THEN 01437 dcosn_dsp = dcosn_dsp 01438 ELSE 01439 IF ( MOD(j/2,2) == 0 ) THEN 01440 dcosn_dsp = dcosn_dsp + choose ( n, j ) * REAL(j,dp) * s**(j-1) * c**i 01441 ELSE 01442 dcosn_dsp = dcosn_dsp - choose ( n, j ) * REAL(j,dp) * s**(j-1) * c**i 01443 END IF 01444 END IF 01445 END DO 01446 END IF 01447 01448 END FUNCTION dcosn_dsp 01449 01450 ! ***************************************************************************** 01451 FUNCTION sinn ( n, c, s ) 01452 01453 INTEGER, INTENT(IN) :: n 01454 REAL(KIND=dp), INTENT(IN) :: c, s 01455 REAL(KIND=dp) :: sinn 01456 01457 INTEGER :: i, j 01458 01459 sinn = 0.0_dp 01460 01461 IF( ABS(s) < EPSILON ( 1.0_dp ) .OR. n==0) THEN 01462 sinn = 0.0_dp 01463 ELSE IF( ABS(c) < EPSILON ( 1.0_dp ) ) THEN 01464 IF ( MOD(n,2) == 0 ) THEN 01465 sinn = 0.0_dp 01466 ELSE 01467 sinn = s*(-1.0_dp)**INT((n-1)/2) 01468 END IF 01469 ELSE 01470 DO i = n-1, 0, -2 01471 IF (i==0) THEN 01472 j = n 01473 IF ( MOD(j/2,2) == 0 ) THEN 01474 sinn = sinn + choose ( n, j ) * s**j 01475 ELSE 01476 sinn = sinn - choose ( n, j ) * s**j 01477 END IF 01478 ELSE 01479 j = n - i 01480 IF ( MOD(j/2,2) == 0 ) THEN 01481 sinn = sinn + choose ( n, j ) * c**i * s**j 01482 ELSE 01483 sinn = sinn - choose ( n, j ) * c**i * s**j 01484 END IF 01485 END IF 01486 END DO 01487 END IF 01488 01489 END FUNCTION sinn 01490 01491 ! ***************************************************************************** 01492 FUNCTION dsinn_dcp ( n, c, s ) 01493 01494 INTEGER, INTENT(IN) :: n 01495 REAL(KIND=dp), INTENT(IN) :: c, s 01496 REAL(KIND=dp) :: dsinn_dcp 01497 01498 INTEGER :: i, j 01499 01500 dsinn_dcp = 0.0_dp 01501 01502 IF( c < EPSILON ( 1.0_dp ) .OR. s < EPSILON ( 1.0_dp ) ) THEN 01503 dsinn_dcp = 0.0_dp 01504 ELSE 01505 DO i = n-1, 0, -2 01506 IF(i == 0) THEN 01507 dsinn_dcp = dsinn_dcp 01508 ELSE 01509 j = n - i 01510 IF ( MOD(j/2,2) == 0 ) THEN 01511 dsinn_dcp = dsinn_dcp + choose ( n, j ) * REAL(i,dp) * c**(i-1) * s**j 01512 ELSE 01513 dsinn_dcp = dsinn_dcp - choose ( n, j ) * REAL(i,dp) * c**(i-1) * s**j 01514 END IF 01515 END IF 01516 END DO 01517 END IF 01518 01519 END FUNCTION dsinn_dcp 01520 01521 ! ***************************************************************************** 01522 FUNCTION dsinn_dsp ( n, c, s ) 01523 01524 INTEGER, INTENT(IN) :: n 01525 REAL(KIND=dp), INTENT(IN) :: c, s 01526 REAL(KIND=dp) :: dsinn_dsp 01527 01528 INTEGER :: i, j 01529 01530 dsinn_dsp = 0.0_dp 01531 01532 IF( c < EPSILON ( 1.0_dp ) .OR. s < EPSILON ( 1.0_dp ) ) THEN 01533 dsinn_dsp = 0.0_dp 01534 ELSE 01535 DO i = n-1, 0, -2 01536 j = n - i 01537 IF (j == 0) THEN 01538 dsinn_dsp = dsinn_dsp 01539 ELSE 01540 IF ( MOD(j/2,2) == 0 ) THEN 01541 dsinn_dsp = dsinn_dsp + choose ( n, j ) * c**i * REAL(j,dp) * s**(j-1) 01542 ELSE 01543 dsinn_dsp = dsinn_dsp - choose ( n, j ) * c**i * REAL(j,dp) * s**(j-1) 01544 END IF 01545 END IF 01546 END DO 01547 END IF 01548 01549 END FUNCTION dsinn_dsp 01550 01551 ! ***************************************************************************** 01552 SUBROUTINE drvy_lm ( r, dy, l, m ) 01553 ! 01554 ! Real Spherical Harmonics 01555 ! _ _ 01556 ! | [(2l+1)(l-|m|)!] |1/2 m cos(m p) m>=0 01557 ! Y_lm ( t, p ) = |---------------------| P_l(cos(t)) 01558 ! |[2Pi(1+d_m0)(l+|m|)!]| sin(|m| p) m<0 01559 ! 01560 ! Input: r == (x,y,z) : normalised x^2 + y^2 + z^2 = 1 01561 ! 01562 REAL(KIND=dp), DIMENSION(:, :), 01563 INTENT(IN) :: r 01564 REAL(KIND=dp), DIMENSION(:, :), 01565 INTENT(OUT) :: dy 01566 INTEGER, INTENT(IN) :: l, m 01567 01568 CHARACTER(len=*), PARAMETER :: routineN = 'drvy_lm', 01569 routineP = moduleN//':'//routineN 01570 01571 INTEGER :: i 01572 REAL(KIND=dp) :: cp, dcp_dx, dcp_dy, dplm, 01573 dsp_dx, dsp_dy, lmm, lpm, pf, 01574 plm, rxy, sp, t, z 01575 01576 SELECT CASE ( l ) 01577 CASE ( : -1 ) 01578 CALL stop_program(routineN,moduleN,__LINE__,"Negative l value") 01579 CASE ( 0 ) 01580 pf = SQRT ( 1.0_dp / ( 4.0_dp * pi ) ) 01581 IF ( m /= 0 ) CALL stop_program(routineN,moduleN,__LINE__,& 01582 "l = 0 and m value out of bounds") 01583 dy(1:3,:) = 0.0_dp 01584 CASE ( 1 ) 01585 SELECT CASE ( m ) 01586 CASE DEFAULT 01587 CALL stop_program(routineN,moduleN,__LINE__,"l = 1 and m value out of bounds") 01588 CASE ( 1 ) 01589 pf = SQRT ( 3.0_dp / ( 4.0_dp * pi ) ) 01590 dy (1, : ) = pf 01591 dy (2, : ) = 0.0_dp 01592 dy (3, : ) = 0.0_dp 01593 CASE ( 0 ) 01594 pf = SQRT ( 3.0_dp / ( 4.0_dp * pi ) ) 01595 dy (1, : ) = 0.0_dp 01596 dy (2, : ) = 0.0_dp 01597 dy (3, : ) = pf 01598 CASE ( -1 ) 01599 pf = SQRT ( 3.0_dp / ( 4.0_dp * pi ) ) 01600 dy (1, : ) = 0.0_dp 01601 dy (2, : ) = pf 01602 dy (3, : ) = 0.0_dp 01603 END SELECT 01604 CASE ( 2 ) 01605 SELECT CASE ( m ) 01606 CASE DEFAULT 01607 CALL stop_program(routineN,moduleN,__LINE__,"l = 2 and m value out of bounds") 01608 CASE ( 2 ) 01609 pf = SQRT ( 15.0_dp / ( 16.0_dp * pi ) ) 01610 dy (1, : ) = pf * 2.0_dp * r ( 1, : ) 01611 dy (2, : ) = -pf * 2.0_dp * r ( 2, : ) 01612 dy (3, : ) = 0.0_dp 01613 CASE ( 1 ) 01614 pf = SQRT ( 15.0_dp / ( 4.0_dp * pi ) ) 01615 dy (1, : ) = pf * r ( 3, : ) 01616 dy (2, : ) = 0.0_dp 01617 dy (3, : ) = pf * r ( 1, : ) 01618 CASE ( 0 ) 01619 pf = SQRT ( 5.0_dp / ( 16.0_dp * pi ) ) 01620 dy (1, : ) = 0.0_dp 01621 dy (2, : ) = 0.0_dp 01622 dy (3, : ) = pf * 6.0_dp * r ( 3, : ) 01623 CASE ( -1 ) 01624 pf = SQRT ( 15.0_dp / ( 4.0_dp * pi ) ) 01625 dy (1, : ) = 0.0_dp 01626 dy (2, : ) = pf * r ( 3, : ) 01627 dy (3, : ) = pf * r ( 2, : ) 01628 CASE ( -2 ) 01629 pf = SQRT ( 15.0_dp / ( 16.0_dp * pi ) ) 01630 dy (1, : ) = pf * 2.0_dp * r ( 2, : ) 01631 dy (2, : ) = pf * 2.0_dp * r ( 1, : ) 01632 dy (3, : ) = 0.0_dp 01633 END SELECT 01634 CASE ( 3 ) 01635 SELECT CASE ( m ) 01636 CASE DEFAULT 01637 CALL stop_program(routineN,moduleN,__LINE__,"l = 3 and m value out of bounds") 01638 CASE ( 3 ) 01639 pf = SQRT ( 35.0_dp / ( 32.0_dp * pi ) ) 01640 dy (1, : ) = pf * 3.0_dp * ( r ( 1, : )*r ( 1, : ) - r ( 2, : )* r ( 2, : )) 01641 dy (2, : ) = -pf * 6.0_dp * r ( 1, : ) * r ( 2, : ) 01642 dy (3, : ) = 0.0_dp 01643 CASE ( 2 ) 01644 pf = SQRT ( 105.0_dp / ( 16.0_dp * pi ) ) 01645 dy (1, : ) = pf * 2.0_dp * r ( 3, : ) * r ( 1, : ) 01646 dy (2, : ) = -pf * 2.0_dp * r ( 3, : ) * r ( 2, : ) 01647 dy (3, : ) = pf * ( r ( 1, : )**2 - r ( 2, : )**2 ) 01648 CASE ( 1 ) 01649 pf = SQRT ( 21.0_dp / ( 32.0_dp * pi ) ) 01650 dy (1, : ) = pf * ( 5.0_dp * r ( 3, : ) * r ( 3, : ) - 1.0_dp ) 01651 dy (2, : ) = 0.0_dp 01652 dy (3, : ) = pf * 10.0_dp * r ( 1, : ) * r ( 3, : ) 01653 CASE ( 0 ) 01654 pf = SQRT ( 7.0_dp / ( 16.0_dp * pi ) ) 01655 dy (1, : ) = 0.0_dp 01656 dy (2, : ) = 0.0_dp 01657 dy (3, : ) = pf * ( 15.0_dp * r ( 3, : ) * r ( 3, : ) - 3.0_dp) 01658 CASE ( -1 ) 01659 pf = SQRT ( 21.0_dp / ( 32.0_dp * pi ) ) 01660 dy (1, : ) = 0.0_dp 01661 dy (2, : ) = pf * ( 5.0_dp * r ( 3, : ) * r ( 3, : ) - 1.0_dp ) 01662 dy (3, : ) = pf * 10.0_dp * r ( 2, : ) * r ( 3, : ) 01663 CASE ( -2 ) 01664 pf = SQRT ( 105.0_dp / ( 16.0_dp * pi ) ) 01665 dy (1, : ) = pf * 2.0_dp * r ( 2, : ) * r ( 3, : ) 01666 dy (2, : ) = pf * 2.0_dp * r ( 1, : ) * r ( 3, : ) 01667 dy (3, : ) = pf * 2.0_dp * r ( 1, : ) * r ( 2, : ) 01668 CASE ( -3 ) 01669 pf = SQRT ( 35.0_dp / ( 32.0_dp * pi ) ) 01670 dy (1, : ) = pf * 6.0_dp * r ( 1, : ) * r ( 2, : ) 01671 dy (2, : ) = pf * 3.0_dp * ( r ( 1, : ) * r ( 1, : ) - r ( 2, : ) * r ( 2, : )) 01672 dy (3, : ) = 0.0_dp 01673 END SELECT 01674 CASE DEFAULT 01675 IF ( m < -l .OR. m > l ) CALL stop_program(routineN,moduleN,__LINE__,& 01676 "m value out of bounds") 01677 lpm = fac ( l + ABS ( m ) ) 01678 lmm = fac ( l - ABS ( m ) ) 01679 IF ( m == 0 ) THEN 01680 t = 4.0_dp * pi 01681 ELSE 01682 t = 2.0_dp * pi 01683 END IF 01684 IF ( ABS ( lpm ) < EPSILON ( 1.0_dp ) ) THEN 01685 pf = REAL ( 2*l+1,KIND=dp) / t 01686 ELSE 01687 pf = ( REAL ( 2*l+1,KIND=dp) * lmm ) / ( t * lpm ) 01688 ENDIF 01689 pf = SQRT ( pf ) 01690 DO i = 1, SIZE ( r ,2 ) 01691 z = r ( 3, i ) 01692 dplm = dlegendre ( z, l, ABS ( m ) ) 01693 plm = legendre ( z, l, m ) 01694 01695 IF ( m == 0 ) THEN 01696 dy (1, i ) = 0.0_dp 01697 dy (2, i ) = 0.0_dp 01698 dy (3, i ) = pf * dplm 01699 ELSE 01700 rxy = SQRT ( r ( 1, i ) ** 2 + r ( 2, i ) ** 2 ) 01701 IF ( rxy < EPSILON ( 1.0_dp ) ) THEN 01702 dy (1:3, i ) = 0.0_dp 01703 ELSE 01704 cp = r ( 1, i ) / rxy 01705 sp = r ( 2, i ) / rxy 01706 dcp_dx = (rxy-2.0_dp*r ( 1, i ) ** 2) / (rxy*rxy) 01707 dcp_dy = -2.0_dp*r ( 1, i ) * r ( 2, i ) / (rxy*rxy) 01708 dsp_dy = (rxy-2.0_dp*r ( 2, i ) ** 2) / (rxy*rxy) 01709 dsp_dx = -2.0_dp*r ( 1, i ) * r ( 2, i ) / (rxy*rxy) 01710 IF ( m > 0 ) THEN 01711 dy (1, i ) = pf * plm * (dcosn_dcp ( m, cp, sp ) * dcp_dx +& 01712 dcosn_dsp( m, cp, sp ) * dsp_dx) 01713 dy (2, i ) = pf * plm * (dcosn_dcp ( m, cp, sp ) * dcp_dy +& 01714 dcosn_dsp( m, cp, sp ) * dsp_dy) 01715 dy (3, i ) = pf * dplm * cosn ( m, cp, sp ) 01716 ELSE 01717 dy (1, i ) = pf * plm * (dsinn_dcp ( -m, cp, sp ) * dcp_dx +& 01718 dsinn_dsp( -m, cp, sp ) * dsp_dx) 01719 dy (2, i ) = pf * plm * (dsinn_dcp ( -m, cp, sp ) * dcp_dy +& 01720 dsinn_dsp( -m, cp, sp ) * dsp_dy) 01721 dy (3, i ) = pf * dplm * sinn ( -m , cp, sp ) 01722 END IF 01723 END IF 01724 END IF 01725 END DO 01726 END SELECT 01727 01728 END SUBROUTINE drvy_lm 01729 01730 END MODULE spherical_harmonics
1.7.3