%  2D waves emitted by an oscillating obstacle and/or moving
%            obstacle in a  dispersive media
%  O. Thual, version 02/11/16
%  Rossby 05/12/02 
function obstacle2d
clear all
close all
global ratio flaclip tmax
global disprelation
% Screen demo
disp('Hello : oscillating obstacle in dispersive media')
while 1
	close all
	choix=input('1: sound 2: Standrew, 3: ship, 4: rossby 9:return   ?')
	choixdef=4;
	% 1: sound 2: Standrew, 3: ship
	if isempty(choix)
		choix=choixdef
	end
	% choix=2 : St andrew cross
	if choix==2 
		disprelation='dispstandrew'
		wstandrew
	end
	% choix=3 : ship
	if choix==3
		disprelation='dispship'
		wship
	end
	% choix=4 : rossby
	if choix==4
		disprelation='disprossby'
		wrossby
	end
	% choix=1 : sound
	if choix==1
		disprelation='dispsound'
		wsound
	end
	% choix=9 : return
	if choix==9 
		disp('Good Bye !')
		return
	end
end
%
% FUNCTIONS
%
%
% << inidomain >>
%
function FK=inidomain
global n omega0 x0 y0 l1 l2 
global px0 py0
global K1 K2
global x y X Y
% PARAMETERS TO VARY
% visu : =1 for visualizations, =0 no visualization 
visu=0;
% p : dimension n=2^p
p=7; 
% sigma2 : variance of the gaussian obstacle
sigma2=.01; 
% PARAMETERS FOR TIME AND SPACE
% omega0 : forcing frequency
omega0 =2*pi;  
% nomber of points in one direction
n=2^p;
% dimensions of the periodic box
l1=10; l2=10;
% position of the obstacle
x0=px0*l1 ; y0=py0*l2;
% initialisation of somme arrays
nh=n/2;
m=[-nh:nh-1]; 
dx=l1/n ; dy=l2/n;
x=m*dx; y=m*dy; 
[X,Y]= meshgrid(x,y);
k1b=2*pi/l1 ; k2b=2*pi/l2 ; 
k1=m*k1b ; k2=m*k2b;
[K1,K2]=meshgrid(k1,k2); 
% SHAPE OF THE FORCING
% forcing shape as a gaussian in Physical Space
la=.5/sigma2;
F=exp(-la*((X-x0).^2+(Y-y0).^2));
% /fig1/ plot of the forcing shape in Pysical Space
iflag=visu;
if iflag==1
	figure 
	surf(F); 
	hold on
	title(' Forcing shape Phys. Space','FontSize',20)
	pause
end
% transfer to the Spectral Space
FK=fftshift(fft2(F));
% /fig2/ plot of the forcing in Spectral Space
iflag=visu;
if iflag==1
	figure
	surf(abs(FK));
	hold on
	title(' modulus of forcing in Sp. Space','FontSize',20)
	pause
end
%
% << shapewf >>
%
function Uclip=shapewf(FK)
global disprelation
global flaclip
% visu : =1 for visualizations, =0 no visualization 
visu=0;
% dispersion relation 
B=eval(disprelation);
% Computation 
E= FK./B;
% /fig3/ plot of the integrand
iflag=visu;
if iflag==1 
	figure
	surf(abs(E));
	hold on
	title(' modulus of Integrand in Sp. Space','FontSize',20)
	pause
end
% COMPUTING OT THE WAVES
% Back from Spectral Space to Physical Space
U=ifft2(ifftshift(E));
%  /fig4/ plot of the wave shape at t=0 
iflag=visu;
if iflag==1 
	figure
	surf(real(U))
	hold on
	title(' real part wave shape at t=0 Phys. Space ','FontSize',20)
	pause
end
% CLIPPING OF THE WAVES
if flaclip==1
	Uclip=clipping(U);
else
	Uclip=U;
end
%
% << clipping >> 
%
function Uclip=clipping(U)
global n x0 y0 l1 l2
% plcip : pourcentage of clipped domain
pclip=0.1;
% clipping indices 
lclip=round(n*pclip);
ix0=round(n/2+1+x0/l1*n);
jy0=round(n/2+1+y0/l2*n); 
ilow=ix0-lclip ;
ihigh=ix0+lclip-1;
jlow=jy0-lclip ;
jhigh=jy0+lclip-1;
% clipping
Uclip=U;
for ix=ilow:ihigh
	for jy=jlow:jhigh
		Uclip(jy,ix)=0;
	end
end
% /fig5/ plot of the clipped wave at t=0
iflag=0;
if iflag==1 
	figure
	%	contour(real(Uclip))
	surf(real(Uclip))
	hold on
	title(' Train onde clippe espace Physique ','FontSize',20)
	pause
