function so = average_raw(si, whichway,bounds,npts)
%SAXS/AVERAGE_RAW Average SAXS image intensities- Spectra-by-spectra
%   AVERAGE(S, 'x',points) returns a SAXS image object S with image
%   columns collapsed and replaced by the average intensity in that column.
%   It uses the original 2D data and bins all the pixels into their
%   respective positions to create a 1D data set.
%
%   In other words, AVERAGE_raw returns the one-dimensional curve of average
%   image intensity versus whatever the x coordinate is, usually momentum
%   transfer (image radius).
%
%   This routine draws heavily on a old routine from Karsten Joensen 
%   which takes the 2D image, the center position and the calibrations
%   and then calculates the r,phi coordinates for each pixel.
%   These routines then step through each point in the 2D image and throws
%   it into the correct "bin" in the 1D spectrum.
%   This approach is not elegant but it works within a reasonable amount of
%   time. It could be sped up by making these routines into C routines and
%   incorporating them as mex files..will be done later
%
% As of September 25,2006 we have included an si.error_norm_fact
% which is multiplied onto the relative error calculated by
% assuming normally distributed pixel values. If there are
% actually 100 photons per pixel-value=1 then the
% si.error_norm_fact will be 1/sqrt(100)


if isempty(si.reduction_state)
    si.reduction_state='000000000000';
end
so=si;
if nargin == 1
    whichway = 'x';  % default
    npts=200; %default
elseif nargin == 2
    npts=200; %default
elseif ~ (isequal(whichway, 'x') || isequal(whichway, 'y'))
    error('The second argument must be ''x'' or ''y''.')
end

error_norm_fact=si.error_norm_fact;

%setting up parameters to calculate both average and relative error


stemp = si;
if isempty(si.livetime)
    stemp.livetime=1;
end

%This hack is included to fix a bug when the user has not entered any calibrations.
if isempty(si.pixelcal)
    si.kcal=1.;
    si.pixelcal=[1.0 1.0];
    si.koffset=0;
    si.wavelength=[];
else
    %this is a correction to allow for WAXS 
    %si.kcal becomes deltaq_at_0*si.pixelcal
    %please note this does not change the so.kcal value
    si.kcal=tan(2*asin(si.kcal*si.wavelength/4/pi))/2*4*pi/si.wavelength;
