% Dominic Grenier décembre 2003
% a,b : dimensions de la section du guide d'ondes [m]
% d : longueur d'observation en terme de lambda_g
% epsr : constante diélectrique du matériau dans le guide d'ondes
% fo : fréquence d'opération [Hz]
% mode : TE:0  TM:1
% m,n : indices du mode
% Nx, Ny, Nz : nombre de points pour l'afichage des graphiques

a=0.02286;
b=0.01016;
d=2;
epsr=1;
fo=21e9;

mode=1;
m=2; n=1;

Nx=21; Ny=11; Nz=31;
%choix des plans de coupe 0/1 [z=0 x=x_o y=y_o]
plan=[1 1 1];

[gamma,vpz,etag]=guidepar(fo,mode,m,n,a,b,epsr);
d=d*2*pi/imag(gamma);

x=0:a/(Nx-1):a;
y=0:b/(Ny-1):b;
z=0:d/(Nz-1):d;
[Z,X,Y]=meshgrid(z,x,y);

[Ex,Ey,Ez,Hx,Hy,Hz]=guideEH(fo,mode,m,n,a,b,epsr,X,Y,Z);

%test pour s'assurer que le mode existe
if ~isempty(Ex)
clf reset;
    
%mesh(X(:,1,:),Y(:,1,:),reshape(sqrt(Ex(:,1,:).^2+Ey(:,1,:).^2),Nx,Ny));

%les 3 plans possibles
nplans=sum(plan);
k=0;
for i=1:3
    if plan(i) %plan demandé
    k=k+1;
	if i==1
        Ix=1:Nx; Iy=1:Ny; Iz=1;
        titre='Plan xy'; axex='x'; axey='y';
	elseif i==2
    	Ix=(Nx-1)/2; Iy=1:Ny; Iz=1:Nz;
        titre='Plan yz'; axex='z'; axey='y';
	else
    	Ix=1:Nx; Iy=(Ny-1)/2; Iz=1:Nz;
        titre='Plan xz'; axex='z'; axey='x';
	end
	Xs=squeeze(X(Ix,Iz,Iy)); Ys=squeeze(Y(Ix,Iz,Iy)); Zs=squeeze(Z(Ix,Iz,Iy));
	Exs=squeeze(Ex(Ix,Iz,Iy)); Eys=squeeze(Ey(Ix,Iz,Iy)); Ezs=squeeze(Ez(Ix,Iz,Iy)); 
	Hxs=squeeze(Hx(Ix,Iz,Iy)); Hys=squeeze(Hy(Ix,Iz,Iy)); Hzs=squeeze(Hz(Ix,Iz,Iy)); 
	if i==1
    	Exx=Exs; Eyy=Eys; Ezz=Ezs; Hxx=Hxs; Hyy=Hys; Hzz=Hzs; Xx=Xs; Yy=Ys; Zz=Zs;
	elseif i==2
    	Exx=Ezs; Eyy=Eys; Ezz=-Exs; Hxx=Hzs; Hyy=Hys; Hzz=-Hxs; Xx=Zs; Yy=Ys; Zz=-Xs;
	else
    	Exx=Ezs'; Eyy=Exs'; Ezz=Eys'; Hxx=Hzs'; Hyy=Hxs'; Hzz=Hys'; Xx=Zs'; Yy=Xs'; Zz=Ys';
	end
    
	%OK figures : mod et dir (E, H) plan choisi
	cmap=hot;
    figure(1);
	subplot(nplans,1,k), surfaces = pcolor(Xx,Yy,sqrt(Exx.^2+Eyy.^2+Ezz.^2));
    title([titre,',E']); xlabel(axex); ylabel(axey);
	set(surfaces,'FaceColor','interp','EdgeColor','none');
	set(surfaces,'AmbientStrength',.6);
	colormap(cmap(end:-1:8,:));
	axis tight;
	hold on;
    quiver(Xx,Yy,Exx,Eyy,0.8,'k');
    hold off;
    figure(2);
	subplot(nplans,1,k), surfaces = pcolor(Xx,Yy,sqrt(Hxx.^2+Hyy.^2+Hzz.^2));
    title([titre,',H']); xlabel(axex); ylabel(axey);
	set(surfaces,'FaceColor','interp','EdgeColor','none');
	set(surfaces,'AmbientStrength',.6);
	colormap(cmap(end:-1:8,[3 2 1]));
	axis tight;
	hold on;
	quiver(Xx,Yy,Hxx,Hyy,0.8,'k');
	hold off;

