function huyghens
clear all
close all
global alpha beta kappa
global g
global theta
%
disp(' Trace de rayons par la methode de Huyghens')
% Important parmeters 
alpha=.01; beta=.01;          % slope and amplitude of slode oscill
nx=100; ny=100;              % number of space point for bathy
dt=1;                       % time step of Huyghens method
nlines=23;                   % number of lines
ntot=11;					 % number of circle on a phase line
k0=1;                        % wave number offshore
% Memory of fine parameters
% alpha, beta, dt, nlines, ntot,  k0 0.05   0.05   0.5   25   10   3 
% alpha, beta, dt, nlines, ntot,  k0  0.05   0.05   0.5   25   10   3 
% alpha, beta, dt, nlines, ntot,  k0  0.05   0.05   0.5   25   10   3 
% alpha, beta, dt, nlines, ntot,  k0  0.01   0.01   1   23   9   1 
disp(sprintf('alpha, beta, dt, nlines, ntot,  k0'))
disp(sprintf('  %0.5g ',alpha, beta,dt,nlines,ntot,k0))
% other parameters
g=10; 
% Bathymetry
xmax=20; xmin=-1; 
[x,dx]=divideit(xmax,xmin,nx);
ymax=10; ymin=-ymax;
[y,dy]=divideit(ymax,ymin,ny);
kappa=1.5*pi/ymax;
[X,Y]=meshgrid(x,y);
[H,H_x,H_y]=bathy(X,Y);
% 3D Visualisation of bathy
iflag=0;
if iflag==1
	figure(3)
	surf(X,Y,-H)
	xlabel('x','Fontsize',30)
	ylabel('y','Fontsize',30)
	rotate3d
	figure(1)
end
% 2D visu of bathy
contour(X,Y,H)
colormap(gray)
axis([-1 19 -10 10])
xlabel('x','Fontsize',30)
ylabel('y','Fontsize',30)
grid
% theta for the circles
ntheta=10; [theta,dtheta]=divideit(pi/2,3*pi/2,ntheta);
% Initial phase lines
ncir=ntot-1;
ycir=divideit(ymin,ymax,ncir);
x0=.9*xmax; xcir=0*ycir+x0;
kcirx=0*xcir-k0; kciry=0*xcir; kcir=kcirx*0+k0;
% Test
iflag=0;
if iflag==1
	test
end
% Iteration
% < movie
nmax=nlines;
n=0;
%   movie > 
hold on
for i=1:nlines
	[xcir,ycir,kcirx,kciry]=phaseline(xcir,ycir,kcirx,kciry,dt,ntot);
	% < movie
	if n==0 
		M=moviein(nmax);
	end
	n=n+1;
	M(:,n)=getframe; 
	filew(gcf,n);
	%   movie >
end
n
% < movie
if n==nmax
	movie(M,2,4)
end
%   movie >
return
%
%% Functions
%
% << bathy >>
%
function [h,h_x,h_y]=bathy(x,y)
global alpha beta kappa
kappay=kappa*y;
h=alpha*x+beta*cos(kappay);
h_x=0*h+alpha;
h_y=-beta*kappa*sin(kappay);
return
%
% << circles >>
%
function circles(x,y,cphi,dt,ntot);
for n=1:ntot
	xn=x(n); yn=y(n); cphin=cphi(n);
	cdt=cphin*dt;
	plot(xn,yn,'.')
	onecircle(xn,yn,cdt);
end
return
%
% << divideit >>
%
function [x,dx]=divideit(xmax,xmin,nx)
% dividing an interval into equal subintervals
dx=(xmax-xmin)/nx;
x=xmin:dx:xmax;
%
% << dispers >>
%
function [omega,omega_h,cphi]=dispers(k,h)
global g
gk=g*k; kh=k.*h; sh=sinh(kh); ch=cosh(kh); th=tanh(kh);
gkth=gk.*th;
omega=sqrt(gkth);
omega_h=.5*sqrt(gk).*k./(sqrt(sh.*ch).*ch);
cphi=omega./k;
return
%
% << onecircl >>
%
function onecircle(x,y,r);
global theta
xc=x+r*cos(theta);
yc=y+r*sin(theta);
plot(xc,yc)
return
%
% << phaseline >>
%
function [xcir,ycir,kcirx,kciry]=phaseline(xcir,ycir,kcirx,kciry,dt,ntot);
k2cir=kcirx.^2+kciry.^2;
kcir=sqrt(k2cir);
[hcir,hcir_x,hcir_y]=bathy(xcir,ycir);
[omegacir, omegacir_h, cphicir]=dispers(kcir,hcir);
plot(xcir,ycir,'r','Linewidth',2)
circles(xcir,ycir,cphicir,dt,ntot);
% Propagation
dtomegask2=dt*omegacir./k2cir;
dtomegacir_h=dt*omegacir_h;
xcir=xcir+kcirx.*dtomegask2;
ycir=ycir+kciry.*dtomegask2;
kcirx=kcirx-hcir_x.*dtomegacir_h;
kciry=kciry-hcir_y.*dtomegacir_h;
return
%
% << test >>
%
function test
figure(2)
ktest=.1:.01:10; htest=0*ktest+1;
[omega,omega_h,cphi]=dispers(ktest,htest);
plot(ktest,cphi)
hold on
plot(ktest,omega)
hold on
plot(ktest,omega_h)
% plot(ktest,sqrt(g*ktest.*tanh(ktest*htest)
figure(1)
% 
% << filew >>
%
% writing image on  files 
function filew(gcf,n)
[a,map]=capture(gcf);
filename=sprintf('image%2.2d.jpg',n)
imwrite(a,map,filename,'jpg')

