function demo13d
% **********************************************************
% *** Last changed 2012-11-10 by Johan Helsing             *
% *** Exterior Neumann problem for Helmholtz               *
% *** RCIP mehtod and regularized combined field equation  *
% *** Field evaluation                                     *
% **********************************************************
  close all
  format long
  format compact
%
% *** user specified quantities ****************************
  itmax=50          % maximum number of GMRES iterations   *
  theta=pi/2        % corner opening angle                 *
  nsub=52;          % number of levels in recursion        *
  premature=46;     % interruption level in reconstruction *
  zsource=0.3+0.1i; % source point                         *
  omega=10          % wave number                          *
  npan=56           % number of panels on coarse mesh      *
  dlim=1.1          % closeness parameter                  *
%  Ng=300;           % field evaluation at Ng^2 points      *
  Ng=100;           % field evaluation at Ng^2 points      *
% **********************************************************
  T16=Tinit16;
  W16=Winit16;
  [IP,IPW]=IPWinit(T16,W16);
  PbcBig=PbcBiginit(IP);  
  PWbcBig=PbcBiginit(IPW);    
  R0=initializeR(W16,T16,theta,PWbcBig,PbcBig);
%
% *** panel breakpoints and interval lengths in parameter ***
  sinter=linspace(0,1,npan+1);
  sinterdiff=ones(npan,1)/npan;
