format long
syms p;
kh=0.01;
cd=6100;cs=3300;m=cs/cd;
y1=-4*p*((p^2+1)*m^2-1)^(1/2)*tan(kh/2*p);
y2=(p^2-1)^2*tan(kh/2*((p^2+1)*m^2-1)^(1/2));
y=y1-y2;
ezplot(y1,[5000,10000]);
grid on
hold on
pause
ezplot(y2,[5000,10000]);
zoom on
pause;
[xx,yy]=ginput(3);  %從圖上獲取曲線零點所對應的p值點
c1=sqrt(xx(1)^2+1)  %求出所對應的c/cs值
c2=sqrt(xx(2)^2+1)
c3=sqrt(xx(3)^2+1)
[x(1),fval,exitflag]=fzero(inline(y1-y2),xx(1),[]);
[x(2),fval,exitflag]=fzero(inline(y1-y2),xx(2),[]);
[x(3),fval,exitflag]=fzero(inline(y1-y2),xx(3),[]);
c11=sqrt(x(1)^2+1)
c22=sqrt(x(2)^2+1)
c33=sqrt(x(3)^2+1)
%檢驗結果
p1=sqrt(c11^2-1);
y11=-4*p1*((p1^2+1)*m^2-1)^(1/2)*tan(kh/2*p1);
y21=(p1^2-1)^2*tan(kh/2*((p1^2+1)*m^2-1)^(1/2));
y01=y11-y21
p2=sqrt(c22^2-1);
y12=-4*p2*((p2^2+1)*m^2-1)^(1/2)*tan(kh/2*p2);
y22=(p2^2-1)^2*tan(kh/2*((p2^2+1)*m^2-1)^(1/2));
y02=y12-y22
p3=sqrt(c33^2-1);
y13=-4*p3*((p3^2+1)*m^2-1)^(1/2)*tan(kh/2*p3);
y23=(p3^2-1)^2*tan(kh/2*((p3^2+1)*m^2-1)^(1/2));
y03=y13-y23
%作圖
%a=jieguo01;x=a(:,1);y0=a(:,2);y1=a(:,3);y2=a(:,4);y3=a(:,5);y4=a(:,6);
%y5=a(:,7);y6=a(:,8);y7=a(:,9);plot(x,y0); hold on;grid on;
%plot(x,y1);plot(x,y2);plot(x,y3);plot(x,y4);plot(x,y5);plot(x,y6);plot(x,y7);