real transmission

chunkie/+chnk/+axissymhelm2d/green.m

@@ -32,6 +32,13 @@
     asym_tables = chnk.axissymhelm2d.load_asym_tables();
 end
 
+persistent asym_tables_sub   
+if isempty(asym_tables_sub)
+    asym_tables_sub = chnk.axissymhelm2d.load_asym_tables_sub();
+end
+
+htables = asym_tables;
+
 [~, ns] = size(src);
 [~, nt] = size(targ);
 
@@ -45,12 +52,19 @@
     ifun = 2;
 end
 
-if ifdiff
+if (ifdiff == 1)
     if abs(zi) > eps
         error('AXISSYMHELM2D.green: Difference of greens function only supported for real wave numbers\n');
     end
     ifun = 3;
 end
+if (ifdiff == 2)
+    if abs(zi) > eps
+        error('AXISSYMHELM2D.green: Difference of greens function only supported for real wave numbers\n');
+    end
+    htables = asym_tables_sub;
+    ifun = 5;
+end
 
 over2pi = 1/2/pi;
 kabs = abs(k);
@@ -61,7 +75,7 @@
 r  = (rt + origin(1))*kabs;
 dz = dz*kabs;
 dr = (rs-rt)*kabs;
-dout = chnk.axissymhelm2d.helm_axi_all(r, dr, dz, asym_tables,ifun);
+dout = chnk.axissymhelm2d.helm_axi_all(r, dr, dz, htables,ifun);
 int = dout.int;
 intdz = dout.intdz;
 intdq = dout.intdq;

chunkie/+chnk/+axissymhelm2d/helm_axi_all.m

@@ -5,10 +5,10 @@
 
     int   = zeros(size(alphs));
 
-    iflag_rk = (ifun == 1 || ifun == 4);
+    iflag_rk = (ifun == 1 || ifun == 4 || ifun == 5);
     iflag_ik = (ifun == 2 || ifun == 4);
-    iflag_dk = (ifun == 3 || ifun == 4);
-
+    iflag_dk = (ifun == 3 || ifun == 4 || ifun == 5);
+    
     if (iflag_rk)
     rk = cell(6);
     for ii=1:6
@@ -186,6 +186,10 @@
         [dsdiff] = proc_kerns(rs,drs,dzs,dk);
         varargout{1} = dsdiff;
     end
+    if(ifun==5)
+        [dsdiff] = proc_kerns(rs,drs,dzs,dk);
+        varargout{1} = dsdiff;
+    end
 
 end
 

chunkie/+chnk/+axissymhelm2d/helm_axi_smooth.m

@@ -68,7 +68,7 @@
     s{6} = kernsdaa;
     s{4} = kernsdkk;
 
-    if (ifun == 1 || ifun==4)
+    if (ifun == 1 || ifun==4 || ifun==5)
         sout.rk = s;
     end
     if (ifun ==2)
@@ -106,6 +106,36 @@
         sout.dk = dk;
     end
 
+
+    if (ifun == 5)
+        r0t = 0;
+        efac = exp(1i*r0t.*rts)./rts.*wle;
+        r2   = rts.^2;
+        k2   = 1i*r0t.*rts;
+        kern2_v = sum(asymint_v(rts,r0t,efac),1);
+        kern2dk_v = sum(asymintdk_v(rts,r0t,efac),1);
+        kern2da_v = sum(asymintda_v(xle,rts,r0t,efac),1);
+        kern2daa_v = sum(asymintdaa_v(xle,rts,r0t,efac,k2,r2),1);
+        kern2dak_v = sum(asymintdak_v(xle,rts,r0t,efac),1);
+        kern2dkk_v = sum(asymintdkk_v(xle,rts,r0t,efac),1);
+
+        ik = {};
+        ik{1} = kern2_v;
+        ik{2} = kern2da_v;
+        ik{3} = 0*kern2dk_v;
+        ik{6} = kern2daa_v;
+        ik{5} = 0*kern2dak_v;
+        ik{4} = 0*kern2dkk_v;
+        sout.k0 = ik;
+
+        dk = {};
+        for ii=1:6
+            dk{ii} = sout.rk{ii}-sout.k0{ii};
+        end
+        sout.dk = dk;
+    end
+
+
 end
 
 function [val] = asymint_v(r,k,efac)

