function [f,c,s,dotsp,fm]=bift_search(alpha,dmax,fprior,M,A,B,ysum,searchparam,verboseparam,x,y,sd,xf)
% this function performs a search for the optimum solution to the 
% regularisation criteria

dotsp=0;
ite=1;
f=fprior;
ntot=max(size(f));
Beye=B.*eye(ntot);

while (ite<=searchparam.maxit && abs(1-dotsp)>searchparam.dotsptol && ...
        searchparam.omega>searchparam.omegamin) || (ite<=searchparam.minit)
    %)
    %     C*************************************************************************
    %     C     Calculation of p(r) for given diameter and alpha
    %     C*************************************************************************
    %     
    if ite ~=1        
        for jj=2:ntot-1
            M(jj)=(f(jj-1)+f(jj+1))/2 ;
        end
        M(1)=f(2)/2;
        M(ntot)=f(ntot-1)/2;

        %tic
%         fold=f;
%         for ii=1:ntot
%             fsumi=0;
%             for jj=1:ntot
%                 fsumi=fsumi+B(ii,jj)*f(jj);
%             end
%             fsumi=fsumi-B(ii,ii)*f(ii);
%             fx=(alpha*M(ii)+ysum(ii)-fsumi)/(alpha+B(ii,ii));
%             f(ii)=(1.-searchparam.omega)*f(ii)+searchparam.omega*fx;
%         end
        %toc
%         tic
        fsumi=f*(B-Beye);
        fx=(alpha*M+ysum-fsumi)./(alpha+sum(Beye));
        fold=f;
        f=(1-searchparam.omega)*f+searchparam.omega*fx;
%         toc
    end
    
    ite=ite+1;
    dotsp=-1;
    
    s=sum(-(f-M).^2);
    gradsi=(-2*(f-M));
    wgrads=sum(gradsi.^2);
    gradci=(2*(f*B-ysum));
    wgradc=sum(gradci.^2);
    dotsp=sum(gradsi.*gradci);
    while (dotsp<0 && alpha<10*max(sum(Beye)) && ite>2 && searchparam.omega>searchparam.omegamin)
        
        %this should take care of unstable scenarios solutions where the optimization is
        %going too fast
        searchparam.omega=searchparam.omega/searchparam.omegareduction;
        f=(1-searchparam.omega)*fold+searchparam.omega*fx;
        
        s=sum(-(f-M).^2);
        gradsi=(-2*(f-M));
        wgrads=sum(gradsi.^2);
        gradci=(2*(f*B-ysum));
        wgradc=sum(gradci.^2);
        dotsp=sum(gradsi.*gradci);
    end

    %this should take care of unstable scenarios solutions where the optimization is
    %going too fast
    %                      
    fm=f*A';
    c=mean((y-fm).^2./sd);           
    
    if ((wgrads==0) || (wgradc ==0)) 
        dotsp=1
    else
        wgrads=sqrt(wgrads);
        wgradc=sqrt(wgradc);
        dotsp=dotsp/(wgrads*wgradc);
        if rem(ite,verboseparam.viewstep)==0 && verboseparam.on==1
            verbosefig=findobj('Tag','Verbosefigure');
            if isempty(verbosefig)
                verbosefig=figure;
                set(verbosefig,'Tag','Verbosefigure')
                set(verbosefig,'Name','Intermediate results from optimization')
            else
                figure(verbosefig)
            end
            subplot(2,1,1)
            plot(xf,M,xf,f)
            title([' iter=',num2str(ite),' dmax=',num2str(dmax),' alpha=',num2str(alpha),' convergence=',num2str(dotsp),' omega=',num2str(searchparam.omega)])
            xlabel('r in A')
            ylabel('p(r)')
            subplot(2,1,2)
            plot(x,y,x,fm)
            xlabel('q (1/A)')
            ylabel('I (a.u.)')
            drawnow
        end
    end
end
