function demo4
  close all
  format long
  format compact
%
% *** user specified quantities **********
  lambda=0.999   % parameter lambda
  theta=pi/2     % corner opening angle
  npan=10;       % number of coarse panels
  evec=1;        % a unit vector 
% ****************************************
  T16=Tinit16;
  W16=Winit16;
  [IP,IPW]=IPWinit(T16,W16);
  T16p1=T16+1;
  for nsub=1:100
    [rhofin_direct,wzpfin]=direct(nsub,lambda,theta,npan,evec,W16,T16);
    rhofin_RCIP=RCIP(nsub,lambda,theta,npan,evec,IP,IPW,W16,T16,T16p1);
    myerr=myL2diff(rhofin_RCIP,rhofin_direct,abs(wzpfin));
%
    set(0,'DefaultAxesFontSize',16)
    set(0,'DefaultAxesLineWidth',1) 
    figure(1)
    semilogy(nsub,myerr,'*')
    axis([0 100 1e-16 1e0])
    title('Difference between direct and reconstructed \rho_{fin}','FontSize',16)
    xlabel('Number of subdivisions {\itn}_{sub}','FontSize',16)
    ylabel('Relative L^2 difference','FontSize',16)  
    axis square
    grid on  
    hold on
  end

  function myerr=myL2diff(f,g,ds)
  myerr=sqrt((ds'*(f-g).^2)/(ds'*f.^2));

  function rhofin=RCIP(nsub,lambda,theta,npan,evec,IP,IPW,W16,T16,T16p1)
%
  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)));
