diff options
Diffstat (limited to 'Dragon/src/MOCCHR.f')
| -rw-r--r-- | Dragon/src/MOCCHR.f | 158 |
1 files changed, 158 insertions, 0 deletions
diff --git a/Dragon/src/MOCCHR.f b/Dragon/src/MOCCHR.f new file mode 100644 index 0000000..dad21fb --- /dev/null +++ b/Dragon/src/MOCCHR.f @@ -0,0 +1,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 |