C
```      SUBROUTINE SOURCE(lin,lou,MPRINT,MWAVE,MTYPE,p,h,AMSOUR,phsour)
C
SAVE T,SF,AUX
DIMENSION T(3,3),SF(3),IPAR(4),PAR(6),p(3),h(3)
COMMON/SOUR/ROS
C
C     RADIATION PATTERNS OF THE SOURCE
C
C     FOR MWAVE.NE.0, COMPUTE THE RADIATION PATTERNS
C
c
PI=4.*ATAN(1.)
IF(MWAVE.NE.0) GO TO 1
if(mtype.eq.0)return
WRITE(lou,100) IPAR,PAR
go to 10
1 continue
if(mtype.eq.0)then
amsour=1.
phsour=0.
return
end if
sqsl=0.
do 9 i=1,3
sqsl=sqsl+p(i)*p(i)
9 continue
10 continue
C
C
GO TO (11,21,31,41),MTYPE
C
C
C.....SINGLE FORCE
C
11 CONTINUE
IF(MWAVE.NE.0)GO TO 12
AUX=PAR(2)/(4.*PI*ROS)
sfaz=PAR(1)
sfdc=PAR(3)
sf(1)=cos(sfaz)*cos(sfdc)
sf(2)=sin(sfaz)*cos(sfdc)
sf(3)=sin(sfdc)
RETURN
C
12 continue
amsour=0.
do 13 i=1,3
amsour=amsour+h(i)*sf(i)
13 continue
AMSOUR=AUX*sqsl*amsour
GO TO 51
C
C.....DOUBLE COUPLE
C
21 CONTINUE
IF(MWAVE.NE.0)GO TO 22
SM=PAR(2)
FS=PAR(3)
D=PAR(1)
AL=PAR(4)
SNF=SIN(FS)
CSF=COS(FS)
SND=SIN(D)
CSD=COS(D)
SNL=SIN(AL)
CSL=COS(AL)
A=CSL*CSF+CSD*SNL*SNF
B=CSL*SNF-CSD*SNL*CSF
C=SNF*SND
D=SNL*SND
T(1,1)=-2.*A*C
T(1,2)=SND*(A*CSF-B*SNF)
T(2,1)=T(1,2)
T(1,3)=-(A*CSD-C*D)
T(3,1)=T(1,3)
T(2,2)=2.*SND*csf*b
T(2,3)=-(B*CSD+D*SND*CSF)
T(3,2)=T(2,3)
T(3,3)=2.*D*CSD
IF(MPRINT.GE.1)WRITE(lou,102)T
AUX=SM/(4.*PI*ROS)
RETURN
C
22 continue
AMSOUR=0.
DO 23 I=1,3
DO 23 J=1,3
AMSOUR=AMSOUR+h(I)*T(I,J)*p(J)
23 CONTINUE
AMSOUR=AUX*sqsl*AMSOUR
GO TO 51
C
C.....EXPLOSIVE (IMPLOSIVE) SOURCE
C
31 CONTINUE
IF(MWAVE.NE.0)GO TO 32
AUX=4.*PI*ROS
AUX=PAR(1)*FLOAT(IPAR(1))/AUX
RETURN
32 continue
AMSOUR=h(1)*p(1)+h(2)*p(2)+h(3)*p(3)
AMSOUR=AUX*sqsl*AMSOUR
GO TO 51
C
41 CONTINUE
IF(MWAVE.NE.0)GO TO 22
T(1,1)=PAR(1)
T(1,2)=PAR(2)
T(2,1)=T(1,2)
T(1,3)=PAR(3)
T(3,1)=T(1,3)
T(2,2)=PAR(4)
T(2,3)=PAR(5)
T(3,2)=T(2,3)
T(3,3)=PAR(6)
IF(MPRINT.GE.1)WRITE(lou,102)T
AUX=1./(4.*PI*ROS)
RETURN
C
51 PHSOUR=0.
50 IF(AMSOUR)52,53,53
52 AMSOUR=-AMSOUR
PHSOUR=PHSOUR-PI
C
100 FORMAT(4I5,6F10.5)
101 FORMAT(8F10.5)
102 FORMAT(//2X,3(3E15.5/))
C
53 RETURN
END
C                                                                 ```