%
%--------------------------------------------------------------------------
% function [map]=FT_read_virgo_map(filename)
%
% A function for Matlab which reads a mirror surface map in the standard
% virgo format.
%
% filename: name of virgo file
% map:      the surface map (variable must be in the standard virgo
%           structure)
%
% Charlotte Bond, Andreas Freise  04.09.2008
%--------------------------------------------------------------------------
%

function [map] = FT_read_virgo_map(filename)

  baseid='read_virgo_map';
    
  % Open file for reading
  [FID, result]=fopen(filename,'r');
    
  if (FID==-1)
    msgid=[baseid,':fileopen'];
    error(msgid,result);
  end
  
  name=filename;
  
  % Standard virgo format maps have type 'phase' (they store surface
  % heights)
  type=0;

  % In virgo fromat maps both (reflected and transmitted) light fields are 
  % affected 
  field=0;
  
  % In the standard virgo format the measurements are given in nanometers
  scaling=1e-09;
  
  % Now read the data
  [M,count]=fscanf(FID,'%g ');
  
  if (count<8)
    msgid=[baseid,':fileread'];
    result='could not read map header';
    error(msgid,result);
  end
  
  % The first two numbers are the grid dimensions in unit steps
  cols=M(1);
  rows=M(2);
  
  % The next two numbers give the optical centre of the grid in coordinates
  x0=M(3);
  y0=M(4);
  
  % The next two numbers are zeros and are not used
  
  % The next number is the zero point number (for all points not on the
  % optical surface) given in nanometers
  %zero=M(7);  
  % sometimes this number is not set properly so we simply take the first data point
  % assuming that with a circular mirror the value always corresponds to the
  % zero point, see below
  
  % The next number is the step size given in meters.  The grid is square
  % so xstep and ystep are the same
  xstep=M(8);
  ystep=M(8);
  
  if(length(M)<9)
    % read dummy line of '=================='
    tmp=fgets(FID);
    [M_tmp,count_tmp]=fscanf(FID,'%g ');
    M=[M;M_tmp];
    count=count+count_tmp;
  end

  % set zero point, see above
  zero=M(10); 

  
  % Close the file
  fclose(FID);
    
  if (count<x0*y0)
    msgid=[baseid,':fileread'];
    result='number of data points not consistent with header information';
    error(msgid,result);
  end
  
  % Then the data follows.  Sometime the number before the data is 0.0 in 
  % which case we have to ignore this point.  If it is the first data point
  % it should be equal to the zero point.  In the case where the first 
  % number is the first data point, the count of the data differs for 
  % different maps.  If the last number is 0.0 then the last data point is
  % at M(count-1).  Otherwise it is M(count).
  
  if M(count)==0.0
      data_count=count-1;
  else
      data_count=count;
  end
  
  % We skip the first point if it is not data
  if M(9)==zero
      tmp_data=M(9:data_count);
  else
      tmp_data=M(10:count);
  end
  
  % Reshape the data into a 2D matrix:
  % - the reshape creates a matrix from the vector.
  % - however, it has x and y exchanged, therefor we take
  %   the transpose. 
  % - further, the virgo map format starts
  %   in the top left corner. More intuitive for using the
  %   data in Matlab or Finesse would be to start in the
  %   lower left corner. Thus we apply 'flipud'
  
  
  [tmp_map_idx]=find(abs(tmp_data)<zero);
  [tmp_idx]=find(abs(tmp_data(tmp_map_idx))>0);
  tmp_data_end=tmp_map_idx(tmp_idx(numel(tmp_idx)));

  % find the extra 0 valued points at the end of the map files (outside the 
  % mirror surface) and set them to the value of the zero point height
  if (zero~=0)
      for i=tmp_data_end+1:numel(tmp_data)
          if (tmp_data(i)~=zero)
              if (i>=tmp_data_end+(1.5*rows));
                tmp_data(i)=zero;
              end
          end
      end
  end
  
  data=flipud(reshape(tmp_data,cols,rows)');
  
  
  % In the standard virgo format the data not on the surface is set to
  % an arbitrary 'zero point' value. We need to set all the points which 
  % are the same numerical value as the zero point to 'not a number'
  zerodata=find(data==zero);
  data(zerodata)=nan;
  
  % Create map structure from the data
  map=FT_new_surface_map(name,data,type,field,x0,y0,xstep,ystep,scaling);
  
