function [angl,spill,plunge,surge,hbig,mdsi] = modsi(hrms,T,hh,slope,Vmax,alpha0)
%
% [angl,spill,plunge,surge,hbig,mdsi] = modsi(hrms,T,hh,slope,Vmax,alpha0);
%
% Calculates modified surf index (see Joint Surf Manual), used for
% prediction of surf zone conditions and decisions on which types
% of landing craft may be used.  
%
% inputs:	hrms - maximum hrms at breaking (meters)
%		T - wave period in seconds
%		hh - depth at breaker location (meters)
%		slope - slope of bottom at breaker location
%		Vmax - maximum longshore current (m/s, >0)
%		alpha0 - initial offshore wave angle (radians)
%
% outputs:	angl - breaker angle (degrees)
%		spill,plunge,surge - percentage of breakers that are spilling,
%			plunging, surging, respectively
%		hbig - set to 1 if parameters lead to modified surf index
%			that is beyond table values
%		mdsi - modified surf index value calculated
%
% Also loads modification tables (stored in modtabs.mat): 
%  		bptab = breaker period modification table
%  		sptab = spilling breaker mod. table
%  		sutab = surging breaker mod. table
%  		watab = wave angle mod. table
%  		wtab = wind mod. table
%
%
% output:	mdsi - modified surf index
%
% Adapted from Commander Naval Surface Force NCSCSURF.F fortran program,
% subroutines "percent", "mdsrf1", and "mdsrf2".
%
% Mark Orzech, 5/98



%%%%%%%% I. Calculate percentages of breaking waves %%%%%%%%%%%%%%%%%%%%

% Calculate breaking wave angle using Snell's law

C0 = 9.81*T/(2*pi);
C = sqrt(9.81*hh);
angl = asin((C/C0)*sin(alpha0));
angl = 180*angl/pi;

angl = abs(5*round(angl/5));		% need wave angle to be >0 and multiple of 5

% Conversions and constants

hrms = hrms*3.28;		% convert m --> ft
hh = hh*3.28;	
Vm = Vmax/0.51308;		% convert from m/s --> knots

gamma = 0.4;
dgam = hh*gamma;

% Integration of breaking wave probability distribution by
% trapezoidal rule

% needed parameters

hbig = 0;
param = 980*(slope*T)^2;
dgam2 = dgam^2;
hrms2 = hrms^2;

% Use 50 integration intervals for each part of integration

numint = 50;
hlow = 0;
h = 0;
pbold = 0;
tot = 0;

% i=1 for surging breakers
% i=2 for surging and plunging breakers
% i=3 for all breakers

% Loop over three types of breakers

for i=1:3
 hhigh = 0.127*param;
 if i==2
  hhigh = 1.013*param;
 elseif i==3
  hhigh = 5.0*hrms;
 end
 
 % Integration wave height interval width
 delt = (hhigh-hlow)/numint;
 
 % Numerical integration loop
 
 for n=1:numint
  h = h+delt;
  h2 = h^2;
  temp = -h2/hrms2;
  ph = 0;		% ph is Rayleigh distribution
  fac12 = 1;		% fac12 is breaking wave weighting function
  
  if temp>=-40		% prevent exp overflow for ph and fac1
  
   ph = ((2*h)/hrms2)*exp(temp);
   fac1 = 1;
   temp = -h2/dgam2;
   if temp>=-40
    fac1 = 1-exp(temp);
   end
   fac2 = (hrms2*hrms2)/(dgam2*dgam2);
   fac12 = fac1*fac2;
   if fac12>1
    fac12 = 1;
   end
  end  % if temp 
  
  pb = fac12*ph;
  temp = (delt/2)*(pbold+pb);
  tot = tot+temp;
  pbold = pb;
  
 end % n=1:numint
 
 hlow = hhigh;
 frac(i) = tot;
 
end % i=1:3


% Percent of breakers in each category (1=surge, 2=plunge, 3=spill)

typ(1) = frac(1);
typ(2) = frac(2)-frac(1);
typ(3) = frac(3)-frac(2);

brk = frac(3);

tot = sum(typ);
tmp = 10*typ/tot;

surge = 10*round(10*typ(1)/tot);		% Round to nearest 10%
plunge = 10*round(10*typ(2)/tot);
spill = 10*round(10*typ(3)/tot);

% Make sure total is 100%

if surge+plunge+spill < 100
 tmp(find(tmp==min(tmp))) = tmp(find(tmp==min(tmp)))+0.5;
elseif surge+plunge+spill > 100
 tmp(find(tmp==max(tmp))) = tmp(find(tmp==max(tmp)))-0.5;
end

if surge+plunge+spill ~= 100
 surge = 10*round(tmp(1));		
 plunge = 10*round(tmp(2));
 spill = 10*round(tmp(3));
end

%%%%%%%%% II. Calculate modified surf index %%%%%%%%%%%%%%%%%%%%%%%%%%%

load modtabs;		% Loads modification table variables:

% Wave height modification term

hs = hrms*sqrt(2);	% significant wave ht (feet)

if round(2*hs)<=10
 hsi = 1 + (round(2*hs));	% hs index (1-11) for tables
else
 hsi = 11;
 hbig = 1;
end


% Period modification

if round(T)<8
% disp('Modified Surf Index:  Wave period too small to calculate (need Tmin = 7.5 sec)');
 mdsi = nan;
 return;
elseif round(T)<=17
 Ti = 18-round(T);
 Tmod = bptab(Ti,hsi);
else
% disp('Modified Surf Index:  Wave period too large to calculate (need Tmax = 17.49 sec)');
 mdsi = nan;
 return;
end
 

% Calculate breaker modification term

if spill > surge

 % Spilling breaker modification

 spi = 11 - spill/10;
 brmod = sptab(spi,hsi);

else

 % Surging breaker modification (no mods for plunging)

 sui = 11 - surge/10;
 brmod = sutab(sui,hsi);
end


% Wave angle modification

if angl<=40
 angli = 9 - (angl/5);
else
 angli = 1;
 hbig = 1;
end
wamod = watab(angli,hsi);


% Littoral current modification (just 3 times max current, knots)

lcmod = Vm * 3;	
if lcmod > 9
 hbig = 1;
 lcmod = 9;
end

% Sum modifications to calculate modified surf index


mdsi = hs + Tmod + brmod + wamod + lcmod;

return;
