## Copyright (C) 2010 Kurnia Wano, Levente Torok
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 2 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program; If not, see <http://www.gnu.org/licenses/>.

## -*- texinfo -*-
##
## @deftypefn {Function File} savevtk ( @var{X}, @var{Y}, @var{Z}, @var{filename} )
##  savevtk Save a 3-D scalar array in VTK format.
##  savevtk( array, filename  ) saves a 3-D array of any size to
##  filename in VTK format. This file format is used by Mayavi2 for example.
##
##  Example:
##    n = 30;
##    X = zeros(n,n,n);
##    for x = 1:n
##      for y = 1:n
##        for z = 1:n
##            X(x, y, z) = 1/sqrt ( x*x + y*y + z*z );
##        endfor
##      endfor
##    endfor
##    X = X * 200 / max(max(max(X)));
##    savevtk( X, "spherical.vtk");
##

##
## Author: Kurnia Wano, Levente Torok <address@hidden>
## Created: 2010-08-02 matlab version
## Updates: 2010-11-03 octave adoptatoin
##

function savevtk(array, filename)
 
    dims = size (array)
    if ( size (dims) != 3)
        error ("Save Vtk requires a 3 dimensional array");
    endif
    [nx, ny, nz] = size (array);

    [fid, msg] = fopen (filename, 'wt');
    if (fid==-1)
        error ("Cannot open file for saving file. %s", msg);
    endif

    fprintf (fid, '# vtk DataFile Version 2.0\n');
    fprintf (fid, 'Comment goes here\n');
    fprintf (fid, 'ASCII\n');
    fprintf (fid, '\n');
    fprintf (fid, 'DATASET STRUCTURED_POINTS\n');
    fprintf (fid, 'DIMENSIONS    %d   %d   %d\n', nx, ny, nz);
    fprintf (fid, '\n');
    fprintf (fid, 'ORIGIN    0.000   0.000   0.000\n');
    fprintf (fid, 'SPACING    1.000   1.000   1.000\n');
    fprintf (fid, '\n');
    fprintf (fid, 'POINT_DATA   %d\n', nx*ny*nz);
    fprintf (fid, 'SCALARS scalars double\n');
    fprintf (fid, 'LOOKUP_TABLE default\n');
    fprintf (fid, '\n');
    for a=1:nz
        for b=1:ny
            for c=1:nx
                fprintf(fid, '%d ', array(c,b,a));
            endfor
            fprintf(fid, '\n');
        endfor
    endfor
    fclose (fid);

endfunction