%   [ r, p, rms_error] = Guess_PCF(data, q_range, D, N, r, p)
%
%   Generate a Pair Correlation Function by the Guess and Check Method
%
%       Inputs -
%                   I versus q -->   data(:,1)= q; data(:,2) = I;
%                   q_range    -->  Restrict q_range(1)<=q<=q_range(2)
%                   D       --->  Maximum correlation length.
%                   N --> Number of segments to use if generating a new
%                   pair correlation function
%                   r,p -->  p(r) --> r(j) and p(j) of an initial guess for
%                   p(r)
%                     note r(1)=0; r(N+1)=D; p(1)=p(N+1)=0!
%
%       Outputs - 
%                   Pair Correlation function p(r) --> r(j) and p(j)
%                   RMS_Error --> Root Mean square error difference between
%                       theoretical I versus Q and measured I versus q.
%

function [r,p,rms_error] = Guess_PCF(data_in, q_range, D, N, r, p)

% Cut the data just to the allowed q_range;
temp = length(data_in); jmin = 1;
while ((jmin<temp)&(data_in(jmin,1)<q_range(1)))
    jmin=jmin+1;
end
jmax=jmin;
while ((jmax<temp)&(data_in(jmax,1)<q_range(2)) )
    jmax = jmax+1;
end
data = data_in([jmin:jmax],:);
Nd = size(data,1);

% Generate the initial guess of p.
if (nargin<5)
r = [0:N]/N * D;
p = [0 , triang(N-1)' , 0];
dp = [ 0 , ones(1,N-1), 0]*0.1;
else
    N = length(r)-1;
    dp = [ 0 , ones(1,N-1), 0]*0.1;
end


resbest = err_fit(data,r,p);

for k=1:30
for j=2:N
    pnew = p;
    pnew(j) = pnew(j) + dp(j)*pnew(j);
    [res, Ic] = err_fit(data, r, pnew);

    if (res<resbest)
        p = pnew;
        resbest=res;
        dp(j) = dp(j)*1.1;
    else
        dp(j)= -0.75*dp(j);
    end
end

    
   fprintf(1,'\n Iteration = %d , RMS Error = %f.', k, resbest);
   rms_error = resbest;
   p = p / (sum(p) * (r(2)-r(1)) );

   figure(1); clf;
   plot(data(:,1), data(:,2),'b'); hold on
   plot(data(:,1), Ic,'r');
   xlabel('q (A^{-1})');
   ylabel('I (Pollacks)');
   figure(2); clf;
   plot(r, p); hold on
   plot(r,p + p.*dp,'r');
   xlabel('r - Angstroms');
   ylabel('p(r) - A^{-1}');
   drawnow
end


function [res,Ic] = err_fit(data, r, p)

    q = data(:,1);
    Io = data(:,2);
    Ic = zeros(size(q));
    N = length(r)-1;
    
    for j=1:N
        u = 0.5*(r(j+1)+r(j));
        Ic = Ic + (0.5*(p(j+1)+p(j)) * (r(j+1)-r(j))/u) * sin(q*u);
    end
   
    temp1= sum(Io.*Io);
    temp2= sum(Io.*Ic);
    temp3 = sum(Ic.*Ic);
    
    res = sqrt( (temp1 - temp2*temp2/temp3) / length(Io) ) / sqrt( temp1/length(Io));
    alpha = temp2/temp3;
    Ic = alpha*Ic;