function [results,h]=biftcalc(aname,x,y,sd,sigy,NoGui,biftfig,h,nmax,alphamin,alphamax,nalpha,ndiam,dmax0,dstep,ntot0,mtot,searchparam,verboseparam)
% C***********************************************************************
% 
% results=biftcalc(x,y,sd,NoGui,biftfig,h,nmax,nalpha,ndiam,dmax0,ntot0,searchparam)
%
% 
%
% This program is part of a small squite of programs that calculate the IFT of SAS data
% using a Bayesian a approach developed
% by Dr. S. Hansen, of the Danish Agricultural University.
%
% This part of the program does the actual grid-evaluation (over alpha and
% dmax) and then calculates the various resultant parameters with errorbars
%
% Details in Hansen, S. (2000) J.Appl.Chryst 33, 1415-1421
%
% It is a matlab revision of a fortran program provided kindly by 
% Dr. S. Hansen.
%
% C***********************************************************************

format40='Calc   Dmax  Alpha     Chi-sq.   Constr.   Check    Sum     log(Evidence)\r';
alpharel=exp(-log(alphamax/alphamin)/(nalpha-1));
if isnan(alpharel)
    alpharel=1;
end
nof=0;
%pre-dimensioning for speed
ftot=zeros(nmax-1,nalpha*ceil(ndiam));
for icount=1:ndiam
    dmax=dmax0+dstep*(icount-1);
    ntot=ceil((dmax/dmax0)*ntot0);
    dfx=dmax/ntot;
    
    [M,xf]=bift_prior(ntot,y,dfx);
    [A]=bift_trans(x,xf);
    
    ysum=(y./sd)*A;
    
    onetemp=1+0*(1:ntot);
    SDMAT=sd'*onetemp;
    B=A'*(A./SDMAT);
    clear SDMAT 
    clear onetemp
    
    %This calculates the special determinant matrix...outside alpha-loop
    %for speed
    Bsize=max(size(B));
    %creating matrix with "next-to-diagonal"=1 and diagonal=0
    uoffset=circshift(eye(Bsize,Bsize),[1,0])+circshift(eye(size(B)),[0,1]);
    uoffset(1,Bsize)=0;
    uoffset(Bsize,1)=0;
    
    %this calculates the intensity close to 0 for normalization
    fm=(A*M')';
    summ=sum((y(1:3)./sd(1:3)))/sum(fm(1:3)./sd(1:3));
    
    M=summ.*M;
    fprior=1.0005*M;
    
    sigma=(1:ntot)*0+1;
    evmax=-1.e10;
    searchparam.omega=searchparam.omegaorig;
    alpha=alphamax;
    wgrads=1;
    %disp (['start variation of alpha',num2str(alpha)])
    fprintf(1,format40)
    
    % c**************************************************************
    % c     Start variation of alpha
    % c**************************************************************
    for alphacount=1:nalpha
        
        nof=nof+1;
        if mod(nof,10)==1 && ~NoGui
            delete(h)
            h=figstatus(['Finished ',num2str(nof), ' of ', num2str(nalpha*ndiam), ' calculations.'],biftfig);
        end
        alpha=alphamax*alpharel^(alphacount-1);
        [f,c,s,dotsp,fm]=bift_search(alpha,dmax,fprior,M,A,B,ysum,searchparam,verboseparam,x,y,sd,xf);
        fprior=f;
        %         C************************************************************************
        %         C     Estimate of p(r) written to file
        %         C     and the probability for this particular solution is calculated
        %         C****fv********************************************************************
        %         
        
        alpha=alpha/(2.); % this is necessary because the original program had an error of factor two in the definition of the optimizing function
        evidence=bift_calc_evidence(alpha,s,c,B,uoffset,mtot);
        
        if(evidence > evmax) 
            evmax=evidence;
            chi2=c;
        end 
        
        ftot(1:ntot,nof)=f(1:ntot)'*dfx;
        %             ftot(ntot+1:nmax-4,nof)=0.;
        ftot(nmax-1,nof)=dmax;
        ftot(nmax-2,nof)=c;
        ftot(nmax-3,nof)=alpha;
        
        dot(nof)=dotsp;
        prob(nof)=evidence;
        evimatrix(nof,1:4)=[nof,dmax,alpha*2,evidence];
        
        alpha=alpha*2.;
        format161= ' %4i %6.1f %8.2e %8.2e %8.2e %8.2e %8.2e %9.2e\n';
        fprintf(1,format161, [nof,dmax,alpha,c,-s,dotsp,summ,evidence]);
        dotsp=0;
    end
    
    %     
    %  End of variation of alpha
    % 
    
end  

% c***********************************************
% c     End of variation of dmax (and alpha)
% c***********************************************


evimax=max(prob);
evimatrix(:,4)=evimatrix(:,4)-evimax; %multiplying evidence by a factor (normalizing to max=1)
prob=exp(prob-evimax);
summodel=sum(prob);
totprob=summodel;

% C*****************************************************************************
% C     Output estimate of p(r) - including errors
% C*****************************************************************************
if ~NoGui
    delete(h)
    h=figstatus('Calculating final parameters',biftfig);
end
xmean=ftot(1:ntot,:)*prob'/totprob;
std=ftot(1:ntot,:).^2*prob'-2*(ftot(1:ntot,:)*prob').*xmean+sum(prob)*xmean.^2;
std=sqrt(std/totprob)/dfx;

