      SUBROUTINE CONVOL(FUN1,FUN2,NPC)
      PARAMETER (NMAX=4096)
      DIMENSION Y1(NMAX+2),Y2(NMAX+2)
      CHARACTER*64 NA,STATUS
      REAL*8 TPI
      EXTERNAL FUN1,FUN2
      DATA EPSILON/1.D-5/
      DATA TPI/6.283185308D0/
C *** FIND THE BEST SET OF POINTS TO REPRESENT FUN1(X), AND FUN2(X)
C *** NOTE THAT FUN1 AND FUN2 ALREADY CONTAIN A POWER OF R
      RMAX=40.
      H=RMAX/NMAX
      DO I=1,NMAX
        R=H*(I-1)
        Y1(I)=FUN1(R)/(2*TPI) 
        Y2(I)=FUN2(R) 
      ENDDO                                                     
      CALL SINFT(Y1,NMAX)
C *** MAKE THIS INTO A 3D SPHERICALLY SYMMETRIC TRANSFORM OF FUN/R
      AKH=TPI/(2*RMAX)
      NA='FKFUN1N.OUT'
      OPEN(1,FILE=NA)
      NA='FKFUN1A.OUT'
      OPEN(2,FILE=NA)
      DO I=2,NMAX
        AK=(I-1)*AKH
        Y1(I)=2*TPI*H*Y1(I)/AK
        T1A=1/(AK*AK)
        WRITE(1,'(2G16.4)')AK,Y1(I)
        WRITE(2,'(2G16.4)')AK,T1A
C *** CORRECTING THE 1/R TRANSFORM
C       Y1(I)=T1A
      ENDDO
      CLOSE(1)
      CLOSE(2)
      CALL SINFT(Y2,NMAX) 
      NA='FKFUN2N.OUT'
      OPEN(1,FILE=NA)
      NA='FKFUN2A.OUT'
      OPEN(2,FILE=NA)
      DO I=2,NMAX
        AK=(I-1)*AKH
        Y2(I)=2*TPI*H*Y2(I)/AK
        T2A=2*TPI*(SIN(10*AK)-10*AK*COS(10*AK))/AK**3
        WRITE(1,'(2G16.4)')AK,Y2(I)
        WRITE(2,'(2G16.4)')AK,T2A
C       Y2(I)=T2A
      ENDDO
      CLOSE(1)
      CLOSE(2)
C *** FORM AND PLOT THE PRODUCT FORM
      NA='FKPRODN.OUT'
      OPEN(1,FILE=NA)
      NA='FKPRODA.OUT'
      OPEN(2,FILE=NA)
      AKH=TPI/(2*RMAX)
      DO I=2,NMAX
        AK=(I-1)*AKH
        Y1(I)=Y1(I)*Y2(I)*AK
        WRITE(1,'(2G16.4)')AK,Y1(I)
        AK10=AK*10
        ANAL=0
        ANAL=(2*TPI)*(SIN(AK10)-AK10*COS(AK10))/(AK**4)
        WRITE(2,'(2G16.4)')AK,ANAL
C       Y1(I)=ANAL
      ENDDO
      CLOSE(1)
      CLOSE(2)
      Y1(1)=2*Y1(2)-Y1(3)
      CALL SINFT(Y1,NMAX)
C *** PLOT THE CONVOLUTION
      NA='FRCONVN.OUT'
      OPEN(1,FILE=NA)
      NA='FRCONVA.OUT'
      OPEN(2,FILE=NA)
      DO I=2,NMAX
        R=(I-1)*H
        Y1(I)=Y1(I)*2*AKH/(R*TPI*TPI)
        WRITE(1,'(2G16.4)')R,Y1(I)
        IF(R.LT.10D0)THEN
          PHI=50-R*R/6
        ELSE
          PHI=1000/(3*R)
        ENDIF
        WRITE(2,'(2G16.4)')R,PHI
      ENDDO
      CLOSE(1)
      CLOSE(2)
      RETURN
      END
