function waves2d(lambda)
%  1D waves emitted by an oscillating obstacle and/or moving
%            obstacle in a  dispersive media
%  O. Thual, version 02/11/14
%
% 
clear all
close all
% EXTERNAL LOOP
for para=0:.001:1
	% TUNING PARAMETERS
	%
	% alpha, beta : parameter of the dispersion relation 
	alpha=1; beta=para;
	% p : dimension n=2^p
	p=7; 
	% sigma2 : variance of the gaussian obstacle
	sigma2=.005; 
	%  eps : dissipation 
	eps=.0001;
	% omega0 : forcing frequency
	omega0 = 2*pi;  
	% visu : =1 for visualizations, =0 no visualization 
	visu=1;
	%
	% PARAMETERS FOR THE DOMAIN DEFINITION
	% nomber of points in one direction
	n=2^p;
	% dimensions of the periodic box
	l1=2*pi;
	% initialisation of somme arrays
	nh=n/2;
	m=[-nh:nh-1]; 
	dx=l1/n ; 
	X=m*dx;
	k1b=2*pi/l1 ;
	K1=m*k1b ; 
	% SHAPE OF THE FORCING
	% position of the obstacle
	x0=0 ;; 
	% forcing shape as a gaussian in Physical Space
	la=.5/sigma2;
	F=exp(-la*(X-x0).^2);
	% /fig1/ plot of the forcing shape in Pysical Space
	iflag=visu
	iflag=0
	if iflag==1
		figure 
		plot(F); 
		hold on
		title(' Forcing shape Phys. Space','FontSize',20)
		pause
	end
	% transfer to the Spectral Space
	FK=fftshift(fft(F));
	% /fig2/ plot of the forcing in Spectral Space
	iflag=visu
	iflag=0
	if iflag==1
		figure
		plot(real(FK));
		hold on
		title(' Real part forcing in Sp. Sapace','FontSize',20)
		pause
	end
	% DISPERSION RELATION
	% dispersion function B
	B=-i*eps-omega0+alpha*K1-beta*K1.^3; 
	% Computation 
	E= FK./B;
	% /fig3/ plot of the integrand
	iflag=visu
	iflag=0
	if iflag==1 
		figure
		plot(real(E));
		hold on
		title(' Integrand in Sp. Space','FontSize',20)
		pause
	end
	% COMPUTING OT THE WAVES
	% Back from Spectral Space to Physical Space
	U=ifft(fftshift(E));
	%  /fig4/ plot of the wave shape at t=0 
	iflag=visu;
	iflag=0
	if iflag==1 
		figure
		plot(real(U))
		hold on
		title(' wave shape at t=0 Phys. Space ','FontSize',20)
		pause
	end
	% CLIPPING OF THE WAVES
	% pclip : percentage of clipping points 
	flaclip=0
	Uclip=U;
	if flaclip==1
		pclip=.05;
		lclip=round(n*pclip);
		ix0=round(n/2+1+x0/l1*n)
		ilow=ix0-lclip 
		ihigh=ix0+lclip-1
		% clipping
		for ix=ilow:ihigh
			Uclip(ix)=0;
		end
	end
	% /fig5/ plot of the clipped wave at t=0
	iflag=0;
	if iflag==1 
		figure
		plot(real(Uclip))
		hold on
		title(' Train onde clippe espace Physique ','FontSize',20)
		pause
	end
	% /fig6/ time animation of the clipped wave
	flag=1;
	if flag==1
		figure
		for t=0:.1:10
			scal=.1;
			UR=real(scal*t+Uclip*exp(-i*omega0*t));
			plot(UR)
			hold on
			%		axis([0 140 -.2 .2])
			%		drawnow
			%		hold off
		end
		drawnow
	end
	%
end

