function [polar_grid,x_grid,y_grid] = radial_grid ( matfile, ...
                                                    out_matfile, ...
						    num_coastlines, ...
						    decimate_coastline_factor )
% RADIAL_GRID:  Constructs radial grid.
%
% Using the data from matfile, construct a radial grid.  All
% that the user need do is click once on the screen to pinpoint the center of
% the grid.  All necessary information to do the gridding is saved in the
% file specified by "out_matfile".
%
% Usage:
%   [rgrid,xgrid,ygrid] = radial_grid ( matfile, 
%                                       out_matfile, 
%                                       num_coastlines,
%                                       decimate_coastline_factor );
%
% Input:
%   matfile:  Matlab *.mat file, contains location of data to grid.  The lat
%             and lon are needed.
%   out_matfile:  Name of *.mat file to save polar_grid, x_grid, and y_grid to.
%   num_coastlines:  Number of coastlines to actually use.  The first one is
%                the big one, so for testing purposes just use 1.
%   decimate_coastline_factor:  Thins out the coastline vector by this much.
%                The coastlines are really big, so taking every 10th point
%                isn't such a bad idea.
%
% Output:  ALL output arguments are optional.
%   polar_grid:  Grid with land areas masked out as -Infs.  Valid grid cells
%                are flagged with +Infs.
%   x_grid, y_grid:  Has x and y locations of valid cells in polar_grid.
%
% If the program fails, a null matrix is returned for polar_grid.
%
% John Evans (jevans@sakhalin.er.usgs.gov)
%

if ( nargin < 2 )
  help radial_grid;
end

if ( nargin < 3 )
  num_coastlines = 1;
end
if ( nargin < 4 )
  decimate_coastline_factor = 10;
end


fig = figure;


eval ( sprintf ( 'load %s', matfile ) );


% Load the coastline file.
load wgom_half.ll;
plot( wgom_half(:,1), wgom_half(:,2), '-c' );
set(gca, 'aspectratio', [1 cos(43*pi/180)] );
hold on


% plot the data, uncomment all this stuff when done debugging
hold on
if ( ~exist('lat') | ~exist('lon') )
  fprintf (1, 'No coastline variables lat or lon in matfile.\n' );
  polar_grid = [];
  return;
else
  plot ( lon, lat, '-c' );
  plot ( lon, lat, '*g' );
end


% get the center of the grid from the mouse
[c_lon, c_lat] = ginput(1);
plot(c_lon, c_lat, 'oc');
centered_lon = lon - c_lon;
centered_lat = lat - c_lat;

% convert the data coords into polar
data_r = sqrt ( centered_lon.^2 + centered_lat.^2 );
data_theta = atan2( centered_lat, centered_lon );
ifind = find(data_theta<0);
data_theta(ifind) = data_theta(ifind) + 2*pi;


% get first radius from the mouse
r1 = min(data_r) * 0.95;
r2 = max(data_r) * 1.05;
theta1 = min(data_theta)*0.95;
theta2 = max(data_theta)*1.05;



% dr should be about 3km.  Need to translate this into degrees lat and lon.
dr = 0.0362075;

