function JSV1984PlusAmp(inputfile)       % S.W. Rienstra 16-08-2021, based on "Acoustic Radiation From A Semi-Infinite Annular Duct In A Uniform Subsonic Mean Flow", J. Sound Vibration 1984, 94(2), 267-288
                                         % used to make the figures in "Acoustic Energy Balances for Sound Radiated from Duct Exit with Mean Flow", International Journal of Aeroacoustics, Aug. 2022.
   %clear all;
   opengl('save', 'software')
   warning off MATLAB:MKDIR:DirectoryExists
   format long
   format compact
   casenm = 'FW';                 % first part of name
   zwart  = [0, 0, 0];            % black
   blauw  = [0, 0, 1];            % blue
   rood   = [1, 0, 0];            % red
   paars  = [0.75, 0, 0.75];      % purple
   bruin  = [0.85, 0.325, 0.098]; % brown
   dkgeel = [0.929, 0.694, 0.125];% dark yellow
   groen  = [0, 0.5, 0];          % green
   lgroen = [0.466, 0.674, 0.188];% lighter green
   pprsze = [0, 0, 6, 4];         % paperpositionsize 6x4 cm
   pltfnt = 6;                    % plotaxis fontsize
   fnext  = '.eps';               % filename extension eps
   tstart = tic;
   if nargin == 1
      display(['inputfile = ', inputfile]);
      fid = fopen(inputfile,'r');
      eval(fscanf(fid,'%c'));
      fclose(fid);
   else
    %% input parameters
    test     = 0;   % =1: plot integration contours
    pwlbl    = 0;   % =0: % minimize labels in power plots
    scenar   = 1;   % scenario 1: omar=[om1,om2],      R    , Pt,     D(xi),    Pf,     PHi,     PHo,     Pw,     plotD, savedat, Kuttacmp, powD
                    % scenario 2: omar=[om1:0.01:om2], R(om), Pt(om), D(om,xi), Pf(om), PHi(om), PHo(om), Pw(om), plotR, plotP,   Kuttacmp, powD
                    % scenario 3: omar=[om1:0.01:om2], R(om), Pt(om),           Pf(om), PHi(om), PHo(om), Pw(om), plotR, plotP,   Kuttacmp         Note: Pf is determined indirectly.
                    % scenario 4: omar=[om1,om2],       Am,   Pt,     D(xi),    Pf,     PHi,     PHo,     Pw,     plotD, savedat,           powD
                    % scenario 5: omar=[om1:0.01:om2],  Am,   Pt(om), D(om,xi), Pf(om), PHi(om), PHo(om), Pw(om), plotD, plotP,             powD
                    % scenario 6: omar=[om1:0.01:om2],  Am,   Pt(om),           Pf(om), PHi(om), PHo(om), Pw(om),        plotP                     Note: Pf is determined indirectly.
    KuttaCmp = 0;   % compares gm=1 with gm=0 in 1 picture per mode
    %
    h     =  0.90;   % hub
    M     =  0.20;   % Mach number
    gm    =  1.0;    % Kutta parameter (unless M=0 or KuttaCmp =1)
    %
    switch scenar
     case 1,     omar  = [0.001                                               ];                 % typical frequency for D-plot
     case {2,3}, omar  = [.001:0.0010:10. ];                 % typical frequency for R-plot and powerplot
     case 4,     omar  = [ 3              ];                 % typical frequency for D-plot
     case {5,6}, omar  = [0.010:0.010:10. ];                 % typical frequency for powerplot
    end
    omfln =   10;                 % omfn  if scenar = 1 or 4
    Nx    =  100;                 % number of xi-angles per lobe (because of GL integration)
    m     =   0;                  % circumferential order
    mumx  =   6;                  % number of alpha's; should be at least muco, nr of cut-on modes + 1; scenar=4: also number of amplitudes
    mudir =   1;                  % only with scenar = 1: mudir=0: plot dirs of all mu; mudir=n plot n-th mode (not efficient: all modes are calculated)
    muR   =  [1,1,2,2];           % plot R_{muR,nuR}
    nuR   =  [1,2,1,2];           %
    klrR  =  {zwart, blauw, rood, paars, bruin, dkgeel};   % colors R plot
    klrD  =  {zwart, blauw, rood, paars, bruin, dkgeel};   % colors D plot
    klPt  =  zwart;
    klPf  =  rood;
    klPHi =  lgroen;
    klPHo =  dkgeel;
    klPw  =  blauw;
    shph  =  [-1, 0,-1, 0]*2*pi;     % shift of phi c.q. phi1 if KuttaCmp=1
    shph0 =  [-1,-1,-1,-1]*2*pi;     % shift of phi0 if KuttaCmp=1
    Edif  =   65;                 % minimal difference in dB between Emax and Emin in polar plot, plotD
    d     =  0.1 ;                % parameters of deformed integration contour
    q     =  6  ;                 % typically 4 or 6
    %% end input parameters
   end
   Nom   =  length(omar);

   if ismember(scenar,[1,4])
      omfn = omfln    ;            % file name is based on omfln
   else
      omfn = omar(Nom);            % file name is based on max(omega)
   end
   if ismember(scenar,[4,5,6])
      Am   = ones(Nom,mumx);
     %Am(1:Nom,1) = 1 ;           % non-zero modal amplitudes
     %Am(1:Nom,2) = 1 ;
     %Am(1:Nom,3) = 1 ;
     %Am(1:Nom,4) = 1 ;
     %Am(1:Nom,5) = 1 ;
     %Am(1:Nom,6) = 1 ;
     %Am(1:Nom,7) = 1 ;
     %Am(1:Nom,8) = 0 ;
      casenm = strcat(casenm(end-1),'A');
   end
   if KuttaCmp && ismember(scenar,[4,5,6])
      display(['scenar = 4,5,6 and KuttaCmp not compatible']); return;
   end
   NR   = length(muR);          % number of R-plots in one picture
   bt   = sqrt(1-M^2);          % beta
   %
   if ~ismember(scenar,[4,5,6])
     Pt   = zeros(mumx,Nom);
     Px   = zeros(mumx,Nom);
     Pf   = zeros(mumx,Nom);
     PHi  = zeros(mumx,Nom);
     PHo  = zeros(mumx,Nom);
     Pw   = zeros(mumx,Nom);
     Pt0  = zeros(mumx,Nom);
     Px0  = zeros(mumx,Nom);
     Pf0  = zeros(mumx,Nom);
     PHi0 = zeros(mumx,Nom);
     PHo0 = zeros(mumx,Nom);
     Pw0  = zeros(mumx,Nom);
     Pt1  = zeros(mumx,Nom);
     Px1  = zeros(mumx,Nom);
     Pf1  = zeros(mumx,Nom);
     PHi1 = zeros(mumx,Nom);
     PHo1 = zeros(mumx,Nom);
     Pw1  = zeros(mumx,Nom);
   else
     PtA  = zeros(1,Nom);
     PfA  = zeros(1,Nom);
     PHiA = zeros(1,Nom);
     PHoA = zeros(1,Nom);
     PwA  = zeros(1,Nom);
     Bm   = zeros(Nom,mumx);    % reflected modal amplitudes
   end
   plotP   = 0;
   plotR   = 0;
   plotD   = 0;
   powD    = 0;
   savedat = 0;
   switch scenar
    case 1, plotD=1; powD =1; savedat= 1 ; if Nom >5, display(['too big omar'])  ; end
    case 2, plotR=1; plotP=1; powD   = 1 ; if Nom <5, display(['too small omar']); end
    case 3, plotR=1; plotP=1;              if Nom <5, display(['too small omar']); end
    case 4, plotD=1; powD =1; savedat= 1 ; if Nom >5, display(['too big omar'])  ; end
    case 5,          plotP=1; powD = 1   ; if Nom <5, display(['too small omar']); end
    case 6,          plotP=1;              if Nom <5, display(['too small omar']); end
    otherwise, display(['wrong scenario']); return;
   end
   % savedat  : save R-data en D-data
   % plotR    : plot R(omega)
   % plotD    : plot D(xi)
   % plotP    : plot Pt(omega), Pf(omega), PHo(omega), PHi(omega), Pw(omega)
   % powD     : Pf(omega) directly determined in scenar =1, 2, 4, 5. Indirectly if =3, 6.
   %            plotD implies powD directly; not-powD implies powD from other powers & energy conservation.
   %
   % not used here, but in case a fixed xi-distribution is required: xi = repmat(pi*[0:Nx-1]./(Nx-1),Nom,1);
   %
   %%%%%%%%%%%%%%
   % Boundingbox correction
   dx1 = -5;
   dx2 =  5;
   dy1 = -5;
   dy2 =  5;
   %%%%%%%%%%%%%%
   % max error polyline simplification
   dperr = 1e-6;
   %%%%%%%%%%%%%%

   %%%%%%
   if h==0, hstr = 'h0.0'; else hstr = num2str(h,'h%0.6g'); end;
   if M==0, Mstr = 'M0.0'; else Mstr = num2str(M,'M%0.6g'); end;
   switch myceil(omfn)
   case  3,    casenm = strcat(casenm,'w',num2str(omfn,'%3.0f'), 'm',num2str(m,'%1i'),hstr,Mstr);
   case  2,    casenm = strcat(casenm,'w',num2str(omfn,'%2.0f.'),'m',num2str(m,'%1i'),hstr,Mstr);
   case  1,    casenm = strcat(casenm,'w',num2str(omfn,'%1.1f'), 'm',num2str(m,'%1i'),hstr,Mstr);
   case  0,    casenm = strcat(casenm,'w',num2str(omfn,'%0.1f'), 'm',num2str(m,'%1i'),hstr,Mstr);
   case -1,    casenm = strcat(casenm,'w',num2str(omfn,'%0.2f'), 'm',num2str(m,'%1i'),hstr,Mstr);
   case -2,    casenm = strcat(casenm,'w',num2str(omfn,'%0.2f'), 'm',num2str(m,'%1i'),hstr,Mstr);
   case -3,    casenm = strcat(casenm,'w',num2str(omfn,'%0.3f'), 'm',num2str(m,'%1i'),hstr,Mstr);
   otherwise   casenm = strcat(casenm,'w',num2str(omfn,'%0.4f'), 'm',num2str(m,'%1i'),hstr,Mstr);
   end
   fprintf('case = %s\n',casenm);

   dirpad  = strcat('.\data\',casenm,'\');
   dirpad1 = strcat('.\\data\\',casenm,'\\');
   dirpad2 = strcat('.\\\\data\\\\',casenm,'\\\\');
   mkdir(dirpad);

   epsR      = strcat(dirpad,casenm,'R',fnext) ;
   epsphi    = strcat(dirpad,casenm,'phi',fnext);
   intcont   = strcat(dirpad,casenm,'intcnt',fnext);
   for j=1:Nom
     if ~ismember(scenar,[4,5,6])
       epsEmmu(j) = string(strcat(dirpad,casenm,'Emmu',num2str(j,'%1i'),fnext)) ;
     else
       epsEmmu(j) = string(strcat(dirpad,casenm,'EmAm',num2str(j,'%1i'),fnext)) ;
     end
   end
   if ~ismember(scenar,[4,5,6])
    for mu=1:mumx
       epsP(mu) = string(strcat(dirpad,casenm,'P',num2str(mu,'%1i'),fnext)) ;
    end
   else
       epsPA    = string(strcat(dirpad,casenm,'PA',fnext)) ;
   end
   Rdatfile  = strcat(dirpad,casenm,'R.dat') ;
   if ~ismember(scenar,[4,5,6])
     Ddatfile  = strcat(dirpad,casenm,'Dmu.dat') ;
   else
     Ddatfile  = strcat(dirpad,casenm,'DAm.dat') ;
   end
   Pdatfile  = strcat(dirpad,casenm,'P.dat') ;
   inpfile   = strcat(dirpad,casenm,'.input') ;
   %%%%%%
   % determine zeros
   % plot(aj0,BesEig(m,aj0,h),'-k',alph,zeros(1,mumx),'xr')
   % axis([-0.01*alph(mumx),alph(mumx)*1.01,-1,1])
   %
   alph  = findalphas(m,h,mumx);          % array of alpha(m,mu), length mumx
    Jm   =  J(m,alph);
   dJm   = dJ(m,alph);
   % U_{m mu}(1)
   if h ==0
     if m ~=0
        Nm = 1./sqrt((1-(m./alph).^2)/2.*Jm.^2);
     else
        Nm = [sqrt(2),1./sqrt(Jm(2:mumx).^2/2)];
     end
     Um = Nm.*Jm;
   else
     alphh = alph*h;
     dJmh  = dJ(m,alphh);
     if m ~=0
        Ym   =  Y(m,alph );
       dYm   = dY(m,alph );
       dYmh  = dY(m,alphh);
       Nm    = alph*pi/sqrt(2)./sqrt((1-(m./alph).^2)./(dJm.^2+dYm.^2) - (1-(m./alphh).^2)./(dJmh.^2+dYmh.^2));
       taum  = atan(dJm./dYm);
       Um    = Nm.*(cos(taum).*Jm-sin(taum).*Ym);
     else
       Ym    = [0, Y(m,alph(2:mumx) )];
       dYm   = [0,dY(m,alph(2:mumx) )];
       dYmh  = [0,dY(m,alphh(2:mumx))];
       Nm    = [sqrt(2/(1-h^2)),alph(2:mumx)*pi/sqrt(2)./sqrt(1./(dJm(2:mumx).^2+dYm(2:mumx).^2) - 1./(dJmh(2:mumx).^2+dYmh(2:mumx).^2))];
       taum  = [0,atan(dJm(2:mumx)./dYm(2:mumx))];
       Um    = [sqrt(2/(1-h^2)),Nm(2:mumx).*(cos(taum(2:mumx)).*Jm(2:mumx)-sin(taum(2:mumx)).*Ym(2:mumx))];
     end
   end
   %
   if isempty(find(bt*alph/max(omar)>1,1,'first'))                   % all modes are cut-on, no first cut-off
       display([sprintf('increase mumx: mumx= %1i < muco',mumx)])    % check if mumx is big enough compared to muco
       return;
   else
       display([sprintf('muco = %1i ... %1i',find(bt*alph/omar(1)>1,1,'first'),find(bt*alph/omar(Nom)>1,1,'first'))]);
   end
   if test
     testKm(M,d,q,intcont,dx1,dy1,dx2,dy2,pprsze,pltfnt);
   end
   %
   % solution per j-th frequency
   for j = 1:Nom
     omega = omar(j);
     display([sprintf('omega = %8.15f', omega)])
     kappa = omega/bt;
     Ombar = vomb(alph/kappa);                 % array of \bar{\Omega}_{m,mu}, length mumx. Note: keep omega>0, because Ombar -> -i*infty if omega=0, except if m=0,mu=1
     mco   = find(alph/kappa>1,1,'first');   % mco-1 = number of cut-on modes for omega(j)
     if mco == 1 && mudir==0
        display('no cut-on modes. mudir := 1')
        mudir = 1;
     end

     if M ~= 0                             % prepare some functions of argument 1/M
        MM   = 1/M +i*1e-16;
        KtlM = Ktilde(MM,kappa,m,h);
        KpM  = Ksplit( 1,MM,kappa,m,h,1/M,d,q);
        KmM  = Ksplit(-1,MM,kappa,m,h,1/M,d,q);
        dKpM = Ksplit( 2,MM,kappa,m,h,1/M,d,q);
        StM  = omega/M;
        dk   = dK(m,StM);
        di   = dI(m,StM);
        iM   = I(m,StM);
        kM   = K(m,StM);
        if h ~=0
           Sth  = omega*h/M;
           dkh  = dK(m,Sth);
           dih  = dI(m,Sth);
           dki  = dkh*di-dih*dk;
        end
     end
     % Ktldmin
     for mu=1:mumx
         Ktldmin(mu) = Ksplit(-1,Ombar(mu),kappa,m,h,1,d,q);
         if M~= 0
           Cmmu(mu)    = i/2*(1-M)*Um(mu)*Ktldmin(mu)*KpM*(1+Ombar(mu))/(1-M*Ombar(mu));
         end
     end
     %% reflection matrix
     if KuttaCmp & M~=0
       for mu=1:mumx
         for nu=1:mumx
           R0 = -1/4*Um(mu)*Um(nu);
           R0 = R0*(1+Ombar(mu))*(1+Ombar(nu))/Ombar(nu)/(Ombar(mu)+Ombar(nu));
           R0 = R0*Ktldmin(mu)*Ktldmin(nu);
           R1 = R0;
           R0 = R0*(1+M*Ombar(nu))/(1-M*Ombar(mu));
           Rm0(j,mu,nu) = R0;
           Rm1(j,mu,nu) = R1;
         end
       end
     else
        for mu=1:mumx
          for nu=1:mumx
            R0 = -1/4*Um(mu)*Um(nu);
            R0 = R0*(1+Ombar(mu))*(1+Ombar(nu))/Ombar(nu)/(Ombar(mu)+Ombar(nu));
            R0 = R0*(1-M*gm*Ombar(mu)-M*(gm-1)*Ombar(nu))/(1-M*Ombar(mu));
           %if mu==1 && nu==1
           %display([sprintf('%+8.15f%+8.15fi',real(Ktldmin(mu)),imag(Ktldmin(mu))),',  ',sprintf('%8.15f%+8.15fi',real(Ombar(mu)),imag(Ombar(mu)))])
           %display([sqrt(2)/Um(mu)])
           %end
            R0 = R0*Ktldmin(mu)*Ktldmin(nu);
            Rm(j,mu,nu) = R0;
          end
        end
        if ismember(scenar,[4,5,6])
          PwB  = 0;
          for nu=1:mumx
            for mu=1:mumx
              Bm(j,nu) = Bm(j,nu)+Am(j,mu)*Rm(j,mu,nu);
            end
          end
        end
     end
     %
     if powD % angles are necessary for plotD and for Pf
       % directional function make angles based on lobes. Each side lobe has 1+Nx+1 points, main lobe (1+Nx+1)+(1+Nx+1)
       xii(1)    = 0;
       if mco == 1+[m==0]
        xii(Nx+2)     = pi;
        [xiaux, waux] = legzo(Nx,0,pi);
        xic(1,1:Nx)   = fliplr(xiaux);
        wic(1,1:Nx)   = fliplr(waux);
        for jj = 1:Nx  %1:Nx
           xii(1+jj) = xic(1,jj);
        end
       else
        for nc=1:mco-1-[m==0]
          xii(nc*(Nx+1)+1) = acos(Ombar(nc+[m==0]));
          [xiaux, waux ]   = legzo(Nx,xii((nc-1)*(Nx+1)+1),xii(nc*(Nx+1)+1));
          xic(nc,1:Nx)     = fliplr(xiaux);
          wic(nc,1:Nx)     = fliplr(waux);
          for jj = 1:Nx
              xii((nc-1)*(Nx+1)+1+jj) = xic(nc,jj);
          end
        end
        xii((mco-[m==0])*(Nx+1)+1) = pi;
        [xiaux, waux ]       = legzo(Nx,xii((mco-1-[m==0])*(Nx+1)+1),xii((mco-[m==0])*(Nx+1)+1));
        xic(mco-[m==0],1:Nx) = fliplr(xiaux);
        wic(mco-[m==0],1:Nx) = fliplr(waux);
        for jj = 1:Nx
          xii((mco-1-[m==0])*(Nx+1)+1+jj) = xic(mco-[m==0],jj);
        end
       end
       NNx(j) = length(xii(:));  % total length of array xi at frequency j = (nr.lobes x (Nx+1)) + 1
     end
     %
     if plotD  % prepare for plot
        muco(j) = mco;
        xi( j,1:NNx(j)) = xii(1:NNx(j));
        xib(j,1:NNx(j)) = acos(bt*cos(xii(1:NNx(j)))./sqrt(1-(M*cos(xii(1:NNx(j)))).^2));  % physical directions, corresponding with xi
        E0(1:NNx(j))    = sqrt(1-(M*cos(xii(1:NNx(j)))).^2);  % factor before D because of transformation from xi to xib
     end
     %
     for mu=1:mumx
       % power Pt, PHi, PHo, Pw
       if KuttaCmp & M~=0
        if mu>=mco
          Ommu = abs(Ombar(mu));
          Pt0(mu,j) = -2*pi*bt^4*Ommu/(1+M^2*Ommu^2)^2;
          Pt1(mu,j) = -2*pi*bt^4*Ommu/(1+M^2*Ommu^2)^2;
          Pt0(mu,j) = Pt0(mu,j)*(imag(Rm0(j,mu,mu))*(1-M^2*Ommu^2)+real(Rm0(j,mu,mu))*2*M*Ommu);
          Pt1(mu,j) = Pt1(mu,j)*(imag(Rm1(j,mu,mu))*(1-M^2*Ommu^2)+real(Rm1(j,mu,mu))*2*M*Ommu);
        else
          Pt0(mu,j) = pi*bt^4*Ombar(mu)/(1-M*Ombar(mu))^2;
          Pt1(mu,j) = Pt0(mu,j);
        end
        for nu=1:mco-1
           Pt0(mu,j) = Pt0(mu,j) - pi*bt^4*abs(Rm0(j,mu,nu))^2*Ombar(nu)/(1+M*Ombar(nu))^2;
           Pt1(mu,j) = Pt1(mu,j) - pi*bt^4*abs(Rm1(j,mu,nu))^2*Ombar(nu)/(1+M*Ombar(nu))^2;
        end
        PHi0(mu,j) = 0;
        PHo0(mu,j) = 0;
        Pw0(mu,j)  = 0;
        if h ==0
           PHi1(mu,j) = pi/2*M*abs(Cmmu(mu))^2*((m^2+StM^2)*(iM/di)^2-StM^2);
        else
           PHi1(mu,j) = pi/2*M*abs(Cmmu(mu))^2*((m^2+StM^2)*((dkh*iM-dih*kM)/dki)^2-StM^2-(m^2+Sth^2)/(Sth*dki)^2);
        end
        PHo1(mu,j) = pi/2*M*abs(Cmmu(mu))^2*(StM^2-(m^2+StM^2)*(kM/dk)^2);
        Pw1(mu,j)  = pi/2*bt^2*(1-M)^2*abs(Um(mu)*(1+Ombar(mu))/(1-M*Ombar(mu))*Ktldmin(mu)/KmM)^2;
        Pw1(mu,j)  = Pw1(mu,j)*real(KtlM*(M/(1-M*Ombar(mu)) - M/(1-M) - dKpM));
       else
        if ~ismember(scenar,[4,5,6])
          if mu>=mco
            Ommu = abs(Ombar(mu));
            Pt(mu,j) = -2*pi*bt^4*Ommu/(1+M^2*Ommu^2)^2;
            Pt(mu,j) = Pt(mu,j)*(imag(Rm(j,mu,mu))*(1-M^2*Ommu^2)+real(Rm(j,mu,mu))*2*M*Ommu);
          else
            Pt(mu,j) = pi*bt^4*Ombar(mu)/(1-M*Ombar(mu))^2;
          end
          for nu=1:mco-1
             Pt(mu,j) = Pt(mu,j) - pi*bt^4*abs(Rm(j,mu,nu))^2*Ombar(nu)/(1+M*Ombar(nu))^2;
          end
          if gm~=0 && M~=0
             if h ==0
                   PHi(mu,j) = pi/2*M*abs(gm*Cmmu(mu))^2*((m^2+StM^2)*(iM/di)^2-StM^2);
             else
                PHi(mu,j) = pi/2*M*abs(gm*Cmmu(mu))^2*((m^2+StM^2)*((dkh*iM-dih*kM)/dki)^2-StM^2-(m^2+Sth^2)/(Sth*dki)^2);
             end
             PHo(mu,j) = pi/2*M*abs(gm*Cmmu(mu))^2*(StM^2-(m^2+StM^2)*(kM/dk)^2);
             Pw(mu,j)  = pi/2*bt^2*(1-M)^2*abs(Um(mu)*(1+Ombar(mu))/(1-M*Ombar(mu))*Ktldmin(mu)/KmM)^2;
             Pw(mu,j)  = Pw(mu,j)*real(KtlM*(conj(gm)*M/(1-M*Ombar(mu)) - abs(gm)^2*M/(1-M) - abs(gm)^2*dKpM));
          else
             PHi(mu,j) = 0;
             PHo(mu,j) = 0;
             Pw(mu,j)  = 0;
          end
        else
          if mu>=mco
            Ommu   = abs(Ombar(mu));
            Ptaux  = 2*pi*bt^4*Ommu/(1+M^2*Ommu^2)^2;
            Ptaux  = Ptaux*(imag(Am(j,mu)*conj(Bm(j,mu)))*(1-M^2*Ommu^2)-real(Am(j,mu)*conj(Bm(j,mu))*2*M*Ommu));
            PtA(j) = PtA(j) + Ptaux;
          else
            Ptaux  = pi*bt^4*Ombar(mu);
            Ptaux  = Ptaux*((abs(Am(j,mu))/(1-M*Ombar(mu)))^2 - (abs(Bm(j,mu))/(1+M*Ombar(mu)))^2);
            PtA(j) = PtA(j) + Ptaux;
          end
          if gm~=0 && M~=0
             PHiA(j) = PHiA(j) + Am(j,mu)*Cmmu(mu);
             PwA(j) = PwA(j) + Am(j,mu)*Um(mu)*(1+Ombar(mu))/(1-M*Ombar(mu))*Ktldmin(mu)/KmM*(M/(1-M*Ombar(mu)) - gm*M/(1-M) -gm*dKpM);
             PwB = PwB + conj(Am(j,mu)*Um(mu)*(1+Ombar(mu))/(1-M*Ombar(mu))*Ktldmin(mu)/KmM);
             if mu == mumx
               PHiA(j) = pi/2*M*abs(gm)^2*abs(PHiA(j))^2;
               PHoA(j) = PHiA(j)*(StM^2-(m^2+StM^2)*(kM/dk)^2);
               if h ==0
                   PHiA(j) = PHiA(j)*((m^2+StM^2)*(iM/di)^2-StM^2);
               else
                   PHiA(j) = PHiA(j)*((m^2+StM^2)*((dkh*iM-dih*kM)/dki)^2-StM^2-(m^2+Sth^2)/(Sth*dki)^2);
               end
               PwA(j) = pi/2*bt^2*(1-M)^2*real(PwA(j)*PwB*conj(gm)*KtlM);
             end
          else
             PHiA(j) = 0;
             PHoA(j) = 0;
             PwA(j)  = 0;
          end
        end
       end % if KuttaCmp
       if powD  % directivity function
         if KuttaCmp & M~=0
          if m==0 & mu==1
             Dir0(1) = sqrt(2)/2*i*kappa^2*sqrt(1-h^2);  % xi=0. Note that Stationaire Phase assumes xi<>0, xi<>pi
             Dir1(1) = Dir0(1);
          else
             Dir0(1) = 0     ;
             Dir1(1) = 0     ;
          end
          if m==0                                             % xi=pi. Note that Stationaire Phase assumes xi<>0, xi<>pi
             Diraux       = i*kappa^2/2*Um(mu)*Ktldmin(mu)*(1+Ombar(mu))/(1-M*Ombar(mu))*Ksplit(1,-1,kappa,m,h,1,d,q);
             Dir0(NNx(j)) = Diraux*((1+M)/(Ombar(mu)+1));
             Dir1(NNx(j)) = Diraux*((1+M)/(Ombar(mu)+1)-M);
          else
             Dir0(NNx(j)) = 0     ;
             Dir1(NNx(j)) = 0     ;
          end
          for nxi = 2: NNx(j)-1
              if mod(nxi,Nx+1)==1
                if floor(nxi/(Nx+1))==mu-[m==0]
                    Diraux0 = 2*kappa^2*Ombar(mu)/(pi*alph(mu)*Um(mu)*dH2(m,alph(mu)));
                    Diraux1 = Diraux0;
                else
                    Diraux0 = 0;
                    Diraux1 = Diraux0;
                end
              else
                snxi = sin(xii(nxi));
                csxi = cos(xii(nxi));
                Diraux0 = Ksplit(1,csxi,kappa,m,h,1,d,q);
                Diraux1 = Diraux0                         ;
                Diraux0 = Diraux0.*(1-csxi)./snxi./dH2(m,kappa*snxi);
                Diraux1 = Diraux0                                      ;
                Diraux0 = Diraux0.*((1-M*csxi)./(Ombar(mu)-csxi));
                Diraux1 = Diraux1.*((1-M*csxi)./(Ombar(mu)-csxi)-M);
                Diraux0 = Diraux0.*i^m*kappa/(2*pi)*Um(mu)*Ktldmin(mu)*(1+Ombar(mu))/(1-M*Ombar(mu));
                Diraux1 = Diraux1.*i^m*kappa/(2*pi)*Um(mu)*Ktldmin(mu)*(1+Ombar(mu))/(1-M*Ombar(mu));
              end
              Dir0(nxi) = Diraux0; % this is D_{m,mu}, gm=0
              Dir1(nxi) = Diraux1; % this is D_{m,mu}, gm=1
          end
          % power Pf
          Pf0(mu,j) = 0;
          Pf1(mu,j) = 0;
          for nc=1:mco-[m==0]
            for nxi = 1:Nx
              Dm0(nc,nxi)  = Dir0((nc-1)*(Nx+1)+1+nxi);      % this is D_{m,mu} with angles xi
              Dm1(nc,nxi)  = Dir1((nc-1)*(Nx+1)+1+nxi);      % this is D_{m,mu} with angles xi
              PwIntg0(nc, nxi) = sin(xic(nc,nxi))./((1-M*cos(xic(nc,nxi))).^2).*abs(Dm0(nc,nxi)).^2;
              PwIntg1(nc, nxi) = sin(xic(nc,nxi))./((1-M*cos(xic(nc,nxi))).^2).*abs(Dm1(nc,nxi)).^2;
            end
            Pf0(mu,j) = Pf0(mu,j)+pi*bt^4/kappa^2*sum(reshape(wic(nc,1:Nx),1,Nx).*reshape(PwIntg0(nc,1:Nx),1,Nx));
            Pf1(mu,j) = Pf1(mu,j)+pi*bt^4/kappa^2*sum(reshape(wic(nc,1:Nx),1,Nx).*reshape(PwIntg1(nc,1:Nx),1,Nx));
          end
          Px0(mu,j) = Pf0(mu,j) + PHi0(mu,j) + PHo0(mu,j) + Pw0(mu,j);  % = Pf0
          Px1(mu,j) = Pf1(mu,j) + PHi1(mu,j) + PHo1(mu,j) + Pw1(mu,j);
         else  % KuttaCmp
          if m==0 & mu==1
             Dir(1) = sqrt(2)/2*i*kappa^2*sqrt(1-h^2);  % xi=0. Note that Stationaire Phase assumes xi<>0, xi<>pi
          else
             Dir(1) = 0 ;
          end
          if m==0                                           % xi=pi. Note that Stationaire Phase assumes xi<>0, xi<>pi
             Dir(NNx(j)) = i*kappa^2/2*Um(mu)*Ktldmin(mu)*(1+Ombar(mu))/(1-M*Ombar(mu))*Ksplit(1,-1,kappa,m,h,1,d,q)*((1+M)/(Ombar(mu)+1)-M*gm);
          else
             Dir(NNx(j)) = 0 ;
          end
          for nxi = 2: NNx(j)-1
              if mod(nxi,Nx+1)==1
                if floor(nxi/(Nx+1))==mu-[m==0]
                    Diraux = 2*i^m*kappa^2*Ombar(mu)/(pi*alph(mu)*Um(mu)*dH2(m,alph(mu)));
                else
                    Diraux = 0;
                end
              else
                snxi = sin(xii(nxi));
                csxi = cos(xii(nxi));
                Diraux = Ksplit(1,csxi,kappa,m,h,1,d,q);
                Diraux = Diraux*(1-csxi)/snxi/dH2(m,kappa*snxi);
                Diraux = Diraux*((1-M*csxi)/(Ombar(mu)-csxi)-M*gm);
                Diraux = Diraux*i^m*kappa/(2*pi)*Um(mu)*Ktldmin(mu)*(1+Ombar(mu))/(1-M*Ombar(mu));
              end
              Dir(nxi) = Diraux; % this is D_{m,mu}
          end
          % Pf
          if ~ismember(scenar,[4,5,6])
            Pf(mu,j) = 0;
               %%%%%%
               %if mu==mudir
               %  figure(20);
               %  clf
               %  hold on
               %  grid on
               %end
               %%%%%%
            for nc=1:mco-[m==0]
              for nxi = 1:Nx
                Dm(nc,nxi)  = Dir((nc-1)*(Nx+1)+1+nxi);      % this is D_{m,mu} with angles xi
                PwIntg(nc, nxi) = sin(xic(nc,nxi))./((1-M*cos(xic(nc,nxi))).^2).*abs(Dm(nc,nxi)).^2;
              end
               %%%%%%
               %if mu == mudir
               %  plot(xic(nc,2:Nx)*180/pi,PwIntg(nc,2:Nx)./(xic(nc,2:Nx).^(2*m+1))./((pi-xic(nc,2:Nx)).^(2*m+1)),'-k');
               %end
               %%%%%%
              Pf(mu,j) = Pf(mu,j) +pi*bt^4/kappa^2*sum(reshape(wic(nc,1:Nx),1,Nx).*reshape(PwIntg(nc,1:Nx),1,Nx));
            end
            Px(mu,j) = Pf(mu,j) + PHi(mu,j) + PHo(mu,j) + Pw(mu,j);
          else
            if mu==1
               DirA = zeros(1,NNx(j));
            end
            DirA = DirA + Am(j,mu)*Dir;
            if mu == mumx
               for nc=1:mco-[m==0]
                 for nxi = 1:Nx
                   DmA(nc,nxi)  = DirA((nc-1)*(Nx+1)+1+nxi);      % this is D_{m,mu} with angles xi
                   PwAIntg(nc, nxi) = sin(xic(nc,nxi))./((1-M*cos(xic(nc,nxi))).^2).*abs(DmA(nc,nxi)).^2;
                 end
                 PfA(j) = PfA(j) + pi*bt^4/kappa^2*sum(reshape(wic(nc,1:Nx),1,Nx).*reshape(PwAIntg(nc,1:Nx),1,Nx));
               end
            end
            PxA(j) = PfA(j) + PHiA(j) + PHoA(j) + PwA(j);
          end
         end  % end KuttaCmp
         if plotD  % prepare for plot
            if KuttaCmp & M~=0
               for nxi =1:NNx(j)
                   Dr0(mu,j,nxi) = E0(nxi)*Dir0(nxi);
                   Dr1(mu,j,nxi) = E0(nxi)*Dir1(nxi);
                   Em0(mu,j,nxi) = 20*log10(max(1e-16,abs(E0(nxi)*Dir0(nxi))));  % this is 20*log10(|E_{m,mu}|) with angles xib
                   Em1(mu,j,nxi) = 20*log10(max(1e-16,abs(E0(nxi)*Dir1(nxi))));
               end
            else
               % prepare for plot
               if ~ismember(scenar,[4,5,6])
                 for nxi =1:NNx(j)
                     Dr(mu,j,nxi) = E0(nxi)*Dir(nxi);
                     Em(mu,j,nxi) = 20*log10(max(1e-16,abs(E0(nxi)*Dir(nxi))));  % this is 20*log10(|E_{m,mu}|) with angles xib
                 end
               else
                 if mu==mumx
                   for nxi =1:NNx(j)
                       DrA(j,nxi) = E0(nxi)*DirA(nxi);
                       EmA(nxi) = 20*log10(max(1e-16,abs(E0(nxi)*DirA(nxi))));  % this is 20*log10(|E_{m,mu}|) with angles xib
                   end
                 end
               end
            end
         end % if plotD
       end  % if powD
     end  % mu
   end   % end j
   %
   if savedat
     if ~ismember(scenar,[4,5,6])
      % store Rm in dat-file
      fid = fopen(Rdatfile,'w');
      for j = 1:Nom
        fprintf(fid,'omega = %2.8f\n',omar(j));
        fprintf(fid,'nu/mu');
        for mu=1:mumx
            if mu==1
               fprintf(fid,blanks(14));
            else
               fprintf(fid,blanks(32));
            end
            fprintf(fid,'%2i',mu);
        end
        fprintf(fid,'\n');
        if KuttaCmp && M~=0
          for nu=1:mumx
            fprintf(fid,'%2i,  ',nu);
            for mu=1:mumx
                fprintf(fid,'%+1.12f%+1.12fi',real(Rm1(j,mu,nu)),imag(Rm1(j,mu,nu)));
                if mu<mumx
                   fprintf(fid,',  ');
                end
            end
            fprintf(fid,'\n');
          end
          fprintf(fid,'\n');
          for nu=1:mumx
            fprintf(fid,'%2i,  ',nu);
            for mu=1:mumx
                fprintf(fid,'%+1.12f%+1.12fi',real(Rm0(j,mu,nu)),imag(Rm0(j,mu,nu)));
                if mu<mumx
                   fprintf(fid,',  ');
                end
            end
            fprintf(fid,'\n');
          end;
        else  %KuttaCmp
          for nu=1:mumx
            fprintf(fid,'%2i,  ',nu);
            for mu=1:mumx
                fprintf(fid,'%+1.12f%+1.12fi',real(Rm(j,mu,nu)),imag(Rm(j,mu,nu)));
                if mu<mumx
                   fprintf(fid,',  ');
                end
            end
            fprintf(fid,'\n');
          end
        end
        if j<Nom
          fprintf(fid,'\n');
        end
      end % j
      fclose(fid);
     end % ~=4
     if plotD
       % store Dir in dat-file
       fid = fopen(Ddatfile,'w');
       for j=1:Nom
        if ~ismember(scenar,[4,5,6])
         for mu = 1:mumx
          for nxi=1:NNx(j)
           if KuttaCmp & M~=0
              fprintf(fid,'%3i,  %8.4f, (%+17.12f %+17.12fi ),  %+17.12f; (%+17.12f %+17.12fi ),  %+17.12f\n',nxi,xib(j,nxi)*180/pi,real(Dr1(mu,j,nxi)),imag(Dr1(mu,j,nxi)), abs(Dr1(mu,j,nxi)), real(Dr0(mu,j,nxi)),imag(Dr0(mu,j,nxi)),abs(Dr0(mu,j,nxi)));
           else
              fprintf(fid,'%3i,  %8.4f, (%+17.12f %+17.12fi ),  %+17.12f\n',nxi,xib(j,nxi)*180/pi,real(Dr(mu,j,nxi)),imag(Dr(mu,j,nxi)),abs(Dr(mu,j,nxi)));
           end
          end
          if j<Nom || mu<mumx
            fprintf(fid,'\n');
          end
         end
        else
         for nxi=1:NNx(j)
           fprintf(fid,'%3i,  %8.4f, (%+17.12f %+17.12fi ),  %+17.12f\n',nxi,xib(j,nxi)*180/pi,real(DrA(j,nxi)),imag(DrA(j,nxi)),abs(DrA(j,nxi)));
         end
         if j<Nom
           fprintf(fid,'\n');
         end
        end
       end
       fclose(fid);
     end
   end % end savedat
   %
   if plotR             % plot Rm(omega)
     if KuttaCmp && M~=0
         movegui(figure(1),'north'), clf
         for nR = 1:NR
           absR0 = abs(Rm0(:,muR(nR),nuR(nR)));
           absR1 = abs(Rm1(:,muR(nR),nuR(nR)));
           ps = dpsimplify([omar.',absR1], dperr);    plot(ps(:,1),ps(:,2),'-', 'Color',cell2mat(klrR(1+mod(nR-1,length(klrR))))); hold on
           ps = dpsimplify([omar.',absR0], dperr);    plot(ps(:,1),ps(:,2),'--','Color',cell2mat(klrR(1+mod(nR-1,length(klrR)))),'Linewidth',0.2);
           xlabel('\omega','Rotation',0);
           ylabel('|R|','Rotation',0);
           grid on
         end
         movegui(figure(2),'south'), clf
         for nR = 1:NR
           angR1 = unwrap(angle(Rm1(:,muR(nR),nuR(nR))))+shph(nR);
           angR0 = unwrap(angle(Rm0(:,muR(nR),nuR(nR))))+shph0(nR);
           ps = dpsimplify([omar.',angR1/pi], dperr); plot(ps(:,1),ps(:,2),'-' ,'Color',cell2mat(klrR(1+mod(nR-1,length(klrR))))); hold on
           ps = dpsimplify([omar.',angR0/pi], dperr); plot(ps(:,1),ps(:,2),'--','Color',cell2mat(klrR(1+mod(nR-1,length(klrR)))),'Linewidth',0.2);
           xlabel('\omega','Rotation',0);
           ylabel('\phi/\pi','Rotation',0);
           grid on
         end
     else
         movegui(figure(1),'north'), clf
         for nR = 1:NR
           absR = abs(Rm(:,muR(nR),nuR(nR)));
           ps = dpsimplify([omar.',absR], dperr);     plot(ps(:,1),ps(:,2),'-', 'Color',cell2mat(klrR(1+mod(nR-1,length(klrR)))));
           xlabel('\omega','Rotation',0);
           ylabel('|R|','Rotation',0);
           grid on
           hold on
         end
         movegui(figure(2),'south'), clf
         for nR = 1:NR
           angR = unwrap(angle(Rm(:,muR(nR),nuR(nR))))+shph(nR);
           ps = dpsimplify([omar.',angR/pi], dperr);  plot(ps(:,1),ps(:,2),'-', 'Color',cell2mat(klrR(1+mod(nR-1,length(klrR)))));
           xlabel('\omega','Rotation',0);
           ylabel('\phi/\pi','Rotation',0);
           grid on
           hold on
         end
     end
     figure(1)
     pax = gca;
     pax.FontSize = pltfnt;
     set(gcf,'color',[1 1 1])
     set(gcf,'inverthardcopy','off')
     set(gcf,'papertype','a4letter')
     set(gcf,'paperunits','centimeters')
     set(gcf,'paperposition',pprsze)
     set(gcf,'clipping','off')
     eval(['print -depsc -cmyk ' epsR]);
     fixepsbbox(epsR,dx1,dy1,dx2,dy2)
     figure(2)
     pax = gca;
     pax.FontSize = pltfnt;
     set(gcf,'color',[1 1 1])
     set(gcf,'inverthardcopy','off')
     set(gcf,'papertype','a4letter')
     set(gcf,'paperunits','centimeters')
     set(gcf,'paperposition',pprsze)
     set(gcf,'clipping','off')
     eval(['print -depsc -cmyk ' epsphi]);
     fixepsbbox(epsphi,dx1,dy1,dx2,dy2)
   end % end plotR
   %
   if plotD     % plot directivity. if powD then powerD instead of directivity
     if KuttaCmp & M~=0
        for j=1:Nom
         if mudir==0
           Emax = ceil(max(max(max([Em0(1:mco-1,j,1:NNx(j)),Em1(1:mco-1,j,1:NNx(j))])))/5)*5;
         else
           Emax = ceil(max(max([Em0(mudir,j,1:NNx(j)),Em1(mudir,j,1:NNx(j))]))/5)*5;
         end
         Emin = Emax - Edif;
         movegui(figure(3+j),'north'), clf
         if mudir == 0
           for mu=1:mco-1  % we plot only cut-on modes
             Em0(mu,j,find(Em0(mu,j,1:NNx(j))<=Emin)) = Emin;
             Em1(mu,j,find(Em1(mu,j,1:NNx(j))<=Emin)) = Emin;
             polarplot( reshape(xib(j,1:NNx(j)),[1,NNx(j)]), reshape(Em1(mu,j,1:NNx(j)),[1,NNx(j)]),'-', 'Color',cell2mat(klrD(1+mod(mu-1,length(klrD)))));  hold on  %Note: the field itself is not always symmetric (m even/odd) but abs(Dir) is
             polarplot(-reshape(xib(j,1:NNx(j)),[1,NNx(j)]), reshape(Em1(mu,j,1:NNx(j)),[1,NNx(j)]),'-', 'Color',cell2mat(klrD(1+mod(mu-1,length(klrD)))));
             polarplot( reshape(xib(j,1:NNx(j)),[1,NNx(j)]), reshape(Em0(mu,j,1:NNx(j)),[1,NNx(j)]),'--','Color',cell2mat(klrD(1+mod(mu-1,length(klrD)))));
             polarplot(-reshape(xib(j,1:NNx(j)),[1,NNx(j)]), reshape(Em0(mu,j,1:NNx(j)),[1,NNx(j)]),'--','Color',cell2mat(klrD(1+mod(mu-1,length(klrD)))));
           end
         else
           Em0(mudir,j,find(Em0(mudir,j,1:NNx(j))<=Emin)) = Emin;
           Em1(mudir,j,find(Em1(mudir,j,1:NNx(j))<=Emin)) = Emin;
           polarplot( reshape(xib(j,1:NNx(j)),[1,NNx(j)]), reshape(Em1(mudir,j,1:NNx(j)),[1,NNx(j)]),'-', 'Color',cell2mat(klrD(1+mod(j-1,length(klrD)))));  hold on  %Note: the field itself is not always symmetric (m even/odd) but abs(Dir) is
           polarplot(-reshape(xib(j,1:NNx(j)),[1,NNx(j)]), reshape(Em1(mudir,j,1:NNx(j)),[1,NNx(j)]),'-', 'Color',cell2mat(klrD(1+mod(j-1,length(klrD)))));
           polarplot( reshape(xib(j,1:NNx(j)),[1,NNx(j)]), reshape(Em0(mudir,j,1:NNx(j)),[1,NNx(j)]),'--','Color',cell2mat(klrD(1+mod(j-1,length(klrD)))));
           polarplot(-reshape(xib(j,1:NNx(j)),[1,NNx(j)]), reshape(Em0(mudir,j,1:NNx(j)),[1,NNx(j)]),'--','Color',cell2mat(klrD(1+mod(j-1,length(klrD)))));
         end
         hold off
         pax = gca      ;
         pax.ThetaLim = [-180 180];
         pax.RLim = [Emin, Emax];
         pax.FontSize = pltfnt;
         set(gcf,'color',[1 1 1])
         set(gcf,'inverthardcopy','off')
         set(gcf,'papertype','a4letter')
         set(gcf,'paperunits','centimeters')
         set(gcf,'paperposition',pprsze)
         set(gcf,'clipping','off')
         eval(['print -depsc -cmyk ' convertStringsToChars(epsEmmu(j))]);
         fixepsbbox(epsEmmu(j),dx1,dy1,dx2,dy2)
        end
     else
        for j=1:Nom
         if ~ismember(scenar,[4,5,6])
           if mudir==0
             Emax = ceil(max(max(Em(1:mco-1,j,1:NNx(j))))/5)*5;
           else
             Emax = ceil(max(Em(mudir,j,1:NNx(j)))/5)*5;
           end
         else
           Emax = ceil(max(EmA(1:NNx(j)))/5)*5;
         end
         Emin = Emax - Edif;
         movegui(figure(3+j),'north'), clf
         if ~ismember(scenar,[4,5,6])
           if mudir ==0
             for mu=1:mco-1 % we plot only cut-on modes
               Em(mu,j,find(Em(mu,j,1:NNx(j))<=Emin)) = Emin;
               polarplot( reshape(xib(j,1:NNx(j)),[1,NNx(j)]), reshape(Em(mu,j,1:NNx(j)),[1,NNx(j)]),'-','Color',cell2mat(klrD(1+mod(mu-1,length(klrD)))));  hold on  %Note: the field itself is not always symmetric (m even/odd) but abs(Dir) is
               polarplot(-reshape(xib(j,1:NNx(j)),[1,NNx(j)]), reshape(Em(mu,j,1:NNx(j)),[1,NNx(j)]),'-','Color',cell2mat(klrD(1+mod(mu-1,length(klrD)))));
             end
           else
             Em(mudir,j,find(Em(mudir,j,1:NNx(j))<=Emin)) = Emin;
             polarplot( reshape(xib(j,1:NNx(j)),[1,NNx(j)]), reshape(Em(mudir,j,1:NNx(j)),[1,NNx(j)]),'-','Color',cell2mat(klrD(1+mod(mudir-1,length(klrD)))));  hold on  %Note: the field itself is not always symmetric (m even/odd) but abs(Dir) is
             polarplot(-reshape(xib(j,1:NNx(j)),[1,NNx(j)]), reshape(Em(mudir,j,1:NNx(j)),[1,NNx(j)]),'-','Color',cell2mat(klrD(1+mod(mudir-1,length(klrD)))));
           end
         else
            EmA(find(EmA(1:NNx(j))<=Emin)) = Emin;
            polarplot( reshape(xib(j,1:NNx(j)),[1,NNx(j)]), reshape(EmA(1:NNx(j)),[1,NNx(j)]),'-','Color',cell2mat(klrD(1+mod(j-1,length(klrD)))));  hold on  %Note: the field itself is not always symmetric (m even/odd) but abs(Dir) is
            polarplot(-reshape(xib(j,1:NNx(j)),[1,NNx(j)]), reshape(EmA(1:NNx(j)),[1,NNx(j)]),'-','Color',cell2mat(klrD(1+mod(j-1,length(klrD)))));
         end
         hold off
         pax = gca      ;
         pax.ThetaLim = [-180 180];
         pax.RLim = [Emin, Emax];
         pax.FontSize = pltfnt;
         set(gcf,'color',[1 1 1])
         set(gcf,'inverthardcopy','off')
         set(gcf,'papertype','a4letter')
         set(gcf,'paperunits','centimeters')
         set(gcf,'paperposition',pprsze)
         set(gcf,'clipping','off')
         eval(['print -depsc -cmyk ' convertStringsToChars(epsEmmu(j))]);
         fixepsbbox(epsEmmu(j),dx1,dy1,dx2,dy2)
        end
     end
   end % end plotD
   % end direct
   if plotP           % plotP = plot power. powD = both for Pt, Pf, PHi, PHo, and Pw. Otherwise: only Pt, PHi, PHo, Pw (xi-integration takes much time)
        if KuttaCmp & M~=0
           if powD
             for mu=1:mumx
               movegui(figure(2+mu),'east'), clf
               ps = dpsimplify([omar.',Pt1(mu,:).'],  dperr); plot(ps(:,1),ps(:,2),'-' ,'Color',         klPt  ); hold on
               ps = dpsimplify([omar.',Pt0(mu,:).'],  dperr); plot(ps(:,1),ps(:,2),'--','Color',         klPt   ,'Linewidth',0.2);
               ps = dpsimplify([omar.',PHi1(mu,:).'], dperr); plot(ps(:,1),ps(:,2),'-' ,'Color',         klPHi );
               ps = dpsimplify([omar.',PHo1(mu,:).'], dperr); plot(ps(:,1),ps(:,2),'-' ,'Color',         klPHo );
               ps = dpsimplify([omar.',Pw1(mu,:).'],  dperr); plot(ps(:,1),ps(:,2),'-' ,'Color',         klPw  );
               ps = dpsimplify([omar.',Pf1(mu,:).'],  dperr); plot(ps(:,1),ps(:,2),'-' ,'Color',         klPf  );
               ps = dpsimplify([omar.',Pf0(mu,:).'],  dperr); plot(ps(:,1),ps(:,2),'--','Color',         klPf   ,'Linewidth',0.2);
               if pwlbl
                title(sprintf('Pt, Pf, PHi, PHo, Pw, m = %1i, mu = %1i', m, mu),'Fontweight','normal')
                xlabel('\omega','Rotation',0);
                ylabel('Power','Rotation',90);
               end
               grid on
               pax = gca;
               pax.FontSize = pltfnt;
               ymax = max([max(Pt(mu,:)),max(Pf(mu,:)),max(PHi(mu,:)),max(PHi(mu,:)),max(Pw(mu,:))]);
               ymin = min([min(Pt(mu,:)),min(Pf(mu,:)),min(PHi(mu,:)),min(PHi(mu,:)),min(Pw(mu,:))]);
               dify = (ymax-ymin)*0.05;
               pax.YLim = [ymin-dify, ymax+dify];
               set(gcf,'color',[1 1 1])
               set(gcf,'inverthardcopy','off')
               set(gcf,'papertype','a4letter')
               set(gcf,'paperunits','centimeters')
               set(gcf,'paperposition',pprsze)
               set(gcf,'clipping','off')
               eval(['print -depsc -cmyk ' convertStringsToChars(epsP(mu))]);
               fixepsbbox(convertStringsToChars(epsP(mu)),dx1,dy1,dx2,dy2)
             end
           else
             for mu=1:mumx
               movegui(figure(2+mu),'east'), clf
               ps = dpsimplify([omar.',Pt1(mu,:).'],  dperr); plot(ps(:,1),ps(:,2),'-' ,'Color',         klPt  ); hold on
               ps = dpsimplify([omar.',Pt0(mu,:).'],  dperr); plot(ps(:,1),ps(:,2),'--','Color',         klPt   ,'Linewidth',0.2);
               ps = dpsimplify([omar.',PHi1(mu,:).'], dperr); plot(ps(:,1),ps(:,2),'-' ,'Color',         klPHi );
               ps = dpsimplify([omar.',PHo1(mu,:).'], dperr); plot(ps(:,1),ps(:,2),'-' ,'Color',         klPHo );
               ps = dpsimplify([omar.',Pw1(mu,:).'],  dperr); plot(ps(:,1),ps(:,2),'-' ,'Color',         klPw  );
               Pf1(mu,:) = Pt1(mu,:) - PHi1(mu,:) - PHo1(mu,:) - Pw1(mu,:);                                        % Pf wordt niet echt berekend, dus
               Pf0(mu,:) = Pt0(mu,:)                                      ;                                        % niet te gebruiken als check.
               ps = dpsimplify([omar.',Pf1(mu,:).'],  dperr); plot(ps(:,1),ps(:,2),'-' ,'Color',         klPf  );
               ps = dpsimplify([omar.',Pf0(mu,:).'],  dperr); plot(ps(:,1),ps(:,2),'--','Color',         klPf   ,'Linewidth',0.2);
               if pwlbl
                title(sprintf('Pt, Pf, PHi, PHo, Pw, m = %1i, mu = %1i', m, mu),'Fontweight','normal')
                xlabel('\omega','Rotation',0);
                ylabel('Power','Rotation',90);
               end
               grid on
               pax = gca;
               pax.FontSize = pltfnt;
               ymax = max([max(Pt(mu,:)),max(Pf(mu,:)),max(PHi(mu,:)),max(PHi(mu,:)),max(Pw(mu,:))]);
               ymin = min([min(Pt(mu,:)),min(Pf(mu,:)),min(PHi(mu,:)),min(PHi(mu,:)),min(Pw(mu,:))]);
               dify = (ymax-ymin)*0.05;
               pax.YLim = [ymin-dify, ymax+dify];
               set(gcf,'color',[1 1 1])
               set(gcf,'inverthardcopy','off')
               set(gcf,'papertype','a4letter')
               set(gcf,'paperunits','centimeters')
               set(gcf,'paperposition',pprsze)
               set(gcf,'clipping','off')
               eval(['print -depsc -cmyk ' convertStringsToChars(epsP(mu))]);
               fixepsbbox(convertStringsToChars(epsP(mu)),dx1,dy1,dx2,dy2)
             end
           end
           hold off
        else
           if powD
            if scenar == 2
             for mu=1:mumx
               movegui(figure(2+mu),'east'), clf
               ps = dpsimplify([omar.',Pt(mu,:).'],  dperr); plot(ps(:,1),ps(:,2),'-', 'Color',         klPt ); hold on
               ps = dpsimplify([omar.',Pf(mu,:).'],  dperr); plot(ps(:,1),ps(:,2),'-', 'Color',         klPf );
               if M~=0
                 ps = dpsimplify([omar.',PHi(mu,:).'], dperr); plot(ps(:,1),ps(:,2),'-', 'Color',         klPHi );
                 ps = dpsimplify([omar.',PHo(mu,:).'], dperr); plot(ps(:,1),ps(:,2),'-', 'Color',         klPHo );
                 ps = dpsimplify([omar.',Pw(mu,:).'],  dperr); plot(ps(:,1),ps(:,2),'-', 'Color',         klPw  );
               end
               if pwlbl
                if M~=0
                  title(sprintf('Pt, Pf, PHi, PHo, Pw, m = %1i, mu = %1i', m, mu),'Fontweight','normal')
                else
                  title(sprintf('Pt, Pf, m = %1i, mu = %1i', m, mu),'Fontweight','normal')
                end
                xlabel('\omega','Rotation',0);
                ylabel('Power','Rotation',90);
               end
               grid on
               pax = gca;
               pax.FontSize = pltfnt;
               ymax = max([max(Pt(mu,:)),max(Pf(mu,:)),max(PHi(mu,:)),max(PHi(mu,:)),max(Pw(mu,:))]);
               ymin = min([min(Pt(mu,:)),min(Pf(mu,:)),min(PHi(mu,:)),min(PHi(mu,:)),min(Pw(mu,:))]);
               dify = (ymax-ymin)*0.05;
               pax.YLim = [ymin-dify, ymax+dify];
               set(gcf,'color',[1 1 1])
               set(gcf,'inverthardcopy','off')
               set(gcf,'papertype','a4letter')
               set(gcf,'paperunits','centimeters')
               set(gcf,'paperposition',pprsze)
               set(gcf,'clipping','off')
               eval(['print -depsc -cmyk ' convertStringsToChars(epsP(mu))]);
               fixepsbbox(convertStringsToChars(epsP(mu)),dx1,dy1,dx2,dy2)
             end
            else
             movegui(figure(2),'east'), clf
             ps = dpsimplify([omar.',PtA(:)],  dperr); plot(ps(:,1),ps(:,2),'-', 'Color',         klPt ); hold on
             ps = dpsimplify([omar.',PfA(:)],  dperr); plot(ps(:,1),ps(:,2),'-', 'Color',         klPf );
             if M~=0
               ps = dpsimplify([omar.',PHiA(:)], dperr); plot(ps(:,1),ps(:,2),'-', 'Color',         klPHi );
               ps = dpsimplify([omar.',PHoA(:)], dperr); plot(ps(:,1),ps(:,2),'-', 'Color',         klPHo );
               ps = dpsimplify([omar.',PwA(:)],  dperr); plot(ps(:,1),ps(:,2),'-', 'Color',         klPw  );
             end
             if pwlbl
              if M~=0
                title(sprintf('Pt, Pf, PHi, PHo, Pw, m = %1i, mu = %1i', m, mu),'Fontweight','normal')
              else
                title(sprintf('Pt, Pf, m = %1i, mu = %1i', m, mu),'Fontweight','normal')
              end
              xlabel('\omega','Rotation',0);
              ylabel('Power','Rotation',90);
             end
             grid on
             pax = gca;
             pax.FontSize = pltfnt;
             ymax = max([max(PtA(:)),max(PfA(:)),max(PHiA(:)),max(PHiA(:)),max(PwA(:))]);
             ymin = min([min(PtA(:)),min(PfA(:)),min(PHiA(:)),min(PHiA(:)),min(PwA(:))]);
             dify = (ymax-ymin)*0.05;
             pax.YLim = [ymin-dify, ymax+dify];
             set(gcf,'color',[1 1 1])
             set(gcf,'inverthardcopy','off')
             set(gcf,'papertype','a4letter')
             set(gcf,'paperunits','centimeters')
             set(gcf,'paperposition',pprsze)
             set(gcf,'clipping','off')
             eval(['print -depsc -cmyk ' convertStringsToChars(epsPA)]);
             fixepsbbox(convertStringsToChars(epsPA),dx1,dy1,dx2,dy2)
            end
           else
            if scenar == 3
             for mu=1:mumx
               movegui(figure(2+mu),'east'), clf
               ps = dpsimplify([omar.',Pt(mu,:).'],  dperr); plot(ps(:,1),ps(:,2),'-', 'Color',         klPt ); hold on
               Pf(mu,:) = Pt(mu,:) - PHi(mu,:) - PHo(mu,:) - Pw(mu,:);                                        % Pf wordt niet echt berekend, dus niet gebruiken als check
               ps = dpsimplify([omar.',Pf(mu,:).'],  dperr); plot(ps(:,1),ps(:,2),'-', 'Color',         klPf );
               if M~=0
                 ps = dpsimplify([omar.',PHi(mu,:).'], dperr); plot(ps(:,1),ps(:,2),'-', 'Color',         klPHi );
                 ps = dpsimplify([omar.',PHo(mu,:).'], dperr); plot(ps(:,1),ps(:,2),'-', 'Color',         klPHo );
                 ps = dpsimplify([omar.',Pw(mu,:).'],  dperr); plot(ps(:,1),ps(:,2),'-', 'Color',         klPw  );
               end
               if pwlbl
                if M~=0
                  title(sprintf('Pt, Pf, PHi, PHo, Pw, m = %1i, mu = %1i', m, mu),'Fontweight','normal')
                else
                  title(sprintf('Pt, Pf, m = %1i, mu = %1i', m, mu),'Fontweight','normal')
                end
                xlabel('\omega','Rotation',0);
                ylabel('Power','Rotation',90);
               end
               grid on
               pax = gca;
               pax.FontSize = pltfnt;
               ymax = max([max(Pt(mu,:)),max(Pf(mu,:)),max(PHi(mu,:)),max(PHi(mu,:)),max(Pw(mu,:))]);
               ymin = min([min(Pt(mu,:)),min(Pf(mu,:)),min(PHi(mu,:)),min(PHi(mu,:)),min(Pw(mu,:))]);
               dify = (ymax-ymin)*0.05;
               pax.YLim = [ymin-dify, ymax+dify];
               set(gcf,'color',[1 1 1])
               set(gcf,'inverthardcopy','off')
               set(gcf,'papertype','a4letter')
               set(gcf,'paperunits','centimeters')
               set(gcf,'paperposition',pprsze)
               set(gcf,'clipping','off')
               eval(['print -depsc -cmyk ' convertStringsToChars(epsP(mu))]);
               fixepsbbox(convertStringsToChars(epsP(mu)),dx1,dy1,dx2,dy2)
             end
            else
               movegui(figure(2),'east'), clf
               ps = dpsimplify([omar.',PtA(:)],  dperr); plot(ps(:,1),ps(:,2),'-', 'Color',         klPt ); hold on
               PfA(:) = PtA(:) - PHiA(:) - PHoA(:) - PwA(:);                                        % Pf wordt niet echt berekend, dus niet gebruiken als check
               ps = dpsimplify([omar.',PfA(:)],  dperr); plot(ps(:,1),ps(:,2),'-', 'Color',         klPf );
               if M~=0
                 ps = dpsimplify([omar.',PHiA(:)], dperr); plot(ps(:,1),ps(:,2),'-', 'Color',         klPHi );
                 ps = dpsimplify([omar.',PHoA(:)], dperr); plot(ps(:,1),ps(:,2),'-', 'Color',         klPHo );
                 ps = dpsimplify([omar.',PwA(:)],  dperr); plot(ps(:,1),ps(:,2),'-', 'Color',         klPw  );
               end
               if pwlbl
                if M~=0
                  title(sprintf('Pt, Pf, PHi, PHo, Pw, m = %1i, mu = %1i', m, mu),'Fontweight','normal')
                else
                  title(sprintf('Pt, Pf, m = %1i, mu = %1i', m, mu),'Fontweight','normal')
                end
                xlabel('\omega','Rotation',0);
                ylabel('Power','Rotation',90);
               end
               grid on
               pax = gca;
               pax.FontSize = pltfnt;
               ymax = max([max(PtA(:)),max(PfA(:)),max(PHiA(:)),max(PHiA(:)),max(PwA(:))]);
               ymin = min([min(PtA(:)),min(PfA(:)),min(PHiA(:)),min(PHiA(:)),min(PwA(:))]);
               dify = (ymax-ymin)*0.05;
               pax.YLim = [ymin-dify, ymax+dify];
               set(gcf,'color',[1 1 1])
               set(gcf,'inverthardcopy','off')
               set(gcf,'papertype','a4letter')
               set(gcf,'paperunits','centimeters')
               set(gcf,'paperposition',pprsze)
               set(gcf,'clipping','off')
               eval(['print -depsc -cmyk ' convertStringsToChars(epsPA)]);
               fixepsbbox(convertStringsToChars(epsPA),dx1,dy1,dx2,dy2)
            end
           end
           hold off
        end
   else           % scenar = 1 or 4
     fid = fopen(Pdatfile,'w');
     fprintf(fid,'co-fr =  % 1.15f',bt*alph(1)); fprintf(fid,', %1.15f',bt*alph(2:mumx)); fprintf(fid,'\n');
     for j=1:Nom
       fprintf(fid,'omega = % 12.15f,  muco = %1i\n', omar(j), muco(j));
       if KuttaCmp & M~=0
           fprintf(fid,'mu    =          %1i', 1); for mu=2:mumx, fprintf(fid,blanks(19));fprintf(fid,'%1i', mu); end; fprintf(fid,'\n');
          myprintf(fid,'Pt1   =  ', Pt1(1,j) ); for mu=2:mumx, myprintf(fid,'  ',Pt1(mu,j) ); end; fprintf(fid,'\n');
          myprintf(fid,'Px1   =  ', Px1(1,j) ); for mu=2:mumx, myprintf(fid,'  ',Px1(mu,j) ); end; fprintf(fid,'\n');
          myprintf(fid,'Pf1   =  ', Pf1(1,j) ); for mu=2:mumx, myprintf(fid,'  ',Pf1(mu,j) ); end; fprintf(fid,'\n');
          myprintf(fid,'PHi1  =  ', PHi1(1,j)); for mu=2:mumx, myprintf(fid,'  ',PHi1(mu,j)); end; fprintf(fid,'\n');
          myprintf(fid,'PHo1  =  ', PHo1(1,j)); for mu=2:mumx, myprintf(fid,'  ',PHo1(mu,j)); end; fprintf(fid,'\n');
          myprintf(fid,'Pw1   =  ', Pw1(1,j) ); for mu=2:mumx, myprintf(fid,'  ',Pw1(mu,j) ); end; fprintf(fid,'\n');
          myprintf(fid,'Pt0   =  ', Pt0(1,j) ); for mu=2:mumx, myprintf(fid,'  ',Pt0(mu,j) ); end; fprintf(fid,'\n');
          myprintf(fid,'Px0   =  ', Px0(1,j) ); for mu=2:mumx, myprintf(fid,'  ',Px0(mu,j) ); end; fprintf(fid,'\n');
          myprintf(fid,'Pf0   =  ', Pf0(1,j) ); for mu=2:mumx, myprintf(fid,'  ',Pf0(mu,j) ); end; fprintf(fid,'\n');
       else
          if scenar ~=4
            fprintf(fid,'mu    =          %1i', 1); for mu=2:mumx, fprintf(fid,blanks(19));fprintf(fid,'%1i', mu); end; fprintf(fid,'\n');
           myprintf(fid,'Pt    =  ', Pt(1,j) ); for mu=2:mumx, myprintf(fid,'  ', Pt(mu,j) ); end; fprintf(fid,'\n');
           myprintf(fid,'Px    =  ', Px(1,j) ); for mu=2:mumx, myprintf(fid,'  ', Px(mu,j) ); end; fprintf(fid,'\n');
           myprintf(fid,'Pf    =  ', Pf(1,j) ); for mu=2:mumx, myprintf(fid,'  ', Pf(mu,j) ); end; fprintf(fid,'\n');
           myprintf(fid,'PHi   =  ', PHi(1,j)); for mu=2:mumx, myprintf(fid,'  ', PHi(mu,j)); end; fprintf(fid,'\n');
           myprintf(fid,'PHo   =  ', PHo(1,j)); for mu=2:mumx, myprintf(fid,'  ', PHo(mu,j)); end; fprintf(fid,'\n');
           myprintf(fid,'Pw    =  ', Pw(1,j) ); for mu=2:mumx, myprintf(fid,'  ', Pw(mu,j) ); end; fprintf(fid,'\n');
          else
           myprintf(fid,'PtA   =  ', PtA(j)  ); fprintf(fid,'\n');
           myprintf(fid,'PxA   =  ', PxA(j)  ); fprintf(fid,'\n');
           myprintf(fid,'PfA   =  ', PfA(j)  ); fprintf(fid,'\n');
           myprintf(fid,'PHiA  =  ', PHiA(j) ); fprintf(fid,'\n');
           myprintf(fid,'PHoA  =  ', PHoA(j) ); fprintf(fid,'\n');
           myprintf(fid,'PwA   =  ', PwA(j)  ); fprintf(fid,'\n');
          end
       end
       if j<Nom, fprintf(fid,'\n'); end
     end
     fclose(fid);
     if Nom>1, fprintf('\n'); end;
     fprintf('co-fr =  % 1.15f',bt*alph(1)); fprintf(', %1.15f',bt*alph(2:mumx)); fprintf('\n');
     for j=1:Nom
       fprintf('omega = % 12.15f,  muco = %1i\n', omar(j), muco(j));
       if KuttaCmp & M~=0
           fprintf('mu    =          %1i', 1);  for mu=2:mumx, fprintf(blanks(19));fprintf('%1i', mu); end; fprintf('\n');
          myprintf(1,'Pt1   =  ', Pt1(1,j) ); for mu=2:mumx, myprintf(1,'  ', Pt1(mu,j) ); end; fprintf('\n');
          myprintf(1,'Px1   =  ', Px1(1,j) ); for mu=2:mumx, myprintf(1,'  ', Px1(mu,j) ); end; fprintf('\n');
          myprintf(1,'Pf1   =  ', Pf1(1,j) ); for mu=2:mumx, myprintf(1,'  ', Pf1(mu,j) ); end; fprintf('\n');
          myprintf(1,'PHi1  =  ', PHi1(1,j)); for mu=2:mumx, myprintf(1,'  ', PHi1(mu,j)); end; fprintf('\n');
          myprintf(1,'PHo1  =  ', PHo1(1,j)); for mu=2:mumx, myprintf(1,'  ', PHo1(mu,j)); end; fprintf('\n');
          myprintf(1,'Pw1   =  ', Pw1(1,j) ); for mu=2:mumx, myprintf(1,'  ', Pw1(mu,j) ); end; fprintf('\n');
          myprintf(1,'Pt0   =  ', Pt0(1,j) ); for mu=2:mumx, myprintf(1,'  ', Pt0(mu,j) ); end; fprintf('\n');
          myprintf(1,'Px0   =  ', Px0(1,j) ); for mu=2:mumx, myprintf(1,'  ', Px0(mu,j) ); end; fprintf('\n');
          myprintf(1,'Pf0   =  ', Pf0(1,j) ); for mu=2:mumx, myprintf(1,'  ', Pf0(mu,j) ); end; fprintf('\n');
       else
        if ~ismember(scenar,[4,5,6])
           fprintf(1,'mu    =          %1i', 1); for mu=2:mumx, fprintf(blanks(19));fprintf('%1i', mu); end; fprintf('\n');
          myprintf(1,'Pt    =  ', Pt(1,j) ); for mu=2:mumx, myprintf(1,'  ', Pt(mu,j) ); end; fprintf('\n');
          myprintf(1,'Px    =  ', Px(1,j) ); for mu=2:mumx, myprintf(1,'  ', Px(mu,j) ); end; fprintf('\n');
          myprintf(1,'Pf    =  ', Pf(1,j) ); for mu=2:mumx, myprintf(1,'  ', Pf(mu,j) ); end; fprintf('\n');
          myprintf(1,'PHi   =  ', PHi(1,j)); for mu=2:mumx, myprintf(1,'  ', PHi(mu,j)); end; fprintf('\n');
          myprintf(1,'PHo   =  ', PHo(1,j)); for mu=2:mumx, myprintf(1,'  ', PHo(mu,j)); end; fprintf('\n');
          myprintf(1,'Pw    =  ', Pw(1,j) ); for mu=2:mumx, myprintf(1,'  ', Pw(mu,j) ); end; fprintf('\n');
        else
          myprintf(1,'PtA   =  ', PtA(j)  ); fprintf('\n');
          myprintf(1,'PxA   =  ', PxA(j)  ); fprintf('\n');
          myprintf(1,'PfA   =  ', PfA(j)  ); fprintf('\n');
          myprintf(1,'PHiA  =  ', PHiA(j) ); fprintf('\n');
          myprintf(1,'PHoA  =  ', PHoA(j) ); fprintf('\n');
          myprintf(1,'PwA   =  ', PwA(j)  ); fprintf('\n');
        end
       end
       if j<Nom, fprintf('\n'); end
     end
     fprintf('\n');
   end
   %
   % save the used input parameters
   %
   fid = fopen(inpfile,'w');
   tend = toc(tstart);
   fprintf(fid,'%%time     =  %2.4f sec\n',tend);
   fprintf(fid,'%%co-fr    =  %1.15f',bt*alph(1)); fprintf(fid,', %1.15f',bt*alph(2:mumx)); fprintf(fid,'\n');
   fprintf(fid,'test      =  %0i;\n',test);
   fprintf(fid,'pwlbl     =  %0i;\n',pwlbl);
   fprintf(fid,'scenar    =  %0i;\n',scenar);
   fprintf(fid,'KuttaCmp  =  %0i;\n',KuttaCmp);
   fprintf(fid,'h         =  %0.15g;\n',h);
   fprintf(fid,'M         =  %0.15g;\n',M);
   fprintf(fid,'gm        =  %0.15g;\n',gm);
   switch scenar
    case {1,4}
   fprintf(fid,['omar      =  [' sprintf(repmat('%4.15g, ', 1, length(omar)-1), omar(1:end-1)) sprintf('%4.15g];\n', omar(end))]);
    case {2,3,5,6}
   fprintf(fid,'omar      =  [%5.3f:%4.3f:%5.3f];\n',omar(1),(omar(end)-omar(1))/(length(omar)-1),omar(end));
   end
   fprintf(fid,'omfln     =  %0.15g;\n',omfn);
   fprintf(fid,'Nx        =  %0i;\n',Nx);
   fprintf(fid,'m         =  %0i;\n',m);
   fprintf(fid,'mumx      =  %0i;\n',mumx);
   fprintf(fid,'mudir     =  %0i;\n',mudir);
   fprintf(fid,['muR       =  [' sprintf(repmat('%1i, ', 1, length(muR)-1), muR(1:end-1)) sprintf('%1i];\n', muR(end))]);
   fprintf(fid,['nuR       =  [' sprintf(repmat('%1i, ', 1, length(nuR)-1), nuR(1:end-1)) sprintf('%1i];\n', nuR(end))]);
      klrRstring = '{';
      for cj=1:length(klrR)-1
          klrRstring = strcat(klrRstring, '[', strip(sprintf(' %.15g',cell2mat(klrR(cj)))),'], ');
      end
      klrRstring = strcat(klrRstring, '[', strip(sprintf(' %.15g',cell2mat(klrR(length(klrR))))),']};\n');
   fprintf(fid,['klrR      =  ' klrRstring]);
      klrDstring = '{';
      for cj=1:length(klrD)-1
          klrDstring = strcat(klrDstring, '[', strip(sprintf(' %.15g',cell2mat(klrD(cj)))),'], ');
      end
      klrDstring = strcat(klrDstring, '[', strip(sprintf(' %.15g',cell2mat(klrR(length(klrD))))),']};\n');
   fprintf(fid,['klrD      =  ' klrDstring]);
   fprintf(fid,['klPt      =  [' strip(sprintf(' %.15g',klPt))  '];\n']);
   fprintf(fid,['klPf      =  [' strip(sprintf(' %.15g',klPf))  '];\n']);
   fprintf(fid,['klPHi     =  [' strip(sprintf(' %.15g',klPHi)) '];\n']);
   fprintf(fid,['klPHo     =  [' strip(sprintf(' %.15g',klPHo)) '];\n']);
   fprintf(fid,['klPw      =  [' strip(sprintf(' %.15g',klPw))  '];\n']);
   fprintf(fid,['shph      =  [' sprintf(repmat('%1i, ', 1, length(shph)-1), round(shph(1:end-1)/2/pi)) sprintf('%1i]*2*pi;\n', round(shph(end)/2/pi))]);
   fprintf(fid,['shph0     =  [' sprintf(repmat('%1i, ', 1, length(shph0)-1), round(shph0(1:end-1)/2/pi)) sprintf('%1i]*2*pi;\n', round(shph0(end)/2/pi))]);
   fprintf(fid,'Edif      =  %0i;\n',Edif);
   fprintf(fid,'d         =  %4.2f;\n',d);
   fprintf(fid,'q         =  %0.2g;\n',q);
   fprintf(fid,'Nom       =  %1i;\n',Nom);
   if ismember(scenar,[4,5,6])
     if imag(Am)==0
       for mu = 1:mumx, fprintf(fid,[ sprintf('Am(1:%1i,%1i)',Nom,mu) sprintf(' =  [') sprintf(repmat('%+2.15g, ', 1, Nom-1), Am(1:Nom-1,mu)) sprintf('%+2.15g];\n', Am(Nom,mu))]); end
     else
       for mu = 1:mumx, fprintf(fid,  sprintf('Am(1:%1i,%1i) = [',Nom,mu));
           for j=1:Nom-1
            fprintf(fid,sprintf('%+2.15g%+2.15g*i, ', real(Am(j,mu)),imag(Am(j,mu))));
           end
           fprintf(fid,sprintf('%+2.15g%+2.15g*i];\n', real(Am(Nom,mu)), imag(Am(Nom,mu)))); end
     end
   end
   fclose(fid);
   %
   fprintf('time  = %2.4f sec\n',tend);
end  % end main program
%
%% auxiliary functions
function y = findalphas(n,hh,P)
    j0    = [0.01:0.1:n+P*ceil(pi/(1-hh)/10)*10];             % fine grid initial values. could be sharpened for h=0 but seems OK
    jm0   = j0(find(abs(diff(BesEig(j0)>=0))));               % find the sign changes >0 and add to initial values Alpha0
    alpha = fsolve(@(x) BesEig(x),jm0,optimoptions('fsolve','Display','off'));
    if n==0
        alpha = [0,alpha];                                    % add 0 (plane wave) if m=0
    end
    y       = alpha(1:P);
    function y = BesEig(x)
       if hh==0
        switch n
        case {1,2,3}
         y = dJ(n,x)              ;                           % beter voor m=1, 2, 3
        otherwise
         y = dJ(n,x).*sqrt(pi*x/2);                           % beter voor m=0, 4,
        end
       else
         switch n
         case {0,1,2}
           y = (dJ(n,x).*dY(n,x*hh) - dY(n,x).*dJ(n,x*hh))                 ;   % m = 0 1 2     5
         otherwise
           y = (dJ(n,x).*dY(n,x*hh) - dY(n,x).*dJ(n,x*hh)).*x*sqrt(hh)*2/pi;   % m = 0     3 4 5 ....
         end
       end
    end
end
function y = Ktilde(t,k,n,hh)            % note symmetry: Ktilde(-t,k,n,hh,btt) = Ktilde(t,k,n,hh,btt)
         y = LL(n,k*v(t),hh);            % k = kappa = omega/bt
end
function y = dKtilde(t,k,n,hh)           % note symmetry: Ktilde(-t,k,n,hh,btt) = Ktilde(t,k,n,hh,btt)
         y = -k*t.*dLL(n,k*v(t),hh)./v(t);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y = Ksplit(option,t,k,n,hh,t0,dd,qq)    % k = kappa, t0 is reference contour bump
         % option =  1: K_+
         % option = -1: K_-
         % option =  2: dK_+(t)/dt / K_+(t) logaritmische afgeleide van K_+
         % option = -2: dK_-(t)/dt / K_-(t) logaritmische afgeleide van K_-
         % let op: t=0 singulier. Let ook op de integratiecontour als Im(t)>0. Als t toch onder de contour, dan effectief in C^-
         % k = kappa = omega/bt wordt meegeven
    switch option
    case -1 % K-
         if imag(t)<0 | (imag(t)<1E-12 & real(t)>0)    % t \in \IC^-
           y = exp(Y1());
         else                                       % t \in \IC^+
           y = Ktilde(t,k,n,hh)./exp(-Y1());
         end
    case  1 % K+
         if imag(t)<0 | (imag(t)<1E-12 & real(t)>0)    % t \in \IC^-
           y = Ktilde(t,k,n,hh)./exp(Y1());
         else                                       % t \in \IC^+
           y = exp(-Y1());
         end
    case -2 % K-'/K-
         if imag(t)<0 | (imag(t)<1E-12 & real(t)>0)    % t \in \IC^-
            y = Y2();
         else                                       % t \in \IC^+
            y = dKtilde(t,k,n,hh)./Ktilde(t,k,n,hh) +Y2();
         end
    case  2 % K+'/K+
         if imag(t)<0 | (imag(t)<1E-12 & real(t)>0)    % t \in \IC^-
           y = dKtilde(t,k,n,hh)./Ktilde(t,k,n,hh) - Y2();
         else                                       % t \in \IC^+
           y = -Y2() ;
         end
    otherwise
         display(['Wrong Ksplit option'])
    end

    function y = Y1()    % integral in Ksplit
            y = t*(integral(@(s) Kintg(s),0,1,'AbsTol',1e-12,'RelTol',1e-8 ))*i/pi ; % exp(integral(@(s) Kintg(s),0,1,'AbsTol',1e-6,'RelTol',1e-3)*i.*t/pi);  % voor als je nauwkeurigheid wilt sturen
    end
    function y = Y2()    % integral in d/dt Ksplit
            y = (integral(@(s) dKintg(s),0,1,'AbsTol',1e-12,'RelTol',1e-8 ))*i/pi ; % exp(integral(@(s) dKintg(s),0,1,'AbsTol',1e-6,'RelTol',1e-3)*i.*t/pi);  % voor als je nauwkeurigheid wilt sturen
    end
    function y = Kintg(s)
             s1 = [s~=0].*s;
             tau = xi(zt(s1));
             y = [s==0].*0 + [s~=0].*log(LL(n,k*v(tau),hh))./(tau.^2-t^2).*dxi(zt(s1)).*dzt(s1);   % LL is better than L because of the Bessel functions
    end
    function y = dKintg(s)
             s1 = [s~=0].*s;
             tau = xi(zt(s1));
             y = [s==0].*0 + [s~=0].*log(LL(n,k*v(tau),hh)).*(tau.^2+t^2)./((tau.^2-t^2).^2).*dxi(zt(s1)).*dzt(s1);   % LL is better than L because of the Bessel functions
    end
    function y = xi(u)
             v = u/t0;
             y = u + i*dd*qq.*v./(qq-1+v.^qq);
    end
    function y = dxi(u)
             v = (u/t0).^qq;
             y = 1 + i*dd/t0*qq*(qq-1)*(1-v)./((qq-1+v).^2);
    end
    function y = zt(s)
             y = s./((1-s).^2);
    end
    function y = dzt(s)
             y = (1+s)./((1-s).^3);
    end
end

function testKm(M,dd,qq,epsfn,ddx1,ddy1,ddx2,ddy2,pprsz,ppltfnt)
    s = [0.0 :0.001:  1.   ];
    s = s(2:length(s));                                                                                 % s(1: ???
    movegui(figure(10),'southeast'), clf
    t0 = 1; plot(real(xi(zt(s))),imag(xi(zt(s))),'-k')   % integration path
    hold on
    plot(1,0,'.k','Markersize',10);
    axis([0,4,-0.5*dd,2*dd]);
    grid on
    if M~=0
      t0 = 1/M; plot(real(xi(zt(s))),imag(xi(zt(s))),'-r')
      plot(1/M,0,'.b','Markersize',10);
      axis([0,3/M,-0.5*dd,2*dd]);
    end
    %title('Integration path','Fontweight','normal')
    hold off
    drawnow
    pax = gca;
    pax.FontSize = ppltfnt;
    set(gcf,'color',[1 1 1])
    set(gcf,'inverthardcopy','off')
    set(gcf,'papertype','a4letter')
    set(gcf,'paperunits','centimeters')
    set(gcf,'paperposition',pprsz)
    set(gcf,'clipping','off')
    eval(['print -depsc -cmyk ' epsfn]);
    fixepsbbox(epsfn,ddx1,ddy1,ddx2,ddy2)
    function y = xi(u)
             v = u/t0;
             y = u + i*dd*qq.*v./(qq-1+v.^qq);
    end
    function y = dxi(u)
             v = (u/t0).^qq;
             y = 1 + i*dd/t0*qq*(qq-1)*(1-v)./((qq-1+v).^2);
    end
    function y = zt(s)
             y = s./((1-s).^2);
    end
    function y = dzt(s)
             y = (1+s)./((1-s).^3);
    end
end
%
function y = LL(n,x,hh) % the same for all h, incl. h=0
         t = i*x                            ;
         if hh ==0
            y = -2*exp(-t+abs(real(t))).*dK(n,t,1).*dI(n,t,1);   % Merk op: -2 K' I'
         else
            y = -2*exp(-t+abs(real(t))).*dK(n,t,1).*dI(n,t,1) + 2*exp(-2*t+t*hh+abs(real(t))*hh).*dK(n,t,1).^2.*dI(n,t*hh,1)./dK(n,t*hh,1);
         end
%         l4 =  (-16*n^4+184*n^2-45)/128;
%         l2 =  (4*n^2-3)/8;
%         t  = i*x;
%         t2 = t.^2;
%         y  = (1 + (l2 + l4./t2)./t2)./t;
end
function y = dLL(n,x,hh)
         t = i*x    ;
         if hh==0
            y1 = LL(n,x,0)./t;
            y2 = (1+(n./t).^2).*(dI(n,t).*K(n,t)+I(n,t).*dK(n,t));
            y  = -2*i*(y1+y2);
         else
            y1 = (d2K(n,t)./dK(n,t)).*LL(n,x,hh);
            y2 = (1+(n./t).^2)./t               ;
            y3 = (dK(n,t)./dK(n,t*hh)).^2.*(1+(n./t./hh).^2)./t;
            y  = 2*i*(y1+y2-y3);
         end
end

function y = J(n,x)
  y = besselj(n,x);
end
function y = dJ(n,x,scale)
  switch nargin
  case 2
    if n==0
       y = -besselj(1,x);
    else
       y = (besselj(n-1,x)-besselj(n+1,x))/2;
    end
  case 3
    if n==0
       y = -besselj(1,x,1);
    else
       y = (besselj(n-1,x,1)-besselj(n+1,x,1))/2;
    end
  end
end
function y = d2J(n,x)
         y = -dJ(n,x)./x - (1-(n./x).^2).*J(n,x);
end
function y = Y(n,x)
  y = bessely(n,x);
end
function y = dY(n,x)
  if n==0
     y = -bessely(1,x);
  else
     y = (bessely(n-1,x)-bessely(n+1,x))/2;
  end
end
function y = d2Y(n,x)
         y = -dY(n,x)./x - (1-(n./x).^2).*Y(n,x);
end
function y = H1(n,x)
  y = besselh(n,1,x);
end
function y = dH1(n,x)
  if n==0
     y = -besselh(1,1,x);
  else
     y = (besselh(n-1,1,x)-besselh(n+1,1,x))/2;
  end
end
function y = d2H1(n,x)
         y = -dH1(n,x)./x - (1-(n./x).^2).*H1(n,x);
end
function y = H2(n,x)
  y = besselh(n,2,x,1);
end
function y = dH2(n,x)
  if n==0
     y = -besselh(1,2,x);
  else
     y = (besselh(n-1,2,x)-besselh(n+1,2,x))/2;
  end
end
function y = d2H2(n,x)
         y = -dH2(n,x)./x - (1-(n./x).^2).*H2(n,x);
end
function y = K(n,x)
         y =besselk(n,x);
end
function y = dK(n,x,scale)
      switch nargin
      case 2
        if n==0
           y = -besselk(1,x);
        else
           y = -(besselk(n-1,x)+besselk(n+1,x))/2;
        end
      case 3
        if n==0
           y = -besselk(1,x,1);
        else
           y = -(besselk(n-1,x,1)+besselk(n+1,x,1))/2;
        end
      end
end
function y = d2K(n,x)
         y = -dK(n,x)./x + (1+(n./x).^2).*K(n,x);
end
function y = I(n,x)
         y =besseli(n,x);
end
function y = dI(n,x,scale)
     switch nargin
     case 2
       if n==0
           y = besseli(1,x);
       else
           y = (besseli(n-1,x)+besseli(n+1,x))/2;
       end
     case 3
       if n==0
           y = besseli(1,x,1);
       else
           y = (besseli(n-1,x,1)+besseli(n+1,x,1))/2;
       end
     end
end
function y = d2I(n,x)
         y = -dI(n,x)./x + (1+(n./x).^2).*I(n,x);
end
function y = vp(t)
         y = sqrt(1-t);      % principal value
end
function y = vm(t)
         y = sqrt(1+t);      % principal value
end
function y = v(t)            % y = vp(t).*vm(t);
         y = sqrt(1-t.^2);   % principal value
end
function y = vomb(x)         % auxiliary functie for Ombar
         y = [1*(abs(x)<1)-i*(abs(x)>1)].*sqrt(abs(1-x.^2));
end
function myprintf(fid,s,x)
         if abs(x)>1e-9
            switch myceil(x)
            case  4,  fprintf(fid,[s,'%+17.12f'],x);
            case  3,  fprintf(fid,[s,'%+17.13f'],x);
            case  2,  fprintf(fid,[s,'%+17.14f'],x);
            otherwise fprintf(fid,[s,'%+17.15f'],x);
            end
         else
            fprintf(fid,[s,'%+17.11e'],x);
         end
end
function y=myceil(x)
         y = ceil(log10(x))+[x==0.1]+[x==1]+[x==10]+[x==100]+[x==1000]+[x==10000];  % number of digits of x
end
%
%% routine used for GL integration, fetched from elsewhere
function [x,w]=legzo(n, a, b)
%       =========================================================
%       Purpose : Compute the zeros of Legendre polynomial Pn(x)
%                 in the interval [a,b], and the corresponding
%                 weighting coefficients for Gauss-Legendre
%                 integration
%       Input :   n    --- Order of the Legendre polynomial
%                 a    --- Lower boundary (optional)
%                 b    --- Upper boundary (optional)
%       Output:   x(n) --- Zeros of the Legendre polynomial
%                 w(n) --- Corresponding weighting coefficients
%       =========================================================
  if nargin == 1
      a = -1;
      b =  1;
  end;
  x = zeros(1, n);
  w = zeros(1, n);
  m = (n+1)/2;
  h = b-a;

  for ii=1:m
      z = cos(pi*(ii-.25)/(n+.5)); % Initial estimate.
      z1 = z+1;
      while abs(z-z1)>eps
          p1 = 1;
          p2 = 0;
          for jj = 1:n
              p3 = p2;
              p2 = p1;
              p1 = ((2*jj-1)*z*p2-(jj-1)*p3)/jj; % The Legendre polynomial.
          end
          pp = n*(z*p1-p2)/(z^2-1); % The L.P. derivative.
          z1 = z;
          z = z1-p1/pp;
      end
      x(ii) = z; % Build up the abscissas.
      x(n+1-ii) = -z;
      w(ii) = h/((1-z^2)*(pp^2)); % Build up the weights.
      w(n+1-ii) = w(ii);
  end

  if a ~= -1 || b ~= 1
      x = (x+1)*(h/2) + a;
  end
end

%% routine fetched from elsewhere, to deplete the excess of plot points
       % dperr = 1e-6;
       % ps = dpsimplify([x(1:N).',y(1:N).'],dperr);   %% LET OP: verticale kolommen
       % plot(ps(:,1),ps(:,2),'r','LineWidth',2);
function [ps,ix] = dpsimplify(p,tol)
   % Recursive Douglas-Peucker Polyline Simplification, Simplify
   %
   % [ps,ix] = dpsimplify(p,tol)
   %
   % dpsimplify uses the recursive Douglas-Peucker line simplification
   % algorithm to reduce the number of vertices in a piecewise linear curve
   % according to a specified tolerance. The algorithm is also know as
   % Iterative Endpoint Fit. It works also for polylines and polygons
   % in higher dimensions.
   %
   % In case of nans (missing vertex coordinates) dpsimplify assumes that
   % nans separate polylines. As such, dpsimplify treats each line
   % separately.
   %
   % For additional information on the algorithm follow this link
   % http://en.wikipedia.org/wiki/Ramer-Douglas-Peucker_algorithm
   %
   % Input arguments
   %
   %     p     polyline n*d matrix with n vertices in d
   %           dimensions.
   %     tol   tolerance (maximal euclidean distance allowed
   %           between the new line and a vertex)
   %
   % Output arguments
   %
   %     ps    simplified line
   %     ix    linear index of the vertices retained in p (ps = p(ix))
   %
   % Examples
   %
   % 1. Simplify line
   %
   %     tol    = 1;
   %     x      = 1:0.1:8*pi;
   %     y      = sin(x) + randn(size(x))*0.1;
   %     p      = [x' y'];
   %     ps     = dpsimplify(p,tol);
   %
   %     plot(p(:,1),p(:,2),'k')
   %     hold on
   %     plot(ps(:,1),ps(:,2),'r','LineWidth',2);
   %     legend('original polyline','simplified')
   %
   % 2. Reduce polyline so that only knickpoints remain by
   %    choosing a very low tolerance
   %
   %     p = [(1:10)' [1 2 3 2 4 6 7 8 5 2]'];
   %     p2 = dpsimplify(p,eps);
   %     plot(p(:,1),p(:,2),'k+--')
   %     hold on
   %     plot(p2(:,1),p2(:,2),'ro','MarkerSize',10);
   %     legend('original line','knickpoints')
   %
   % 3. Simplify a 3d-curve
   %
   %     x = sin(1:0.01:20)';
   %     y = cos(1:0.01:20)';
   %     z = x.*y.*(1:0.01:20)';
   %     ps = dpsimplify([x y z],0.1);
   %     plot3(x,y,z);
   %     hold on
   %     plot3(ps(:,1),ps(:,2),ps(:,3),'k*-');
   %
   % Author: Wolfgang Schwanghart, 13. July, 2010.
   % w.schwanghart[at]unibas.ch
   if nargin == 0
       help dpsimplify
       return
   end
   error(nargchk(2, 2, nargin))
   % error checking
   if ~isscalar(tol) || tol<0;
       error('tol must be a positive scalar')
   end
   % nr of dimensions
   nrvertices    = size(p,1);
   dims    = size(p,2);
   % anonymous function for starting point and end point comparision
   % using a relative tolerance test
   compare = @(a,b) abs(a-b)/max(abs(a),abs(b)) <= eps;
   % what happens, when there are NaNs?
   % NaNs divide polylines.
   Inan      = any(isnan(p),2);
   % any NaN at all?
   Inanp     = any(Inan);
   % if there is only one vertex
   if nrvertices == 1 || isempty(p);
       ps = p;
       ix = 1;
   % if there are two
   elseif nrvertices == 2 && ~Inanp;
       % when the line has no vertices (except end and start point of the
       % line) check if the distance between both is less than the tolerance.
       % If so, return the center.
       if dims == 2;
           d    = hypot(p(1,1)-p(2,1),p(1,2)-p(2,2));
       else
           d    = sqrt(sum((p(1,:)-p(2,:)).^2));
       end
       if d <= tol;
           ps = sum(p,1)/2;
           ix = 1;
       else
           ps = p;
           ix = [1;2];
       end
   elseif Inanp;
       % case: there are nans in the p array
       % --> find start and end indices of contiguous non-nan data
       Inan = ~Inan;
       sIX = strfind(Inan',[0 1])' + 1;
       eIX = strfind(Inan',[1 0])';
       if Inan(end)==true;
           eIX = [eIX;nrvertices];
       end
       if Inan(1);
           sIX = [1;sIX];
       end
       % calculate length of non-nan components
       lIX = eIX-sIX+1;
       % put each component into a single cell
       c   = mat2cell(p(Inan,:),lIX,dims);
       % now call dpsimplify again inside cellfun.
       if nargout == 2;
           [ps,ix]   = cellfun(@(x) dpsimplify(x,tol),c,'uniformoutput',false);
           ix        = cellfun(@(x,six) x+six-1,ix,num2cell(sIX),'uniformoutput',false);
       else
           ps   = cellfun(@(x) dpsimplify(x,tol),c,'uniformoutput',false);
       end
       % write the data from a cell array back to a matrix
       ps = cellfun(@(x) [x;nan(1,dims)],ps,'uniformoutput',false);
       ps = cell2mat(ps);
       ps(end,:) = [];
       % ix wanted? write ix to a matrix, too.
       if nargout == 2;
           ix = cell2mat(ix);
       end
   else
       % if there are no nans than start the recursive algorithm
       ixe     = size(p,1);
       ixs     = 1;
       % logical vector for the vertices to be retained
       I   = true(ixe,1);
       % call recursive function
       p   = simplifyrec(p,tol,ixs,ixe);
       ps  = p(I,:);
       % if desired return the index of retained vertices
       if nargout == 2;
           ix  = find(I);
       end
   end
   % _________________________________________________________
   function p  = simplifyrec(p,tol,ixs,ixe)
      % check if startpoint and endpoint are the same
      % better comparison needed which included a tolerance eps
      c1 = num2cell(p(ixs,:));
      c2 = num2cell(p(ixe,:));
      % same start and endpoint with tolerance
      sameSE = all(cell2mat(cellfun(compare,c1(:),c2(:),'UniformOutput',false)));
      if sameSE;
          % calculate the shortest distance of all vertices between ixs and
          % ixe to ixs only
          if dims == 2;
              d    = hypot(p(ixs,1)-p(ixs+1:ixe-1,1),p(ixs,2)-p(ixs+1:ixe-1,2));
          else
              d    = sqrt(sum(bsxfun(@minus,p(ixs,:),p(ixs+1:ixe-1,:)).^2,2));
          end
      else
          % calculate shortest distance of all points to the line from ixs to ixe
          % subtract starting point from other locations
          pt = bsxfun(@minus,p(ixs+1:ixe,:),p(ixs,:));
          % end point
          a = pt(end,:)';
          beta = (a' * pt')./(a'*a);
          b    = pt-bsxfun(@times,beta,a)';
          if dims == 2;
              % if line in 2D use the numerical more robust hypot function
              d    = hypot(b(:,1),b(:,2));
          else
              d    = sqrt(sum(b.^2,2));
          end
      end
      % identify maximum distance and get the linear index of its location
      [dmax,ixc] = max(d);
      ixc  = ixs + ixc;
      % if the maximum distance is smaller than the tolerance remove vertices
      % between ixs and ixe
      if dmax <= tol;
          if ixs ~= ixe-1;
              I(ixs+1:ixe-1) = false;
          end
          % if not, call simplifyrec for the segments between ixs and ixc (ixc and ixe)
      else
          p   = simplifyrec(p,tol,ixs,ixc);
          p   = simplifyrec(p,tol,ixc,ixe);
      end
   end
end

%% van elders opgeduikelde routine om fout in eps te corrigeren
function fixepsbbox(filename,dx1,dy1,dx2,dy2)
% function fixepsbbox(filename)
%
% matlab seems to compute a bounding box on eps files which is too
% large in the x-direction
%
% this script fixes the bounding box
% it is 99% stolen from fixeps.m, located here:
% http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=4818&objectType=file
% the only change is that this changes the bbox numbers
% seriously, the only change is the addition of lines 22,23,33,34
%
% boundingbox has form of:
% %%BoundingBox:    x1   y1   x2   y2
% where (x1,y1) is lower-left and (x2,y2) is upper-right
%
% matlab computes x1 too small and x2 too large
% changes lines to:
% %%BoundingBox:    x1+dx1 y1 x2+dx2 y2
% default amount to change bbox - found this fixed my plots just fine
%dx1 = -10;	% amount to move x1
%dx2 =  25;	% amount to move x2
fid = fopen(filename,'r+');
k=0;
while k <2                                                  % 2 locations to replace.
    tline = fgetl(fid);                                     % get one line text
    stridx=strfind(tline,'Box:');
	if isempty(stridx)==0
        len=length(tline);                                  % the original line length
		bb=sscanf(tline(stridx+4:end),'%i');                % read the numbers
		bb(1) = bb(1) + dx1;								% change x1
		bb(2) = bb(2) + dy1;								% change y1
		bb(3) = bb(3) + dx2;								% change x2
		bb(4) = bb(4) + dy2;								% change y2
		bbstr=sprintf('%g %g %g %g',bb);                    % write bb numbers to string
        tline=tline(1:stridx+3);                             % keep the "%%(page)boundingbox" string (with starting '%%')
		spaces(1:len-length(tline)-length(bbstr)-1)=' ';    % add trailing spaces as to overwrite old line completely
		tline=[tline ' ' bbstr spaces];                     % concate numbers and blank spaces to "%%(page)boundingbox"
        fseek(fid,-len-2,'cof');                            % before using fprintf search to correct position
		count = fprintf(fid,'%s',tline);
        fseek(fid,2,'cof');                                 % seek to beginning of line (for windows text file) on
                                                            % for linux: change '2' to '1' I think
        k=k+1;
	end
end
fclose(fid);
end

