function [xfp,yfp,DI]=findDI2D(Nodes,phiundam,phidam,numdivx,numdivy,mode_anal,nu)
%
%  Name: findDI2D.m
%
%  Usage:  [xfp,yfp,DI]=finDI2D(Nodes,phiundam,phidam,numdivx,numdivy,mode_anal,nu)
%
%  This m-file calculated the fractional energies for each mode
%
%   input:	
%  	numdiv=number of subdivisions for the beam
%  	
%
%   output:
%	
%
%  Version SWD970115
%
%  Author(s): Phillip Cornwell, cornwell@lanl.gov
%	modified by Scott W. Doebling, doebling@lanl.gov
%
%  Institution: Los Alamos National Laboratory,
%               ESA-EA
%
%  Copyright (c) 1996-1997, 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
%

betanum=zeros(numdivy,numdivx);
betaden=zeros(numdivy,numdivx);

fprintf(1,'Calculating damage index using 2-D technique...');
  
for i=mode_anal,
  [Efrac,Etot,xloc,yloc]=fracten2(Nodes,phiundam,numdivx,numdivy,i,nu);
  [Efracd,Etotd,xloc,yloc]=fracten2(Nodes,phidam,numdivx,numdivy,i,nu);

  betanum=betanum+Efracd/Etotd;
  betaden=betaden+Efrac/Etot;
end

beta=betanum./betaden;

%
%  check for NAN's
%
nanloc=find(isnan(beta)==1);

if ~isempty(nanloc),
  beta(nanloc)=ones(nanloc);
end
%
%  determine the normalized values of the damage index
%
avebeta=mean(beta(:));
stdbeta=std(beta(:));
DI=(beta-avebeta)/stdbeta;

%
%  determine the x and y location for plotting
%
xvalues=(xloc(1:numdivx)+xloc(2:numdivx+1))/2;
yvalues=(yloc(1:numdivy)+yloc(2:numdivy+1))/2;

xfp=ones(length(yvalues),1)*xvalues;
yfp=yvalues'*ones(1,length(xvalues));


fprintf(1,'done\n')


return