function deltak = mrptstmt(K_u,M_u,phi_d,Freq_d)
%
%
% Function: mrptstmt
% 
% Usage:  deltak = mrptstmt(K_u,M_u,phi_d,Freq_d)
% 
% Description:
%    Computes the perturbation to the global stiffness matrix
%
%

% Version SWD980701
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This matlab source code was originally     %
% developed as part of "DIAMOND" at          %
% Los Alamos National Laboratory. It may     %
% be copied, modified, and distributed in    %
% any form, provided:                        %
%  a) This notice accompanies the files and  %
%     appears near the top of all source     %
%     code files.                            %
%  b) No payment or commercial services are  %
%     received in exchange for the code.     %
%                                            %
% Original copyright is reserved by the      %
% Regents of the University of California,   %
% in addition to Scott W. Doebling, Phillip  %
% J. Cornwell, Erik G. Straser, and Charles  %
% R. Farrar.                                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if iscell(Freq_d),
   for i = 1:length(Freq_d),
      Freq_d_mat(i) = Freq_d{i};
      phi_d_mat(:,i) = phi_d{i}.';
   end
   Freq_d = Freq_d_mat;
   phi_d = phi_d_mat;
else
   phi_d = phi_d.';
end

Lambda_d = diag((Freq_d*2*pi).^2);

dvecs = K_u * phi_d - M_u * phi_d * Lambda_d;

H = inv(dvecs.' * phi_d);

deltak = dvecs * H * dvecs.' * 100;

return