% function [r, p, q, I] = Compute_SAXS(positions, charges)
%
%      Computes the scattering pattern for a molecule.
%
%                   positions(j,:) = [x,y,z] coordinates of j-th scatterer
%                   in molecule.
%                   charges(j) = number of charges/scatterer strength of
%                   j-th scatterer in molecule.
%
%       Output information
%           Pair Correlation Function - p(r) --->  p(j) versus r(j).
%           Scattering Pattern - I(q) ---> I(j) versus q(j)
%

function [r, p, q, I] = Compute_SAXS(positions, charges);

R = positions'; % R(:,j) are the 3-d coordinates of the j-th scattering object.
charge =charges'; % charge(j) is the charge associated with the j-th scattering object.

% Make a histogram of pair-pair distances weighted by charge..

%  Set up the bins for the histogram.
r_min = 0; % Minimum allowed pair correlation length.
r_max = 2*sqrt(3)*max(max(R)) ; % Maximum pair correlation length.
Nbins = 300; % Number of bins in our histogram.
r = ([1:Nbins]-0.5)/Nbins*(r_max-r_min)+r_min;
g = zeros(Nbins,1);

% Treat the scatterers as gaussian balls with a size "sigma" Angstroms.
sigma = 1.0;
sf = -1.0/(4*sigma*sigma);

% Sum all atom-atom pairs.
for j=1:length(charge)
    for k=1:length(charge);

        % Compute distance and weight of this pair of scatterers
        dist = norm(R(:,j)-R(:,k));
        wt = charge(j)*charge(k);

        % Distribute this pair into the histogram bins.
        if (dist>0.1*sigma)
            wt = wt*sigma*sigma;
            for m=1:Nbins
                 d1=r(m)-dist; d2 = r(m)+dist;
                 g(m)=g(m)+wt * (exp(sf*d1*d1) - exp(sf*d2*d2)) * ( r(m)/dist ) ;     
            end
        else
            for m=1:Nbins
                 g(m)=g(m)+wt * exp(sf*r(m)*r(m)) * r(m)* r(m) ;     
             end 
        end
    
    end
end
 
% Normalize the pair correlation function to one.
g = g/ sum(g);

% When we return the pair correlation function we wish to normalize it so 
% Integral(g(r).dr)=1
p = g / (r(2)-r(1));

% Now convolve g(r) with sinc(qr) to determine I as a function of q.
% Initialize the q-range
q_min = 0 ; 
Rg = sqrt(sum(g.* r'.*r')/2); % Natural size of protein 
q_max = 8/Rg;
Nqbins = 500;

q = ([1:Nqbins]-1)/Nqbins * (q_max-q_min) + q_min;
I = zeros(size(q));

for j=1:Nqbins
    for k=1:Nbins
        I(j) = I(j) + g(k) * sinc( q(j) * r(k)/pi);
    end
end