function demo12
% ********************************************************
% *** Last changed 2012-10-25 by Johan Helsing           *
% *** Exterior Neumann problem for Helmholtz             *
% *** RCIP method and analogy with approach for Laplace  *
% *** Initial guess from fixed-point iteration           *
% ********************************************************
  close all
  format long
  format compact
%
% *** user specified quantities **************************
  itmax=800         % maximum number of GMRES iterations *
  theta=pi/2        % corner opening angle               *
  nsub=52;          % number of steps in recursion       *
  zsource=0.3+0.1i; % source point                       *
  ztarg=-0.1+0.2i;  % target point                       *
% ********************************************************
  T16=Tinit16;
  W16=Winit16;
  [IP,IPW]=IPWinit(T16,W16);
  Pbc=Pbcinit(IP);
  PWbc=Pbcinit(IPW);
  R0=initializeR(W16,T16,theta,PWbc,Pbc);
  sidiloc=[1;0.5;0.5;0.5;0.5;1];
  LogCloc=LogCinit(sidiloc,T16,W16,6,96,0);
%
  omegavec=logspace(0,3,1000);
  for omega=omegavec
    omega
    npan=round(0.6*omega+18)
% *** 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;
    awzpcoa=abs(wzpcoa);
    arclength=sum(awzpcoa)
    area=0.5*imag(conj(zcoa).'*wzpcoa)
%
% *** The K_coa^{\circ} matrix is set up ***
    disp('Setup of K_coa^{\circ} starts')
    LogCcoa=LogCinit(sinterdiff,T16,W16,npan,npcoa,1);
    Kcirccoa=-Kadjoperinit(LogCcoa,zcoa,zpcoa,zppcoa,awzpcoa,wcoa, ...
			 nzcoa,omega,npcoa);
    starind=[npcoa-31:npcoa 1:32];
    Kcirccoa(starind,starind)=zeros(64);
%
% *** recursion for the R matrix ***
    disp('Recursion starts')
    R=R0;
    for level=1:nsub
      denom=2^(nsub-level)*npan;
      sloc=[T16/4+0.25;T16/4+0.75;T16/2+1.5]/denom;
      zloc=zfuncloc(sloc,theta);
      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);
      awzploc=wloc.*abs(zploc);
      Kloc=-Kadjoperinit(LogCloc,zloc,zploc,zpploc,awzploc,wloc, ...
			 nzloc,omega,96);
      MAT=eye(96)+Kloc;
      R=SchurBana(PWbc,Pbc,MAT,R);
    end
    Rcoa=speye(npcoa);
    Rcoa(starind,starind)=R;
%
% *** 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]=myGMRESR(Kcirccoa,Rcoa,rhs,npcoa,itmax,eps);
    GMRES_iter_RCIP=it
%
% *** Post processing ***
    rhohatcoa=Rcoa*rhotildecoa;
    ufield=rhohatcoa.'*Starg(ztarg,zcoa,awzpcoa,omega)
    ufiref=besselh(0,omega*abs(zsource-ztarg))
    myerr=abs(ufiref-ufield)
    myerr(myerr<eps)=eps;
    myplot(omega,myerr,GMRES_iter_RCIP)
  end    

  function myplot(omega,myerr,GMRES_iter_RCIP)    
  set(0,'DefaultAxesFontSize',16)
  set(0,'DefaultAxesLineWidth',1) 
  figure(1)
  loglog(omega,GMRES_iter_RCIP,'*')
  axis([1 1e3 1 1e3])
  title('GMRES iterations for full convergence','FontSize',16)
  xlabel('\omega','FontSize',18)
  ylabel('number of iterations ','FontSize',16)  
  grid on  
  hold on
  drawnow
  figure(2)
  loglog(omega,myerr,'*')
  axis([1 1e3 1e-16 1e0])
  title('Analogy with standard approach for Laplace','FontSize',16)
  xlabel('\omega','FontSize',18)
  ylabel('estimated error in {\itU}({\itx})','FontSize',16)  
  grid on  
  hold on
  drawnow
  
  function R=initializeR(W16,T16,theta,PWbc,Pbc)
  myerr=1;
  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;
  Kloc=-M1Ainit(zloc,zploc,zpploc,nzloc,wloc,wzploc,96);
  MAT=eye(96)+Kloc;
  R=inv(MAT(17:80,17:80));
  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 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 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 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'*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]=myGMRESR(A,R,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;                                   
    w=A*(R*V(:,it));
    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*(R*x)-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 M1=LogCinit(pinterdiff,T16,W16,nseg,N,iper)
% *** Corrections to Logarithmic potential log(|tau-z|) ***   
% *** block-tri-diagonal output ***
% *** iper=0,1 (0 is open arc, 1 is closed contour) ***
  [TMP,~]=LGIcompRecR(0,1,T16);
  TMP=diadivR(TMP,W16);
  M1=zeros(N);
  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);
      M1(myind,m)=-log(abs(T16(nj)-T16));
      M1(m,m)=0;
    end
    M1(myind,myind)=M1(myind,myind)+TMP;
  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);          
    [TMP,accept]=LGIcompRecR(-1-alpha,alpha,T16);
    mi=mi(accept);
    for nj=1:16
      M1(mi,myinds(nj))=-log(abs(T16(nj)+1+(1-T16(accept))*alpha));       
    end
    TMP=TMP(accept,:);
    M1(mi,myinds)=M1(mi,myinds)+diadivR(TMP,W16);
  end
% *** subdiagonal blocks (targets to the right) ***
  for k=1:kend
    myinds=k*16-15:k*16;
    kp1=mod(k,nseg)+1;
    mi=kp1*16-15:kp1*16;       
    alpha=pinterdiff(kp1)/pinterdiff(k);               
    [TMP,accept]=LGIcompRecR(1+alpha,alpha,T16);
    mi=mi(accept);
    for nj=1:16
      M1(mi,myinds(nj))=-log(abs(T16(nj)-1-(T16(accept)+1)*alpha));   
    end
    TMP=TMP(accept,:);
    M1(mi,myinds)=M1(mi,myinds)+diadivR(TMP,W16);
  end
  if iper==1
    M1=sparse(M1);
  end
  
  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 [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));