% Gyration3.m just calculates the radius of gyration of a small-angle
% x-ray scattering curve.
% The input should be a two column matrix of q's and scattering strengths.
% For a gyration calculation, early data will be unreasonable because of the 
% beamstop while data with a high-q will no longer fit the expansion.
% The fitting function is I(q) = I(0)*exp(-q^2 Rg^2/3)
% This is the guinier approximation

function [Io,Rg] = Gyration3(data)

%  Find the beginning and the end of the region for a Guiner fit!
X = [data(:,1)];
Y = [data(:,2)];
plot(X,Y);
hold ;
min_q = input('What is the lowest value of Q to start fit at?');
max_q = input('What is the maximum value of Q to terminate fit at?');
min_i = 1;
while(X(min_i)<min_q) min_i=min_i+1; end
max_i = min_i + 1;
while(X(max_i)<max_q) max_i= max_i+1; end

%  Do the fit
A(max_i-min_i+1,2) = 0.0;
B(max_i-min_i+1) = 0.0;
for i=min_i: max_i,
    A(i-min_i+1,1) = 1;
    A(i-min_i+1,2) = X(i)*X(i); 
    B(i-min_i+1) = log(Y(i));
end
C = A'*A;
D = A'*B';
E = inv(C)* D;

% And plot it
maximum_y = max(Y)*1.25
Z(size(X)) = 0.0;
for i=1:size(X),
    Z(i) = exp(E(1) + E(2)*X(i)*X(i));
    if (Z(i)<0) Z(i) = 0.0 ; end
    if (Z(i)>maximum_y) Z(i) = maximum_y ; end
end
plot(X,Z,'c');
hold off ;


Io = exp(E(1));
Rg = (-3.0*E(2))^0.5;