function[s_p,Fsp]= nonstat_transfer(x,f,H,C,cota,TFtype) % Nonstationary integral transform along flowline to predict surface response s_p(x) % (and its Fourier Transform Fsp) due to bed perturbation forcing f. % % Inputs: x - row vector of horizontal distance % H, C, cota - row vectors of background variables (ice thickness, slip ratio and cot(alpha) respectively) % *** NOTE that H must have the same unit as x *** % f - row vector of the basal forcing % TFtype - type of transfer function, either 'TSB' or 'TSC'. (This specifies whether the forcing f is bed-topographic perturbation (b) or bed-slipperiness perturbation (c).) % % Outputs: s_p - surface response, at those positions specified in x % Fsp - Discrete Fourier transform of s_p % (1) set up dx = x(2)-x(1); N = length(x); % (2) dimensional wavenumber axis dk = 2*pi/(N*dx); if rem(N,2) == 0 % N even k = dk * [ 0:1:N/2 , -(N/2-1):1:-1]'; else % N odd k = dk * [ 0:1:(N-1)/2 , -(N-1)/2:1:-1]'; end % (3) compile transfer matrix M = T(k,x)*exp(-ikx), each column of M referring to a position x M = zeros(length(k),N); for j = 1:N % call transfer function TSB or TSC, where input k is dimensionLESS knd = H(j) * k; eval([ 'M(:,j) = ' , TFtype , '(knd,C(j),cota(j));' ]) M(:,j) = M(:,j) .* exp(-i*k*x(j)); end % (4) apply nonstationary filter if TFtype == 'TSC' Fsp = dx * (M * (H.*f)').'; elseif TFtype == 'TSB' Fsp = dx * (M * f').'; end % being careful here, using [.'] to transpose the resulting vector without complex conjugation (because ' means CONJUGATE TRANSPOSE in Matlab) % inverse DFT to find predicted surface topographic perturbation s_p = real(ifft(Fsp)) / dx;