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.EQ.0, READ THE INPUT DATA FOR RADIATION PATTERNS
C     FOR MWAVE.NE.0, COMPUTE THE RADIATION PATTERNS
C
c     READ THE INPUT DATA
c
      PI=4.*ATAN(1.)
      IF(MWAVE.NE.0) GO TO 1
      if(mtype.eq.0)return
      READ(lin,100) IPAR,PAR
      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     COMPUTE THE RADIATION PATTERNS
C
      GO TO (11,21,31),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
   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