test ising

Author: tristangdwl

chunkie/+chnk/+helm2dquas/green.m

@@ -100,8 +100,8 @@
 val_near= zeros(nkappa,nptclose);
 grad_near = zeros(nkappa,nptclose,2);
 hess_near = zeros(nkappa,nptclose,3);
+ls = -l:l;
 if ~isempty(rxclose)
-    ls = -l:l;
     for i = ls
         if ising == 1
             iuse = true(nptclose,1);
@@ -204,9 +204,11 @@
     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);
@@ -217,11 +219,13 @@
     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);
@@ -234,6 +238,7 @@
     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);
@@ -242,6 +247,7 @@
         val(:,isub,:) = val(:,isub,:) - vali;
         grad(:,isub,:) = grad(:,isub,:) - gradi;
         hess(:,isub,:) = hess(:,isub,:) - hessi;
+        end
     end
 
     val = reshape(val,nkappa*ntarg,nsrc);

chunkie/@kernel/helm2dquas.m

@@ -88,42 +88,42 @@
 quas_param.sn = sn;
 
 obj.params.quas_param = quas_param;
-
+obj.params.ising = ising;
 switch lower(type)
 
     case {'s', 'single'}
         obj.type = 's';
-        obj.eval = @(s,t) chnk.helm2dquas.kern(zk, s, t, 's',quas_param,ising);
+        obj.eval = @(s,t) chnk.helm2dquas.kern(zk, s, t, 's',quas_param,[],ising);
         obj.fmm = [];
         obj.sing = 'log';
         if isscalar(kappa)
         obj.splitinfo = [];
         obj.splitinfo.type = {[0 0 0 0],[1 0 0 0]};
         obj.splitinfo.action = {'r','r'};
-        obj.splitinfo.functions = @(s,t) helm2dquas_s_split(zk,s,t,quas_param,ising);
+        obj.splitinfo.functions = @(s,t) helm2dquas_s_split(zk,s,t,quas_param);
         end
 
     case {'d', 'double'}
         obj.type = 'd';
-        obj.eval = @(s,t) chnk.helm2dquas.kern(zk, s, t, 'd',quas_param,ising);
+        obj.eval = @(s,t) chnk.helm2dquas.kern(zk, s, t, 'd',quas_param,[],ising);
         obj.fmm = [];
         obj.sing = 'log';
         if isscalar(kappa)
         obj.splitinfo = [];
         obj.splitinfo.type = {[0 0 0 0],[1 0 0 0],[0 0 -1 0]};
         obj.splitinfo.action = {'r','r','r'};
-        obj.splitinfo.functions = @(s,t) helm2dquas_d_split(zk,s,t,quas_param,ising);
+        obj.splitinfo.functions = @(s,t) helm2dquas_d_split(zk,s,t,quas_param);
         end
 
     case {'sp', 'sprime'}
         obj.type = 'sp';
-        obj.eval = @(s,t) chnk.helm2dquas.kern(zk, s, t, 'sprime',quas_param,ising);
+        obj.eval = @(s,t) chnk.helm2dquas.kern(zk, s, t, 'sprime',quas_param,[],ising);
         obj.fmm = [];
         obj.sing = 'log';
 
     case {'dp', 'dprime'}
         obj.type = 'dp';
-        obj.eval = @(s,t) chnk.helm2dquas.kern(zk, s, t, 'dprime',quas_param,ising);
+        obj.eval = @(s,t) chnk.helm2dquas.kern(zk, s, t, 'dprime',quas_param,[],ising);
         obj.fmm = [];
         obj.sing = 'hs';
 
@@ -141,7 +141,7 @@
         obj.splitinfo = [];
         obj.splitinfo.type = {[0 0 0 0],[1 0 0 0],[0 0 -1 0]};
         obj.splitinfo.action = {'r','r','r'};
-        obj.splitinfo.functions = @(s,t) helm2dquas_c_split(zk,s,t,quas_param,coefs,ising);
+        obj.splitinfo.functions = @(s,t) helm2dquas_c_split(zk,s,t,quas_param,coefs);
         end
 
     case {'cp', 'cprime'}
@@ -159,6 +159,10 @@
         error('Unknown Helmholtz kernel type ''%s''.', type);
 
 end
+if obj.params.ising == 0
+    obj.sing = 'smooth'; 
+    obj.splitinfo = [];
+end
 end
 
 function f = helm2dquas_s_split(zk,s,t,quas_param)

devtools/test/quasiperiodicTest.m

@@ -87,6 +87,25 @@ function quasiperiodicTest0()
 c2trval_assemb(1:2, :) = coefs(1)*dval + coefs(2)*sval;
 c2trval_assemb(3:end, :) = coefs(1)*dpval + coefs(2)*spval;
 assert(norm(c2trval_assemb-c2trval) <  1e-12)
+
+
+skern1a = kernel('hq','s',zk,kappa,d,[],[],0);
+skern1b = kernel('h','s',zk);
+sval1 = skern1a.eval(src,targ) + skern1b.eval(src,targ);
+assert(norm(sval-sval1) <  1e-12)
+
+dkern1a = kernel('hq','d',zk,kappa,d,[],[],0);
+dkern1b = kernel('h','d',zk);
+dval1 = dkern1a.eval(src,targ) + dkern1b.eval(src,targ);
+assert(norm(dval-dval1) <  1e-12)
+
+allkern1a = kernel(@(s,t) chnk.helm2dquas.kern(zk,s,t,'all',quas_param,coefa,0));
+allkern1b = kernel(@(s,t) chnk.helm2d.kern(zk,s,t,'all',coefa));
+allval1a = allkern1a.eval(src,targ);
+allval1b = allkern1b.eval(src,targ); 
+allval1b = reshape(allval1b,1,2,2); allval1b = repmat(allval1b,length(kappa),1,1);
+allval1b = reshape(allval1b,length(kappa)*2,2);
+assert(norm(allval1a+allval1b -allval) <  1e-12)
 end