exercise 9.10

%%  mue : Effective Permeability

MUe = mu0*mur*(1+fm*f./(F0.^2-f^2));    % (9.25a)

%% B contains all betas for each t/a
B=zeros(2,numel(H0)+1,numel(N3));


for s2=1:1:numel(N3)

% n4 low means thin ferrite, n4 large means ferrite completely
% filling up  waveguide cross-section.

   t1=t(s2);
   c1=c(s2);
   d1=d(s2);   % c=d imposed, from question header but leave d in case future use

   beta=Beta_range;  % it was syms but didn't work  

   Beta=[t1;c1];

   for s1=1:1:numel(H0)

      Kfsq=4*pi*f^2*MUe(s1)*eps0*er-beta.^2;
      Kasq=4*pi^2/c0^2*f^2-beta.^2;

% building equation
      Eq1=Kfsq./(MUe(s1)).^2;                                                          % [1]

      Eq2=K2(s1).^2.*beta.^2./((mu0*mur)^2*MUe(s1).^2);             % [2]

      KEq3=abs(Kasq.^.5).*cot(abs(Kasq.^.5)*c1)./(mu0*MUe(s1));    % [3]
      Eq31=KEq3.*abs(Kfsq.^.5).*cot(abs(Kfsq.^.5)*t1);
      Eq32=KEq3.*K2(s1).*beta/mur;

      Eq4=Kasq/mu0^2;                                                                    % [4]

      Eq5=cot(abs(Kasq.^.5)*c1).*tan(abs(Kasq.^.5)*d1);                  % [5]

      KEq6=abs(Kasq.^.5).*tan(abs(Kasq.^.5)*d1)./(mu0*MUe(s1));   % [6]
      Eq61=KEq6.*abs(Kfsq.^.5).*cot(abs(Kfsq.^.5)*t1);
      Eq62=-KEq6.*K2(s1).*beta./mur;

% assembling equation

      S1vor=Eq1+Eq2+Eq31+Eq32+Eq4+Eq5+Eq61+Eq62;

%% solving equation
% When plotting the following
% figure;plot(beta,S1vor);grid on;xlabel('\beta');title('S1vor')
% the asymptote is the sought point [beta(n) H0(n)].
% To capture this point support function intersections.m is particularly handy
% directly using the Y points generated by plot function 

      hf1=figure;
      hf1.Visible='off';
      hp1=plot(beta,S1vor);grid on;xlabel('\beta');title('S1vor');

% capturing asymptote
      [x0,y0]=intersections([beta(1) beta(end)],[0 0],beta,hp1.YData,'robust');

      if numel(x0)<2
          x0=[x0;0];
     end

      if numel(x0)>2
          x0([3:end])=[];
      end

%% logging single point result
      Beta=[Beta x0];

      close all;
% Not closing so many figures causes forced shut down that in turn
% shows up as the following low level error on launching MATLAB again:
%
% Warning: MATLAB previously crashed due to a low-level graphics error. 
% To prevent another crash in this session, MATLAB is using software
% OpenGL instead of using your graphics
% hardware. To save this setting for future sessions, use the opengl('save', 'software') command. 
% For more information, see Resolving Low-Level Graphics Issues. 
% > In matlab.graphics.internal.initialize (line 15) 


      end

% Beta(:,1)=[]; % 1st round left log header empty, but when calculating for
% different [t c] I am putting [t(s2);c(s2)] in each Beta header, to be
% able to tell what fill-out corresponds to each pair of curves.


%% Displaying beta(H0) for particular [t c] fill-up

   B(:,:,s2)=Beta;

   b1=Beta(1,[2:end]);
   b2=Beta(2,[2:end]);
   hf1=figure;plot(H0,b1/k0,H0,b2/k0)
   ax2=gca;
   grid on;
   grid minor;
   ax2.PlotBoxAspectRatioMode='manual';
   ax2.PlotBoxAspectRatio=[1 1 1];
   grid on;
   title('beta(H0)/k0');
   str1=[ ' t/a = ' num2str(t1/a) ];
   text(300,2.4,str1,'FontSize',14)
   xlabel('H_0');ylabel('\beta/H_0');

   legend({' \beta1','\beta2'},'Location','northeastoutside')

   str2=['graph_' num2str(s2) '_chap_09_ex_10']
   saveas(hf1,str2,'jpg')
        
end