%% savexyzmat.m
%%
%% Transform vectors x and y and cartesian matrix z (i.e., z=z(x,y))
%% into a resulting matrix with 3 columns, something like that written
%% to ascii files with 3 columns (like for gnuplot -- but there are no
%% blanks lines here of course).
%%
%% Usage: [xyzmat = ] savexyzkssmat(x, y, z, filename [, 'fastx' | 'fasty']);
%%
%% Filename can be either:
%%	- a file name: then this file is opened, written and closed;
%%	- a stream, e.g. stdout, stderr, or a file descriptor from fopen()
%%	- an empty string '' -- then, do no output (why you call this routine?)
%%
%% Option 'fastx' writes the matrix with x-column running faster:
%%	[ x1 y1 z11; x2 y1 z21; x3 y1 z31; ...; x1 y2 z12; ... ]
%% Option 'fasty' writes the matrix with y-column running faster:
%%	[ x1 y1 z11; x1 y2 z12; x1 y3 z13; ...; x2 y1 z21; ... ]
%% Default is 'fasty'.
%%
%% Example:
%%	x=1:2; y=1:4; [xx,yy]=meshgrid(y,x); % Cartesian: reversed order here!
%%	z=xx.+yy; 
%%	savexyzkssmat(x,y,z, 'out_fastx.dat', 'fastx')
%%	savexyzkssmat(x,y,z, 'out_fasty.dat', 'fasty')
%%
%% Author: Petr Mikulik
%% Version: 1. 6. 2009
%%
%% Compatibility: Octave 3, Matlab
%% History:
%%     1. 6. 2009: Added demo. Tested with Octave 3.1 and Matlab.
%%     7. 6. 2004:  Documentation update. Improved options.
%%     5. 2. 2001: Original version.


function xyz = savexyzkssmat ( x, y, z, filename, fastwhat )

if (nargin < 4 | nargin > 5)
    error('Usage: savexyzkssmat(x, y, z, filename [, opt])');
end

if ~is_vector(x) | ~is_vector(y) | ~is_matrix(z)
    error('x, y must be vectors and z cartesian matrix');
end
xlen = length(x);
ylen = length(y);
if [xlen, ylen] ~= size(z)
    error('cartesian matrix z must be of size length(x) times length(y)');
end

want_fasty = 1;
if nargin==5 && strcmp(fastwhat,'fastx')
    want_fasty = 0;
end

xyz=0;

want_output=1;
if isstr(filename)
    if length(filename)>0 f = fopen(filename,'w');
    else want_output=0;
    end
elseif isstream(filename) % FILE* passed
    f=filename;
end

xlen = length(x);
ylen = length(y);
if want_fasty	% case fasty
    for ix=1:xlen
	% [x(ix)*ones(1,ylen); y; z(ix,:)];
	if want_output
	    fprintf(f, '%g\t%g\t%g\n', [x(ix)*ones(1,ylen); y; z(ix,:)]);
	    fprintf(f, '\n'); % a blank line after each scan
	end
    end
else		% case fastx
    for iy=1:ylen
	% [x; y(iy)*ones(1,xlen); z(:,iy)'];
	if want_output
	    fprintf(f, '%g\t%g\t%g\n', [x; y(iy)*ones(1,xlen); z(:,iy)']);
	    fprintf(f, '\n'); % a blank line after each scan
	end
    end
end

if isstr(filename) & length(filename)>0
    fclose(f);
end

if nargout > 0
    if ~want_fasty
	[yy,xx]=meshgrid(y,x);
	xx=xx(:);
	yy=yy(:);
	zz=z(:);
	xyz = [xx, yy, zz];
    else
	[xx,yy]=meshgrid(x,y);
	xx=xx(:);
	yy=yy(:);
	zz=z(:);
	xyz = [xx, yy, zz]
    end
end


%!demo
%! x=1:2; y=1:4; [xx,yy]=meshgrid(y,x); % Cartesian: reversed order here!
%! z=xx.+yy;
%! savexyzkssmat(x,y,z, 'out_fastx.dat', 'fastx')
%! savexyzkssmat(x,y,z, 'out_fasty.dat', 'fasty')


% eof savexyzkssmat.m