function [lambda_pos,energy,dE,spectrum,corrected_spectrum]=reconstruct_spectra_at677(base_current,polarisation,E0,alpha1,opt)


% [lambda_pos,energy,dE,spectrum,corrected_spectrum]=reconstruct_spectra_at677(base_current,polarisation,E0,alpha1,opt)
% ------------------------------------------------------------------------
% input parameters:
% base_current: current vector
% polarisation: vector wiht the measured polarisation (real + 1i*imag)
% E0: incident beam energy in meV
% alpha1: tilt angle in degree
% opt: field with optional values
% default opt.forcesymm = 1 
% forces real part to be symmetric and imag part to be antisymmetric
% default opt.binning = 1 
% applys the constant binning method for the energy spectrum
% --------------------------------------------------------------------------
% output parameters:
% lambda_pos: lambda_final of the scattered He atoms
% energy: final energy of the scattered He atoms
% dE: energy exchange E_final - E_initial
% spectrum: energy spectrum (without correction)
% corrected_spectrum: energy spectrum corrected by the Jacobian


if nargin<4;
   error('too few parameters');
end

%setting default values if no options are passed
if nargin == 4;
   opt.forcesymm = 0;
   opt.binning = 1;
end

if nargin > 4;
   if ~isfield(opt,'forcesymm');
       opt.forcesymm = 1;
   end
   if ~isfield(opt,'binning');
       opt.binning = 1;
   end
end

N=length(polarisation);
load_chess_parameters;
m = SE_amu * 3.01603;
alpha1=alpha1*pi/180;

% forcing real part to be symmetric and imag part to be antisymmetric
%###########################################################################################################################
if opt.forcesymm==1;
     polarisation=polarisation(:)'; %cause we like to work with a column vector
     real_sig=real(polarisation); imag_sig=imag(polarisation);
     real_corrected=mean([real_sig(1:ceil(length(real_sig)/2));fliplr(real_sig(ceil(length(real_sig)/2):end))],1);
     real_corrected=[real_corrected,fliplr(real_corrected(1:end-1))];
     imag_corrected=mean([imag_sig(1:ceil(length(imag_sig)/2));-1*fliplr(imag_sig(ceil(length(imag_sig)/2):end))],1);
     imag_corrected=[imag_corrected,-1*fliplr(imag_corrected(1:end-1))];
     polarisation=real_corrected+1i*imag_corrected;
     polarisation=polarisation(:)'; %bringing everythihng back into shape
end
%###########################################################################################################################

%win=sigwin.hann(length(base_current),'symmetric');
%window=generate(win);
window=hann_n(length(base_current))';
bs=polarisation;
bs=polarisation.*window; %applying a window function to force the spectrum to go to zero at the edges
%bsf=fftshift(fft(fftshift(bs))/N); %doing the actual fft of the data...
%should just give you an offset
%bsf=fftshift(fft(fftshift(bs))); %doing the actual fft of the data
bsf=fftshift(fft(fftshift(bs))); %doing the actual fft of the data
%what about the issue of the last point being the same as the first and the
%fft window...
%bsf=fftshift(fft(fftshift(bs(1:end-1)))); %doing the actual fft of the data
%bsf=[bsf,NaN];

wavelength_res=1/(2*max(base_current)*((SE_gamma/2/pi)*SE_Beff1*m/SE_h));
N=round(length(base_current)/2);
exp_lambda=linspace(-(N-1)*wavelength_res,(N-1)*wavelength_res,length(base_current)); %wavlength vector proportional to current vector


K0=beamprops('energy',E0,3)*1e10;
V0=SE_hbar*K0/m;
lambda0=SE_h/(m*V0);
%accounts for the fact that we are measuring along a tilt angle... 
%project lambda back according to the tilt angle
alpha=2*pi-alpha1;
%alpha=pi-alpha1;
lambda_shift=(exp_lambda./cos(alpha-pi/2)+lambda0*tan(pi/2-alpha)); %eq. 33 of Andys review


lambda_pos=lambda_shift(lambda_shift>0);
spectrum=bsf;
spectrum=spectrum(lambda_shift>0);
energy=6.2415e21*SE_h^2./(2*m*lambda_pos.^2);
dE=energy-E0;
jacobian=(SE_h/m^2)*(2*energy/m).^(-3/2);
corrected_spectrum=spectrum.*jacobian;


%applying the constant binning method
%###########################################################################################################################
if opt.binning==1;
    %let's put dE and spectrum in a neat ascending order first
    [dE,ind]=sort(dE); corrected_spectrum=corrected_spectrum(ind);
    %constant binning only on creation side...
    %binsize=0.05; %change the binning size (in meV) if you want to
    binsize=60/N;
    Enew=[min(dE):binsize:-0.2];
    for kk=1:length(Enew)-1;
        index = (Enew(kk) <= dE)&(dE <= Enew(kk+1));
        ynew(kk)=mean(corrected_spectrum(index));
    end
    Enew=Enew+binsize/2; Enew(end)=[];
    index=length(dE)-sum(dE>Enew(end));
    Enew=[Enew,dE(index:end)]; ynew=[ynew,corrected_spectrum(index:end)];
    dE=Enew;
    corrected_spectrum=abs(ynew);
end
%###########################################################################################################################