function [] = surfcomp(Ho,slop,ho,T,adeg)
%
%  [] = surfcomp(Ho,slop,ho,T,adeg)
%
%
% An non-interactive function adapted from "SURFER.EXEC" on the mainframe.
% For oc4213, Nearshore Processes, Dr. Thornton
% Laboratory #3
%
% Called by initializing program surfer.m.
%
% Function inputs several wave statistics, then
% calculates cross-shore values of significant wave height,
% longshore current velocity, and percentage of breaking waves.
% Calculations are based on the Thornton-Guza surf model and the
% Longuet-Higgins longshore current model.
% Results are plotted on a single figure, along with some of the
% initial wave statistics.  Additional results, including maximum
% values of significant wave height and longshore current, breaker
% angle, and modified surf index, are displayed on a second figure.
%
% Written by Mark Orzech, 6/98

global J;
vmat=version;

if isempty(Ho)
 Ho = 1;
end
if isempty(slop)
 slop = 0.01;
end
if isempty(ho)
 ho = 50;
end
if isempty(T)
 T = 10;
end
if isempty(adeg)
 adeg = 0;
end

fp = 1/T;    
g = 9.81;		% gravitational constant
Co = g*T/(2*pi);	% deep water wave speed
Lo = g*(T^2)/(2*pi);	% deep water wavelength


% Check for reasonable values:

if (Ho<0.1) | (Ho>5)
 J=figure; set(J,'position',[600 20 560 420]);
 axes; set(gca,'visible','off');
 text(-.1,.5,'Wave height is not within required range.  Please try again.');
elseif (ho<4*Ho)
 J=figure; set(J,'position',[600 20 560 420]);
 axes; set(gca,'visible','off');
 text(-.1,.5,'Water depth is not deep enough for chosen wave height.  Please try again.');
 ho = 50;
elseif (Ho/Lo > 1/7)
 J=figure; set(J,'position',[600 20 560 420]);
 axes; set(gca,'visible','off');
 text(-.1,.5,'Waves too steep.  Please choose a longer period.');
elseif (slop < 0)
 J=figure; set(J,'position',[600 20 560 420]);
 axes; set(gca,'visible','off');
 text(-.1,.5,'Slope must be positive.  Please try a new slope.');
elseif (T<7.5 | T>17.49)
 J=figure; set(J,'position',[600 20 560 420]);
 axes; set(gca,'visible','off');
 text(-.1,.5,'Peak wave period must be between 7.5 and 17.49 sec.  Please try again.');
elseif (abs(adeg)>90)
 J=figure; set(J,'position',[600 20 560 420]);
 axes; set(gca,'visible','off');
 text(-.1,.5,'Initial wave direction must be between -90 and 90.  Please try again.');

 
else	% If above requirements satisfied, continue with calculations

 if (~isempty(J))		% Close error window if one was used
   close(J);
   J=[];
 end
 
a0 = adeg*pi/180;	% convert alpha to radians for calculations
Cf = 0.005;		% friction coefficient
B = 1;			% breaking fraction of wave face
gam = 0.4;  	
dh = ho/100;
h = 0.001:dh:(ho+0.001);
x = h/slop;


% Calculate Hrms, then Hs at each cross-shore location

ro = (Ho^2)*(ho^.5);
a = 23*((g/pi)^.5)*(gam^4)*slop/(15*fp*(B^3));

hrms1 = (a^.2)*(h.^.9);
hrms2 = (ho^(-5.75))-(a/(ro^2.5));

Hrms = hrms1.*((1-(h.^5.75).*hrms2).^(-.2));

Hs = Hrms*sqrt(2);


% Calculate Longshore current velocity

V = 3*pi*(B^3)*fp*(g^.5)*sin(a0)*(Hrms.^6)./(8*Cf*(gam^4)*Co*(h.^4.5));


% Calculate percentage of breaking waves at each location

num = (Hrms./(gam*h)).^2;
Pb = (num.^2).*(1 - ((1+num).^(-1)));
Pb = Pb.*(Pb<1)+(Pb>=1);		% So that maximum of Pb is 1.


% Determine main type of breaker

ep = slop/(sqrt(2*pi*max(Hs)/(g*T^2)));
if ep < 0.4
 typ = 'Spilling';
elseif ep < 2.0
 typ = 'Plunging';
else
 typ = 'Surging';
end


% Plot results

