function [pxxmean,pxx1sig,pyymean,pyy1sig,xm,ym] = md2curv(phi,Nodes,phi1sig,nmc)
%
%
%   md2curv
%
%   Finds error bounds on curvature of mode shapes
%	based on error bounds on mode shapes
%	using a Monte Carlo approach
%
%	Caution: geometry customized for i40 geometry
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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.                                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


phi = phi(1:26,:);

nresp = size(phi,1);
nmod = size(phi,2);

pxx_phi = zeros(nresp * nmod, nmc);
pxy_phi = zeros(nresp * nmod, nmc);
pyy_phi = zeros(nresp * nmod, nmc);

pxxmean = zeros(nresp,nmod);
pxx1sig = zeros(nresp,nmod);
pyymean = zeros(nresp,nmod);
pyy1sig = zeros(nresp,nmod);

%
%  Compute the variables for the gradient computation
%

znode = reshape([1:26]',13,2)';

xm = zeros(size(znode,2),1);
ym = zeros(size(znode,1),1);

for ty = 1:size(znode,1),

	nr = find(Nodes(:,1) == znode(ty,1));
 	ym(ty) = Nodes(nr(1),3);
end


for tx = 1:size(znode,2),

        nr = find(Nodes(:,1) == znode(1,tx));
 	xm(tx) = Nodes(nr(1),2);
       
end

%  Do for each mode
%


for i = 1:nmod,


    disp(['Processing Mode #',num2str(i)])

%
%   Do nmc MC runs
%

    for j = 1:nmc,


	if j == 1,
		mcphi = phi(:,i);
	else 
		mcphi = phi(:,i) + phi1sig(:,i).*randn(size(phi1sig,1),1);
	end

	zphi = reshape(mcphi,13,2)';

	[px,py] = gradient(zphi,xm,ym);
	[pxx,pxy] = gradient(px,xm,ym);
	[pyx,pyy] = gradient(py,xm,ym);

	pxx_phi( (i-1)*nresp+1 : i*nresp , j) = reshape(pxx',26,1);
	pxy_phi( (i-1)*nresp+1 : i*nresp , j) = reshape(pxy',26,1);
	pyy_phi( (i-1)*nresp+1 : i*nresp , j) = reshape(pyy',26,1);


    end

end
  
   pxxmean = reshape(mean(pxx_phi')',26,nmod);
   pxx1sig = reshape(std(pxx_phi')',26,nmod);
   pyymean = reshape(mean(pyy_phi')',26,nmod);
   pyy1sig = reshape(std(pyy_phi')',26,nmod);


return

