% TESTSR_RIMKLTPM  Test all four routines sr_RMKLTp.m, sr_RMKLTm.m, 
%    sr_IMKLTp.m, and sr_IMKLTm.m
%
f1='%8.0f %22.14e\n';
N=40;
a=.5;
for b=[.01 .1 .5 1 3 5 10 15]
%  dig=70
  dig=75
%  dig=80
%  ab1=sr_RMKLTp(dig,N,a,b);
%  ab2=sr_RMKLTm(dig,N,a,b);
  ab3=sr_IMKLTp(dig,N,a,b);
  ab4=sr_IMKLTm(dig,N,a,b);
%  ab1=vpa(ab1,16);
%  ab1=double(ab1);
%  ab2=vpa(ab2,16);
%  ab2=double(ab2);
   ab3=vpa(ab3,16);
   ab3=double(ab3);
   ab4=vpa(ab4,16);
   ab4=double(ab4);
%
% Example 4.3 
%
  abjp=r_jacobi01(N,0,a);
  abjm=r_jacobi01(N,0,-a);
  lnGp=mfun('lnGAMMA',a+b*i);
  gp=imag(lnGp);
  Gp=exp(real(lnGp));
  lnGm=mfun('lnGAMMA',-a-b*i);
  gm=imag(lnGm);
  Gm=exp(real(lnGm));
  for n=3:3:15
%    xw1=gauss(n,ab1);
    xw3=gauss(n,ab3);
%    xw2=gauss(n,ab2);
    xw4=gauss(n,ab4);
%    x1=xw1(:,1); w1=xw1(:,2);
    x3=xw3(:,1); w3=xw3(:,2);
%    x2=xw2(:,1); w2=xw2(:,2);
    x4=xw4(:,1); w4=xw4(:,2);
    xwjp=gauss(n,abjp);
    xwjm=gauss(n,abjm);
    xjp=xwjp(:,1); wjp=xwjp(:,2);
    xjm=xwjm(:,1); wjm=xwjm(:,2);
%    I1=sum(w1.*exp(-2*x1));
%    I2=sum(w2.*exp(-2*x2));
    I1=sum(w3.*exp(-2*x3));
    I2=sum(w4.*exp(-2*x4));
    Ijp=sum(wjp.*exp(-2*xjp));
    Ijm=sum(wjm.*exp(-2*xjm));
%    RI3=Gp*I1+Gm*I2-Gp*Ijm-Gm*Ijp;
    II3=Gp*I1-Gm*I2-Gp*Ijm+Gm*Ijp;
%    fprintf(f1,n,RI3)
    fprintf(f1,n,II3)
  end
end