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

% check that arrays of one element, empty arrays work

s = size(f);
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; 
f = permute(f,p);   % switch dimensions dim and 1

% code for dim=1
N = size(f,1) + 1;
S = size(f); S(1) = 1;
z = zeros(S);
G = [z;f;z;-flip(f,1)];      % odd extension of each column of length n=N-1 to length 2*N
Gh = fft(G,[],1)/(2*N);      % apply four() wrt to dimension 1
b = zeros(size(f));
b(:,:) = 1i*(Gh(2:N,:)-Gh(end:-1:N+2,:));

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