complex u(3,3)
        real k,mu
        common/prop1/h,k,mu,eta
        write (*,*) 'Enter k:'
        read (*,*) k
        h=0.5*k
        mu=1.0
        xc=0.
        yc=0.
        zc=0.
        write (*,*) 'Enter x,y,z:'
        read (*,*) x,y,z
        call gf(x,y,z,xc,yc,zc,u)
        write (*,'(6f10.4)') u(1,1),u(1,2),u(1,3)
        write (*,'(6f10.4)') u(2,1),u(2,2),u(2,3)
        write (*,'(6f10.4)') u(3,1),u(3,2),u(3,3)
        end

        subroutine gf(x,y,z,xc,yc,zc,u)
        real k,mu
        common/prop1/h,k,mu,eta
        complex u(3,3)
        complex uu(3)
        xx=x-xc
        yy=y-yc
        zz=z-zc
        call gfr(xx,yy,zz,uu)
        u(1,3)=uu(1)
        u(2,3)=uu(2)
        u(3,3)=uu(3)
        call gfr(zz,yy,-xx,uu)
        u(1,1)=uu(3)
        u(2,1)=uu(2)
        u(3,1)=-uu(1)
        call gfr(zz,xx,yy,uu)
        u(1,2)=uu(2)
        u(2,2)=uu(3)
        u(3,2)=uu(1)
        return
        end

        subroutine gfr(x,y,z,uu)
        real h,h2,k,k2,mu
        common/prop1/h,k,mu,eta
        complex uu(3)
        complex eh,ek,expa,expb
        complex dxeh,dyeh,dzeh
        complex dxzeh,dyzeh,dzzeh
        complex dxek,dyek,dzek
        complex dxzek,dyzek,dzzek
        complex factor

        pi=3.141593

        r=sqrt(x*x+y*y+z*z)
        r2=r*r
        r3=r2*r
        r4=r2*r2
        write (*,*) 'r=',r,r2,r3,r4
        x2=x*x
        y2=y*y
        z2=z*z
        h2=h*h
        k2=k*k

        gk=1./(4.*pi*mu*k2)
        expa=cexp((0.,-1.)*h*r)
        expb=cexp((0.,-1.)*k*r)
        eh=expa/r
        ek=expb/r

        factor=cmplx(1./r , h)
        dxeh=-(x/r2)*factor*expa
        dyeh=-(y/r2)*factor*expa
        dzeh=-(z/r2)*factor*expa
        factor=cmplx(3./r2-h2 , 3.*h/r)
        dxzeh=(x*z/r3)*factor*expa
        dyzeh=(y*z/r3)*factor*expa
        factor=cmplx(3.*z2/r3-(h2*z2+1.)/r , -h+3.*h*z2/r2)
        dzzeh=(1./r2)*factor*expa

        factor=cmplx(1./r , k)
        dxek=-(x/r2)*factor*expb
        dyek=-(y/r2)*factor*expb
        dzek=-(z/r2)*factor*expb
        factor=cmplx(3./r2-k2 , 3.*k/r)
        dxzek=(x*z/r3)*factor*expb
        dyzek=(y*z/r3)*factor*expb
        factor=((3.*z2-r2)/r)*(1./r+(0.,1.)*k)-k2*z2
        dzzek=(1./r3)*factor*expb

        uu(1)=gk*(dxzek-dxzeh)
        uu(2)=gk*(dyzek-dyzeh)
        uu(3)=gk*(dzzek-dzzeh+k2*ek)
        return
        end