%OK figure 3 et 4 : mod et dir (E, H) 3D
%figure(3);
%i=1:2:Nx;
%j=1:1:Ny;
%k=1:1:Nz;
%Cx=X(i,k,j);
%Cy=Y(i,k,j);
%Cz=Z(i,k,j);
%quiver3(Cz,Cx,Cy,Ez(i,k,j),Ex(i,k,j),Ey(i,k,j),0.5,'r');
%axis equal,tight;
%figure(4);
%quiver3(Cz,Cx,Cy,Hz(i,k,j),Hx(i,k,j),Hy(i,k,j),0.5,'b');
%axis equal,tight;

%OK figure 5 : lignes E et H plan choisi
%figure(5);
%xs=[Xx(2:end-1,1)' Xx(2:end-1,1)' Xx(1,2:end-1) Xx(end-1,2:end-1)];
%ys=[Yy(2:end-1,1)' Yy(2:end-1,end-1)' Yy(1,2:end-1) Yy(1,2:end-1)];
%linese=streamline(Xx',Yy',Exx',Eyy',xs,ys);
%set(linese,'Color','red');
%hold on;
%linesh=streamline(Xx',Yy',Hxx',Hyy',xs,ys);
%set(linesh,'Color','blue');
%hold off;

%figure 6 et 7 : lignes (E, H) 3D
%figure(6);
%xs=[X(2,2:end-1,1) X(2,2:end-1,1) X(end-1,2:end-1,1) X(end-1,2:end-1,1)];
%ys=[Y(1,2:end-1,2) Y(1,2:end-1,end-1) Y(1,2:end-1,2) Y(1,2:end-1,end-1)];
%zs=[Z(1,2:end-1,1) Z(1,2:end-1,1) Z(1,2:end-1,1) Z(1,2:end-1,1)];
%xs=[X(2:end-1,1,1)' X(2:end-1,1,1)' squeeze(X(2,1,2:end-1))' squeeze(X(end-1,1,2:end-1))'];
%ys=[Y(2:end-1,1,2)' Y(2:end-1,1,end-1)' squeeze(Y(1,1,2:end-1))' squeeze(Y(1,1,2:end-1))'];
%zs=kron(ones(1,length(xs)),Z(1,2:end-1,1));
%xs=kron(ones(1,length(Z(1,2:end-1,1))),xs);
%ys=kron(ones(1,length(Z(1,2:end-1,1))),ys);
%linese=streamline(Z,X,Y,Ez,Ex,Ey,zs,xs,ys);
%set(linese,'Color','red');
%figure(7);
%linesh=streamline(Z,X,Y,Hz,Hx,Hy,zs,xs,ys);
%set(linesh,'Color','blue');

%figure 8 et 9 : mod E et H 3D coupes divers plans
%cmap=hot;
%figure(8);
%Cx=[a/2];
%Cy=[b/5 4*b/5];
%Cz=[0 d/4 d/2 d];
%surfaces=slice(Z,X,Y,sqrt(Ex.^2+Ey.^2+Ex.^2),Cz,Cx,Cy);
%set(surfaces,'FaceColor','interp','EdgeColor','none');
%set(surfaces,'AmbientStrength',.6);
%colormap(cmap(end:-1:8,:));
%axis equal,tight;
%figure(9);
%surfaces=slice(Z,X,Y,sqrt(Hx.^2+Hy.^2+Hx.^2),Cz,Cx,Cy);
%set(surfaces,'FaceColor','interp','EdgeColor','none');
%set(surfaces,'AmbientStrength',.6);
%colormap(cmap(end:-1:8,[3 2 1]));
%axis equal,tight;

%figure 10 : E ou H 3D avec cones
%figure(10);
%i=1:2:Nx;
%j=1:2:Ny;
%k=1:2:Nz;
%Cx=X(i,k,j);
%Cy=Y(i,k,j);
%Cz=Z(i,k,j);
%hcones = coneplot(Z,X,Y,Ez,Ex,Ey,Cz,Cx,Cy,5);
%set(hcones,'FaceColor','red','EdgeColor','none');
%clf reset;
%hcones = coneplot(Z,X,Y,Hz,Hx,Hy,Cz,Cx,Cy,5);
%set(hcones,'FaceColor','blue','EdgeColor','none');
%axis equal,tight; view(30,40);
%camproj perspective; camzoom(1.5);
%camlight right; lighting phong;
%set(hcones,'DiffuseStrength',.8);
	end
end
end