end
% 
% << animwaves2 >>
%
function animwaves2(Uclip)
global n nmax M
global dt tmax leg omega0
global dc cmin cmax
global x y X Y
global drawl
% jpg =1 for jpeg file creation
jpg=0;
for t=0:dt:tmax
	clear UR
	UR=real(Uclip*exp(-i*omega0*t));
	contour(x,y,UR,cmin:dc:cmax)
	% line for standrew
	if drawl==2
		hold on
		drawcrossline
		hold off
	end
	% line for ship
	if drawl==3
		hold on
		drawbetaline
		hold off
	end
	% line for rossby
	if drawl==4
		hold on
		drawsingline
		hold off
	end
	title(leg,'FontSize',20)
	URmin=min(min(UR));
	URmax=max(max(UR));
	tbounds=[t URmin URmax cmin cmax dc]
	if n==0
		pause
		M=moviein(n);
	else
		M(:,n)=getframe;
		if jpg==1 
			filew(gcf,n)
		end
	end
	n=n+1
end
if n==nmax
	newplot
	movie(M,2,3)
end
% 
% << animwaves >>
%
function animwaves(Uclip)
global dt tmax leg omega0
global dc cmin cmax
global x y X Y
% jpg =1 for jpeg file creation
jpg=0;
figure
%hold off
nmax=tmax/dt; 
n=0
for t=0:dt:tmax
	clear UR
	UR=real(Uclip*exp(-i*omega0*t));
	contour(x,y,UR,cmin:dc:cmax)
	title(leg,'FontSize',20)
	URmin=min(min(UR));
	URmax=max(max(UR));
	tbounds=[t URmin URmax cmin cmax dc]
	if n==0
		pause
		M=moviein(n);
		n=n+1;
	end
	M(:,n)=getframe;
	if jpg==1
		filew(gcf,n)
	end
	n=n+1;
end
movie(M,2,3)
% 
% << animsurf >>
%
function animsurf(Uclip)
global dt tmax leg omega0
global dc cmin cmax
global x y X Y
global l1 l2 X Y
figure
newplot
%	axis([-l1/2 l1/2 -l2/2 l2/2 -.05 .05])
%	shading interp;                       
%	colormap(pink); 	
for t=0:dt:tmax
	UR=real(Uclip*exp(-i*omega0*t));
	surfl(X,Y,UR)
	%	set(h2,'EraseMode','normal')
	%	contour(x,y,UR,cmin:dc:cmax)
	title(leg,'FontSize',20)
	URmin=min(min(UR));
	URmax=max(max(UR));
	tbounds=[t URmin URmax cmin cmax dc]
	pause(0)
end
%
% FUNCTION DEPENDING OF THE DISPERSION
%
% SAINT ANDREW CROSS
%
% << wstandrew >> 
%
function wstandrew
% StAndrew : ratio=.4:.1:1, eps=.6-.5*ratio;
% ratio : omega0 = ratio * N 
global eps ratio omega0 flaclip tmax
global dt dc cmin cmax
global px0 py0
global leg
global n nmax M
global drawl
drawl=2
% position of the obstacle
px0=0 ; py0=0;
% maximum value for contour
cmax=.1;
cmin=-cmax;
dc=cmax/10;
% timestep for animation
dt=.1; 
% Initialisation of the forcing
FK=inidomain;
% infinite loop of the demo
% demos for wstandrew
while 1 
	choice=input(' 1: demo1 (default), 2: demo2 ? 9 to quit  ')
	% default choice
	if isempty(choice)==1
		choice=1
	end
	% 9 : quit
	if choice == 9
		break
	end
	% choice 1 : ratio
	if choice==1
		while 1
			ratio=input(' ratio ?  (defaut=.7)  9 to quit ?   ')
			if isempty(ratio)==1
				ratio=.7
			end
			% 9 : quit
			if ratio == 9
				break
			end
			% animation with no clip
			flaclip=0 ; 
			tmax=1/ratio;
			Uclip=shapewf(FK);
			leg=sprintf('omega0/N = %0.2g     dissipation=%0.2g',ratio,eps);
			animwaves(Uclip)
			
		end
	end
	% choice 2 : loop on ratio
	if choice ==2
		tmax=1;
		nlop=3;	rmin=.4;dr=.2 ; rmax=rmin+(nlop-1)*dr;
		nmax=nlop*(round(tmax/dt)+1)
		n=0;
		for ratio=rmin:dr:rmax
			% animation with clip
			flaclip=1 ; 
			Uclip=shapewf(FK);
			leg=sprintf('omega0/N = %0.2g     dissipation=%0.2g',ratio,eps);
			animwaves2(Uclip)
			%	whos global
		end
	end