end
warning off MATLAB:divideByZero
switch whichway
    case 'x'
        if str2num(si.reduction_state(12))==0 % if the beamstopmask is not to be used
            stemp.bsmask_pixels=1+0.*si.raw_pixels;
        end
        
        if str2num(si.reduction_state(8))==0 % if the mask is not to be used
            stemp.mask_pixels=1+0.*si.raw_pixels;
        end

        % this is kept as a reminder that pixel by pixel flatfield division could be used
        stemp.flatfield_pixels=1+0*si.raw_pixels;   %...but is not    
        flatmask=(stemp.bsmask_pixels.*stemp.mask_pixels)./stemp.flatfield_pixels; 
       
        %if readoutnoise
        if str2num(si.reduction_state(9))==1    %this means that readoutnoise should be subtracted !!!
            pixels=prepare_pixels(si.readoutnoise,flatmask,si.reduction_state);
            [xreadoutnoise,yreadoutnoiseave,yreadoutnoiserelerr,yreadoutnoiserelerr2]=...
                azimbinning(pixels.*flatmask,bounds(1),bounds(1)+bounds(3),npts,...
                si.xaxisbintype,bounds(2),bounds(2)+bounds(4),si.raw_center(2),si.raw_center(1),...
                si.kcal/si.pixelcal(2),si.kcal/si.pixelcal(1),si.wavelength,error_norm_fact);
        else
            yreadoutnoiseave=0;
            yreadoutnoiserelerr=0;
            yreadoutnoiserelerr2=0;
        end

        % now we are ready look at some data.
        pixels=prepare_pixels(si.raw_pixels,flatmask,si.reduction_state); %maybe do readout_subtraction and zinger_removal
        [x,ysampleave,ysamplerelerr,ysamplerelerr2]=azimbinning(pixels.*flatmask,bounds(1),bounds(1)+bounds(3),npts,...
            si.xaxisbintype,bounds(2),bounds(2)+bounds(4),si.raw_center(2),si.raw_center(1),...
            si.kcal/si.pixelcal(2),si.kcal/si.pixelcal(1),si.wavelength,error_norm_fact);
        warning off MATLAB:divideByZero
        ysamplerelerr=sqrt((ysamplerelerr.*ysampleave).^2+(yreadoutnoiseave.*yreadoutnoiserelerr).^2)./abs(ysampleave-yreadoutnoiseave);
        ysamplerelerr2=sqrt((ysamplerelerr2.*ysampleave).^2+(yreadoutnoiseave.*yreadoutnoiserelerr2).^2)./abs(ysampleave-yreadoutnoiseave);
        ysampleave=ysampleave-yreadoutnoiseave;

        if (str2num(si.reduction_state(5))==0 || isempty(stemp.transfact))  % if sample transmission is not to be used
            stemp.transfact=1;
        end
        
        if str2num(si.reduction_state(7))==1 %if flatfield correction average is to be used.
            stemp.flatfield_pixels=1+0*si.raw_pixels; % not implemented at this time
        end
        
        if str2num(si.reduction_state(4))==1 % if empty holder subtraction is to be used.
            empty_corr=1;
            pixels=prepare_pixels(si.empty_pixels,flatmask,si.reduction_state); %maybe do readout_subtraction and zinger_removal
            [xempty,yemptyave,yemptyrelerr,yemptyrelerr2]=azimbinning(pixels.*flatmask,bounds(1),bounds(1)+bounds(3),npts,...
                si.xaxisbintype,bounds(2),bounds(2)+bounds(4),si.raw_center(2),si.raw_center(1),...
                si.kcal/si.pixelcal(2),si.kcal/si.pixelcal(1),si.wavelength,error_norm_fact);
            warning off MATLAB:divideByZero
            yemptyrelerr=sqrt((yemptyrelerr.*yemptyave).^2+(yreadoutnoiseave.*yreadoutnoiserelerr).^2)./abs(yemptyave-yreadoutnoiseave);
            yemptyrelerr2=sqrt((yemptyrelerr2.*yemptyave).^2+(yreadoutnoiseave.*yreadoutnoiserelerr2).^2)./abs(yemptyave-yreadoutnoiseave);
            yemptyave=yemptyave-yreadoutnoiseave;
        else
            empty_corr=0;
            yemptyave=0;
            yemptyrelerr=0; % this does not take into account the error on the transmission number
            yemptyrelerr2=0; % this does not take into account the error on the transmission number
            stemp.empty_transfact=1; 
            stemp.empty_livetime=1;
        end
        
        if str2num(si.reduction_state(11))==1 % if solvent subtraction is to be used.
            pixels=prepare_pixels(si.solvent_pixels,flatmask,si.reduction_state); %check  do readout_subtraction and zinger_removal
            [xsolvent,ysolventave,ysolventrelerr,ysolventrelerr2]=azimbinning(pixels.*flatmask,bounds(1),bounds(1)+bounds(3),npts,...
                si.xaxisbintype,bounds(2),bounds(2)+bounds(4),si.raw_center(2),si.raw_center(1),...
                si.kcal/si.pixelcal(2),si.kcal/si.pixelcal(1),si.wavelength,error_norm_fact);
            warning off MATLAB:divideByZero
            ysolventrelerr=sqrt((ysolventrelerr.*ysolventave).^2+(yreadoutnoiseave.*yreadoutnoiserelerr).^2)./abs(ysolventave-yreadoutnoiseave);
            ysolventrelerr2=sqrt((ysolventrelerr2.*ysolventave).^2+(yreadoutnoiseave.*yreadoutnoiserelerr2).^2)./abs(ysolventave-yreadoutnoiseave);
            ysolventave=ysolventave-yreadoutnoiseave;
        else
            ysolventave=0;
            ysolventrelerr=0; % this does not take into account the error on the transmission number
            ysolventrelerr2=0; % this does not take into account the error on the transmission number
            stemp.solvent_transfact=1; 
            stemp.solvent_livetime=1;
        end
        
        if str2num(si.reduction_state(3))==1 % if darkcurrent correction average is to be used.
            pixels=prepare_pixels(si.darkcurrent_pixels,flatmask,si.reduction_state); %maybe do readout_subtraction and zinger_removal
            [xdarkcurrent,ydarkcurrentave,ydarkcurrentrelerr,ydarkcurrentrelerr2]=...
                azimbinning(pixels.*flatmask,bounds(1),bounds(1)+bounds(3),npts,...
                si.xaxisbintype,bounds(2),bounds(2)+bounds(4),si.raw_center(2),si.raw_center(1),...
                si.kcal/si.pixelcal(2),si.kcal/si.pixelcal(1),si.wavelength,error_norm_fact);
            warning off MATLAB:divideByZero
            ydarkcurrentrelerr=sqrt((ydarkcurrentrelerr.*ydarkcurrentave).^2+(yreadoutnoiseave.*yreadoutnoiserelerr).^2)./abs(ydarkcurrentave-yreadoutnoiseave);
            ydarkcurrentrelerr2=sqrt((ydarkcurrentrelerr2.*ysolventave).^2+(yreadoutnoiseave.*yreadoutnoiserelerr2).^2)./abs(ydarkcurrentave-yreadoutnoiseave);
            ydarkcurrentave=ydarkcurrentave-yreadoutnoiseave;
        else
            ydarkcurrentave=0;
            ydarkcurrentrelerr=0; % this does not take into account the error on the transmission number
            ydarkcurrentrelerr2=0; % this does not take into account the error on the transmission number
            stemp.darkcurrent_livetime=1;
        end
                
        if   ~isempty(si.flatfield_filename)&& (str2num(si.reduction_state(6))==1 || str2num(si.reduction_state(7))==1) % if the flatfield is to be used 
            
            %maxflatfield=mean(mean(si.flatfield_pixels(floor(si.raw_center(1)-50):floor(si.raw_center(1)+50),floor(si.raw_center(2)-50):floor(si.raw_center(2)+50))));
            maxflatfield=ffvalue(si,3); % maxflatfield is actually not the maximum but rather an "average" defined in ffvalue.m
            flatmask=flatmask/maxflatfield;
            pixels=prepare_pixels(si.flatfield_pixels,flatmask,si.reduction_state); %maybe do readout_subtraction and zinger_removal
            [xflatfield,yflatfieldave,yflatfieldrelerr,yflatfieldrelerr2]=...
                azimbinning(pixels.*flatmask,bounds(1),bounds(1)+bounds(3),npts,...
                si.xaxisbintype,bounds(2),bounds(2)+bounds(4),si.raw_center(2),si.raw_center(1),...
                si.kcal/si.pixelcal(2),si.kcal/si.pixelcal(1),si.wavelength,error_norm_fact);
                yflatfieldave=yflatfieldave-yreadoutnoiseave;

            if str2num(si.reduction_state(7))==1      %if flatfield correction average is to be used
                %here we try to fit in the average plot to a 2nd order
                %polynomial
                P=polyfit(xflatfield(find(~isnan(yflatfieldave))),yflatfieldave(find(~isnan(yflatfieldave)))',2); % fit to a second order polynomial
                yflatfieldold=yflatfieldave;
                yflatfieldave=polyval(P,xflatfield)';
            end
            warning off MATLAB:divideByZero
            ysampleave=ysampleave./yflatfieldave;
            yemptyave=yemptyave./yflatfieldave;
            ysolventave=ysolventave./yflatfieldave;
            ydarkcurrentave=ydarkcurrentave./yflatfieldave;
            
            ysamplerelerr=ysamplerelerr./yflatfieldave;
            yemptyrelerr=yemptyrelerr./yflatfieldave;
            ysolventrelerr=ysolventrelerr./yflatfieldave;
            ydarkcurrentrelerr=ydarkcurrentrelerr./yflatfieldave;
            ysamplerelerr2=ysamplerelerr2./yflatfieldave;
            yemptyrelerr2=yemptyrelerr2./yflatfieldave;
            ydarkcurrentrelerr2=ydarkcurrentrelerr2./yflatfieldave;
        end
        % no considerations are taken as to how the flatfield division
        % affects the errorbars...this should be fixed at a later stage.
        
        %
        % OK so now we have everything needed to do the correct
        % subtractions: yave is the average corrected value per pixel in
        % the whole scan...to find the value per pixel per second one has
        % to divide by livetime.

        % This is the case with sample and background only
        if isempty(si.solvent_filename)
            warning off MATLAB:divideByZero
            yave=(ysampleave-(stemp.livetime/stemp.darkcurrent_livetime)*ydarkcurrentave)/stemp.transfact-...
                 empty_corr*(yemptyave-(stemp.empty_livetime/stemp.darkcurrent_livetime)*ydarkcurrentave)*(stemp.livetime/stemp.empty_livetime/stemp.empty_transfact);
            warning off MATLAB:divideByZeroa
            yrelerr=sqrt(((ysampleave.*ysamplerelerr).^2+(stemp.livetime/stemp.darkcurrent_livetime).^2*(ydarkcurrentave.*ydarkcurrentrelerr).^2)/stemp.transfact.^2+...
                empty_corr*((yemptyave.*yemptyrelerr).^2+(stemp.empty_livetime/stemp.darkcurrent_livetime).^2*(ydarkcurrentave.*ydarkcurrentrelerr).^2)*...
                (stemp.livetime/stemp.empty_livetime/stemp.empty_transfact).^2)./yave;
            warning off MATLAB:divideByZero
            yrelerr2=sqrt(((ysampleave.*ysamplerelerr2).^2+(stemp.livetime/stemp.darkcurrent_livetime).^2*(ydarkcurrentave.*ydarkcurrentrelerr2).^2)/stemp.transfact.^2+...
                empty_corr*((yemptyave.*yemptyrelerr2).^2+(stemp.empty_livetime/stemp.darkcurrent_livetime).^2*(ydarkcurrentave.*ydarkcurrentrelerr2).^2)*...
                (stemp.livetime/stemp.empty_livetime/stemp.empty_transfact).^2)./yave;

            % This is the case with sample and background and solvent
        else 

            warning off MATLAB:divideByZero
            %subtracting darkcurrent from each one
            ysampleavedc=(ysampleave-(stemp.livetime/stemp.darkcurrent_livetime)*ydarkcurrentave);
            ysolventavedc=(ysolventave-(stemp.solvent_livetime/stemp.darkcurrent_livetime)*ydarkcurrentave);
            yemptyavedc=(yemptyave-(stemp.empty_livetime/stemp.darkcurrent_livetime)*ydarkcurrentave);
            yave=ysampleavedc/stemp.transfact-...
                ysolventavedc/stemp.solvent_transfact*(stemp.livetime/stemp.solvent_livetime)*(stemp.sample_thickness/stemp.solvent_thickness)-...
                yemptyavedc/stemp.empty_transfact*(stemp.livetime/stemp.empty_livetime)*(1-stemp.sample_thickness/stemp.solvent_thickness);
            warning off MATLAB:divideByZero
            ysampleerrdc=(ysampleave.*ysamplerelerr).^2+(stemp.livetime/stemp.darkcurrent_livetime).^2*(ydarkcurrentave.*ydarkcurrentrelerr).^2;
            ysolventerrdc=(ysolventave.*ysolventrelerr).^2+(stemp.solvent_livetime/stemp.darkcurrent_livetime).^2*(ydarkcurrentave.*ydarkcurrentrelerr).^2;
            yemptyerrdc=(yemptyave.*yemptyrelerr).^2+(stemp.empty_livetime/stemp.darkcurrent_livetime).^2*(ydarkcurrentave.*ydarkcurrentrelerr).^2;
            yerr=ysampleerrdc/stemp.transfact.^2+...
                ysolventerrdc/stemp.solvent_transfact.^2*(stemp.livetime/stemp.solvent_livetime)*(stemp.sample_thickness/stemp.solvent_thickness)^2+...
                yemptyerrdc/stemp.empty_transfact.^2*(stemp.livetime/stemp.empty_livetime)*(1-stemp.sample_thickness/stemp.solvent_thickness)^2;
            yrelerr=sqrt(yerr)./yave;
            warning off MATLAB:divideByZero
            ysampleerr2dc=(ysampleave.*ysamplerelerr2).^2+(stemp.livetime/stemp.darkcurrent_livetime).^2*(ydarkcurrentave.*ydarkcurrentrelerr).^2;
            ysolventerr2dc=(ysolventave.*ysolventrelerr2).^2+(stemp.solvent_livetime/stemp.darkcurrent_livetime).^2*(ydarkcurrentave.*ydarkcurrentrelerr).^2;
            yemptyerr2dc=(yemptyave.*yemptyrelerr2).^2+(stemp.empty_livetime/stemp.darkcurrent_livetime).^2*(ydarkcurrentave.*ydarkcurrentrelerr).^2;
            yerr2=ysampleerr2dc/stemp.transfact^2+...
                ysolventerr2dc/stemp.solvent_transfact^2*(stemp.livetime/stemp.solvent_livetime)*(stemp.sample_thickness/stemp.solvent_thickness)^2+...
                yemptyerr2dc/stemp.empty_transfact^2*(stemp.livetime/stemp.empty_livetime)*(1-stemp.sample_thickness/stemp.solvent_thickness)^2;
            yrelerr2=sqrt(yerr2)./yave;
        end
        
        if str2num(si.reduction_state(1))==0 % if sample thickness is to besed
            stemp.sample_thickness=1;
        end
        
        if str2num(si.reduction_state(2))==0 % if absolute intensity scaling factor is to be used
            stemp.abs_int_fact=1;
        end
        
        % now make an intensity correction to bring the units into
        % intensity per solid angle (mr^2)
        if isempty(si.wavelength)
            solid_angle=1;
        else
            solid_angle=(si.kcal^2/si.pixelcal(1)/si.pixelcal(2))*(si.wavelength/2/pi)^2*...
                cos(2*asin(x'*si.wavelength/4/pi)).^3*1E6;
        end
        
        yave=(yave/stemp.sample_thickness*stemp.abs_int_fact)./solid_angle;
        
        switch so.type
            case 'xy'
                so.type = 'x-ave';
                event = 'average vs x';
            case 'rth'
                so.type = 'r-ave';
                event = 'average vs r';
        end
        so.pixels=yave';
        so.relerr=yrelerr';
        so.xaxisfull=x;
    case 'y'
        % we use a different approach here due to problems with discretization 
        % issues at odd angles 
        % here we use the pol2cart method of the original SAXSGUI versions.
        % 
        if str2num(si.reduction_state(12))==0 % if the beamstopmask is not to be used
            stemp.bsmask_pixels=1+0.*si.raw_pixels;
        end
        
        %For mask
        if str2num(si.reduction_state(8))==0 % if the mask is not to be used
            stemp.mask_pixels=1+0.*si.raw_pixels;
        end
        
        % this is kept as a reminder that pixel by pixel flatfield division could be used
        stemp.flatfield_pixels=1+0*si.raw_pixels;   %...but is not
        flatmask=(stemp.bsmask_pixels.*stemp.mask_pixels)./stemp.flatfield_pixels;
        
        %we calculate the intensity average versur azimuthal angel
        
        st=si; %transferring the actual data into temporary saxs structure
        st.pixels=prepare_pixels(st.raw_pixels,flatmask,si.reduction_state).*flatmask; %transferring the 2D raw data into the temporary saxs strucutre
        st=calc_azim_int(st,bounds,npts); %calculating the average vs azimuthal angle
        
        
        %this following calculation is valid for all following
        %averages...but needs to be done at least once.
        len = length(st.pixels);
        x = st.xaxis;
        y = st.yaxis;
        x = [x(1) : (x(2) - x(1)) / (len - 1) : x(2)];
        y = [y(1) : (y(2) - y(1)) / (len - 1) : y(2)];
        
        % calibrated
        if stemp.kcal
            UnitsPerPixel=stemp.kcal/mean(st.pixelcal);
        else
            UnitsPerPixel=1;
        end
        
        dx=x(2)-x(1);
        dy=y(2)-y(1);
        xmin=x(1);
        ymin=y(1);
        xmax=x(len);
        ymax=y(len);
        
        NumOfPix=pi*(xmax^2-xmin^2)/UnitsPerPixel^2*(dy)/360;
        
        % now we have what we need for calculating the average and errors
        % on the actual data.
        ysampleave=st.pixels;
        ysamplerelerr=1./sqrt(NumOfPix.*st.pixels);
        
     
        if (str2num(si.reduction_state(5))==0 || isempty(stemp.transfact))% if sample transmission is not to be used
            stemp.transfact=1;
        end

        
        if str2num(si.reduction_state(4))==1 % if empty holder subtraction is to be used.
            empty_corr=1;
            st=si; %transferring the actual data into temporary saxs structure
            st.pixels=prepare_pixels(st.empty_pixels,flatmask,si.reduction_state).*flatmask; %transferring the 2D empty-holder data into the temporary saxs strucutre
            st=calc_azim_int(st,bounds,npts); %calculating the average vs azimuthal angle
            yemptyave=st.pixels;
            yemptyrelerr=1./sqrt(NumOfPix.*st.pixels)*error_norm_fact;
        else
            empty_corr=0;
            yemptyave=0;
            yemptyrelerr=0; % this does not take into account the error on the transmission number
            stemp.empty_transfact=1; 
            stemp.empty_livetime=1;
        end
        
        if str2num(si.reduction_state(11))==1 % if solvent holder subtraction is to be used.
            st=si; %transferring the actual data into temporary saxs structure
            st.pixels=prepare_pixels(st.solvent_pixels,flatmask,si.reduction_state).*flatmask; %transferring the 2D empty-holder data into the temporary saxs strucutre
            st=calc_azim_int(st,bounds,npts); %calculating the average vs azimuthal angle
            ysolventave=st.pixels;
            ysolventrelerr=1./sqrt(NumOfPix.*st.pixels)*error_norm_fact;
        else
            ysolventave=0;
            ysolventrelerr=0; % this does not take into account the error on the transmission number
            stemp.solvent_transfact=1; 
            stemp.solvent_livetime=1;
            stemp.solvent_thickness=1;
        end

        if str2num(si.reduction_state(3))==1 % if darkcurrent correction average is to be used.
            st=si; %transferring the actual data into temporary saxs structure
            st.pixels=prepare_pixels(st.darkcurrent_pixels,flatmask,si.reduction_state).*flatmask; %transferring the 2D empty-holder data into the temperature saxs strucutre
            st=calc_azim_int(st,bounds,npts); %calculating the average vs azimuthal angle
            ydarkcurrentave=st.pixels;
            ydarkcurrentrelerr=1./sqrt(NumOfPix.*st.pixels)*error_norm_fact;
        else
            ydarkcurrentave=0;
            ydarkcurrentrelerr=0; % this does not take into account the error on the transmission number
            stemp.darkcurrent_livetime=1;
        end
        
        if   ~isempty(si.flatfield_filename)& (str2num(si.reduction_state(6))==1 | str2num(si.reduction_state(7))==1) % if the flatfield is to be used 
            
            %maxflatfield=mean(mean(si.flatfield_pixels(floor(si.raw_center(1)-50):floor(si.raw_center(1)+50),floor(si.raw_center(2)-50):floor(si.raw_center(2)+50))));
            maxflatfield=ffvalue(si,3) %maxflatfield is actually not the maximum but rather an "average" defined in ffnormalize.m
            flatmask=flatmask/maxflatfield;
            st=si; %transferring the actual data into temporary saxs structure
            st.pixels=prepare_pixels(st.flatfield_pixels,flatmask,si.reduction_state).*flatmask; %transferring the 2D empty-holder data into the temperature saxs strucutre
            st=calc_azim_int(st,bounds,npts); %calculating the average vs azimuthal angle
            xflatfield=y;
            yflatfieldave=st.pixels;
            yflatfieldrelerr=1./sqrt(NumOfPix.*st.pixels)*error_norm_fact;
            if str2num(si.reduction_state(7))==1      %if flatfield correction average is to be used
                %here we try to fit in the average plot to a 2nd order
                %polynomial
                P=polyfit(xflatfield(find(~isnan(yflatfieldave))),yflatfieldave(find(~isnan(yflatfieldave)))',2); % fit to a seconorder polynomial
                yflatfieldold=yflatfieldave;
                yflatfieldave=polyval(P,xflatfield)';
            end
            ysampleave=ysampleave./yflatfieldave;
            yemptyave=yemptyave./yflatfieldave;
            ydarkcurrentave=ydarkcurrentave./yflatfieldave;
            ysamplerelerr=ysamplerelerr./yflatfieldave;
            yemptyrelerr=yemptyrelerr./yflatfieldave;
            ysolventrelerr=ysolventrelerr./yflatfieldave;
            ydarkcurrentrelerr=ydarkcurrentrelerr./yflatfieldave;
        end
        % no considerations are taken as to have the flatfield division
        % affects the errorbars...this should be fixed at a later stage.
        
        % OK so now we have everything needed to do the correct
        % subtractions: yave is the average corrected value per pixel in
        % the whole scan...to find the value per pixel per second one has
        % to divide by livetime.
        if isempty(si.solvent_filename)
            warning off MATLAB:divideByZero
            yave=(ysampleave-(stemp.livetime/stemp.darkcurrent_livetime)*ydarkcurrentave)/stemp.transfact-...
                empty_corr*(yemptyave-(stemp.empty_livetime/stemp.darkcurrent_livetime)*ydarkcurrentave)*(stemp.livetime/stemp.empty_livetime/stemp.empty_transfact);
            warning off MATLAB:divideByZero
            yrelerr=sqrt(((ysampleave.*ysamplerelerr).^2+(stemp.livetime/stemp.darkcurrent_livetime).^2*(ydarkcurrentave.*ydarkcurrentrelerr).^2)/stemp.transfact.^2+...
                empty_corr*((yemptyave.*yemptyrelerr).^2+(stemp.empty_livetime/stemp.darkcurrent_livetime).^2*(ydarkcurrentave.*ydarkcurrentrelerr).^2)*...
                (stemp.livetime/stemp.empty_livetime/stemp.empty_transfact).^2)./yave;
            warning off MATLAB:divideByZero
        else

            warning off MATLAB:divideByZero
            %subtracting darkcurrent from each one
            ysampleavedc=(ysampleave-(stemp.livetime/stemp.darkcurrent_livetime)*ydarkcurrentave);
            ysolventavedc=(ysolventave-(stemp.solvent_livetime/stemp.darkcurrent_livetime)*ydarkcurrentave);
            yemptyavedc=(yemptyave-(stemp.empty_livetime/stemp.darkcurrent_livetime)*ydarkcurrentave);
            yave=ysampleavedc/stemp.transfact-...
                ysolventavedc/stemp.solvent_transfact*(stemp.livetime/stemp.solvent_livetime)*(stemp.sample_thickness/stemp.solvent_thickness)-...
                yemptyavedc/stemp.empty_transfact*(stemp.livetime/stemp.empty_livetime)*(1-stemp.sample_thickness/stemp.solvent_thickness);
            warning off MATLAB:divideByZero
            ysampleerrdc=(ysampleave.*ysamplerelerr).^2+(stemp.livetime/stemp.darkcurrent_livetime).^2*(ydarkcurrentave.*ydarkcurrentrelerr).^2;
            ysolventerrdc=(ysolventave.*ysolventrelerr).^2+(stemp.solvent_livetime/stemp.darkcurrent_livetime).^2*(ydarkcurrentave.*ydarkcurrentrelerr).^2;
            yemptyerrdc=(yemptyave.*yemptyrelerr).^2+(stemp.empty_livetime/stemp.darkcurrent_livetime).^2*(ydarkcurrentave.*ydarkcurrentrelerr).^2;
            yerr=ysampleerrdc/stemp.transfact^2+...
                ysolventerrdc/stemp.solvent_transfact^2*(stemp.livetime/stemp.solvent_livetime)*(stemp.sample_thickness/stemp.solvent_thickness)^2+...
                yemptyerrdc/stemp.empty_transfact^2*(stemp.livetime/stemp.empty_livetime)*(1-stemp.sample_thickness/stemp.solvent_thickness)^2;
            yrelerr=sqrt(yerr)./yave;
        end
        
        if str2num(si.reduction_state(1))==0 % if sample thickness is to be used
            stemp.sample_thickness=1;
        end
        
        if str2num(si.reduction_state(2))==0 % if absolute intensity scaling factor is to be used
            stemp.abs_int_fact=1;
        end
        
        yave=yave/stemp.sample_thickness*stemp.abs_int_fact;
        
        so.type = 'th-ave';
        event = 'average vs theta';
        so.pixels=yave;
        so.relerr=yrelerr;
        so.xaxisfull=y;       
end

warning on MATLAB:divideByZero
so.history = strvcat(so.history, event);
so.xaxis=[bounds(1) bounds(1)+bounds(3)];
so.yaxis=[bounds(2) bounds(2)+bounds(4)];

function pixelsout=prepare_pixels(pixelsin,flatmask,reduction_state)
% this routine makes sure that the data to be operated on is corrected
% prior to later binning

if str2num(reduction_state(10))==1    %this means that data should be "zinger removed"
    pixelsin=zinger_removal(pixelsin,flatmask);
end
pixelsout=pixelsin;






