function export2DLegacyVTK(XX,YY,fields,field_names,filename)

% Exports to VTK a field given at evenly spaced points.
% 
% XX, YY ------ are matrices (obtained from meshgrid) that contain the (x,y) coordinates
%
% fields ------ is a 3D matrix. The first and second indices indicate the 2D coordinates, while the third
%               index is used to stack together the fields. Each field is assumed to be a scalar, so if you
%               want to print a vector, you should give it every component separately.
% field_names - is a cell array of strings, each entry being the name of the respective fields. Used to 
%               assign a proper name to the field into the VTK file.
% filename ---- destination VTK file. You may want to put the ".vtk" extension.

Nx = size(XX,2);
Ny = size(XX,1);

x_min = XX(1,1);
x_max = XX(1,end);
y_min = YY(1,1);
y_max = YY(end,1);

dx = (x_max-x_min)/Nx;
dy = (y_max-y_min)/Ny;

Nfields = size(fields,3);

if (numel(field_names) ~= Nfields)
  error('Attention! Field names cell array "field_names" does not contain the same number of entries as the "fields" matrix. ABORTING')
end

% Open file in writing mode
fid = fopen(filename, 'w');

fprintf(fid,'# vtk DataFile Version 2.0\n')
fprintf(fid,'Created with export2DLegacyVTK.m\n')
fprintf(fid,'ASCII\n')
fprintf(fid,'DATASET STRUCTURED_POINTS\n')
fprintf(fid,'DIMENSIONS %d %d %d\n', Nx, Ny, 1)
fprintf(fid,'SPACING %e %e %e\n', dx, dy, 0.0)
fprintf(fid,'ORIGIN %e %e %e\n', x_min, y_min, 0.0)
fprintf(fid,'POINT_DATA %d\n', Nx*Ny)

for iii = 1:Nfields

  fprintf(fid,'SCALARS %s float 1\n', field_names{iii})
  fprintf(fid,'LOOKUP_TABLE default\n')

  for j = 1:Ny
    for i = 1:Nx

      fprintf(fid, '%e ', fields(j,i,iii))

    end
  end

  fprintf(fid, '\n') % new line

end

% % This is not needed:
% fprintf(fid, 'METADATA\n')
% fprintf(fid, 'INFORMATION 0\n')

fclose(fid);

return

close all
clear
clc

% Create domain
x = linspace(0,1,100);
y = linspace(0,0.5,100);

[XX,YY] = meshgrid(x,y);

% Create fields
field1 = sin(pi*XX);
field2 = sin(2*pi*XX).*cos(4*pi*YY);
field3 = sqrt(XX) + YY.^2;

% Assemble fields into a 3D matrix
all_fields(:,:,1) = field1;
all_fields(:,:,2) = field2;
all_fields(:,:,3) = field3;

field_names = {'field1', 'field2', 'field3'};

% Export fields to VTK
export2DLegacyVTK(XX,YY,all_fields,field_names,'output.vtk')