% function y=imp1(x)      %thesis pic 1,the relationship of w and beta
%  n=[1 3 1.5];
%  l1=0.13;
%  l2=0.265;
% %  xx=5:2:30;
% % for   w=xx
% global w
%      k=sqrt((n.^2)*(w.^2) -x^2) ;% w=w/c
%      A=exp(i*k(2)*l1)*(i*sin(k(3)*l2)*(k(2)^2+k(3)^2)/(2*k(2)*k(3))+cos(k(3)*l2));
%      e=real(A)+((real(A))^2-1)^0.5;
%      B=i*exp(-i*k(2)*l1)*sin(k(3)*l2)*(k(2)^2-k(3)^2)/(2*k(2)*k(3));
%      y=k(2)*besselj(k(1),-1)*(e-A+B)+i*k(1)*besselj(k(1),0)*(e-A-B);
% 
%      %  end   

% % %      f=inline(['k(2)*besselj(k(1),-1)*(e-A+B)+i*k(1)*besselj(k(1),0)*(e-A-B)'],'beta')
% % %    ezplot(' k(2)*besselj(k(1),-1)*(e+A+B)+i*k(1)*besselj(k(1),0)*(e-A-B)')


% % % y=[];xx=0:.05:10;x0=0;h=optimset;h.Display='off';    %循環畫圖例子177
% % % for x=xx
% % %     f=inline(['w.*exp(w)-'num2str(x)],'w');
% % %     y1=fsolve(f,x0,h);x0=y1;y=[y,y1];
% % % end
% % % plot(xx,y)


% function y=imp1(x)                         %write imp&imp1 in the same file cannot get the answer
%  n=[1 3 1.5];
%  l1=0.13;
%  l2=0.265;
%  y=[];xx=5:.05:33;x0=0.3;h=optimset;h.Display='off';
%  for w=xx
%      k=sqrt((n.^2)*(w.^2) -x^2) ;% w=w/c
%      A=exp(i*k(2)*l1)*(i*sin(k(3)*l2)*(k(2)^2+k(3)^2)/(2*k(2)*k(3))+cos(k(3)*l2));
%      e=real(A)+((real(A))^2-1)^0.5;
%      B=i*exp(-i*k(2)*l1)*sin(k(3)*l2)*(k(2)^2-k(3)^2)/(2*k(2)*k(3));
%      p=k(2)*besselj(k(1),-1)*(e-A+B)+i*k(1)*besselj(k(1),0)*(e-A-B);
%      y1=fsolve(p,x0,h);x0=y1;y=[y,y1];
%  end
%  plot(xx,y)
 

function y=imp1(x) %write imp&imp1 in the same file cannot get the answer
% y=x^2+1                   %this example show fsolve can not be used in compex answer
 syms x;
 n=[1 3 1.5];
 l1=0.13;
 l2=0.265;
     k=sqrt((n.^2)*(6.^2) -x^2) ;%w=w/c
     A=exp(i*k(2)*l1)*(i*sin(k(3)*l2)*(k(2)^2+k(3)^2)/(2*k(2)*k(3))+cos(k(3)*l2));
     e=real(A)+((real(A))^2-1)^0.5;
     B=i*exp(-i*k(2)*l1)*sin(k(3)*l2)*(k(2)^2-k(3)^2)/(2*k(2)*k(3));
     y=[k(2)*besselj(k(1),-1)*(e-A+B)+i*k(1)*besselj(k(1),0)*(e-A-B)];  
% y=0.5;

% function y=imp1(x)  % as an example,this can run well
% syms x k;      %用solve時必須注明這個！
% k=x+2;
% a=k*2;
% y=cos(a*x)+x;

%%%solve函數和fsolve語言上不同。
% S=solve( 'k(2)*besselj(k(1),-1)*(e-A+B)+i*k(1)*besselj(k(1),0)*(e-A-B)=0','k=sqrt(x^2)',...  %try to use 'solve'
%     'A=exp(i*k(2)*l1)*(i*sin(k(3)*l2)*(k(2)^2+k(3)^2)/(2*k(2)*k(3))+cos(k(3)*l2))',...
%     ' e=real(A)+((real(A))^2-1)^0.5',...
%     'B=i*exp(-i*k(2)*l1)*sin(k(3)*l2)*(k(2)^2-k(3)^2)/(2*k(2)*k(3))')
% %%用solve時顯示explicit solution could not be found