clear all

Omega=linspace(0.4+0.001*sqrt(2),3.5,2*1801);
EP=linspace(0,3,2001);

for ie=1:length(EP);
    e=EP(ie);
    for io=1:length(Omega);
        om=Omega(io);
        o1=om+e;
        o2=om-e;
        c1=cos(pi*o1);s1=sin(pi*o1);
        c2=cos(pi*o2);s2=sin(pi*o2);
        PHI1=[c1, s1/o1;-o1*s1,c1];
        PHI2=[c2, s2/o2;-o2*s2,c2];
        PHI=PHI1*PHI2;
        sigma=max(abs(eig(PHI)));
        if sigma<1.0000001;
            lambda(io,ie)=-10;
        else;
        lambda(io,ie)=log(sigma)/(2*pi);
        end; 
    end;
end;
    v=linspace(0.0001,2,41);
    figure(32)
    contourf(Omega,EP,lambda',v);colormap jet;colorbar
    xlabel('\omega');ylabel('\epsilon')
    xlim([0 3]);ylim([0 3])
    axis('square')
        o1=2;t1=2*pi/(4*o1);
        o2=1/2;t2=2*pi/(4*o2);
        
        c1=cos(t1*o1);s1=sin(t1*o1);
        c2=cos(t2*o2);s2=sin(t2*o2);
        PHI1=[c1, s1/o1;-o1*s1,c1];
        PHI2=[c2, s2/o2;-o2*s2,c2];
        PHI=PHI1*PHI2;
        sigma=(abs(eig(PHI)));