function growthrates = growthParameters2GrowthRates(gp, deltaSt)
% function growthrates = growthParameters2GrowthRates(gp, deltaSt)
%
% A function thats takes a structure gp with fields, s_max, s_min and
% theta. Each of the fields is a Nx1 vector where N denotes the number of
% regions for that organ. So for an organ with 12 regions we would have 
% 
% gp = 
%    s_max: [12x1 double]
%    s_min: [12x1 double]
%    theta: [12x1 double]
%
% The other input is T, this denotes the number of hours it takes for the
% organ to double in size (default T = 20).
%
% Returns the doubling times and the anisotropy rates for each region.
%
% See also: getDoublingTimeColours, getAnisotropyColours
%
% Dr. A. I. Hanna (2007)
if nargin<2
    deltaSt = 20;
end

num_regions = size(gp,3);
growthrates(num_regions) =initGrowthRates;


for i=1:num_regions
    growth_tensor = gp(:,:,i);
    [Smaj, Smin, theta] = gtlib_growthTensor2Params(growth_tensor);
    gr = initGrowthRates;
    
    gr.Smaj = Smaj;
    gr.Smin = Smin;
    gr.Sarea = Smaj + Smin;
    
    gr.Rmaj = log(Smaj)/(deltaSt*log(2));
    gr.Rmin = log(Smin)/(deltaSt*log(2));
    gr.Rarea = gr.Rmaj + gr.Rmin;

    gr.Kmaj = log(Smaj)/(deltaSt);
    gr.Kmin = log(Smin)/(deltaSt);
    gr.Karea = gr.Kmaj + gr.Kmin;
    
    kmin = gr.Kmin;
    kmin(kmin==0) = realmin;
    gr.Kanis = gr.Kmaj/kmin;
    
    % check for zero growth between stages and set the doubling time to
    % realmax
    if gr.Rmin==0
        gr.Dmin = realmax;
    else
        gr.Dmin = 1/gr.Rmin;
    end
    if gr.Rmaj==0
        gr.Dmaj = realmax;
    else
        gr.Dmaj = 1/gr.Rmaj;
    end
    if gr.Rarea==0
        gr.Darea = realmax;
    else
        gr.Darea = 1/gr.Rarea;
    end


    gr.theta = theta;
    growthrates(i) = gr;
end
return;
%%
function gr = initGrowthRates
gr.Smaj = 0;
gr.Smin = 0;
gr.Sarea = 0;

gr.Kmaj = 0;
gr.Kmin = 0;
gr.Karea = 0;
gr.Kanis = 0;

gr.Rmaj = 0;
gr.Rmin = 0;
gr.Rarea = 0;

gr.Dmaj = 0;
gr.Dmin = 0;
gr.Darea = 0;

gr.theta = 0;

return