function ffs=polk(kx1,dkx,nkx,rho,KY,ff);

[m,nky]=size(ff);


i_l=8;
i_r=8192;
xintp=intp_coef(i_l,i_r);

kxx=(rho(:).^2*ones(1,nky))-(ones(m,1)*KY).^2;
kxx=kxx.*(kxx > 0); kxx=sqrt(kxx);

  dkx0=dkx/i_r;

  ffs=zeros(nkx,nky);

  for i=1:m;
i
   fkx=kxx(i,:)-kx1;
   ikx=floor(fkx./(dkx*ones(1,nky)));
   dif=fkx-(dkx*ikx);
   ikx=ikx+1;
   ikx=ones(i_l,1)*ikx+[-floor(i_l/2)+1:floor(i_l/2)]'*ones(1,nky);
   ikx=max(ikx,ones(i_l,nky));
   ikx=min(ikx,nkx*ones(i_l,nky));
   ikx=ikx+nkx*ones(i_l,1)*[0:nky-1];
   ikx0=round(dif/dkx0)+1;

   ffs(ikx)=ffs(ikx)+ ((ones(i_l,1)*ff(i,:)).*xintp(:,ikx0));

  end;