function demo11b
% ********************************************************
% *** Last changed 2012-11-10 by Johan Helsing           *
% *** Exterior Dirichlet problem for Helmholtz           *
% *** RCIP method and combined field integral equation   *
% *** Field evaluation                                   *
% ********************************************************
  close all
  format long
  format compact
%
% *** user specified quantities ****************************
  itmax=100         % maximum number of GMRES iterations   *
  theta=pi/2        % corner opening angle                 *
  nsub=52;          % number of levels in recursion        *
  premature=48;     % interruption level in reconstruction *
  zsource=0.3+0.1i; % source point                         *
  omega=10          % wave number                          *
  npan=46           % number of panels on coarse mesh      *
  dlim=0.7          % 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);
  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);
%
% *** 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 K_coa^{\circ} matrix is set up ***
  disp('Setup of K_coa^{\circ} starts')
  LogCcoa=LogCinit(sinterdiff,T16,W16,npan,npcoa,1);
  Kopercoa=Koperinit(LogCcoa,zcoa,zpcoa,zppcoa,wzpcoa,wcoa,omega,npcoa);
  Sopercoa=Soperinit(LogCcoa,zcoa,zpcoa,awzpcoa,sinterdiff,omega,npcoa);
  clear LogCcoa
  Kcirccoa=Kopercoa-0.5i*omega*Sopercoa;
  clear Kopercoa Sopercoa
  starind=[npcoa-31:npcoa 1:32];
  Kcirccoa(starind,starind)=zeros(64);
%
% *** Recursion for the R matrix ***
  disp('Recursion starts')
  R=R0;
  Rstor=zeros(64,64,nsub);
  Klocstor=zeros(96,96,nsub);
  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;
    wzploc=wloc.*zploc;
    awzploc=abs(wzploc);
    sidiloc=[1;0.5;0.5;0.5;0.5;1]/denom;
    Kloc=Koperinit(LogCloc,zloc,zploc,zpploc,wzploc,wloc,omega,96);
    Sloc=Soperinit(LogCloc,zloc,zploc,awzploc,sidiloc,omega,96);
    Klocstor(:,:,level)=Kloc-0.5i*omega*Sloc;
    MAT=eye(96)+Kloc-0.5i*omega*Sloc;
    Rstor(:,:,level)=R;
    R=SchurBana(PWbc,Pbc,MAT,R);
  end
  Rcoa=speye(npcoa);
  Rcoa(starind,starind)=R;
%
% *** Solving main linear system ***
  disp('Solving main linear system starts')
  g=besselh(0,omega*abs(zsource-zcoa));     
  rhs=2*g;
  [rhotildecoa,it]=myGMRESR(Kcirccoa,Rcoa,rhs,npcoa,itmax,eps);
  clear Kcirccoa
  GMRES_iter_RCIP=it
%
% *** Post processing ***
  disp('Post processing starts')
%
% *** rhofin locally ***
  nse=nsub-premature+1;
  rhotilde=rhotildecoa(starind);
  rhofinloc=zeros(16*(4+2*nse),1);
  Pbc=PbcinitBig(IP);
  for k=nsub:-1:premature
    Kcirc=Klocstor(:,:,k);
    Kcirc(17:80,17:80)=zeros(64);
    MAT=eye(96)+Kcirc;
    MAT(17:80,17:80)=inv(Rstor(:,:,k));
    myind1L=(nsub-k)*16+(1:16);
    myind2L=(nsub+k-2*premature+5)*16+(1:16);
    tmp=Pbc*rhotilde;
    rhotilde=tmp-Kcirc*(MAT\tmp);
    rhofinloc([myind1L myind2L])=rhotilde([1:16 81:96]);
    rhotilde=rhotilde(17:80);
  end
  rhofinloc(nse*16+1:(4+nse)*16)=Rstor(:,:,premature)*rhotilde;
  clear Klocstor Rstor
%
% *** rhofin globally ***
  rhofin=zeros(16*(npan+2*nse),1);
  rhofin(1:16*(2+nse))=rhofinloc(16*(2+nse)+1:16*(4+2*nse));
  rhofin(16*(2+nse)+1:16*(2+nse)+(npan-4)*16)=rhotildecoa(33:(npan-2)*16);
  rhofin(16*(2+nse)+(npan-4)*16+1:16*(npan+2*nse))=rhofinloc(1:16*(2+nse));
%
  [zfin,wzpfin,nzfin,zinterfin,npanfin]=finegrid(nse,theta,npan,W16,T16);  
%
% ***  
  rhohatcoa=Rcoa*rhotildecoa;
  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=rhohatcoa.'*(Ktarg(zg(k),zcoa,wzpcoa,omega) ...
		     -0.5i*omega*Starg(zg(k),zcoa,awzpcoa,omega));
    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);
      rhohat=rhohatcoa(myind);
      d=min(abs(zsc-ztarg))/panlen(kk);
      if d>dlim
	Unum=Unum+rhohat.'*(Ktarg(ztarg,zsc,wzp,omega) ...
			    -0.5i*omega*Starg(ztarg,zsc,awzp,omega));
      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+rhohat.'*(KtargClose(LogC,CauC,ztarg,zsc,wzp,omega) ...
			    -0.5i*omega*StargClose(LogC,ztarg,zsc,awzp,omega));
      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);
      rhohat=rhofin(myind);
      d=min(abs(zsc-ztarg))/sum(awzp);
      if d>dlim
	Unum=Unum+rhohat.'*(Ktarg(ztarg,zsc,wzp,omega) ...
			    -0.5i*omega*Starg(ztarg,zsc,awzp,omega));
      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+rhohat.'*(KtargClose(LogC,CauC,ztarg,zsc,wzp,omega) ...
			    -0.5i*omega*StargClose(LogC,ztarg,zsc,awzp,omega));
      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('Dirichlet: 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 [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=M1Rinit(z,zp,zpp,w,wzp,N)
% *** double layer potential ***   
  M1=zeros(N);
  for m=1:N
    M1(:,m)=imag(wzp(m)./(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
    S(:,m)=besselh(0,omega*abs(z-z(m)))*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 K=Koperinit(LogC,z,zp,zpp,wzp,w,omega,N)
  K=zeros(N);
  for m=1:N
    tmp=omega*abs(z-z(m));
    K(:,m)=tmp.*besselh(1,tmp).*imag(wzp(m)./(z-z(m)));
    K(m,m)=1i/pi*imag(zpp(m)/zp(m))*w(m);  
  end
  myind=find(LogC);
  K(myind)=K(myind)+2i/pi*(real(K(myind)).*LogC(myind));
  K=1i/2*K;

  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 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];
  wzploc=wloc.*zploc;
  Kloc=-M1Rinit(zloc,zploc,zpploc,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 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=diadivR(A,d)
  [~,np]=size(A);
  for k=1:np
    A(:,k)=A(:,k)/d(k);
  end

  function A=diamultL(d,A)
  [~,np]=size(A);
  for k=1:np
    A(:,k)=d.*A(:,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 [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 [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=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 Pbc=Pbcinit(IP);
  Pbc=zeros(64,32);
  Pbc( 1:32, 1:16)=IP;
  Pbc(33:64,17:32)=flipud(fliplr(IP));
  
  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;