%
  PWbc=zeros(96,64);
  PWbc( 1:16, 1:16)=eye(16);
  PWbc(17:48,17:32)=IPW;
  PWbc(49:96,33:64)=flipud(fliplr(PWbc(1:48,1:32)));
  check_PwTP=norm(PWbc'*Pbc-eye(64))
%
% *** 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) ;
  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)));
  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;
  arclength=sum(abs(wzpcoa))
  area=0.5*imag(conj(zcoa).'*wzpcoa)
%
% *** The K_coa^{\circ} matrix is set up ***
  Kcirccoa=M1Ainit(zcoa,zpcoa,zppcoa,nzcoa,wcoa,wzpcoa,npcoa);
  starind=[npcoa-31:npcoa 1:32];
  Kcirccoa(starind,starind)=zeros(64);
%
% *** recursion for the R matrix ***
  Rstor=zeros(64,64,nsub);
  Klocstor=zeros(96,96,nsub);
  for level=1:nsub
    sloc=[T16p1/4;T16/4+0.75;T16/2+1.5]/2^(nsub-level)/npan;
    wloc=[W16/4;W16/4;W16/2]/2^(nsub-level)/npan;
    zloc=zfunc(sloc,theta);
    zploc=zpfunc(sloc,theta);
    zpploc=zppfunc(sloc,theta);
    zloc=[conj(flipud(zloc));zloc];
    zploc=[-conj(flipud(zploc));zploc];
    zpploc=[conj(flipud(zpploc));zpploc];
    wloc=[flipud(wloc);wloc];
    nzloc=-1i*zploc./abs(zploc);
    wzploc=wloc.*zploc;
    Kloc=M1Ainit(zloc,zploc,zpploc,nzloc,wloc,wzploc,96);
    Klocstor(:,:,level)=Kloc;  
    MAT=eye(96)+lambda*Kloc;
    if level==1
      R=inv(MAT(17:80,17:80));
    end
    Rstor(:,:,level)=R;  
    R=SchurBana(PWbc,Pbc,MAT,R);
  end
  Rcoa=eye(npcoa);
  Rcoa(starind,starind)=R;
%
% *** Solving main linear system ***
  g=2*real(conj(evec)*nzcoa);     
  [rhotildecoa,it]=myGMRESR(lambda*Kcirccoa,Rcoa,lambda*g, ...
			    abs(wzpcoa),npcoa,100,eps);
  GMRES_iter_RCIP=it
%
% *** Post processing ***
  rhohatcoa=Rcoa*rhotildecoa;
  zeta=real(conj(evec)*zcoa.*abs(wzpcoa));
  qcoa=rhohatcoa'*zeta
%
% *** Reconstruction of rhofin on Gamma^{\star} ***
  rhotilde=rhotildecoa(starind);
  rhofinloc=zeros(16*(4+2*nsub),1);
  for k=nsub:-1:1
    lmbKcirc=lambda*Klocstor(:,:,k);
    lmbKcirc(17:80,17:80)=zeros(64);
    MAT=eye(96)+lmbKcirc;
    MAT(17:80,17:80)=inv(Rstor(:,:,k));
    myind1L=(nsub-k)*16+(1:16);
    myind2L=(nsub+k+3)*16+(1:16);
    tmp=Pbc*rhotilde;
    rhotilde=tmp-lmbKcirc*(MAT\tmp);
    rhofinloc([myind1L myind2L])=rhotilde([1:16 81:96]);
    rhotilde=rhotilde(17:80);
  end
  rhofinloc(nsub*16+1:(4+nsub)*16)=Rstor(:,:,1)*rhotilde;
%
% *** rhofin globally ***
  rhofin=zeros(16*(npan+2*nsub),1);
  rhofin(1:16*(2+nsub))=rhofinloc(16*(2+nsub)+1:16*(4+2*nsub));
  rhofin(16*(2+nsub)+1:16*(2+nsub)+(npan-4)*16)=rhotildecoa(33:(npan-2)*16);
  rhofin(16*(2+nsub)+(npan-4)*16+1:16*(npan+2*nsub))=rhofinloc(1:16*(2+nsub));

  function [rhofin,wzpfin]=direct(nsub,lambda,theta,npan,evec,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) ;
  zpfin=zpfunc(sfin,theta);
  zppfin=zppfunc(sfin,theta);
% *** some extra presicion gained from symmetry ***
  zfin(npfin/2+1:npfin)=conj(flipud(zfin(1:npfin/2)));
  zpfin(npfin/2+1:npfin)=-conj(flipud(zpfin(1:npfin/2)));
  zppfin(npfin/2+1:npfin)=conj(flipud(zppfin(1:npfin/2)));
% *************************************************
  nzfin=-1i*zpfin./abs(zpfin);
  wzpfin=wfin.*zpfin;
%
% *** The K_fin matrix is set up ***
  Kfin=M1Ainit(zfin,zpfin,zppfin,nzfin,wfin,wzpfin,npfin);
%
% *** Solving main linear system ***
  g=2*real(conj(evec)*nzfin);     
  [rhofin,it]=myGMRES(lambda*Kfin,lambda*g,abs(wzpfin),npfin,100,eps);
  GMRES_iter_direct=it

  function A=SchurBana(PW,P,K,A)
  starL=17:80;
  circL=[1:16 81:96];
  starS=17:48;
  circS=[1:16 49:64];
  VA=K(circL,starL)*A;
  PTA=PW(starL,starS)'*A;
  PTAU=PTA*K(starL,circL);
  DVAUI=inv(K(circL,circL)-VA*K(starL,circL));
  DVAUIVAP=DVAUI*(VA*P(starL,starS));
  A(starS,starS)=PTA*P(starL,starS)+PTAU*DVAUIVAP;
  A(circS,circS)=DVAUI;
  A(circS,starS)=-DVAUIVAP;
  A(starS,circS)=-PTAU*DVAUI;
  
  function [x,it]=myGMRES(A,b,ds,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;                                   
    w=A*V(:,it)+V(:,it).'*ds;
    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);
      true_residual=norm(x+A*x+x.'*ds-b)/bnrm2
      break;
    end
  end

  function [x,it]=myGMRESR(A,R,b,ds,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;                                   
    tmp=R*V(:,it);
    w=A*tmp+tmp'*ds;
    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);
      tmp=R*x;
      true_residual=norm(x+A*tmp+tmp'*ds-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 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 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 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=IP.*(W32*(1./W16)');