function ssi_out = regrid_ssi(ssi_in,wvl_edges_in,wvl_edges_out)

% Tim Kruschke
% mail: tkruschke[at]geomar.de
%
% Function to regrid spectral solar irradiance. Prepared in context of
% ROMIC-SOLIC project and CMIP6 solar forcing, where SATIRE and NRL data
% are considered.
% 
% Input:
% ssi_in: spectral solar irradiance (dimensions wavelengths x time)
% wvl_edges_in: edges of original wavelength bins (dimensions wavelengths x 2)
% wvl_edges_out: edges of target wavelength bins (dimensions wavelengths x 2)
%
% Output:
% ssi_out: spectral solar irradiance fro specified target wavelength bins (bin means)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ssi_out=zeros(size(wvl_edges_out,1),size(ssi_in,2));

for i=1:size(wvl_edges_out,1);
    index_in=find(wvl_edges_in(:,2)>wvl_edges_out(i,1) & wvl_edges_in(:,1)<wvl_edges_out(i,2));
    if ~isempty(index_in);
        if wvl_edges_in(index_in(1),1)>wvl_edges_out(i,1); % Linear interpolation for part between low edge of target bin
            tmp_ssi=ssi_in(index_in(1)-1,:)+ ...           %  and first original bin found inside
                (ssi_in(index_in(1),:)-ssi_in(index_in(1)-1,:))/ ...
                (wvl_edges_in(index_in(1),1)-wvl_edges_in(index_in(1)-1,2))* ...
                ((wvl_edges_in(index_in(1),1)+wvl_edges_out(i,1))/2-wvl_edges_in(index_in(1)-1,2));
            ssi_out(i,:)=ssi_out(i,:)+tmp_ssi*(wvl_edges_in(index_in(1),1)-wvl_edges_out(i,1));
        end
        for ii=1:length(index_in);
            lowmarg=(wvl_edges_in(index_in(ii),1)<wvl_edges_out(i,1));
            upmarg=(wvl_edges_in(index_in(ii),2)>wvl_edges_out(i,2));
            if ~lowmarg && ~upmarg;
                ssi_out(i,:)=ssi_out(i,:)+ssi_in(index_in(ii),:)* ...
                    (wvl_edges_in(index_in(ii),2)-wvl_edges_in(index_in(ii),1));
            elseif lowmarg && upmarg;
                ssi_out(i,:)=ssi_out(i,:)+ssi_in(index_in(ii),:)*(wvl_edges_out(i,2)-wvl_edges_out(i,1));
            elseif lowmarg
                ssi_out(i,:)=ssi_out(i,:)+ssi_in(index_in(ii),:)* ...
                    (wvl_edges_in(index_in(ii),2)-wvl_edges_out(i,1));
            elseif upmarg
                ssi_out(i,:)=ssi_out(i,:)+ssi_in(index_in(ii),:)* ...
                    (wvl_edges_out(i,2)-wvl_edges_in(index_in(ii),1));
                            
            else
                error('Logical error in function!')                
            end
            if ii~=length(index_in);
                if wvl_edges_in(index_in(ii),2)~=wvl_edges_in(index_in(ii)+1,1); % Linear interpolation for gap between original bins
                    tmp_ssi=(ssi_in(index_in(ii),:)+ssi_in(index_in(ii+1),:))/2;
                    ssi_out(i,:)=ssi_out(i,:)+tmp_ssi* ...
                        (wvl_edges_in(index_in(ii+1),1)-wvl_edges_in(index_in(ii),2));
                end
            end
        end
        if wvl_edges_in(index_in(end),2)<wvl_edges_out(i,2); % Linear interpolation for part between high edge of target bin
            tmp_ssi=ssi_in(index_in(end),:)+ ...             %  and last original bin found inside
                (ssi_in(index_in(end)+1,:)-ssi_in(index_in(end),:))/ ...
                (wvl_edges_in(index_in(end)+1,1)-wvl_edges_in(index_in(end),2))* ...
                ((wvl_edges_in(index_in(end),2)+wvl_edges_out(i,2))/2-wvl_edges_in(index_in(end),2));
            ssi_out(i,:)=ssi_out(i,:)+tmp_ssi*(wvl_edges_out(i,2)-wvl_edges_in(index_in(end),2));
        end
    else % Linear interpolation between two original bins onto center of empty target bin
        index_below=find(wvl_edges_in(:,2)<wvl_edges_out(i,1),1,'last');
        index_above=find(wvl_edges_in(:,1)>wvl_edges_out(i,2),1,'first');
        tmp_ssi=ssi_in(index_below,:)+ ...
            (ssi_in(index_above,:)-ssi_in(index_below,:))/ ...
            (wvl_edges_in(index_above,1)-wvl_edges_in(index_below,2))* ...
            ((wvl_edges_out(i,1)+wvl_edges_out(i,2))/2-wvl_edges_in(index_below,2));
        ssi_out(i,:)=tmp_ssi*(wvl_edges_out(i,2)-wvl_edges_out(i,1));
    end
    ssi_out(i,:)=ssi_out(i,:)/(wvl_edges_out(i,2)-wvl_edges_out(i,1));
end

return