function f = isintr(b,dim)
% f = isintr(b)
% f = isintr(b,dim)      % apply isintr for dimension dim of array
% 
% Inverse Discrete Sine Transform
% given coefficients b=[b_1,...,b_{N-1}] of sine series sum_{k=1...N-1} b_k sin(k*pi*x)
% find values f=[f_1,...,f_{N-1}] at x=1/N,...,(N-1)/N

s = size(b);
if nargin==1
  dim = find(s>1,1);
  if isempty(dim)
    dim = 1;
  end
end

p = 1:length(s);
p(dim)=1; p(1)=dim; 
b = permute(b,p);   % switch dimensions dim and 1

% code for dim=1
N = size(b,1) + 1;
S = size(b); S(1) = 1;
z = zeros(S);
Gh = [z;b;z;-flip(b,1)]/(2i);   % use sin(t) = (exp(it)-exp(-it))/(2i)
G = ifft(Gh,2*N,1)*(2*N);       % apply ifour() to each column
f = zeros(size(b));
f(:,:) = G(2:N,:);

% switch dimensions back
f = permute(f,p);
