% Gyration2.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.
% In this case, we allow for Rg^2 and Rg^4 moments, but no higher
% This is a vast improvement on Gyration1.m which only fits the first moment


function [Io,Rg,Rg4] = Gyration2(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,3) = 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);
    A(i-min_i+1,3) = X(i)*X(i)*X(i)*X(i);
    B(i-min_i+1) = 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) = E(1) + E(2)*X(i)*X(i) + E(3)*(X(i))^4.0;
    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 = E(1)
Rg = (-3.0*E(2)/E(1))^0.5;
Rg4 = (18.0*E(3))^0.25;