function lyngby_write_sdt(filename, X, siz, arg1, ...
    arg2, arg3, arg4, arg5, arg6, arg7);

% lyngby_write_sdt     - Write a volume/datamatrix to SDT/SPR file
%
%        function lyngby_write_sdt(filename, X, siz, ...
%                'PropertyName', 'PropertyValue');
%                
%        Input:    filename  Filename 
%                  X         Volume or datamatrix
%                  siz       x,y,z,[t] Number of pixels and timepoints
%
%        Property: VoxelDim  {[0.001 0.0001 0.0001 1]} Voxel size
%                            in meters and time intervals in
%                            seconds ('interval' in the SPR)
%                  Precision [ 8 | {16} | 32 | 64 ] Precision
%                            (bit per voxel).
%                  Datatype  [ float | {signed} ] 
%                            If float is off then the default
%                            precision is 16
%                  Endian    [ {big} | little ] Byte ordering
%                            within the word 
%                  SdtOrient [ {ax} | cor | sag ] Orientation of
%                            the slices
%                  Fov       [ {voxelDim * siz} ] Field of view
%                  Origin    [ {-fov/2+voxelDim/2} ] Position of
%                            the center of the first voxel 
%                  Scale     [ {0} | 1 ] Scaling of pixel
%                            values. Does not apply in connection
%                            with 'float'.
% 
%       This functions write Stimulate SDT (data) and SPR (header)
%       files. It is not able to handle write the COMPLEX datatype. 
%
%       Ref: Stimulate, http://www.cmrr.umn.edu/stimulate/
%            http://www.cmrr.umn.edu/stimulate/stimUsersGuide/
%
%       See also LYNGBY, LYNGBY_READ_SDT, LYNGBY_WRITE_ANA,
%                LYNGBY_WRITE_VAPET. 
%
% $Id: lyngby_write_sdt.m,v 1.6 2003/03/04 15:38:20 fnielsen Exp $


    if ~isstr(filename)
      error('''filename'' argument should be a string')
    end

    if prod(size(X)) ~= prod(siz)
      error(sprintf('X and siz do not correspond. V: %d, siz: %d', ...
	  prod(size(X)), prod(siz)));
    end

    % Default properties
    numDim = length(siz);
    voxelDim = [ 0.001 0.001 0.001 1];
    precision = 16;
    bScale = 0;
    datatype = 'signed';
    endian = 'big';
    sdtOrient = 'ax';         % 'sag', cor, ax
    origin = [];
    fov = [];
    
    % Parse Properties
    n = 1;
    while n <= nargin-3
      eval(sprintf('arg = lower(arg%d);', n));
      if strcmp(arg, 'voxeldim')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	if length(arg) == length(siz)
	  voxelDim = arg;
	else
	  error(sprintf(['Argument with ''VoxelDim'' should be ' ...
		' the same length as siz. VoxelDim: %dx%d, siz: %dx%d'], ...
		size(arg), size(siz)));
	  end
      elseif strcmp(arg, 'precision')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	if isreal(arg) 
	  precision = arg;
	else
	  error('The argument with ''Precision'' should be a number.'); 
	end
      elseif strcmp(arg, 'datatype')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	if strcmp(lower(arg), 'float') 
	  datatype = 'float';
	  precision = 32;
	else
	  if strcmp(lower(arg), 'signed') 
	    datatype = 'signed';
	  %elseif strcmp(lower(arg), 'unsigned') 
	  %datatype = 'unsigned';
	  else
	    error(['The argument with ''Datatype'' should be a string: '...
		  '''float'', ''signed'' or ''unsigned''']); 
	  end
	  if ~(precision == 8 | precision == 16)
	    precision = 16;
	  end
	end
      elseif strcmp(arg, 'endian')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	if strcmp(lower(arg), 'big') 
	  endian = 'big';
	elseif strcmp(lower(arg), 'little')
	  endian = 'little';
	else
	  error(['The argument with ''Endian'' should be a string: '...
		  '''big'' or½ ''little'' ']); 
	end
      elseif strcmp(arg, 'scale')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	if arg
	  bScale = 1;
	else
	  bScale = 0;
	end
	
      elseif strcmp(arg, 'fov')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	if length(arg) == length(siz)
	  fov = arg;
	else
	  error(sprintf(['Argument with ''Fov'' should be ' ...
		' the same length as siz. Fov: %dx%d, siz: %dx%d'], ...
		size(arg), size(siz)));
	end
      
      elseif strcmp(arg, 'origin')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	if length(arg) == length(siz)
	  origin = arg;
	else
	  error(sprintf(['Argument with ''Origin'' should be ' ...
		' the same length as siz. Origin: %dx%d, siz: %dx%d'], ...
		size(arg), size(siz)));
	end
      
      elseif strcmp(arg, 'sdtorient')
	n = n + 1;
	eval(sprintf('arg = lower(arg%d);', n));
	if isstr(arg)
	  if strcmp(arg, 'ax')
	    sdtOrient = 'ax';
	  elseif strcmp(arg, 'sag')
	    sdtOrient = 'sag';
	  elseif strcmp(arg, 'cor')
	    sdtOrient = 'cor';
	  else 
	    error(['Argument to ''stdOrient'' should be ' ...
		  '''ax'', ''sag'' or ''cor''']);
	  end
	else
	  error('Argument to ''stdOrient'' should be a string');
	end

      else
	error(sprintf('Invalid property: %s', arg));
      end
      n = n + 1;
    end

    if strcmp(datatype, 'float')
      if precision ~= 32 
	error('When using float the precision should be 32.');
      end
    else
      if precision ~= 8 &  precision ~= 16 & precision ~= 32
	error(['When not using float the precision should ' ...
	      'be 8, 16 or 32']);
      end
    end
    
    scale = 1;
    if bScale & ~strcmp(datatype, 'float')
      maxabs = max(max(abs(X)));
      scale = (2^(precision-1)-1) / maxabs;
      X = scale * X; 
    end

    voxelDim = voxelDim(1:length(siz));
    
    if isempty(fov)
      fov = voxelDim .* siz;
    end
    
    if isempty(origin)
      origin = - fov/2 + voxelDim/2;
    end
   
    
    % Construct filename
    index = findstr(filename, '.sdt');
    if ~isempty(index) 
      if length(filename) - length('.sdt') + 1 == index
	filename = filename(1:index-1);
      end
    end
    index = findstr(filename, '.spr');
    if ~isempty(index) 
      if length(filename) - length('.spr') + 1 == index
	filename = filename(1:index-1);
      end
    end
    sdtFilename = [ filename '.sdt'];
    sprFilename = [ filename '.spr'];
    
    
    % Open SDT file
    if strcmp(endian, 'big')
      fid = fopen(sdtFilename, 'w', 'ieee-be');
    else
      fid = fopen(sdtFilename, 'w', 'ieee-le');
    end
    if fid == -1
      error(sprintf('Could not open file: %s', sdtFilename));
    end
    
    % Write SDT file
    if precision == 8
      count = fwrite(fid, X, 'char');
    elseif precision == 16
      count = fwrite(fid, X, 'int16');
    elseif precision == 32 & strcmp(datatype, 'float')
      count = fwrite(fid, X, 'int32');
    else
      count = fwrite(fid, X, 'float32');
    end
    if count ~= prod(size(X))
      error('Could not write all data')
    end
      
    fclose(fid);
    
    % Open SPR file
    fid = fopen(sprFilename, 'w');
    if fid == -1
      error(sprintf('Could not open file: %s', sdtFilename));
    end
    
    fprintf(fid, 'numDim:         %d\n', numDim);

    fprintf(fid, 'dim:            ');
    fprintf(fid, '%d ', siz);
    fprintf(fid, '\n');
  
    fprintf(fid, 'origin:         ');
    fprintf(fid, '%f ', origin(1:numDim));
    fprintf(fid, '\n');

    fprintf(fid, 'fov:            ');
    fprintf(fid, '%f ', fov(1:numDim));
    fprintf(fid, '\n');
    
    fprintf(fid, 'interval:       ');
    fprintf(fid, '%d ', voxelDim(1:numDim));
    fprintf(fid, '\n');

    fprintf(fid, 'dataType:       ');    
    if precision == 8
      fprintf(fid, 'BYTE');
    elseif precision == 16
      fprintf(fid, 'WORD');
    elseif precision == 32 & strcmp(datatype, 'signed')
      fprintf(fid, 'LWORD');
    else
      fprintf(fid, 'REAL');
    end
    fprintf(fid, '\n');
    
    fprintf(fid, 'displayRange:   %d %d\n', min(min(X)), ...
	max(max(X))); 

    fprintf(fid, 'sdtOrient:      ');
    fprintf(fid, '%s\n', sdtOrient);
    
    if bScale
      if precision == 16
        fprintf(fid, 'Real2WordScale: %f\n', scale);
      end
    end
    
    fclose(fid);













