%
%------------------------------------------------------------------
% function [map]=FT_read_surface_map(filename)
%
% Function for Matlab, reads a mirror surface map as written
% by FT_write_surface_map.
% 
% filename: name of file to be written
% map:      the surface map (variable must be a 'map' structure,
%           see FT_new_surface_map)
%
% Andreas Freise  09.07.2008
%------------------------------------------------------------------
%

% TODO: the reading of the header could be more robust with fgetl and
% textscan

function [map] = FT_read_surface_map(filename)
  
  baseid='read_surface_map';
    
  % open file for reading
  [FID, result]=fopen(filename,'r');
  
  if (FID==-1)
    msgid=[baseid,':fileopen'];
    error(msgid,result);
  end
  
  % find the first line of the header
  found=FT_find_text_in_file(FID,'% Surface map');

  if (~found)
    msgid=[baseid,':fileread'];
    result='Could not find map data header in map file';
    error(msgid,result);
  end
  
  % start reading the header. Some general remarks:
  % - a '%%' in the fscanf string means a '%' in the file
  % - the line break is not read properly by fscanf, therefore
  %   each fscanf is followed by a fgetl(FID) which simple reads
  %   in the rest of the current text line 
  % - when multiple numbers are to be read fscanf works, if multiple
  %   strings are to be read it does not. In that case I use
  %   textscan. Textscan returns 'cells' rather than the string
  %   itself. Thus we have to convert them into strings with the
  %   cell2mat command.
  
  name=fscanf(FID,'%% Name: %s');
  fgetl(FID);
  
  % this command must be done without the leading '%%'. I don't know why
  tstr=fgetl(FID);
  [stype]=textscan(tstr,'%% Type: %s %s');
  %fgetl(FID);
  
  smaptype=cell2mat(stype{1});
  smapfield=cell2mat(stype{2});
  
  type=0;
  if (strcmp(smaptype,'phase'))
    type=0;
  elseif (strcmp(smaptype,'absorption'))
    type=1;
  elseif (strcmp(smaptype,'reflectivity'))
    type=2;
  end
  
  field=0;
  if (strcmp(smapfield,'both'))
    field=0;
  elseif (strcmp(smapfield,'reflection'))
    field=1;
  elseif (strcmp(smapfield,'transmission'))
    field=2;
  end
    
  [rc,count]=fscanf(FID,'%% Size: %d %d');
  fgetl(FID);
  rows=rc(1);
  cols=rc(2);

  [xy]=fscanf(FID,'%% Optical center (x,y): %g %g');
  fgetl(FID);
  x0=xy(1);
  y0=xy(2);
  
  [xystep]=fscanf(FID,'%% Step size (x,y): %g %g');
  fgetl(FID);
  xstep=xystep(1);
  ystep=xystep(2);

  scaling=fscanf(FID,'%% Scaling: %g');
  fgetl(FID);

  % finished reading the header, now reading the data.

  % read the numerical data of the map as one long vector
  tmp_data=fscanf(FID,'%g ');
  
  % reshape data into a 2D matrix
  data=reshape(tmp_data,cols,rows)';
  %data=tmp_data;
  
  % close the file 
  fclose(FID);
  
  % create map structure from the data
  map=FT_new_surface_map(name,data, type, field, x0, y0, xstep, ystep, scaling);
  