end
% 
% << drawcrossline >>
%
function drawcrossline
global x y x0 y0
global ratio
beta=asin(ratio);
tbeta=tan(beta)
ybeta=y0+tbeta*(x-x0);
plot(x,ybeta,'.k')
ybeta=y0-tbeta*(x-x0);
ybeta=y0-tbeta*(x-x0);
plot(x,ybeta,'.k')
%
% << dispstandrew >>
%
function  B=dispstandrew
global eps  omega0  ratio 
global K1 K2
%  eps : dissipation 
eps=.3-.5*(ratio-.7); 
eps=.2-.5*(ratio-.7); 
% wave number
N=omega0/ratio;
THETA=atan2(K2,K1);
% dispersion function B
iflab=2;
if iflab==1
	omegae=omega0-i*eps
	B=-K2*(omegae^2*K2.^2-(N^2-omegae^2)*K1.^2)./omegae;
end
if iflab==2
	B=-i*eps-omega0+N*abs(cos(atan2(K2,K1))); 
end
%
% SHIP WAVES
%
%
% << wship >> 
%
function wship
global V
global dt tmax eps omega0 flaclip
global dc cmin cmax
global px0 py0
global leg
global n nmax M
global drawl
drawl=3;
% position of x0=px0*l1
px0=1/3 ; py0=0 ;
% maximum value for contour camx=1
cmax=1;
cmin=-cmax;
dc=(cmax-cmin)/21;
% time step
dt=.01;
% Initialisation of the forcing
FK=inidomain;
% infinite loop of the demo
% demos for wsship
while 1 
	choice=input(' 1: demo1 (default), 2: demo2 ? 9 to quit  ')
	% default choice
	if isempty(choice)==1
		choice=1
	end
	% 9 : quit
	if choice == 9
		break
	end
	% choice 1 : V
	if choice==1
		while 1
			V=input(' Velocity ?  default=.5,  9 to quit  ')
			if isempty(V)==1
				V=.5
			end
			% 9 : quit
			if V == 9
				break
			end
			% animation with no clip
			flaclip=0 ; tmax=0;
			Uclip=shapewf(FK);
			leg=sprintf('V= %0.2g     dissipation=%0.2g',V,eps);
			% animation and drawing the lines
			animwaves(Uclip) 
			hold on
			drawbetaline
			% 		
			% loop on ship velocity
			%for V=.28:.1:.28
			%for V=.26:.01:.30
			
		end
	end
	% choice 2 : multiple V
	if choice==2
		nlop=9 ;
		Vmin=.25; dV=.05 ; Vmax=Vmin+(nlop-1)*dV;
		nmax=nlop
		n=0;
		hold off 
		for V=Vmin:dV:Vmax
			%for V=.2:.1:2
			flaclip=0 ; tmax=0;
			Uclip=shapewf(FK);
			leg=sprintf('V= %0.2g     dissipation=%0.2g',V,eps);
			% animation and drawing the lines
			animwaves2(Uclip) ;
			% 
			%	whos global
		end
	end
end
% 
% << drawbetaline >>
%
function drawbetaline
global x y x0 y0
beta=asin(1/3);
tbeta=tan(beta)
ybeta=y0+tbeta*(x-x0);
plot(x,ybeta)
ybeta=y0-tbeta*(x-x0);
ybeta=y0-tbeta*(x-x0);
plot(x,ybeta)
%
% << dispship >>
%
function  B=dispship
global eps  V
global K1 K2
%  eps : dissipation .01
eps=.01;
eps=.01+.15*(V-.25)
% dispersion function B
K=sqrt(K1.^2+K2.^2); 
B=-i*eps-sqrt(K)-K1*V; 
%
% ROSSBY
%
%
% << wrossby >> 
%
function wrossby
global V
global a angl beta C
global dt tmax eps omega0 flaclip
global dc cmin cmax
global px0 py0
global leg
global n nmax M
global drawl
drawl=4;
% position of x0=px0*l1
px0=1/3 ; py0=0 ;
px0=0 ; py0=0 ;
% maximum value for contour camx=1
cmax=.004;
cmin=-cmax;
dc=(cmax-cmin)/11;
% time step
dt=.01;
% Initialisation of the forcing
FK=inidomain;
% Dispersion relation parameters
a=0;
angl=a*2*pi/360;
C=10; 
beta=1;
% infinite loop of the demo
% demos 
while 1 
	choice=input(' 1: demo1 (default), 2: demo2 ? 9 to quit  ')
	% default choice
	if isempty(choice)==1
		choice=1
	end
	% 9 : quit
	if choice == 9
		break
	end
	% choice 1 : V
	if choice==1
		while 1
			V=input(' Velocity ?  default=.01,  9 to quit  ')
			if isempty(V)==1
				V=.01
			end
			% 9 : quit
			if V == 9
				break
			end
			% animation with no clip
			flaclip=0 ; tmax=0;
			Uclip=shapewf(FK);
			leg=sprintf('V= %0.2g a=%0.3g  C=%0.2g bet=%0.2g eps=%0.2g  ',V,a,C,beta,eps);
			% animation and drawing the lines
			animwaves(Uclip)
			pause
			hold on 
			drawsingline
			% 					
		end
	end
	% choice 2 : multiple V
	if choice==2
		nlop=37 ;
		% Dispersion relation parameters
		a=0;
		angl=a*2*pi/360;
		C=10; 
		beta=1;
		V=.01
		% New choices
		a=30;
		%
		Vmin=.001; dV=.001 ; Vmax=Vmin+(nlop-1)*dV;
		Cmin=10; dC=10 ; Cmax=Cmin+(nlop-1)*dC;
		amin=0; da=10; amax=amin+(nlop-1)*da;
		nmax=nlop
		n=0;
		hold off 
		%for V=Vmin:dV:Vmax
		% for C=Cmin:dC:Cmax
		for a=amin:da:amax
			angl=a*2*pi/360;
			flaclip=0 ; tmax=0;
			Uclip=shapewf(FK);
			leg=sprintf('V= %0.2g a=%0.3g  C=%0.2g bet=%0.2g eps=%0.2g  ',V,a,C,beta,eps);
			% animation and drawing the lines
			animwaves2(Uclip) ;
			% 
			%	whos global
		end
	end