H=figure(2); clf; set(H,'position',[20 20 560 420]); set(H,'numbertitle','off');
set(H,'name','KEY WAVE STATISTICS');
subplot(4,1,1); 
plot(x,100*Pb); axis([0 max(x) 0 100]); ylabel('Percent');
title(['alpha_0 = ' num2str(adeg) '^o, H_{so} = ' num2str(Ho) ' m, slope = ' num2str(slop) ', h_o = ' num2str(ho) ' m, T_p = ' num2str(T) ' sec'],'fontsize',10);
text(0.72*max(x),90,'Percent Breaking','fontsize',10);
if vmat(1)=='5'
  set(gca,'XTickLabel',[]);
else
  set(gca,'xticklabel',[]);
end

subplot(4,1,2);
plot(x,Hs); 
if max(Hs)<3
 ymax=3.6;
else
 ymax=1.2*max(Hs);
end
axis([0 max(x) 0 ymax]); ylabel('meters');
if vmat(1)=='5'
  set(gca,'XTickLabel',[]);
else
  set(gca,'xticklabel',[]);
end

text(0.93*max(x),0.9*ymax,'Hs','fontsize',10);

subplot(4,1,3);
plot(x,V); 
if max(abs(V))==0
 ymax=1;
elseif max(abs(V))<3
 ymax=3;
else
 ymax=1.05*max(abs(V));
end
axis([0 max(x) -ymax ymax]);
hold on; plot([0:(max(x)/100):max(x)],[zeros(1,length(0:(max(x)/100):max(x)))],':'); 
ylabel('m/s');
text(0.7*max(x),0.8*ymax,'Longshore Current','fontsize',10);
text(0.05*max(x),0.8*ymax,'To Right','fontsize',10);
text(0.05*max(x),-0.8*ymax,'To Left','fontsize',10);
if vmat(1)=='5'
  set(gca,'XTickLabel',[]);
else
  set(gca,'xticklabel',[]);
end

subplot(4,1,4);
plot(x,h); axis([0 max(x) -0.1*max(h) 1.1*max(h)]);
text(0.75*max(x),0.02*max(h),'Water Depth, h','fontsize',10); 
ylabel('meters'); xlabel('Distance Offshore (m)');
set(gca,'ydir','reverse');
orient tall

H=figure(3); clf; set(H,'position',[600 20 560 420]); set(H,'numbertitle','off');
set(H,'name','RESULTS OF ANALYSIS');
axes; set(gca,'visible','off');

text(-.1,.9,'INPUT PARAMETERS:');
text(-.1,.8,['alpha_o=',num2str(adeg),'^o;  H_{so}=', num2str(Ho),' m;  slope=',num2str(slop),';  h_o=',num2str(ho),' m;  T_p=',num2str(T),' sec'],'fontsize',12);

if sum(Pb<1)
 text(-.1,.7,['All waves have broken shoreward of x = ' num2str(x(min(find(Pb<1)))) ' m'],'fontsize',16);
else
 text(-.1,.7,['All waves have broken seaward of x = ' num2str(max(x)) ' m, h = ' num2str(max(h)-0.001) ' m'],'fontsize',16);
end
 text(-.1,.6,['Maximum breaker height = ' num2str(max(Hs)) ' m at x = ' num2str(x(find(Hs==max(Hs)))) ' m'],'fontsize',16);

if max(V)~=0
 text(-.1,.5,['Maximum longshore current = ' num2str(V(find(abs(V)==max(abs(V))))) ' m/s at x = ' num2str(x(find(abs(V)==max(abs(V))))) ' m'],'fontsize',16);
else
 text(-.1,.5,'Maximum longshore current = 0 m/s','fontsize',16);
end

% Calculate modified surf index (call modsi subroutine)

hh = h(find(Hrms==max(Hrms)));

[angl,spill,plunge,surge,hbig,mdsi] = modsi(max(Hrms),T,hh,slop,max(abs(V)),a0);

if imag(angl)~=0
 text(-.1,.4,'Breaker angle not calculated; initial depth too shallow','fontsize',16);
else
 text(-.1,.4,['Breaker angle is ' num2str(angl) ' degrees'],'fontsize',16);
end

% Display results for breaker type

text(-.1,.3,['Breaker type:  ' num2str(spill) '% Spilling,   ' num2str(plunge) '% Plunging,   ' num2str(surge) '% Surging'],'fontsize',16);

if round(T)<8
 text(-.1,.2,'Modified Surf Index:  Wave period too small to calculate','fontsize',16);
elseif round(T)<=17
% continue;
else
 text(-.1,.2,'Modified Surf Index:  Wave period too large to calculate','fontsize',16);
end

if hbig==0
 text(-.1,.1,['Modified surf index = ' num2str(mdsi)],'fontsize',16);
else
 text(-.1,.1,['Modified surf index is greater than ' num2str(mdsi)],'fontsize',16);
end

%text(-.1,.5,'******************************************************');

end 	%% if

return;
