% Script: sincfcn
% Programmer: t. shores
% Last Rev: 7/24/01
% Description:
% a package of functions written to evaluate sincs, S(k,h)(z)'s and 
% truncated sinc series C(v,h)(x), as well as their derivatives, reliably. 
% Observe that Octave already has a built-in sinc function.
% Matlab users will have to create separate function files for each
% function defined below and keep them in a common directory, and also
% delete any occurrence of endfunction.

% Matlab users: you will have to unpack each function in this file into 
% a separate function file.  Octave users: simply run the script.
%
% (for Octave) not a function file:
1;

%

function retval = arrayf(fname,arr)
% usage: arrayf(fname,arr)
% description: this returns the matrix resulting from replacing
% each coordinate of arr by its evaluation by the function with
% name fname.

% local variables
% ii,jj: index variables
% nr,nc: number of rows, cols

[nr,nc] = size(arr);
retval = zeros(nr,nc);
for ii = 1:nr
  for jj = 1:nc
    retval(ii,jj) = feval(fname,arr(ii,jj));
  end;
end;

endfunction

%

function idn = Idn(ndx,n)
%usage: Idn(ndx,n)
% description: Creates the necessary sinc identity matrices I^(i), with exponent i=-1,0,1,2.  Be warned that this proc is of a split personality.  For i=0,1,2, the matrices of Lund and Bowers [LB, p.106] are used, while for i=-1, Stenger [St,p.176] is used...-1 is in working order, but very slow.  One could improve this by doing table lookup for the necessary values of Si(k).

% Programmer: t. shores	    Last Rev: 4/6/96

% local variables:
% idn: temporary storage for returned matrix
% k,m: loop indices
% inputs:
% ndx: exponent of requested identity
% n: dimension of the requested identity
% outputs:
% return floating point value of the identity.

idn = eye(n);
c = zeros(n,1);
if (ndx == 1)
  for k = 2:n 
     c(k,1) = (-1)^k/(k-1);
  end
  idn = toeplitz(c,-c);
elseif (ndx == -1)
% use sinc integration on intervals [k,k+1] here
% assuming value 0 at endpoints, so do c(2,1) separately
  c(2,1) = 0.5894898722360839;
  M = 60;
  h = pi/sqrt(2*M);
  nodes = h*[-M:M];
  enodes = exp(nodes);
  for k = 3:n
    a = k-2;
    b = k-1;
    inodes = (b*enodes + a) ./ (enodes + 1);
    c(k,1) = c(k-1,1)+h*sum((sin(pi*inodes) ./ (pi*inodes)) .* (inodes - a) .* (b - inodes));
  end
  idn = toeplitz(c,-c);
  idn = 0.5*ones(n,n)+idn;
elseif (ndx == 2)
  c(1,1) = -pi^2/3;
  for k = 2:n
      c(k,1) = -2*(-1)^(k-1)/(k-1)^2;
  end
  idn = toeplitz(c,c);
end

endfunction

%

function retval = sincp(x)
% usage: y = sincp(x)
% description:
% this computes sinc'(x) reliably, x a matrix.

% local variables:
% qq,qq: row,column size of x
% kk: loop variable
% epsl: to 4th power is approximate 0..adjust with precision
% xx: temporary float
% inputs:
% x: argument
% outputs: sinc'(x)

epsl = 0.0001;
[pp,qq] = size(x);
retval = zeros([pp,qq]);
for kk = 1:pp
  for ll = 1:qq
    xx = pi*x(kk,ll);
    if abs(xx) > epsl
      retval(kk,ll) = pi*(xx*cos(xx)-sin(xx))/(xx*xx);
    else
      retval(kk,ll) = pi*xx*(xx*xx/10.0-1.0)/3.0;
    end
  end
end

endfunction

%

function retval = sincpp(x)
% usage: y = sincpp(x)
%description:
% this computes sinc''(x) reliably, x a vector.

% local variables:
% nn: length(x)
% kk: loop variable
% epsl: to 4th power is approximate 0..adjust with precision
% xx: temporary float
% inputs:
% x: argument
% outputs: sinc''(x)

