real*8 k,h,mu,eta
        common/prop1/h,k,mu,eta
        complex uu(2,2),tt(2,2),sum(2)
        complex t(2),u(2)
        real*8 x,y,z
        write (*,*) 'Enter wave number k:'
        read (*,*) k
        h=k/2.
        mu=1.0

10      write (*,*) 'Enter x,y:'
        read (*,*) x,y

        z=0.d0
        call gfr(x,y,z)
        go to 10
11      stop
        end

        subroutine gfr(x,y,z)
        real*8 h,k,mu,eta
        common/prop1/h,k,mu,eta
        complex*16 uu(3),ss(3,3)
        complex*16 eh,ek,expa
        complex*16 expb,dxeh,dyeh,dzeh
        complex*16 dxzeh,dyzeh,dzzeh
        complex*16 dxxzeh,dyyzeh,dzzzeh,dxyzeh,dxzzeh,dyzzeh
        complex*16 dxek,dyek,dzek
        complex*16 dxzek,dyzek,dzzek
        complex*16 dxxzek,dyyzek,dzzzek,dxyzek,dxzzek,dyzzek
        complex*16 factor,f1,f2

        complex*16 dyzphi,dyzzpsi,dypsi
        complex*16 dzzzek1
        complex*16 uqq


        real*8 x,y,z,r,r2,r3,r4,x2,y2,z2,h2,k2,pi,gk
        pi=3.1415926535d0

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

        gk=1./(4.*pi*mu*k2)
c        expa=cexpd((0.,-1.)*h*r)
c        expb=cexpd((0.,-1.)*k*r)
        expa=dcmplx(dcos(h*r),-dsin(h*r))
        expb=dcmplx(dcos(k*r),-dsin(k*r))

        eh=expa/r
        ek=expb/r

        factor=dcmplx(1./r , h)
        dxeh=-(x/r2)*factor*expa
        dyeh=-(y/r2)*factor*expa
        dzeh=-(z/r2)*factor*expa
        factor=dcmplx(3./r2-h2 , 3.*h/r)
        dxzeh=(x*z/r3)*factor*expa
        dyzeh=(y*z/r3)*factor*expa
        factor=dcmplx(3.*z2/r3-(h2*z2+1.)/r , -h+3.*h*z2/r2)
        dzzeh=(1./r2)*factor*expa
        f1=-(15.*x*z/r3)*(1./r+(0.,1.)*h*x)+(3.*z/r)*(1./r+(0.,1.)*h)
        f2=(h2*x2*z/r)*(6./r+(0.,1.)*h)-h2*z
        dxxzeh=(f1+f2)/r3*expa
        f1=-(15.*y*z/r3)*(1./r+(0.,1.)*h*y)+(3.*z/r)*(1./r+(0.,1.)*h)
        f2=(h2*y2*z/r)*(6./r+(0.,1.)*h)-h2*z
        dyyzeh=(f1+f2)*expa/r3
        f1=-(15.*z/r3)*(1./r+(0.,1.)*h*z2)+(9.*z/r)*(1./r+(0.,1.)*h)
        f2=(h2*z2*z/r)*(6./r+(0.,1.)*h)-3.*h2*z
        dzzzeh=(f1+f2)/r3*expa
        factor=-(15./r2)*(1./r+(0.,1.)*h)+h2*(6./r+(0.,1.)*h)
        dxyzeh=(x*y*z/r4)*factor*expa
        f1=((3.*r2-15.*z2)/r3)*(1./r+(0.,1.)*h)
        f2=(h2*z2/r)*(6./r+(0.,1.)*h)-h2
        dxzzeh=(x/r3)*(f1+f2)*expa
        dyzzeh=(y/r3)*(f1+f2)*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
        f1=-(15.*x*z/r3)*(1./r+(0.,1.)*k*x)+(3.*z/r)*(1./r+(0.,1.)*k)
        f2=(k2*x2*z/r)*(6./r+(0.,1.)*k)-k2*z
        dxxzek=(f1+f2)/r3*expb
        f1=-(15.*y*z/r3)*(1./r+(0.,1.)*k*y)+(3.*z/r)*(1./r+(0.,1.)*k)
        f2=(k2*y2*z/r)*(6./r+(0.,1.)*k)-k2*z
        dyyzek=(f1+f2)*expb/r3
        f1=-(15.*z/r3)*(1./r+(0.,1.)*k*z2)+(9.*z/r)*(1./r+(0.,1.)*k)
        f2=(k2*z2*z/r)*(6./r+(0.,1.)*k)-3.*k2*z
        dzzzek=(f1+f2)/r3*expb
        factor=-(15./r2)*(1./r+(0.,1.)*k)+k2*(6./r+(0.,1.)*k)
        dxyzek=(x*y*z/r4)*factor*expb
        f1=((3.*r2-15.*z2)/r3)*(1./r+(0.,1.)*k)
        f2=(k2*z2/r)*(6./r+(0.,1.)*k)-k2
        dxzzek=(x/r3)*(f1+f2)*expb
        dyzzek=(y/r3)*(f1+f2)*expb

        uu(1)=gk*(dxzek-dxzeh)
        uu(2)=gk*(dyzek-dyzeh)
        uu(3)=gk*(dzzek-dzzeh+k2*ek)
        uqq=r*uu(3)
        write (*,'(2f15.10)') uqq

        ss(1,1)=gk*(k2*dzeh+2.*dyyzeh+2.*dzzzeh+2.*dxxzek)
        ss(1,2)=2.*gk*(dxyzek-dxyzeh)
        ss(1,3)=gk*(2.*(dxzzek-dxzzeh)+k2*dxek)
        ss(2,2)=gk*(k2*dzeh+2.*dxxzeh+2.*dzzzeh+2.*dyyzek)
        ss(2,3)=gk*(2.*(dxzzek-dxzzeh)+k2*dyek)
        ss(3,3)=gk*(k2*(dzeh+2.*dzek)+2.*dxxzeh+2.*dyyzeh+2.*dzzzek)
        ss(2,1)=ss(1,2)
        ss(3,1)=ss(1,3)
        ss(3,2)=ss(2,3)
        return
        end