function [mmom,canmin]=smmom_Enufin(dig,N,nu,c)
%   MMOM=SMMOM_ENUFIN(DIG,N,NU,C) computes in DIG-digit
%   arithmetic the first 2N modified moments, relative
%   to monic shifted Jacobi01 polynomials with Jacobi
%   parameters alpha=0, beta=nu-1, transformed from  
%   the interval [0,1] to [0,C] of the subrange exponential 
%   integral weight function w(x)=E_nu(x) on [0,C].     

tic
disp('Be sure to load the correct file ab_powexp190_nu_c.txt') 
disp('on line 20 of the routine') 
digits(dig); nu=vpa(nu); c=vpa(c); s0=vpa(0); 
s1=vpa(1);
mmom=vpa(zeros(2*N,1)); can=vpa(zeros(N)); canmin=zeros(2*N,1);
P=vpa(zeros(2*N,1)); Q=vpa(zeros(2*N,1));
abl=sr_laguerre(dig,N); xwl=sgauss(dig,N,abl); 
abj=sr_jacobi01(dig,2*N-1,s0,nu-s1); xwj=sgauss(dig,N,abj);
abjc(:,1)=c*abj(:,1); abjc(1,2)=(c^nu)*abj(1,2);
abjc(2:2*N-1,2)=c^2*abj(2:2*N-1,2);
abstar=loadvpa('ab_powexp190_1_5.txt',9,190,40,2); 
abstar=vpa(abstar,dig); xwstar=sgauss(dig,N,abstar); 
for i=1:N
  for j=1:N
    x(i,j)=xwl(i,1)*xwj(j,1); y(i,j)=c+xwl(i,1)*xwstar(j,1);
    P(1,i,j)=1; Q(1,i,j)=1;
    P(2,i,j)=x(i,j)-abjc(1,1); Q(2,i,j)=y(i,j)-abjc(1,1);
    for k=2:2*N-1
      P(k+1,i,j)=(x(i,j)-abjc(k,1))*P(k,i,j)-abjc(k,2)*P(k-1,i,j);
      Q(k+1,i,j)=(y(i,j)-abjc(k,1))*Q(k,i,j)-abjc(k,2)*Q(k-1,i,j);
    end 
  end
end
for k=1:2*N
  u1=s0;
  for i=1:N
    u2=s0;
    for j=1:N
        u2=u2+xwj(j,2)*P(k,i,j)-xwstar(j,2)*Q(k,i,j);
    end
    u1=u1+xwl(i,2)*u2;
    can(i,j)=abs((xwj(j,2)*P(k,i,j)-xwstar(j,2)*Q(k,i,j))/ ...
      (xwj(j,2)*P(k,i,j)));
  end
  canmin(k)=min(min(can)); canmin(k)=double(canmin(k));
  mmom(k)=u1;
end
toc
