% Nonspecular scattered light in a spectrometer in one dimension. 
% Modeled (somewhat unrealistically) as a Fredholm equation of the
% second kind.
alpha = 0.5;
beta = 0.5;
max_counts = 1e3; % Maximum counts per bin in the unblurred spectrum.
Niter = 5; % Number of iterations in the Neumann series

dlambda = 0.001; %wavelength resolution in nm
lambda0 = 115; %low end of measured wavelengths in nm
lambda1 = 125; %high end of wavelength range in nm
lambda = lambda0:dlambda:lambda1;
N = numel(lambda); %number of wavelengths
k = k_array(N, dlambda); %coordinate of fft

% Create scattering wings of GDF
width = 1.0; %width of scattering core in nm
f_scat = exp(-abs(k*width/(pi)));
scat = real(ifftn(f_scat));

figure(1)
hold off
plot(lambda, scat);
hold on
xlabel('\lambda, nm')
ylabel('scattering function, 1/nm')
%orient tall
%orient landscape
print('scattering_fig1.pdf', '-dpdfwrite')

% Create spectrum with one spectral line (lyman alpha)
linewidth = 0.02; %gaussian (thermal) line width, nm
center = 121.6; %line center, nm
height = max_counts; %counts per bin
spectrum = height.*exp(-(lambda - center).^2./(2.*linewidth.^2));


% Convolve with instrument GDF
f_spec = fftn(spectrum);
f_spec_meas = f_spec .* (alpha + beta .* f_scat);
spec_meas = real(ifftn(f_spec_meas));
spec_meas = poissrnd(spec_meas, 1, N); %add Poisson noise
f_spec_meas = fftn(spec_meas); %Since we use f_spec_meas for
   %our data analysis, it is imperative to put the noise into it!

figure(2)
hold off
plot(lambda, spectrum, 'g-', lambda, spec_meas,'r-')
hold on
xlabel('\lambda, nm')
ylabel('spectrum, counts per bin')
legend('true spectrum','measured spectrum')
%orient tall
%orient landscape
print('scattering_fig2.pdf', '-dpdfwrite')


% Attempt to recover spectrum from measured (f_spec_meas),
% using Neumann series:
f_recovered = f_spec_meas;
for n=1:Niter
   f_recovered = f_spec_meas/alpha - ...
      ( beta * f_scat .* f_recovered) / alpha;
end
recovered = real( ifftn(f_recovered) );

figure(3)
hold off
plot(lambda, spectrum, 'g-', lambda, ...
   spec_meas,'r-', lambda, recovered,'b-')
hold on
axis([center - 10*linewidth, center + 10*linewidth,...
   min(recovered), max(recovered)])
xlabel('\lambda, nm')
ylabel('spectrum, counts per bin')
legend('true spectrum','measured spectrum','recovered spectrum')
%orient tall
%orient landscape
print('scattering_fig3.pdf', '-dpdfwrite')


