function demo13c
% **********************************************************
% *** Last changed 2012-11-09 by Johan Helsing             *
% *** Exterior Neumann problem for Helmholtz               *
% *** RCIP mehtod and regularized combined field equation  *
% *** Local regularization used for the Cauchy operator    *
% **********************************************************
  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 steps in recursion
  zsource=0.3+0.1i; % source point
  ztarg=-0.1+0.2i;  % target point
% ***************************************************
  T16=Tinit16;
  W16=Winit16;
  diff16=diffinit(T16);
  [IP,IPW]=IPWinit(T16,W16);
  PbcBig=PbcBiginit(IP);  
  PWbcBig=PbcBiginit(IPW);    
  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+48)
% *** 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)
    disp('Setup of coarse matrices starts')
%
% *** The coarse matrices are set up ***
    LogCcoa=LogCinit(sinterdiff,T16,W16,npan,npcoa,1);
    CauCcoa=CauCinit(zcoa,zinter,wzpcoa,nzcoa,wcoa,sinterdiff,diff16,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);
    disp('Recursion starts')
%
% *** recursion for the R matrix ***
    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;
      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);
      wzploc=wloc.*zploc;
      awzploc=abs(wzploc);
      sidiloc=[1;0.5;0.5;0.5;0.5;1]/denom;
      if level==nsub
	sidilocX=[1;1;0.5;0.5;0.5;0.5;1;1]/denom;	
	slocX=[T16/4+0.25;T16/4+0.75;T16/2+1.5;T16/2+2.5]/denom;      
	sinterlocX=[0;0.5;1;2;3]/denom;      
	wlocX=[W16/2;W16/2;W16/4;W16/4;W16/4;W16/4;W16/2;W16/2]/denom;
      else
	sidilocX=[2;1;0.5;0.5;0.5;0.5;1;2]/denom;		
	slocX=[T16/4+0.25;T16/4+0.75;T16/2+1.5;T16+3]/denom;            
	sinterlocX=[0;0.5;1;2;4]/denom;            
	wlocX=[W16;W16/2;W16/4;W16/4;W16/4;W16/4;W16/2;W16]/denom;
      end
      zlocX=zfuncloc(slocX,theta);
      zinterlocX=zfuncloc(sinterlocX,theta);
      zinterlocX(5)=[];
      zplocX=zpfuncloc(slocX,theta);
      nzlocX=-1i*zplocX./abs(zplocX);      
      wzplocX=wlocX.*zplocX;
      CauCloc=CauCinit(zlocX,zinterlocX,wzplocX,nzlocX,wlocX,sidilocX, ...
		       diff16,8,128,0);
      CauCloc=CauCloc(17:112,17:112);
      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;
      MAT=eye(288)+ABCloc;
      if level==1
	R=inv(MAT(starL,starL));
      end
      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 ***
    ufiref=besselh(0,omega*abs(zsource-ztarg))
    rho1hatcoa=R1G*rhotildecoa;
    B2rho1hatcoa=B2circcoa*rho1hatcoa;
    C2rho1hatcoa=C2circcoa*rho1hatcoa;
    rho2hatcoa=R2G*rhotildecoa+R5214G*B2rho1hatcoa+R8217G*C2rho1hatcoa;
    ufield=Starg(ztarg,zcoa,awzpcoa,omega).'*rho1hatcoa+ ...
	   1i*(Ktarg(ztarg,zcoa,wzpcoa,omega).'*rho2hatcoa)
    myerr=abs(ufiref-ufield)
    myerr(myerr<eps)=eps;
    myplot(omega,myerr,GMRES_iter_RCIP)    
    clear Acirccoa B1circcoa B2circcoa C1circcoa C2circcoa
  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 1e2])
  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('Regularized field integral equation','FontSize',16)
  xlabel('\omega','FontSize',18)
  ylabel('estimated error in {\itU}({\itx})','FontSize',16)  
  grid on  
  hold on
  drawnow

  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 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 M1=CauCinit(z,zinter,wzp,nz,w,pinterdiff,diff16,nseg,np,iper)
% *** Corrections to Cauchy operator with Local regulariaztion ***
  M1=zeros(np);
  if iper==0
    kstart=2;
    kend=nseg-1;
  else
    kstart=1;
    kend=nseg;
  end
  for k=kstart:kend
    myindt=(k-1)*16+1:k*16;
    if k==1
      z1=zinter(nseg);
      z2=zinter(3);
      myinds=[(nseg-1)*16+1:nseg*16 1:32];
    elseif k==nseg
      z1=zinter(k-1);
      z2=zinter(2);
      myinds=[(nseg-2)*16+1:nseg*16 1:16];    
    else
      z1=zinter(k-1);
      z2=zinter(k+2);
      myinds=[(k-2)*16+1:(k+1)*16];  
    end
    M0=zeros(16,48);
    for j=1:48
      m=myinds(j);
      M0(:,j)=wzp(m)./(z(m)-z(myindt))/pi/1i;
    end
    for j=1:16
      M0(j,16+j)=0;
      m=myindt(j);
      logterm=log((z2-z(m))/(z1-z(m)))/pi/1i;
      if real(logterm)<-0.25
	logterm=logterm+1;
      elseif real(logterm)>0.25
	logterm=logterm-1;
      else
	logterm
	error('Strange')
      end
      M1(m,m)=logterm-sum(M0(j,:));
    end
  end
  for k=1:nseg
    pdif=pinterdiff(k)/2;
    myind=k*16-15:k*16;
    M1(myind,myind)=M1(myind,myind)+diamultL(w(myind)/pdif/pi/1i,diff16);
    for j=1:16
      M1(myind,myind(j))=conj(nz(myind(j)))*nz(myind).*M1(myind,myind(j));
    end
    M1(myind,myind)=2i*imag(M1(myind,myind));
  end
  if iper==1
    M1=sparse(M1);
  end

  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 [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 diffmat=diffinit(TA)
  A=ones(16);
  B=zeros(16);
  for k=2:16
    A(:,k)=A(:,k-1).*TA;
    B(:,k)=(k-1)*TA.^(k-2);  
  end
  diffmat=B/A;
  
  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;