added code
58 changed files+2798−0
.DS_Storeadded+0−0View file
Binary file not shown.
wfp_1Dspringscattering/.DS_Storeadded+0−0View file
Binary file not shown.
wfp_1Dspringscattering/alphaEvolution_spectral.madded+51−0View file
@@ -0,0 +1,51 @@
1+function [anp1,bnp1] = alphaEvolution_spectral(an,bn,sn,zeroLoc,K,dt,W,snIdxNow,p,q,p0,q0,sn_hat)
2+% ALPHAEVOLUTION_SPECTRAL computes the Fourier coefficients to spectral
3+% accuracy
4+%
5+% [anp1,bnp1] = alphaEvolution_spectral(an,bn,sn,zeroLoc,K,dt,W,snIdxNow,...
6+% p,q,p0,q0,s,N,tol)
7+% returns Fourier coefficients at the new time t = n*dt, anp1,...
8+% and time derivativesof anp1, bnp1.
9+%
10+% Inputs:
11+% snIdxNow : the current time index
12+% an: Fourier coefficients at time t = tInitial + (n-1)*dt
13+% bn: 1st time derivative of an
14+% sn: matrix of density grid functions sn(time, source number)
15+% zeroLoc: location of the zero wave number k
16+% K: vector of fourier wave numbers
17+% N: number of Fourier modes
18+% tol: error tolerance for finufft
19+% dt: time-step
20+% W: W*dt = delta where [0,delta] is the support of phit
21+% s: [s(1),...,s(M)] location of the M spring scatterers
22+% p,q,p0,q0: integrals needed in the evaluation of Fourier coefficients
23+% (see notes)
24+% expMat: expMat = exp(1i*k*(s')) needed in Fourier coefficient evaluation
25+
26+if nargin == 0, test_alphaEvaluation_spectral; return; end
27+
28+%% Treatment of u_history
29+k = K';
30+
31+% update k values not equal to zero
32+i = 0:W;
33+h = sum(p(:,i+1).*sn_hat,2); % sum rows
34+g = sum(q(:,i+1).*sn_hat,2);
35+h = h*dt; g = g*dt;
36+
37+% evolution equations k ~= 0
38+anp1 = an.*cos(k*dt) + bn.*(1./k).*sin(k*dt) + h;
39+bnp1 = -k.*an.*sin(k*dt) + bn.*cos(k*dt) + g;
40+
41+% update k = 0 frequency
42+sig0 = sum(sn(snIdxNow-i,:),2); % sum rows
43+h0 = sum(p0(i+1).*sig0);
44+g0 = sum(q0(i+1).*sig0);
45+h0 = h0*dt; g0 = g0*dt;
46+
47+% evolution equations for k = 0 case
48+anp1(zeroLoc) = an(zeroLoc) + bn(zeroLoc)*dt + h0;
49+bnp1(zeroLoc) = bn(zeroLoc) + g0;
50+
51+end
wfp_1Dspringscattering/alphaRBCproj.madded+147−0View file
@@ -0,0 +1,147 @@
1+function [an,bn,info] = alphaRBCproj(an,bn,phi,phip)
2+% ALPHARBCPROJ apply projector to Fourier coeff rep for 1D free space BCs
3+%
4+% [an,bn] = alphaRBCproj(an,bn,phi,phip) applies RBCs to the alpha_n and beta_n
5+% vectors (set of Fourier coefficients of u and v=u_t, respectively).
6+% phi is a length-w vector giving a (roll-up) window evaluated on
7+% the Nyquist grid, and phip the evaluation of its derivative (scaled as if
8+% the grid were unit sized). w sets the tolerance, and phi and phip could
9+% be as in the history split. The cost is O(N log N), where
10+% N=length(an)=length(bn) is the number of coeffs (N/2 the max frequency).
11+%
12+% Usage:
13+% The user's computational domain is Omega = (-pi+2wh,pi-2wh), where
14+% h = 2pi/N is the Nyquist grid spacing implied by the number of coefficients,
15+% for the periodic domain [-pi,pi). The user simulation must live in Omega.
16+% [-pi,-pi+2wh] and [pi-2wh,pi) are padding
17+% regions needed to apply the RBCs, and waves in them will be destroyed.
18+% Thus the user sacrifices 4w of the N points to the RBCs, in order that
19+% the Fourier series well approximates the free-space solution in Omega.
20+% To maintain outgoing Fourier coefficients, this RBC projection must be
21+% reapplied in tRBC time units, for wh <= tRBC <= 2wh.
22+% The lower limit is needed so that traveling waves have cleared the
23+% windowing region (otherwise the window hits them multiple times and
24+% aliasing error occurs). The upper limit is so that waves do not wrap around
25+% into the other transition region, since the assumption is made that the
26+% waves in transition region (eg (pi-2wh,pi-wh)) are out-going.
27+% Even without this assumption, spectral differentiation would be needed
28+% to recover v in the padding regions, but waves would be hit >1 times by
29+% the window, causing aliasing error. (For t > 3wh,
30+% plain old periodic pollution would also occur back in Omega).
31+% The phi and phip should be as in the history split (as given by window.jl).
32+%
33+% [an,bn,info] = alphaRBCproj(an,bn,phi,phip) also prints & returns debug info.
34+%
35+% When called with no arguments, a self-test is done.
36+%
37+% The method is:
38+% 1) evaluate u and v on the h grid, via the FFT from the coeffs an and bn.
39+% 2) multiplication of u by phi-windows in the padding regions (which
40+% broadens the spectral bandwidth, hence gam is needed),
41+% 3) building v=u_t=+-u_x on the grid points in the padded regions to give
42+% outgoing 1D wave equation solutions there.
43+% 4) inverse FFT to the output coefficients an, bn.
44+%
45+% Notes:
46+% i) N even for now.
47+% ii) The total padding region could easily be shrunk from 4w to 3w points.
48+% The former gives a bit more flexibility in the projection interval.
49+% iii) check Fourier series defn, sign, prefactor, may not matter.
50+
51+% Barnett 9/8/24
52+if nargin==0, test_alphaRBCproj; return; end
53+verb = (nargout>2); % verbosity
54+N = length(an); assert(length(bn)==N)
55+w = length(phi); assert(length(phip)==w)
56+assert(4*w<N)
57+h = 2*pi/N;
58+
59+% 1)
60+u0 = fft(an(:)); v0 = fft(bn(:)); % x on [0,2pi), so padding is in middle
61+
62+% 2)
63+j = N/2-2*w+(1:4*w)'; % indices to bleach out (combine R and L padding)
64+phi=phi(:); phip=phip(:); % force col vecs
65+phij = [flipud(phi); zeros(2*w,1); phi]; % R then L padding
66+pou = ones(N,1); pou(j) = phij; info.pou=pou; % save the POU, full [0,2pi) grid
67+u = u0 .* pou;
68+v = v0 .* pou; % the naive step 3, alone would cause O(1/w) reflection.
69+% This assumes eg, waves in (pi-2wh,pi-wh) are R-going
70+
71+% 3) apply prod rule, overwrite v = (phi.u0)_x = phi.u0_x + phi'.u0 in padding,
72+% noting that phi.u0_x is already accounted for in v.
73+phipj = (1/h) * [flipud(phip); zeros(2*w,1); phip]; % takes d/dx or -d/dx
74+v(j) = v(j) + phipj.*u0(j); % add term +- phi'.u0 in prod rule
75+
76+% 4)
77+an = reshape(ifft(u),size(an)); bn = reshape(ifft(v),size(bn));
78+
79+
80+%%%%%% helpers for test
81+function [u,v,xg] = show_uv(an,bn) % plot Fourier rep in real space on Nyq grid
82+N = length(an);
83+xg = -pi + (0:N-1)/N*2*pi;
84+u = fftshift(fft(an)); v = fftshift(fft(bn)); % shift since grid starts -pi
85+plot(xg, real([u;v]), '.-')
86+a=axis; a(1:2)=[-pi,pi]; axis(a); drawnow
87+
88+function [an, bn] = propexact(t, kg, an0, bn0) % evol by t the coeffs on k grid
89+ak=abs(kg); s=sin(ak*t); c=cos(ak*t); % (taken from tryEFrep.jl)
90+sok=s./ak; sok(kg==0) = t; % note all act on vectors
91+an = c.*an0 + sok.*bn0;
92+bn = -s.*ak.*an0 + c.*bn0;
93+
94+%%%%%%
95+function test_alphaRBCproj % test can evolve a source-free WE in free space
96+addpath ./utils
97+verb = 1; % 1,2,.. for debug tests
98+tol = 1e-12;
99+gam = 0.5; % fractional exceeding of Nyquist of the sigma
100+
101+theta = log(1/tol); % set up window on regular grid in Nour's way
102+w = ceil(2*theta/(pi*gam));
103+phipfun = window(1.0,theta,w); % hack dt=1 (not so happy window needs dt)
104+g = (1:w) - 0.5; % not sure of off-by-one shifts here
105+phip = phipfun(g)'; % col vec
106+phi = blending(phipfun, g, tol);
107+if verb>1
108+ figure; plot(g, [phi, phip], '-+'); axis tight; xlabel('gridpoints');
109+ legend('\phi','\phi'''); drawnow
110+ % check that phip is really deriv of phi on the same grid...
111+ foldphipe = [phip;flipud(-phip)]; % append flipped so can do periodic diff
112+ foldphip = (pi/w) * perispecdiff([phi;flipud(phi)]); % 2pi/2w rescales t
113+ fprintf('max err in phip on grid: %.3g\n', norm(foldphip-foldphipe,inf))
114+end
115+
116+N = 300; assert(mod(N,2)==0) % grid size N even
117+x0 = 0.8; % test initial condition center loc
118+k0 = N/sqrt(8*theta); % k-width of Gaussian which hits tol by N/2
119+kg = [0:N/2-1, -N/2:-1]; % freq k grid, fft ordering (may differ from Nour)
120+c0 = 2/(sqrt(2*pi)*k0); % Gaussian height 2 (so splits into packets of 1)
121+an = c0 * exp(-0.5*(kg/k0).^2) .* exp(1i*x0*kg); % bump at x0
122+bn = 0*an;
123+h = 2*pi/N; % implied grid (used inside proj only, plus for plotting)
124+tRBC = 1.0*w*h; % may check fails either side of allowed range [1,2] :)
125+fprintf('%.3g of N is in Omega; tRBC=%.3g\n',1-4*w/N,tRBC)
126+T = 5*pi;
127+nt = ceil(T/tRBC)+1; tg=tRBC*(0:nt-1); uxt = nan(N,nt);
128+info.pou = nan(1,N);
129+if verb, figure; end
130+for n=1:nt, t = tg(n); % t step
131+ if verb, subplot(2,1,1); [u,v,xg] = show_uv(an,bn);
132+ if verb>1, hold on; plot(xg,fftshift(info.pou),'k-'); hold off; end
133+ title(sprintf('u and v: t=%6.2f',t));
134+ vline([-1 1]*(pi-2*w*h)); vline([-1 1]*(pi-w*h)); % POU on/off pts
135+ uxt(:,n) = u(:);
136+ subplot(2,1,2); imagesc(xg,tg,log10(abs(uxt))'); axis xy;
137+ caxis([-12 0]); xlabel x; ylabel t; colorbar; title('log_{10}|u(x,t)|')
138+ vline([-1 1]*(pi-2*w*h)); vline([-1 1]*(pi-w*h)); % POU on/off pts
139+ drawnow;
140+ end
141+ [an,bn,info] = alphaRBCproj(an,bn,phi,phip); % comment to switch off
142+ [an, bn] = propexact(tRBC, kg, an, bn); % we plot *after* prop
143+ pause(.5); % anim already too slow anyway :(
144+end
145+
146+
147+
wfp_1Dspringscattering/checkMatchingGrids.madded+14−0View file
@@ -0,0 +1,14 @@
1+function checkMatchingGrids(t_cell,x_cell,numResolutions)
2+
3+tGridErr = 0; xGridErr = 0;
4+for m = 2:numResolutions
5+ tGridErr = max(tGridErr,maxNorm(t_cell{m} - t_cell{1}));
6+ xGridErr = max(xGridErr,maxNorm(x_cell{m} - x_cell{1}));
7+end
8+
9+tol = 1e-12;
10+if(tGridErr>tol || xGridErr>tol)
11+ warning('Time or space grids do not match to report self-convergence');
12+end
13+
14+end
wfp_1Dspringscattering/data/.DS_Storeadded+0−0View file
Binary file not shown.
wfp_1Dspringscattering/fig/.DS_Storeadded+0−0View file
Binary file not shown.
wfp_1Dspringscattering/fixDtBasedOnTypNeighborsNum.madded+35−0View file
@@ -0,0 +1,35 @@
1+function [dt,typNumOfNeighbors] = fixDtBasedOnTypNeighborsNum(maxNumNeighbors,src_dmn,W,M,dt)
2+% FIXDTBASEDONTYPNEIGHBORSNUM adjusts dt to make sure num of local
3+% neighbors is less than a maximum value
4+%
5+% INPUTS:
6+% maxNumNeighbors: maximum number of neighbors
7+% [as,bs]: domain in which the sources are located
8+% W: W*dt is the compact support of the window function
9+% M: number of spring scatterers
10+% dt: time step (first estimate)
11+% N: 2pi/dt: number of Fourier modes
12+% tFinal: final time.
13+%
14+% OUTPUTS:
15+% dt: adjusted timestep
16+% N: adjusted number of Fourier modes
17+% typNumOfNeighbors: typical number of neighbors expected on average for
18+% each source.
19+
20+as = src_dmn(1); bs = src_dmn(2);
21+
22+if(isnan(maxNumNeighbors))
23+ maxNumNeighbors = ceil(sqrt(M*log(2*W*M))); % choose to balance otherContrib and finufft
24+end
25+
26+typNumOfNeighbors = ceil(2*W*dt*M/(bs - as));
27+
28+if(typNumOfNeighbors>maxNumNeighbors)
29+ dt = (bs - as)*maxNumNeighbors/(2*W*M);
30+ fprintf('time step changed to dt = %1.2e to ensure maximum number of neighbors %d\n',dt,maxNumNeighbors);
31+ typNumOfNeighbors = maxNumNeighbors;
32+end
33+
34+
35+end
wfp_1Dspringscattering/get_ExactSol.madded+64−0View file
@@ -0,0 +1,64 @@
1+function U = get_ExactSol(x,t,mu,t0,M,s)
2+% GETEXACTSOL gives the exact solutions for a Gaussian densities
3+%
4+% U = get_ExactSol(x,t,mu,t0,M,s,addProj)
5+% Compute an exact solution for testing by choosing Gaussian densities
6+% with at locations 'x' (may be an input vector) and time 't' (one fixed
7+% time value). The Gaussians take the form
8+% sig(t) = exp(-mu*(t - t0)^2) where 'mu' and 't0' are vectors of size 'M'
9+% where 'M' is the number of sources and 's' represents the source
10+% locations s = [s(1),...,s(M)]
11+%
12+% If addProj == 1, the exact solution computed is the one corresponding to
13+% the solution of the free-space problem
14+% If addProj == 0, the exact solution is the one corresponding to the
15+% solution of the periodic problem
16+
17+if nargin == 0, test_getExactSol; return; end
18+
19+Ng = length(x);
20+U = zeros(Ng,1);
21+
22+sqrtpi = sqrt(pi);
23+sqrtmu = sqrt(mu);
24+
25+%%% Without the images
26+for i = 1:Ng % for each x value
27+ for m = 1:M % for each source
28+ if(t>abs(x(i) - s(m))) % check needed for integration
29+ A = t - abs(x(i) - s(m));
30+ U(i) = U(i) + (sqrtpi/(4*sqrtmu(m)))*(erf(sqrtmu(m)*t0(m)) - erf(sqrtmu(m)*(t0(m) - A)));
31+ end
32+ end
33+end
34+
35+end % end function
36+
37+function test_getExactSol
38+ax = -pi; bx =pi; Nx = 50;
39+as = -pi; bs = pi;
40+dt = 1e-1;
41+t = 0:dt:4*pi; Nt = length(t);
42+M = 100; Ns = M;
43+x = linspace(ax,bx,Nx);
44+s = linspace(as,bs,Ns); s = s(randperm(Ns,M));
45+mu = linspace(10,50,M);
46+t0 = linspace(0.5,t(end)-0.5,M);
47+
48+u = zeros(Nx,Nt);
49+for n = 1:Nt
50+ u(:,n) = get_ExactSol(x(:),t(n),mu,t0,M,s);
51+end
52+
53+u = real(u);
54+
55+imagesc(x,t,log10(abs(u))'); axis xy;
56+% hold on; plot(s,0,'|r','MarkerSize',10); hold off;
57+a=axis; a(1:2)=[-pi,pi]; axis(a);
58+xlabel x; ylabel t; colorbar; title('u(x,t)');
59+
60+% plot(x,U,'LineWidth',2);
61+% title(sprintf('solution evaluated at t = %1.2f',t));
62+% xlabel('x'); ylabel('u'); axis tight;
63+
64+end % end the test function
wfp_1Dspringscattering/get_InterpNodes.madded+151−0View file
@@ -0,0 +1,151 @@
1+function [interpNodes,interpNodesIdx] = get_InterpNodes(auxNode, dt, order, endTime)
2+% GETINTERPNODES retrieves the interpolation nodes on the time-grid to
3+% interpolate unknown function values at auxiliary nodes
4+%
5+% [interpNodes,interpNodesIdx] = get_InterpNodes(auxNode, dt, order, endTime)
6+% returns the interpolation nodes on the time-grid and their
7+% corresponding indices to interpolate function values at a given
8+% auxiliary node (off the time-grid). The function ensures that
9+% interpolation nodes do not go beyong a given endTime
10+%
11+% Inputs:
12+% auxNode: given value that may not coincide with the uniform time-grid
13+% dt: time-step
14+% order: equals the number of interpolation nodes used
15+% endTime: ensure interpolation nodes < endTime
16+%
17+% Outputs:
18+% interpNodes: a vector of size order that carries the interpolation nodes
19+% needed to determine function values at auxNode
20+% interpNodesIdx: a vector of size order that carries the interpolation
21+% nodes indices on the time-grid
22+
23+
24+%%% This is the one-sided treatment (towards the left)
25+% rightGridPoint_idx = ceil(glnode/dt);
26+% rightGridPoint = rightGridPoint_idx*dt; % find the nearest node on the uniform grid
27+% interpNodes = linspace(rightGridPoint - (order - 1)*dt,rightGridPoint,order)';
28+% interpNodesIdx = linspace((rightGridPoint_idx - (order - 1)),rightGridPoint_idx,order)';
29+
30+if nargin==0, test_getInterpNodes; return; end
31+
32+%%% This is the centered/skewed treatment
33+m = order; % set m to equal the order
34+
35+tau_idx = round(auxNode/dt); % index of nearest node on time-grid
36+tau = tau_idx*dt; % nearest node on time-grid
37+endTime_idx = round(endTime/dt); % index of endTime
38+k = endTime_idx - tau_idx; % how far is nearest node from endTime
39+
40+if(mod(m,2) == 0) % if m is even
41+ if(k>=((m-2)/2)) % centered interpolation treatment
42+
43+ if(auxNode>tau)
44+ % rightward skew: add interpolation point to the right end
45+ startingNode = tau - ((m-2)/2)*dt;
46+ startingIdx = tau_idx - ((m-2)/2);
47+
48+ % if last node is outside domain: add interpolation point to
49+ % the left end
50+ if((tau+((m/2)*dt))>endTime)
51+ startingNode = tau - (m/2)*dt;
52+ startingIdx = tau_idx - (m/2);
53+ end
54+
55+ else % leftward skew: add interpolation point to the left end
56+ startingNode = tau - (m/2)*dt;
57+ startingIdx = tau_idx - (m/2);
58+ end
59+ else
60+ % skewed treatment: when centered interpolation is not possible
61+ % because interpolation nodes become larger than endTime
62+ startingNode = tau - (m-k-1)*dt;
63+ startingIdx = tau_idx - (m-k-1);
64+ end
65+else % if m is odd
66+ if(k>=((m-1)/2)) % centered interpolation
67+ startingNode = (tau - ((m-1)/2)*dt);
68+ startingIdx = tau_idx - ((m-1)/2);
69+ else % skewed treatment
70+ startingNode = tau - (m-k-1)*dt;
71+ startingIdx = tau_idx - (m-k-1);
72+ end
73+end
74+
75+% interpolation Nodes: startingNode, startingNode + dt,...
76+% ...,startingNode + (m-1)*dt
77+
78+interpNodes = (startingNode:dt:(startingNode + ((m-1)*dt)))';
79+interpNodesIdx = (startingIdx:(startingIdx + (m-1)))';
80+
81+end
82+
83+function test_getInterpNodes
84+clf;
85+
86+startTime = 0; endTime = 1;
87+t1 = .25; t2 = 1;
88+N = 16;
89+order = 5;
90+saveToFile = 0; % set to '1'
91+saveFile = '/Users/nalhassanieh/Desktop';
92+
93+tn = linspace(startTime,endTime,N);
94+dt = tn(2) - tn(1);
95+glnodes = lgwt(N,t1,t2);
96+
97+% plot integration domain and GL nodes
98+figure(1);
99+hAxes = axes('NextPlot','add',...
100+ 'DataAspectRatio',[1 1 1],...
101+ 'XLim',[startTime endTime],...
102+ 'YLim',[0 eps],...
103+ 'Color','none');
104+p1 = plot(glnodes,0,'bo','MarkerSize',10);
105+p2 = plot(tn,0,'r|','MarkerSize',10);
106+p3 = plot([t1;t2],0,'kx','MarkerSize',10);
107+l_temp = [p1(1);p2(1);p3(1)];
108+% legend(l_temp,'GL nodes','time-grid','endpoints','interpreter','latex','location','north');
109+% title('GL nodes for [t_1,t_2] on time-grid');
110+set(gca,'XTick',[t1 t2], 'YTick', []);
111+xticklabels({'$t_1$','$t_2$'})
112+set(gca,'TickLabelInterpreter','latex');
113+
114+if(saveToFile == 1)
115+ saveas(gcf,sprintf('%s/GLnodesOnTimeGrid',saveFile),'epsc');
116+end
117+
118+pos = {'right','center','left'};
119+
120+cnt = 1;
121+for glnodeNum = [2,8,15]
122+ auxNode = glnodes(glnodeNum);
123+ [interpNodes,~] = get_InterpNodes(auxNode, dt, order, endTime);
124+
125+ % plot one GL node and interpolation nodes
126+ figure;
127+ hAxes = axes('NextPlot','add',...
128+ 'DataAspectRatio',[1 1 1],...
129+ 'XLim',[startTime endTime],...
130+ 'YLim',[0 eps],...
131+ 'Color','none');
132+ p1 = plot(glnodes(glnodeNum),0,'bo','MarkerSize',10);
133+ p2 = plot(interpNodes,0,'*','color',[0.4660 0.6740 0.1880],'MarkerSize',10);
134+ p3 = plot(tn,0,'r|','MarkerSize',10);
135+ p4 = plot([t1;t2],0,'kx','MarkerSize',10);
136+ l_temp = [p1(1);p2(1);p3(1);p4(1)];
137+ % legend(l_temp,'GL node','interpolation nodes','time-grid','endpoints','interpreter','latex','location','north');
138+ % title(sprintf('%s GL node interpolation',pos{cnt}));
139+ set(gca,'XTick',[t1 t2], 'YTick', []);
140+ xticklabels({'$t_1$','$t_2$'})
141+ set(gca,'TickLabelInterpreter','latex');
142+
143+ if(saveToFile == 1)
144+ saveas(gcf,sprintf('%s/exampleGLInterp_%s',saveFile,pos{cnt}),'epsc');
145+ end
146+
147+ pause(0.1);
148+ cnt = cnt + 1;
149+end
150+
151+end
wfp_1Dspringscattering/get_g.madded+11−0View file
@@ -0,0 +1,11 @@
1+function gn = get_g(t,dt,s,dataParam,M,beta,solnType)
2+
3+if(strcmp(solnType,'ms'))
4+ U = get_ExactSol(s,(t + dt),dataParam.mu,dataParam.t0,M,s); % get the manufactured solution at xval
5+ gn = -dataParam.sig(t + dt) - beta.*U;
6+else
7+ uin = get_uIncident(s,(t + dt),dataParam.mu,dataParam.t0);
8+ gn = beta.*uin;
9+end
10+
11+end
wfp_1Dspringscattering/get_gvalue.madded+15−0View file
@@ -0,0 +1,15 @@
1+function g = get_gvalue(t,s,j,sig,beta,M,mu,t0,addProj)
2+% GET_GVALUE get the data at a given time t
3+%
4+% g = get_gvalue(t,s,j,sig,beta,M,mu,t0) returns the data g for each spring
5+% scatterer condition at a time 't'. 's' is a vector containing the source
6+% locations s = [s(1),...,s(M)], where 'M' is the number of sources. 'sig'
7+% is a cell that holds manufactured density values (Gaussians) where
8+% sig{j} = @(t) exp(mu(j)*(t - t0(j))). 'beta' is a factor of the spring
9+% constant, and 'j' represents the index of the current spring scatterer
10+
11+xval = s(j); % location at the current source
12+U = get_ExactSol(xval,t,mu,t0,M,s,addProj); % get the manufactured solution at xval
13+g = -sig(t,j) - beta*U;
14+
15+end
wfp_1Dspringscattering/get_gvalueFromIncomingField.madded+18−0View file
@@ -0,0 +1,18 @@
1+function g = get_gvalueFromIncomingField(t,s,j,beta,mu,t0)
2+% GET_GVALUEFROMINCOMINGFIELD get data from incoming field
3+%
4+% g = get_gvalueFromIncomingField(t,s,j,beta) returns the data g for each spring
5+% scatterer condition at a time 't'. 's' is a vector containing the source
6+% locations s = [s(1),...,s(M)], where 'M' is the number of sources.
7+% 'beta' is a factor of the spring constant, and 'j' represents the index
8+% of the current spring scatterer. 'mu' and 't0' are parameters needed for
9+% the incident wave.
10+
11+x = s(j);
12+uin = get_uIncident(x,t,mu,t0);
13+g = beta*uin;
14+
15+% testing
16+% g = beta*cos(10*(s(j) - t - t0));
17+
18+end
wfp_1Dspringscattering/get_historyContribution.madded+38−0View file
@@ -0,0 +1,38 @@
1+function [anp1,bnp1,historySum,tRBC_now] = get_historyContribution(t,an,bn,sn,zeroLoc,K,...
2+ dt,tRBC, W, snIdxNow, p,q,p0,q0,sn_hat,tRBC_now,phi_tn,phit_tn,s,tol)
3+% GET_HISTORYCONTRIBUTION a function to get the history contribution in BIE
4+%
5+% INPUTs:
6+% t: current time
7+% an,bn: Fourier coefficients
8+% sn,sn_hat: density grid function and transform
9+% zeroLoc: location of the zero frequency
10+% K: grid of frequency values
11+% dt: time-step
12+% tRBC,tRBC_now: time-step for radiation boundary conditions and tRBC_now
13+% is parameter to indicate which value on the time grid is being used at
14+% the current step.
15+% dt; time-step
16+% W: width of the window function
17+% snIdxNow: index of sn at the first time-step
18+% p,q,p0,q0: the integrals p and q needed for the history evaluation
19+% (see notes). p0 and q0 correspond to the zero-frequency case.
20+% phi_tn,phit_tn: blending and window function for free-space projection
21+% s:location of sources
22+% tol: error tol
23+% addProj: button to add free-space projection
24+%
25+% OUTPUTs:
26+% anp1,bnp1: Fourier coefficients at the new time-step
27+% historySum: sum needed in the history evaluation on the right side
28+% tRBC_now: the updated time to apply RBC
29+
30+[anp1,bnp1] = alphaEvolution_spectral(an,bn,sn,zeroLoc,K,dt,W,snIdxNow,p,q,p0,q0,sn_hat);
31+if(t>=tRBC_now)
32+ [anp1,bnp1,~] = alphaRBCproj(anp1,bnp1,phi_tn,phit_tn);
33+ tRBC_now = tRBC_now + tRBC;
34+end
35+
36+historySum = finufft1d2(s,-1,tol,anp1,struct('modeord',1));
37+
38+end
wfp_1Dspringscattering/get_lmmwts.madded+98−0View file
@@ -0,0 +1,98 @@
1+function [wts,LebesgueConstant] = get_lmmwts(xpts, npts, nqloc,rlen,dt,delta,winData,computeLebesgueCst)
2+% GET_LMMWTS gives the weights for the LMM-style treatment of [0,dt]
3+% integrals or [0,L] integrals where L<dt
4+%
5+% wts = get_lmmwts(xpts, npts, nqloc,rlen,phit,blendingTol,dt)
6+% performs barycentric interpolation using 'xpts' (a vector of size 'npts')
7+% as interpolation nodes, then applies GL quadrature with 'nqloc'
8+% quadrature nodes to evaluate the integral over
9+% [xpts(npts-1),xpts(npts-1) + rlen].
10+% The weights are modulated by the blending function needed for the windowing
11+% in the analytic split of the solution. Here 'phi_wts' are cheb weights
12+% to compute the blending function.
13+% 'chebApproxInfo': a struct with fields domain [ax,bx] over which cheb
14+% weights are computed; ninters, number of subintervals, and nord,
15+% order of polynomial interpolation or number of cheb nodes.
16+% 'computeLebesgueConstant' is a button to compute Lebesgue constant if
17+% needed (set to one).
18+%
19+% Note: function rewritten from Fortran function provided by Leslie
20+% Greengard
21+
22+if nargin==0, test_getlmmwts; return; end
23+
24+rintmatloc = zeros(nqloc, npts); % Allocate memory for integration matrix
25+
26+% Obtain Gauss-Legendre nodes and weights -> on [-1,1] interval
27+gl = glwt_prep(nqloc);
28+glnodest = gl.x0;
29+glweightst = gl.w0;
30+
31+% precompute the barycentric weights
32+X=repmat(xpts,1,npts);
33+w = 1./prod(-X+X.'+eye(npts),1);
34+
35+% Loop over intervals (xpts(1),..,xpts(npts))
36+aa = xpts(npts-1); % here the interval of integration [aa,bb] = [xpts(npts-1),xpt(npts)]
37+
38+% GL nodes adjusted to interval
39+xt = aa + rlen*(glnodest + 1.0)/2.0;
40+phi = generalwindow((dt - xt)/delta,winData);
41+phimult = (1 - phi);
42+
43+glwloc = rlen*glweightst.*phimult;
44+
45+tol = 1.0e-20;
46+for iquad = 1:nqloc % loop over each GL node
47+ rr = prod((xt(iquad) - xpts));
48+ rintmatloc(iquad, :) = rr*w./(xt(iquad) - xpts');
49+ rintmatloc(iquad,(abs(xt(iquad) - xpts) < tol)) = 1;
50+end
51+
52+% Compute the Lebesgue constant
53+if(computeLebesgueCst == 1)
54+ LebesgueConstant = norm(rintmatloc,inf);
55+else
56+ LebesgueConstant = NaN;
57+end
58+
59+% Compute weights: GL quadrature applied to each lj
60+wts = glwloc'*rintmatloc;
61+
62+end % end of get_lmmwts function
63+
64+function test_getlmmwts
65+clf;
66+ax = 0; bx = 1;
67+h = 1e-4;
68+xpts = ax:h:bx;
69+rlen = h;
70+tol = 1e-12; % error tolerance
71+gam = .5; % pad gam*Nyquist for window
72+theta = log(1/tol);
73+W = ceil(2*theta/(pi*gam));
74+[phit1,phitt1] = window(theta,W);
75+chebApproxInfo = struct('ninters',12,'nord',5,'domain',[-36,100]);
76+[phi_wts,~,~]= prepWindowChebWts(phit1,phitt1,chebApproxInfo,tol);
77+
78+% npts = 7;
79+% nqloc = npts;
80+% [wts,LebesgueConstant] = get_lmmwts(xpts, npts, nqloc,rlen,phi_wts,h,chebApproxInfo);
81+% bar(wts);
82+% title(sprintf('Order %2d LMM weights for grid points\n with separation %1.2e',npts,h));
83+
84+p = 8;
85+for i = 1:p
86+ npts = 2*i;
87+ nqloc = npts;
88+ [wts,LebesgueConstant] = get_lmmwts(xpts, npts, nqloc,rlen,phi_wts,h,chebApproxInfo);
89+
90+ subplot((p/2),2,i)
91+ bar(wts);
92+ str = '\Lambda';
93+ title(sprintf('order $%2d$, $%s_{%d} = %1.1f$',npts,str,npts,LebesgueConstant),'Interpreter','latex');
94+end
95+
96+sgtitle(sprintf('LMM weights for grid points\n with separation %1.2e',h))
97+saveas(gcf,'/Users/nalhassanieh/Desktop/stabilityResults/LMM_wts2','epsc');
98+end
wfp_1Dspringscattering/get_solnAndErr.madded+63−0View file
@@ -0,0 +1,63 @@
1+function [u,ue,tn_sol,n_sol,tSOL_now,sn_sol] = get_solnAndErr(u,ue,tn_sol,n_sol,...
2+ tSOL_now,x,t,snIdxNow,sn,an,interpNodesIdx_sol,...
3+ interpAndGLWeights_sol,interpShift_sol,s,M,delta,tol,tSOL,solnType,dataParam,fixSolCnt,n,tgh,sn_sol)
4+% GET_SOLNANDERR to compute the solution and error (if available)
5+%
6+% INPUTs
7+% u,ue: computed and exact solutions
8+% tn_sol,n_sol,tSOL_now: solution time grid and index, tSOL_now is a
9+% variable to track the time of the solution with respect to the actual
10+% time grid
11+% x: spatial grid for testing
12+% t,tn: current time, and time grid
13+% snIdxNow: index of sn at the first time-step
14+% sn: the density grid functions
15+% an: Fourier Coefficients
16+% interpNodesIdx_sol,interpAndGLWeights_sol,interpShift_sol:
17+% interpolation parameters needed for the computation of the local part
18+% of the solution.
19+% s,M: source location and total num of sources
20+% tol: error tolerance
21+% mu,t0: parameters needed for the computation of the manufactured
22+% solution for error computation. These variables are held in the
23+% dataParam struct for ease of passing in functions and switching between
24+% test solutions.
25+% W,delta: width of the window function: W*dt = delta
26+% h_temp: paramter needed for free-space projection
27+% tSOL: time-step for the test solution grid
28+% addProj: button to get free-space solution
29+% solnType: manufactured or true
30+% writerObj: object to save movie
31+% plotMovie: button to plot movie
32+% saveMovie: button to save movie
33+%
34+% OUTPUTs
35+% u,ue: computed and exact solutions
36+% tn_sol,n_sol,tSOL_now: solution time grid and index, tSOL_now is a
37+% variable to track the time of the solution with respect to the actual
38+% time grid
39+
40+
41+mu = dataParam.mu;
42+t0 = dataParam.t0;
43+
44+% if(t>=tSOL_now)
45+if(mod(n-tgh,fixSolCnt) == 0 && t>=0)
46+
47+ sn_sol(:,(n_sol + 1)) = sn(end,:)';
48+
49+ u(:,(n_sol+1)) = get_uSol(x,t,snIdxNow,sn,an,interpNodesIdx_sol,...
50+ interpAndGLWeights_sol,interpShift_sol,s,M,delta,tol);
51+
52+ if(strcmp(solnType,'ms'))
53+ ue(:,(n_sol+1)) = get_ExactSol(x,t,mu,t0,M,s);
54+ else
55+ ue(:,(n_sol+1)) = get_uIncident(x,t,mu,t0);
56+ end
57+
58+ tn_sol(n_sol+1) = t;
59+ n_sol = n_sol + 1;
60+ tSOL_now = tSOL_now + tSOL;
61+end
62+
63+end
wfp_1Dspringscattering/get_sourceGrid.madded+16−0View file
@@ -0,0 +1,16 @@
1+function s = get_sourceGrid(src_dmn,ds,M,uniform_sgrid)
2+
3+as = src_dmn(1); bs = src_dmn(2);
4+
5+if(uniform_sgrid == 1)
6+ Ns = M;
7+ s_grid = linspace(as,bs,Ns);
8+ s = s_grid';
9+else
10+ Ns = ceil((bs - as)/ds);
11+ s_grid = linspace(as,bs,Ns); %s = s_grid';
12+ s_idx = randperm(Ns,M);
13+ s = s_grid(s_idx); s = s';
14+end
15+
16+end
wfp_1Dspringscattering/get_springConstants.madded+9−0View file
@@ -0,0 +1,9 @@
1+function beta = get_springConstants(betaMax,M,uniform_beta)
2+
3+if(uniform_beta == 1)
4+ beta = betaMax*ones(M,1);
5+else
6+ beta = 0.1 + (betaMax-0.1)*rand(M,1); % constant for the spring
7+end
8+
9+end
wfp_1Dspringscattering/get_uHistory.madded+28−0View file
@@ -0,0 +1,28 @@
1+function uh = get_uHistory(an,x,tol)
2+% GET_UHISTORY computes the history part of the solution
3+%
4+% uh = get_uHistory(an,x,K)
5+% computes the history part of the numerical solution, given Fourier
6+% modes 'an', 'x' grid points
7+
8+%%% Using Finufft (NU to NU)
9+% uh = (1/(2*pi))*finufft1d3(K,an,-1,tol,x);
10+
11+%%% Using Finufft (U to NU)
12+opts.modeord = 1;
13+uh = (1/(2*pi))*finufft1d2(x,-1,tol,an,opts);
14+
15+% Did not use fft because number of test points is small. We are using here
16+% N frequencies to compute Nx values where Nx<N. The fft operator truncates
17+% at Nx frequencies, which produces incorrect results.
18+
19+%%% Standard if missing Finufft
20+% uh = zeros(size(x));
21+% cnt = 1;
22+% for k = K
23+% uh = uh + an(cnt)*exp(-1i*k*x);
24+% cnt = cnt + 1;
25+% end
26+% uh = (1/(2*pi))*uh;
27+
28+end
wfp_1Dspringscattering/get_uIncident.madded+24−0View file
@@ -0,0 +1,24 @@
1+function uin = get_uIncident(x,t,mu,t0)
2+% GET_UINCIDENT returns an incident wave
3+%
4+% Inputs:
5+% x,t: space and time variables
6+% mu,t0: for the Gaussian as seen in formula below
7+
8+if nargin == 0, test_get_uIncident; return; end
9+
10+uin = exp(-mu.*(x - t - t0).^2);
11+
12+end
13+
14+function test_get_uIncident
15+ax = -pi; bx = pi; Nx = 100;
16+x = linspace(ax,bx,Nx);
17+t = 3;
18+mu = 30;
19+t0 = -3;
20+uin = get_uIncident(x,t,mu,t0);
21+plot(x,uin,'LineWidth',2);
22+xlabel x; ylabel u;
23+title 'plot of u_in';
24+end
wfp_1Dspringscattering/get_uIncidentInfo.madded+28−0View file
@@ -0,0 +1,28 @@
1+function [h0,dataParam] = get_uIncidentInfo(mu)
2+% GET_UINCIDENTINFO prepares information needed for an incident wave of the
3+% form uin = exp(-mu*(x - t - t0).^2);
4+%
5+% INPUT:
6+% tFinal: final time
7+% tol: error tolerance for the full problem (to compute starting dt)
8+% gam: pad gam*Nyquist for window function (to compute starting dt)
9+%
10+% OUTPUT:
11+% mu,t0 as indicated above
12+% Nt0: number of initial time steps to resolve incident pulse
13+% h0: minimum step size to resolve incident pulse
14+% dataParam: is struct to hold mu and t0 (for ease of switching between
15+% types of test solutions
16+
17+t0 = -3;
18+
19+% Choose h0 such that the densities are resolved
20+N0 = 30; % number of grid points per standard deviation (sd)
21+bumpWidth = sqrt(log(1/eps)./mu);
22+h0_vec = bumpWidth./N0;
23+h0 = min(h0_vec);
24+
25+dataParam.mu = mu;
26+dataParam.t0 = t0;
27+
28+end
wfp_1Dspringscattering/get_uLocal.madded+50−0View file
@@ -0,0 +1,50 @@
1+function ul = get_uLocal(x,t,snIdxNow,sn,interpNodesIdx_sol,interpAndGLWeights_sol,interpShift_sol,s,M,delta)
2+% GET_ULOCAL computes the local part of the solution
3+%
4+% ul = get_uLocal(x,t,snIdxNow,sn,interpNodesIdx_sol,...
5+% interpAndGLWeights_sol,interpShift_sol,s,M,delta)
6+% returns the local part of the solution at each value of the vector 'x' at
7+% time 't'. The function takes the unifrom time-grid 'tn', 'delta' ([0,delta]
8+% is the compact support of the window), 'M' source locations in the vector
9+% s = [s(1),...,s(M)], the computed density values sn.
10+% snIdxNow: index of the density at the current time step.
11+% interpNodesIdx_sol: interpolation nodes used to evaluate the density at
12+% GL nodes
13+% interpAndGLWeights_sol: interpolation/integration weights needed to
14+% evaluate the local integral
15+% interpShift_sol: initial shift of interpolation nodes (to undo the shift)
16+
17+% Pick interpolation points between [startTime,endTime], where startTime
18+% and endTime lie on the time-grid
19+Nx = length(x);
20+ ul = zeros(Nx,1);
21+ for i = 1:Nx % for each x value
22+ for j = 1:M % for each source
23+ L = abs(x(i) - s(j));
24+ if(L<delta && (t - L)>0)
25+ snj = sn(:,j);
26+ % perform the local integration using GL and barycentric
27+ % interpolation
28+ shift = (snIdxNow - 1) + interpShift_sol;
29+
30+ shiftedInterpNodesIdx = interpNodesIdx_sol{i,j} + shift;
31+ snIdx = shiftedInterpNodesIdx + 1; %snIdx(snIdx<=0) = 1;
32+ snvec = sum(interpAndGLWeights_sol{i,j}.*snj(snIdx),2);
33+
34+ I = sum(snvec);
35+ ul(i) = ul(i) + I;
36+ end
37+ end
38+ ul(i) = 0.5*ul(i);
39+ end
40+
41+end
42+
43+% Note replaced for loops by matvecs. Original:
44+% snvec = zeros(P,1);
45+% for k = 1:P % for each GL node calculate interpolated value at GL node
46+% shiftedInterpNodesIdx = interpNodesIdx_sol{i,j}(k,:) + shift;
47+% snIdx = shiftedInterpNodesIdx + 1; snIdx(snIdx<=0) = 1;
48+% snvec(k) = interpAndGLWeights_sol{i,j}(k,:)*(sn{j}(snIdx));
49+% end
50+% I = sum(snvec);
wfp_1Dspringscattering/get_uSol.madded+30−0View file
@@ -0,0 +1,30 @@
1+function u = get_uSol(x,t,snIdxNow,sn,an,interpNodesIdx_sol,interpAndGLWeights_sol,interpShift_sol,s,M,delta,tol)
2+% GET_USOL gets the solution by forming local and history parts
3+%
4+% u = get_uSol(x,t,snIdxNow,sn,an,K,interpNodesIdx_sol,...
5+% interpAndGLWeights_sol,interpShift_sol,s,M,delta)
6+% returns the solution u by evaluating the local and the history parts.
7+% The function requires the spatial grid 'x', the current time 't', the
8+% unifrom time grid 'tn', the Fourier coefficients 'an', the frequency
9+% grid 'K', 'delta' the width of the window function used in the
10+% local/history split, the location of sources 's' s = [s(1),...,s(M)]
11+% where 'M' is the number of sources, 'sn' is the computed density.
12+%
13+% Other input values:
14+% snIdxNow: index of the density at the current time step.
15+% interpNodesIdx_sol: interpolation nodes used to evaluate the density at
16+% GL nodes
17+% interpAndGLWeights_sol: interpolation/integration weights needed to
18+% evaluate the local integral
19+% interpShift_sol: initial shift of interpolation nodes (to undo the shift)
20+
21+% get the local part of the solution
22+ul = get_uLocal(x,t,snIdxNow,sn,interpNodesIdx_sol,interpAndGLWeights_sol,interpShift_sol,s,M,delta);
23+
24+% get the history part of the solution
25+uh = get_uHistory(an,x,tol);
26+
27+% compute the solution
28+u = uh + ul;
29+
30+end
wfp_1Dspringscattering/get_xgrid.madded+23−0View file
@@ -0,0 +1,23 @@
1+function x = get_xgrid(ax,bx,src_dmn,Nx,h0,W)
2+% GET_XGRID a function to get the spatial grid for plotting
3+%
4+% INPUTS:
5+% [ax,bx]: spatial domain
6+% [as,bs]: interval where sources are placed
7+% Nx: number of grid points in the spatial discretization
8+% h0: intial time-step size
9+% W: width of the window function
10+% addProj: button to turn on free-space projection
11+% getSelfConvergence: button to turn on self-convergence study
12+
13+as = src_dmn(1); bs = src_dmn(2);
14+
15+h_temp = h0;
16+
17+xmin = (ax + 2*h_temp*W); xmax = (bx - 2*h_temp*W);
18+x = linspace(xmin,xmax,Nx)';
19+if(xmin>as || xmax<bs)
20+ error('Sources are outside allowed region; reduce the time-step or narrow down the domain where sources live');
21+end
22+
23+end
wfp_1Dspringscattering/main.madded+231−0View file
@@ -0,0 +1,231 @@
1+startMatlabFile
2+
3+%% Set file directories to save data, figures, tables and movies
4+addpath ./utils;
5+dataFile = './data';
6+figFile = './fig';
7+
8+%% Testing Buttons
9+plotOpt = 0; % '1' to plot sol at final time, '2' heat maps
10+saveWorkspaceOpt = 0;
11+addErrResults = 0;
12+savePlot = 0;
13+logScale = 0; % '1' to use log scale in heat maps
14+printTimeStep = 0; % '1' to print out each time step
15+evalSol = 1;
16+
17+%% Problem parameters
18+%%% order and test solution
19+order = 2; % order of accuracy = interpolation
20+numResolutions = 3; % num of grid resolutions
21+solnType = 'ms'; % type of solution 'ms' or 'true'
22+
23+%%% domain and final time
24+ax = -pi; bx = pi; % Domain
25+tFinal = 3*pi; % final time
26+tsStop = -1; % stop time stepping at tsStop. ('-1' continue to tFinal)
27+
28+%%% Choose number of spatial and temporal set for solution computation
29+Nx = 10; % size of spatial grid for testing and plotting
30+Nt_sol = 10; % size of the time grid for testing and plotting
31+tSOL = tFinal/Nt_sol; % time-step for the solution evaluation
32+
33+%%% sources
34+M = 10; % number of sources
35+maxNumNeighbors = 10; % set max number of neighbors
36+src_dmn = [-1,1];
37+ds = 1e-4; % set min distance between sources
38+uniform_sgrid = 0;
39+s = get_sourceGrid(src_dmn,ds,M,uniform_sgrid);
40+
41+%%% spring constants
42+betaMax = 3; uniform_beta = 0;
43+beta = get_springConstants(betaMax,M,uniform_beta);
44+
45+%%% incident pulse parameter
46+mu = 30;
47+
48+%%% window
49+tol = 1e-12; % error tolerance
50+[winData,W] = setup_generalwindow(tol);
51+
52+%% Compute some values for the window function and GL
53+P = W; % GL nodes
54+gl = glwt_prep(P); % GL nodes and weights on [-1,1]
55+
56+%% Manufactured density values for testing or incident wave for true scattering
57+if(strcmp(solnType,'ms'));[h0,dataParam] = manufacturedSolution(M,tFinal);
58+else; [h0,dataParam] = get_uIncidentInfo(mu); end
59+
60+% fix initial time step
61+[h0,typNumOfNeighbors] = fixDtBasedOnTypNeighborsNum(maxNumNeighbors,src_dmn,W,M,h0);
62+Nt0 = ceil(tFinal/h0); if(mod(Nt0,2) == 1); Nt0 = Nt0 + 1; end
63+
64+% get the spatial grid
65+x = get_xgrid(ax,bx,src_dmn,Nx,h0,W);
66+
67+%% Print title
68+printTitle(tFinal, tol, W, P, M, typNumOfNeighbors, order, solnType);
69+prechar = sprintf('%s%d_%d_mn%d_',solnType,betaMax,dataParam.mu(1),maxNumNeighbors); % char/string to precede title of fig or tab
70+
71+%% Grid resolution study
72+if(evalSol == 1)
73+ err = zeros(1,numResolutions);
74+ h = zeros(1,numResolutions);
75+ ug = cell(1,numResolutions);
76+ t_cell = cell(1,numResolutions);
77+ x_cell = cell(1,numResolutions);
78+ sn_cell = cell(1,numResolutions);
79+end
80+
81+for m = 1:numResolutions
82+ N = Nt0*(2^(m-1));
83+
84+ dt = 2*pi/N;
85+ Nt = ceil(tFinal/dt);
86+ dt = tFinal/Nt; % timestep
87+ h_RBC = 2*pi/N;
88+
89+ tgh = W + order;
90+ tInitial = -tgh*dt; % initial time
91+ Ntg = Nt + 1 + tgh; % total number of time-steps
92+ tn = linspace(tInitial ,tFinal, Ntg)'; % time grid
93+ it1 = tgh + 1; it2 = Ntg; % interior time indices (start from zero)
94+ It = it1:it2; % time indices starting from zero
95+ delta = dt*W;
96+
97+ %%% time parameter for self convergence study
98+ if(m == 1); fixSolCnt = ceil((Nt + 1)./Nt_sol);
99+ else; fixSolCnt = 2*fixSolCnt; end
100+ Nt_sol = ceil((Ntg-tgh)/fixSolCnt);
101+
102+ % frequency paramenters
103+ K0 = N/2; % max wave number
104+ K = [0:(K0-1),(-K0:-1)]; % frequency set
105+ zeroLoc = 1; % where k = 0 is located
106+
107+ % set up tRBC
108+ tRBC = W*h_RBC; % choose tFinal to be a multiple of pi
109+
110+ %% Allocate space for time-stepping
111+ % let sn hold the density function values at any time tn
112+ timeLevels = tgh + 2;
113+ sn = zeros(timeLevels,M);% store tgh + current time + next time
114+ snIdxNow = tgh+1; % index of sn at the current time
115+
116+ % Let an hold the values of the Fourier coefficients, and bn
117+ % correspond to time derivative of an, needed in the history treatment
118+ an = zeros(N,1);
119+ bn = an; anp1 = an; bnp1 = bn;
120+
121+ % Allocate space for computed solution for heat maps
122+ u = zeros(Nx,Nt_sol); % u(x,t) to generate heat map
123+ ue = zeros(Nx,Nt_sol);
124+ tn_sol = zeros(1,Nt_sol);
125+ sn_sol = zeros(M,Nt_sol);
126+
127+ %% Prepare for time-stepping
128+ t = dt; % this represents the current time
129+
130+ tic
131+ % prepare nodes and weights needed for the evaluation of self local
132+ % integrals and other local integrals
133+ [implicitMat,A_sp,deltaInteractions,dtInteractions] = ...
134+ prepNodesAndWeights(tn(It),dt,beta,gl,s,src_dmn,snIdxNow,P,M,W,order,...
135+ timeLevels,typNumOfNeighbors,winData);
136+
137+ % prepare the nodes and weights for the solution local evaluation
138+ if(evalSol == 1)
139+ [interpNodesIdx_sol,interpAndGLWeights_sol,interpShift_sol] = prepForLocalSolEval(x,dt,gl,s,M,P,W,order,winData);
140+ end
141+
142+ % prepare integrals needed in the evaluation of Fourier coefficients
143+ [p,q,p0,q0,sn_hat,phi_tn,phit_tn] = prepValuesForFourierCoefs(sn,snIdxNow,dt,K,N,W,gl,s,tol,winData);
144+
145+ tm.prep = toc;
146+ %% Time Stepping
147+ tRBC_now = 0; tSOL_now = tSOL; n_sol = 1; totalSteps = 1;
148+ tic
149+ for n = It(2:(end-1)) % note: t(n) = tInitial + (n-1)*dt = (-tgh + (n-1))*dt
150+
151+ % get the history contribution
152+ [anp1,bnp1,historySum,tRBC_now] = get_historyContribution(t,an,bn,sn,zeroLoc,K,...
153+ dt, tRBC, W, snIdxNow, p,q,p0,q0,sn_hat,tRBC_now,phi_tn,phit_tn,s,tol);
154+
155+ % Get the gn for the spring scattering conditions
156+ gn = get_g(t,dt,s,dataParam,M,beta,solnType);
157+
158+ % Perform a sparse matvec to evaluate integrals
159+ sn_v = reshape(sn,[],1);
160+ RHS = gn + (beta/2).*((A_sp*sn_v)+ (1/pi)*historySum);
161+
162+ % solve the implicit system
163+ sn(end,:) = implicitMat\RHS;
164+
165+ % Take the real part of the density
166+ sn = real(sn);
167+
168+ % new time
169+ t = tn(n+1);
170+
171+ if(evalSol == 1)
172+ [u,ue,tn_sol,n_sol,tSOL_now,sn_sol] = get_solnAndErr(u,ue,tn_sol,n_sol,...
173+ tSOL_now,x,t,snIdxNow,sn,anp1,interpNodesIdx_sol,...
174+ interpAndGLWeights_sol,interpShift_sol,s,M,delta,tol,...
175+ tSOL,solnType,dataParam,fixSolCnt,n,tgh,sn_sol);
176+ end
177+
178+ % update values for the next time step
179+ [an,bn,sn,sn_hat] = prepForNextStep(anp1,bnp1,sn,sn_hat,snIdxNow,tgh,s,N,tol);
180+
181+ if(printTimeStep == 1)
182+ fprintf('t = %1.4e\n',t);
183+ end
184+
185+ if(totalSteps == tsStop)
186+ break;
187+ end
188+ totalSteps = totalSteps + 1;
189+ end
190+
191+ timeSteppingTime = toc;
192+ tm.ts = timeSteppingTime/totalSteps;
193+
194+ %% Error Analysis
195+ if(evalSol == 1)
196+ u = real(u);
197+ ug{m} = u;
198+ h(m) = dt;
199+ t_cell{m} = tn_sol;
200+ x_cell{m} = x';
201+ sn_cell{m} = sn_sol;
202+
203+ % Check the error at all time
204+ if(strcmp(solnType,'ms'))
205+ err(m) = max(max(abs(u-ue))); % max-norm error
206+ rate = printOutResults(m,t,Nt,K0,dt,err,h);
207+ end
208+
209+ % Plot and save into a file if needed
210+ if(plotOpt == 1 && strcmp(solnType,'ms'))
211+ plotSolnAndError(t,x,u(:,end),ue,figFile,prechar,order,tFinal,savePlot);
212+ pause(.1);
213+ elseif(plotOpt == 2)
214+ plotHeatMap(u,ue,x,tn_sol,s,figFile,order,tFinal,prechar,savePlot,solnType,logScale,dataFile);
215+ pause(1);
216+ end
217+ end
218+end
219+
220+checkMatchingGrids(t_cell,x_cell,numResolutions);
221+
222+if(numResolutions>2)
223+ [sc_err,sc_rate] = selfConvergence(ug,h,numResolutions);
224+elseif(numResolutions==2)
225+ sc_err = maxNorm(ug{2} - ug{1});
226+ fprintf('self convergence error = %1.2e\n',sc_err);
227+end
228+
229+if(saveWorkspaceOpt == 1 || addErrResults == 1)
230+ saveWorkspace;
231+end
wfp_1Dspringscattering/manufacturedSolution.madded+56−0View file
@@ -0,0 +1,56 @@
1+function [h0,dataParam] = manufacturedSolution(M,tFinal)
2+% MANUFACTUREDSOLUTION analytic density functions
3+%
4+% [sig,Nx0,h0,t0,mu] = manufacturedSolution(M,tFinal,s,beta)
5+% returns in sig a cell array of M density function handles
6+%
7+% Input values:
8+% M: number of spring scatterers (sources)
9+% tFinal: final time
10+% tol: error tolerance for the full problem (to compute starting dt)
11+% gam: pad gam*Nyquist for window function (to compute starting dt)
12+%
13+% Outputs:
14+% sig: function handle in terms of time that produces a vector of density
15+% values evaluated at a specific time t:
16+% [sig_1(t),...,sig_M(t)]
17+% Nx0, h0: Number of grid points needed to ensure sig is resolved on the
18+% time grid, as well as the maximum step-size.
19+% t0, mu : (to document)
20+% dataParam: is a struct involving t0,mu and sig.
21+%
22+% Notes:
23+% Take the densities to be Gaussians for testing.
24+
25+override = 0; % turn on '1' to avoid restrictions on chosen test function
26+
27+mu = linspace(40,50,M)';
28+t0 = linspace(1,3,M)';
29+
30+% choose t0 such that sig(0) = 0.01*eps
31+% t0 = sqrt(log(amp/(0.01*eps))./mu);
32+
33+%%% Working resolution using standard deviation
34+N0 = 20; % number of grid points per standard deviation (sd)
35+bumpWidth = sqrt(log(1/eps)./mu);
36+h0_vec = bumpWidth./N0;
37+h0 = min(h0_vec);
38+
39+%%% Resolution computed from Fourier analysis
40+% omega0 = 2*sqrt(mu).*sqrt(log((1/tol)*sqrt(pi./mu)));
41+% h0 = min(pi*(1 - gam)./omega0);
42+% Nt0 = ceil(tFinal/h0); % or Nt0 = ceil(2*pi/h0);
43+
44+% ensure that the forcing function is bandlimited within the computational time domain
45+if(any(2*t0>tFinal)&&(override == 0))
46+ error('Forcing inappropriate for testing. For tFinal=%1.2f\nChoose other mu or tFinal',tFinal);
47+end
48+
49+% Define the density functions to be Gaussians
50+sig = @(t) exp(-mu.*((t - t0).^2));
51+
52+dataParam.sig = sig;
53+dataParam.t0 = t0;
54+dataParam.mu = mu;
55+
56+end
wfp_1Dspringscattering/maxNorm.madded+5−0View file
@@ -0,0 +1,5 @@
1+function Nrm = maxNorm(v)
2+
3+Nrm = max(max(abs(v)));
4+
5+end
wfp_1Dspringscattering/mergeWeights.madded+26−0View file
@@ -0,0 +1,26 @@
1+function [wv_out] = mergeWeights(A,w,Wth)
2+% MERGEWEIGHTS a function to add weights at similar nodes to avoid double
3+% sums
4+%
5+% Inputs:
6+% A: some matrix A with some repeated elements
7+% w: weights corresponding to each element of matrix A
8+% Outputs:
9+% v: the collapsed matrix A with no repeated elements
10+% wv: the sum of weights of repeated elements of A. Each such sum
11+% corresponds to the unique value of the element saved in v.
12+
13+[m,n] = size(A);
14+v = unique(A); % create a vector with unique A elements
15+wv = zeros(length(v),1);
16+for i = 1:m
17+ for j = 1:n
18+ idx = find(v == A(i,j));
19+ wv(idx) = wv(idx) + w(i,j);
20+ end
21+end
22+
23+wv_out = zeros(Wth,1);
24+wv_out(v(1):(v(1) + length(v)-1)) = wv;
25+
26+end
wfp_1Dspringscattering/plotHeatMap.madded+75−0View file
@@ -0,0 +1,75 @@
1+function plotHeatMap(u,ue,x,tn,s,figFile,order,tFinal,prechar,savePlot,solnType,logScale,dataFile)
2+% PLOTHEATMAP plots a heatmap of a given solution u(x,t)
3+%
4+% plotHeatMap(u,ue,x,tn,figFile,order,tFinal,prechar,savePlot)
5+% returns a heatmap of u with respect to x and t given discrete 'x' grid
6+% and 'tn' time grid. 'figFile' is the file where the figure will be
7+% saved, 'order' and 'tFinal' are needed to record the order of accuracy
8+% and the final time in the file name. 'prechar' is a character used to
9+% precede the file name if needed for special cases. 'savePlot' = 1 to
10+% save plot in the given figFile. 'solnType' type of testing solution used
11+% ms (manufactured solution), true (true scattering example).
12+
13+M = length(s);
14+
15+colormap(jet(256));
16+
17+figure(1);
18+if(strcmp(solnType,'ms'))
19+ if(logScale == 1)
20+ U = log10(abs(u))';
21+ climInt = [-12 0];
22+ solTitle = 'log_{10}|u(x,t)|';
23+ else
24+ U = abs(u)';
25+ climInt = [0 0.5];
26+ solTitle = '|u(x,t)|';
27+ end
28+ imagesc(x,tn,U); hold on;
29+ plot(s,0,'|r','MarkerSize',10); hold off;
30+ axis xy; clim(climInt);
31+ xlabel x; ylabel t; colorbar; title(sprintf('M = %d, order = %d, $%s$',M,order,solTitle));
32+else % in this case ue = u_incident
33+ if(logScale == 1)
34+ U = log10(abs(u + ue))';
35+ climInt = [-12 0];
36+ solTitle = 'log_{10}|u_{tot}(x,t)|';
37+ else
38+ % U = abs(u + ue)';
39+ % climInt = [0 0.05];
40+ % solTitle = '|u_{tot}(x,t)|';
41+
42+ U = (u + ue)';
43+ climInt = [-0.005 0.005];
44+ solTitle = 'u_{tot}(x,t)';
45+ end
46+ imagesc(x,tn,U); hold on;
47+ plot(s,0,'|r','MarkerSize',10);
48+ % xline(s,'--r');
49+ hold off;
50+ axis xy; clim(climInt);
51+ xlabel x; ylabel t; colorbar; title(sprintf('M = %d, order = %d, $%s$',M,order,solTitle));
52+end
53+
54+if(savePlot == 1)
55+ fileName = sprintf('%s/%sDSheatmap%d_%d_%d',figFile,prechar,order,round(tFinal),M);
56+ print(fileName,'-depsc');
57+
58+ fileName2 = sprintf('%s/%sDSheatmapDATA%d_%d_%d',dataFile,prechar,order,round(tFinal),M);
59+ save(fileName2);
60+end
61+
62+% plot the heatmap of the error
63+if(strcmp(solnType,'ms'))
64+ figure(2);
65+ imagesc(x,tn,log10(abs(u-ue))'); axis xy; clim([-12 0]); hold on;
66+ plot(s,0,'|r','MarkerSize',10); hold off;
67+ xlabel x; ylabel t; colorbar; title(sprintf('M = %d, order = %d, log_{10}|u(x,t)-u_e(x,t)|',M,order));
68+
69+ if(savePlot == 1)
70+ fileName = sprintf('%s/%sDSheatmap%d_%d_%d',figFile,prechar,order,round(tFinal),M);
71+ print(fileName,'-depsc');
72+ end
73+end
74+
75+end
wfp_1Dspringscattering/plotSolnAndError.madded+35−0View file
@@ -0,0 +1,35 @@
1+function plotSolnAndError(t,x,un,ue,figFile,prechar,order,tFinal,savePlot)
2+% PLOTSOLNANDERROR plots the solution and error
3+%
4+% plotSolnAndError(t,tn,sn,src,sige,figFile,prechar,order,tFinal,plotOption)
5+% Inputs:
6+% tn: uniform time grid
7+% sn: density grid function
8+% sige: exact density function
9+% figFile: file to store figures
10+% prechar: character to precede name of figure
11+
12+figure(1);
13+plot(x,un,'-o','LineWidth',2); hold on
14+plot(x,ue,'-x','LineWidth',2);hold off
15+title(sprintf('Computed vs exact at t = %1.1f',t));
16+axis tight;
17+legend('u','u_e','Location','best');
18+xlabel('x'); ylabel('y');
19+set(gca,'fontSize',16);
20+if(savePlot == 1)
21+ saveas(gcf,sprintf('%s/%sDScomputed%d_%d',figFile,prechar,order,round(tFinal)),'png');
22+end
23+
24+figure(2);
25+plot(x,un-ue,'LineWidth',2);
26+title(sprintf('Error t = %1.1f',t));
27+axis tight;
28+legend('error','Location','best');
29+xlabel('t'); ylabel('y');
30+set(gca,'fontSize',16);
31+if(savePlot == 1)
32+ saveas(gcf,sprintf('%s/%sDSerror%d_%d',figFile,prechar,order,round(tFinal)),'png');
33+end
34+
35+end
wfp_1Dspringscattering/prepForLocalSolEval.madded+75−0View file
@@ -0,0 +1,75 @@
1+function [interpNodesIdx_sol,interpAndGLWeights_sol,interpShift_sol] = prepForLocalSolEval(x,dt,gl,s,M,P,W,order,winData)
2+% PREPFORLOCALSOLEVAL precomputes parameters for local solution evaluation
3+%
4+% [interpNodesIdx_sol,interpAndGLWeights_sol,interpShift_sol] = ...
5+% prepForLocalSolEval(x,dt,delta,gl,s,M,phi_wts,P,W,order,chebApproxInfo)
6+% prepares interpolation nodes and weights needed for the evaluation
7+% of the local part of the solution. The function takes 'x' the
8+% spatial grid, 'dt' the time-step, 'delta = W*dt' the width of the window
9+% function 'gl' GL nodes and weights on [-1,1], 's' a vector
10+% carrying the location of the sources, 'M' the number of sources, the
11+% blending function cheb weights 'phi_wts' and the number of GL nodes 'P'.
12+%
13+% NOTE:'chebApproxInfo': a struct with fields domain [ax,bx] over which
14+% cheb weights are computed; ninters, number of subintervals, and nord,
15+% order of polynomial interpolation or number of cheb nodes
16+%
17+% OUTPUT:
18+% interpNodesIdx_sol: indices of interpolation nodes
19+% interpAndGLWeights_sol: weights needed to perform the integration of
20+% the local part of the soluiton
21+% intepShift_sol: how many grid points are the interpolation points
22+% shifted. We undo the shift when we use these nodes in the
23+% time-stepping loop.
24+
25+delta = W*dt;
26+
27+Nx = length(x); % number of points in the spatial grid
28+
29+% allocate space for variables
30+interpNodesIdx_sol = cell(Nx,M); interpNodesIdx_sol(:) = {zeros(P,order)};
31+interpAndGLWeights_sol = cell(Nx,M); interpAndGLWeights_sol(:) = {zeros(P,order)};
32+
33+interpShift_sol = -(W-1); % shift index of the interpolation nodes (nodes shifted forward by delta)
34+
35+for i = 1:Nx % for each x vlue
36+ for j = 1:M % for each source
37+ L = abs(x(i) - s(j));
38+ if(L<delta)
39+
40+ % prepare the GL nodes eta needed for window factor
41+ [eta,~] = glwt(L,delta,gl);
42+
43+ % Get the blending function at eta
44+ phiL_to_delta_sol = generalwindow((eta./delta),winData);
45+
46+ % actual interval of integration
47+ % t1 = t - delta;
48+ % t2 = t - L;
49+
50+ % interval of integration (will be shifted later)
51+ t1 = 0;
52+ t2 = delta - L;
53+
54+ % Get the GL nodes and weights in the interval of integration
55+ [glnodes,glwts] = glwt(t1,t2,gl);
56+
57+ % update the glwts by the window factor
58+ glwts = glwts.*(1 - phiL_to_delta_sol(end:-1:1));
59+
60+ % perform Barycentric interpolation at the GL nodes and get
61+ % weights
62+ interpWeights_sol = zeros(P,order);
63+ endTime = t2;
64+ for k = 1:P
65+ [interpNodes,interpNodesIdx_sol{i,j}(k,:)] = get_InterpNodes(glnodes(k), dt, order,endTime);
66+ interpWeights_sol(k,:) = barycentricInterp(interpNodes, order, glnodes(k));
67+ end
68+
69+ % Multiply the interpolation weights by the GL weights
70+ interpAndGLWeights_sol{i,j} = interpWeights_sol.*glwts;
71+ end
72+ end
73+end
74+
75+end
wfp_1Dspringscattering/prepForNextStep.madded+27−0View file
@@ -0,0 +1,27 @@
1+function [an,bn,sn,sn_hat] = prepForNextStep(anp1,bnp1,sn,sn_hat,snIdxNow,tgh,s,N,tol)
2+% PREPFORNEXTSTEP prepares for the next step in the the time-stepping
3+%
4+% INPUTs:
5+% anp1,bnp1: Fourier coefficients at new step and time derivative
6+% sn, sn_hat: density grid function and transfrom
7+% snIdxNow: index of density at the first time step
8+% tgh: number of ghost points in time
9+% s: location of sources
10+% N: number of Fourier coefficients
11+% tol: error tolerance
12+%
13+% OUTPUTs:
14+% an,bn: Fourier coefficients at the old step
15+% sn,sn_hat: density grid function and transfrom
16+
17+% update values for the next time step
18+an = anp1; bn = bnp1;
19+
20+% Update the density values for the next time-step
21+for kk = 1:(tgh + 1)
22+ sn(kk,:) = sn((kk+1),:);
23+end
24+sn_hat1 = finufft1d1(s,sn(snIdxNow,:)',1,tol,N,struct('modeord',1));
25+sn_hat = [sn_hat1,sn_hat(:,1:(end-1))];
26+
27+end
wfp_1Dspringscattering/prepNodesAndWeights.madded+182−0View file
@@ -0,0 +1,182 @@
1+function [implicitMat,A_sp,deltaInteractions,dtInteractions] = ...
2+ prepNodesAndWeights(tn,dt,beta,gl,s,src_dmn,snIdxNow,P,M,W,order,...
3+ timeLevels,typNumOfNeighbors,winData)
4+% PREPNODESANDWEIGHTS to prepare nodes and weights for the numerical
5+% evaluation of integrals
6+%
7+% INPUTS:
8+% tn,dt: time grid (array) with time step dt
9+% phi_wts: array of weights to evaluate the blending function using
10+% Chebyshev polynomial
11+% beta: an array containing spring constants
12+% gl: GL nodes gl.x0 and weights gl.w0 over the interval [-1,1];
13+% s: a vector holding the locations of the spring sources
14+% snIdxNow: index of density grid function at the current step. For
15+% example, the density functions are arranged in the following order
16+% -tgh, ..., -1,0,1, where tgh corresponds to temporal ghost points,
17+% snIdxNow corresponds to the index of 0.
18+% P: number of GL nodes
19+% M: number of sources
20+% W: width of the window function (i.e. support of the window function up
21+% to tol)
22+% chebApproxInfo: information needed for the Chebyshev polynomial
23+% approximation of window and blending functions
24+% timeLevels: the number of time levels stored for the density functions;
25+% this is equal to W + order to account for centered barycentric
26+% interpolation at the left end of the time axis
27+% typNumOfNeighbors: estimated average typical number of neighbors for
28+% spring sources
29+%
30+% OUTPUTS:
31+% implicitMat: implicit sparse matrix to solve for the density values at
32+% the new time
33+% A_sp: a sparse matrix to evaluate integrals numerically using a
34+% combination of barycentric interpolation and Gauss-Legendre quadrature
35+
36+% domain of sources
37+as = src_dmn(1); bs = src_dmn(2);
38+
39+% support of the window function
40+delta = W*dt;
41+
42+% prepare the implicit matrix
43+implicitMat = speye(M);
44+
45+% get the LMM weights for the treatment of selfLocal1
46+selfLocal1_wts = get_lmmwts(tn(1:order), order, order,dt,dt,delta,winData,0);
47+coef = -(1 + (selfLocal1_wts(end).*(beta./2))); % coefficient of the sigma at (n+1)
48+implicitMat = coef.*implicitMat;
49+
50+% GL nodes and weights for selfLocal2 and blending function on [dt,delta]
51+[eta,~] = glwt(dt,delta,gl);
52+phidt_to_delta = generalwindow((eta./delta),winData);
53+
54+% GL nodes and weights for the selfLocal2 implementation and
55+% interpolation weights
56+[selfLocal2_glnodes,selfLocal2_glwts] = glwt(0,(delta-dt),gl);
57+selfLocal2_interpWeights_temp = zeros(P,order);
58+interpNodesIdx_temp = zeros(P,order);
59+for kk = 1:P % for each GL node calculate interpolated value at GL node
60+ [interpNodes,interpNodesIdx_temp(kk,:)] = get_InterpNodes(selfLocal2_glnodes(kk), dt, order,(delta-dt));
61+ selfLocal2_interpWeights_temp(kk,:) = barycentricInterp(interpNodes, order, selfLocal2_glnodes(kk));
62+end
63+% multiply the weights by GL weights and window factor
64+selfLocal2_interpWeights_temp = selfLocal2_interpWeights_temp.*(1 - phidt_to_delta(end:-1:1)).*selfLocal2_glwts;
65+
66+% shift the nodes to fit the position of each node
67+interpNodesIdx_temp = interpNodesIdx_temp+ snIdxNow + 1 - W;
68+
69+% merge the weights of similar interpolation nodes
70+[selfLocal2_interpWeights] ...
71+ = mergeWeights(interpNodesIdx_temp,selfLocal2_interpWeights_temp,timeLevels);
72+
73+% create bins to avoid checks over the full set of sources
74+% each box is delta big, so check sources in the current box and other
75+% surrounding boxes
76+nboxes = ceil((bs - as)/delta); % number of boxes
77+snew = (s - as)./(bs - as); % transform s so it is on [0,1] for assign function
78+[ioffst, ibox, isradr,icnt] = assign(nboxes, snew, M); % assign sources to boxes
79+
80+% Fill in a sparse matrix with weights to evaluate within the time loop
81+% The product of the sparse matrix with the density vector evaluates
82+% numerical integrals needed in the solver
83+nz = M*timeLevels*typNumOfNeighbors;
84+i_sp = ones(1,nz);
85+j_sp = ones(1,nz);
86+a_sp = zeros(1,nz);
87+
88+cnt = 1;
89+dtInteractions = zeros(M,1); % set up counts for delta interactions for testing
90+deltaInteractions = zeros(M,1); % set up counts for dt interactions for testing
91+for j = 1:M
92+ adr = ibox(j);
93+ sourcesNearby_left = []; sourcesNearby_right = [];
94+ sourcesNearby_center = isradr(ioffst(adr):(ioffst(adr) + icnt(adr) - 1));
95+ if(adr>1)
96+ sourcesNearby_left = [isradr(ioffst(adr-1):(ioffst(adr-1) + icnt(adr-1) - 1))];
97+ end
98+ if(adr<nboxes)
99+ sourcesNearby_right = [isradr(ioffst(adr+1):(ioffst(adr+1) + icnt(adr+1) - 1))];
100+ end
101+ sourcesNearby = [sourcesNearby_center; sourcesNearby_left; sourcesNearby_right];
102+ for l = sourcesNearby'
103+ % prepare the indices of the sparse matrix
104+ lInd = ((l-1)*timeLevels + 1):((l-1)*timeLevels + timeLevels);
105+ Idx = ((cnt - 1)*timeLevels+1):((cnt - 1)*timeLevels+timeLevels);
106+
107+ if(l~=j)
108+ % local evaluation from other sources
109+ L = abs(s(j) - s(l));
110+ if(L<delta && L>dt)
111+ % get GL nodes and weights for the integral involving other
112+ % neighboring sources, but does not include the density at
113+ % the new step
114+ [otherLocal2_glnodes,otherLocal2_glwts] = glwt(0,(delta - L),gl);
115+ otherLocal2_interpWeights_temp = zeros(P,order);
116+ otherInterpNodesIdx_temp = zeros(P,order);
117+ for kk = 1:P % for each GL node calculate interpolated value at GL node
118+ [otherInterpNodes,otherInterpNodesIdx_temp(kk,:)] = get_InterpNodes(otherLocal2_glnodes(kk), dt, order,(delta - L));
119+ otherLocal2_interpWeights_temp(kk,:) = barycentricInterp(otherInterpNodes, order, otherLocal2_glnodes(kk));
120+ end
121+ % prepare the values of phi to evaluate the local
122+ % contribution from each of the other sources
123+ [eta,~] = glwt(L,delta,gl);
124+ phiL_to_delta = generalwindow((eta./delta),winData);
125+
126+ % multiply the interpolation weights with the blending
127+ % function
128+ otherLocal2_interpWeights_temp = otherLocal2_interpWeights_temp.*(1 - phiL_to_delta(end:-1:1)).*otherLocal2_glwts;
129+
130+ % shift the indices;
131+ otherInterpNodesIdx_temp = otherInterpNodesIdx_temp + snIdxNow + 1 - W;
132+
133+ [otherLocal2_interpWeights] ...
134+ = mergeWeights(otherInterpNodesIdx_temp,otherLocal2_interpWeights_temp,timeLevels);
135+
136+ i_sp(Idx) = j;
137+ j_sp(Idx) = lInd;
138+ a_sp(Idx) = otherLocal2_interpWeights;
139+
140+ cnt = cnt + 1 ;
141+ deltaInteractions(j) = deltaInteractions(j) + 1;
142+
143+ elseif(L<dt)
144+
145+ otherLocal1_wts = get_lmmwts(tn(1:order), order, order,(dt - L),dt,delta,winData,0);
146+ coef_other = (-beta(j)/2)*otherLocal1_wts(end); % coefficient of the sigma at (n+1)
147+ implicitMat(j,l) = coef_other;
148+
149+ i_sp(Idx) = j;
150+ j_sp(Idx) = lInd;
151+ a_sp(Idx) = selfLocal2_interpWeights;
152+
153+ % Adding the local 1 contribution
154+ endTerm = Idx(end)-1; % Index refering tothe current step
155+ a_sp((endTerm+1-(order-1)):endTerm) = a_sp((endTerm+1-(order-1)):endTerm) + otherLocal1_wts((end-(order-1)):end-1);
156+ cnt = cnt + 1 ;
157+
158+ dtInteractions(j) = dtInteractions(j) + 1;
159+ end
160+ elseif(j == l)
161+ % Adding the self local 2 contributions
162+ i_sp(Idx) = j;
163+ j_sp(Idx) = lInd;
164+ a_sp(Idx) = selfLocal2_interpWeights;
165+
166+ % Adding the self local 1 contributions
167+ endTerm = Idx(end)-1; % Index refering to the current step
168+ a_sp((endTerm+1-(order-1)):endTerm) = a_sp((endTerm+1-(order-1)):endTerm) + selfLocal1_wts((end-(order-1)):end-1);
169+ cnt = cnt + 1 ;
170+ end
171+
172+ end % end of l loop
173+end % end of j loop
174+
175+% decompose the matrix to speed up solution in the time-stepping loop
176+implicitMat = decomposition(implicitMat);
177+
178+A_sp = sparse(i_sp,j_sp,a_sp,M,(M*timeLevels));
179+
180+clear i_sp j_sp a_sp;
181+
182+end % end of function
wfp_1Dspringscattering/prepValuesForFourierCoefs.madded+61−0View file
@@ -0,0 +1,61 @@
1+function [p,q,p0,q0,sn_hat,phi_tn,phit_tn] = prepValuesForFourierCoefs(sn,snIdxNow,dt,K,N,W,gl,s,tol,winData)
2+% PREPVALUESFORFOURIERCOEFS prepares parameters for the evaluation of the
3+% history part of the solution.
4+%
5+% OUTPUT:
6+% the integrals p and q needed for the history evaluation (see notes). p0
7+% and q0 correspond to the zero-frequency case.
8+% sn_hat = \sum_{j = 1}^M \sigma_j(\tau)e^{ikx_j} at initial tau =
9+% (0:W)*dt
10+%
11+% INPUT:
12+% 'dt' timestep, 'K' vector of frequencies, 'W' width
13+% of the window function, 'gl' GL nodes and weights on [-1,1], the
14+% window function and its derivative 'phit' and 'phitt'; 's' is a vector
15+% containing source locations.
16+% 'chebApproxInfo': a struct with fields domain [ax,bx] over which
17+% cheb weights are computed; ninters, number of subintervals, and nord,
18+% order of polynomial interpolation or number of cheb nodes
19+%
20+% NOTE: 'expMat' and 'expMat2'
21+% correspond to matrices involving exponential evaluations needed in
22+% history treatment (needed in earlier version of the code)
23+
24+% Get GL nodes and weights on the interval [0,dt]
25+[tau,w] = glwt(0,dt,gl);
26+
27+delta = W*dt;
28+
29+k = K';
30+tau = tau';
31+lK = length(K);
32+
33+% allocate space for outputs
34+p = zeros(lK,(W+1)); q = zeros(lK,(W+1));
35+p0 = zeros(W+1,1); q0 = p0;
36+
37+for i = 0:W
38+ gam = i*dt;
39+
40+ winArg = (tau + gam)./delta;
41+ [~,phit,phitt] = generalwindow(winArg,winData);
42+ phit = (1/delta)*phit;
43+ phitt = (1/delta^2)*phitt;
44+
45+ I = (2*cos(k*(tau+gam)).*phit+(1./k).*(sin(k*(tau+gam)).*phitt));
46+ p(:,i+1) = (((1./k).*sin(k*(dt - tau))).*I)*w;
47+ q(:,i+1) = ((cos(k*(dt - tau))).*I)*w;
48+
49+ I0 = (2*phit + (tau+gam).*phitt);
50+ p0(i+1) = sum(((dt - tau).*I0)'.*w);
51+ q0(i+1) = sum(I0'.*w);
52+end
53+
54+sn_hat = finufft1d1(s,sn(snIdxNow-(0:W),:)',1,tol,N,struct('modeord',1));
55+
56+%%% free-space parameters
57+phi_grid = (1:W) - 0.5;
58+[phi_tn1,phit_tn1,~] = generalwindow(phi_grid./W,winData);
59+phi_tn = phi_tn1'; phit_tn = (1/W)*phit_tn1';
60+
61+end
wfp_1Dspringscattering/printOutResults.madded+16−0View file
@@ -0,0 +1,16 @@
1+function [rate] = printOutResults(m,t,Nt,K0,dt,err,ht)
2+% PRINTOUTRESULTS prints out error results and rates in the command window
3+%
4+% printOutResults(m,t,Nt,K0,dt,err,L2err,ht)
5+
6+rate = 0;
7+fprintf('t=%8.2e: Nt=%4d K = %3d dt=%9.3e max-err=%8.2e',t,Nt,K0,dt,err(m));
8+if( m>1 )
9+ rate = log(err(m-1)/err(m))/log(ht(m-1)/ht(m));
10+ fprintf(' <strong>inf-rate=%4.2f</strong>',rate);
11+ fprintf('\n');
12+else
13+ fprintf('\n');
14+end
15+
16+end
wfp_1Dspringscattering/printTitle.madded+21−0View file
@@ -0,0 +1,21 @@
1+function printTitle(tFinal, tol, W, P, M, typNumOfNeighbors, order, solnType)
2+% PRINTTITLE print a title for the driver file indicating important
3+% problem parameters
4+%
5+% printTitle(tFinal, tol, W, P, gam, order, src)
6+% Displays important problem paramters in the command window
7+%
8+% Input values:
9+% tFinal: final time
10+% tol: tolerance for the window function
11+% W: W*dt is the support of window function
12+% P: number of GL nodes usually set to W
13+% gam: pad gam*Nyquist for window
14+% order: order of accuracy
15+% M: number of spring scatterers
16+% solnType: type of testing solution
17+
18+fprintf('wave: t = %3.0f, tol = %1.0e, W = %d, P = %d, M = %d, typNumNeighbors = %d, order = %d, solnType = %s\n',...
19+ tFinal, tol, W, P, M, typNumOfNeighbors, order, solnType);
20+
21+end
wfp_1Dspringscattering/saveWorkspace.madded+36−0View file
@@ -0,0 +1,36 @@
1+%% determine the file name
2+if (uniform_sgrid == 1)
3+ str = 'uni'; else; str = 'rand'; end
4+
5+% Save only the selected variables
6+fileName = sprintf('%s/WFPworkspaceO%d_T%d_%sM%d_b%d',dataFile,order,round(tFinal),str,M,betaMax);
7+
8+if (saveWorkspaceOpt == 1)
9+ % all variables in the workspace
10+ vars = who;
11+
12+ % variables to exclude
13+ exclude = {'A_sp','It','RHS','deltaInteractions',...
14+ 'dtInteractions','evalSol','gl','gn','historySum',...
15+ 'implicitMat','interpAndGLWeights_sol','interpNodesIdx_sol',...
16+ 'interpShift_sol','it1','it2','p','p0','phi_tn','phit_tn',...
17+ 'printTimeStep','q','q0','saveWorkspaceOpt',...
18+ 'sn_v','timeSteppingTime','uniform_beta','uniform_sgrid',...
19+ 'winData'};
20+
21+ % Find variables to save
22+ vars_to_save = setdiff(vars, exclude);
23+ save(fileName, vars_to_save{:});
24+ fprintf('saved workspace to %s\n',matlabDataFile);
25+
26+elseif(addErrResults == 1)
27+ matlabDataFile =strcat(fileName,'.mat');
28+ if(isfile(matlabDataFile))
29+ save(fileName, 'sc_err', '-append');
30+ fprintf('added error to the data file in %s\n',matlabDataFile);
31+ else
32+ fprintf('Did not add the error to the data file %s\n',matlabDataFile);
33+ end
34+end
35+
36+
wfp_1Dspringscattering/selfConvergence.madded+58−0View file
@@ -0,0 +1,58 @@
1+function [err,rate] = selfConvergence(ug,h,M)
2+
3+if nargin == 0, test_selfConvergence; return; end
4+
5+if(M<3)
6+ error('need at least 3 resolutions to perform self-convergence study');
7+end
8+
9+m = M-2;
10+
11+Dmp1 = ug{m+1} - ug{m};
12+Dmp2 = ug{m+2} - ug{m+1};
13+
14+NDmp1 = maxNorm(Dmp1);
15+NDmp2 = maxNorm(Dmp2);
16+ratio = NDmp1/NDmp2;
17+r = h(m)/h(m+1);
18+rate = log(ratio)/log(r);
19+
20+C = Dmp2/abs((h(m+1)^(rate)) - (h(m+2)^(rate)));
21+
22+err = zeros(1,M);
23+for k=1:M
24+ uerr = C*(h(k)^(rate));
25+ err(k) = maxNorm(uerr);
26+end
27+
28+fprintf('self convergence rate = %1.2e ', rate);
29+fmt=['errors =' repmat(' %1.2e;',1,numel(err)) '\n'];
30+fprintf(fmt,err);
31+end
32+
33+%%% test function
34+function test_selfConvergence
35+
36+ue = @(x,t) cos(x + 3).*exp(x + 1).*log(1 + x).*sin(2*t - 2).*exp(-t);
37+
38+h0 = 1e-3;
39+h(1) = h0;
40+h(2) = h0/2;
41+h(3) = h0/4;
42+h(4) = h0/8;
43+
44+x = -pi:h0:pi;
45+t = linspace(0,1,length(x));
46+
47+rate = 1.5;
48+numRes = 4;
49+ug = cell(numRes,1);
50+C = (1 + (100-1).*rand(size(x)));
51+for k = 1:numRes
52+ ug{k} = ue(x,t) + C*(h(k)^(rate)) + C*(h(k)^(rate+3));
53+end
54+
55+selfConvergence(ug,h,numRes);
56+
57+end
58+
wfp_1Dspringscattering/startMatlabFile.madded+5−0View file
@@ -0,0 +1,5 @@
1+clear; clf; set(0,'DefaultLineLineWidth',2);
2+set(groot,'defaultAxesFontSize',16);
3+set(groot, 'defaultAxesTickLabelInterpreter','latex');
4+set(groot, 'defaultLegendInterpreter','latex');
5+set(groot, 'defaultTextInterpreter', 'latex');
wfp_1Dspringscattering/utils/.DS_Storeadded+0−0View file
Binary file not shown.
wfp_1Dspringscattering/utils/antiderivmat_1d.madded+37−0View file
@@ -0,0 +1,37 @@
1+function L = antiderivmat_1d(s)
2+% ANTIDERIVMAT_1D matrix from nodes in 1D to antiderivatives at nodes
3+%
4+% L = antiderivmat_1d(s) returns (N-1)*N matrix taking values on s,
5+% a list of N nodes, to their integrals from the first node to each of the
6+% other N-1 nodes in turn. Ie, the antiderivative with constant chosen so that
7+% its value at the first node would be zero. Nodes must be spaced in a sensible
8+% way. length(s) should not exceed around 30 for stability reasons.
9+%
10+% Notes: 1) should match Leslie's get_integratemat().
11+% 2) Computed in Helsing style, with centering and scaling for stability.
12+
13+% Barnett 8/13/24.
14+
15+if nargin==0, test_antiderivmat_1d; return; end
16+
17+cen = (max(s)+min(s))/2; hwid = (max(s)-min(s))/2; % affine map s to [-1,1]
18+s = (s-cen)/hwid;
19+n = numel(s); s = s(:); % col vec
20+V = ones(n); for j=2:n, V(:,j) = V(:,j-1).*s; end % Vandermonde (polyval) mat
21+U = diag(s)*V*diag(1./(1:n)); % mat evaluating an antideriv of poly
22+L = (V'\U')'; % backwards-stable way to solve for it (Helsing)
23+L = L(2:end,:) - L(1,:); % adjust const to zero at first node
24+L = L*hwid; % unscale
25+
26+%%%%%%
27+function test_antiderivmat_1d
28+off = 4.3; sc = 1.7; % test centering and scaling
29+x = sc*linspace(-1,1,16)' + off; % nodes
30+%x = sc*gauss(16) + off; % nodes
31+f = @(x) sin(0.8*x + 0.7); % the antiderivative
32+fp = @(x) 0.8*cos(0.8*x + 0.7); % the input func
33+Fex = f(x(2:end))-f(x(1)); % exact ans at nodes 2...N, col vec
34+L = antiderivmat_1d(x);
35+F = L * fp(x); % hit L against vec of func values
36+fprintf('max abs err for antideriv on nodes : %.3g\n',max(abs(F - Fex)))
37+fprintf('[mat inf-norm = %.3g; max element size = %.3g]\n',norm(L,inf),max(abs(L(:))))
wfp_1Dspringscattering/utils/assign.madded+77−0View file
@@ -0,0 +1,77 @@
1+function [ioffst, ibox, isradr,icnt] = assign(nboxes, xat, natoms)
2+% ASSIGN Assigns sources and targets to boxes.
3+%
4+% INPUTS:
5+% nboxes - Number of boxes
6+% xat - Positions of sources
7+% natoms - Number of sources
8+%
9+% OUTPUTS:
10+% ioffst - Offsets for sources in each box
11+% ibox - Addresses of the boxes containing each target
12+% center - Centers of each box
13+if nargin == 0, test_assign; return; end
14+
15+% xat =(xat - xat(1))./(xat(end) - xat(1));
16+
17+% Initialize variables
18+h = 1 / nboxes;
19+icnt = zeros(nboxes, 1);
20+ibox = zeros(natoms, 1);
21+isradr = zeros(natoms, 1);
22+
23+% Find box in which each source lies and increment counter
24+for j = 1:natoms
25+ ixh = floor(xat(j)/h);
26+ if ixh >= nboxes
27+ ixh = nboxes-1;
28+ elseif ixh <= 0
29+ ixh = 0;
30+ end
31+ iadr = ixh + 1;
32+ icnt(iadr) = icnt(iadr) + 1;
33+ ibox(j) = iadr;
34+end
35+
36+% Compute ioffst array
37+ioffst = zeros(nboxes, 1);
38+ioffst(1) = 1;
39+for j = 2:nboxes
40+ ioffst(j) = ioffst(j-1) + icnt(j-1);
41+end
42+
43+% Reset icnt for recording source locations in the original array
44+icnt(:) = 0;
45+for j = 1:natoms
46+ iadr = ibox(j);
47+ indx = ioffst(iadr) + icnt(iadr);
48+ isradr(indx) = j;
49+ icnt(iadr) = icnt(iadr) + 1;
50+end
51+
52+end
53+
54+function test_assign
55+as = -1; bs = 0;
56+Ns = 10;
57+s = linspace(as,bs,Ns);
58+s_idx = randperm(Ns,Ns);
59+s = s(s_idx);
60+
61+figure(1);
62+plot(s,0,'*b');
63+
64+nboxes = 5;
65+
66+snew = (s - as)./(bs - as);
67+
68+[ioffst, ibox, isradr,icnt] = assign(nboxes, snew, Ns);
69+
70+plot(s,0,'b*',s(ioffst),0,'r|');
71+
72+box = 3
73+start = ioffst(box)
74+isradr(ioffst(box):(ioffst(box) + icnt(box) - 1))
75+
76+pause;
77+end
wfp_1Dspringscattering/utils/barycentricInterp.madded+38−0View file
@@ -0,0 +1,38 @@
1+function interpWeights = barycentricInterp(interpNodes, npts, xval)
2+% BARYCENTRICINTERP returns the barycentric interpolation weights given a
3+% set of interpolation nodes
4+%
5+% interpWeights = barycentricInterp(interpNodes, npts, xval)
6+% returns barycentric interpolation weights given 'interpNodes', a vector
7+% of size 'npts' containing interpolation nodes, and 'xval', a value where
8+% the barycentric interpolating polynomial is evaluated
9+
10+% if nargin == 0, test_barycentricInterp; return; end
11+
12+% Precompute barycentric interpolation factors
13+X=repmat(interpNodes,1,npts);
14+w = 1./prod(-X+X.'+eye(npts),1);
15+
16+tol = 1.0e-20;
17+xvalminusInterpNodes = xval - interpNodes';
18+rr = prod(xvalminusInterpNodes);
19+
20+interpWeights = rr*w./(xvalminusInterpNodes);
21+interpWeights(abs(xval - interpNodes')<tol) = 1;
22+
23+end
24+
25+% function test_barycentricInterp
26+% clf;
27+% ax = 0; bx = 1;
28+% npts = 6;
29+% interpNodes = linspace(ax,bx,npts);
30+%
31+% xval = (bx-ax)*rand + ax;
32+% interpWeights = barycentricInterp(interpNodes, npts, xval);
33+%
34+% bar(interpWeights);
35+% title(sprintf(['Interpolation weights\n for %d equidistant interpolation...' ...
36+% ' nodes\n in [%1.0f,%1.0f] evaluated at x = %1.2f'],npts,ax,bx,xval));
37+%
38+% end
wfp_1Dspringscattering/utils/blending.madded+32−0View file
@@ -0,0 +1,32 @@
1+function phi = blending(phit, tau, tol)
2+% BLENDING retrieves the blending function
3+%
4+% phi = blending(phit, tau, tol) gives the value of blending function, the
5+% antiderivative of the window function 'phit' at a given value 'tau' to
6+% an error tolerance 'tol'
7+
8+if nargin == 0, test_blending; return; end
9+
10+phi = zeros(length(tau),1);
11+for i = 1:length(tau)
12+ % can replace quadgk by integral (not sure which is better)
13+ phi(i) = quadgk(phit,0, tau(i),'RelTol',tol,'AbsTol',tol);
14+end
15+
16+phi = real(phi);
17+end
18+
19+function test_blending
20+tFinal = 1;
21+Nt = 100;
22+tol = 1e-12;
23+b = log(1/tol);
24+gam = 0.5;
25+w = ceil(2*b/(pi*gam));
26+t = linspace(0,tFinal,Nt); dt = t(2) - t(1);
27+phit = window(dt,b,w);
28+phi = blending(phit, t, tol);
29+plot(t,phi, 'LineWidth',2);
30+title('blending function \phi');
31+xlabel('x'); ylabel('y');
32+end
wfp_1Dspringscattering/utils/chebEval.madded+72−0View file
@@ -0,0 +1,72 @@
1+function y = chebEval(x,wts,chebApproxInfo)
2+% CHEBEval retruns the approximation of a function at x
3+%
4+% y = chebApprox(x,ax,h,wts,ninters,nord) approximates a function f
5+% (included in the wts argument) using chebyshev polynomial expansion on
6+% subintervals of a domain that starts at ax
7+%
8+% INPUTS:
9+% x: value in the domain where f is to be approximated
10+% wts: chebyshev weights multiplied by function values in all
11+% subintervals, wts(nord,ninters)
12+% chebApproxInfo.domain: domain [ax,bx] overwhich wts are computed.
13+% chebApproxInfo.ninters: number of subintervals
14+% chebApproxInfo.nord: order of polynomial interpolation; number of cheb
15+% nodes
16+% NOTE: code adjusted from code provided by Leslie Greengard.
17+
18+if nargin == 0, test_chebEval; return; end
19+
20+domain = chebApproxInfo.domain; ax = domain(1); bx = domain(2);
21+ninters = chebApproxInfo.ninters;
22+nord = chebApproxInfo.nord;
23+
24+Nx = length(x);
25+y = zeros(Nx,1);
26+
27+h = (bx-ax)/ninters;
28+
29+for ix = 1:Nx
30+ % get index of left node in current interval
31+ iint = min((floor((x(ix)-ax)/h) + 1),ninters);
32+
33+ % get interval endpoints [a,b]
34+ a = ax + (iint-1)*h; b = ax + iint*h;
35+
36+ % transform to u in [-1,1]
37+ u = (2*x(ix) - a - b)/(b-a);
38+
39+ % Evaluate the cheb polynomial
40+ y(ix) = wts(1,iint);
41+ for J = 2:nord
42+ TJ = cos((J - 1) * acos(u)); % Chebyshev polynomial T_J(X)
43+ y(ix) = y(ix) + TJ * wts(J,iint); % Accumulate the value
44+ end
45+end
46+end
47+
48+function test_chebEval
49+ninters = 12; % number of subintervals
50+nord = 4; % number of cheb nodes in each interval
51+
52+% define the function that we wish to approximate
53+fun = @(x) log(x + 3).*exp(-2*x).*(x.^2).*cos(20*x);
54+
55+% define the domain over which f is defined
56+ax = -2; bx = -1;
57+chebApproxInfo.domain = [ax,bx];
58+chebApproxInfo.ninters = ninters;
59+chebApproxInfo.nord = nord;
60+
61+% table of cheb weights times function values
62+wts = mktab_wts(fun, chebApproxInfo);
63+
64+Nx = 120;
65+x = linspace(ax,bx,Nx)';
66+y = chebEval(x,wts,chebApproxInfo);
67+
68+plot(x,fun(x),'xb',x,y,'or','LineWidth',2);
69+xlabel x; ylabel y; legend('true','approx');
70+title(sprintf('Approximation of f(x) = x^2e^{-2x}cos(20x)log(x + 3) via\n order %d Chebyshev polynomials\n on %d subintervals',nord,ninters))
71+set(gca,'fontsize',15);
72+end
wfp_1Dspringscattering/utils/chexfcdir.madded+41−0View file
@@ -0,0 +1,41 @@
1+function TEXP = chexfcdir(F, N)
2+% CHEXFCDIR computes coefficients of Chebyshev expansion
3+% of function tabulated at classical Chebyshev nodes.
4+%
5+% INPUT:
6+% F = function values at Chebyshev nodes in increasing order
7+% on the interval (-1,1)
8+% N = number of Chebyshev nodes
9+%
10+% OUTPUT:
11+% TEXP = array of expansion coefficients of (N-1)st degree
12+% interpolant
13+%
14+% NOTE: code adjusted from code provided by Leslie Greengard.
15+
16+if nargin == 0, test_chexfcdir; return; end
17+
18+F = F(:); % ensure the input is a column vector
19+
20+FAC = 2/N; % Factor to normalize
21+TEXP = zeros(N, 1); % Initialize TEXP array
22+
23+J = (1:N)'; % cheb node index
24+for I = 1:N % loop for each cheb coefficient
25+ TEXP(I) = sum(FAC * F(J).*cos((J - 0.5)*(I - 1)*pi./N));
26+ TEXP(I) = -TEXP(I)*(-1)^I; % Apply the sign change
27+end
28+TEXP(1) = TEXP(1) / 2; % Adjust the first coefficient
29+
30+end
31+
32+function test_chexfcdir
33+N = 16;
34+[x, ~, ~, ~, ~, ~] = chnodc(-1, 1, N);
35+
36+F = cos(x - 3).*exp(-2*x)';
37+
38+TEXP = chexfcdir(F, N);
39+disp('Array of expansion coefficients of (N-1)st degree interpolant');
40+disp(TEXP);
41+end
wfp_1Dspringscattering/utils/chnodc.madded+62−0View file
@@ -0,0 +1,62 @@
1+function [CHPTS, SINCH, U, V, U1, V1] = chnodc(A, B, M)
2+% CHNODC Constructs Chebyshev nodes and mapping coefficients.
3+%
4+% [CHPTS, SINCH, U, V, U1, V1] = CHNODC(A, B, M)
5+%
6+% This function calculates classical Chebyshev nodes for the
7+% interval [A, B] and the coefficients for linear mappings
8+% between the intervals [-1, 1] and [A, B].
9+%
10+% INPUT PARAMETERS:
11+% A - Lower bound of the interval.
12+% B - Upper bound of the interval.
13+% M - Number of Chebyshev nodes to generate.
14+%
15+% OUTPUT PARAMETERS:
16+% CHPTS - Array of Chebyshev nodes in the interval [A, B].
17+% SINCH - Array containing U * sin(theta) for the Chebyshev points.
18+% U - Coefficient for mapping from [-1, 1] to [A, B].
19+% V - Coefficient for mapping from [-1, 1] to [A, B].
20+% U1 - Coefficient for mapping from [A, B] to [-1, 1].
21+% V1 - Coefficient for mapping from [A, B] to [-1, 1].
22+%
23+% The kth Chebyshev node is computed using:
24+% CHPTS(k) = U * cos((2k - 1) * pi / (2M)) + V
25+%
26+% Example:
27+% A = -1;
28+% B = 1;
29+% M = 5; % Number of Chebyshev nodes
30+% [CHPTS, SINCH, U, V, U1, V1] = chnodc(A, B, M);
31+%
32+% NOTE: code adjusted from code provided by Leslie Greengard.
33+
34+if(nargin == 0), test_chnodc; return; end
35+
36+% Construct the scaling parameters
37+U = (B - A) / 2;
38+V = (B + A) / 2;
39+U1 = 2 / (B - A);
40+V1 = 1 - U1 * B;
41+
42+% Preallocate the arrays for Chebyshev nodes and sin values
43+CHPTS = zeros(1, M);
44+SINCH = zeros(1, M);
45+
46+% Construct the Chebyshev nodes and corresponding SIN array
47+K = 1:M;
48+CHPTS(M - K + 1) = U * cos((2 * K - 1) * pi / (2 * M)) + V;
49+SINCH(M - K + 1) = U * sin((2 * K - 1) * pi / (2 * M));
50+
51+end
52+
53+function test_chnodc
54+clf;
55+
56+A = -2;
57+B = 1;
58+M = 20; % Number of Chebyshev nodes
59+[CHPTS, SINCH, U, V, U1, V1] = chnodc(A, B, M);
60+
61+plot(CHPTS,0,'*r','markerSize',10);
62+end
wfp_1Dspringscattering/utils/generalwindow.madded+86−0View file
@@ -0,0 +1,86 @@
1+function [y,yp,ypp] = generalwindow(x,data)
2+%
3+% [y,yp,ypp] = generalwindow(tol,x,data)
4+%
5+% evaluates blending window on [0,1] which decays to O(tol) at 0^+,
6+% and 1-O(tol) at 1^-. y is vals, yp are derivs, ypp are 2nd derivs.
7+% yp has support in [0,1],
8+% x can be an array, and y, yp, ypp will be the same size as x.
9+%
10+% Note: does not return function handles!
11+%
12+% To test, call without arguments.
13+
14+if nargin==0, test_generalwindow; return; end
15+
16+tol = data.tol;
17+gam = data.gam;
18+
19+% handle values x<0 and x>1...
20+y = 0*x; yp= 0*x; ypp = 0*x; % also makes outputs correct sizes
21+y(x>=1.0) = 1.0;
22+% remaining x's actually need to evaluate something (not 0 or 1)...
23+jj = x >0.0 & x < 1.0; % indices to eval
24+
25+if isempty(data.info) % slow evaluation
26+ theta = log(1/tol); % set up window
27+ W = ceil(2*theta/(pi*gam));
28+ t = W*x(jj); % rescaled ordinates
29+ [phit1,phitt1] = window(theta,W); % get func handles for [0,W]
30+ yp(jj) = W*phit1(t); % make domain of x be [0,1]. w/ Jacobian
31+ ypp(jj) = W*W*phitt1(t);
32+ y(jj) = blending(phit1,t,tol); % must be scalar input
33+
34+else % use data fast eval
35+ t = x(jj); % rescaled ordinates
36+ y(jj) = chebEval(t,data.wei,data.info);
37+ yp(jj) = chebEval(t,data.weip,data.info);
38+ ypp(jj) = chebEval(t,data.weipp,data.info);
39+end
40+
41+end
42+
43+ %%%%%%%%%%%%%
44+function test_generalwindow
45+verb = 1;
46+tol = 1e-12;
47+% setup the data struct
48+data = struct('tol',tol,'gam',0.5,'info',[]);
49+
50+x = -1:1e-2:2.0;
51+[y,yp,ypp] = generalwindow(x,data);
52+if verb,
53+
54+ %figure; plot(x,[y;yp;ypp],'o-'); legend('y','yp','ypp');
55+
56+ %%% testing here
57+ delta = 0.5;
58+ t = x*delta;
59+ [y_new,~,~] = generalwindow((t./delta),data);
60+ figure;
61+ plot(x,y,'o-'); xline(1); hold on;
62+ plot(t,y_new,'r*'); xline(delta)
63+end
64+
65+n=1e3; z = rand(n,1);
66+tic;
67+[Y,Yp,Ypp] = generalwindow(z,data);
68+fprintf("throughput of generalwindow slow = %.3g pts/sec\n",n/toc)
69+
70+% test the fast cheb eval...
71+data = setup_generalwindow(tol);
72+%data
73+%data.info
74+[yf,ypf,yppf] = generalwindow(x,data);
75+
76+disp("test fast vs slow eval...")
77+norm(y-yf,inf)
78+norm(yp-ypf,inf)/norm(yp,inf)
79+norm(ypp-yppf,inf)/norm(ypp,inf)
80+
81+n=1e5; z = rand(n,1);
82+tic;
83+[Y,Yp,Ypp] = generalwindow(z,data);
84+fprintf("throughput of generalwindow cheb = %.3g pts/sec\n",n/toc)
85+end
86+
wfp_1Dspringscattering/utils/glwt.madded+17−0View file
@@ -0,0 +1,17 @@
1+function [x,w] = glwt(a,b,gl)
2+%GLWT convert GL nodes and weights over [-1,1] to ones over [a,b]
3+%
4+%[x,w] = glwt(a,b,x0,w0) returns GL nodes and weights x and w over [a,b]
5+% given GL nodes and weights x0,w0 over [-1,-1]
6+
7+% get nodes and weights
8+x0 = gl.x0;
9+w0 = gl.w0;
10+
11+% Linear map from[-1,1] to [a,b]
12+x=(a*(1-x0)+b*(1+x0))/2;
13+
14+% Compute the weights
15+w=(b-a)*w0;
16+
17+end
wfp_1Dspringscattering/utils/glwt_prep.madded+64−0View file
@@ -0,0 +1,64 @@
1+function gl = glwt_prep(N)
2+% GLWT_PREP to get GL nodes and weights on [-1,1]
3+%
4+% [y,w]=glwt_prep(N)
5+% returns GL nodes 'y' and GL weights 'w' given a number of GL nodes 'N'
6+%
7+% Notes:
8+% This script is for computing definite integrals using Legendre-Gauss
9+% Quadrature. Computes the Legendre-Gauss nodes and weights on an interval
10+% [a,b] with truncation order N
11+%
12+% Suppose you have a continuous function f(x) which is defined on [a,b]
13+% which you can evaluate at any x in [a,b]. Simply evaluate it at all of
14+% the values contained in the x vector to obtain a vector f. Then compute
15+% the definite integral using sum(f.*w);
16+%
17+% Written by Greg von Winckel - 02/25/2004
18+
19+N=N-1;
20+N1=N+1; N2=N+2;
21+
22+xu=linspace(-1,1,N1)';
23+
24+% Initial guess
25+y=cos((2*(0:N)'+1)*pi/(2*N+2))+(0.27/N1)*sin(pi*xu*N/N2);
26+
27+% Legendre-Gauss Vandermonde Matrix
28+L=zeros(N1,N2);
29+
30+% Derivative of LGVM
31+Lp=zeros(N1,N2);
32+
33+% Compute the zeros of the N+1 Legendre Polynomial
34+% using the recursion relation and the Newton-Raphson method
35+
36+y0=2;
37+
38+% Iterate until new points are uniformly within epsilon of old points
39+while max(abs(y-y0))>eps
40+
41+ L(:,1)=1;
42+ Lp(:,1)=0;
43+
44+ L(:,2)=y;
45+ Lp(:,2)=1;
46+
47+ for k=2:N1
48+ L(:,k+1)=( (2*k-1)*y.*L(:,k)-(k-1)*L(:,k-1) )/k;
49+ end
50+
51+ Lp=(N2)*( L(:,N1)-y.*L(:,N2) )./(1-y.^2);
52+
53+ y0=y;
54+ y=y0-L(:,N2)./Lp; % Newton's iteration
55+
56+end
57+
58+w= 1./((1-y.^2).*Lp.^2)*(N2/N1)^2;
59+
60+
61+gl.x0 = y;
62+gl.w0 = w;
63+
64+end
wfp_1Dspringscattering/utils/hline.madded+107−0View file
@@ -0,0 +1,107 @@
1+function hhh=hline(y,in1,in2)
2+% function h=hline(y, linetype, label)
3+%
4+% Draws a horizontal line on the current axes at the location specified by 'y'. Optional arguments are
5+% 'linetype' (default is 'r:') and 'label', which applies a text label to the graph near the line. The
6+% label appears in the same color as the line.
7+%
8+% The line is held on the current axes, and after plotting the line, the function returns the axes to
9+% its prior hold state.
10+%
11+% The HandleVisibility property of the line object is set to "off", so not only does it not appear on
12+% legends, but it is not findable by using findobj. Specifying an output argument causes the function to
13+% return a handle to the line, so it can be manipulated or deleted. Also, the HandleVisibility can be
14+% overridden by setting the root's ShowHiddenHandles property to on.
15+%
16+% h = hline(42,'g','The Answer')
17+%
18+% returns a handle to a green horizontal line on the current axes at y=42, and creates a text object on
19+% the current axes, close to the line, which reads "The Answer".
20+%
21+% hline also supports vector inputs to draw multiple lines at once. For example,
22+%
23+% hline([4 8 12],{'g','r','b'},{'l1','lab2','LABELC'})
24+%
25+% draws three lines with the appropriate labels and colors.
26+%
27+% By Brandon Kuczenski for Kensington Labs.
28+% brandon_kuczenski@kensingtonlabs.com
29+% 8 November 2001
30+
31+if length(y)>1 % vector input
32+ for I=1:length(y)
33+ switch nargin
34+ case 1
35+ linetype='r:';
36+ label='';
37+ case 2
38+ if ~iscell(in1)
39+ in1={in1};
40+ end
41+ if I>length(in1)
42+ linetype=in1{end};
43+ else
44+ linetype=in1{I};
45+ end
46+ label='';
47+ case 3
48+ if ~iscell(in1)
49+ in1={in1};
50+ end
51+ if ~iscell(in2)
52+ in2={in2};
53+ end
54+ if I>length(in1)
55+ linetype=in1{end};
56+ else
57+ linetype=in1{I};
58+ end
59+ if I>length(in2)
60+ label=in2{end};
61+ else
62+ label=in2{I};
63+ end
64+ end
65+ h(I)=hline(y(I),linetype,label);
66+ end
67+else
68+ switch nargin
69+ case 1
70+ linetype='r:';
71+ label='';
72+ case 2
73+ linetype=in1;
74+ label='';
75+ case 3
76+ linetype=in1;
77+ label=in2;
78+ end
79+
80+
81+
82+
83+ g=ishold(gca);
84+ hold on
85+
86+ x=get(gca,'xlim');
87+ h=plot(x,[y y],linetype);
88+ if ~isempty(label)
89+ yy=get(gca,'ylim');
90+ yrange=yy(2)-yy(1);
91+ yunit=(y-yy(1))/yrange;
92+ if yunit<0.2
93+ text(x(1)+0.02*(x(2)-x(1)),y+0.02*yrange,label,'color',get(h,'color'))
94+ else
95+ text(x(1)+0.02*(x(2)-x(1)),y-0.02*yrange,label,'color',get(h,'color'))
96+ end
97+ end
98+
99+ if g==0
100+ hold off
101+ end
102+ set(h,'tag','hline','handlevisibility','off') % this last part is so that it doesn't show up on legends
103+end % else
104+
105+if nargout
106+ hhh=h;
107+end
wfp_1Dspringscattering/utils/interpmat_1d.madded+38−0View file
@@ -0,0 +1,38 @@
1+function L = interpmat_1d(t,s)
2+% INTERPMAT_1D interpolation matrix from nodes in 1D to any target nodes
3+%
4+% L = interpmat_1d(t,s) returns interpolation matrix taking values on nodes s
5+% (a list of nodes) to target nodes t. It assumes smooth functions.
6+% length(s) should be kept small, eg, 30 or less.
7+%
8+% Run without arguments does a self test (see test code for usage example).
9+%
10+% Notes: Computed in Helsing style, with centering and scaling for stability.
11+
12+% Barnett 7/17/16. Auto-centering & scaling for stability 12/23/21.
13+
14+if nargin==0, test_interpmat_1d; return; end
15+
16+cen = (max(s)+min(s))/2; hwid = (max(s)-min(s))/2; % affine s to [-1,1]
17+s = (s-cen)/hwid;
18+t = (t-cen)/hwid;
19+
20+p = numel(s); q = numel(t); s = s(:); t = t(:); % all col vecs
21+n = p; % set the polynomial order we go up to
22+V = ones(p,n); for j=2:n, V(:,j) = V(:,j-1).*s; end % polyval matrix on nodes
23+R = ones(q,n); for j=2:n, R(:,j) = R(:,j-1).*t; end % polyval matrix on targs
24+L = (V'\R')'; % backwards-stable way to do it (Helsing) See corners/interpdemo.m
25+
26+%%%%%%
27+function test_interpmat_1d
28+off = 4.3; % test centering and scaling
29+sc = 1.7;
30+x = sc*linspace(-1,1,16)' + off;
31+f = @(x) sin(x + 0.7);
32+data = f(x); % func on smooth (src) nodes
33+t = sc*(2*rand(1000,1) - 1) + off; % cover same interval as the x lie
34+uex = f(t); % col vec
35+L = interpmat_1d(t,x);
36+u = L * data;
37+fprintf('max abs err for interp in [a,b] : %.3g\n',max(abs(u - uex)))
38+fprintf('interp mat inf-norm = %.3g; max element size = %.3g\n',norm(L,inf),max(abs(L(:))))
wfp_1Dspringscattering/utils/lgwt.madded+59−0View file
@@ -0,0 +1,59 @@
1+function [x,w]=lgwt(N,a,b)
2+
3+% lgwt.m
4+%
5+% This script is for computing definite integrals using Legendre-Gauss
6+% Quadrature. Computes the Legendre-Gauss nodes and weights on an interval
7+% [a,b] with truncation order N
8+%
9+% Suppose you have a continuous function f(x) which is defined on [a,b]
10+% which you can evaluate at any x in [a,b]. Simply evaluate it at all of
11+% the values contained in the x vector to obtain a vector f. Then compute
12+% the definite integral using sum(f.*w);
13+%
14+% Written by Greg von Winckel - 02/25/2004
15+N=N-1;
16+N1=N+1; N2=N+2;
17+
18+xu=linspace(-1,1,N1)';
19+
20+% Initial guess
21+y=cos((2*(0:N)'+1)*pi/(2*N+2))+(0.27/N1)*sin(pi*xu*N/N2);
22+
23+% Legendre-Gauss Vandermonde Matrix
24+L=zeros(N1,N2);
25+
26+% Derivative of LGVM
27+Lp=zeros(N1,N2);
28+
29+% Compute the zeros of the N+1 Legendre Polynomial
30+% using the recursion relation and the Newton-Raphson method
31+
32+y0=2;
33+
34+% Iterate until new points are uniformly within epsilon of old points
35+while max(abs(y-y0))>eps
36+
37+
38+ L(:,1)=1;
39+ Lp(:,1)=0;
40+
41+ L(:,2)=y;
42+ Lp(:,2)=1;
43+
44+ for k=2:N1
45+ L(:,k+1)=( (2*k-1)*y.*L(:,k)-(k-1)*L(:,k-1) )/k;
46+ end
47+
48+ Lp=(N2)*( L(:,N1)-y.*L(:,N2) )./(1-y.^2);
49+
50+ y0=y;
51+ y=y0-L(:,N2)./Lp; % Newton's iteration
52+
53+end
54+
55+% Linear map from [-1,1] to [a,b]
56+x=(a*(1-y)+b*(1+y))/2;
57+
58+% Compute the weights
59+w=(b-a)./((1-y.^2).*Lp.^2)*(N2/N1)^2;
wfp_1Dspringscattering/utils/mktab_wts.madded+57−0View file
@@ -0,0 +1,57 @@
1+function wts = mktab_wts(fun, chebApproxInfo)
2+% MKTAB_WTS make a table of function values over cheb nodes in subintervals
3+% of domain [ax,bx]
4+%
5+% INPUTs:
6+% fun: function handle
7+% chebApproxInfo.domain: domain [ax,bx] overwhich wts are computed.
8+% chebApproxInfo.ninters: number of subintervals
9+% chebApproxInfo.nord: order of polynomial interpolation; number of cheb
10+% nodes
11+%
12+% OUTPUT
13+% wts(norder,ninters): table of cheb weights in each of
14+% the subintervals where cheb poly approximation is applied
15+%
16+% NOTE: code adjusted from code provided by Leslie Greengard.
17+
18+if(nargin == 0), test_mktab_wts; return; end
19+
20+nord = chebApproxInfo.nord;
21+ninters = chebApproxInfo.ninters;
22+domain = chebApproxInfo.domain;
23+ax = domain(1); bx = domain(2);
24+
25+% Define function
26+% fun = @(x) ...;
27+
28+% Preallocate the output matrix
29+wts = zeros(nord, ninters);
30+
31+h = (bx - ax)/ninters;
32+
33+for iint = 1:ninters
34+ A = ax + (iint - 1)*h;
35+ B = ax + iint*h;
36+
37+ % Get Chebyshev nodes and sine values
38+ [CHPTS,~,~,~,~,~] = chnodc(A, B, nord);
39+ wts(:, iint) = chexfcdir(fun(CHPTS), nord);
40+end
41+end
42+
43+function test_mktab_wts
44+ninters = 10; % Number of intervals
45+nord = 5; % Order for Chebyshev nodes
46+fun = @(x) cos(x - 3).*exp(-2*x);
47+
48+ax = -1; bx = 1;
49+chebApproxInfo.nord = nord;
50+chebApproxInfo.ninters = ninters;
51+chebApproxInfo.domain = [ax,bx];
52+
53+wtsTab = mktab_wts(fun, chebApproxInfo);
54+
55+disp('Function at cheb nodes in subintervals:');
56+disp(wtsTab);
57+end
wfp_1Dspringscattering/utils/perispecdiff.madded+25−0View file
@@ -0,0 +1,25 @@
1+function g = perispecdiff(f)
2+% PERISPECDIFF - use FFT to take periodic spectral differentiation of vector
3+%
4+% g = perispecdiff(f) returns g the derivative of the spectral interpolant
5+% of f, which is assumed to be the values of a smooth 2pi-periodic function
6+% at the N gridpoints 2.pi.j/N, for j=1,..,N (or any translation of such
7+% points). Can be row or col vec, and output is same shape.
8+%
9+% Without arguments, does a self-test.
10+
11+% Barnett 2/18/14
12+if nargin==0, test_perispecdiff; return; end
13+N = numel(f);
14+if mod(N,2)==0 % even
15+ g = ifft(fft(f(:)).*[0 1i*(1:N/2-1) 0 1i*(-N/2+1:-1)].');
16+else
17+ g = ifft(fft(f(:)).*[0 1i*(1:(N-1)/2) 1i*((1-N)/2:-1)].');
18+end
19+g = reshape(g,size(f));
20+
21+%%%%%%
22+function test_perispecdiff
23+N = 50; tj = 2*pi/N*(1:N)';
24+f = sin(3*tj); fp = 3*cos(3*tj); % trial periodic function & its deriv
25+norm(fp-perispecdiff(f))
wfp_1Dspringscattering/utils/setup_generalwindow.madded+37−0View file
@@ -0,0 +1,37 @@
1+function [data,W] = setup_generalwindow(tol)
2+
3+% [data,W] = SETUP_GENERALWINDOW(tol)
4+% set up fast cheb eval for window
5+%
6+% See: generalwindow.m
7+
8+gam = 0.5; data.gam = 0.5;
9+data.tol = tol;
10+data.info = {};
11+
12+chebApproxInfo.domain = [0,1];
13+chebApproxInfo.ninters = max(1,ceil(0.2*log(1/tol))); % *** hack to make sure accurate to tol!;
14+chebApproxInfo.nord = 16;
15+
16+% table of cheb weights times function values
17+data.wei = mktab_wts(@(x) yfun(x,data), chebApproxInfo);
18+data.weip = mktab_wts(@(x) ypfun(x,data), chebApproxInfo);
19+data.weipp = mktab_wts(@(x) yppfun(x,data), chebApproxInfo);
20+data.info = chebApproxInfo;
21+
22+
23+% Some more information needed
24+theta = log(1/tol); % set up window
25+W = ceil(2*theta/(pi*gam));
26+
27+end
28+
29+function y = yfun(x,data)
30+ [y,~,~] = generalwindow(x,data);
31+end
32+function yp = ypfun(x,data)
33+ [~,yp,~] = generalwindow(x,data);
34+end
35+function ypp = yppfun(x,data)
36+ [~,~,ypp] = generalwindow(x,data);
37+end
wfp_1Dspringscattering/utils/vline.madded+107−0View file
@@ -0,0 +1,107 @@
1+function hhh=vline(x,in1,in2)
2+% function h=vline(x, linetype, label)
3+%
4+% Draws a vertical line on the current axes at the location specified by 'x'. Optional arguments are
5+% 'linetype' (default is 'r:') and 'label', which applies a text label to the graph near the line. The
6+% label appears in the same color as the line.
7+%
8+% The line is held on the current axes, and after plotting the line, the function returns the axes to
9+% its prior hold state.
10+%
11+% The HandleVisibility property of the line object is set to "off", so not only does it not appear on
12+% legends, but it is not findable by using findobj. Specifying an output argument causes the function to
13+% return a handle to the line, so it can be manipulated or deleted. Also, the HandleVisibility can be
14+% overridden by setting the root's ShowHiddenHandles property to on.
15+%
16+% h = vline(42,'g','The Answer')
17+%
18+% returns a handle to a green vertical line on the current axes at x=42, and creates a text object on
19+% the current axes, close to the line, which reads "The Answer".
20+%
21+% vline also supports vector inputs to draw multiple lines at once. For example,
22+%
23+% vline([4 8 12],{'g','r','b'},{'l1','lab2','LABELC'})
24+%
25+% draws three lines with the appropriate labels and colors.
26+%
27+% By Brandon Kuczenski for Kensington Labs.
28+% brandon_kuczenski@kensingtonlabs.com
29+% 8 November 2001
30+
31+if length(x)>1 % vector input
32+ for I=1:length(x)
33+ switch nargin
34+ case 1
35+ linetype='r:';
36+ label='';
37+ case 2
38+ if ~iscell(in1)
39+ in1={in1};
40+ end
41+ if I>length(in1)
42+ linetype=in1{end};
43+ else
44+ linetype=in1{I};
45+ end
46+ label='';
47+ case 3
48+ if ~iscell(in1)
49+ in1={in1};
50+ end
51+ if ~iscell(in2)
52+ in2={in2};
53+ end
54+ if I>length(in1)
55+ linetype=in1{end};
56+ else
57+ linetype=in1{I};
58+ end
59+ if I>length(in2)
60+ label=in2{end};
61+ else
62+ label=in2{I};
63+ end
64+ end
65+ h(I)=vline(x(I),linetype,label);
66+ end
67+else
68+ switch nargin
69+ case 1
70+ linetype='r:';
71+ label='';
72+ case 2
73+ linetype=in1;
74+ label='';
75+ case 3
76+ linetype=in1;
77+ label=in2;
78+ end
79+
80+
81+
82+
83+ g=ishold(gca);
84+ hold on
85+
86+ y=get(gca,'ylim');
87+ h=plot([x x],y,linetype);
88+ if length(label)
89+ xx=get(gca,'xlim');
90+ xrange=xx(2)-xx(1);
91+ xunit=(x-xx(1))/xrange;
92+ if xunit<0.8
93+ text(x+0.01*xrange,y(1)+0.1*(y(2)-y(1)),label,'color',get(h,'color'))
94+ else
95+ text(x-.05*xrange,y(1)+0.1*(y(2)-y(1)),label,'color',get(h,'color'))
96+ end
97+ end
98+
99+ if g==0
100+ hold off
101+ end
102+ set(h,'tag','vline','handlevisibility','off')
103+end % else
104+
105+if nargout
106+ hhh=h;
107+end
wfp_1Dspringscattering/utils/window.madded+20−0View file
@@ -0,0 +1,20 @@
1+function [phit,phitt] = window(b,w)
2+% WINDOW1 generates the Kaiser Bessel window function and its time
3+% derivative
4+%
5+% [phit,phitt] = window(b,w)
6+% returns the Kaiser Bessel window function phit and its time derivative
7+% phitt whose support is w.
8+%
9+% Input values:
10+% b = log(1/tol), where tol phit(0) = phit(delta) = tol
11+% w : the width of the support of phit
12+%
13+% Note: fuNction rewritten from a Julia code provided by Alex Barnett
14+
15+cen = w/2.0; % center of the window function
16+prefac = (1/w) * b/(sinh(b));
17+phit = @(x) prefac*besseli(0,b*sqrt(1-((2/w)*(x-cen)).^2));
18+phitt = @(x) prefac*besseli(1,b*sqrt(1-((2/w)*(x-cen)).^2)).*(-4*b*(x - cen))./((w^2)*sqrt(1-((2/w)*(x-cen)).^2));
19+
20+end