Update shape of vector valued quasiperiodic kernels

chunkie/+chnk/+helm2d/kern.m

@@ -40,7 +40,7 @@
 %                       [coef(1,1)*D coef(1,2)*S; coef(2,1)*D' coef(2,2)*S']
 %                type == 'c2trans' returns the combined field, and the 
 %                          normal derivative of the combined field
-%                        [coef(1)*D + coef(2)*S; coef(1)*D' + coef(2)*S']
+%                        [coef(1,1)*D + coef(1,2)*S; coef(2,1)*D' + coef(2,2)*S']
 %                type == 'trans_rep' returns the potential corresponding
 %                           to the transmission representation
 %                        [coef(1)*D coef(2)*S]
@@ -250,6 +250,7 @@
 if strcmpi(type,'c2trans') 
   coef = ones(2,1);
   if(nargin == 5); coef = varargin{1}; end
+  if numel(coef) == 2, coef = repmat(coef(:).',2,1); end
   targnorm = targinfo.n(:,:);
   srcnorm = srcinfo.n(:,:);
 
@@ -272,8 +273,8 @@
 
   submat = zeros(2*nt,ns);
 
-  submat(1:2:2*nt,:) = coef(1)*submatd + coef(2)*submats;
-  submat(2:2:2*nt,:) = coef(1)*submatdp + coef(2)*submatsp;
+  submat(1:2:2*nt,:) = coef(1,1)*submatd + coef(1,2)*submats;
+  submat(2:2:2*nt,:) = coef(2,1)*submatdp + coef(2,2)*submatsp;
 
 end
 

chunkie/+chnk/+helm2dquas/green.m

@@ -1,5 +1,5 @@
-function [val,grad,hess] = green(src,targ,zk,kappa,d,sn,l)
-%CHNK.HELM2DQUAS.GREEN evaluate the quasiperiodic helmholtz green's function
+function [val,grad,hess] = green(src,targ,zk,kappa,d,sn,l,ising)
+%CHNK.HELM2DQUAS.GREEN evaluate the quasiperiodic Helmholtz Green's function
 % for the given sources and targets
 %
 % Input:
@@ -8,12 +8,14 @@
 %   d - period
 %   sn - precomputed lattice sum integrals 
 %       (see chnk.helm2dquas.latticecoefs)
+%   l - number of periodic copies computed explicitly
+%   ising - if set to 0, only include the periodic copies. If set to 1,
+%       include the free-space part
 %
 % see also CHNK.HELM2DQUAS.KERN
 [~,nsrc] = size(src);
 [~,ntarg] = size(targ);
 
-
 xs = repmat(src(1,:),ntarg,1);
 ys = repmat(src(2,:),ntarg,1);
 
@@ -40,8 +42,7 @@
 
 npt = size(r,1);
 
-ythresh = d/2;
-% ythresh = d/10;
+ythresh = 2*d/2;
 iclose = abs(ry) < ythresh;
 ifar = ~iclose;
 
@@ -52,6 +53,9 @@
 ryclose = ry(iclose);
 rclose = r(iclose);
 
+nxclose = nx(iclose);
+nptclose = size(rxclose, 1);
+
 nkappa = length(kappa);
 
 val = zeros(nkappa,npt,1);
@@ -76,55 +80,57 @@
 % beta = sqrt((xi_m.^2-zk^2));
 beta = sqrt(1i*(xi_m-zk)).*sqrt(-1i*(xi_m+zk));
 
-% fhat = exp(-beta.*sqrt(ryfar.^2))./beta.*exp(1i*xi_m.*rxfar)/2;
 fhat = exp(-beta.*sqrt(ryfar.^2) + 1i*xi_m.*rxfar)./(2*beta);
 val(:,ifar,:) = sum(fhat,3)/(d);
 if nargout > 1
 grad(:,ifar,1) = sum(1i*xi_m.*fhat,3)/d;
 grad(:,ifar,2) = sum(-beta.*(sqrt(ryfar.^2)./ryfar).*fhat,3)/d;
-% grad(:,ifar,:) =[gx,gy];
 end
 
 if nargout >2
 hess(:,ifar,1) = sum(-xi_m.^2.*fhat,3)/d;
 hess(:,ifar,2) = sum(-1i*xi_m.*beta.*(sqrt(ryfar.^2)./ryfar).*fhat,3)/d;
 hess(:,ifar,3) = sum((beta.*(sqrt(ryfar.^2)./ryfar)).^2.*fhat,3)/d;
-% hess(:,ifar,:) = [hxx,hxy,hyy];
 end
 
-
 end
 
 alpha = (exp(1i*kappa(:)*d));
 
-val_near=0;
-grad_near = zeros(1,1,2);
-hess_near = zeros(1,1,3);
+val_near= zeros(nkappa,nptclose);
+grad_near = zeros(nkappa,nptclose,2);
+hess_near = zeros(nkappa,nptclose,3);
+ls = -l:l;
 if ~isempty(rxclose)
-    for i = -l:l
+    for i = ls
+        if ising == 1
+            iuse = true(nptclose,1);
+        else
+            iuse = nxclose ~= -i;
+        end
+
         rxi = rxclose - i*d;
         if nargout>2
         [vali,gradi,hessi] = chnk.helm2d.green(zk,[0;0],[rxi.';ryclose.']);
         vali = reshape(vali,1,[],1);
         gradi = reshape(gradi,1,[],2);
         hessi = reshape(hessi,1,[],3);
-        val_near = val_near + vali.*alpha.^i;
-        grad_near = grad_near + gradi.*alpha.^i;
-        hess_near = hess_near + hessi.*alpha.^i;
+        val_near(:,iuse) = val_near(:,iuse) + vali(:,iuse).*alpha.^i;
+        grad_near(:,iuse,:) = grad_near(:,iuse,:) + gradi(:,iuse,:).*alpha.^i;
+        hess_near(:,iuse,:) = hess_near(:,iuse,:) + hessi(:,iuse,:).*alpha.^i;
         elseif nargout > 1
         [vali,gradi] = chnk.helm2d.green(zk,[0;0],[rxi.';ryclose.']);
         vali = reshape(vali,1,[],1);
         gradi = reshape(gradi,1,[],2);
-        val_near = val_near + vali.*alpha.^i;
-        grad_near = grad_near + gradi.*alpha.^i;
+        val_near(:,iuse) = val_near(:,iuse) + vali(:,iuse).*alpha.^i;
+        grad_near(:,iuse,:) = grad_near(:,iuse,:) + gradi(:,iuse,:).*alpha.^i;
         else
         vali = chnk.helm2d.green(zk,[0;0],[rxi.';ryclose.']);
         vali = reshape(vali,1,[],1);
-        val_near = val_near + vali.*alpha.^i;
+        val_near(:,iuse) = val_near(:,iuse) + vali(:,iuse).*alpha.^i;
         end
     end
-
-    % sn = sn.';
+
     N = size(sn,2)-1;
     ns = (0:N);
     ns_use = (0:N+2);
@@ -149,9 +155,6 @@
     val_far = 0.25*1i*Js(:,:,1).*sn(:,:,1) + 0.5*1i*sum(sn(:,:,2:end).*Js(:,:,2:end-2).*cs(:,:,2:end),3);
     val(:,iclose) = val_near+val_far;
 
-
-
-
     if nargout >1
         DJs = cat(3,-Js(:,:,2),.5*(Js(:,:,1:end-3)-Js(:,:,3:end-1)))*zk;
         ss = (eipn-1./eipn)/2i;
@@ -188,21 +191,67 @@
 
         hess_far = cat(3,hess_far_xx, hess_far_xy, hess_far_yy);
 
-
         hess(:,iclose,:) = hess_near + hess_far;
     end
 end
 
 quasi_phase = exp(1i*kappa(:)*nx(:).'*d);
 
-val = reshape(quasi_phase.*val,nkappa*ntarg,nsrc);
-if nargout>1
-grad = reshape(quasi_phase.*grad,nkappa*ntarg,nsrc,2);
-end
-if nargout>2
-hess = reshape(quasi_phase.*hess,nkappa*ntarg,nsrc,3);
-end
-% end
 
-end
+if nargout == 1
+    val = quasi_phase.*val;
 
+    if ising == 0
+        isub = (abs(nx(:)) > max(ls)) | ifar;
+
+        if any(isub)
+        vali = chnk.helm2d.green(zk,[0;0],[rx(isub).' + nx(isub).'*d;ry(isub).']);
+        vali = reshape(vali,1,[],1);
+        val(:,isub,:) = val(:,isub,:) - vali;
+        end
+    end
+
+    val = reshape(val,nkappa*ntarg,nsrc);
+elseif nargout == 2
+    val = quasi_phase.*val;
+    grad = quasi_phase.*grad;
+
+    if ising == 0
+        isub = (abs(nx(:)) > max(ls)) | ifar;
+
+        if any(isub)
+        [vali, gradi] = chnk.helm2d.green(zk,[0;0],[rx(isub).' + nx(isub).'*d;ry(isub).']);
+        vali = reshape(vali,1,[],1);
+        gradi = reshape(gradi,1,[],2);
+        val(:,isub,:) = val(:,isub,:) - vali;
+        grad(:,isub,:) = grad(:,isub,:) - gradi;
+        end
+    end
+
+    val = reshape(val,nkappa*ntarg,nsrc);
+    grad = reshape(grad,nkappa*ntarg,nsrc,2);
+elseif nargout == 3
+    val = quasi_phase.*val;
+    grad = quasi_phase.*grad;
+    hess = quasi_phase.*hess;
+
+    if ising == 0
+        isub = (abs(nx(:)) > max(ls)) | ifar;
+
+        if any(isub)
+        [vali, gradi, hessi] = chnk.helm2d.green(zk,[0;0],[rx(isub).' + nx(isub).'*d;ry(isub).']);
+        vali = reshape(vali,1,[],1);
+        gradi = reshape(gradi,1,[],2);
+        hessi = reshape(hessi,1,[],3);
+
+        val(:,isub,:) = val(:,isub,:) - vali;
+        grad(:,isub,:) = grad(:,isub,:) - gradi;
+        hess(:,isub,:) = hess(:,isub,:) - hessi;
+        end
+    end
+
+    val = reshape(val,nkappa*ntarg,nsrc);
+    grad = reshape(grad,nkappa*ntarg,nsrc,2);
+    hess = reshape(hess,nkappa*ntarg,nsrc,3);
+end
+end