
format compact

b=4;
A=100;

figure(1);clf;
figure(2);clf;
figure(3);clf;

lw2=2;

w=linspace(1,10,90*2+1);
wp=linspace(1,10,90*2+1);
[W,WP]=meshgrid(w,wp);
wop=real(sqrt(w.^2-2*b^2));
Dd=(w.^2-w.^2).^2+4*(b^2)*w.^2;
Dop=(w.^2-wop.^2).^2+4*(b^2)*wop.^2;
D=(W.^2-WP.^2).^2+4*(b^2)*WP.^2;

figure(1);
surf(w,wp,log(1./(D+A)));
hold on;
jjj=96:length(w);
plot3(w(jjj),wop(jjj),(log(1./(Dop(jjj)+A))),'r','LineWidth',lw2);
plot3(w,w,log(1./(Dd+A)),'LineWidth',lw2);
plot3([1 1]*8,[1 10],-5*[1 1],'r--');
plot3([1 10],[1 1]*8,-5*[1 1],'b--');
hold off;
xlabel('\omega');ylabel('\omega''');view(2);shading interp;
axis([min(w) max(w) min(wp) max(wp)])

figure(3);
jj=141;
plot(w,log(1./D(:,jj)),'r');
hold on
m1=min(log(1./D(:,jj)));
m1=-9;
M1=max(log(1./D(:,jj)));
plot(wop(jj)*[1 1],[m1 M1],'k');
hold off
title(['\omega=',num2str(wp(jj))])
axis([1 10 -9.0 -8.0])

figure(2);
ii=141;
plot(wp,log(1./D(ii,:)));
hold on
m2=min(log(1./D(ii,:)));
m2=-9;
M2=max(log(1./D(ii,:)));
plot(w(ii)*[1 1],[m2 M2],'k');
hold off
title(['\omega''=',num2str(w(ii))])
axis([1 10 -9.0 -8.0])

