function [Efrac,Etot,xloc,yloc]=fracten2(Nodes,phi,numdivx,numdivy,imode,nu)
%
%  Name:  fracten2.m
%
%  Usage:  
%
%  This m-file calculated the fractional energies for each mode
%
%   input:	
%  	numdiv=number of subdivisions for the beam
%  	
%
%   output:
%	
%	Etot(i)=total strain energy for mode i
% 	Efrac(i,j)=energy of the jth element of the ith mode
%
%  Version PJC960531
%
%  Author(s): Phillip Cornwell, cornwell@lanl.gov
%
%  Institution: Los Alamos National Laboratory,
%               ESA-EA
%
%  Copyright (c) 1996, The Regents of the University
%       of California. (See top level of code for
%       full copyright statement.)
%

%
% Define the number of subdivisions for the numerical integration
%

nn=3;
x=Nodes(:,2);
y=Nodes(:,3);
z=phi(:,imode);
dx=(max(x)-min(x))/(numdivx*nn);
dy=(max(y)-min(y))/(numdivy*nn);
xi=min(x):dx:max(x);
yi=(min(y):dy:max(y))';

[XI,YI,ZI]=griddata(x,y,z,xi,yi);

xloc=min(x):dx*nn:max(x);
yloc=min(y):dy*nn:max(y);

%
% find derivatives
%
[px,py]=gradient(ZI,dx,dy);
[pxx,pxy]=gradient(px,dx,dy);
[pyx,pyy]=gradient(py,dx,dy);

%
% Determine the integrand
%
fx=pxx.^2+pyy.^2+2*nu*pxx.*pyy+2*(1-nu)*pxy.^2;

%
%  Do the integration
%
for i=1:numdivx
   for j=1:numdivy
	xvalues=xi((i-1)*nn+1:i*nn+1);
	yvalues=yi((j-1)*nn+1:j*nn+1);
	zvalues=fx((j-1)*nn+1:j*nn+1,(i-1)*nn+1:i*nn+1);
        for k=1:nn+1
	   temp(k)=trapz(yvalues,zvalues(:,k));
	end
	Efrac(j,i)=trapz(xvalues,temp);
   end
end

Etot=sum(sum(Efrac));
