
%Jan 15 10 modelSens adapted from lqnSolve

function [h,H]=modelSens(param,controlled,parFit,parInteger)

[htemp] = modelCalc(param, controlled);
%controlled is row vector, rest are columns
h=[controlled,htemp']';  %column vector, pad for controlled variables
nn = size(param); %number of parameters including those not being fitted
%if (size(nn)>1) n = nn(1); else n = nn; end
n = nn(1);
mm = size(h);  %number of results including padding
%if (size(mm)>1) m = mm(1); else m = mm; end
m = mm(1);
%begin sensitivity computation
H = zeros(m,n);

for k = 1:n %derivatives wrto param k
    if(parFit(k)>0)
        if(parInteger(k)>0) %for integer increments
           inc = 1;
        end;
        if(parInteger(k)==0)
            inc = 0.01*param(k);
        end;        
        param(k) = param(k) + inc;
        [hpert] = modelCalc(param, controlled);
        hpert = [controlled,hpert']'; % padded
        param(k) = param(k) - inc;
        for j = 1:m
            H(j,k)=(hpert(j) - h(j))/inc; %includes rows j of zero values for fixed and non-fitted parameters
        end
    end;
end;
%H has empty columns for the parameters not flagged by 1's in selectpar