% Zoe Boekelheide
% Math 164: Scientific Computing
% Harvey Mudd College
% Started April 8, 2003
% Finished May 2, 2003

% Code to make surface gravity waves break. Units are such that
% the wavelength of the wave is 2*pi, and gravity g = 1. The density of 
% water is also taken to be 1.

% Things that can be varied to change the wave physically:
%   y(x)  --  this changes the shape of the surface
%   Phi(x, y)  --  this changes the initial velocity of particles in the wave (stability)
%   p_s  -- this changes the surface pressure on the wave (actually p_s is contained
%   in timestep3.m if you want to change it)


% function makewave(a) takes as an input 2*pi*a
% suggested input values are 0.2, 0.4, 0.8, 1.6, 2.4
% but the program will try to accommodate for other values

function makewave(a);

% Define initial wave function to look pretty!
num_points=60;
num_steps = 40;
delta_x = 2*pi/num_points;
x=0:delta_x:2*pi-delta_x;
y=a*(2/1.6)*((sech(x-0.95*pi)).^2 - 0.5 - 0.3*cos(x+0.15));

% Define Phi to solve Laplace's equation. lambda can be whatever. In order for
% Phi to be periodic, lambda should be 1. But to add a instability, we can
% make Phi slightly unperiodic and make lambda something close to 1 that is
% not exactly 1. It turns out Phi will still end up matching up very closely
% even when lambda is not exactly 1.
lambda=0.75;
A=1;
B=1;
Phi=exp(lambda*y).*(A*sin(lambda*x)+B*cos(lambda*x));


%Finally do some time stepping


% Master arrays hold the values for x, y, and Phi at all time steps

Master_X = zeros(num_points, num_steps);
Master_Y = zeros(num_points, num_steps);
Master_Phi = zeros(num_points, num_steps);

Master_X(:, 1) = x';
Master_Y(:, 1) = y';
Master_Phi(:, 1) = Phi';

% this sequence specifies what time step to use; dependent
% on the amplitude (assuming the original Phi and y(x))
if((2.39 < a) & (a < 2.41))
    h=.0000247
elseif((1.59 < a) & (a < 1.61))
    h=.0001357
elseif((0.79 < a) & (a < 0.81))
    h=.0005589
elseif((0.39 < a) & (a < 0.41))
    h=.000852
elseif ((0.19 < a) & (a < 0.21))
    h=.000865
else
h=0.000907*exp(-a^2/1.4085^2)+ 1.4984e-05
end


surfability = 0;

for t_step = 1:num_steps-1

%TIMESTEP!
new_crappy_vector = timestep3(num_points, x, y, Phi, h);

for x_step = 1:num_points
Master_X(x_step, t_step+1) = new_crappy_vector(x_step);
Master_Y(x_step, t_step+1) = new_crappy_vector(x_step + num_points);
Master_Phi(x_step, t_step+1) = new_crappy_vector(x_step + 2*num_points);
end

% This is a smoothing function that runs every timestep.
for x_step = 1:num_points
        if (x_step > 2 & x_step < num_steps -2) 
		    x(x_step) = (Master_X(mod(x_step-1-1, num_points)+1, t_step+1) + Master_X(x_step, t_step+1) + Master_X(mod(x_step+1-1, num_points)+1, t_step+1))/3; 
        else
            x(x_step) = Master_X(x_step, t_step+1);
        end
        y(x_step) = (Master_Y(mod(x_step-1-1, num_points)+1, t_step+1) + Master_Y(x_step, t_step+1) + Master_Y(mod(x_step+1-1, num_points)+1, t_step+1))/3;
		Phi(x_step) = (Master_Phi(mod(x_step-1-1, num_points)+1, t_step+1) + Master_Phi(x_step, t_step+1) + Master_Phi(mod(x_step+1-1, num_points)+1, t_step+1))/3;             
end
% Calculate surfability for one timestep
for i = 2:num_points-1
if ( (atan((y(i+1)-y(i))/(x(i+1)-x(i))) < .611) & (atan((y(i+1)-y(i))/(x(i+1)-x(i))) > 0) & ((y(i+1)-y(i))/(x(i+1)-x(i)) > (y(i)-y(i-1))/(x(i)-x(i-1))))
surfability = surfability + sqrt((y(i+1)-y(i)).^2 + (x(i+1)-x(i)).^2)*h;
end
end

end


for i=1:num_steps/4
subplot(num_steps/4,1,i); plot(Master_X(:,4*i), Master_Y(:,4*i));
end

surfability

% function will output the value it chose for h, and the surfability of the wavefunction