% 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, Io , std_Io] = gyration_l(data)

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

end

plot(X,Y);
disp('You can now play with the graph until you get bored.');
disp('The graph data is in X and Y.');
disp('You probably just want to adjust the range a bit.');
disp('When you are happy, type "return".');

keyboard ;
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) = data(i,1) * data(i,2);
end
C = C';
[B , correl_B, phi] = wlin(A,C,w);


% And plot it
Z = B(1) + B(2)*X;
plot(X,Z,'c');

err_Y(size(Y)) = 0.0;
for i=1:size(Y)
     err_Y(i) = (phi/w(i))^0.5;
end
phi
Z = Y + err_Y;
plot(X,Z,'r');
Z = Y - err_Y ;
plot(X,Z,'g');


hold off ;

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