chunkie/+chnk/+axissymhelm2d/kern.m

@@ -1,4 +1,4 @@
-function submat = kern(zk, srcinfo, targinfo, origin, type, varargin)
+function submat = kern(zks, srcinfo, targinfo, origin, type, varargin)
 %CHNK.AXISSYMHELM2D.KERN axissymmetric Helmholtz layer potential kernels in 2D
 % 
 % Syntax: submat = chnk.axissymhelm2d.kern(zk,srcinfo,targingo,type,htables)
@@ -61,6 +61,13 @@
 [~, ns] = size(src);
 [~, nt] = size(targ);
 
+if (numel(zks) == 1)
+    zk = zks(1);
+else
+    zk1 = zks(1);
+    zk2 = zks(2);
+end
+
 if strcmpi(type, 'd')
     srcnorm = srcinfo.n;
     [~, grad] = chnk.axissymhelm2d.green(zk, src, targ, origin);
@@ -183,6 +190,31 @@
               - hess(:,:,6).*nxsrc.*nytarg + hess(:,:,3).*nysrc.*nytarg;
 end
 
+if strcmpi(type, 'dprime_re_diff')
+
+  targnorm = targinfo.n;
+  srcnorm = srcinfo.n;
+  ifdiff = 2;
+  [~,~,hess] = chnk.axissymhelm2d.green(zk1, src, targ, origin,ifdiff);
+  nxsrc = repmat(srcnorm(1,:),nt,1);
+  nysrc = repmat(srcnorm(2,:),nt,1);
+  nxtarg = repmat((targnorm(1,:)).',1,ns);
+  nytarg = repmat((targnorm(2,:)).',1,ns);
+  submat1 = hess(:,:,4).*nxsrc.*nxtarg - hess(:,:,5).*nysrc.*nxtarg ... 
+              - hess(:,:,6).*nxsrc.*nytarg + hess(:,:,3).*nysrc.*nytarg;
+  ifdiff = 2;
+  [~,~,hess] = chnk.axissymhelm2d.green(zk2, src, targ, origin,ifdiff);
+  nxsrc = repmat(srcnorm(1,:),nt,1);
+  nysrc = repmat(srcnorm(2,:),nt,1);
+  nxtarg = repmat((targnorm(1,:)).',1,ns);
+  nytarg = repmat((targnorm(2,:)).',1,ns);
+  submat2 = hess(:,:,4).*nxsrc.*nxtarg - hess(:,:,5).*nysrc.*nxtarg ... 
+              - hess(:,:,6).*nxsrc.*nytarg + hess(:,:,3).*nysrc.*nytarg;
+
+  submat = submat1-submat2;
+
+end
+
 
 
 if strcmpi(type, 'dprimediff')

chunkie/+chnk/+axissymhelm2d/load_asym_tables_sub.m

@@ -0,0 +1,95 @@
+function asym_tables = load_asym_tables_sub()
+%CHNK.AXISSYMHELM2D.load_asym_tables loads the precomputed tables
+% for axissymmetric Helmholtz kernels
+    p = mfilename('fullpath');
+    dirname = dir([p ,'.m']).folder;
+    fname = [dirname '/asym_helm_data2.mat'];
+    load(fname,'allvs');
+    msizes = size(allvs);
+    ncheb = msizes(1);
+    nks = msizes(5);
+    nas = msizes(4);
+    nkers = msizes(3);
+    xcheb = (0:(ncheb-1))/(ncheb-1)*pi;
+
+    [N,X]=meshgrid(0:(ncheb-1),xcheb);
+    Tc2v = cos(N.*X);
+    Tv2c = inv(Tc2v);
+    Tfull = kron(Tv2c,Tv2c);
+    Tc2vd = N.*sin(N.*X)./sqrt(1-cos(X).^2);
+    Tc2vd(1,:) = (0:(ncheb-1)).^2;
+    Tc2vd(end,:) = (-1).^(1:ncheb).*(0:(ncheb-1)).^2;
+    Tc2cd = Tv2c*Tc2vd;
+
+    xcheb = cos(xcheb)';
+
+
+    allda = allvs;
+    alldk = allvs;
+    alldaa = allvs;
+    alldak = allvs;
+    alldkk = allvs;
+
+    for kk=1:nks
+        for jj=1:nas
+            for ii=1:nkers
+%                vmat = squeeze(allvs(:,:,ii,jj,kk));
+                if (ii >=1)
+                    vmat = squeeze(allvs(:,:,ii,jj,kk));
+                else
+                    vmat = squeeze(allvs(:,:,ii,jj,kk));
+                    vmat = reshape(Tfull*vmat(:),[12,12]);
+                    allvs(:,:,ii,jj,kk) = vmat;
+                end
+                ia = jj;
+                vda = (2)^(ia-1)*20*Tc2cd*vmat;
+                vdaa= (2)^(ia-1)*20*Tc2cd*vda;
+                vdk = 2/(pi)*Tc2cd*vmat.';
+                vdkk= (2)/pi*Tc2cd*vdk;
+                vdak = 2/(pi)*Tc2cd*vda.';
+                allda(:,:,ii,jj,kk) = vda;
+                alldaa(:,:,ii,jj,kk)= vdaa;
+                alldk(:,:,ii,jj,kk) = vdk.';
+                alldak(:,:,ii,jj,kk)= vdak.';
+                alldkk(:,:,ii,jj,kk)= vdkk.';
+            end
+        end
+    end    
+
+    asym_tables =[];
+    asym_tables.allvs = allvs;
+    asym_tables.allda = allda;
+    asym_tables.alldaa= alldaa;
+    asym_tables.alldk = alldk;
+    asym_tables.alldak= alldak;
+    asym_tables.alldkk= alldkk;
+    asym_tables.ncheb = ncheb;
+
+    nlege = 500;
+    [xlege,wlege,~,~] = lege.exps(nlege);
+    xlege = (pi*(xlege+1)/2);
+    wlege = wlege*pi/2;
+    asym_tables.xlege = xlege;
+    asym_tables.wlege = wlege;
+
+   	nlege = 100;
+    [xlege,wlege,~,~] = lege.exps(nlege);
+    xlege = (pi*(xlege+1)/2);
+    wlege = wlege*pi/2;
+    asym_tables.xlege_mid = xlege;
+    asym_tables.wlege_mid = wlege;
+
+  	nlege = 50;
+    [xlege,wlege,~,~] = lege.exps(nlege);
+    xlege = (pi*(xlege+1)/2);
+    wlege = wlege*pi/2;
+    asym_tables.xlege_midnear = xlege;
+    asym_tables.wlege_midnear = wlege;
+
+  	nlege = 20;
+    [xlege,wlege,~,~] = lege.exps(nlege);
+    xlege = (pi*(xlege+1)/2);
+    wlege = wlege*pi/2;
+    asym_tables.xlege_midnearnear = xlege;
+    asym_tables.wlege_midnearnear = wlege;
+end

chunkie/+chnk/+axissymhelm2d/test1.m

@@ -0,0 +1,96 @@
+k = 10.1;
+src = [1;2];
+targ = [1.01;2];
+origin = [0;0];
+ifdiff = 2;
+[val, grad, hess] = chnk.axissymhelm2d.green(k, src, targ, origin, ifdiff);
+
+dx = targ-src;
+dr = dx(1);
+dz = dx(2);
+rs = src(1);
+
+r0   = sqrt(rs.^2+(rs+dr).^2+dz.^2);
+alph = (dr.^2+dz.^2)./r0.^2;
+alph = 1-alph;
+frt = @(x) sqrt(1-alph*cos(x));
+fex = @(x) exp(1i*k*r0*frt(x));
+fdi = @(x) (fex(x)-1)./frt(x)/r0;
+
+qint = integral(fdi,0,2*pi,'AbsTol',1E-12);
+darea = rs;
+qint = qint*darea/(4*pi);
+
+dh = 1E-5;
+
+vph = gfunc(src,targ.*[1+dh;1],k);
+vmh = gfunc(src,targ.*[1-dh;1],k);
+vdt1 = (vph-vmh)/(2*dh*targ(1))*darea
+
+grad(1) - vdt1
+
+vph = gfunc(src.*[1+dh;1],targ,k);
+vmh = gfunc(src.*[1-dh;1],targ,k);
+vds1 = (vph-vmh)/(2*dh*src(1))*darea
+
+grad(2) - vds1
+
+vph = gfunc(src,targ.*[1;1+dh],k);
+vmh = gfunc(src,targ.*[1;1-dh],k);
+vdt2 = (vph-vmh)/(2*dh*targ(2))*darea
+
+grad(3) - vdt2
+
+ds = zeros(4,4);
+
+svec= [src;targ];
+for ii=1:4
+    for jj=1:4
+        s = svec;
+        s(ii) = s(ii)*(1+dh);
+        s(jj) = s(jj)*(1+dh);
+        sv = s(1:2);
+        tv = s(3:4);
+        vpp = gfunc(sv,tv,k);
+        s = svec;
+        s(ii) = s(ii)*(1-dh);
+        s(jj) = s(jj)*(1+dh);
+        sv = s(1:2);
+        tv = s(3:4);
+        vmp = gfunc(sv,tv,k); 
+        s = svec;
+        s(ii) = s(ii)*(1+dh);
+        s(jj) = s(jj)*(1-dh);
+        sv = s(1:2);
+        tv = s(3:4);
+        vpm = gfunc(sv,tv,k); 
+        s = svec;
+        s(ii) = s(ii)*(1-dh);
+        s(jj) = s(jj)*(1-dh);
+        sv = s(1:2);
+        tv = s(3:4);
+        vmm = gfunc(sv,tv,k); 
+        ds(ii,jj) = (vpp-vmp+vmm-vpm)/(4*dh^2*svec(ii)*svec(jj));
+    end
+end
+
+
+function [val] = gfunc(s_in,t_in,k)
+
+dx = t_in-s_in;
+dr = dx(1);
+dz = dx(2);
+rs = s_in(1);
+
+r0   = sqrt(rs.^2+(rs+dr).^2+dz.^2);
+alph = (dr.^2+dz.^2)./r0.^2;
+alph = 1-alph;
+frt = @(x) sqrt(1-alph*cos(x));
+fex = @(x) exp(1i*k*r0*frt(x));
+fdi = @(x) (fex(x)-1)./frt(x)/r0;
+
+qint = integral(fdi,0,2*pi,'AbsTol',1E-12);
+darea = 1;
+val = qint*darea/(4*pi);
+
+end

chunkie/@kernel/axissymhelm2ddiff.m

@@ -27,8 +27,12 @@
 zr1 = real(zks(1)); zi1 = imag(zks(1));
 zr2 = real(zks(2)); zi2 = imag(zks(2));
 
+ifreal = 0;
 if zi1 ~=0 || zr2 ~=0 || zr1 ~= zi2
-    error('Wave numbers must be of the form zk, 1i*zk with zk real');
+    ifreal = 1;
+    if (zi1~= 0 || zi2 ~= 0)
+        error('Wave numbers must be of the form zk, 1i*zk with zk real');
+    end
 end
 
 obj = kernel();
@@ -89,11 +93,19 @@
             error('Coefs must be of form [c c]');   
         end
         obj.type = 'dp';
+        if(~ifreal)
         obj.eval = @(s,t) coefs(1)*chnk.axissymhelm2d.kern(real(zks(1)),  ... 
                                                 s, t, [0,0], 'dprimediff');
         obj.shifted_eval = @(s,t,o) coefs(1)*chnk.axissymhelm2d.kern(real(zks(1)),  ... 
                                                 s, t, o, 'dprimediff');
         obj.fmm = [];
+        else
+        obj.eval = @(s,t) coefs(1)*chnk.axissymhelm2d.kern(zks,  ... 
+                                                s, t, [0,0], 'dprime_re_diff');
+        obj.shifted_eval = @(s,t,o) coefs(1)*chnk.axissymhelm2d.kern(zks,  ... 
+                                                s, t, o, 'dprime_re_diff');
+        obj.fmm = [];
+        end
         if ( abs(coefs(1)-coefs(2)) < eps )
             obj.sing = 'log';
         else