%
% *** discretization points and weights ***  
  npcoa=16*npan
  scoa=zeros(npcoa,1);
  wcoa=zeros(npcoa,1);
  for k=1:npan
    myind=k*16-15:k*16;
    sdif=sinterdiff(k)/2;
    scoa(myind)=(sinter(k)+sinter(k+1))/2+sdif*T16;
    wcoa(myind)=W16*sdif;
  end
  zcoa=zfunc(scoa,theta) ;
  zinter=zfunc(sinter',theta);
  zpcoa=zpfunc(scoa,theta);
  zppcoa=zppfunc(scoa,theta);
% *** some extra presicion gained from symmetry ***
  zcoa(npcoa/2+1:npcoa)=conj(flipud(zcoa(1:npcoa/2)));
  itmp=round(npan/2+1);
  zinter(itmp:npan+1)=conj(flipud(zinter(1:npan-itmp+2)));
  zpcoa(npcoa/2+1:npcoa)=-conj(flipud(zpcoa(1:npcoa/2)));
  zppcoa(npcoa/2+1:npcoa)=conj(flipud(zppcoa(1:npcoa/2)));
% *************************************************
  nzcoa=-1i*zpcoa./abs(zpcoa);
  wzpcoa=wcoa.*zpcoa;
  awzpcoa=abs(wzpcoa);
  arclength=sum(awzpcoa)
  area=0.5*imag(conj(zcoa).'*wzpcoa)
%
% *** The coarse matrices are set up ***
  disp('Setup of coarse matrices starts')
  [LogCcoa,CauCcoa]=LogCauCinit(zcoa,nzcoa,wzpcoa,zinter,sinterdiff,T16, ...
				W16,npan,npcoa,1);
  Acirccoa=-Kadjoperinit(LogCcoa,zcoa,zpcoa,zppcoa,awzpcoa,wcoa, ...
			 nzcoa,omega,npcoa);
  B2circcoa=SoperinitI(LogCcoa,zcoa,zpcoa,awzpcoa,sinterdiff,omega,npcoa);
  B1circcoa=-1i*B1operinit(LogCcoa,zcoa,zpcoa,awzpcoa,sinterdiff, ...
			   nzcoa,omega,npcoa);
  C1circcoa=-1i*SXoperinit(LogCcoa,CauCcoa,zcoa,zpcoa,nzcoa,awzpcoa, ...
			   omega,npan,npcoa);
  C2circcoa=SXoperinitI(LogCcoa,CauCcoa,zcoa,zpcoa,nzcoa,awzpcoa,omega, ...
			npan,npcoa);
  clear LogCcoa CauCcoa
  starind=[npcoa-31:npcoa 1:32];
  Acirccoa(starind,starind)=zeros(64);
  B1circcoa(starind,starind)=zeros(64);
  B2circcoa(starind,starind)=zeros(64);
  C1circcoa(starind,starind)=zeros(64);
  C2circcoa(starind,starind)=zeros(64);
%
% *** recursion for the R matrix ***
  disp('Recursion starts')
  R=R0;
  Rstor=zeros(192,192,nsub);
  ABClocstor=zeros(288,288,nsub);
  starL=[17:80 113:176 209:272];
  for level=1:nsub
    denom=2^(nsub-level)*npan;
    sloc=[T16/4+0.25;T16/4+0.75;T16/2+1.5]/denom;
    sinterloc=[0;0.5;1;2]/denom;
    zloc=zfuncloc(sloc,theta);
    zinterloc=zfuncloc(sinterloc,theta);
    zinterloc(4)=[];
    zploc=zpfuncloc(sloc,theta);
    zpploc=zppfuncloc(sloc,theta);
    wloc=[W16/2;W16/4;W16/4;W16/4;W16/4;W16/2]/denom;
    nzloc=-1i*zploc./abs(zploc);
    wzploc=wloc.*zploc;
    awzploc=wloc.*abs(zploc);
    sidiloc=[1;0.5;0.5;0.5;0.5;1]/denom;
    [LogCloc,CauCloc]=LogCauCinit(zloc,nzloc,wzploc,zinterloc,sidiloc, ...
				  T16,W16,6,96,0);
    Aloc=-Kadjoperinit(LogCloc,zloc,zploc,zpploc,awzploc,wloc,nzloc,omega,96);
    B1loc=-1i*B1operinit(LogCloc,zloc,zploc,awzploc,sidiloc,nzloc,omega,96);
    B2loc=SoperinitI(LogCloc,zloc,zploc,awzploc,sidiloc,omega,96);
    C1loc=-1i*SXoperinit(LogCloc,CauCloc,zloc,zploc,nzloc,awzploc,omega,6,96);
    C2loc=SXoperinitI(LogCloc,CauCloc,zloc,zploc,nzloc,awzploc,omega,6,96);
    ABCloc=zeros(288);
    ABCloc(  1: 96,  1: 96)= Aloc;
    ABCloc(  1: 96, 97:192)= B1loc;
    ABCloc(  1: 96,193:288)= C1loc;
    ABCloc( 97:192,  1: 96)=-B2loc;
    ABCloc(193:288,  1: 96)=-C2loc;
    ABClocstor(:,:,level)=ABCloc;
    MAT=eye(288)+ABCloc;
    Rstor(:,:,level)=R;
    R=SchurBana(PWbcBig,PbcBig,MAT,R);
  end
  R1=R(  1: 64,  1: 64);
  R2=R( 65:128,  1: 64);    
  R3=R(129:192,  1: 64);    
  R4=R(  1: 64, 65:128);
  R5=R( 65:128, 65:128);    
  R6=R(129:192, 65:128);        
  R7=R(  1: 64,129:192);
  R8=R( 65:128,129:192);    
  R9=R(129:192,129:192);        
%    
  R14=R1\R4;
  R17=R1\R7;
  R5214=R5-R2*R14;
  R8217=R8-R2*R17;
  R6314=R6-R3*R14;
  R9317=R9-R3*R17;
%
  R1G=speye(npcoa);
  R1G(starind,starind)=R1;
  R2G=sparse(npcoa);
  R2G(starind,starind)=R2;
  R3G=sparse(npcoa);
  R3G(starind,starind)=R3;
%
  R14G=sparse(npcoa);
  R14G(starind,starind)=R14;
  R17G=sparse(npcoa);
  R17G(starind,starind)=R17;
%
  R5214G=speye(npcoa);
  R5214G(starind,starind)=R5214;
  R6314G=sparse(npcoa);
  R6314G(starind,starind)=R6314;
  R8217G=sparse(npcoa);
  R8217G(starind,starind)=R8217;
  R9317G=speye(npcoa);
  R9317G(starind,starind)=R9317;
%
% *** Solving main linear system ***
  disp('Solving main linear system starts')
  tmp=omega*abs(zsource-zcoa);
  g=tmp.*besselh(1,tmp).*real(nzcoa./(zsource-zcoa));
  rhs=-2*g;
  [rhotildecoa,it]=myGMRESR5(Acirccoa,B1circcoa,B2circcoa,C1circcoa, ...
			     C2circcoa,R1G,R2G, R3G,R14G,R17G,R5214G, ...
			     R6314G,R8217G,R9317G,rhs,npcoa,itmax,eps);
  GMRES_iter_RCIP=it
%
% *** Post processing starts ***
%
  rho1hatcoa=R1G*rhotildecoa;  
  B2rho1hatcoa=B2circcoa*rho1hatcoa;
  C2rho1hatcoa=C2circcoa*rho1hatcoa;
  rho2hatcoa=R2G*rhotildecoa+R5214G*B2rho1hatcoa+R8217G*C2rho1hatcoa;
%
% *** rhofin locally ***
  nse=nsub-premature+1;
  rhotildecoa1=rhotildecoa-R14G*B2rho1hatcoa-R17G*C2rho1hatcoa;
  rhotildeBB=zeros(192,1);
  rhotildeBB(  1: 64)=rhotildecoa1(starind);
  rhotildeBB( 65:128)=B2rho1hatcoa(starind);  
  rhotildeBB(129:192)=C2rho1hatcoa(starind);    
  pts=16*(4+2*nse);
  pts2=16*(2+nse);
  rhofinlocBB=zeros(3*pts,1);
  PbcBB=PbcBiginitBig(IP);
  starL=[17:80 113:176 209:272];
  circL=[1:16 81:112 177:208 273:288];
  for k=nsub:-1:premature
    Kcirc=ABClocstor(:,:,k);
    Kcirc(starL,starL)=zeros(192);
    MAT=eye(288)+Kcirc;
    MAT(starL,starL)=inv(Rstor(:,:,k));
    myind1L=(nsub-k)*16+(1:16);
    myind2L=(nsub+k-2*premature+5)*16+(1:16);
    myind12L=[myind1L myind2L];
    myindBB=[myind12L myind12L+pts myind12L+2*pts];
    tmp=PbcBB*rhotildeBB;
    rhotildeBB=tmp-Kcirc*(MAT\tmp);
    rhofinlocBB(myindBB)=rhotildeBB(circL);
    rhotildeBB=rhotildeBB(starL);
  end
  myind12L=nse*16+1:(4+nse)*16;
  myindBB=[myind12L myind12L+pts myind12L+2*pts];  
  rhofinlocBB(myindBB)=Rstor(:,:,premature)*rhotildeBB;
  clear Klocstor Rstor
%
% *** rhofin globally ***
  rhofin1=zeros(16*(npan+2*nse),1);
  rhofin1(1:pts2)=rhofinlocBB(pts2+1:pts);
  rhofin1(pts2+1:pts2+(npan-4)*16)=rhotildecoa1(33:(npan-2)*16);
  rhofin1(pts2+(npan-4)*16+1:16*(npan+2*nse))=rhofinlocBB(1:pts2);
%
  rhofin2=zeros(16*(npan+2*nse),1);
  rhofin2(1:pts2)=rhofinlocBB(pts+(pts2+1:pts));
  rhofin2(pts2+1:pts2+(npan-4)*16)=B2rho1hatcoa(33:(npan-2)*16);
  rhofin2(pts2+(npan-4)*16+1:16*(npan+2*nse))=rhofinlocBB(pts+(1:pts2));
%
  [zfin,wzpfin,nzfin,zinterfin,npanfin]=finegrid(nse,theta,npan,W16,T16);  
%
% ***  
  Ngtot=Ng*Ng;
  xg=linspace(-0.1,1.1,Ng);
  yg=linspace(-0.53,0.53,Ng);
  zg=zeros(Ngtot,1);
  for k=1:Ng
    zg((k-1)*Ng+(1:Ng))=xg(k)+1i*yg;
  end
  Uref=besselh(0,omega*abs(zsource-zg));
%
% *** Initial screening ***
  disp('Screening starts')
  mylist=zeros(Ngtot,1); 
% *** 0 => no evaluation
% *** 1 => regular evaluation on coarse grid
% *** 2 => some special close evaluation on coarse grid
% *** 3 => some special close evaluation on fine grid
  cornpan=[1 npan];
  panlen=zeros(npan,1);
  for k=1:npan
    panlen(k)=sum(awzpcoa((k-1)*16+1:k*16));
  end
  for k=1:Ngtot
    inout=Ineval(zg(k),theta);
    if inout
      mylist(k)=1;
      for kk=1:npan
	myind=(kk-1)*16+1:kk*16;
	d=min(abs(zcoa(myind)-zg(k)))/panlen(kk);
	if d<dlim && mylist(k)<=2
	  mylist(k)=2;
	  if isempty(find(kk==cornpan))==0
	    mylist(k)=3;
	  end
	end
      end
    end
  end
  mylist1=find(mylist==1)';
  mylist2=find(mylist==2)';
  mylist3=find(mylist==3)';
  evalpoints=length(mylist1)+length(mylist2)+length(mylist3)
  a1=eps*ones(Ngtot,1);
%
  disp('List1 starts')
  for k=mylist1
    Unum=Starg(zg(k),zcoa,awzpcoa,omega).'*rho1hatcoa+ ...
	   1i*(Ktarg(zg(k),zcoa,wzpcoa,omega).'*rho2hatcoa);
    a1(k)=abs(Uref(k)-Unum);  
  end
%
  disp('List2 starts')    
  for k=mylist2
    ztarg=zg(k);
    Unum=0;
    for kk=1:npan
      myind=(kk-1)*16+1:kk*16;
      zsc=zcoa(myind);
      wzp=wzpcoa(myind);
      awzp=abs(wzp);
      rho1hat=rho1hatcoa(myind);
      rho2hat=rho2hatcoa(myind);      
      d=min(abs(zsc-ztarg))/panlen(kk);
      if d>dlim
	Unum=Unum+Starg(ztarg,zsc,awzp,omega).'*rho1hat+ ...
	      1i*(Ktarg(ztarg,zsc,wzp,omega).'*rho2hat);
      else
	a=zinter(kk);
	b=zinter(kk+1);
	zmid=b+a;
	zdif=b-a;     
	cc=zdif/2;
	nz=nzcoa(myind);
	ztgtr=(2*ztarg-zmid)/zdif;
	zsctr=(2*zsc-zmid)/zdif;
	LogC=LGIcompRecFAS(ztgtr,zsctr,ztarg,zsc,nz,awzp,cc);	  
	ztup=b-ztarg;
	ztlo=a-ztarg;
	CauC=M1IcompRecFS(ztgtr,ztup,ztlo,zsctr,wzp./(zsc-ztarg));
	Unum=Unum+StargClose(LogC,ztarg,zsc,awzp,omega).'*rho1hat+ ...
	      1i*(KtargClose(LogC,CauC,ztarg,zsc,wzp,omega).'*rho2hat);
      end
    end
    a1(k)=abs(Uref(k)-Unum);
  end
%
  disp('List3 starts')    
  for k=mylist3
    ztarg=zg(k);
    Unum=0;
    for kk=1:npanfin
      myind=(kk-1)*16+1:kk*16;
      zsc=zfin(myind);
      wzp=wzpfin(myind);
      awzp=abs(wzp);
      rho1hat=rhofin1(myind);
      rho2hat=rhofin2(myind);      
      d=min(abs(zsc-ztarg))/sum(awzp);
      if d>dlim
	Unum=Unum+Starg(ztarg,zsc,awzp,omega).'*rho1hat+ ...
	      1i*(Ktarg(ztarg,zsc,wzp,omega).'*rho2hat);
      else
	a=zinterfin(kk);
	b=zinterfin(kk+1);
	zmid=b+a;
	zdif=b-a;     
	cc=zdif/2;
	nz=nzfin(myind);
	ztgtr=(2*ztarg-zmid)/zdif;
	zsctr=(2*zsc-zmid)/zdif;
	LogC=LGIcompRecFAS(ztgtr,zsctr,ztarg,zsc,nz,awzp,cc);	  
	ztup=b-ztarg;
	ztlo=a-ztarg;
	CauC=M1IcompRecFS(ztgtr,ztup,ztlo,zsctr,wzp./(zsc-ztarg));
	Unum=Unum+StargClose(LogC,ztarg,zsc,awzp,omega).'*rho1hat+ ...
	      1i*(KtargClose(LogC,CauC,ztarg,zsc,wzp,omega).'*rho2hat);
      end
    end
    a1(k)=abs(Uref(k)-Unum);
  end
%
  a1(a1<eps)=eps;
  F1=zeros(Ng);
  for k=1:Ng
    F1(1:Ng,k)=log10(a1((k-1)*Ng+(1:Ng)));
  end
%
  set(0,'DefaultAxesFontSize',16)
  set(0,'DefaultAxesLineWidth',1) 
  figure(1)
  imagesc(xg,yg,F1);      
  colormap(flipud(pink))
  axis xy
  colorbar
  hold on
  title('Neumann: Log_{10} of estimated error in {\itU}({\itr})')
  xlabel('{\itx}','FontSize',18)
  ylabel('{\ity}','FontSize',18)
  plot(real(zsource),imag(zsource),'b*')
  axis([-0.1 1.1 -0.4 0.4])
  axis equal
  
  function R=initializeR(W16,T16,theta,PWbc,Pbc)
  myerr=1;
  starL=[17:80 113:176 209:272];
  sloc=[T16/4+0.25;T16/4+0.75;T16/2+1.5];
  zloc=sloc*exp(-1i*theta/2);
  zploc=ones(48,1)*exp(-1i*theta/2);
  zloc=[conj(flipud(zloc));zloc];
  zploc=[-conj(flipud(zploc));zploc];
  zpploc=zeros(96,1);
  wloc=[W16/2;W16/4;W16/4;W16/4;W16/4;W16/2];
  nzloc=-1i*zploc./abs(zploc);
  wzploc=wloc.*zploc;
  sinterloc=[0;0.5;1;2];
  zinterloc=sinterloc*exp(-1i*theta/2);
  zinterloc=[conj(flipud(zinterloc));zinterloc];  
  zinterloc(4)=[];
  Aloc=-M1Ainit(zloc,zploc,zpploc,nzloc,wloc,wzploc,96);
  Cloc=CauCinit0(zloc,wzploc,nzloc,zinterloc,T16,W16,6,96);
  ABCloc=zeros(288);
  ABCloc(  1: 96,  1: 96)= Aloc;
  ABCloc(  1: 96,193:288)= Cloc;
  ABCloc(193:288,  1: 96)=-1i*Cloc;
  MAT=eye(288)+ABCloc;
  R=inv(MAT(starL,starL));
  iter=0;
  while myerr>eps
    Rold=R;
    R=SchurBana(PWbc,Pbc,MAT,R);
    myerr=norm(R-Rold,'fro')/norm(R,'fro');
    iter=iter+1;
  end
  Init_iter=iter

  function [zfin,wzpfin,nzfin,zinterfin,npanfin]=finegrid(nsub, ...
						  theta,npan,W16,T16)  
% *** panel breakpoints and interval lengths in parameter ***
  npanfin=npan+2*nsub;
  sinterfin=zeros(npanfin+1,1);
  sinterfin(1:npan+1)=linspace(0,1,npan+1);
  sinterfindiff=ones(npan+2*nsub,1)/npan;
  for k=1:nsub
    sinterfin(3:end)=sinterfin(2:end-1);
    sinterfin(2)=(sinterfin(1)+sinterfin(2))/2;   
    sinterfindiff(3:end)=sinterfindiff(2:end-1);
    sinterfindiff(2)=sinterfindiff(1)/2;
    sinterfindiff(1)=sinterfindiff(1)/2;   
  end
  sinterfin(end-nsub:end)=1-flipud(sinterfin(1:nsub+1));
  sinterfindiff(end-nsub-1:end)=flipud(sinterfindiff(1:nsub+2));
%
% *** discretization points and weights ***  
  npfin=16*npanfin;
  sfin=zeros(npfin,1);
  wfin=zeros(npfin,1);
  for k=1:npanfin
    myind=k*16-15:k*16;
    sdif=sinterfindiff(k)/2;
    sfin(myind)=(sinterfin(k)+sinterfin(k+1))/2+sdif*T16;
    wfin(myind)=W16*sdif;
  end
  zfin=zfunc(sfin,theta) ;
  zinterfin=zfunc(sinterfin,theta);
  zpfin=zpfunc(sfin,theta);
% *** some extra presicion gained from symmetry ***
  zfin(npfin/2+1:npfin)=conj(flipud(zfin(1:npfin/2)));
  itmp=round(npanfin/2+1);
  zinterfin(itmp:npanfin+1)=conj(flipud(zinterfin(1:npanfin-itmp+2)));
  zpfin(npfin/2+1:npfin)=-conj(flipud(zpfin(1:npfin/2)));
% *************************************************
  nzfin=-1i*zpfin./abs(zpfin);
  wzpfin=wfin.*zpfin;

  function iout=Ineval(zg,theta)
  iout=0;
  thetag=abs(angle(zg));
  if thetag>theta/2
    iout=1;
  else
    sg=0.5+thetag/theta;
    if abs(zg)>sin(pi*sg)
      iout=1;
    end
  end
  
  function M1=M1Ainit(z,zp,zpp,nz,w,wzp,N)
% *** adjoint of double layer potential ***   
  M1=zeros(N);
  for m=1:N
    M1(:,m)=abs(wzp(m))*real(nz./(z(m)-z));
    M1(m,m)=-w(m)*imag(zpp(m)/zp(m))/2;      
  end
  M1=M1/pi;
  
  function S=Soperinit(LogC,z,zp,awzp,sinterdiff,omega,N)
  S=zeros(N);
  for m=1:N
    tmp=omega*abs(z-z(m));
    S(:,m)=besselh(0,tmp)*awzp(m);
  end
  myind=find(LogC);
  S(myind)=S(myind)+2i/pi*(real(S(myind)).*LogC(myind));
  tmp=1+2i/pi*(log(omega/2)-psi(1));
  for m=1:N
    dsd=sinterdiff(fix((m-1)/16)+1)/2;
    S(m,m)=(2i/pi*(LogC(m,m)+log(dsd*abs(zp(m))))+tmp)*awzp(m);
  end
  S=1i/2*S;

  function S=SoperinitI(LogC,z,zp,awzp,sinterdiff,omega,N)
  S=zeros(N);
  for m=1:N
    tmp=1i*omega*abs(z-z(m));
    S(:,m)=besselh(0,tmp)*awzp(m);
  end
  [myindi,myindj]=find(LogC);
  tmp=1i*omega*abs(z(myindi)-z(myindj));
  TMP=besselj(0,tmp).*awzp(myindj);
  myind=find(LogC);
  S(myind)=S(myind)+2i/pi*(TMP.*LogC(myind));
  tmp=1+2i/pi*(log(1i*omega/2)-psi(1));
  for m=1:N
    dsd=sinterdiff(fix((m-1)/16)+1)/2;
    S(m,m)=(2i/pi*(LogC(m,m)+log(dsd*abs(zp(m))))+tmp)*awzp(m);
  end
  S=1i/2*S;

  function SX=SXoperinit(LogC,CauC,z,zp,nz,awzp,omega,nseg,N)
  SX=zeros(N);
  for m=1:N
    tmp=omega*abs(z-z(m));
    SX(:,m)=tmp.*besselh(1,tmp).*imag(nz./(z-z(m)))*awzp(m);
    SX(m,m)=0;
  end
  myind=find(LogC);
  SX(myind)=SX(myind)+2i/pi*(real(SX(myind)).*LogC(myind));
  myind=find(CauC);
  SX(myind)=SX(myind)+CauC(myind);  
  SX=1i/2*SX;

  function SX=SXoperinitI(LogC,CauC,z,zp,nz,awzp,omega,nseg,N)
  SX=zeros(N);
  for m=1:N
    tmp=1i*omega*abs(z-z(m));
    SX(:,m)=tmp.*besselh(1,tmp).*imag(nz./(z-z(m)))*awzp(m);
    SX(m,m)=0;
  end
  myind=find(LogC);
  [mi,mj]=find(LogC);
  mij=find(mi==mj);
  tmp=1i*omega*abs(z(mi)-z(mj));
  TMP=tmp.*besselj(1,tmp).*imag(nz(mi)./(z(mi)-z(mj))).*awzp(mj);
  TMP(mij)=0;
  SX(myind)=SX(myind)+2i/pi*(TMP.*LogC(myind));  
  myind=find(CauC);
  SX(myind)=SX(myind)+CauC(myind);
  SX=1i/2*SX;
  
  function B1=B1operinit(LogC,z,zp,awzp,sinterdiff,nz,omega,N)  
  B1=omega^2*Soperinit(LogC,z,zp,awzp,sinterdiff,omega,N);
  for m=1:N
    B1(:,m)=B1(:,m).*real(nz*conj(nz(m)));
  end

  function Kadj=Kadjoperinit(LogC,z,zp,zpp,awzp,w,nz,omega,N)
  Kadj=zeros(N);
  for m=1:N
    tmp=omega*abs(z-z(m));
    Kadj(:,m)=tmp.*besselh(1,tmp).*real(nz./(z(m)-z))*awzp(m);
    Kadj(m,m)=1i/pi*imag(zpp(m)/zp(m))*w(m);        
  end
  myind=find(LogC);
  Kadj(myind)=Kadj(myind)+2i/pi*(real(Kadj(myind)).*LogC(myind));
  Kadj=1i/2*Kadj;
  
  function u=Starg(ztarg,z,awzp,omega)
  u=1i/4*besselh(0,omega*abs(ztarg-z)).*awzp;  
  
  function u=Ktarg(ztarg,z,wzp,omega)
  tmp=omega*abs(z-ztarg);
  u=1i/4*tmp.*besselh(1,tmp).*imag(wzp./(ztarg-z));

  function u=StargClose(LogC,ztarg,z,awzp,omega)
  u=besselh(0,omega*abs(ztarg-z)).*awzp;
  u=u+2i/pi*(real(u).*LogC);
  u=1i/4*u;
  
  function u=KtargClose(LogC,CauC,ztarg,z,wzp,omega)
  tmp=omega*abs(z-ztarg);
  u=tmp.*besselh(1,tmp).*imag(wzp./(ztarg-z));
  u=u+2i/pi*(real(u).*LogC+CauC);  
  u=1i/4*u;
  
  function A=SchurBana(PW,P,K,A)
  starL=[17:80 113:176 209:272];
  circL=[1:16 81:112 177:208 273:288];
  starS=[17:48 81:112 145:176];
  circS=[1:16 49:80 113:144 177:192];
  VA=K(circL,starL)*A;
  PTA=PW'*A;
  PTAU=PTA*K(starL,circL);
  DVAUI=inv(K(circL,circL)-VA*K(starL,circL));
  DVAUIVAP=DVAUI*(VA*P);
  A(starS,starS)=PTA*P+PTAU*DVAUIVAP;
  A(circS,circS)=DVAUI;
  A(circS,starS)=-DVAUIVAP;
  A(starS,circS)=-PTAU*DVAUI;
  
  function [x,it]=myGMRESR5(A,B1,B2,C1,C2,R1,R2,R3,R14,R17,R5214,R6314, ...
			    R8217,R9317,b,n,m,tol);
% *** GMRES with low-threshold stagnation control ***
  V=zeros(n,m+1);
  H=zeros(m);
  cs=zeros(m,1);
  sn=zeros(m,1);
  bnrm2=norm(b);
  V(:,1)=b/bnrm2;
  s=bnrm2*eye(m+1,1);
  for it = 1:m                                  
    it1=it+1;                                   
    x=V(:,it);
    R1x=R1*x;
    B2R1x=B2*R1x;
    C2R1x=C2*R1x;
    w=A*R1x-R14*B2R1x-R17*C2R1x+B1*(R2*x+R5214*B2R1x+R8217*C2R1x)+ ...
      C1*(R3*x+R6314*B2R1x+R9317*C2R1x);
    for k = 1:it
      H(k,it)=V(:,k)'*w;
      w=w-H(k,it)*V(:,k);
    end
    H(it,it)=H(it,it)+1;
    wnrm2=norm(w);
    V(:,it1)=w/wnrm2;
    for k=1:it-1                                
      temp     = cs(k)*H(k,it)+sn(k)*H(k+1,it);
      H(k+1,it)=-sn(k)*H(k,it)+cs(k)*H(k+1,it);
      H(k,it)  = temp;
    end
    [cs(it),sn(it)]=rotmat(H(it,it),wnrm2);     
    H(it,it)= cs(it)*H(it,it)+sn(it)*wnrm2;
    s(it1) =-sn(it)*s(it);                      
    s(it)  = cs(it)*s(it);                         
    myerr=abs(s(it1))/bnrm2;
    if (myerr<=tol)|(it==m)                     
      predicted_residual=myerr
      y=triu(H(1:it,1:it))\s(1:it);             
      x=fliplr(V(:,1:it))*flipud(y);
      R1x=R1*x;
      B2R1x=B2*R1x;
      C2R1x=C2*R1x;
      w=-R14*B2R1x-R17*C2R1x+A*R1x+B1*(R2*x+R5214*B2R1x+R8217*C2R1x)+ ...
	C1*(R3*x+R6314*B2R1x+R9317*C2R1x);
      true_residual=norm(x+w-b)/bnrm2
      break;
    end
  end

  function [c,s]=rotmat(a,b)
  if  b==0
    c=1;
    s=0;
  elseif abs(b)>abs(a)
    temp=a/b;
    s=1/sqrt(1+temp^2);
    c=temp*s;
  else
    temp=b/a;
    c=1/sqrt(1+temp^2);
    s=temp*c;
  end

  function zout=zfunc(s,theta)
  zout=sin(pi*s).*exp(1i*theta*(s-0.5));

  function zpout=zpfunc(s,theta)
  zpout=(pi*cos(pi*s)+1i*theta*sin(pi*s)).*exp(1i*theta*(s-0.5));
  
  function zppout=zppfunc(s,theta)
  zppout=(2i*pi*theta*cos(pi*s)-(theta^2+pi^2)*sin(pi*s)).* ...
	 exp(1i*theta*(s-0.5));

  function zout=zfuncloc(s,theta)
  zout=zfunc(s,theta);
  zout=[conj(flipud(zout));zout];
  
  function zpout=zpfuncloc(s,theta)
  zpout=zpfunc(s,theta);
  zpout=[-conj(flipud(zpout));zpout];
  
  function zppout=zppfuncloc(s,theta)
  zppout=zppfunc(s,theta);
  zppout=[conj(flipud(zppout));zppout];

  function [ML,M1]=LogCauCinit(z,nz,wzp,zinter,pinterdiff,T16,W16,nseg,N,iper)
% *** Local panelwise evaluation ***
% *** iper=0,1 (0 is open arc, 1 is closed contour) ***
  [TMP,~]=LGIcompRecR(0,1,T16);
  TMPL=diadivR(TMP,W16);
  ML=zeros(N);
  M1=zeros(N);
  M0=zeros(16);   
  if iper==1
    kstart=1;
    kend=nseg;  
  else
    kstart=2;
    kend=nseg-1;  
  end
% *** central blocks ***
  for k=1:nseg
    myind=k*16-15:k*16;
    for nj=1:16
      m=myind(nj);
      ML(myind,m)=-log(abs(T16(nj)-T16));
      ML(m,m)=0;
    end
    ML(myind,myind)=ML(myind,myind)+TMPL;
%
    z1=zinter(k);
    z2=zinter(k+1);
    zmid=z2+z1;
    zdif=z2-z1;     
    zs=(2*z(myind)-zmid)/zdif;
    zt=(2*z(myind)-zmid)/zdif;
    ztup=z2-z(myind);
    ztlo=z1-z(myind);
    [TMPM,~]=M1IcompRecC(zt,ztup,ztlo,zs,1);     
    for m=1:16
     M0(:,m)=wzp(myind(m))./(z(myind(m))-z(myind))/pi/1i;
     M0(m,m)=0;
    end
    M1(myind,myind)=TMPM/pi/1i-M0;
    for m=myind
      M1(myind,m)=conj(nz(m))*nz(myind).*M1(myind,m);
    end
  end
% *** superdiagonal blocks (targets to the left) ***
  for k=kstart:nseg
    myinds=k*16-15:k*16;
    km1=mod(k-2,nseg)+1;
    mi=km1*16-15:km1*16;       
    alpha=pinterdiff(km1)/pinterdiff(k);          
    [TMPL,~]=LGIcompRecR(-1-alpha,alpha,T16);
    for nj=1:16
      ML(mi,myinds(nj))=-log(abs(T16(nj)+1+(1-T16)*alpha));       
    end
    ML(mi,myinds)=ML(mi,myinds)+diadivR(TMPL,W16);
%  
    z1=zinter(k);
    z2=zinter(k+1);
    zmid=z2+z1;
    zdif=z2-z1;     
    zs=(2*z(myinds)-zmid)/zdif;
    zt=(2*z(mi)-zmid)/zdif;          
    ztup=z2-z(mi);
    ztlo=z1-z(mi);
    [TMPM,~]=M1IcompRecC(zt,ztup,ztlo,zs,0);     
    for m=1:16
     M0(:,m)=wzp(myinds(m))./(z(myinds(m))-z(mi))/pi/1i;
    end
    M1(mi,myinds)=TMPM/pi/1i-M0;
    for m=myinds
      M1(mi,m)=conj(nz(m))*nz(mi).*M1(mi,m);
    end
  end
% *** subdiagonal blocks (targets to the right) ***
  for k=1:kend
    myinds=k*16-15:k*16;
    myinds=k*16-15:k*16;
    kp1=mod(k,nseg)+1;
    mi=kp1*16-15:kp1*16;       
    alpha=pinterdiff(kp1)/pinterdiff(k);               
    [TMPL,~]=LGIcompRecR(1+alpha,alpha,T16);
    for nj=1:16
      ML(mi,myinds(nj))=-log(abs(T16(nj)-1-(T16+1)*alpha));   
    end
    ML(mi,myinds)=ML(mi,myinds)+diadivR(TMPL,W16);
%
    z1=zinter(k);
    z2=zinter(k+1);
    zmid=z2+z1;
    zdif=z2-z1;     
    zs=(2*z(myinds)-zmid)/zdif;
    zt=(2*z(mi)-zmid)/zdif;
    ztup=z2-z(mi);
    ztlo=z1-z(mi);
    [TMPM,~]=M1IcompRecC(zt,ztup,ztlo,zs,0);     
    for m=1:16
     M0(:,m)=wzp(myinds(m))./(z(myinds(m))-z(mi))/pi/1i;
    end
    M1(mi,myinds)=TMPM/pi/1i-M0;
    for m=myinds
      M1(mi,m)=conj(nz(m))*nz(mi).*M1(mi,m);
    end
  end
  M1=2i*imag(M1);
  if iper==1
    ML=sparse(ML);
    M1=sparse(M1);
  end
  
  function M1=CauCinit0(z,wzp,nz,zinter,T16,W16,nseg,N)
% *** Local panelwise evaluation ***
  M1=zeros(N);
  for m=1:N
    M1(:,m)=wzp(m)*conj(nz(m))*nz./(z(m)-z)/pi/1i;
    M1(m,m)=0;
  end
% *** central blocks ***
  for k=1:nseg
    myind=k*16-15:k*16;
    z1=zinter(k);
    z2=zinter(k+1);
    zmid=z2+z1;
    zdif=z2-z1;     
    zs=(2*z(myind)-zmid)/zdif;
    zt=(2*z(myind)-zmid)/zdif;
    ztup=z2-z(myind);
    ztlo=z1-z(myind);
    [TMPM,~]=M1IcompRecC(zt,ztup,ztlo,zs,1);     
    M1(myind,myind)=TMPM/pi/1i;
    for m=myind
      M1(myind,m)=conj(nz(m))*nz(myind).*M1(myind,m);
    end
  end
% *** superdiagonal blocks (targets to the left) ***
  kstart=2;
  for k=kstart:nseg
    myinds=k*16-15:k*16;
    if k==1
      mi=nseg*16-15:nseg*16;       
    else
      mi=(k-1)*16-15:(k-1)*16;
    end
    z1=zinter(k);
    z2=zinter(k+1);
    zmid=z2+z1;
    zdif=z2-z1;     
    zs=(2*z(myinds)-zmid)/zdif;
    zt=(2*z(mi)-zmid)/zdif;          
    ztup=z2-z(mi);
    ztlo=z1-z(mi);
    [TMPM,~]=M1IcompRecC(zt,ztup,ztlo,zs,0);     
    M1(mi,myinds)=TMPM/pi/1i;
    for m=myinds
      M1(mi,m)=conj(nz(m))*nz(mi).*M1(mi,m);
    end
  end
% *** subdiagonal blocks (targets to the right) ***
  kend=nseg-1;
  for k=1:kend
    myinds=k*16-15:k*16;
    if k==nseg
      mi=1*16-15:1*16;       
    else
      mi=(k+1)*16-15:(k+1)*16;
    end
    z1=zinter(k);
    z2=zinter(k+1);
    zmid=z2+z1;
    zdif=z2-z1;     
    zs=(2*z(myinds)-zmid)/zdif;
    zt=(2*z(mi)-zmid)/zdif;
    ztup=z2-z(mi);
    ztlo=z1-z(mi);
    [TMPM,~]=M1IcompRecC(zt,ztup,ztlo,zs,0);     
    M1(mi,myinds)=TMPM/pi/1i;
    for m=myinds
      M1(mi,m)=conj(nz(m))*nz(mi).*M1(mi,m);
    end
  end
  M1=1i*imag(M1);

  function A=diamultL(d,A)
  [~,np]=size(A);
  for k=1:np
    A(:,k)=d.*A(:,k);
  end

  function A=diadivR(A,d)
  [~,np]=size(A);
  for k=1:np
    A(:,k)=A(:,k)/d(k);
  end

  function M1IV=M1IcompRecFS(zt,ztup,ztlo,zs,kern)
% *** zt is target vector in transformed plane ***
  A=ones(16);
  for k=2:16
    A(:,k)=A(:,k-1).*zs;
  end
  p=zeros(16,1);
  c=((1-(-1).^(1:16))./(1:16))';
  p(1)=log(ztup/ztlo);
  for k=2:16
    p(k)=zt*p(k-1)+c(k-1);
  end
  M1IV=real(((A.')\p-kern)/1i);

  function [LGIV,accept]=LGIcompRecFAS(ztgtr,zsctr,ztg,zsc,nz,awzp,cc);
% *** absolute values, transformed ztarg and z ***
  A=ones(16);
  for k=2:16
    A(:,k)=A(:,k-1).*zsctr;
  end
  p=zeros(17,1);
  q=zeros(16,1);
  c=((1-(-1).^(1:16))./(1:16))';
  p(1)=log(1-ztgtr)-log(-1-ztgtr);
  p111=log(1-ztgtr)+log(-1-ztgtr);
  for k=2:17
    p(k)=ztgtr*p(k-1)+c(k-1);
  end
  q(1:2:15)=(p111-p(2:2:16))./(1:2:15)';
  q(2:2:16)=(p(1)-p(3:2:17))./(2:2:16)';
  LGIV=imag(1i*log(cc)+cc*conj(nz).*((A.')\q./awzp));
  LGIV=LGIV-log(abs(zsc-ztg));
  
  function [M1IV,accept]=M1IcompRecC(zt,ztup,ztlo,zs,cauch)
% *** zt is target vector in transformed plane ***
  A=ones(16);
  for k=2:16
    A(:,k)=A(:,k-1).*zs;
  end
  M1IV=zeros(16);
  accept=1:16;
  accept=accept(abs(zt)<2);
  p=zeros(16,1);
  c=((1-(-1).^(1:16))./(1:16))';
  for j=1:16
    p(1)=log(ztup(j)/ztlo(j));
    if cauch==1
      if imag(p(1))<-pi/2
	p(1)=p(1)+pi*i;
      elseif imag(p(1))>pi/2
	p(1)=p(1)-pi*i;       
      else
	error('Konstigt')
      end
    end
    for k=2:16
      p(k)=zt(j)*p(k-1)+c(k-1);
    end
    M1IV(j,:)=p.'/A;
  end
  
  function [LGIV,accept]=LGIcompRecR(trans,mscale,T16);
% *** T is target vector, sources on canonical interval ***
  LGIV=zeros(16);
  A=ones(16);
  for k=2:16
    A(:,k)=A(:,k-1).*T16;
  end
  accept=1:16;
  T=trans+mscale*T16;
  accept=accept(abs(T)<2);
  p=zeros(17,1);
  q=zeros(16,1);
  c=((1-(-1).^(1:16))./(1:16))';
  for j=1:16
    p(1)=log(abs((1-T(j))/(1+T(j))));
    p111=log(abs(1-T(j)^2));
    for k=2:17
      p(k)=T(j)*p(k-1)+c(k-1);
    end
    q(1:2:15)=(p111-p(2:2:16))./(1:2:15)';
    q(2:2:16)=(p(1)-p(3:2:17))./(2:2:16)';
    LGIV(j,:)=q.'/A;
  end
  
  function T=Tinit16
% *** 16-point Gauss-Legendre nodes ***  
  T=zeros(16,1);
  T( 1)=-0.989400934991649932596154173450332627;
  T( 2)=-0.944575023073232576077988415534608345;
  T( 3)=-0.865631202387831743880467897712393132;
  T( 4)=-0.755404408355003033895101194847442268;
  T( 5)=-0.617876244402643748446671764048791019;
  T( 6)=-0.458016777657227386342419442983577574;
  T( 7)=-0.281603550779258913230460501460496106;
  T( 8)=-0.095012509837637440185319335424958063;
  T( 9)= 0.095012509837637440185319335424958063;
  T(10)= 0.281603550779258913230460501460496106;
  T(11)= 0.458016777657227386342419442983577574;
  T(12)= 0.617876244402643748446671764048791019;
  T(13)= 0.755404408355003033895101194847442268;
  T(14)= 0.865631202387831743880467897712393132;
  T(15)= 0.944575023073232576077988415534608345;
  T(16)= 0.989400934991649932596154173450332627;

  function W=Winit16
% *** 16-point Gauss-Legendre weights ***  
  W=zeros(16,1); 
  W( 1)= 0.027152459411754094851780572456018104;
  W( 2)= 0.062253523938647892862843836994377694;
  W( 3)= 0.095158511682492784809925107602246226;
  W( 4)= 0.124628971255533872052476282192016420;
  W( 5)= 0.149595988816576732081501730547478549;
  W( 6)= 0.169156519395002538189312079030359962;
  W( 7)= 0.182603415044923588866763667969219939;
  W( 8)= 0.189450610455068496285396723208283105;
  W( 9)= 0.189450610455068496285396723208283105;
  W(10)= 0.182603415044923588866763667969219939;
  W(11)= 0.169156519395002538189312079030359962;
  W(12)= 0.149595988816576732081501730547478549;
  W(13)= 0.124628971255533872052476282192016420;
  W(14)= 0.095158511682492784809925107602246226;
  W(15)= 0.062253523938647892862843836994377694;
  W(16)= 0.027152459411754094851780572456018104;

  function [IP,IPW]=IPWinit(T16,W16)
  A=ones(16);
  AA=ones(32,16);
  T32=[(T16-1)/2;(T16+1)/2];
  W32=[W16;W16]/2;
  for k=2:16
    A(:,k)=A(:,k-1).*T16;
    AA(:,k)=AA(:,k-1).*T32;         
  end
  IP=AA/A;
  IPW=diadivR(diamultL(W32,IP),W16);

  function Pbc=Pbcinit(IP);
  Pbc=zeros(64,32);
  Pbc( 1:32, 1:16)=IP;
  Pbc(33:64,17:32)=flipud(fliplr(IP));

  function PbcBig=PbcBiginit(IP);  
  Pbc=Pbcinit(IP);
  PbcBig=zeros(192,96);
  PbcBig(  1: 64, 1:32)=Pbc;
  PbcBig( 65:128,33:64)=Pbc;
  PbcBig(129:192,65:96)=Pbc;
  
  function Pbc=PbcinitBig(IP);
  Pbc=zeros(96,64);
  Pbc( 1:16, 1:16)=eye(16);
  Pbc(17:48,17:32)=IP;
  Pbc(49:96,33:64)=flipud(fliplr(Pbc(1:48,1:32)));
  
  function PbcBB=PbcBiginitBig(IP);  
  Pbc=PbcinitBig(IP);
  PbcBB=zeros(288,192);
  PbcBB(  1: 96,  1: 64)=Pbc;
  PbcBB( 97:192, 65:128)=Pbc;
  PbcBB(193:288,129:192)=Pbc;