1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
|
*DECK MOCCHR
SUBROUTINE MOCCHR(NDIM,NANI,NFUNL,NMU,XMUANG,PHI1,PHI2,R)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Generates all spherical harmonics R(L,M) at point (XMUANG,PHI1).
*
*Copyright:
* Copyright (C) 2002 Ecole Polytechnique de Montreal
* This library is free software; you can redistribute it and/or
* modify it under the terms of the GNU Lesser General Public
* License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version
*
*Author(s): R. Roy
*
*Parameters: input
* NDIM number of dimensions for the geometry.
* NANI scattering anisotropy (=0 for isotropic scattering).
* NFUNL number of spherical harmonics per polar angle.
* NMU order of the polar quadrature in 2D and 1 in 3D.
* XMUANG cosines of the different tracking polar angles
* (polar quadrature in 2D and tracking angles in 3D).
* PHI1 first cosine of the tracking azimuthal angle.
* PHI2 second cosine of the tracking azimuthal angle.
*
*Parameters: output
* R spherical harmonics.
*
*Comments:
* for 0 <= L <= NANI (and for -L <= M <= L)
* and for all angles (for 1 <= IMU <= NMU).
* Definition of spherical harmonics R(L,M):
* R(L, M)= FACT(L,M)*P(L,M)*COS(M*PHI1) for 0 < M <= L
* R(L, 0)= P(L,0)
* R(L,-M)= FACT(L,M)*P(L,M)*SIN(M*PHI1) for 0 < M <= L
* where FACT(L,M)= SQRT( 2*(L-M)!/(L+M)! )
*
*-----------------------------------------------------------------------
*
IMPLICIT NONE
*----
* SUBROUTINE ARGUMENTS
*----
INTEGER NDIM,NANI,NFUNL,NMU
REAL XMUANG(NMU),R(NMU,NFUNL)
DOUBLE PRECISION PHI1,PHI2
*----
* LOCAL VARIABLES
*----
INTEGER IMU,L,M,LPM,LMMP1,IND,NSELEC
LOGICAL LROK
DOUBLE PRECISION DPHI,DCOP2,DSIP2,DL00,DL10,DL01,DL11,DM0,DM1,
> DM2,DCOS,DSIN,DCOT,DFAC,DZERO,DONE
DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: RWORK
PARAMETER ( DZERO= 0.0D0, DONE= 1.0D0 )
*
* INDEX FOR MATRIX 'RWORK'
IND(L,M)= L*(L+1) + M + 1
*
***** FIRST, COMPUTES ALL SPHERICAL HARMONICS
*
* GENERATES ALL LEGENDRE POLYNOMIALS P(L,M) AT POINT XMU
* FOR 0 <= L <= NANI (AND FOR 0 <= M <= L)
* USING THE RECURENCE RELATIONS:
* P(L+1,M)= ((2*L+1)*XMU*P(L,M)-(L+M)*P(L-1,M))/(L-M+1)
* P(L,M+1)= 2*M*XMU*P(L,M)/SQRT(1-XMU**2)-(L+M)*(L-M+1)*P(L,M-1)
*
* ESTABLISH SIMPLE CASES
ALLOCATE(RWORK(NMU,(NANI+1)*(NANI+1)))
IF( NANI.GE.0 )THEN
DO 10 IMU= 1, NMU
RWORK(IMU,IND(0,0))= DONE
10 CONTINUE
ELSE
CALL XABORT('MOCCHR: THE FIRST ARGUMENT MUST NON NEGATIVE')
ENDIF
*
IF( NANI.GE.1 )THEN
DPHI = SIGN(ACOS(PHI1),PHI2)
DO 40 IMU= 1, NMU
DCOS= DBLE(XMUANG(IMU))
DSIN= SQRT( DONE - DCOS *DCOS )
RWORK(IMU,IND(1,-1))= DSIN*PHI2
RWORK(IMU,IND(1, 0))= DCOS
RWORK(IMU,IND(1, 1))= DSIN*PHI1
*
IF( NANI.GE.2 )THEN
* RECURENCE PLG(IND(L,0)),PLG(IND(L,1) FOR 2 <= L <= NANI
DCOT= DCOS/DSIN
DL00= DONE
DL10= DCOS
DL01= DZERO
DL11= DSIN
DO 30 L= 1, NANI-1
* IF M=1, L=L+1 THEN L+M=L+2, L-M+1=L+1
LPM= L + 2
LMMP1= L + 1
DFAC= (DONE+DONE)/DBLE(LPM*LMMP1)
DM0=(DBLE(2*L+1)*DCOS*DL10 - DBLE(L)*DL00)/DBLE(L+1)
DM1=(DBLE(2*L+1)*DCOS*DL11 - DBLE(L+1)*DL01)/DBLE(L)
* ESTABLISH RELATIONS FOR L=L+1, ABS(M)<=1
RWORK(IMU,IND(L+1,-1))= SQRT(DFAC)*DM1*PHI2
RWORK(IMU,IND(L+1, 0))= DM0
RWORK(IMU,IND(L+1, 1))= SQRT(DFAC)*DM1*PHI1
DL00= DL10
DL01= DL11
DL10= DM0
DL11= DM1
* RECURENCE PLG(IND(L,M)) FOR 2 <= M <= L
DO 20 M= 1, L
* HERE DM0=PLG(L+1,0), DM1=PLG(L+1,1)
* ESTABLISH RELATIONS FOR L=L+1, ABS(M)>1
DFAC= DFAC/DBLE((LMMP1-1)*(LPM+1))
DCOP2= COS(DBLE(M+1)*DPHI)
DSIP2= SIN(DBLE(M+1)*DPHI)
DM2= DBLE(2*M)*DCOT*DM1 - DBLE(LPM*LMMP1)*DM0
RWORK(IMU,IND(L+1,-(M+1)))= SQRT(DFAC)*DM2*DSIP2
RWORK(IMU,IND(L+1, M+1 ))= SQRT(DFAC)*DM2*DCOP2
DM0= DM1
DM1= DM2
LPM= LPM + 1
LMMP1= LMMP1 - 1
20 CONTINUE
30 CONTINUE
ENDIF
40 CONTINUE
ENDIF
*
***** SELECTS THE GOOD SPHERICAL HARMONICS RWORK(L,M) FUNCTIONS
* FOR NDIM=1(SLAB),2(TWO-D RECT),3(THREE-D).
* COMPRESSES RWORK INTO R.
*
NSELEC= 0
LROK= .FALSE.
DO 80 L= 0, NANI
DO 70 M= -L, L
IF( NDIM.EQ.1 )THEN
LROK= M.EQ.0
ELSEIF( NDIM.EQ.2 )THEN
LROK= MOD(L+M,2).EQ.0
ELSEIF( NDIM.EQ.3 )THEN
LROK= .TRUE.
ENDIF
IF( LROK )THEN
NSELEC= NSELEC+1
DO 50 IMU= 1, NMU
R(IMU,NSELEC)= REAL(RWORK(IMU,IND(L,M)))
50 CONTINUE
ENDIF
70 CONTINUE
80 CONTINUE
IF(NSELEC.NE.NFUNL) CALL XABORT('MOCCHR: INVALID NSELEC')
DEALLOCATE(RWORK)
*
RETURN
END
|
