%
%--------------------------------------------------------------------------
% function [map]=FT_read_zygo_map(filename, filetype)
%
% A function for Matlab which reads a mirror surface map in the zygo
% format.  The function reads in two different file formats, either ascii
% (.asc) or .xyz.  The function returns a mirror map in the standard
% structure, rescaling the surface height data into nanometres.
%
% filename: name of file
% filetype: 1:  xyz
%           2:  ascii
%
% map:      the surface map
%
% Part of the Simtools package, http://www.gwoptics.org/simtools
% 27.10.2010    Charlotte Bond
%--------------------------------------------------------------------------
%

% Description: Reads a mirror surface map in Zygo interferometer format
% Keywords: zygo, mirror, map, read

function [map] = FT_read_zygo_map(filename,filetype)

  baseid='read_zygo_map';
    
  % Open file for reading
  [FID, result]=fopen(filename,'r');
    
  if (FID==-1)
    msgid=[baseid,':fileopen'];
    error(msgid,result);
  end
  
  % Value of surface height data (for ascii files) which indicates a map
  % element outside the mirror surface.
  zygo_data_lim=2147483640;
  
  name=FT_remove_file_extension(filename);
  
  % Standard maps have type 'phase' (they store surface
  % heights)
  type=0;
  
  % Both (reflected and transmitted) light fields are affected 
  field=0;
  
  % Scaling for nanometers
  scaling=1e-09;
 
  % Read header
  for i=1:2
    [M_tmp,count_tmp]=fscanf(FID,'%g');
    tmp=fgets(FID);
  end
  
  % Set values from header
  [M,count]=fscanf(FID,'%g');
  x0=M(7);
  y0=M(8);
  cols=M(9);
  rows=M(10);
  tmp=fgets(FID);
  
  % Ascii files contain intensity data aswell
  if filetype==2
      Icols=M(3);
      Irows=M(4);
  end
  
  % Read rest of header
  for i=1:2
    [M_tmp,count_tmp]=fscanf(FID,'%g');
    tmp=fgets(FID);
  end
  
  % Set scaling factors from header
  [M,count]=fscanf(FID,'%g');
  tmp=fgets(FID);
  
  % S: interferometric scale factor
  % lambda: wavelength
  % O: obliquity factor
  S=M(2);
  lambda=M(3);
  O=M(5);
  
  xstep = M(7);
  ystep = M(7);
  
  [M,count]=fscanf(FID,'%g');
  tmp=fgets(FID);
  
  % Resolution of phase data points, 1 or 0
  PhaseRes=M(11);
  
  if PhaseRes==0
    R=4096;
  elseif PhaseRes==1
    R=32768;
  end
  
  if filetype==1
      % .xyz files give phase data in microns
      Hscale=1e-6;
  elseif filetype==2
      % .asc files give phase data in internal units.  To convert to m use
      % Hscale factor
	  Hscale=S*O*lambda/R;
  end
  
  % Read rest of header
  [M_tmp,count_tmp]=fscanf(FID,'%g');
  tmp=fgets(FID);
  
  if filetype==2
      [M_tmp,count_tmp]=fscanf(FID,'%g');
      tmp=fgets(FID);
  end
 

  if filetype==1
      % .xyz file
      count=1;
      data=linspace(0,0,cols*rows);
      % Timer settings
      total_runs=cols*rows;
      elapsed_time=0;
      tic;
        
      % Read in data
      while count <= cols*rows
       
            [M_tmp,count_tmp]=fscanf(FID,'%g ');
            % Skip data if only two numbers on line (x and y coordinates).
            if count_tmp==2
                data(count)=nan;
                count=count+1;
            % End of data if end of numerical data
            elseif count_tmp==0
                count=cols*rows+1;
            % Otherwise read in data
            else
                points=floor(count_tmp/3);
                for p=1:points
                    data(count)=M_tmp(3*p);
                    count=count+1;
                end
          
                if points~=(count_tmp/3)
                    data(count)=nan;
                    count=count+1;
                end
            end
            tmp=fgets(FID);
            
            % Timer
            if (rem(count,2000)==0)  
                newline=0;
                [elapsed_time]=FT_print_progress(elapsed_time,count,total_runs,0);
                tic;
            end
      
      end
      [elapsed_time]=FT_print_progress(elapsed_time,count,total_runs,0);
      fprintf('\n');
  
  % Ascii file
  else
      
    % Read in intensity data
    [M,count]=fscanf(FID,'%g ');
    Idata=M(1:count);
    
    % Reshape intensity data
    Idata=rot90(reshape(Idata,Icols,Irows));
    tmp=fgets(FID);
    
    % Read in phase data
    [M,count]=fscanf(FID,'%g');
    data=M(1:count);
    
    % Set invalid elements to nan
    idx=find(data>=zygo_data_lim);
    data(idx)=nan;
    
  end
    
  % Reshape phase data and rotate it to corespond to (0,0) in bottom left
  % hand corner
  data=rot90(reshape(data,cols,rows));
  
  % Rescale data to nanometers
  data=data*Hscale/scaling;
  
  % Create new map structure with given data
  map=FT_new_surface_map(name,data,type,field,x0,y0,xstep,ystep,scaling);
  

end