function beta=findBeta(Etot,Efrac,Etotd,Efracd,mode_anal)
%
%  Name:  findBeta.m
%
%  Usage:  beta=findBeta(Etot,Efrac,Etotd,Efracd,mode_anal)
%
%  Calculates the raw damage index for each subdivision of 1 bar element
%       over all the modes
%
%   input:
%	Etot(i)=int(0,L,Phipp^2)=total strain energy for undamaged ith mode 
%	Efrac(i,j)=energy of the jth element of the undamaged ith mode
%	Etotd(i)=int(0,L,Phipp^2)=total strain energy for damaged ith mode 
%	Efracd(i,j)=energy of the jth element of the damaged ith mode
%
%   output:
%	betanum=((Efracd(1,:)/Etotd(1)+1)+(Efracd(2,:)/Etotd+1) +...)
%	betaden=((Efrac(1,:)/Etot(1)+1)+(Efrac(2,:)/Etot+1) +...)
%  beta = betanum ./ betaden
%
%   notes:
%    	work needs to be done on the best way to combine multiple modes
%	and whether or not to include the +1
%
%  Version SWD970403
%
%  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, The Regents of the University
%       of California. (See top level of code for
%       full copyright statement.)
%

%
%  define the +1 as a variable so it is easy to change
%
ii=1;

betanum=zeros(1,size(Efrac,2));
betaden=zeros(1,size(betanum));
for i=mode_anal
   betanum=betanum+Efracd(i,:)/Etotd(i)+ii;
   betaden=betaden+Efrac(i,:)/Etot(i)+ii;
end
 
beta=betanum ./ betaden;

% alternate approach
%
% beta=zeros(size(Efrac,1))
% for i=mode_anal
%    beta=beta+(Efracd(i,:)/Etotd(i)+ii)./(Efrac(i,:)/Etot(i)+ii)
% end
%

%
% calculate the normalized variable
%
avebeta=mean(beta);
stdbeta=std(beta);
DI=(beta-avebeta)/stdbeta;

return
