function S=m_loadsir(filename,varargin)
%
% function S=loadsir(filename,<options>)
%
% function to load a sir file and optionally display
% returns a structure S that includes data array and other information
%
% inputs: 
%  filename - name of file to read
%  options  - 'show' with default projection in currrent figure
%           - 'nolatlon' do not create lat/lon arrays, overrides show
%           - 'noset' do not set projection, overrides show
% 
% Note: precise SIR projections are imperfectly represented in M_map.
% M_map converted SIR files are adequate for plotting, but use SIR
% projections for precise pixel level operations

% written by D.G. Long 21 May 2020

global MAP_PROJECTION

% set run options
show=0;
nolatlon=0;
noset=0;
if nargin > 1
  for k=1:nargin-1
    arg=lower(varargin{k});
    if strcmp(arg,'show')
      show=1;
    elseif strcmp(arg,'nolatlon')
      nolatlon=1;
    elseif strcmp(arg,'noset')
      noset=1;
    end
  end
end

S=[]; % initialize empty output structure
[S.im S.head]=loadsir(filename); % read file

% decode projection information from header
S.iopt=sirheadvalue('iopt',S.head);
nsx=sirheadvalue('nsx',S.head);
nsy=sirheadvalue('nsy',S.head);
S.sh=0; % southern hemisphere azimuth when set to 1
switch S.iopt
  case -1 % image only (no projection information)
    S.projname='None';
    S.proj='None';
  case 0  % rectalinear
    S.projname='Equidistant Cylindrical';
    S.proj='equidistant';
    [x,y]=pix2latlon([1,nsx],[1,nsy],S.head);
    S.plon=[min(x),max(x)];
    S.plat=[min(y),max(y)];
    S.sph='normal';
  case 1  % Lambert on sphere 
    S.projname='Lambert Conformal Conic';
    S.proj='lambert';
    [x,y]=pix2latlon([1,nsx],[1,nsy],S.head);
    % convert to +/- 180 longitude range
    x(x>180)=x(x>180)-360;
    x(x<-180)=x(x<-180)+360;
    S.plon=[min(x),max(x)];
    S.plat=[min(y),max(y)];
    S.clon=mean(S.plon);
    S.ell='sphere';
  case 2  % Lambert on elipsoid
    S.projname='Lambert Conformal Conic';
    S.proj='lambert';
    [x,y]=pix2latlon([1,1,nsx,nsx],[1,nsy,1,nsy],S.head);
    % convert to +/- 180 longitude range
    x(x>180)=x(x>180)-360;
    x(x<-180)=x(x<-180)+360;
    S.plon=[min(x),max(x)];
    S.plat=[min(y),max(y)];
    S.clon=mean(S.plon);
    S.ell='wgs72';  % approx 
  case 5  % Polar stereographic
    S.projname='Polar Stereographic';
    S.proj='stereographic';
    S.plon=0; % sirheadvalue('xdeg',S.head);
    S.plat=90*sign(sirheadvalue('ydeg',S.head));
    [x,y]=pix2latlon(1,nsy/2,S.head);
    S.rad=90-y*sign(sirheadvalue('ydeg',S.head));
    S.rot=0; %-sirheadvalue('xdeg',S.head);
    if S.plat<0, S.sh=1; end;
  case 8  % EASE2 NH
    S.projname='Azimuthal Equal-Area (EASE2 NH)';
    S.proj='azimuthal equal-area';
    S.plon=0;
    S.plat=90;
    S.rad=90;
  case 9  % EASE2 SH
    S.projname='Azimuthal Equal-Area (EASE2 SH)';
    S.proj='azimuthal equal-area';
    S.plon=0;
    S.plat=-90;
    S.rad=90;
    S.sh=1;
  case 10 % EASE2 Global
    S.projname='Cylindrical Equal-Area (EASE2 Global)';
    S.proj='cylindrical';
    [x,y]=pix2latlon([1,nsx],[1,nsy],S.head);
    S.plon=[min(x),max(x)];
    S.plat=[min(y),max(y)];
    S.sph='wgs84';
  case 11 % EASE1 NH
    S.projname='Azimuthal Equal-Area (EASE1 NH)';
    S.proj='azimuthal equal-area';
    S.plon=0;
    S.plat=90;
    S.rad=90;
  case 12 % EASE1 SH
    S.projname='Azimuthal Equal-Area (EASE1 SH)';
    S.proj='azimuthal equal-area';
    S.plon=0;
    S.plat=-90;
    S.rad=90;
    S.sh=1;
  case 13 % EASE1 Global
    S.projname='Cylindrical Equal-Area (EASE1 global)';
    S.proj='cylindrical';
    S.plon=[-180,180];
    S.plat=[-90,90];
    S.sph='wgs84';
  otherwise
    S.projname='M_map unsupported projection: unknown';
    S.proj='None';
    disp('*** invalid sir projection');
end
S.title=strtrim(sirheadvalue('title',S.head));

%fprintf('Filename: %s\n Projection: %s\n',filename,S.projname);
%disp(S.plon);
%disp(S.plat);
%disp(S.clon);
%disp(S.rad);
%disp(S.ell);
      
if ~noset
  % set M_map projection based on SIR file projection information
  switch S.iopt
    case 0  % rectalinear lat/lon grid (type of cylindrical)
      m_proj(S.proj,'lon',S.plon,'lat',S.plat,'sph',S.sph);
    case 2 % Lambert
      m_proj(S.proj,'lon',S.plon,'lat',S.plat,'clon',S.clon,'ell',S.ell);
    case 5 % Polar Stereographic (type of azimuthal)
      m_proj(S.proj,'lon',S.plon,'lat',S.plat,'rad',S.rad,'circle');
    case 8  % EASE2 NH (type of azimuthal)
      m_proj(S.proj,'lon',S.plon,'lat',S.plat,'rad',S.rad,'circle');
    case 9  % EASE2 SH (type of azimuthal)
      m_proj(S.proj,'lon',S.plon,'lat',S.plat,'rad',S.rad,'circle');
    case 10  % EASE2 global (type of cylindrical)
      m_proj(S.proj,'lon',S.plon,'lat',S.plat,'sph',S.sph);
    case 11  % EASE1 NH (type of azimuthal)
      m_proj(S.proj,'lon',S.plon,'lat',S.plat,'rad',S.rad,'circle');
    case 12  % EASE1 SH (type of azimuthal)
      m_proj(S.proj,'lon',S.plon,'lat',S.plat,'rad',S.rad,'circle');
    case 13  % EASE1 global (type of cylindrical)
      m_proj(S.proj,'lon',S.plon,'lat',S.plat,'sph',S.sph);
    otherwise
      fprintf('*** Cannot set M_map projection: %s\n',S.proj);
  end
end

if ~nolatlon
  % compute lat/lon arrays
  [m n]=size(S.im);
  [X,Y]=meshgrid(1:n,1:m);
  [S.lon,S.lat]=pix2latlon(X,Y,S.head);
  % convert to +/- 180 longitude range
  S.lon(S.lon>180)=S.lon(S.lon>180)-360;
  S.lon(S.lon<-180)=S.lon(S.lon<-180)+360;

  if show
    % use M_map to display image using default projection
    m_pcolor(S.lon,S.lat,S.im); shading flat;
    % add grid
    if S.sh
      % Note: for southern hemisphere, need to set xaxislocation to top
      m_grid('tickdir','out','xaxislocation','top');
    else
      m_grid('tickdir','out');
    end
    % add color bar
    colorbar;
    % add title
    title(S.title);
  end
end