function fft2D(N,tinit,tfin,dt,f,h,c)
% N est le nombre de subdivisions en espace dans chaque direction. 
% La largeur des subdivisions est dx
% tinit est le temps initial
% tfin est le temps final
% dt est le pas de temps
% f le type de fonction pour u0
% C'est une fonction paire centrŽe en (0.5,0.5) ˆ choisir parmi
%    'aa' tensorisŽ de la fonction creneau
%    'bb' tensorisŽ de la fonction chapeau
%    'ff' tensorisŽ de la fonction cos^2
%    'ra' fonction radiale creneau (marche plutot mal)
%    'rf' fonction radiale cos^2
% Les symŽtries sont utilisŽes pour simplifier les formules
% h parametre (demi-largeur du support)
% Ce programme est autonome et calcule la solution des ondes
% pour l'opŽrateur en espace Au = - u" + 4*c*pi^2 . 
% avec conditions aux limites periodiques, a l'aide de la
% transformation de fft2.m

% Preparation des proprietes de la figure 1
figure(1)
WHITEBG(1,'k')
set(1,'Renderer','zbuffer')

% La partition en espace
dx = 1/N ;
m = 1/2 ;
x1 = [0:dx:1-dx] ;
x2 = [0:dx:1-dx] ;
xmi = max(0,m-h) ;
xma = min(1,m+h) ;
if f == 'ra'
	for j = 1:N ;
		for k = 1:N
			y(j,k) = ( (x1(j)-m)^2 + (x2(k)-m)^2 < h^2 ) ;
		end
	end
elseif f == 'rf'
	for j = 1:N ;
		for k = 1:N
			r(j,k) = sqrt((x1(j)-m)^2 + (x2(k)-m)^2) ;
			y(j,k) = ( (x1(j)-m)^2 + (x2(k)-m)^2 < h^2 ) ;
			y(j,k) = y(j,k) * cos( pi*r(j,k)/(2*h) )^2 ;
		end
	end
elseif f == 'aa'
	y1 = (xmi<x1 & x1<xma) ;
	y2 = (xmi<x2 & x2<xma) ;
	y = y1' * y2 ;
elseif f == 'bb'
	y1 = (xmi<x1 & x1<xma) .* max (0*x1 , 1 - abs( (x1-m)/h ) ) ;
	y2 = (xmi<x2 & x2<xma) .* max (0*x2 , 1 - abs( (x2-m)/h ) ) ;
	y = y1' * y2 ;
elseif f == 'ff'
	y1 = (xmi<x1 & x1<xma) .* cos( pi*(x1-m)/(2*h) ).^2 ;
	y2 = (xmi<x2 & x2<xma) .* cos( pi*(x2-m)/(2*h) ).^2;
	y = y1' * y2 ;
end
z = fft2(y) ;
FC = z/(N*N) ;
for t = tinit:dt:tfin
	zzt = 0*y ;
	for j = 1:N/2 
		for k = 1:N/2
			jm1 = j-1 ;
			km1 = k-1 ;
			if jm1 == 0, c1 = 1; else c1 = 2 ; end
			if km1 == 0, c2 = 1; else c2 = 2 ; end
			zzt = zzt + FC(j,k) * cos(2*pi*jm1*x1)' * cos(2*pi*km1*x2)...
			 * c1 * c2 * cos(2*pi*sqrt(jm1*jm1+km1*km1+c)*t) ;
		end
	end
	zr = zzt ;
	zmax = max( max(abs(zr)) ) ;
	contourf(x1,x2,zr,[-zmax:zmax/5:-zmax/5 -zmax/15 zmax/15 zmax/5:zmax/5:zmax]) ;
	caxis([-zmax zmax])
	colorbar
	title (['t  =  ',num2str(t),'    &    zmax  =  ',num2str(zmax,'%1.3e')])
	pause(0.05)
end

% Apres des experiences avec c=0, on pourra essayer
% fft2D(50,0,2,0.2,'rf',0.1,-0.4)
% fft2D(50,0.45,1,0.05,'rf',0.1,-2)
