% Gyration.m calculates the radius of gyration of a small-angle
% x-ray scattering curve using the Guinier method.
% The input should be a two column matrix of q's and scattering strengths.
% The data is plotted to allow the user to pick which range of q-values 
% should be fitted.  Clearly the very low Q-data is meaningless because of the
% beam-stop, while the very high Q-data will not fit the Guinier approximation % well.  

% The variance of the data is assumed to go as 1/q, which is correct for 
% independent errors at each scattering angle. 

% The fitting function is I(q) = I(0)*exp(-q^2 Rg^2/3)
% This is the guinier approximation

% The curve-fit is performed by wlin.m.

function [Rg, std_Rg] = gyration_l(data)

%  Find the beginning and the end of the region for a Guiner fit!
X = [data(:,1)];
Y = [data(:,2)];
maximum_Y = max(Y);
for i=1:size(Y)
   if (Y(i)>(0.0002*maximum_Y)) 
   Y(i) = log(Y(i)); 
   else 
   Y(i) = log(0.0002*maximum_Y) ;
   end
   X(i) = X(i)*X(i);

end

plot(X,Y);
hold on ;
min_q = input('What is the lowest value of Q^2 to start fit at?');
max_q = input('What is the maximum value of Q^2 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;
C(max_i-min_i+1) = 0.0;
w(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); 
    C(i-min_i+1) = Y(i);
    w(i-min_i+1) = X(i)^0.5 * Y(i);
end
C = C';
[B , correl_B, phi] = wlin(A,C,w);

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

Rg = (-3.0*B(2))^0.5;
std_Rg = 1.5*sqrt(correl_B(2,2))/Rg ;