function [out1,out2,savepts] = hatano(mstruct,in1,in2,object,direction,savepts)

%HATANO  Hatano Asymmetrical Equal Area Pseudocylindrical Projection
%
%  This is an equal area projection.  Scale is true along 40 deg,
%  42 min N and 38 deg, 27 min S, and is constant along any parallel
%  but generally not between pairs of parallels equidistant
%  from the Equator.  It is free of distortion along the central
%  meridian at 40 deg, 42 min N and 38 deg, 27 min S.  This
%  projection is not conformal or equidistant.
%
%  This projection was presented by Masataka Hatano in 1972.

%  Copyright 1996-2003 The MathWorks, Inc.
% $Revision: 1.9.4.1 $
%  Written by:  E. Byrns, E. Brown


if nargin == 1                  %  Set the default structure entries
	  mstruct.trimlat = angledim([ -90  90],'degrees',mstruct.angleunits);
      mstruct.trimlon = angledim([-180 180],'degrees',mstruct.angleunits);
	  mstruct.mapparallels = angledim([3827 4042],'dms',mstruct.angleunits);
      mstruct.nparallels   = 0;
	  mstruct.fixedorient  = [];
	  out1 = mstruct;          return
elseif nargin ~= 5 & nargin ~= 6
      error('Incorrect number of arguments')
end


units  = mstruct.angleunits;
aspect = mstruct.aspect;
radius = rsphere('authalic',mstruct.geoid);
origin  = angledim(mstruct.origin,units,'radians');
trimlat = angledim(mstruct.flatlimit,units,'radians');
trimlon = angledim(mstruct.flonlimit,units,'radians');
scalefactor = mstruct.scalefactor;
falseeasting = mstruct.falseeasting;
falsenorthing = mstruct.falsenorthing;

%  Adjust the origin latitude to the auxiliary sphere

origin(1) = geod2aut(origin(1),mstruct.geoid,'radians');
trimlat   = geod2aut(trimlat  ,mstruct.geoid,'radians');

switch direction
case 'forward'

     lat  = angledim(in1,units,'radians');
     lat  = geod2aut(lat,mstruct.geoid,'radians');
     long = angledim(in2,units,'radians');

     [lat,long] = rotatem(lat,long,origin,direction);   %  Rotate to new origin
     [lat,long,clipped] = clipdata(lat,long,object);    %  Clip at the date line
     [lat,long,trimmed] = trimdata(lat,trimlat,long,trimlon,object);

%  Construct the structure of altered points

     savepts.trimmed = trimmed;
     savepts.clipped = clipped;

%  Determine the location of positive and negative entries
%  Necessary since this projection treats northern and southern
%  latitudes differently.

     indx1 = find(lat >= 0);    indx2 = find(lat < 0);

%  Pick up NaN place holders in the output data.

	 x = long;  y = lat;

%  Projection transformation
%  Convergence is selected to ensure successful testing of forward
%  and inverse points (the hard point set) using TMAPCALC

     convergence = 1E-10;     maxsteps = 100;
     steps = 1;               thetanew = lat;
	 converged = 0;

     while ~converged & steps <= maxsteps
          steps = steps + 1;
          thetaold = thetanew;
		  delnorth = -(thetaold(indx1) + sin(thetaold(indx1)) - ...
		               2.67595*sin(lat(indx1))) ./ ...
		                       (1 + cos(thetaold(indx1)));

		  delsouth = -(thetaold(indx2) + sin(thetaold(indx2)) - ...
		               2.43763*sin(lat(indx2))) ./ ...
		                       (1 + cos(thetaold(indx2)));

          if max(abs(delnorth(:))) <= convergence & ...
		     max(abs(delsouth(:))) <= convergence
		            converged = 1;
		  else
		            thetanew(indx1) = thetaold(indx1) + delnorth;
		            thetanew(indx2) = thetaold(indx2) + delsouth;
	      end
     end
     thetanew = thetanew / 2;

     x(indx1) = 0.85 * radius * long(indx1) .* cos(thetanew(indx1));
     x(indx2) = 0.85 * radius * long(indx2) .* cos(thetanew(indx2));
     y(indx1) = 1.75859 * radius * sin(thetanew(indx1));
     y(indx2) = 1.93052 * radius * sin(thetanew(indx2));

%  Apply scale factor, false easting, northing

	x = x*scalefactor+falseeasting;
	y = y*scalefactor+falsenorthing;

%  Set the output variables

     switch  aspect
	    case 'normal',         out1 = x;      out2 = y;
	    case 'transverse',	   out1 = y;      out2 = -x;
        otherwise,             error('Unrecognized aspect string')
     end


case 'inverse'

     switch  aspect
	    case 'normal',         x = in1;    y = in2;
	    case 'transverse',	   x = -in2;   y = in1;
        otherwise,             error('Unrecognized aspect string')
     end

%  Apply scale factor, false easting, northing

	x = (x-falseeasting)/scalefactor;
	y = (y-falsenorthing)/scalefactor;

%  Pick up NaN place holders in the output data.

	 long = x;  lat = y;

%  Determine the location of positive and negative entries
%  Necessary since this projection treats northern and southern
%  latitudes differently.

     indx1 = find(y >= 0);    indx2 = find(y < 0);

%  Northern latitudes

     theta = asin(y(indx1) ./ (1.75859 * radius));
	 lat(indx1)  = asin((2*theta +sin(2*theta))/2.67595);
	 long(indx1) = x(indx1) ./(0.85*radius*cos(theta));

%  Southern Latitudes

     theta = asin(y(indx2) ./ (1.93052 * radius));
	 lat(indx2)  = asin((2*theta +sin(2*theta))/2.43763);
	 long(indx2) = x(indx2) ./(0.85*radius*cos(theta));

%  Undo trims and clips

     [lat,long] = undotrim(lat,long,savepts.trimmed,object);
     [lat,long] = undoclip(lat,long,savepts.clipped,object);

%  Now rotate data back to baseline coordinate system

     [lat,long] = rotatem(lat,long,origin,direction);
	 lat = aut2geod(lat,mstruct.geoid,'radians');

%  Store the output data and convert to the proper units

     out1 = angledim(lat, 'radians', units);
     out2 = angledim(long,'radians', units);

otherwise
     error('Unrecognized direction string')
end

%  Some operations on NaNs produce NaN + NaNi.  However operations
%  outside the map may product complex results and we don't want
%  to destroy this indicator.

indx = find(isnan(out1) | isnan(out2));
out1(indx) = NaN;   out2(indx) = NaN;