% AR_PZ_SURF.m
%
% P. Flandrin, Feb. 15, 2006
%
% calcule et trace le module de la transformee en z
% d'un modele AR(2) reel pour differentes valeurs
% du module des poles (complexes conjugues)

[th,r] = meshgrid((0:1:360)*pi/180,0:.025:1) ;
[X,Y] = pol2cart(th,r) ;
[th2,r2] = meshgrid((0:5:360)*pi/180,0:.1:1) ;
[X2,Y2] = pol2cart(th2,r2) ;
[th3,r3] = meshgrid((0:5:360)*pi/180,0:1:1) ;
[X3,Y3] = pol2cart(th3,r3) ;
Z = X+i*Y ; Z2 = X2+i*Y2 ;

M = 25 ;

p = [.5 .6 .7 .75 .8 .85 .9 .94 .96 .98 .99] ; % modules des poles

for k = 1: length(p)
    
    pause
    f = Z.^2./(Z+(1+i)*p(k)/sqrt(2))./(Z+(-1+i)*p(k)/sqrt(2)) ; % module de la TZ
    surf(X,Y,min(M+abs(f),2*M)) ; view(-15,30) ; colormap(copper)
    hold on
    mesh(X2,Y2,zeros(size(X2)))
    mesh(X3,Y3,M+zeros(size(X3)))
    plot3([-p(k)/sqrt(2) -p(k)/sqrt(2)],[-p(k)/sqrt(2) -p(k)/sqrt(2)],[0 2*M],'r','LineWidth',6)
    plot3([-1/sqrt(2) 0],[-1/sqrt(2) 0],[0 0],'b','LineWidth',6)
    plot3([p(k)/sqrt(2) p(k)/sqrt(2)],[-p(k)/sqrt(2) -p(k)/sqrt(2)],[0 2*M],'r','LineWidth',6)
    plot3([1/sqrt(2) 0],[-1/sqrt(2) 0],[0 0],'b','LineWidth',6)
    plot3([-p(k)/sqrt(2) -p(k)/sqrt(2)],[-p(k)/sqrt(2) -p(k)/sqrt(2)],[0 0],'ro','LineWidth',6)
    plot3([p(k)/sqrt(2) p(k)/sqrt(2)],[-p(k)/sqrt(2) -p(k)/sqrt(2)],[0 0],'ro','LineWidth',6)
    plot3([-p(k)/sqrt(2) -p(k)/sqrt(2)],[-p(k)/sqrt(2) -p(k)/sqrt(2)],[M M],'ro','LineWidth',6)
    plot3([p(k)/sqrt(2) p(k)/sqrt(2)],[-p(k)/sqrt(2) -p(k)/sqrt(2)],[M M],'ro','LineWidth',6)
    plot3([0 0],[0 0],[0 2*M],'g','LineWidth',6)
    plot3([0 0],[1 -1],[0 0],'g','LineWidth',6)
    plot3([-1 1],[0 0],[0 0],'g','LineWidth',6)
    set(gca,'ZTick',[])
    text(0,1,fix(7*M/4),['module = ',num2str(p(k))],'FontSize',14)
    xlabel('partie reelle')
    ylabel('partie imaginaire')
    zlabel('module de la TZ')
    hold off
    
end