From 1967c48eb82e451813af97d63c60dfb75fe1d837 Mon Sep 17 00:00:00 2001 From: jghoskins <32715900+jghoskins@users.noreply.github.com> Date: Fri, 24 May 2024 10:27:41 -0500 Subject: [PATCH] updated kernels --- .../+chnk/+flex2d/besseldiff_etc_pscoefs.m | 32 + chunkie/+chnk/+flex2d/even_pseval.m | 17 + chunkie/+chnk/+flex2d/fmm.m | 121 ++ chunkie/+chnk/+flex2d/green.m | 52 + chunkie/+chnk/+flex2d/helmdiffgreen.m | 178 +++ chunkie/+chnk/+flex2d/kern.m | 1322 +++++++++++++++++ .../+chnk/+flex2d/kernel_replacement_test.m | 44 + .../smooth_part_derivatives.m | 0 chunkie/+chnk/+flex2d/test_helmdiffgreen.m | 80 + chunkie/+chnk/+flex2d/test_ps_besseldiff.m | 16 + .../transmission_helper.m | 0 .../+chnk/+helm2d/thirdforth_derivatives.m | 178 --- 12 files changed, 1862 insertions(+), 178 deletions(-) create mode 100644 chunkie/+chnk/+flex2d/besseldiff_etc_pscoefs.m create mode 100644 chunkie/+chnk/+flex2d/even_pseval.m create mode 100644 chunkie/+chnk/+flex2d/fmm.m create mode 100644 chunkie/+chnk/+flex2d/green.m create mode 100644 chunkie/+chnk/+flex2d/helmdiffgreen.m create mode 100644 chunkie/+chnk/+flex2d/kern.m create mode 100644 chunkie/+chnk/+flex2d/kernel_replacement_test.m rename chunkie/+chnk/{+helm2d => +flex2d}/smooth_part_derivatives.m (100%) create mode 100644 chunkie/+chnk/+flex2d/test_helmdiffgreen.m create mode 100644 chunkie/+chnk/+flex2d/test_ps_besseldiff.m rename chunkie/+chnk/{+helm2d => +flex2d}/transmission_helper.m (100%) delete mode 100644 chunkie/+chnk/+helm2d/thirdforth_derivatives.m diff --git a/chunkie/+chnk/+flex2d/besseldiff_etc_pscoefs.m b/chunkie/+chnk/+flex2d/besseldiff_etc_pscoefs.m new file mode 100644 index 0000000..7d91028 --- /dev/null +++ b/chunkie/+chnk/+flex2d/besseldiff_etc_pscoefs.m @@ -0,0 +1,32 @@ +function [cf1,cf2] = besseldiff_etc_pscoefs(nterms) +% powerseries coefficients for J_0 - 1 and the other series in +% definition of Y_0 +% +% these are both series with the terms z^2,z^4,z^6,..z^(2*nterms) +% +% cf1 are power series coeffs of even terms (excluding constant) +% for J_0(z) - 1 +% +% cf2 are power series coeffs for the similar series +% z^2/4 - (1+1/2)*(z^2/4)^2/(2!)^2 + (1+1/2+1/3)*(z^2/4)^3/(3!)^2 - ... +% +% which appears in the power series for H_0(z) + +sum1 = 0; +fac1 = 1; +pow1 = 1; + +cf1 = zeros(nterms,1); +cf2 = zeros(nterms,1); +sgn = 1; + +for j = 1:nterms + fac1 = fac1/(j*j); + sum1 = sum1 + 1.0/j; + pow1 = pow1*0.25; + cf1(j) = -sgn*pow1*fac1; + cf2(j) = sgn*sum1*pow1*fac1; + sgn = -sgn; +end + +end \ No newline at end of file diff --git a/chunkie/+chnk/+flex2d/even_pseval.m b/chunkie/+chnk/+flex2d/even_pseval.m new file mode 100644 index 0000000..70be417 --- /dev/null +++ b/chunkie/+chnk/+flex2d/even_pseval.m @@ -0,0 +1,17 @@ +function [v] = even_pseval(cf,x) +% evaluates an even power series with no constant term +% with coefficients given in cf. +% +% x can be an array + +x2 = x.*x; + +xp = x2; + +v = zeros(size(x)); + +nterms = length(cf); +for j = 1:nterms + v = v + cf(j)*xp; + xp = xp.*x2; +end diff --git a/chunkie/+chnk/+flex2d/fmm.m b/chunkie/+chnk/+flex2d/fmm.m new file mode 100644 index 0000000..b62111a --- /dev/null +++ b/chunkie/+chnk/+flex2d/fmm.m @@ -0,0 +1,121 @@ +function varargout = fmm(eps, zk, srcinfo, targinfo, type, sigma, varargin) +%CHNK.HELM2D.FMM Fast multipole methods for evaluating Helmholtz layer +%potentials, their gradients, and Hessians. +% +% Syntax: +% pot = chnk.helm2d.fmm(eps, zk, srcinfo, targinfo, type, sigma, varargin) +% [pot, grad] = chnk.helm2d.fmm(eps, zk, srcinfo, targinfo, type, sigma, varargin) +% [pot, grad, hess] = chnk.helm2d.fmm(eps, zk, srcinfo, targinfo, type, sigma, varargin) +% +% Let x be targets and y be sources for these formulas, with +% n_x and n_y the corresponding unit normals at those points +% (if defined). +% +% Kernels based on G(x,y) = i/4 H_0^{(1)}(zk |x-y|) +% +% D(x,y) = \nabla_{n_y} G(x,y) +% S(x,y) = G(x,y) +% +% Note: the density should be scaled by the weights +% +% Input: +% eps - precision requested +% zk - complex number, Helmholtz wave number +% srcinfo - description of sources in ptinfo struct format, i.e. +% ptinfo.r - positions (2,:) array +% ptinfo.d - first derivative in underlying +% parameterization (2,:) +% ptinfo.d2 - second derivative in underlying +% parameterization (2,:) +% ptinfo.n - normal (2,:) array +% +% targinfo - description of targets in ptinfo struct format, +% type - string, determines kernel type +% type == 'd', double layer kernel D +% type == 's', single layer kernel S +% type == 'c', combined layer kernel coefs(1)*D + coefs(2)*S +% type = 'sprime' normal derivative of single layer, +% note that targinfo must contain targ.n +% type = 'dprime' normal derivative of double layer, +% note that targinfo must contain targ.n +% sigma - density +% varargin{1} - coef in the combined layer formula, otherwise +% does nothing +% +% Optional Output: +% pot - potential/neumann data corresponding to the kernel at the target locations +% grad - gradient at target locations +% hess - hessian at target locations + +if ( nargout == 0 ) + warning('CHUNKIE:helm2d:fmm:empty', ... + 'Nothing to compute in HELM2D.FMM. Returning empty array.'); + return +end + +srcuse = []; +srcuse.sources = srcinfo.r(1:2,:); +switch lower(type) + case {'s', 'sprime'} + srcuse.charges = sigma(:).'; + case {'d', 'dprime'} + srcuse.dipstr = sigma(:).'; + srcuse.dipvec = srcinfo.n(1:2,:); + case 'c' + coefs = varargin{1}; + srcuse.charges = coefs(2)*sigma(:).'; + srcuse.dipstr = coefs(1)*sigma(:).'; + srcuse.dipvec = srcinfo.n(1:2,:); +end + +if ( isstruct(targinfo) ) + targuse = targinfo.r(:,:); +else + targuse = targinfo; +end + +pg = 0; +pgt = min(nargout, 2); +U = hfmm2d(eps, zk, srcuse, pg, targuse, pgt); + +% Assign potentials +if ( nargout > 0 ) + switch lower(type) + case {'s', 'd', 'c'} + varargout{1} = U.pottarg.'; + case {'sprime', 'dprime'} + if ( ~isfield(targinfo, 'n') ) + error('CHUNKIE:helm2d:fmm:normals', ... + 'Targets require normal info when evaluating Helmholtz kernel ''%s''.', type); + end + varargout{1} = ( U.gradtarg(1,:).*targinfo.n(1,:) + ... + U.gradtarg(2,:).*targinfo.n(2,:) ).'; + otherwise + error('CHUNKIE:lap2d:fmm:pot', ... + 'Potentials not supported for Helmholtz kernel ''%s''.', type); + end +end + +% Assign gradients +if ( nargout > 1 ) + switch lower(type) + case {'s', 'd', 'c'} + varargout{2} = U.gradtarg; + otherwise + error('CHUNKIE:helm2d:fmm:grad', ... + 'Gradients not supported for Helmholtz kernel ''%s''.', type); + end +end + +% Assign Hessians +if ( nargout > 2 ) + switch lower(type) + case {'s', 'd', 'c'} + varargout{3} = U.hesstarg; + otherwise + error('CHUNKIE:helm2d:fmm:hess', ... + 'Hessians not supported for Helmholtz kernel ''%s''.', type); + end +end + +end diff --git a/chunkie/+chnk/+flex2d/green.m b/chunkie/+chnk/+flex2d/green.m new file mode 100644 index 0000000..04fdb2f --- /dev/null +++ b/chunkie/+chnk/+flex2d/green.m @@ -0,0 +1,52 @@ +function [val,grad,hess] = green(k,src,targ) +%CHNK.HELM2D.GREEN evaluate the helmholtz green's function +% for the given sources and targets + +[~,ns] = size(src); +[~,nt] = size(targ); + +xs = repmat(src(1,:),nt,1); +ys = repmat(src(2,:),nt,1); + +xt = repmat(targ(1,:).',1,ns); +yt = repmat(targ(2,:).',1,ns); + +rx = xt-xs; +ry = yt-ys; + +rx2 = rx.*rx; +ry2 = ry.*ry; + +r2 = rx2+ry2; + +r = sqrt(r2); + + +if nargout > 0 + h0 = besselh(0,1,k*r); + val = 0.25*1i*h0; +end + +[m,n] = size(xs); + +if nargout > 1 + h1 = besselh(1,1,k*r); + grad = zeros(m,n,2); + + grad(:,:,1) = -1i*k*0.25*h1.*rx./r; + grad(:,:,2) = -1i*k*0.25*h1.*ry./r; + +end + +if nargout > 2 + + hess = zeros(m,n,3); + + r3 = r.^3; + + h2 = 2*h1./(k*r)-h0; + + hess(:,:,1) = 0.25*1i*k*((rx-ry).*(rx+ry).*h1./r3 - k*rx2.*h0./r2); + hess(:,:,2) = 0.25*1i*k*k*rx.*ry.*h2./r2; + hess(:,:,3) = 0.25*1i*k*((ry-rx).*(rx+ry).*h1./r3 - k*ry2.*h0./r2); +end diff --git a/chunkie/+chnk/+flex2d/helmdiffgreen.m b/chunkie/+chnk/+flex2d/helmdiffgreen.m new file mode 100644 index 0000000..cf9489b --- /dev/null +++ b/chunkie/+chnk/+flex2d/helmdiffgreen.m @@ -0,0 +1,178 @@ +function [val,grad,hess,der3,der4] = helmdiffgreen(k,src,targ) +%HELMDIFFGREEN evaluate the difference of the +% Helmholtz Green function and the Laplace Green function +% for the given sources and targets, i.e. +% +% G(x,y) = i/4 H_0^(1)(k|x-y|) + 1/(2 pi) log(|x-y|) +% +% where H_0^(1) is the principal branch of the Hankel function +% of the first kind. This routine avoids numerical cancellation +% when |k||x-y| is small. +% +% - grad(:,:,1) has G_{x1}, grad(:,:,2) has G_{x2} +% - hess(:,:,1) has G_{x1x1}, hess(:,:,2) has G_{x1x2}, +% hess(:,:,3) has G_{x2x2} +% - der3 has the third derivatives in the order G_{x1x1x1}, G_{x1x1x2}, +% G_{x1x2x2}, G_{x2x2x2} +% - der4 has the fourth derivatives in the order G_{x1x1x1x1}, +% G_{x1x1x1x2}, G_{x1x1x2x2}, G_{x1x2x2x2}, G_{x2x2x2x2} +% +% derivatives are on the *target* variables + +[~,ns] = size(src); +[~,nt] = size(targ); + +xs = repmat(src(1,:),nt,1); +ys = repmat(src(2,:),nt,1); + +xt = repmat(targ(1,:).',1,ns); +yt = repmat(targ(2,:).',1,ns); + +dx = xt-xs; +dy = yt-ys; + +dx2 = dx.*dx; +dx3 = dx2.*dx; +dx4 = dx3.*dx; +dy2 = dy.*dy; +dy3 = dy2.*dy; +dy4 = dy3.*dy; + +r2 = dx2 + dy2; +r = sqrt(r2); +rm1 = 1./r; +rm2 = rm1.*rm1; +rm3 = rm1.*rm2; +rm4 = rm1.*rm3; +rm5 = rm1.*rm4; + +% get value and r derivatives + +[g0,g1,g2,g3,g4] = diff_h0log_and_rders(k,r); + +% evaluate potential and derivatives + +if nargout > 0 + val = g0; +end +if nargout > 1 + grad(:,:,1) = dx.*g1.*rm1; + grad(:,:,2) = dy.*g1.*rm1; +end +if nargout > 2 + hess(:,:,1) = dx2.*g2.*rm2+g1.*(1.*rm1-dx2.*rm3); + hess(:,:,2) = dx.*dy.*(g2.*rm2-g1.*rm3); + hess(:,:,3) = dy2.*g2.*rm2+g1.*(1.*rm1-dy2.*rm3); +end +if nargout > 3 + der3(:,:,1) = (dx3.*g3+3*dy2.*dx.*(g2.*rm1-g1.*rm2)).*rm3; + der3(:,:,2) = dx2.*dy.*(g3.*rm3-3*(g2.*rm4-g1.*rm5)) + ... + dy.*(g2.*rm2-g1.*rm3); + der3(:,:,3) = dx.*dy2.*(g3.*rm3-3*(g2.*rm4-g1.*rm5)) + ... + dx.*(g2.*rm2-g1.*rm3); + der3(:,:,4) = (dy3.*g3+3*dx2.*dy.*(g2.*rm1-g1.*rm2)).*rm3; +end + +if nargout > 4 + der4(:,:,1) = (dx4.*(g4-6*g3.*rm1+15*(g2.*rm2-g1.*rm3))).*rm4 + ... + (6*dx2.*(g3-3*(g2.*rm1-g1.*rm2))).*rm3 + ... + (3*(g2-g1.*rm1)).*rm2; + der4(:,:,2) = (dx3.*dy.*(g4-6*g3.*rm1+15*(g2.*rm2-g1.*rm3))).*rm4 + ... + (3*dx.*dy.*(g3-3*(g2.*rm1-g1.*rm2))).*rm3; + der4(:,:,3) = dx2.*dy2.*(g4-6*g3.*rm1+15*g2.*rm2-15*g1.*rm3).*rm4 + ... + g3.*rm1 - 2*g2.*rm2 + 2*g1.*rm3; + der4(:,:,4) = dx.*dy3.*(g4-6*g3.*rm1+15*(g2.*rm2-g1.*rm3)).*rm4 + ... + 3*dx.*dy.*(g3-3*(g2.*rm1-g1.*rm2)).*rm3; + der4(:,:,5) = dy4.*(g4-6*g3.*rm1+15*(g2.*rm2-g1.*rm3)).*rm4 + ... + 6*dy2.*(g3-3*(g2.*rm1-g1.*rm2)).*rm3 + ... + 3*(g2-g1.*rm1).*rm2; +end + +end + +function [g0,g1,g2,g3,g4] = diff_h0log_and_rders(k,r) + +g0 = zeros(size(r)); +g1 = zeros(size(r)); +g2 = zeros(size(r)); +g3 = zeros(size(r)); +g4 = zeros(size(r)); + +io4 = 1i*0.25; +o2p = 1/(2*pi); + +isus = abs(k)*r < 1; +%isus = false(size(r)); + +% straightforward formulas for sufficiently large + +rnot = r(~isus); + +kr = k*rnot; + +h0 = besselh(0,1,kr); +h1 = besselh(1,1,kr); + +rm1 = 1./rnot; +rm2 = rm1.*rm1; +rm3 = rm1.*rm2; +rm4 = rm1.*rm3; + +g0(~isus) = io4*h0 + o2p*log(rnot); +g1(~isus) = -k*io4*h1 + o2p*rm1; +g2(~isus) = -k*k*io4*h0 + k*io4*h1.*rm1 - o2p*rm2; +g3(~isus) = k*k*io4*h0.*rm1 + io4*k*(k*k-2*rm2).*h1 + 2*o2p*rm3; +g4(~isus) = k*io4*(3*rm2-k*k).*(2*h1.*rm1-k*h0) - 6*o2p*rm4; + +% manually cancel when small + +rsus = r(isus); +rm1 = 1./rsus; +rm2 = rm1.*rm1; +rm3 = rm1.*rm2; +rm4 = rm1.*rm3; + +gam = 0.57721566490153286060651209; +nterms = 14; +const1 = (io4+(log(2)-gam-log(k))*o2p); + +% relevant parts of hankel function represented as power series +[cf1,cf2] = chnk.flex2d.besseldiff_etc_pscoefs(nterms); +kpow = (k*k).^(1:nterms); +cf1 = cf1(:).*kpow(:); cf2 = cf2(:).*kpow(:); + +logr = log(rsus); + +j0m1 = chnk.flex2d.even_pseval(cf1,rsus); +f = chnk.flex2d.even_pseval(cf2,rsus); + +% differentiate power series to get derivatives +fac = 2*(1:nterms); +cf1 = cf1.*fac(:); cf2 = cf2.*fac(:); +j0m1d1 = chnk.flex2d.even_pseval(cf1,rsus).*rm1; +fd1 = chnk.flex2d.even_pseval(cf2,rsus).*rm1; +cf1 = cf1.*(fac(:)-1); cf2 = cf2.*(fac(:)-1); +j0m1d2 = chnk.flex2d.even_pseval(cf1,rsus).*rm2; + +fd2 = chnk.flex2d.even_pseval(cf2,rsus).*rm2; +cf1 = cf1(:).*(fac(:)-2); cf1 = cf1(2:end); +cf2 = cf2(:).*(fac(:)-2); cf2 = cf2(2:end); +j0m1d3 = chnk.flex2d.even_pseval(cf1,rsus).*rm1; + +fd3 = chnk.flex2d.even_pseval(cf2,rsus).*rm1; +fac = fac(1:end-1); +cf1 = cf1(:).*(fac(:)-1); cf2 = cf2(:).*(fac(:)-1); +j0m1d4 = chnk.flex2d.even_pseval(cf1,rsus).*rm2; +fd4 = chnk.flex2d.even_pseval(cf2,rsus).*rm2; + +% combine to get derivative of i/4 H + log/(2*pi) +g0(isus) = const1*(j0m1+1) - o2p*(f + logr.*j0m1); +g1(isus) = const1*j0m1d1 - o2p*(fd1 + logr.*j0m1d1 + j0m1.*rm1); +g2(isus) = const1*j0m1d2 - o2p*(fd2 + logr.*j0m1d2 + 2*j0m1d1.*rm1 - j0m1.*rm2); +g3(isus) = const1*j0m1d3 - o2p*(fd3 + logr.*j0m1d3 + 3*j0m1d2.*rm1 - ... + 3*j0m1d1.*rm2 + 2*j0m1.*rm3); +g4(isus) = const1*j0m1d4 - o2p*(fd4 + logr.*j0m1d4 + 4*j0m1d3.*rm1 - ... + 6*j0m1d2.*rm2 + 8*j0m1d1.*rm3 - 6*j0m1.*rm4); + +end + diff --git a/chunkie/+chnk/+flex2d/kern.m b/chunkie/+chnk/+flex2d/kern.m new file mode 100644 index 0000000..0893d6f --- /dev/null +++ b/chunkie/+chnk/+flex2d/kern.m @@ -0,0 +1,1322 @@ +function submat= kern(zk,srcinfo,targinfo,type,varargin) +%CHNK.HELM2D.KERN standard Helmholtz layer potential kernels in 2D +% +% Syntax: submat = chnk.heml2d.kern(zk,srcinfo,targingo,type,varargin) +% +% Let x be targets and y be sources for these formulas, with +% n_x and n_y the corresponding unit normals at those points +% (if defined). Note that the normal information is obtained +% by taking the perpendicular to the provided tangential deriviative +% info and normalizing +% +% Kernels based on G(x,y) = i/4 H_0^{(1)}(zk |x-y|) +% +% D(x,y) = \nabla_{n_y} G(x,y) +% S(x,y) = G(x,y) +% S'(x,y) = \nabla_{n_x} G(x,y) +% D'(x,y) = \nabla_{n_x} \nabla_{n_y} G(x,y) +% +% Input: +% zk - complex number, Helmholtz wave number +% srcinfo - description of sources in ptinfo struct format, i.e. +% ptinfo.r - positions (2,:) array +% ptinfo.d - first derivative in underlying +% parameterization (2,:) +% ptinfo.d2 - second derivative in underlying +% parameterization (2,:) +% targinfo - description of targets in ptinfo struct format, +% if info not relevant (d/d2) it doesn't need to +% be provided. sprime requires tangent info in +% targinfo.d +% type - string, determines kernel type +% type == 'd', double layer kernel D +% type == 's', single layer kernel S +% type == 'sprime', normal derivative of single +% layer S' +% type == 'dprime', normal derivative of double layer D' +% type == 'c', combined layer kernel coef(1) D + coef(2) S +% type == 'stau', tangential derivative of single layer +% type == 'all', returns all four layer potentials, +% [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'] +% type == 'trans_rep' returns the potential corresponding +% to the transmission representation +% [coef(1)*D coef(2)*S] +% type == 'trans_rep_prime' returns the normal derivative +% corresponding to the transmission representation +% [coef(1)*D' coef(2)*S'] +% type == 'trans_rep_grad' returns the gradient corresponding +% to the transmission representation +% [coef(1)*d_x D coef(2)*d_x S; +% coef(1)*d_y D coef(2)*d_y S] +% +% varargin{1} - coef: length 2 array in the combined layer +% formula, 2x2 matrix for all kernels +% otherwise does nothing +% +% Output: +% submat - the evaluation of the selected kernel for the +% provided sources and targets. the number of +% rows equals the number of targets and the +% number of columns equals the number of sources +% +% see also CHNK.HELM2D.GREEN + +src = srcinfo.r; +targ = targinfo.r; + +[~,ns] = size(src); +[~,nt] = size(targ); + +% double layer +if strcmpi(type,'d') + srcnorm = srcinfo.n; + [~,grad] = chnk.flex2d.green(zk,src,targ); + nx = repmat(srcnorm(1,:),nt,1); + ny = repmat(srcnorm(2,:),nt,1); + submat = -(grad(:,:,1).*nx + grad(:,:,2).*ny); +end + +% normal derivative of single layer +if strcmpi(type,'sprime') + targnorm = targinfo.n; + [~,grad] = chnk.flex2d.green(zk,src,targ); + nx = repmat((targnorm(1,:)).',1,ns); + ny = repmat((targnorm(2,:)).',1,ns); + submat = (grad(:,:,1).*nx + grad(:,:,2).*ny); +end + + +% Tangential derivative of single layer +if strcmpi(type,'stau') + targtan = targinfo.d; + [~,grad] = chnk.flex2d.green(zk,src,targ); + dx = repmat((targtan(1,:)).',1,ns); + dy = repmat((targtan(2,:)).',1,ns); + ds = sqrt(dx.*dx+dy.*dy); + submat = (grad(:,:,1).*dx + grad(:,:,2).*dy)./ds; +end + +% single layer +if strcmpi(type,'s') + submat = chnk.flex2d.green(zk,src,targ); +end + +% normal derivative of double layer +if strcmpi(type,'dprime') + targnorm = targinfo.n; + srcnorm = srcinfo.n; + [~,~,hess] = chnk.flex2d.green(zk,src,targ); + nxsrc = repmat(srcnorm(1,:),nt,1); + nysrc = repmat(srcnorm(2,:),nt,1); + nxtarg = repmat((targnorm(1,:)).',1,ns); + nytarg = repmat((targnorm(2,:)).',1,ns); + submat = -(hess(:,:,1).*nxsrc.*nxtarg + hess(:,:,2).*(nysrc.*nxtarg+nxsrc.*nytarg)... + + hess(:,:,3).*nysrc.*nytarg); +end + +% Combined field +if strcmpi(type,'c') + srcnorm = srcinfo.n; + coef = ones(2,1); + if(nargin == 5); coef = varargin{1}; end + [submats,grad] = chnk.flex2d.green(zk,src,targ); + nx = repmat(srcnorm(1,:),nt,1); + ny = repmat(srcnorm(2,:),nt,1); + submatd = -(grad(:,:,1).*nx + grad(:,:,2).*ny); + submat = coef(1)*submatd + coef(2)*submats; +end + +% normal derivative of combined field +if strcmpi(type,'cprime') + coef = ones(2,1); + if(nargin == 5); coef = varargin{1}; end; + targnorm = targinfo.n; + srcnorm = srcinfo.n; + + submat = zeros(nt,ns); + + % Get gradient and hessian info + [submats,grad,hess] = chnk.flex2d.green(zk,src,targ); + + nxsrc = repmat(srcnorm(1,:),nt,1); + nysrc = repmat(srcnorm(2,:),nt,1); + nxtarg = repmat((targnorm(1,:)).',1,ns); + nytarg = repmat((targnorm(2,:)).',1,ns); + + % D' + submatdp = -(hess(:,:,1).*nxsrc.*nxtarg ... + + hess(:,:,2).*(nysrc.*nxtarg+nxsrc.*nytarg)... + + hess(:,:,3).*nysrc.*nytarg); + % S' + submatsp = (grad(:,:,1).*nxtarg + grad(:,:,2).*nytarg); + + submat = coef(1)*submatdp + coef(2)*submatsp; +end + + +% Dirichlet and neumann data corresponding to combined field +if strcmpi(type,'c2trans') + coef = ones(2,1); + if(nargin == 5); coef = varargin{1}; end; + targnorm = targinfo.n; + srcnorm = srcinfo.n; + + % Get gradient and hessian info + [submats,grad,hess] = chnk.flex2d.green(zk,src,targ); + + nxsrc = repmat(srcnorm(1,:),nt,1); + nysrc = repmat(srcnorm(2,:),nt,1); + nxtarg = repmat((targnorm(1,:)).',1,ns); + nytarg = repmat((targnorm(2,:)).',1,ns); + + % D' + submatdp = -(hess(:,:,1).*nxsrc.*nxtarg ... + + hess(:,:,2).*(nysrc.*nxtarg+nxsrc.*nytarg)... + + hess(:,:,3).*nysrc.*nytarg); + % S' + submatsp = (grad(:,:,1).*nxtarg + grad(:,:,2).*nytarg); + % D + submatd = -(grad(:,:,1).*nxsrc + grad(:,:,2).*nysrc); + + 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; + +end + + +% all kernels, [c11 D, c12 S; c21 D', c22 S'] +if strcmpi(type,'all') + + targnorm = targinfo.n; + srcnorm = srcinfo.n; + cc = varargin{1}; + + submat = zeros(2*nt,2*ns); + + % Get gradient and hessian info + [submats,grad,hess] = chnk.flex2d.green(zk,src,targ); + + nxsrc = repmat(srcnorm(1,:),nt,1); + nysrc = repmat(srcnorm(2,:),nt,1); + nxtarg = repmat((targnorm(1,:)).',1,ns); + nytarg = repmat((targnorm(2,:)).',1,ns); + + submatdp = -(hess(:,:,1).*nxsrc.*nxtarg ... + + hess(:,:,2).*(nysrc.*nxtarg+nxsrc.*nytarg)... + + hess(:,:,3).*nysrc.*nytarg); + % S' + submatsp = (grad(:,:,1).*nxtarg + grad(:,:,2).*nytarg); + % D + submatd = -(grad(:,:,1).*nxsrc + grad(:,:,2).*nysrc); + + + submat(1:2:2*nt,1:2:2*ns) = submatd*cc(1,1); + submat(1:2:2*nt,2:2:2*ns) = submats*cc(1,2); + submat(2:2:2*nt,1:2:2*ns) = submatdp*cc(2,1); + submat(2:2:2*nt,2:2:2*ns) = submatsp*cc(2,2); +end + + +% clamped plate kernel for modified biharmonic problem, K21 Kernel is +% handled separately. +if strcmpi(type, 'clamped-plate') + srcnorm = srcinfo.n; + srctang = srcinfo.d; + targnorm = targinfo.n; + + [hess, third, forth] = chnk.flex2d.thirdforth_derivatives(zk, src, targ); % Hankel part + nx = repmat(srcnorm(1,:),nt,1); + ny = repmat(srcnorm(2,:),nt,1); + + nxtarg = repmat((targnorm(1,:)).',1,ns); + nytarg = repmat((targnorm(2,:)).',1,ns); + + zkimag = (1i)*zk; + [hessK, thirdK, forthK] = chnk.flex2d.thirdforth_derivatives(zkimag, src, targ); % modified bessel K part + + dx = repmat(srctang(1,:),nt,1); + dy = repmat(srctang(2,:),nt,1); + + ds = sqrt(dx.*dx+dy.*dy); + + taux = dx./ds; + tauy = dy./ds; + + Kxx = 1/(2*zk^2).*(third(:, :, 1).*(nx.*nx.*nx) + third(:, :, 2).*(3*nx.*nx.*ny) +... + third(:, :, 4).*(3*nx.*ny.*ny) + third(:, :, 6).*(ny.*ny.*ny)) - ... + 1/(2*zk^2).*(thirdK(:, :, 1).*(nx.*nx.*nx) + thirdK(:, :, 2).*(3*nx.*nx.*ny) +... + thirdK(:, :, 4).*(3*nx.*ny.*ny) + thirdK(:, :, 6).*(ny.*ny.*ny)) + ... + (3/(2*zk^2).*(third(:, :, 1).*(nx.*taux.*taux) + third(:, :, 2).*(2*nx.*taux.*tauy + ny.*taux.*taux) +... + third(:, :, 4).*(nx.*tauy.*tauy + 2*ny.*taux.*tauy) + third(:, :, 6).*(ny.*tauy.*tauy))-... + 3/(2*zk^2).*(thirdK(:, :, 1).*(nx.*taux.*taux) + thirdK(:, :, 2).*(2*nx.*taux.*tauy + ny.*taux.*taux) +... + thirdK(:, :, 4).*(nx.*tauy.*tauy + 2*ny.*taux.*tauy) + thirdK(:, :, 6).*(ny.*tauy.*tauy))); % G_{ny ny ny} + 3G_{ny tauy tauy} + + Kxy = -(1/(2*zk^2).*(hess(:, :, 1).*(nx.*nx) + hess(:, :, 2).*(2*nx.*ny) + hess(:, :, 3).*(ny.*ny))-... + 1/(2*zk^2).*(hessK(:, :, 1).*(nx.*nx) + hessK(:, :, 2).*(2*nx.*ny) + hessK(:, :, 3).*(ny.*ny)))+... + (1/(2*zk^2).*(hess(:, :, 1).*(taux.*taux) + hess(:, :, 2).*(2*taux.*tauy) + hess(:, :, 3).*(tauy.*tauy))-... + 1/(2*zk^2).*(hessK(:, :, 1).*(taux.*taux) + hessK(:, :, 2).*(2*taux.*tauy) + hessK(:, :, 3).*(tauy.*tauy))); % -G_{ny ny} + G_{tauy tauy} + + Kyx = 0; + Kyy = 1/(2*zk^2).*(third(:,:, 1).*(nx.*nx.*nxtarg) +third(:, :, 2).*(nx.*nx.*nytarg + 2*nx.*ny.*nxtarg) + third(:, :, 4).*(2*nx.*ny.*nytarg + ny.*ny.*nxtarg)+... + third(:, :,6).*(ny.*ny.*nytarg)) - 1/(2*zk^2).*(thirdK(:,:, 1).*(nx.*nx.*nxtarg) +thirdK(:, :, 2).*(nx.*nx.*nytarg + 2*nx.*ny.*nxtarg) + thirdK(:, :, 4).*(2*nx.*ny.*nytarg + ny.*ny.*nxtarg)+... + thirdK(:, :,6).*(ny.*ny.*nytarg)) - 1/(2*zk^2).*(third(:,:, 1).*(taux.*taux.*nxtarg) +third(:, :, 2).*(taux.*taux.*nytarg + 2*taux.*tauy.*nxtarg) + third(:, :, 4).*(2*taux.*tauy.*nytarg + tauy.*tauy.*nxtarg)+... + third(:, :,6).*(tauy.*tauy.*nytarg)) + 1/(2*zk^2).*(thirdK(:,:, 1).*(taux.*taux.*nxtarg) +thirdK(:, :, 2).*(taux.*taux.*nytarg + 2*taux.*tauy.*nxtarg) + thirdK(:, :, 4).*(2*taux.*tauy.*nytarg + tauy.*tauy.*nxtarg)+... + thirdK(:, :,6).*(tauy.*tauy.*nytarg)); + + submat = zeros(2*nt,2*ns); + + submat(1:2:end,1:2:end) = Kxx; + submat(1:2:end,2:2:end) = Kxy; + + submat(2:2:end,1:2:end) = Kyx; + submat(2:2:end,2:2:end) = Kyy; +end + +if strcmpi(type, 'clamped-plate K21') + srcnorm = srcinfo.n; + srctang = srcinfo.d; + targnorm = targinfo.n; + R = 0.001; + + [~, ~, forth] = chnk.flex2d.thirdforth_derivatives(zk, src, targ); % Hankel part + nx = repmat(srcnorm(1,:),nt,1); + ny = repmat(srcnorm(2,:),nt,1); + + nxtarg = repmat((targnorm(1,:)).',1,ns); + nytarg = repmat((targnorm(2,:)).',1,ns); + + zkimag = (1i)*zk; + [~, ~, forthK] = chnk.flex2d.thirdforth_derivatives(zkimag, src, targ); % modified bessel K part + + dx = repmat(srctang(1,:),nt,1); + dy = repmat(srctang(2,:),nt,1); + + ds = sqrt(dx.*dx+dy.*dy); + + taux = dx./ds; + tauy = dy./ds; + + rx = targ(1,:).' - src(1,:); + ry = targ(2,:).' - src(2,:); + r2 = rx.^2 + ry.^2; + dists = sqrt(r2); + inds = and((dists < R),(dists > 0)); + + + submat = 1/(2*zk^2).*(forth(:, :, 1).*(nx.*nx.*nx.*nxtarg) + forth(:, :, 2).*(nx.*nx.*nx.*nytarg + 3*nx.*nx.*ny.*nxtarg) + ... + forth(:, :, 4).*(3*nx.*nx.*ny.*nytarg + 3*nx.*ny.*ny.*nxtarg) + forth(:, :, 6).*(3*nx.*ny.*ny.*nytarg +ny.*ny.*ny.*nxtarg)+... + forth(:, :, 8).*(ny.*ny.*ny.*nytarg)) - 1/(2*zk^2).*(forthK(:, :, 1).*(nx.*nx.*nx.*nxtarg) + forthK(:, :, 2).*(nx.*nx.*nx.*nytarg + 3*nx.*nx.*ny.*nxtarg) + ... + forthK(:, :, 4).*(3*nx.*nx.*ny.*nytarg + 3*nx.*ny.*ny.*nxtarg) + forthK(:, :, 6).*(3*nx.*ny.*ny.*nytarg +ny.*ny.*ny.*nxtarg)+... + forthK(:, :, 8).*(ny.*ny.*ny.*nytarg)) + ... + (3/(2*zk^2).*(forth(:, :, 1).*(nx.*taux.*taux.*nxtarg)+ forth(:, :, 2).*(nx.*taux.*taux.*nytarg + 2*nx.*taux.*tauy.*nxtarg + ny.*taux.*taux.*nxtarg) +... + forth(:, :, 4).*(2*nx.*taux.*tauy.*nytarg + ny.*taux.*taux.*nytarg + nx.*tauy.*tauy.*nxtarg + 2*ny.*taux.*tauy.*nxtarg) + ... + forth(:, :, 6).*(nx.*tauy.*tauy.*nytarg +2*ny.*taux.*tauy.*nytarg + ny.*tauy.*tauy.*nxtarg) + forth(:, :, 8).*(ny.*tauy.*tauy.*nytarg))-... + 3/(2*zk^2).*(forthK(:, :, 1).*(nx.*taux.*taux.*nxtarg)+ forthK(:, :, 2).*(nx.*taux.*taux.*nytarg + 2*nx.*taux.*tauy.*nxtarg + ny.*taux.*taux.*nxtarg) +... + forthK(:, :, 4).*(2*nx.*taux.*tauy.*nytarg + ny.*taux.*taux.*nytarg + nx.*tauy.*tauy.*nxtarg + 2*ny.*taux.*tauy.*nxtarg) + ... + forthK(:, :, 6).*(nx.*tauy.*tauy.*nytarg +2*ny.*taux.*tauy.*nytarg + ny.*tauy.*tauy.*nxtarg) + forthK(:, :, 8).*(ny.*tauy.*tauy.*nytarg))); + % + if any(inds(:)) % if statement for evalauation points that are closed to the boundary + temp = 1i*imag(submat(inds)); + + nx = nx(inds); ny = ny(inds); + + nxtarg = nxtarg(inds); + nytarg = nytarg(inds); + + dx = dx(inds); + dy = dy(inds); + + + ds = sqrt(dx.*dx+dy.*dy); + + + taux = dx./ds; % normalization + tauy = dy./ds; + + + rx = rx(inds); + ry = ry(inds); + r2 = r2(inds); + + + normaldot = nx.*nxtarg + ny.*nytarg; + + + rn = rx.*nx + ry.*ny; + + rntarg = rx.*nxtarg + ry.*nytarg; + + rtau = rx.*taux + ry.*tauy; + + + taudotnsrc = taux.*nx + tauy.*ny; + taudotntarg = taux.*nxtarg + tauy.*nytarg; + + [smoothfirst, smoothsecond] = chnk.flex2d.smooth_part_derivatives(zk, rx, ry); + + eulerconstantpart = (smoothfirst(:, :, 1)).*(nx.*nx.*nx.*nxtarg) +... + (smoothfirst(:, :, 2)).*(nx.*nx.*nx.*nytarg + 3*nx.*nx.*ny.*nxtarg) + ... + (smoothfirst(:, :, 3)).*(3*nx.*nx.*ny.*nytarg + 3*nx.*ny.*ny.*nxtarg) +... + (smoothfirst(:, :, 4)).*(3*nx.*ny.*ny.*nytarg +ny.*ny.*ny.*nxtarg)+... + (smoothfirst(:, :, 5)).*(ny.*ny.*ny.*nytarg) + ... + 3.*((smoothfirst(:, :, 1)).*(nx.*taux.*taux.*nxtarg) +... + (smoothfirst(:, :, 2)).*(nx.*taux.*taux.*nytarg + 2*nx.*taux.*tauy.*nxtarg + ny.*taux.*taux.*nxtarg) +... + (smoothfirst(:, :, 3)).*(2*nx.*taux.*tauy.*nytarg + ny.*taux.*taux.*nytarg + nx.*tauy.*tauy.*nxtarg + 2*ny.*taux.*tauy.*nxtarg) + ... + (smoothfirst(:, :, 4)).*(nx.*tauy.*tauy.*nytarg +2*ny.*taux.*tauy.*nytarg + ny.*tauy.*tauy.*nxtarg) +... + (smoothfirst(:, :, 5)).*(ny.*tauy.*tauy.*nytarg)); + + puresmoothpart = (smoothsecond(:, :, 1)).*(nx.*nx.*nx.*nxtarg) +... + (smoothsecond(:, :, 2)).*(nx.*nx.*nx.*nytarg + 3*nx.*nx.*ny.*nxtarg) + ... + (smoothsecond(:, :, 3)).*(3*nx.*nx.*ny.*nytarg + 3*nx.*ny.*ny.*nxtarg) +... + (smoothsecond(:, :, 4)).*(3*nx.*ny.*ny.*nytarg +ny.*ny.*ny.*nxtarg)+... + (smoothsecond(:, :, 5)).*(ny.*ny.*ny.*nytarg) + ... + 3.*((smoothsecond(:, :, 1)).*(nx.*taux.*taux.*nxtarg) +... + (smoothsecond(:, :, 2)).*(nx.*taux.*taux.*nytarg + 2*nx.*taux.*tauy.*nxtarg + ny.*taux.*taux.*nxtarg) +... + (smoothsecond(:, :, 3)).*(2*nx.*taux.*tauy.*nytarg + ny.*taux.*taux.*nytarg + nx.*tauy.*tauy.*nxtarg + 2*ny.*taux.*tauy.*nxtarg) + ... + (smoothsecond(:, :, 4)).*(nx.*tauy.*tauy.*nytarg +2*ny.*taux.*tauy.*nytarg + ny.*tauy.*tauy.*nxtarg) +... + (smoothsecond(:, :, 5)).*(ny.*tauy.*tauy.*nytarg)); + + submat(inds) = -(3/(4*pi)).*(normaldot)./r2 + (3/(2*pi)).*(rntarg.*rn)./(r2.^2) +... + (3/(2*pi)).*((rn.^2).*normaldot)./(r2.^2) - (2/pi).*((rn.^3).*rntarg)./(r2.^3) + ... + 3.*(-1/(2*pi).*(taudotntarg.*taudotnsrc)./r2 + (1/pi).*(rntarg.*rtau.*taudotnsrc)./(r2.^2) + ... + (1/(2*pi)).*(normaldot.*(rtau.^2))./(r2.^2) + (1/pi).*rtau.*rn.*taudotntarg./(r2.^2) -... + (2/pi).*(rntarg.*rn.*(rtau.^2))./(r2.^3) - (1/(4*pi)).*normaldot./(r2) + (1/(2*pi)).*(rntarg.*rn)./(r2.^2)) + ... + temp + eulerconstantpart + puresmoothpart; + + + end + + +end + +% free plate kernel for the modified biharmonic problem (K11 with no +% hilbert transform subtraction, K12 kernel, K22 with no curvature part) K21 is +% handled in a separate type. +if strcmpi(type, 'free plate first part') + srcnorm = srcinfo.n; + targnorm = targinfo.n; + targtang = targinfo.d; + coefs = varargin{1}; + + + + + [~, ~, hess, third, ~] = chnk.flex2d.helmdiffgreen(zk, src, targ); % Hankel part + nx = repmat(srcnorm(1,:),nt,1); + ny = repmat(srcnorm(2,:),nt,1); + + nxtarg = repmat((targnorm(1,:)).',1,ns); + nytarg = repmat((targnorm(2,:)).',1,ns); + + zkimag = (1i)*zk; + [~, ~, hessK, thirdK, ~] = chnk.flex2d.helmdiffgreen(zkimag, src, targ); % modified bessel K part + + + + dx1 = repmat((targtang(1,:)).',1,ns); + dy1 = repmat((targtang(2,:)).',1,ns); + + + + ds1 = sqrt(dx1.*dx1+dy1.*dy1); + + + + tauxtarg = dx1./ds1; + tauytarg = dy1./ds1; + + + + + K11 = -(1/(2*zk^2).*(third(:, :, 1).*(nxtarg.*nxtarg.*nx) + third(:, :, 2).*(nxtarg.*nxtarg.*ny + 2*nxtarg.*nytarg.*nx) +... + third(:, :, 3).*(2*nxtarg.*nytarg.*ny + nytarg.*nytarg.*nx) +... + third(:, :, 4).*(nytarg.*nytarg.*ny)) - ... + 1/(2*zk^2).*(thirdK(:, :, 1).*(nxtarg.*nxtarg.*nx) + thirdK(:, :, 2).*(nxtarg.*nxtarg.*ny + 2*nxtarg.*nytarg.*nx) +... + thirdK(:, :, 3).*(2*nxtarg.*nytarg.*ny + nytarg.*nytarg.*nx) +... + thirdK(:, :, 4).*(nytarg.*nytarg.*ny))) - ... + coefs(1).*(1/(2*zk^2).*(third(:, :, 1).*(tauxtarg.*tauxtarg.*nx) + third(:, :, 2).*(tauxtarg.*tauxtarg.*ny + 2*tauxtarg.*tauytarg.*nx) +... + third(:, :, 3).*(2*tauxtarg.*tauytarg.*ny + tauytarg.*tauytarg.*nx) +... + third(:, :, 4).*(tauytarg.*tauytarg.*ny)) - ... + 1/(2*zk^2).*(thirdK(:, :, 1).*(tauxtarg.*tauxtarg.*nx) + thirdK(:, :, 2).*(tauxtarg.*tauxtarg.*ny + 2*tauxtarg.*tauytarg.*nx) +... + thirdK(:, :, 3).*(2*tauxtarg.*tauytarg.*ny + tauytarg.*tauytarg.*nx) +... + thirdK(:, :, 4).*(tauytarg.*tauytarg.*ny))); % first kernel with no hilbert transforms (G_{nx nx ny + nu G_{taux taux ny}). + + + + + K12 = 1/(2*zk^2).*(hess(:, :, 1).*nxtarg.*nxtarg + hess(:, :, 2).*(2*nxtarg.*nytarg) + hess(:, :, 3).*nytarg.*nytarg)-... + 1/(2*zk^2).*(hessK(:, :, 1).*nxtarg.*nxtarg + hessK(:, :, 2).*(2*nxtarg.*nytarg) + hessK(:, :, 3).*nytarg.*nytarg)+... + coefs(1)/(2*zk^2).*(hess(:, :, 1).*tauxtarg.*tauxtarg + hess(:, :, 2).*(2*tauxtarg.*tauytarg) + hess(:, :, 3).*tauytarg.*tauytarg)-... + coefs(1)/(2*zk^2).*(hessK(:, :, 1).*tauxtarg.*tauxtarg + hessK(:, :, 2).*(2*tauxtarg.*tauytarg) + ... + hessK(:, :, 3).*tauytarg.*tauytarg) ; % G_{nx nx} + nu G_{taux taux} + + K21 = 0; + + + + + K22 = 1./(2*zk^2).*(third(:, :, 1).*(nxtarg.*nxtarg.*nxtarg) + third(:, :, 2).*(3*nxtarg.*nxtarg.*nytarg) +... + third(:, :, 3).*(3*nxtarg.*nytarg.*nytarg) + third(:, :, 4).*(nytarg.*nytarg.*nytarg)) - ... + 1./(2*zk^2).*(thirdK(:, :, 1).*(nxtarg.*nxtarg.*nxtarg) + thirdK(:, :, 2).*(3*nxtarg.*nxtarg.*nytarg)+... + thirdK(:, :, 3).*(3*nxtarg.*nytarg.*nytarg) + thirdK(:, :, 4).*(nytarg.*nytarg.*nytarg)) +... + (2-coefs(1))/(2*zk^2).*(third(:, :, 1).*(tauxtarg.*tauxtarg.*nxtarg) + third(:, :, 2).*(tauxtarg.*tauxtarg.*nytarg + 2*tauxtarg.*tauytarg.*nxtarg) +... + third(:, :, 3).*(2*tauxtarg.*tauytarg.*nytarg + tauytarg.*tauytarg.*nxtarg) +... + + third(:, :, 4).*(tauytarg.*tauytarg.*nytarg)) - ... + (2-coefs(1))/(2*zk^2).*(thirdK(:, :, 1).*(tauxtarg.*tauxtarg.*nxtarg) + thirdK(:, :, 2).*(tauxtarg.*tauxtarg.*nytarg + 2*tauxtarg.*tauytarg.*nxtarg) +... + thirdK(:, :, 3).*(2*tauxtarg.*tauytarg.*nytarg + tauytarg.*tauytarg.*nxtarg) +... + + thirdK(:, :, 4).*(tauytarg.*tauytarg.*nytarg)); % G_{nx nx nx} + (2-nu) G_{taux taux nx} + + submat = zeros(2*nt,2*ns); + + submat(1:2:end,1:2:end) = K11; + submat(1:2:end,2:2:end) = K12; + + submat(2:2:end,1:2:end) = K21; + submat(2:2:end,2:2:end) = K22; + +end + + +% free plate kernel K21 for the modified biharmonic problem. This part +% handles the singularity subtraction and swap the evaluation to its +% asymptotic expansions if the targets and sources are close. + +if strcmpi(type, 'free plate K21 first part') + srcnorm = srcinfo.n; + srctang = srcinfo.d; + targnorm = targinfo.n; + targtang = targinfo.d; + coefs = varargin{1}; + %R = 0.001; % radius of kernel replacement + + + + [~, ~,~, ~, forth] = chnk.flex2d.helmdiffgreen(zk, src, targ); % Hankel part + nx = repmat(srcnorm(1,:),nt,1); + ny = repmat(srcnorm(2,:),nt,1); + + nxtarg = repmat((targnorm(1,:)).',1,ns); + nytarg = repmat((targnorm(2,:)).',1,ns); + + zkimag = (1i)*zk; + [~, ~, ~,~, forthK] = chnk.flex2d.helmdiffgreen(zkimag, src, targ); % modified bessel K part + + + + + + dx = repmat(srctang(1,:),nt,1); + dy = repmat(srctang(2,:),nt,1); + + dx1 = repmat((targtang(1,:)).',1,ns); + dy1 = repmat((targtang(2,:)).',1,ns); + + + ds = sqrt(dx.*dx+dy.*dy); + ds1 = sqrt(dx1.*dx1+dy1.*dy1); + + taux = dx./ds; % normalization + tauy = dy./ds; + + tauxtarg = dx1./ds1; + tauytarg = dy1./ds1; + + rx = targ(1,:).' - src(1,:); + ry = targ(2,:).' - src(2,:); + r2 = rx.^2 + ry.^2; + dists = sqrt(r2); + %inds = and((dists <= R),(dists > 0)); + + submat = -(1/(2*zk^2).*(forth(:, :, 1).*(nxtarg.*nxtarg.*nxtarg.*nx) + forth(:, :, 2).*(nxtarg.*nxtarg.*nxtarg.*ny + 3*nxtarg.*nxtarg.*nytarg.*nx) + ... + forth(:, :, 3).*(3*nxtarg.*nxtarg.*nytarg.*ny + 3*nxtarg.*nytarg.*nytarg.*nx) +... + forth(:, :, 4).*(3*nxtarg.*nytarg.*nytarg.*ny +nytarg.*nytarg.*nytarg.*nx)+... + forth(:, :, 5).*(nytarg.*nytarg.*nytarg.*ny)) -... + 1/(2*zk^2).*(forthK(:, :, 1).*(nxtarg.*nxtarg.*nxtarg.*nx) + forthK(:, :, 2).*(nxtarg.*nxtarg.*nxtarg.*ny + 3*nxtarg.*nxtarg.*nytarg.*nx) + ... + forthK(:, :, 3).*(3*nxtarg.*nxtarg.*nytarg.*ny + 3*nxtarg.*nytarg.*nytarg.*nx) +... + forthK(:, :, 4).*(3*nxtarg.*nytarg.*nytarg.*ny +nytarg.*nytarg.*nytarg.*nx)+... + forthK(:, :, 5).*(nytarg.*nytarg.*nytarg.*ny))) + (3/(4*pi)).*(nxtarg.*nx + nytarg.*ny)./r2 - ... + ((2-coefs(1))/(2*zk^2).*(forth(:, :, 1).*(tauxtarg.*tauxtarg.*nxtarg.*nx) + forth(:, :, 2).*(tauxtarg.*tauxtarg.*nxtarg.*ny + tauxtarg.*tauxtarg.*nytarg.*nx + 2*tauxtarg.*tauytarg.*nxtarg.*nx)+... + forth(:, :, 3).*(tauxtarg.*tauxtarg.*nytarg.*ny + 2*tauxtarg.*tauytarg.*nxtarg.*ny + tauytarg.*tauytarg.*nxtarg.*nx + 2*tauxtarg.*tauytarg.*nytarg.*nx)+... + forth(:, :, 4).*(tauytarg.*tauytarg.*nxtarg.*ny + 2*tauxtarg.*tauytarg.*nytarg.*ny + tauytarg.*tauytarg.*nytarg.*nx) +... + forth(:, :, 5).*(tauytarg.*tauytarg.*nytarg.*ny)) -... + (2-coefs(1))/(2*zk^2).*(forthK(:, :, 1).*(tauxtarg.*tauxtarg.*nxtarg.*nx) + forthK(:, :, 2).*(tauxtarg.*tauxtarg.*nxtarg.*ny + tauxtarg.*tauxtarg.*nytarg.*nx + 2*tauxtarg.*tauytarg.*nxtarg.*nx)+... + forthK(:, :, 3).*(tauxtarg.*tauxtarg.*nytarg.*ny + 2*tauxtarg.*tauytarg.*nxtarg.*ny + tauytarg.*tauytarg.*nxtarg.*nx + 2*tauxtarg.*tauytarg.*nytarg.*nx)+... + forthK(:, :, 4).*(tauytarg.*tauytarg.*nxtarg.*ny + 2*tauxtarg.*tauytarg.*nytarg.*ny + tauytarg.*tauytarg.*nytarg.*nx) +... + forthK(:, :, 5).*(tauytarg.*tauytarg.*nytarg.*ny))) +... + (2-coefs(1))/(4*pi).*(nxtarg.*nx + nytarg.*ny)./r2 - ... + ((2-coefs(1))/(2*pi)).*(nxtarg.*nx + nytarg.*ny).*(rx.*tauxtarg + ry.*tauytarg).*(rx.*tauxtarg + ry.*tauytarg)./(r2.^2) - ... + ((1+coefs(1))/(4*pi)).*((taux.*tauxtarg + tauy.*tauytarg)./(r2) - 2*(rx.*tauxtarg + ry.*tauytarg).*(rx.*taux + ry.*tauy)./(r2.^2))-... + (3/(4*pi)).*(nxtarg.*nx + nytarg.*ny)./r2 - (2-coefs(1))/(4*pi).*(nxtarg.*nx + nytarg.*ny)./r2 + ... + (2-coefs(1))/(2*pi).*(nxtarg.*nx + nytarg.*ny).*(rx.*tauxtarg + ry.*tauytarg).*(rx.*tauxtarg + ry.*tauytarg)./(r2.^2); + % third kernel with singularity subtraction and H[delta'] + % + % if any(inds(:)) % if statement for evalauation points that are closed to the boundary + % temp = 1i*imag(submat(inds)); + % + % nx = nx(inds); ny = ny(inds); + % + % nxtarg = nxtarg(inds); + % nytarg = nytarg(inds); + % + % dx = dx(inds); + % dy = dy(inds); + % + % dx1 = dx1(inds); + % dy1 = dy1(inds); + % + % ds = sqrt(dx.*dx+dy.*dy); + % ds1 = sqrt(dx1.*dx1+dy1.*dy1); + % + % taux = dx./ds; % normalization + % tauy = dy./ds; + % + % tauxtarg = dx1./ds1; + % tauytarg = dy1./ds1; + % + % rx = rx(inds); + % ry = ry(inds); + % r2 = r2(inds); + % + % tauxtargnsrc = tauxtarg.*nx + tauytarg.*ny; + % tauxtargntarg = tauxtarg.*nxtarg + tauytarg.*nytarg; + % + % normaldot = nx.*nxtarg + ny.*nytarg; + % tangentdot = taux.*tauxtarg + tauy.*tauytarg; + % + % rn = rx.*nx + ry.*ny; + % + % rntarg = rx.*nxtarg + ry.*nytarg; + % + % rtau = rx.*taux + ry.*tauy; + % + % rtautarg = rx.*tauxtarg + ry.*tauytarg; + % + % + % + % [smoothfirst, smoothsecond] = chnk.flex2d.smooth_part_derivatives(zk, rx, ry); + % + % eulerconstantpart = (smoothfirst(:, :, 1)).*(nxtarg.*nxtarg.*nxtarg.*nx) +... + % (smoothfirst(:, :, 2)).*(nxtarg.*nxtarg.*nxtarg.*ny + 3*nxtarg.*nxtarg.*nytarg.*nx) + ... + % (smoothfirst(:, :, 3)).*(3*nxtarg.*nxtarg.*nytarg.*ny + 3*nxtarg.*nytarg.*nytarg.*nx) +... + % (smoothfirst(:, :, 4)).*(3*nxtarg.*nytarg.*nytarg.*ny +nytarg.*nytarg.*nytarg.*nx)+... + % (smoothfirst(:, :, 5)).*(nytarg.*nytarg.*nytarg.*ny)+ ... + % (2-coefs(1)).*((smoothfirst(:, :, 1)).*(tauxtarg.*tauxtarg.*nxtarg.*nx) +... + % (smoothfirst(:, :, 2)).*(tauxtarg.*tauxtarg.*nxtarg.*ny + tauxtarg.*tauxtarg.*nytarg.*nx + 2*tauxtarg.*tauytarg.*nxtarg.*nx)+... + % (smoothfirst(:, :, 3)).*(tauxtarg.*tauxtarg.*nytarg.*ny + 2*tauxtarg.*tauytarg.*nxtarg.*ny + tauytarg.*tauytarg.*nxtarg.*nx + 2*tauxtarg.*tauytarg.*nytarg.*nx)+... + % (smoothfirst(:, :, 4)).*(tauytarg.*tauytarg.*nxtarg.*ny + 2*tauxtarg.*tauytarg.*nytarg.*ny + tauytarg.*tauytarg.*nytarg.*nx) +... + % (smoothfirst(:, :, 5)).*(tauytarg.*tauytarg.*nytarg.*ny)); + % + % puresmoothpart = (smoothsecond(:, :, 1)).*(nxtarg.*nxtarg.*nxtarg.*nx) +... + % (smoothsecond(:, :, 2)).*(nxtarg.*nxtarg.*nxtarg.*ny + 3*nxtarg.*nxtarg.*nytarg.*nx) + ... + % (smoothsecond(:, :, 3)).*(3*nxtarg.*nxtarg.*nytarg.*ny + 3*nxtarg.*nytarg.*nytarg.*nx) +... + % (smoothsecond(:, :, 4)).*(3*nxtarg.*nytarg.*nytarg.*ny +nytarg.*nytarg.*nytarg.*nx)+... + % (smoothsecond(:, :, 5)).*(nytarg.*nytarg.*nytarg.*ny)+ ... + % (2-coefs(1)).*((smoothsecond(:, :, 1)).*(tauxtarg.*tauxtarg.*nxtarg.*nx) +... + % (smoothsecond(:, :, 2)).*(tauxtarg.*tauxtarg.*nxtarg.*ny + tauxtarg.*tauxtarg.*nytarg.*nx + 2*tauxtarg.*tauytarg.*nxtarg.*nx)+... + % (smoothsecond(:, :, 3)).*(tauxtarg.*tauxtarg.*nytarg.*ny + 2*tauxtarg.*tauytarg.*nxtarg.*ny + tauytarg.*tauytarg.*nxtarg.*nx + 2*tauxtarg.*tauytarg.*nytarg.*nx)+... + % (smoothsecond(:, :, 4)).*(tauytarg.*tauytarg.*nxtarg.*ny + 2*tauxtarg.*tauytarg.*nytarg.*ny + tauytarg.*tauytarg.*nytarg.*nx) +... + % (smoothsecond(:, :, 5)).*(tauytarg.*tauytarg.*nytarg.*ny)); + % + % submat(inds) = 3/(2*pi).*(rn.*rntarg./(r2.^2)) + 3/(2*pi).*(rntarg.^2).*(normaldot)./(r2.^2) -... + % (2/pi).*(rn.*(rntarg.^3))./(r2.^3) + ... + % (2-coefs(1)).*(-1/(2*pi).*(tauxtargnsrc).*(tauxtargntarg)./(r2) +... + % 1/(pi).*(rn.*rtautarg.*tauxtargntarg)./(r2.^2) + ... + % 1/(pi).*(rtautarg.*tauxtargnsrc.*rntarg)./(r2.^2) -.... + % 2/(pi).*(rn.*rntarg.*(rtautarg.*rtautarg))./(r2.^3) + ... + % 1/(2*pi).*(rn.*rntarg)./(r2.^2)) - 3/(4*pi).*normaldot./r2 +... + % ((2-coefs(1))/(2*pi)).*((rtautarg.*rtautarg).*(normaldot))./(r2.^2) -... + % ((2-coefs(1))/(4*pi)).*(normaldot)./r2 - ... + % ((1+coefs(1))/(4*pi)).*((tangentdot)./(r2) - 2*(rtautarg).*(rtau)./(r2.^2)) + temp + eulerconstantpart + puresmoothpart; + % end +end + +% kernels in K11 with hilbert transform subtractions. +% (i.e. beta*(G_{nx nx tauy} + 1/4 H + nu*(G_{taux taux tauy} + 1/4 H)) + +if strcmpi(type, 'free plate hilbert subtract') + srcnorm = srcinfo.n; + srctang = srcinfo.d; + targnorm = targinfo.n; + targtang = targinfo.d; + coefs = varargin{1}; + + [~,~,~, third, ~] = chnk.flex2d.helmdiffgreen(zk, src, targ); % Hankel part + nx = repmat(srcnorm(1,:),nt,1); + ny = repmat(srcnorm(2,:),nt,1); + + nxtarg = repmat((targnorm(1,:)).',1,ns); + nytarg = repmat((targnorm(2,:)).',1,ns); + + zkimag = (1i)*zk; + [~,~, ~, thirdK, ~] = chnk.flex2d.helmdiffgreen(zkimag, src, targ); % modified bessel K part + [~,grad] = chnk.lap2d.green(src,targ,true); + + + + + dx = repmat(srctang(1,:),nt,1); + dy = repmat(srctang(2,:),nt,1); + + dx1 = repmat((targtang(1,:)).',1,ns); + dy1 = repmat((targtang(2,:)).',1,ns); + + + ds = sqrt(dx.*dx+dy.*dy); + ds1 = sqrt(dx1.*dx1+dy1.*dy1); + + taux = dx./ds; % normalization + tauy = dy./ds; + + tauxtarg = dx1./ds1; + tauytarg = dy1./ds1; + + + hilb = 2*(grad(:,:,1).*ny - grad(:,:,2).*nx); + + submat = -((1+ coefs(1))/2)*(1./(2*zk^2).*(third(:, :, 1).*(nxtarg.*nxtarg.*taux) + third(:, :, 2).*(nxtarg.*nxtarg.*tauy+ 2*nxtarg.*nytarg.*taux) +... + third(:, :, 3).*(2*nxtarg.*nytarg.*tauy +nytarg.*nytarg.*taux) +... + third(:, :, 4).*(nytarg.*nytarg.*tauy)) - ... + 1./(2*zk^2).*(thirdK(:, :, 1).*(nxtarg.*nxtarg.*taux) + thirdK(:, :, 2).*(nxtarg.*nxtarg.*tauy+ 2*nxtarg.*nytarg.*taux) +... + thirdK(:, :, 3).*(2*nxtarg.*nytarg.*tauy +nytarg.*nytarg.*taux) +... + thirdK(:, :, 4).*(nytarg.*nytarg.*tauy)))+ ((1+ coefs(1))/2).*0.25*hilb -... + ((1+ coefs(1))/2)*coefs(1).*(1./(2*zk^2).*(third(:, :, 1).*(tauxtarg.*tauxtarg.*taux) + third(:, :, 2).*(tauxtarg.*tauxtarg.*tauy + 2*tauxtarg.*tauytarg.*taux) +... + third(:, :, 3).*(2*tauxtarg.*tauytarg.*tauy + tauytarg.*tauytarg.*taux) +... + third(:, :, 4).*(tauytarg.*tauytarg.*tauy)) - ... + 1./(2*zk^2).*(thirdK(:, :, 1).*(tauxtarg.*tauxtarg.*taux) + thirdK(:, :, 2).*(tauxtarg.*tauxtarg.*tauy + 2*tauxtarg.*tauytarg.*taux) +... + thirdK(:, :, 3).*(2*tauxtarg.*tauytarg.*tauy + tauytarg.*tauytarg.*taux) +... + thirdK(:, :, 4).*(tauytarg.*tauytarg.*tauy))) + ((1+ coefs(1))/2)*coefs(1).*0.25*hilb; % hilbert subtraction + +end + +% kernels in K21 that are coupled with Hilbert transforms. +% (i.e beta* G_{nx nx nx tauy} + (2-nu)*beta*G_{taux taux tauy nx}) +if strcmpi(type, 'free plate coupled hilbert') + srctang = srcinfo.d; + targnorm = targinfo.n; + targtang = targinfo.d; + coefs = varargin{1}; + + [~, ~,~, ~, forth] = chnk.flex2d.helmdiffgreen(zk, src, targ); % Hankel part + + + nxtarg = repmat((targnorm(1,:)).',1,ns); + nytarg = repmat((targnorm(2,:)).',1,ns); + + zkimag = (1i)*zk; + [~, ~,~, ~, forthK] = chnk.flex2d.helmdiffgreen(zkimag, src, targ); % modified bessel K part + + + + + + dx = repmat(srctang(1,:),nt,1); + dy = repmat(srctang(2,:),nt,1); + + dx1 = repmat((targtang(1,:)).',1,ns); + dy1 = repmat((targtang(2,:)).',1,ns); + + + ds = sqrt(dx.*dx+dy.*dy); + ds1 = sqrt(dx1.*dx1+dy1.*dy1); + + taux = dx./ds; % normalization + tauy = dy./ds; + + tauxtarg = dx1./ds1; + tauytarg = dy1./ds1; + + + + submat = -((1+ coefs(1))/2).*(1/(2*zk^2).*(forth(:, :, 1).*(nxtarg.*nxtarg.*nxtarg.*taux) + ... + forth(:, :, 2).*(nxtarg.*nxtarg.*nxtarg.*tauy + 3*nxtarg.*nxtarg.*nytarg.*taux) + ... + forth(:, :, 3).*(3*nxtarg.*nxtarg.*nytarg.*tauy + 3*nxtarg.*nytarg.*nytarg.*taux) +... + forth(:, :, 4).*(3*nxtarg.*nytarg.*nytarg.*tauy + nytarg.*nytarg.*nytarg.*taux) +... + forth(:, :, 5).*(nytarg.*nytarg.*nytarg.*tauy)) - ... + 1/(2*zk^2).*(forthK(:, :, 1).*(nxtarg.*nxtarg.*nxtarg.*taux) + ... + forthK(:, :, 2).*(nxtarg.*nxtarg.*nxtarg.*tauy + 3*nxtarg.*nxtarg.*nytarg.*taux) + ... + forthK(:, :, 3).*(3*nxtarg.*nxtarg.*nytarg.*tauy + 3*nxtarg.*nytarg.*nytarg.*taux) +... + forthK(:, :, 4).*(3*nxtarg.*nytarg.*nytarg.*tauy + nytarg.*nytarg.*nytarg.*taux) +... + forthK(:, :, 5).*(nytarg.*nytarg.*nytarg.*tauy))) -... + ((2-coefs(1))/2)*(1+coefs(1)).*(1/(2*zk^2).*(forth(:, :, 1).*(tauxtarg.*tauxtarg.*nxtarg.*taux) + ... + forth(:, :, 2).*(tauxtarg.*tauxtarg.*nxtarg.*tauy + tauxtarg.*tauxtarg.*nytarg.*taux + 2*tauxtarg.*tauytarg.*nxtarg.*taux) + ... + forth(:, :, 3).*(tauxtarg.*tauxtarg.*nytarg.*tauy + 2*tauxtarg.*tauytarg.*nxtarg.*tauy + tauytarg.*tauytarg.*nxtarg.*taux + 2*tauxtarg.*tauytarg.*nytarg.*taux) +... + forth(:, :, 4).*(tauytarg.*tauytarg.*nxtarg.*tauy + 2*tauxtarg.*tauytarg.*nytarg.*tauy + tauytarg.*tauytarg.*nytarg.*taux) +... + forth(:, :, 5).*(tauytarg.*tauytarg.*nytarg.*tauy)) - ... + 1/(2*zk^2).*(forthK(:, :, 1).*(tauxtarg.*tauxtarg.*nxtarg.*taux) + ... + forthK(:, :, 2).*(tauxtarg.*tauxtarg.*nxtarg.*tauy + tauxtarg.*tauxtarg.*nytarg.*taux + 2*tauxtarg.*tauytarg.*nxtarg.*taux) + ... + forthK(:, :, 3).*(tauxtarg.*tauxtarg.*nytarg.*tauy + 2*tauxtarg.*tauytarg.*nxtarg.*tauy + tauytarg.*tauytarg.*nxtarg.*taux + 2*tauxtarg.*tauytarg.*nytarg.*taux) +... + forthK(:, :, 4).*(tauytarg.*tauytarg.*nxtarg.*tauy + 2*tauxtarg.*tauytarg.*nytarg.*tauy + tauytarg.*tauytarg.*nytarg.*taux) +... + forthK(:, :, 5).*(tauytarg.*tauytarg.*nytarg.*tauy))); + +end + +% Updated part in K21 that is not coupled with Hilbert transform. +%(i.e. (1-nu)*(-G_{nx nx ny} + G_{taux taux ny}). + +if strcmpi(type, 'free plate K21 second part') + srcnorm = srcinfo.n; + targnorm = targinfo.n; + targtang = targinfo.d; + coefs = varargin{1}; + + + nx = repmat(srcnorm(1,:),nt,1); + ny = repmat(srcnorm(2,:),nt,1); + + nxtarg = repmat((targnorm(1,:)).',1,ns); + nytarg = repmat((targnorm(2,:)).',1,ns); + + + dx1 = repmat((targtang(1,:)).',1,ns); + dy1 = repmat((targtang(2,:)).',1,ns); + + + + ds1 = sqrt(dx1.*dx1+dy1.*dy1); + + + tauxtarg = dx1./ds1; + tauytarg = dy1./ds1; + + zkimag = (1i)*zk; + + [~,~,~, third, ~] = chnk.flex2d.helmdiffgreen(zk, src, targ); % Hankel part + [~,~,~, thirdK, ~] = chnk.flex2d.helmdiffgreen(zkimag, src, targ); % modified bessel K part + + submat = ((1-coefs(1))/(2*zk^2).*(third(:, :, 1).*(nxtarg.*nxtarg.*nx) + third(:, :, 2).*(nxtarg.*nxtarg.*ny + 2*nxtarg.*nytarg.*nx) +... + third(:, :, 3).*(2*nxtarg.*nytarg.*ny + nytarg.*nytarg.*nx) +... + third(:, :, 4).*(nytarg.*nytarg.*ny)) - ... + (1-coefs(1))/(2*zk^2).*(thirdK(:, :, 1).*(nxtarg.*nxtarg.*nx) + thirdK(:, :, 2).*(nxtarg.*nxtarg.*ny + 2*nxtarg.*nytarg.*nx) +... + thirdK(:, :, 3).*(2*nxtarg.*nytarg.*ny + nytarg.*nytarg.*nx) +... + thirdK(:, :, 4).*(nytarg.*nytarg.*ny))) - ... + ((1-coefs(1))/(2*zk^2).*(third(:, :, 1).*(tauxtarg.*tauxtarg.*nx) + third(:, :, 2).*(tauxtarg.*tauxtarg.*ny + 2*tauxtarg.*tauytarg.*nx) +... + third(:, :, 3).*(2*tauxtarg.*tauytarg.*ny + tauytarg.*tauytarg.*nx) +... + third(:, :, 4).*(tauytarg.*tauytarg.*ny)) - ... + (1-coefs(1))/(2*zk^2).*(thirdK(:, :, 1).*(tauxtarg.*tauxtarg.*nx) + thirdK(:, :, 2).*(tauxtarg.*tauxtarg.*ny + 2*tauxtarg.*tauytarg.*nx) +... + thirdK(:, :, 3).*(2*tauxtarg.*tauytarg.*ny + tauytarg.*tauytarg.*nx) +... + thirdK(:, :, 4).*(tauytarg.*tauytarg.*ny))); % (1-nu)*(- G_{nx nx ny} + G_{ny taux taux}) +end + +% Updated part in K21 that is coupled with Hilbert transform. +%(i.e. (1-nu)*(beta*(G_{taux taux tauy} + 1/4 H) - beta*(G_{nx nx tauy} + 1/4 H)). + +if strcmpi(type, 'free plate K21 hilbert part') + srcnorm = srcinfo.n; + srctang = srcinfo.d; + targnorm = targinfo.n; + targtang = targinfo.d; + coefs = varargin{1}; + + [~, ~, ~,third, ~] = chnk.flex2d.helmdiffgreen(zk, src, targ); % Hankel part + nx = repmat(srcnorm(1,:),nt,1); + ny = repmat(srcnorm(2,:),nt,1); + + nxtarg = repmat((targnorm(1,:)).',1,ns); + nytarg = repmat((targnorm(2,:)).',1,ns); + + zkimag = (1i)*zk; + [~, ~, ~, thirdK, ~] = chnk.flex2d.helmdiffgreen(zkimag, src, targ); % modified bessel K part + [~,grad] = chnk.lap2d.green(src,targ,true); + + + + + dx = repmat(srctang(1,:),nt,1); + dy = repmat(srctang(2,:),nt,1); + + dx1 = repmat((targtang(1,:)).',1,ns); + dy1 = repmat((targtang(2,:)).',1,ns); + + + ds = sqrt(dx.*dx+dy.*dy); + ds1 = sqrt(dx1.*dx1+dy1.*dy1); + + taux = dx./ds; % normalization + tauy = dy./ds; + + tauxtarg = dx1./ds1; + tauytarg = dy1./ds1; + + + hilb = 2*(grad(:,:,1).*ny - grad(:,:,2).*nx); + + submat = (1-coefs(1)).*(-((1+ coefs(1))/2).*(1./(2*zk^2).*(third(:, :, 1).*(tauxtarg.*tauxtarg.*taux) + third(:, :, 2).*(tauxtarg.*tauxtarg.*tauy + 2*tauxtarg.*tauytarg.*taux) +... + third(:, :, 3).*(2*tauxtarg.*tauytarg.*tauy + tauytarg.*tauytarg.*taux) +... + third(:, :, 4).*(tauytarg.*tauytarg.*tauy)) - ... + 1./(2*zk^2).*(thirdK(:, :, 1).*(tauxtarg.*tauxtarg.*taux) + thirdK(:, :, 2).*(tauxtarg.*tauxtarg.*tauy + 2*tauxtarg.*tauytarg.*taux) +... + thirdK(:, :, 3).*(2*tauxtarg.*tauytarg.*tauy + tauytarg.*tauytarg.*taux) +... + thirdK(:, :, 4).*(tauytarg.*tauytarg.*tauy))) + ((1+ coefs(1))/2).*0.25*hilb + ... + ((1+ coefs(1))/2)*(1./(2*zk^2).*(third(:, :, 1).*(nxtarg.*nxtarg.*taux) + third(:, :, 2).*(nxtarg.*nxtarg.*tauy+ 2*nxtarg.*nytarg.*taux) +... + third(:, :, 3).*(2*nxtarg.*nytarg.*tauy +nytarg.*nytarg.*taux) +... + third(:, :, 4).*(nytarg.*nytarg.*tauy)) - ... + 1./(2*zk^2).*(thirdK(:, :, 1).*(nxtarg.*nxtarg.*taux) + thirdK(:, :, 2).*(nxtarg.*nxtarg.*tauy+ 2*nxtarg.*nytarg.*taux) +... + thirdK(:, :, 3).*(2*nxtarg.*nytarg.*tauy +nytarg.*nytarg.*taux) +... + thirdK(:, :, 4).*(nytarg.*nytarg.*tauy))) - ((1+ coefs(1))/2).*0.25*hilb) ; + + + +end + +% Updated part in K22. (i.e. (1-nu)*(G_{taux taux} - G_{nx nx}) + +if strcmpi(type, 'free plate K22 second part') + targnorm = targinfo.n; + targtang = targinfo.d; + coefs = varargin{1}; + + + + [~, ~, hess, ~, ~] = chnk.flex2d.helmdiffgreen(zk, src, targ); % Hankel part + + nxtarg = repmat((targnorm(1,:)).',1,ns); + nytarg = repmat((targnorm(2,:)).',1,ns); + + zkimag = (1i)*zk; + [~, ~, hessK, ~, ~] = chnk.flex2d.helmdiffgreen(zkimag, src, targ); % modified bessel K part + + + + + + dx1 = repmat((targtang(1,:)).',1,ns); + dy1 = repmat((targtang(2,:)).',1,ns); + + + + ds1 = sqrt(dx1.*dx1+dy1.*dy1); + + + tauxtarg = dx1./ds1; + tauytarg = dy1./ds1; + + + submat = (1-coefs(1)).*(1/(2*zk^2).*(hess(:, :, 1).*tauxtarg.*tauxtarg + hess(:, :, 2).*(2*tauxtarg.*tauytarg) + hess(:, :, 3).*tauytarg.*tauytarg)-... + 1/(2*zk^2).*(hessK(:, :, 1).*tauxtarg.*tauxtarg + hessK(:, :, 2).*(2*tauxtarg.*tauytarg) + ... + hessK(:, :, 3).*tauytarg.*tauytarg)-... + (1/(2*zk^2).*(hess(:, :, 1).*nxtarg.*nxtarg + hess(:, :, 2).*(2*nxtarg.*nytarg) + hess(:, :, 3).*nytarg.*nytarg)-... + 1/(2*zk^2).*(hessK(:, :, 1).*nxtarg.*nxtarg + hessK(:, :, 2).*(2*nxtarg.*nytarg) + hessK(:, :, 3).*nytarg.*nytarg))); + + +end + +% K11 kernel of the supported plate with 'Hilbert transform subtractions'. +if strcmpi(type, 'supported plate K11') + srcnorm = srcinfo.n; + srctang = srcinfo.d; + targnorm = targinfo.n; + targtang = targinfo.d; + coefs = varargin{1}; + + + + + [~, third, ~] = chnk.flex2d.thirdforth_derivatives(zk, src, targ); % Hankel part + nx = repmat(srcnorm(1,:),nt,1); + ny = repmat(srcnorm(2,:),nt,1); + + nxtarg = repmat((targnorm(1,:)).',1,ns); + nytarg = repmat((targnorm(2,:)).',1,ns); + + zkimag = (1i)*zk; + [~, thirdK, ~] = chnk.flex2d.thirdforth_derivatives(zkimag, src, targ); % modified bessel K part + [~,grad] = chnk.lap2d.green(src,targ,true); % lap kernels for the adjoint of Hilbert transform + + dx = repmat(srctang(1,:),nt,1); + dy = repmat(srctang(2,:),nt,1); + + dx1 = repmat((targtang(1,:)).',1,ns); + dy1 = repmat((targtang(2,:)).',1,ns); + + + ds = sqrt(dx.*dx+dy.*dy); + ds1 = sqrt(dx1.*dx1+dy1.*dy1); + + taux = dx./ds; % normalization + tauy = dy./ds; + + tauxtarg = dx1./ds1; + tauytarg = dy1./ds1; + + + stau = 2*(grad(:,:,1).*nytarg - grad(:,:,2).*nxtarg); % kernel of adjoint of Hilbert transform + + c = (3-coefs(1))/(1 + coefs(1)); + + submat = -(1/(2*zk^2).*(third(:, :, 1).*(nx.*nx.*tauxtarg) + third(:, :, 2).*(nx.*nx.*tauytarg + 2*nx.*ny.*tauxtarg) +... + third(:, :, 4).*(2*nx.*ny.*tauytarg + ny.*ny.*tauxtarg) +... + third(:, :, 6).*(ny.*ny.*tauytarg)) - ... + 1/(2*zk^2).*(thirdK(:, :, 1).*(nx.*nx.*tauxtarg) + thirdK(:, :, 2).*(nx.*nx.*tauytarg + 2*nx.*ny.*tauxtarg) +... + thirdK(:, :, 4).*(2*nx.*ny.*tauytarg + ny.*ny.*tauxtarg) +... + thirdK(:, :, 6).*(ny.*ny.*tauytarg))) - (0.25*stau) + ... + c.*(-(1/(2*zk^2).*(third(:, :, 1).*(taux.*taux.*tauxtarg) +... + third(:, :, 2).*(taux.*taux.*tauytarg + 2*taux.*tauy.*tauxtarg) +... + third(:, :, 4).*(2*taux.*tauy.*tauytarg + tauy.*tauy.*tauxtarg) +... + third(:, :, 6).*(tauy.*tauy.*tauytarg)) - ... + 1/(2*zk^2).*(thirdK(:, :, 1).*(taux.*taux.*tauxtarg) +... + thirdK(:, :, 2).*(taux.*taux.*tauytarg + 2*taux.*tauy.*tauxtarg) +... + thirdK(:, :, 4).*(2*taux.*tauy.*tauytarg + tauy.*tauy.*tauxtarg) +... + thirdK(:, :, 6).*(tauy.*tauy.*tauytarg))) - (0.25)*stau); + + % (G_{ny ny nx} - 0.25*Stau) + c*(G_{tauy tauy taux}-0.25 Stau)). + + + +end + + +if strcmpi(type, 'test kernel') % test G_{taux taux tauy nx} + srcnorm = srcinfo.n; + srctang = srcinfo.d; + targnorm = targinfo.n; + targtang = targinfo.d; + coefs = varargin{1}; + %R = 0.001; % radius of kernel replacement + + + + [~, ~, forth] = chnk.flex2d.thirdforth_derivatives(zk, src, targ); % Hankel part + nx = repmat(srcnorm(1,:),nt,1); + ny = repmat(srcnorm(2,:),nt,1); + + nxtarg = repmat((targnorm(1,:)).',1,ns); + nytarg = repmat((targnorm(2,:)).',1,ns); + + zkimag = (1i)*zk; + [~, ~, forthK] = chnk.flex2d.thirdforth_derivatives(zkimag, src, targ); % modified bessel K part + + + + + + dx = repmat(srctang(1,:),nt,1); + dy = repmat(srctang(2,:),nt,1); + + dx1 = repmat((targtang(1,:)).',1,ns); + dy1 = repmat((targtang(2,:)).',1,ns); + + + ds = sqrt(dx.*dx+dy.*dy); + ds1 = sqrt(dx1.*dx1+dy1.*dy1); + + taux = dx./ds; % normalization + tauy = dy./ds; + + tauxtarg = dx1./ds1; + tauytarg = dy1./ds1; + + rx = targ(1,:).' - src(1,:); + ry = targ(2,:).' - src(2,:); + r2 = rx.^2 + ry.^2; + dists = sqrt(r2); + %inds = and((dists <= R),(dists > 0)); + + submat = 1/(2*zk^2).*(forth(:, :, 1).*(nxtarg.*nxtarg.*nxtarg.*nx) + forth(:, :, 2).*(nxtarg.*nxtarg.*nxtarg.*ny + 3*nxtarg.*nxtarg.*nytarg.*nx) + ... + forth(:, :, 4).*(3*nxtarg.*nxtarg.*nytarg.*ny + 3*nxtarg.*nytarg.*nytarg.*nx) + forth(:, :, 6).*(3*nxtarg.*nytarg.*nytarg.*ny +nytarg.*nytarg.*nytarg.*nx)+... + forth(:, :, 8).*(nytarg.*nytarg.*nytarg.*ny)) -... + 1/(2*zk^2).*(forthK(:, :, 1).*(nxtarg.*nxtarg.*nxtarg.*nx) + forthK(:, :, 2).*(nxtarg.*nxtarg.*nxtarg.*ny + 3*nxtarg.*nxtarg.*nytarg.*nx) + ... + forthK(:, :, 4).*(3*nxtarg.*nxtarg.*nytarg.*ny + 3*nxtarg.*nytarg.*nytarg.*nx) + forthK(:, :, 6).*(3*nxtarg.*nytarg.*nytarg.*ny +nytarg.*nytarg.*nytarg.*nx)+... + forthK(:, :, 8).*(nytarg.*nytarg.*nytarg.*ny)) + (3/(4*pi)).*(nxtarg.*nx + nytarg.*ny)./r2 + ... + ((2-coefs(1))/(2*zk^2).*(forth(:, :, 1).*(tauxtarg.*tauxtarg.*nxtarg.*nx) + forth(:, :, 2).*(tauxtarg.*tauxtarg.*nxtarg.*ny + tauxtarg.*tauxtarg.*nytarg.*nx + 2*tauxtarg.*tauytarg.*nxtarg.*nx)+... + forth(:, :, 4).*(tauxtarg.*tauxtarg.*nytarg.*ny + 2*tauxtarg.*tauytarg.*nxtarg.*ny + tauytarg.*tauytarg.*nxtarg.*nx + 2*tauxtarg.*tauytarg.*nytarg.*nx)+... + forth(:, :, 6).*(tauytarg.*tauytarg.*nxtarg.*ny + 2*tauxtarg.*tauytarg.*nytarg.*ny + tauytarg.*tauytarg.*nytarg.*nx) +... + forth(:, :, 8).*(tauytarg.*tauytarg.*nytarg.*ny)) -... + (2-coefs(1))/(2*zk^2).*(forthK(:, :, 1).*(tauxtarg.*tauxtarg.*nxtarg.*nx) + forthK(:, :, 2).*(tauxtarg.*tauxtarg.*nxtarg.*ny + tauxtarg.*tauxtarg.*nytarg.*nx + 2*tauxtarg.*tauytarg.*nxtarg.*nx)+... + forthK(:, :, 4).*(tauxtarg.*tauxtarg.*nytarg.*ny + 2*tauxtarg.*tauytarg.*nxtarg.*ny + tauytarg.*tauytarg.*nxtarg.*nx + 2*tauxtarg.*tauytarg.*nytarg.*nx)+... + forthK(:, :, 6).*(tauytarg.*tauytarg.*nxtarg.*ny + 2*tauxtarg.*tauytarg.*nytarg.*ny + tauytarg.*tauytarg.*nytarg.*nx) +... + forthK(:, :, 8).*(tauytarg.*tauytarg.*nytarg.*ny)) +... + (2-coefs(1))/(4*pi).*(nxtarg.*nx + nytarg.*ny)./r2 - ... + ((2-coefs(1))/(2*pi)).*(nxtarg.*nx + nytarg.*ny).*(rx.*tauxtarg + ry.*tauytarg).*(rx.*tauxtarg + ry.*tauytarg)./(r2.^2)) - ... + ((1+coefs(1))/(4*pi)).*((taux.*tauxtarg + tauy.*tauytarg)./(r2) - 2*(rx.*tauxtarg + ry.*tauytarg).*(rx.*taux + ry.*tauy)./(r2.^2))-... + (3/(4*pi)).*(nxtarg.*nx + nytarg.*ny)./r2 - (2-coefs(1))/(4*pi).*(nxtarg.*nx + nytarg.*ny)./r2 + ... + (2-coefs(1))/(2*pi).*(nxtarg.*nx + nytarg.*ny).*(rx.*tauxtarg + ry.*tauytarg).*(rx.*tauxtarg + ry.*tauytarg)./(r2.^2); + % third kernel with singularity subtraction and H[delta'] + + % if any(inds(:)) % if statement for evalauation points that are closed to the boundary + % temp = 1i*imag(submat(inds)); + % + % nx = nx(inds); ny = ny(inds); + % + % nxtarg = nxtarg(inds); + % nytarg = nytarg(inds); + % + % dx = dx(inds); + % dy = dy(inds); + % + % dx1 = dx1(inds); + % dy1 = dy1(inds); + % + % ds = sqrt(dx.*dx+dy.*dy); + % ds1 = sqrt(dx1.*dx1+dy1.*dy1); + % + % taux = dx./ds; % normalization + % tauy = dy./ds; + % + % tauxtarg = dx1./ds1; + % tauytarg = dy1./ds1; + % + % rx = rx(inds); + % ry = ry(inds); + % r2 = r2(inds); + % + % tauxtargnsrc = tauxtarg.*nx + tauytarg.*ny; + % tauxtargntarg = tauxtarg.*nxtarg + tauytarg.*nytarg; + % + % normaldot = nx.*nxtarg + ny.*nytarg; + % tangentdot = taux.*tauxtarg + tauy.*tauytarg; + % + % rn = rx.*nx + ry.*ny; + % + % rntarg = rx.*nxtarg + ry.*nytarg; + % + % rtau = rx.*taux + ry.*tauy; + % + % rtautarg = rx.*tauxtarg + ry.*tauytarg; + % + % + % + % [smoothfirst, smoothsecond] = chnk.flex2d.smooth_part_derivatives(zk, rx, ry); + % + % eulerconstantpart = (smoothfirst(:, :, 1)).*(nxtarg.*nxtarg.*nxtarg.*nx) +... + % (smoothfirst(:, :, 2)).*(nxtarg.*nxtarg.*nxtarg.*ny + 3*nxtarg.*nxtarg.*nytarg.*nx) + ... + % (smoothfirst(:, :, 3)).*(3*nxtarg.*nxtarg.*nytarg.*ny + 3*nxtarg.*nytarg.*nytarg.*nx) +... + % (smoothfirst(:, :, 4)).*(3*nxtarg.*nytarg.*nytarg.*ny +nytarg.*nytarg.*nytarg.*nx)+... + % (smoothfirst(:, :, 5)).*(nytarg.*nytarg.*nytarg.*ny)+ ... + % (2-coefs(1)).*((smoothfirst(:, :, 1)).*(tauxtarg.*tauxtarg.*nxtarg.*nx) +... + % (smoothfirst(:, :, 2)).*(tauxtarg.*tauxtarg.*nxtarg.*ny + tauxtarg.*tauxtarg.*nytarg.*nx + 2*tauxtarg.*tauytarg.*nxtarg.*nx)+... + % (smoothfirst(:, :, 3)).*(tauxtarg.*tauxtarg.*nytarg.*ny + 2*tauxtarg.*tauytarg.*nxtarg.*ny + tauytarg.*tauytarg.*nxtarg.*nx + 2*tauxtarg.*tauytarg.*nytarg.*nx)+... + % (smoothfirst(:, :, 4)).*(tauytarg.*tauytarg.*nxtarg.*ny + 2*tauxtarg.*tauytarg.*nytarg.*ny + tauytarg.*tauytarg.*nytarg.*nx) +... + % (smoothfirst(:, :, 5)).*(tauytarg.*tauytarg.*nytarg.*ny)); + % + % puresmoothpart = (smoothsecond(:, :, 1)).*(nxtarg.*nxtarg.*nxtarg.*nx) +... + % (smoothsecond(:, :, 2)).*(nxtarg.*nxtarg.*nxtarg.*ny + 3*nxtarg.*nxtarg.*nytarg.*nx) + ... + % (smoothsecond(:, :, 3)).*(3*nxtarg.*nxtarg.*nytarg.*ny + 3*nxtarg.*nytarg.*nytarg.*nx) +... + % (smoothsecond(:, :, 4)).*(3*nxtarg.*nytarg.*nytarg.*ny +nytarg.*nytarg.*nytarg.*nx)+... + % (smoothsecond(:, :, 5)).*(nytarg.*nytarg.*nytarg.*ny)+ ... + % (2-coefs(1)).*((smoothsecond(:, :, 1)).*(tauxtarg.*tauxtarg.*nxtarg.*nx) +... + % (smoothsecond(:, :, 2)).*(tauxtarg.*tauxtarg.*nxtarg.*ny + tauxtarg.*tauxtarg.*nytarg.*nx + 2*tauxtarg.*tauytarg.*nxtarg.*nx)+... + % (smoothsecond(:, :, 3)).*(tauxtarg.*tauxtarg.*nytarg.*ny + 2*tauxtarg.*tauytarg.*nxtarg.*ny + tauytarg.*tauytarg.*nxtarg.*nx + 2*tauxtarg.*tauytarg.*nytarg.*nx)+... + % (smoothsecond(:, :, 4)).*(tauytarg.*tauytarg.*nxtarg.*ny + 2*tauxtarg.*tauytarg.*nytarg.*ny + tauytarg.*tauytarg.*nytarg.*nx) +... + % (smoothsecond(:, :, 5)).*(tauytarg.*tauytarg.*nytarg.*ny)); + % + % submat(inds) = 3/(2*pi).*(rn.*rntarg./(r2.^2)) + 3/(2*pi).*(rntarg.^2).*(normaldot)./(r2.^2) -... + % (2/pi).*(rn.*(rntarg.^3))./(r2.^3) + ... + % (2-coefs(1)).*(-1/(2*pi).*(tauxtargnsrc).*(tauxtargntarg)./(r2) +... + % 1/(pi).*(rn.*rtautarg.*tauxtargntarg)./(r2.^2) + ... + % 1/(pi).*(rtautarg.*tauxtargnsrc.*rntarg)./(r2.^2) -.... + % 2/(pi).*(rn.*rntarg.*(rtautarg.*rtautarg))./(r2.^3) + ... + % 1/(2*pi).*(rn.*rntarg)./(r2.^2)) - 3/(4*pi).*normaldot./r2 +... + % ((2-coefs(1))/(2*pi)).*((rtautarg.*rtautarg).*(normaldot))./(r2.^2) -... + % ((2-coefs(1))/(4*pi)).*(normaldot)./r2 - ... + % ((1+coefs(1))/(4*pi)).*((tangentdot)./(r2) - 2*(rtautarg).*(rtau)./(r2.^2)) + temp + eulerconstantpart + puresmoothpart; + %end +end + + + + + +% first-kernel of evaluation for the clamped plate +if strcmpi(type, 'first kernel') + srcnorm = srcinfo.n; + srctang = srcinfo.d; + nx = repmat(srcnorm(1,:),nt,1); + ny = repmat(srcnorm(2,:),nt,1); + dx = repmat(srctang(1,:),nt,1); + dy = repmat(srctang(2,:),nt,1); + ds = sqrt(dx.*dx+dy.*dy); + + taux = dx./ds; + tauy = dy./ds; + + + + [~, third, ~] = chnk.flex2d.thirdforth_derivatives(zk, src, targ); % Hankel part + + zkimag = 1i*zk; + [~, thirdK, ~] = chnk.flex2d.thirdforth_derivatives(zkimag, src, targ); % modified bessel K part + + + + submat = 1/(2*zk^2).*(third(:, :, 1).*(nx.*nx.*nx) + third(:, :, 2).*(3*nx.*nx.*ny) +... + third(:, :, 4).*(3*nx.*ny.*ny) + third(:, :, 6).*(ny.*ny.*ny)) - ... + 1/(2*zk^2).*(thirdK(:, :, 1).*(nx.*nx.*nx) + thirdK(:, :, 2).*(3*nx.*nx.*ny) +... + thirdK(:, :, 4).*(3*nx.*ny.*ny) + thirdK(:, :, 6).*(ny.*ny.*ny)) + ... + (3/(2*zk^2).*(third(:, :, 1).*(nx.*taux.*taux) + third(:, :, 2).*(2*nx.*taux.*tauy + ny.*taux.*taux) +... + third(:, :, 4).*(nx.*tauy.*tauy + 2*ny.*taux.*tauy) + third(:, :, 6).*(ny.*tauy.*tauy))-... + 3/(2*zk^2).*(thirdK(:, :, 1).*(nx.*taux.*taux) + thirdK(:, :, 2).*(2*nx.*taux.*tauy + ny.*taux.*taux) +... + thirdK(:, :, 4).*(nx.*tauy.*tauy + 2*ny.*taux.*tauy) + thirdK(:, :, 6).*(ny.*tauy.*tauy))); + +end + +% second-kernel of evaluation for the clamped plate +if strcmpi(type, 'second kernel') + srcnorm = srcinfo.n; + srctang = srcinfo.d; + nx = repmat(srcnorm(1,:),nt,1); + ny = repmat(srcnorm(2,:),nt,1); + dx = repmat(srctang(1,:),nt,1); + dy = repmat(srctang(2,:),nt,1); + + ds = sqrt(dx.*dx+dy.*dy); + + taux = dx./ds; + tauy = dy./ds; + + [hess, ~, ~] = chnk.flex2d.thirdforth_derivatives(zk, src, targ); % Hankel part + + zkimag = 1i*zk; + [hessK, ~, ~] = chnk.flex2d.thirdforth_derivatives(zkimag, src, targ); % modified bessel K part + + submat = -(1/(2*zk^2).*(hess(:, :, 1).*(nx.*nx) + hess(:, :, 2).*(2*nx.*ny) + hess(:, :, 3).*(ny.*ny))-... + 1/(2*zk^2).*(hessK(:, :, 1).*(nx.*nx) + hessK(:, :, 2).*(2*nx.*ny) + hessK(:, :, 3).*(ny.*ny)))+... + (1/(2*zk^2).*(hess(:, :, 1).*(taux.*taux) + hess(:, :, 2).*(2*taux.*tauy) + hess(:, :, 3).*(tauy.*tauy))-... + 1/(2*zk^2).*(hessK(:, :, 1).*(taux.*taux) + hessK(:, :, 2).*(2*taux.*tauy) + hessK(:, :, 3).*(tauy.*tauy))); % -G_{ny ny} + G_{tauy tauy} +end + + +if strcmpi(type, 'free plate eval first') % G_{ny} + srcnorm = srcinfo.n; + [~,grad] = chnk.flex2d.helmdiffgreen(zk,src,targ); % Hankel part + nx = repmat(srcnorm(1,:),nt,1); + ny = repmat(srcnorm(2,:),nt,1); + + + + zkimag = (1i)*zk; + [~,gradK] = chnk.flex2d.helmdiffgreen(zkimag,src,targ); % modified bessel K part + + + + + submat = (-1/(2*zk^2).*(grad(:, :, 1).*(nx) + grad(:, :, 2).*ny) + ... + 1/(2*zk^2).*(gradK(:, :, 1).*(nx) + gradK(:, :, 2).*ny)); + + +end + +if strcmpi(type, 'free plate eval first hilbert') % G_{tauy} + srctang = srcinfo.d; + coefs = varargin{1}; + + [~,grad] = chnk.flex2d.helmdiffgreen(zk,src,targ); % Hankel part + + + zkimag = (1i)*zk; + [~,gradK] = chnk.flex2d.helmdiffgreen(zkimag,src,targ); % modified bessel K part + + dx = repmat(srctang(1,:),nt,1); + dy = repmat(srctang(2,:),nt,1); + + ds = sqrt(dx.*dx+dy.*dy); + + + taux = dx./ds; % normalization + tauy = dy./ds; + + + + submat = ((1 + coefs(1))/2).*(-1/(2*zk^2).*(grad(:, :, 1).*(taux) + grad(:, :, 2).*tauy) + ... + 1/(2*zk^2).*(gradK(:, :, 1).*(taux) + gradK(:, :, 2).*tauy)); % G_{tauy} +end + + +if strcmpi(type, 'free plate eval second') % G = 1/(2k^2) (i/4 H_0^{1} - 1/2pi K_0) + + [val,~] = chnk.flex2d.helmdiffgreen(zk,src,targ); % Hankel part + + zkimag = (1i)*zk; + [valK,~] = chnk.flex2d.helmdiffgreen(zkimag,src,targ); % modified bessel K part + + submat = 1/(2*zk^2).*val - 1/(2*zk^2).*valK; + +end + + + + + + +% Dirichlet data/potential correpsonding to transmission rep +if strcmpi(type,'trans_rep') + + coef = ones(2,1); + if(nargin == 5); coef = varargin{1}; end; + srcnorm = srcinfo.n; + [submats,grad] = chnk.flex2d.green(zk,src,targ); + nx = repmat(srcnorm(1,:),nt,1); + ny = repmat(srcnorm(2,:),nt,1); + submatd = -(grad(:,:,1).*nx + grad(:,:,2).*ny); + + submat = zeros(1*nt,2*ns); + submat(1:1:1*nt,1:2:2*ns) = coef(1)*submatd; + submat(1:1:1*nt,2:2:2*ns) = coef(2)*submats; +end + +% Neumann data corresponding to transmission rep +if strcmpi(type,'trans_rep_prime') + coef = ones(2,1); + if(nargin == 5); coef = varargin{1}; end; + targnorm = targinfo.n; + srcnorm = srcinfo.n; + + submat = zeros(nt,ns); + + % Get gradient and hessian info + [submats,grad,hess] = chnk.flex2d.green(zk,src,targ); + + nxsrc = repmat(srcnorm(1,:),nt,1); + nysrc = repmat(srcnorm(2,:),nt,1); + nxtarg = repmat((targnorm(1,:)).',1,ns); + nytarg = repmat((targnorm(2,:)).',1,ns); + + % D' + submatdp = -(hess(:,:,1).*nxsrc.*nxtarg ... + + hess(:,:,2).*(nysrc.*nxtarg+nxsrc.*nytarg)... + + hess(:,:,3).*nysrc.*nytarg); + % S' + submatsp = (grad(:,:,1).*nxtarg + grad(:,:,2).*nytarg); + submat = zeros(nt,2*ns) + submat(1:1:1*nt,1:2:2*ns) = coef(1)*submatdp; + submat(1:1:1*nt,2:2:2*ns) = coef(2)*submatsp; +end + + +% Gradient correpsonding to transmission rep +if strcmpi(type,'trans_rep_grad') + coef = ones(2,1); + if(nargin == 5); coef = varargin{1}; end; + + srcnorm = srcinfo.n; + + submat = zeros(nt,ns,6); + % S + [submats,grad,hess] = chnk.flex2d.green(zk,src,targ); + + nxsrc = repmat(srcnorm(1,:),nt,1); + nysrc = repmat(srcnorm(2,:),nt,1); + % D + submatd = -(grad(:,:,1).*nxsrc + grad(:,:,2).*nysrc); + + submat = zeros(2*nt,2*ns); + + submat(1:2:2*nt,1:2:2*ns) = -coef(1)*(hess(:,:,1).*nxsrc + hess(:,:,2).*nysrc); + submat(1:2:2*nt,2:2:2*ns) = coef(2)*grad(:,:,1); + + submat(2:2:2*nt,1:2:2*ns) = -coef(1)*(hess(:,:,2).*nxsrc + hess(:,:,3).*nysrc); + submat(2:2:2*nt,2:2:2*ns) = coef(2)*grad(:,:,2); + + + + + + + + +end + + diff --git a/chunkie/+chnk/+flex2d/kernel_replacement_test.m b/chunkie/+chnk/+flex2d/kernel_replacement_test.m new file mode 100644 index 0000000..993f468 --- /dev/null +++ b/chunkie/+chnk/+flex2d/kernel_replacement_test.m @@ -0,0 +1,44 @@ +clear +clc + + +zk = 6; +nu = 1/3; + +% t1 = 0; +% t2 = pi/10000; + +coefs = [nu;0]; + +srcinfo = []; + +srcinfo.r = [3; 0]; +srcinfo.n = [1; 0]; +srcinfo.d = [0; 1]; + + + +targinfo = []; + + +targinfo.r = [3; 1/1000]; +targinfo.n = [1; 0]; +targinfo.d = [0; 1]; + +fkern = @(s,t) chnk.helm2d.kern(zk, s, t, 'test kernel', coefs); +ans1 = fkern(srcinfo, targinfo); + +fkern1 = @(s,t) chnk.helm2d.kern(zk, s, t, 'free plate K21 first part', coefs); +ans2 = fkern1(srcinfo, targinfo); + + +abs(ans1-ans2) + + + + + + + + + diff --git a/chunkie/+chnk/+helm2d/smooth_part_derivatives.m b/chunkie/+chnk/+flex2d/smooth_part_derivatives.m similarity index 100% rename from chunkie/+chnk/+helm2d/smooth_part_derivatives.m rename to chunkie/+chnk/+flex2d/smooth_part_derivatives.m diff --git a/chunkie/+chnk/+flex2d/test_helmdiffgreen.m b/chunkie/+chnk/+flex2d/test_helmdiffgreen.m new file mode 100644 index 0000000..96e1db1 --- /dev/null +++ b/chunkie/+chnk/+flex2d/test_helmdiffgreen.m @@ -0,0 +1,80 @@ +%test_helmdiffgreen test the difference kernel functions +% using finite differences +% + +clearvars; clc +rng(1234); + +ns = 3; +nt = 4; +k = 1.3; + +src0 = randn(2,ns); +targ0 = randn(2,nt); + +% test derivatives + +[val0,grad0,hess0,der30,der40] = chnk.helm2d.helmdiffgreen(k,src0,targ0); + +for j = 1:5 + h = 10^(-j); + dx = h*[1;0]; + targ1 = targ0 + dx; + [val1,grad1,hess1,der31,der41] = chnk.helm2d.helmdiffgreen(k,src0,targ1); + + errdx = norm(ones(size(val0)) - 2*(val1-val0)/h./(grad0(:,:,1)+grad1(:,:,1))); + errdxx = norm(ones(size(val0)) - 2*(grad1(:,:,1)-grad0(:,:,1))/h./(hess0(:,:,1)+hess1(:,:,1))); + errdxy = norm(ones(size(val0)) - 2*(grad1(:,:,2)-grad0(:,:,2))/h./(hess0(:,:,2)+hess1(:,:,2))); + errdxxx = norm(ones(size(val0)) - 2*(hess1(:,:,1)-hess0(:,:,1))/h./(der30(:,:,1)+der31(:,:,1))); + errdxxy = norm(ones(size(val0)) - 2*(hess1(:,:,2)-hess0(:,:,2))/h./(der30(:,:,2)+der31(:,:,2))); + errdxyy = norm(ones(size(val0)) - 2*(hess1(:,:,3)-hess0(:,:,3))/h./(der30(:,:,3)+der31(:,:,3))); + errdxxxx = norm(ones(size(val0)) - 2*(der31(:,:,1)-der30(:,:,1))/h./(der40(:,:,1)+der41(:,:,1))); + errdxxxy = norm(ones(size(val0)) - 2*(der31(:,:,2)-der30(:,:,2))/h./(der40(:,:,2)+der41(:,:,2))); + errdxxyy = norm(ones(size(val0)) - 2*(der31(:,:,3)-der30(:,:,3))/h./(der40(:,:,3)+der41(:,:,3))); + errdxyyy = norm(ones(size(val0)) - 2*(der31(:,:,4)-der30(:,:,4))/h./(der40(:,:,4)+der41(:,:,4))); + + dx = h*[0;1]; + targ1 = targ0 + dx; + [val1,grad1,hess1,der31,der41] = chnk.helm2d.helmdiffgreen(k,src0,targ1); + + errdy = norm(ones(size(val0)) - 2*(val1-val0)/h./(grad0(:,:,2)+grad1(:,:,2))); + errdyy = norm(ones(size(val0)) - 2*(grad1(:,:,2)-grad0(:,:,2))/h./(hess0(:,:,3)+hess1(:,:,3))); + errdyyy = norm(ones(size(val0)) - 2*(hess1(:,:,3)-hess0(:,:,3))/h./(der30(:,:,4)+der31(:,:,4))); + errdyyyy = norm(ones(size(val0)) - 2*(der31(:,:,4)-der30(:,:,4))/h./(der40(:,:,5)+der41(:,:,5))); + + fprintf('%5.2e : err in dx\n',errdx) + fprintf('%5.2e : err in dy\n',errdy) + fprintf('%5.2e : err in dxx\n',errdxx) + fprintf('%5.2e : err in dxy\n',errdxy) + fprintf('%5.2e : err in dyy\n',errdyy) + fprintf('%5.2e : err in dxxx\n',errdxxx) + fprintf('%5.2e : err in dxxy\n',errdxxy) + fprintf('%5.2e : err in dxyy\n',errdxyy) + fprintf('%5.2e : err in dyyy\n',errdyyy) + fprintf('%5.2e : err in dxxxx\n',errdxxxx) + fprintf('%5.2e : err in dxxxy\n',errdxxxy) + fprintf('%5.2e : err in dxxyy\n',errdxxyy) + fprintf('%5.2e : err in dxyyy\n',errdxyyy) + fprintf('%5.2e : err in dyyyy\n',errdyyyy) +end + +%% +% some high precision calcs for comparison + +% from mathematica +srct = [0;0]; targt = 10^(-4)*[1;1]; +kt = sqrt(2); +valt = -0.03670784159258519723+0.24999999750000000625 *1i; +gradxt = -0.00014535819351409084213-0.00002499999987500000021 *1i; +hessxyt = 0.079577479012801035120+1.249999995833*1e-9*1i; +der3xxyt = 0.000070689659841738687819+0.000012499999916666666823*1i; +der3xxxt = 1591.54965094568054413183407228+0.00003749999983333333359375*1i; +der4xxyyt = 7.95774786149136280529936105345*1e6+0.12499999875000000364583*1i; +[val,grad,hess,der3,der4] = chnk.helm2d.helmdiffgreen(kt,srct,targt); + +abs(val-valt)/abs(valt) +abs(grad(1)-gradxt)/abs(gradxt) +abs(hess(2)-hessxyt)/abs(hessxyt) +abs(der3(1)-der3xxxt)/abs(der3xxxt) +abs(der3(2)-der3xxyt)/abs(der3xxyt) % this derivative appears to be just hard to get +abs(der4(3)-der4xxyyt)/abs(der4xxyyt) diff --git a/chunkie/+chnk/+flex2d/test_ps_besseldiff.m b/chunkie/+chnk/+flex2d/test_ps_besseldiff.m new file mode 100644 index 0000000..bec9af0 --- /dev/null +++ b/chunkie/+chnk/+flex2d/test_ps_besseldiff.m @@ -0,0 +1,16 @@ + +% test bessel0-1 and other ps coefficients + +clearvars; clc + +gam = 0.57721566490153286060651209; + +nterms = 12; +[cf1,cf2] = besseldiff_etc_pscoefs(nterms); + +z = 1.3 + 0.1*1i; +v1 = even_pseval(cf1,z); +abs(v1-(besselj(0,z)-1)) + +v2 = even_pseval(cf2,z); +abs(besselh(0,z)-besselj(0,z)-1i*(2/pi)*((log(z/2)+gam)*besselj(0,z) + v2)) diff --git a/chunkie/+chnk/+helm2d/transmission_helper.m b/chunkie/+chnk/+flex2d/transmission_helper.m similarity index 100% rename from chunkie/+chnk/+helm2d/transmission_helper.m rename to chunkie/+chnk/+flex2d/transmission_helper.m diff --git a/chunkie/+chnk/+helm2d/thirdforth_derivatives.m b/chunkie/+chnk/+helm2d/thirdforth_derivatives.m deleted file mode 100644 index 3342611..0000000 --- a/chunkie/+chnk/+helm2d/thirdforth_derivatives.m +++ /dev/null @@ -1,178 +0,0 @@ -function [hess, third, forth] = thirdforth_derivatives(k,src,targ) -%CHNK.HELM2D.GREEN evaluate the second, third, and forth derivatives of -% helmholtz green's function: i/4 H_0^1(k|x-y|). -% for the given sources and targets - -% the third derivatives are: (listed in order) -% G_{x1x1y1}, G_{x1x1y2}, G_{x1x2y1}, G_{x1x2y2}, G_{x2x2y1}, G_{x2x2y2} -% -% forth derivatives are: -%G_{x1x1x1y1}, G_{x1x1x1y2}, G_{x1x1x2y1}, G_{x1x1x2y2}, G_{x1x2x2y1}, -%G_{x1x2x2y2}, G_{x2x2x2y1}, G_{x2x2x2y2} - -[~,ns] = size(src); -[~,nt] = size(targ); - -xs = repmat(src(1,:),nt,1); -ys = repmat(src(2,:),nt,1); - -xt = repmat(targ(1,:).',1,ns); -yt = repmat(targ(2,:).',1,ns); - -rx = xt-xs; -ry = yt-ys; - -rx2 = rx.*rx; -ry2 = ry.*ry; - -r2 = rx2+ry2; - -r = sqrt(r2); - -[m,n] = size(xs); - -if nargout > 0 - h0 = besselh(0,1,k*r); - hess = zeros(m,n,3); - h1 = besselh(1,1,k*r); - r3 = r.^3; - - h2 = 2*h1./(k*r)-h0; - - hess(:,:,1) = 0.25*1i*k*((rx-ry).*(rx+ry).*h1./r3 - k*rx2.*h0./r2); - hess(:,:,2) = 0.25*1i*k*k*rx.*ry.*h2./r2; - hess(:,:,3) = 0.25*1i*k*((ry-rx).*(rx+ry).*h1./r3 - k*ry2.*h0./r2); -end - - -if nargout > 1 - h3 = besselh(0,1,k*r) - besselh(2, 1, k*r); - h4 = h1 - besselh(3, 1, k*r); - - third = zeros(m,n,6); - - third(:,:,1) = 3i*k*0.25*h1.*(rx.^3)./(r.^5) - 3i*k*0.25*rx.*h1./(r.^3) - ... - 3i*(k^2)*0.25*(rx.^3).*h3./(2*r.^4) + 3i*(k^2)*0.25.*rx.*h3./(2*r.^2) - 1i*(k^3)*0.25*(rx.^3).*h1./(2*r.^3) -... - 1i*(k^3)*0.25.*(rx.^3).*h4./(4*r.^3); - - third(:,:,2) = 3i*k*0.25*(rx.^2).*ry.*h1./(r.^5) - 1i*0.25*k.*ry.*h1./(r.^3) - ... - (1i)*(0.25)*3*(k^2)*(rx.^2).*(ry).*h3./(2*r.^4) + 1i*(k^2)*0.25*ry.*h3./(2*r.^2) - ... - 1i*(k^3)*0.25*(rx.^2).*ry.*h1./(2*r.^3) - 1i*(k^3)*0.25*(rx.^2).*ry.*h4./(4*r.^3); - - third(:,:, 3) = 3i*k*0.25*(rx.^2).*ry.*h1./(r.^5) - 1i*0.25*k.*ry.*h1./(r.^3) - ... - (1i)*(0.25)*3*(k^2)*(rx.^2).*(ry).*h3./(2*r.^4) + 1i*(k^2)*0.25*ry.*h3./(2*r.^2) - ... - 1i*(k^3)*0.25*(rx.^2).*ry.*h1./(2*r.^3) - 1i*(k^3)*0.25*(rx.^2).*ry.*h4./(4*r.^3); - - third(:, :, 4) = -1i*k*0.25*rx.*h1./(r.^3) + 3i*k*0.25.*rx.*(ry.^2).*h1./(r.^5) + ... - 1i*(k^2).*rx.*h3./(8*r.^2) - 3i*(k^2).*rx.*(ry.^2).*h3./(8*r.^4) - ... - 1i*(k^3).*(rx).*(ry.^2).*h1./(8*r.^3) - 1i*(k^3)*rx.*(ry.^2).*h4./(16*r.^3); - - third(:, :, 5) = -1i*k*0.25.*rx.*h1./(r.^3) + 3i*k*0.25.*rx.*(ry.^2).*h1./(r.^5) + ... - 1i*(k^2).*rx.*h3./(8*r.^2) - 3i*(k^2).*rx.*(ry.^2).*h3./(8*r.^4) - ... - 1i*(k^3).*(rx).*(ry.^2).*h1./(8*r.^3) - 1i*(k^3)*rx.*(ry.^2).*h4./(16*r.^3); - - third(:, :, 6) = -3i*k*0.25.*ry.*h1./(r.^3) + 3i*0.25*k*(ry.^3).*h1./(r.^5) + ... - 3i*0.25*(k^2)*ry.*h3./(2*r.^2) - 3i*(k^2)*(0.25)*(ry.^3).*h3./(2*r.^4) - ... - 1i*(0.25)*(k^3).*(ry.^3).*h1./(2*r.^3) - 1i*(k^3)*0.25.*(ry.^3).*h4./(4*r.^3); -end - -if nargout > 2 - h5 = besselh(2, 1, k*r) - besselh(4, 1, k*r); - - forth = zeros(m,n,8); - - forth(:, :, 1) = 0.25*(1i)*(-15)*k.*(rx.^4).*h1./(r.^7) + 0.25*(1i)*18*k*(rx.^2).*h1./(r.^5)-... - 3*(1i)*(0.25)*k*h1./(r.^3) + 0.25*(1i)*(15*k^2).*(rx.^4).*h3./(2*r.^6) - ... - (1i)*(0.25)*9*(k^2)*(rx.^2).*h3./(r.^4) + (1i)*(0.25)*3*(k^2).*h3./(2*r.^2)+... - (1i)*(0.25)*(k^3)*(rx.^4).*h1./(r.^5) + (1i)*(0.25)*(k^3)*(rx.^4).*h4./(2*r.^5) + ... - -(1i)*(0.25)*(k^3)*(rx.^2).*h1./(r.^3) - (1i)*(0.25)*(k^3)*(rx.^2).*h4./(2*r.^3) +... - (1i)*(0.25)*(3*k^3)*(rx.^4).*h1./(2*r.^5) + (1i)*(0.25)*(3*k^3)*(rx.^4).*h4./(4*r.^5) -... - (1i)*(0.25)*(3*k^3)*(rx.^2).*h1./(2*r.^3) - (1i)*(0.25)*(3*k^3)*(rx.^2).*h4./(4*r.^3) +... - (1i)*(0.25)*(0.5)*(k^3)*(rx.^4).*h1./(r.^5) - (1i)*(0.25)*(0.5)*(k^3)*(rx.^2).*h1./(r.^3)-... - (1i)*(0.25)*(0.5)*(k^4)*(rx.^4).*h3./(2*r.^4) + (1i)*(0.25)*(0.5)*(k^3)*(rx.^4).*h4./(2*r.^5) -... - (1i)*(0.25)*(0.5)*(k^3)*(rx.^2).*h4./(2*r.^3) + (1i)*0.25*(0.5)*(k^4)*(-rx.^4).*h3./(4*r.^4) +... - (1i)*0.25*(0.5)*(k^4)*(rx.^4).*h5./(4*r.^4); - - forth(:, :, 2) = (0.25)*(1i)*(-15)*k.*(rx.^3).*(ry).*h1./(r.^7) + 0.25*(1i)*9*k*(rx.*ry).*h1./(r.^5) + ... - (0.25)*(1i)*(15)*(k^2).*(rx.^3).*(ry).*(h3)./(2*r.^6) - (0.25)*(1i)*9*(k^2)*(rx.*ry).*h3./(2*r.^4) -... - ((0.25)*(1i)*(k^3)*(-rx.^3).*(ry).*h1./(r.^5) + 0.25*(1i)*(k^3)*(-rx.^3).*ry.*h4./(2*r.^5)) + ... - (0.25)*(1i)*(3*k^3).*(rx.^3).*ry.*h1./(2*r.^5) + (0.25)*(1i)*(3*k^3).*(rx.^3).*ry.*h4./(4*r.^5) -... - (0.25)*(1i)*(3*k^3)*(rx.*ry).*h1./(2*r.^3) -(0.25)*(1i)*(3*k^3).*(rx.*ry).*h4./(4*r.^3) +... - (0.25)*(1i)*(k^3)*(rx.^3).*ry.*h1./(2*r.^5) - (0.25)*(1i)*(k^4).*(rx.^3).*ry.*h3./(4*r.^4) + ... - (0.25)*(1i)*(k^3)*(rx.^3).*ry.*h4./(4*r.^5) + (0.25)*(1i)*(k^4).*(rx.^3).*(-ry).*h3./(8*r.^4) + ... - (0.25)*(1i)*(k^4)*(rx.^3).*ry.*h5./(8*r.^4); - - forth(:, :, 3) = (0.25)*(1i)*(-15)*k.*(rx.^3).*(ry).*h1./(r.^7) + 0.25*(1i)*9*k*(rx.*ry).*h1./(r.^5) + ... - (0.25)*(1i)*(15)*(k^2).*(rx.^3).*(ry).*(h3)./(2*r.^6) - (0.25)*(1i)*9*(k^2)*(rx.*ry).*h3./(2*r.^4) -... - ((0.25)*(1i)*(k^3)*(-rx.^3).*(ry).*h1./(r.^5) + 0.25*(1i)*(k^3)*(-rx.^3).*ry.*h4./(2*r.^5)) + ... - (0.25)*(1i)*(3*k^3).*(rx.^3).*ry.*h1./(2*r.^5) + (0.25)*(1i)*(3*k^3).*(rx.^3).*ry.*h4./(4*r.^5) -... - (0.25)*(1i)*(3*k^3)*(rx.*ry).*h1./(2*r.^3) -(0.25)*(1i)*(3*k^3).*(rx.*ry).*h4./(4*r.^3) +... - (0.25)*(1i)*(k^3)*(rx.^3).*ry.*h1./(2*r.^5) - (0.25)*(1i)*(k^4).*(rx.^3).*ry.*h3./(4*r.^4) + ... - (0.25)*(1i)*(k^3)*(rx.^3).*ry.*h4./(4*r.^5) + (0.25)*(1i)*(k^4).*(rx.^3).*(-ry).*h3./(8*r.^4) + ... - (0.25)*(1i)*(k^4)*(rx.^3).*ry.*h5./(8*r.^4); - - forth(:, :, 4) = (0.25)*(1i)*(3*k)*(rx.^2).*h1./(r.^5) -(0.25)*(1i)*k.*h1./(r.^3) - ... - (0.25)*(1i)*15*k.*(rx.^2).*(ry.^2).*(h1)./(r.^7) + (0.25)*(1i)*(3*k).*(ry.^2).*h1./(r.^5)-... - (0.25)*(1i)*(3*k^2).*(rx.^2).*h3./(2*r.^4) + (0.25)*(1i)*(k^2).*h3./(2*r.^2) +... - (0.25)*(1i)*15*(k^2).*(rx.^2).*(ry.^2).*h3./(2*r.^6) - (0.25)*(1i)*(3*k^2)*(ry.^2).*h3./(2*r.^4) -... - ((0.25)*(1i)*(k^3).*(-rx.^2).*(ry.^2).*h1./(r.^5) + (0.25)*(1i)*(k^3).*(-rx.^2).*(ry.^2).*h4./(2*r.^5)) +... - (0.25)*(1i)*(3*k^3).*(rx.^2).*(ry.^2).*h1./(2*r.^5) + (0.25)*(1i)*(3*k^3).*(rx.^2).*(ry.^2).*h4./(4*r.^5) - ... - (0.25)*(1i)*(k^3).*(ry.^2).*h1./(2*r.^3) - (0.25)*(1i)*(k^3).*(ry.^2).*h4./(4*r.^3) +... - (0.25)*(1i)*(k^3).*(-rx.^2).*h1./(2*r.^3) + (0.25)*(1i)*(k^3).*(rx.^2).*(ry.^2).*h1./(2*r.^5)-... - (0.25)*(1i)*(k^4).*(rx.^2).*(ry.^2).*h3./(4*r.^4) - (0.25)*(1i)*(k^3).*(rx.^2).*h4./(4*r.^3) + ... - (0.25)*(1i)*(k^3).*(rx.^2).*(ry.^2).*h4./(4*r.^5) -(0.25)*(1i)*(k^4)*(rx.^2).*(ry.^2).*h3./(8*r.^4)+... - (0.25)*(1i)*(k^4).*(rx.^2).*(ry.^2).*h5./(8*r.^4); - - forth(:, :, 5) = (0.25)*(1i)*(3*k)*(rx.^2).*h1./(r.^5) -(0.25)*(1i)*k.*h1./(r.^3) - ... - (0.25)*(1i)*15*k.*(rx.^2).*(ry.^2).*(h1)./(r.^7) + (0.25)*(1i)*(3*k).*(ry.^2).*h1./(r.^5)-... - (0.25)*(1i)*(3*k^2).*(rx.^2).*h3./(2*r.^4) + (0.25)*(1i)*(k^2).*h3./(2*r.^2) +... - (0.25)*(1i)*15*(k^2).*(rx.^2).*(ry.^2).*h3./(2*r.^6) - (0.25)*(1i)*(3*k^2)*(ry.^2).*h3./(2*r.^4) -... - ((0.25)*(1i)*(k^3).*(-rx.^2).*(ry.^2).*h1./(r.^5) + (0.25)*(1i)*(k^3).*(-rx.^2).*(ry.^2).*h4./(2*r.^5)) +... - (0.25)*(1i)*(3*k^3).*(rx.^2).*(ry.^2).*h1./(2*r.^5) + (0.25)*(1i)*(3*k^3).*(rx.^2).*(ry.^2).*h4./(4*r.^5) - ... - (0.25)*(1i)*(k^3).*(ry.^2).*h1./(2*r.^3) - (0.25)*(1i)*(k^3).*(ry.^2).*h4./(4*r.^3) +... - (0.25)*(1i)*(k^3).*(-rx.^2).*h1./(2*r.^3) + (0.25)*(1i)*(k^3).*(rx.^2).*(ry.^2).*h1./(2*r.^5)-... - (0.25)*(1i)*(k^4).*(rx.^2).*(ry.^2).*h3./(4*r.^4) - (0.25)*(1i)*(k^3).*(rx.^2).*h4./(4*r.^3) + ... - (0.25)*(1i)*(k^3).*(rx.^2).*(ry.^2).*h4./(4*r.^5) -(0.25)*(1i)*(k^4)*(rx.^2).*(ry.^2).*h3./(8*r.^4)+... - (0.25)*(1i)*(k^4).*(rx.^2).*(ry.^2).*h5./(8*r.^4); - - forth(:, :, 6) = (-0.25)*(1i)*k*((-9*rx.*ry./(r.^5) + 15*rx.*(ry.^3)./(r.^7)).*h1 + 2*k*rx.*(ry).*h3./(r.^4)-... - k*(-rx./(r.^3) + 3*rx.*(ry.^2)./(r.^5)).*ry.*h3./(2*r) - 4*k*rx.*(ry.^3).*h3./(r.^6) + ... - k*rx.*(ry).*h3./(2*r.^4) - k*rx.*(ry.^3).*h3./(2*r.^6) - (k^2)*rx.*(ry.^3).*h1./(2*r.^5) - (k^2)*rx.*(ry.^3).*h4./(4*r.^5)-... - (k^2)*rx.*(ry.^3).*h1./(r.^5) - (k^2)*rx.*(ry.^3).*h4./(2*r.^5) + 3*k*rx.*ry.*h3./(2*r.^4) -... - 3*k*rx.*(ry.^3).*h3./(2*r.^6) + (k^2)*rx.*ry.*h1./(2*r.^3) + (k^2)*rx.*ry.*h4./(4*r.^3) - (k^2)*rx.*(ry.^3).*h1./(2*r.^5)-... - (k^2)*rx.*(ry.^3).*h4./(4*r.^5) + (k^2)*rx.*(ry).*h1./(2*r.^3) +(k^2)*rx.*ry.*h4./(4*r.^3) - ... - (k^2)*rx.*(ry.^3).*h1./(2*r.^5) - (k^2)*rx.*(ry.^3).*h4./(4*r.^5) + (k^2)*rx.*ry.*h1./(2*r.^3)-... - (k^2)*rx.*(ry.^3).*h1./(2*r.^5) +(k^3)*rx.*(ry.^3).*h3./(4*r.^4) + (k^2)*rx.*(ry).*h4./(4*r.^3)-... - (k^2)*rx.*(ry.^3).*h4./(4*r.^5) - (k^3)*rx.*(-ry.^3).*h3./(8*r.^4) - (k^3)*rx.*(ry.^3).*h5./(8*r.^4)); - - forth(:, :, 7) = (-0.25)*(1i)*k*((-9*rx.*ry./(r.^5) + 15*rx.*(ry.^3)./(r.^7)).*h1 + 2*k*rx.*(ry).*h3./(r.^4)-... - k*(-rx./(r.^3) + 3*rx.*(ry.^2)./(r.^5)).*ry.*h3./(2*r) - 4*k*rx.*(ry.^3).*h3./(r.^6) + ... - k*rx.*(ry).*h3./(2*r.^4) - k*rx.*(ry.^3).*h3./(2*r.^6) - (k^2)*rx.*(ry.^3).*h1./(2*r.^5) - (k^2)*rx.*(ry.^3).*h4./(4*r.^5)-... - (k^2)*rx.*(ry.^3).*h1./(r.^5) - (k^2)*rx.*(ry.^3).*h4./(2*r.^5) + 3*k*rx.*ry.*h3./(2*r.^4) -... - 3*k*rx.*(ry.^3).*h3./(2*r.^6) + (k^2)*rx.*ry.*h1./(2*r.^3) + (k^2)*rx.*ry.*h4./(4*r.^3) - (k^2)*rx.*(ry.^3).*h1./(2*r.^5)-... - (k^2)*rx.*(ry.^3).*h4./(4*r.^5) + (k^2)*rx.*(ry).*h1./(2*r.^3) +(k^2)*rx.*ry.*h4./(4*r.^3) - ... - (k^2)*rx.*(ry.^3).*h1./(2*r.^5) - (k^2)*rx.*(ry.^3).*h4./(4*r.^5) + (k^2)*rx.*ry.*h1./(2*r.^3)-... - (k^2)*rx.*(ry.^3).*h1./(2*r.^5) +(k^3)*rx.*(ry.^3).*h3./(4*r.^4) + (k^2)*rx.*(ry).*h4./(4*r.^3)-... - (k^2)*rx.*(ry.^3).*h4./(4*r.^5) - (k^3)*rx.*(-ry.^3).*h3./(8*r.^4) - (k^3)*rx.*(ry.^3).*h5./(8*r.^4)); - - forth(:, :, 8) = (-0.25)*(1i)*(3*k)*h1./(r.^3) + 0.25*(1i)*(18*k)*(ry.^2).*h1./(r.^5) - (0.25)*(1i)*(15*k).*(ry.^4).*h1./(r.^7) +... - (0.25)*(1i)*(3*k^2).*(h3)./(2*r.^2) - (0.25)*(1i)*(9*k^2).*(ry.^2).*h3./(r.^4) + (0.25)*(1i)*(15*k^2).*(ry.^4).*h3./(2*r.^6) +... - (0.25)*(1i)*(k^3)*(-ry.^2).*h1./(r.^3) - (0.25)*(1i)*(k^3).*(ry.^2).*h4./(2*r.^3) - (0.25)*(1i)*(k^3).*(-ry.^4).*h1./(r.^5) + ... - (0.25)*(1i)*(k^3)*(ry.^4).*h4./(2*r.^5) - (0.25)*(1i)*(3*k^3).*(ry.^2).*h1./((2*r.^3)) - (0.25)*(1i)*(3*k^3).*(ry.^2).*h4./(4*r.^3) +... - (0.25)*(1i)*(3*k^3)*(ry.^4).*h1./(2*r.^5) + (0.25)*(1i)*(3*k^3).*(ry.^4).*h4./(4*r.^5) - (0.25)*(1i)*(k^3).*(ry.^2).*h1./(2*r.^3) + ... - (0.25)*(1i)*(k^3)*(ry.^4).*h1./(2*r.^5) - (0.25)*(1i)*(k^4).*(ry.^4).*h3./(4*r.^4) - (0.25)*(1i)*(k^3).*(ry.^2).*h4./(4*r.^3) + ... - (0.25)*(1i)*(k^3)*(ry.^4).*h4./(4*r.^5) + (0.25)*(1i)*(k^4).*(-ry.^4).*h3./(8*r.^4) + (0.25)*(1i)*(k^4)*(ry.^4).*h5./(8*r.^4); - -end -end - - - - - - - - - - -