epsl = 0.0001;
nn = length(x);
retval = zeros(size(x));
for kk = 1:nn
  xx = pi*x(kk);
  xxx = xx*xx;
  if abs(xx) > epsl
    retval(kk) = pi*pi*(2*(sin(xx)-xx*cos(xx))-xxx*sin(xx))/(xx*xxx);
  else
    retval(kk) = pi*pi*(xxx/10.0-1.0/3.0);
  end
end

endfunction

%

function retval = Skh(k,h,z)
% usage: Skh(k,h,z)
%description:
% this computes Skh(k,h,z)..it assumes h not 0, z a vector

retval = sinc(z/h-k);

endfunction

%

function retval = Skhp(k,h,z)
% usage: Skhp(k,h,z)
%description:
% this computes Skh(k,h,z)..it assumes h not 0, z a vector.

retval = sincp(z/h-k)/h;

endfunction


function retval = Skhpp(k,h,z)
% usage: Skhpp(k,h,z)
%description:
% this computes Skh(k,h,z)..it assumes h not 0, z a vector

retval = sincpp(z/h-k)/h^2;

endfunction

%

function retval = Cvh(m,n,h,v,z)
% usage: Cvh(m,n,h,v,z)
% description:
% computes truncated cardinal function at vector z
% with coefficients in vector v.

% local variables:
% kk: index variable
% inputs:
% m: lower limit of cardinal sum
% n: upper limit of cardinal sum
% h: step size
% v: array of values of v at appropriate nodes
% z: the argument of the cardinal function
% outputs:
% returns the value of cardinal function at z

retval = zeros(size(z));
for kk  = -m:n
  retval = retval + v(kk+m+1)*Skh(kk,h,z);
end

endfunction

%

function retval = Cvhp(m,n,h,v,z)
% usage: Cvhp(m,n,h,v,z)
% description:
% computes derivative of truncated cardinal function
% at vector z with coefficients in vector v.

% local variables:
% kk: index variable
% inputs:
% m: lower limit of cardinal sum
% n: upper limit of cardinal sum
% h: step size
% v: array of values of v at appropriate nodes
% z: the argument of the cardinal function
% outputs:
% returns the value of der cardinal function at z

retval = zeros(size(z));
for kk  = -m:n
  retval = retval + v(kk+m+1)*Skhp(kk,h,z);
end

endfunction

%

function retval = Cvhpp(m,n,h,v,z)
% usage: Cvhpp(m,n,h,v,z)
% description:
% computes second derivative of truncated cardinal function
% at vector z with coefficients in vector v.

% local variables:
% kk: index variable
% inputs:
% m: lower limit of cardinal sum
% n: upper limit of cardinal sum
% h: step size
% v: array of values of v at appropriate nodes
% z: the argument of the cardinal function
% outputs:
% returns the value of der^2 cardinal function at z

retval = zeros(size(z));
for kk  = -m:n
  retval = retval + v(kk+m+1)*Skhpp(kk,h,z);
end

endfunction

%

function retval = capprox(m,n,h,d,gamnodes,catzero,z)
% usage:capprox(m,n,h,d,gamnodes,catzero,z)
% description:
% computes approximate solution to Mueller-Shores problem at z.
% Caution: assumes weight fcn gam of paper exists.
% Also assumes phi and transformation fcns q0, q1.

% local variables:
% kk: index variable
% phiz: phi(z)
% inputs:
% m: lower limit of cardinal sum
% n: upper limit of cardinal sum
% h: step size
% d: array of values of c0 at appropriate nodes
% z: the argument of the cardinal function
% outputs:
% returns the value of der^2 cardinal function at z

retval = zeros(size(z));
% first compute c0(z)
phiz = arrayf('phi',z);
for kk  = -m:n
  retval = retval+(d(kk+m+1)/gamnodes(kk+m+1))*arrayf('gam',z).*...
           Skh(kk,h,phiz);
end
% now use c(z)=c0(z) + q0(z) + catzero*q1(z)
retval = retval + arrayf('q0',z) + catzero*arrayf('q1',z);

endfunction