end
% 
% << drawsingline >>
%
function drawsingline
global x y x0 y0
global angl
tangl=tan(angl)
ybeta=y0+tangl*(x-x0);
plot(x,ybeta)
%
% << disprossby >>
%
function  B=disprossby
global eps  V omega0
global a angl beta C
global K1 K2
%  eps : dissipation .01
V1=V*cos(angl);
V2=V*sin(angl);
eps=.01;
eps=.01;
%C=.000001; 
% dispersion function B
KS=K1.^2+K2.^2;
omega1=omega0*.000001;
Omega=-beta*K1./(KS+C)+K1*V1+K2*V2;
B=-i*eps-omega1-Omega;
%contour(Omega,.5*omega1:omega1/20:1.5*omega1)
ndc=5;
cmin=0+omega1; dc=.01 ; cmax=cmin+ndc*dc;
figure(1)
contour(K1,K2,Omega,[cmin:dc:cmax])
hold on
contour(K1,K2,Omega,[-cmax:dc:cmin],'-.k')
grid
leg=sprintf('V= %0.2g a=%0.3g  C=%0.2g bet=%0.2g eps=%0.2g  ',V,a,C,beta,eps);
title(leg,'FontSize',20)
hold off
%surf(Omega)
figure(2)
%
% SOUND
% 
% << wsound >> 
%
function wsound
global csound
global dt tmax eps omega0 flaclip
global dc cmin cmax
global px0 py0
global leg
global n nmax M
global drawl
drawl=1
% position of x0=px0*l1
px0=0 ; py0=0 ;
% maximum value for contour
cmax=.1;
cmin=-cmax;
dc=(cmax-cmin)/20;
% time step
dt=.1;
% Initialisation of the forcing
FK=inidomain;
% demos for wsound
while 1 
	choice=input(' 1: demo1 (defaut), 2: demo2 ? 9 to quit  ?')
	% default choice
	if isempty(choice)==1
		choice=1
	end
	% 9 : quit
	if choice == 9
		break
	end
	% choice 1 : c
	if choice==1
		while 1
			csound=input(' Velocity csound  ? default=1,  9 to quit  ?')
			if isempty(csound)==1
				csound=1
			end
			% 9 : quit
			if csound == 9
				break
			end
			% animation with no clip
			tmax=1;
			flaclip=0 ; 
			Uclip=shapewf(FK);
			leg=sprintf('c= %0.2g     dissipation=%0.2g',csound,eps);
			% animation and drawing the lines
			animwaves(Uclip) ;
			%	animsurf(Uclip) ;
			hold on 
		end
	end
	if choice==2
		% loop on sound speed
		tmax=1.;
		nlop=3; csmin=1; dcs=.5; csmax=csmin+(nlop-1)*dcs;
		nmax=nlop*(round(tmax/dt)+1);
		n=0;
		hold off
		for csound=csmin:dcs:csmax
			flaclip=0 ; 
			Uclip=shapewf(FK);
			leg=sprintf('c= %0.2g     dissipation=%0.2g',csound,eps);
			% animation and drawing the lines
			animwaves2(Uclip) ;
			hold off
			pause
			%	animsurf(Uclip) ;
			% 
			%	whos global
		end
	end
end
%
% << dispsound >>
%
function  B=dispsound
global eps  csound omega0
global K1 K2
%  eps : dissipation .01
eps=.5;
% dispersion function B
K=sqrt(K1.^2+K2.^2); 
B=-i*eps-omega0+csound*K; 
% 
% << 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')