f=xmean/dfx;
cc=dfx*4*3.14159;
f=f/cc;
sigf=std/cc;

% c******************************************************
% c     Output Rg
% c******************************************************
rmax=inf;
sumtemp=0.;
sum2=0.;
esum=0.;
esum2=0;
sumtemp=sumtemp+.5*f(1);
sumtemp=sumtemp+f(2);
sum2=sum2+.5*f(1)*xf(1).^2;
sum2=sum2+f(2)*xf(2).^2;
esum=esum+.5*sigf(1).^2;
esum=esum+sigf(2).^2;
esum2=esum2+.5*xf(1).^4*sigf(1).^2;
esum2=esum2+xf(2).^4*sigf(2).^2;
for ii=3:ntot
    if(xf(ii) <= rmax) 
        sumtemp=sumtemp+f(ii);
        esum=esum+sigf(ii).^2;
        sum2=sum2+f(ii)*xf(ii).^2;
        esum2=esum2+xf(ii).^4*sigf(ii).^2;
    end
end
const=dfx*4*pi;
rg=sqrt(.5*sum2/sumtemp);

sigmar=0.5*rg*sqrt((esum/sumtemp.^2+esum2/sum2.^2));

% c******************************************************
% c     Output of Diameter
% c******************************************************

diammean=0;
probtot=0;
% dmax is kept in ftot(nmax-1,:)
diammean=ftot(nmax-1,:)*prob'/sum(prob);
sigmad=(diammean-ftot(nmax-1,:)).^2*prob'/sum(prob);
sigmad=sqrt(sigmad);



% 
% c******************************************************
% c     Output fit of data
% c******************************************************

fm=(f'*A'*cc);


format1143=' Qmin     Rg    S_Rg     I(0)         S_I(0)    Dmax   S_Dmax\r';
format1133=('%6.4f %6.1f  %6.2f  %10.3g %10.3g  %6.2f  %6.2f \r');

fprintf(1,format1143)
fprintf(1,format1133,[min(x),rg,sigmar,sumtemp*const,sqrt(esum)*const,diammean,sigmad]);

sizeevimatrix=size(evimatrix);
zmat=reshape(evimatrix,round(sizeevimatrix(1)/ndiam),ndiam,4);
evix=zmat(1,:,2);
eviy=zmat(:,1,3);
eviz=zmat(:,:,4);

results.name=aname;
results.xin=x;
results.yin=y;
results.sigmayin=sigy;
results.yfit=fm;
results.r=xf;
results.pr=f;
results.sigmapr=sigf;
results.evix=evix;
results.eviy=eviy;
results.eviz=eviz;
results.ndmax=ndiam;
results.rg=rg;
results.sigmarg=sigmar;
results.Io=sumtemp*const;
results.sigmaIo=sqrt(esum)*const;
results.DmaxMean=diammean;
results.SigmaDmaxMean=sigmad;
results.DmaxMostLikely=results.evix(floor(find(results.eviz==max(max(results.eviz)))/length(results.evix)));



if ~NoGui
    delete(h)
end
end
