% sine_fit tries to find a reasonable pair correlation function.
% It does so by making assumptions about the pair correlation function.
% First, it assumes that the correlation function goes to zero for r>D.
% Thus, D must be larger than the size of the protein and if there are protein
% protein correlations contained in the data these will not show up.
% Next, it assumes that within the range 0<r<D, r the function is made up of a % sum of the first N harmonics.  So,
%  P(r)/r = Sum ( sin( n pi / D * r)* a_n )

%  It then attempts a least square's fit of these co-efficients from the data.

% for input, it needs a 2 column data matrix with values of q and scattering
% intensities for those q values.
% Also, it needs to know the proposed maximum correlation length, D, and
% the number of harmonics in the fit, m.

function [P] = sine_fit( data, D, m)

N = size(data,1);

% Set up the fitting functions.
   for j=1:m
            temp1 = j*pi*(-1.0)^(j+1)/D ;
            for i=1:N
            A(i,j) = sin(data(i,1)*D)/( (j*pi/D)^2.0 - data(i,1)^2 );
            A(i,j) = temp1 * A(i,j) / data(i,1);
            end
   end

% Set up the weights
  C = data(:,2); w = zeros(N,1);
  for j=1:N
         w(j) = data(j,1); 
  end

% Do the inversion
  [b , std_b] = wlin(A,C,w);

% What the fit looks like in q-space.
F = A*b ;


% plot(data(:,1),C);
% hold on
% plot(data(:,1),F,'c');
% pause

% What the fit looks like in r-space
P = b;
 hold off
 g(200) = 0.0;
 r(200) = 0.0 ;
 for i=1:200
     r(i) = (i-0.5)/200.0*D ;
     g(i) = 0.0 ;
     for j=1:m
        g(i) = g(i) + b(j)*sin(j*pi*r(i)/D);
     end
     g(i) = g(i)*r(i) ;
 end

hold on 
plot(r,g)