8748c66added codenalhassanieh 1function [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
26if nargin == 0, test_alphaEvaluation_spectral; return; end
28%% Treatment of u_history
29k = K';
31% update k values not equal to zero
32i = 0:W;
33h = sum(p(:,i+1).*sn_hat,2); % sum rows
34g = sum(q(:,i+1).*sn_hat,2);
35h = h*dt; g = g*dt;
37% evolution equations k ~= 0
38anp1 = an.*cos(k*dt) + bn.*(1./k).*sin(k*dt) + h;
39bnp1 = -k.*an.*sin(k*dt) + bn.*cos(k*dt) + g;
41% update k = 0 frequency
42sig0 = sum(sn(snIdxNow-i,:),2); % sum rows
43h0 = sum(p0(i+1).*sig0);
44g0 = sum(q0(i+1).*sig0);
45h0 = h0*dt; g0 = g0*dt;
47% evolution equations for k = 0 case
48anp1(zeroLoc) = an(zeroLoc) + bn(zeroLoc)*dt + h0;
49bnp1(zeroLoc) = bn(zeroLoc) + g0;
51end