% GDF of a spectrometer (1D deconvolution). 
% This is done naively, as a result it is extremely sensitive
% to noise. Best results are obtained with > 1e6 peak counts per bin in
% the spectral lines.

max_counts = 1e5; % Maximum counts per bin in the unblurred spectrum.
continuum = max_counts/20; % Continuum intensity, counts per bin

dlambda = 0.05; %wavelength resolution in nm
lambda0 = 152; %low end of measured wavelengths in nm
lambda1 = 160; %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 properly normalized GDF
width = 0.5; %width of scattering core in nm
f_gdf = exp(-abs(k*width/(pi)));
gdf = real(ifftn(f_gdf));

figure(1)
hold off
plot(lambda, gdf);
hold on
xlabel('\lambda, nm')
ylabel('scattering function, 1/nm')
print('deconvolution_naif_fig1.pdf', '-dpdfwrite')

% Create spectrum with two spectral lines (C IV) and continuum.
linewidth = 0.1; %gaussian (thermal) line width, nm
center1 = 154.820; %line center, nm
center2 = 155.077
height = max_counts; %counts per bin
spectrum = height.*exp(-(lambda - center1).^2./(2*linewidth^2))...
   + height.*exp(-(lambda - center2).^2./(2*linewidth^2))...
   + continuum; 


% Simulated observation
f_spec = fftn(spectrum);
f_spec_meas = f_spec .* f_gdf; %Convolve with instrument gdf
spec_meas = real(ifftn(f_spec_meas));
spec_meas = poissrnd(spec_meas); %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')
print('deconvolution_naif_fig2.pdf', '-dpdfwrite')


% Attempt to recover spectrum from measured (f_spec_meas),
% using formal inverse with FFT:
f_recovered = f_spec_meas ./ f_gdf;
recovered = real( ifftn(f_recovered) );

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


