

function [x, y, z] = ss_block(xr, yr, nsub)

% function [x, y, z] = ss_block(xr, yr, nsub)
% x, y, and z are lon, lat, and depth from the Smith/Sandwell
% topography, version 6.2, in the lon range xr, and lat range yr.
% If nsub is omitted, the full resolution is provided.
% Set nsub = 3 to get every third point, or 9 for every 9th point.
%
% This is based on m_tbase by Rich Pawlowizc of UBC.
%
%  Eric Firing
% changed path for SIO. Lisa Sep99.

PATHNAME = '/Users/jhooper/Documents/Work_data/lou_toolbox/sand/topo_8.2.img'

% PATHNAME = 'C:\MATLAB7p1\toolbox\m_map\etopo2_new.raw'




if nargin == 2,

   nsub = 1;

elseif nargin ~= 3,

   error('Number of arguments must be 2 or 3.')

end

if nsub == 3 | nsub == 9,

  % topofile = sprintf([PATHNAME], nsub);
    topofile = sprintf('%s',[PATHNAME]) %mod by R.Smith 10/04/06
  
elseif nsub == 1,

   topofile = [PATHNAME];

else

   error('If specified, nsub must be either 1 or 3 or 9.');

end



p.pts_per_degree = 30/nsub;

p.nlong = 360 * p.pts_per_degree;

p.nlat = 6336/nsub;



% If the left side is negative, shift the entire range by 360 to make

% all longitudes positive.

if xr(1) < 0,

   xr = xr + 360;

   wrapped = 1;

else

   wrapped = 0;

end





[ii,jj] = imgxy2ij(xr, yr, p)

nlongz = diff(ii) + 1;

nlatz =  diff(jj) + 1;

% Note: z is transposed relative to the input array; Matlab

% surface and contour functions want Z(lat, lon).

z = zeros(nlatz, nlongz);



% If the right side wraps past 360 degrees, break longitudes into

% two pieces, the first (ii(1):ii(2)) going up to just befor 360,

% and the second (ii(3):ii(4)) going from 0 to the right limit.

if ii(2) > p.nlong

   ii = [ii(1); p.nlong; 1; ii(2)-p.nlong];

end



% More informative error messages can be added to the following:

if (any(ii<1 | ii>p.nlong))

   error('longitude index out of range');

end

if (any(jj>p.nlat | jj<1))

   error('latitude index out of range');

end



[fid, msg] = fopen(topofile, 'rb', 'b');   % BigEndian file

if fid == -1,

   error(['Can''t open ', topofile, msg]);

end



if length(ii) == 2,      % No crossing of the prime meridian; single read.

   [x, y] = imgij2xy(ii(1):ii(2), jj(2):jj(1), p);

   if wrapped, x = x - 360; end

   iz = 1;

   for jlat = jj(2):jj(1)

      fseek(fid, 2 * (ii(1)-1 + (jlat-1)*p.nlong), 'bof');

      z(iz,:) = fread(fid, [1 nlongz], 'int16');

      iz = iz + 1;

   end

else                     % Crosses the prime meridian; 2 reads per latitude

   [x, y] = imgij2xy([ii(1):ii(2) ii(3):ii(4)], jj(2):jj(1), p);

   nlongz1 = diff(ii([1 2])) + 1;

   nlongz2 = diff(ii([3 4])) + 1;

   jz1 = 1:nlongz1;

   jz2 = nlongz1 + (1:nlongz2);

   if wrapped,  x(jz1) = x(jz1) - 360; end % x(jz2) was unwrapped

                                           % when ii was split in two.  

             

   iz = 1;

   for jlat = jj(2):jj(1)

      fseek(fid, 2 * (ii(3)-1 + (jlat-1)*p.nlong), 'bof');

      z(iz,jz2) = fread(fid, [1 nlongz2], 'int16');

      fseek(fid, 2 * (ii(1)-1 + (jlat-1)*p.nlong), 'bof');

      z(iz,jz1) = fread(fid, [1 nlongz1], 'int16');

      iz = iz + 1;

   end

end



fclose(fid);



%%--------- end of ss_block function ------------------------------



function [ii, jj] = imgxy2ij(x,y, p);



tanlat = tan(y*pi/360);  % tan(lat/2)

jj = round(0.5 + p.nlat/2 - p.pts_per_degree*180/pi * ...
log( (1+tanlat)./(1-tanlat) ) );

     % 0.5 + p.nlat/2 is the location of the equator in "index space";

     % Note that the grid straddles the equator, because it is symmetric

     % and has an even number of latitudes.

ii = round(0.5 + eps + x*p.pts_per_degree);



%%-------------------------------------------------------------------



function [x, y] = imgij2xy(ii, jj, p)

expy = exp((jj-0.5-p.nlat/2)*pi/(p.pts_per_degree*180));

y = -2*(180/pi)*atan( (expy-1)./(expy+1) );

x = (ii-0.5)/p.pts_per_degree;



%%-------------------------------------------------------------------





