| 1 |
module shapes |
| 2 |
implicit none |
| 3 |
PRIVATE |
| 4 |
|
| 5 |
INTEGER, PARAMETER:: CHEBYSHEV_TN = 1 |
| 6 |
INTEGER, PARAMETER:: CHEBYSHEV_UN = 2 |
| 7 |
INTEGER, PARAMETER:: LAGUERRE = 3 |
| 8 |
INTEGER, PARAMETER:: HERMITE = 4 |
| 9 |
|
| 10 |
public :: do_shape_pair |
| 11 |
|
| 12 |
|
| 13 |
SUBROUTINE Get_Associated_Legendre(x, l, m, lmax, plm, dlm) |
| 14 |
|
| 15 |
! Purpose: Compute the associated Legendre functions |
| 16 |
! Plm(x) and their derivatives Plm'(x) |
| 17 |
! Input : x --- Argument of Plm(x) |
| 18 |
! l --- Order of Plm(x), l = 0,1,2,...,n |
| 19 |
! m --- Degree of Plm(x), m = 0,1,2,...,N |
| 20 |
! lmax --- Physical dimension of PLM and DLM |
| 21 |
! Output: PLM(l,m) --- Plm(x) |
| 22 |
! DLM(l,m) --- Plm'(x) |
| 23 |
|
| 24 |
real (kind=8), intent(in) :: x |
| 25 |
integer, intent(in) :: lmax, l, m |
| 26 |
real (kind=8), dimension(0:MM,0:N), intent(inout) :: PLM(0:lmax, 0:m) |
| 27 |
real (kind=8), dimension(0:MM,0:N), intent(inout) :: DLM(0:lmax, 0:m) |
| 28 |
integer :: i, j |
| 29 |
real (kind=8) :: xq, xs |
| 30 |
integer :: ls |
| 31 |
|
| 32 |
! zero out both arrays: |
| 33 |
DO I = 0, m |
| 34 |
DO J = 0, l |
| 35 |
PLM(J,I) = 0.0D0 |
| 36 |
DLM(J,I) = 0.0D0 |
| 37 |
end DO |
| 38 |
end DO |
| 39 |
|
| 40 |
! start with 0,0: |
| 41 |
PLM(0,0) = 1.0D0 |
| 42 |
|
| 43 |
! x = +/- 1 functions are easy: |
| 44 |
IF (abs(X).EQ.1.0D0) THEN |
| 45 |
DO I = 1, m |
| 46 |
PLM(0, I) = X**I |
| 47 |
DLM(0, I) = 0.5D0*I*(I+1.0D0)*X**(I+1) |
| 48 |
end DO |
| 49 |
DO J = 1, m |
| 50 |
DO I = 1, l |
| 51 |
IF (I.EQ.1) THEN |
| 52 |
DLM(I, J) = 1.0D+300 |
| 53 |
ELSE IF (I.EQ.2) THEN |
| 54 |
DLM(I, J) = -0.25D0*(J+2)*(J+1)*J*(J-1)*X**(J+1) |
| 55 |
ENDIF |
| 56 |
end DO |
| 57 |
end DO |
| 58 |
RETURN |
| 59 |
ENDIF |
| 60 |
|
| 61 |
LS = 1 |
| 62 |
IF (abs(X).GT.1.0D0) LS = -1 |
| 63 |
XQ = sqrt(LS*(1.0D0-X*X)) |
| 64 |
XS = LS*(1.0D0-X*X) |
| 65 |
|
| 66 |
DO I = 1, l |
| 67 |
PLM(I, I) = -LS*(2.0D0*I-1.0D0)*XQ*PLM(I-1, I-1) |
| 68 |
enddo |
| 69 |
|
| 70 |
DO I = 0, l |
| 71 |
PLM(I, I+1)=(2.0D0*I+1.0D0)*X*PLM(I, I) |
| 72 |
enddo |
| 73 |
|
| 74 |
DO I = 0, l |
| 75 |
DO J = I+2, m |
| 76 |
PLM(I, J)=((2.0D0*J-1.0D0)*X*PLM(I,J-1) - (I+J-1.0D0)*PLM(I,J-2))/(J-I) |
| 77 |
end DO |
| 78 |
end DO |
| 79 |
|
| 80 |
DLM(0, 0)=0.0D0 |
| 81 |
|
| 82 |
DO J = 1, m |
| 83 |
DLM(0, J)=LS*J*(PLM(0,J-1)-X*PLM(0,J))/XS |
| 84 |
end DO |
| 85 |
|
| 86 |
DO I = 1, l |
| 87 |
DO J = I, m |
| 88 |
DLM(I,J) = LS*I*X*PLM(I, J)/XS + (J+I)*(J-I+1.0D0)/XQ*PLM(I-1, J) |
| 89 |
end DO |
| 90 |
end DO |
| 91 |
|
| 92 |
RETURN |
| 93 |
END SUBROUTINE Get_Associated_Legendre |
| 94 |
|
| 95 |
|
| 96 |
subroutine Get_Orthogonal_Polynomial(x, m, function_type, pl, dpl) |
| 97 |
|
| 98 |
! Purpose: Compute orthogonal polynomials: Tn(x) or Un(x), |
| 99 |
! or Ln(x) or Hn(x), and their derivatives |
| 100 |
! Input : function_type --- Function code |
| 101 |
! =1 for Chebyshev polynomial Tn(x) |
| 102 |
! =2 for Chebyshev polynomial Un(x) |
| 103 |
! =3 for Laguerre polynomial Ln(x) |
| 104 |
! =4 for Hermite polynomial Hn(x) |
| 105 |
! n --- Order of orthogonal polynomials |
| 106 |
! x --- Argument of orthogonal polynomials |
| 107 |
! Output: PL(n) --- Tn(x) or Un(x) or Ln(x) or Hn(x) |
| 108 |
! DPL(n)--- Tn'(x) or Un'(x) or Ln'(x) or Hn'(x) |
| 109 |
|
| 110 |
real(kind=8), intent(in) :: x |
| 111 |
integer, intent(in):: m |
| 112 |
integer, intent(in):: function_type |
| 113 |
real(kind=8), dimension(0:n), intent(inout) :: pl, dpl |
| 114 |
|
| 115 |
real(kind=8) :: a, b, c, y0, y1, dy0, dy1 |
| 116 |
|
| 117 |
A = 2.0D0 |
| 118 |
B = 0.0D0 |
| 119 |
C = 1.0D0 |
| 120 |
Y0 = 1.0D0 |
| 121 |
Y1 = 2.0D0*X |
| 122 |
DY0 = 0.0D0 |
| 123 |
DY1 = 2.0D0 |
| 124 |
PL(0) = 1.0D0 |
| 125 |
PL(1) = 2.0D0*X |
| 126 |
DPL(0) = 0.0D0 |
| 127 |
DPL(1) = 2.0D0 |
| 128 |
IF (function_type.EQ.CHEBYSHEV_TN) THEN |
| 129 |
Y1 = X |
| 130 |
DY1 = 1.0D0 |
| 131 |
PL(1) = X |
| 132 |
DPL(1) = 1.0D0 |
| 133 |
ELSE IF (function_type.EQ.LAGUERRE) THEN |
| 134 |
Y1 = 1.0D0-X |
| 135 |
DY1 = -1.0D0 |
| 136 |
PL(1) = 1.0D0-X |
| 137 |
DPL(1) = -1.0D0 |
| 138 |
ENDIF |
| 139 |
DO K = 2, N |
| 140 |
IF (function_type.EQ.LAGUERRE) THEN |
| 141 |
A = -1.0D0/K |
| 142 |
B = 2.0D0+A |
| 143 |
C = 1.0D0+A |
| 144 |
ELSE IF (function_type.EQ.HERMITE) THEN |
| 145 |
C = 2.0D0*(K-1.0D0) |
| 146 |
ENDIF |
| 147 |
YN = (A*X+B)*Y1-C*Y0 |
| 148 |
DYN = A*Y1+(A*X+B)*DY1-C*DY0 |
| 149 |
PL(K) = YN |
| 150 |
DPL(K) = DYN |
| 151 |
Y0 = Y1 |
| 152 |
Y1 = YN |
| 153 |
DY0 = DY1 |
| 154 |
DY1 = DYN |
| 155 |
end DO |
| 156 |
RETURN |
| 157 |
|
| 158 |
end subroutine Get_Orthogonal_Polynomial |
