function  matrix_computation(funct, params, X, Y_experimental)
%          SSD matrix computation for 3D function
%
%  3/21/08 - uses function with swapped order of parameters
%                - uses X with column-oriented structure
%
% 3/22/08 - working version after bug fixes
%
% 3/26/08  - now it computes SSD to EXPERIMENTAL DATA rather than to Y DATA
%                   produced with the model parameter set. This tells you
%                   what are other parameter sets that will describe you
%                   experimental data as good as the best-fit result.
% 

global X_column

%============= First plot: data presentation ==================
Y=Y_experimental;

%=================== Generate parameter arrays =======================
SSD_matrix_size=[]; %-- collect all dimensions ---
for i=1:length(params.param)
    param_grid(i).range=params.range(i,:);
    param_grid(i).array = [ param_grid(i).range(1): params.grid_step(i) : param_grid(i).range(2)];
    param_grid(i).length = length( param_grid(i).array );   
    SSD_matrix_size=[ SSD_matrix_size param_grid(i).length ];
end


%================ Scan parameter grid and generate SSD matrix ==============

SSD_matrix=zeros(SSD_matrix_size);  %--- make N-D matrix to store SSD values

current_pars=params; %--- parameters to probe
switch length(params.param)

%=================================== 3D ===========================================        
    case 3
        %--- calculate matrix --
        for i=1:param_grid(1).length
		    disp(sprintf('Computation of plane %d out of %d', i,param_grid(1).length));
            for j=1:param_grid(2).length
                for k=1:param_grid(3).length
                    %-- set new parameters ---
                    p1=feval(params.funct{1},param_grid(1).array(i) );
                    p2=feval(params.funct{2},param_grid(2).array(j) ); 
                    p3=feval(params.funct{3},param_grid(3).array(k) );
                    current_pars.param = [ p1; p2; p3 ];
	
                    %--- calculate deviation of function with new parameters
                    %   from ideal data 
                    Ynew=feval(funct, current_pars, X);
                    %SSD_matrix(i,j,k)= sqrt(sum( (Ynew-Y).^2 ));
                    SSD_matrix(i,j,k)= sum( (Ynew-Y).^2 );
              

                end
            end
        end
        save -mat 'SSD_matrix.mat' SSD_matrix;
            


                 

    otherwise
        error(sprintf('I cannot do %d dimensions!',length(params.param)));
end
        