% dA should be about 9 square kilometers.  This is approximately dr^2 in
% terms of "square degrees" (my math profs would twist in their graves over
% this, too bad they're not dead).
dA = 9.7793e-4;

% Now get what dtheta should be.
dtheta = dA / ( (0.5 * (r1 + r2)) * dr );




% construct the grid
% Determine # of discrete values of r, theta.
nr = ceil ((r2 - r1)/dr + eps) + 1;
ntheta = ceil ((theta2 - theta1)/dtheta + eps) + 1;

% vector of radii, theta
rv = [r1:dr:r1+(nr-1)*dr];
thetav = [theta1:dtheta:theta1+(ntheta-1)*dtheta];

[RV,THETAV] = meshgrid(rv,thetav);
polar_grid = 0 * ones(size(RV));



% We want a grid in terms of our original data.  Use this to get a decent
% view as to what the grid actually looks like.
x_grid = RV .* cos(THETAV) + c_lon;
y_grid = RV .* sin(THETAV) + c_lat;
xy_grid = 0 * ones(size(x_grid));



% outline of grid.
% This is necessary to prune the coastlines.  By checking whether or not a
% coastline lies within these values of r and theta, we can quickly discard
% any coastline that lies outside the grid.  This is faster than checking
% all the grid cells to determine if they lie within EVERY SINGLE coastline.
pgrid_bndry = [r1 theta1; r1 theta2; r2 theta2; r2 theta1];


% convert coastline data into polar.  First center around user origin.
centered_xy_coast(:,1) = wgom_half(:,1) - c_lon;
centered_xy_coast(:,2) = wgom_half(:,2) - c_lat;

% compute r and theta, adjusting theta to be between 0 and 2*pi
polar_coast(:,1) = sqrt( centered_xy_coast(:,1).^2 + centered_xy_coast(:,2).^2 );
polar_coast(:,2) = atan2 ( centered_xy_coast(:,2), centered_xy_coast(:,1) );
ifind = find(polar_coast(:,2) < 0);
polar_coast(ifind,2) = polar_coast(ifind,2) + 2*pi;



% for each separate coastline, set all grid points contained within equal to
% -Inf
coast = wgom_half;

% The "breaks" in the coastline files are marked by NaNs.  By counting
% these, we determine how many coastlines there are.
coastline_breaks = find(isnan(coast(:,1)));

% make sure we get the last coastline
n = length(coastline_breaks);
coastline_breaks(n+1) = length(polar_coast(:,1)) + 1;


num_coastlines = min(n, num_coastlines);
endi = 0;
for i = 1:num_coastlines

  % extract the ith coastline
  starti = endi+1;
  endi = coastline_breaks(i);
  
  % extract the ith coast, decimate it
%  disp( sprintf ( 'getting coastline, decimating by %i', ...
%                   decimate_coastline_factor ) );
  coasti = polar_coast(starti:decimate_coastline_factor:endi-1,:);
  
  n = length(coasti(:,1));
  coasti(n+1,:) = coasti(1,:);

  
  % check to see if the coastline intersects the grid boundary.  If not,
  % don't bother checking for whether or not any of the grid cells are
  % contained by the coastline.  This is sort of the inverse problem, but
  % this way if there is a small island just outside of the grid, we don't
  % have to check every friggin' grid cell to see if it's inside the island.
  inside_out_coast = insider ( coasti(:,1), coasti(:,2), ...
                               pgrid_bndry(:,1), pgrid_bndry(:,2) );
  if ( find(inside_out_coast) )

    % need to do each column of the grid.  Can't do it all at once because of
    % memory problems.
    for j = 1:nr
  
      inside_or_out_indices = insider ( RV(:,j), THETAV(:,j), ...
	                                coasti(:,1), coasti(:,2) );
      inside = find ( inside_or_out_indices );
      polar_grid(inside,j) = -Inf * ones(size(inside'));
      xy_grid(inside,j) = -Inf * ones(size(inside'));
      x_grid(inside,j) = -Inf * ones(size(inside'));
      y_grid(inside,j) = -Inf * ones(size(inside'));
  
    end %for j

  end %if
  

  
end  % for i


mesh(x_grid,y_grid,xy_grid);
set(gca,'view',[0 90]);
set(gca, 'aspectratio', [1 cos(43*pi/180)] );



% Set all data cells that are NOT -Inf to Inf.  These are the cells to grid.
% THIS DOESN'T MEAN THAT -INFs ARE CHANGED INTO INFs!  Go back and read it
% again. 
ifind = find(polar_grid ~= -Inf);
polar_grid(ifind) = Inf * ones(size(ifind));


% save grid in *.mat file
eval ( sprintf ( 'save %s polar_grid x_grid y_grid c_lon c_lat r1 r2 theta1 theta2 dr dtheta', out_matfile ) );


return;
