C       revised Nov. 20, 1996 by petra

C       This subroutine calculates the density at points X2
C       using points X.  It is returned in dens1.  The number
C       of variables in X is given by ncol1. The bandwidth is
C       passed in hd1.
C       The program expects to receive R, the interquartile range
C       of the data which it uses in calculating the silverman
C       optimal bandwidth formula for a biweight kernel

	SUBROUTINE CALCDENS (Xm,X2m,n,n2,R,dens1,hd1)
	
C	Variable Declarations

	integer n, n2, j, o
	double precision Xm(n), X2m(n2), hd , R
	double precision dist(n), K(n),dens1(n2)
	double precision arg,sum1,k2,k3,k4 

C----------------------------------------------------------------
C       Calculate the optimal bandwidth by the Silverman formula
C       adjusted for the biweight kernel

c       Silverman rule of thumb bandwidth using interquartile range R

	hd = hd*(R/1.34)*(n**(-1.0/5.0))


	do 2 o=1,n2

c	------------------------------------------------------------
c	Calculate K((X2_o - X(j)) / h) using quartic (biweight) kernel function 


	arg = 0.0
	do 140 j=1,n
          K(j) = 0.0
	  arg =  ((X2m(o)-Xm(j)) / hd) 
          if(abs(arg) .LE. 1.0) then
              K(j) = (15.0/16.0)*(arg**2.0-1.0)**2.0
          end if
140 	continue

C       calculate the density as the sum of the products of the
C       kernels

        sum1 = 0.0
        do 150 j = 1,n
          sum1 = sum1+K(j)
150     continue
       
        dens1(o) = (1/(n*hd))*sum1

2  	continue


	RETURN
	END