1function [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)
24% Get GL nodes and weights on the interval [0,dt]
25[tau,w] = glwt(0,dt,gl);
27delta = W*dt;
29k = K';
30tau = tau';
31lK = length(K);
33% allocate space for outputs
34p = zeros(lK,(W+1)); q = zeros(lK,(W+1));
35p0 = zeros(W+1,1); q0 = p0;
37for i = 0:W
38 gam = i*dt;
40 winArg = (tau + gam)./delta;
41 [~,phit,phitt] = generalwindow(winArg,winData);
42 phit = (1/delta)*phit;
43 phitt = (1/delta^2)*phitt;
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;
49 I0 = (2*phit + (tau+gam).*phitt);
50 p0(i+1) = sum(((dt - tau).*I0)'.*w);
51 q0(i+1) = sum(I0'.*w);
52end
54sn_hat = finufft1d1(s,sn(snIdxNow-(0:W),:)',1,tol,N,struct('modeord',1));
56%%% free-space parameters
57phi_grid = (1:W) - 0.5;
58[phi_tn1,phit_tn1,~] = generalwindow(phi_grid./W,winData);
59phi_tn = phi_tn1'; phit_tn = (1/W)*phit_tn1';
61end