%C$*****************************************************************************
%C$                                                                            *
%C$                                                                            *
%C$    DESCRIPTION:                                                            *
%C$                                                                            *
%C$    This subroutine converts from geodetic latitude and longitude to Polar  *
%C$    Stereographic (X,Y) coordinates for the polar regions.  The equations   *
%C$    are from Snyder, J. P., 1982,  Map Projections Used by the U.S.         *
%C$    Geological Survey, Geological Survey Bulletin 1532, U.S. Government     *
%C$    Printing Office.  See JPL Technical Memorandum 3349-85-101 for further  *
%C$    details.                                                                *
%C$                                                                            *
%C$                                                                            *
%C$    ARGUMENTS:                                                              *
%C$                                                                            *
%C$    Variable    Description                              *
%C$                                                         *
%C$    ALAT        Geodetic Latitude (degrees, +90 to -90)  *
%C$    ALONG       Geodetic Longitude (degrees, 0 to 360)   *
%C$    X           Polar Stereographic X Coordinate (km)    *
%C$    Y           Polar Stereographic Y Coordinate (km)    *
%      ii
%      jj
%C$                                                                            *
%C$                                                                            *
%C$                  Written by C. S. Morris - April 29, 1985                  *
%C$                  Revised by C. S. Morris - December 11, 1985               *
%C$                                                                     	      *
%C$                  Revised by V. J. Troisi - January 1990                    *
%C$                  SGN - provides hemisphere dependency (+/- 1)              *
%C$		    Revised by Xiaoming Li - October 1996                     *
%C$		    Corrected equation for RHO                                *
%
%      Modified to Matlab version and vectorized by Yingtao Ma; 120918
%C$*****************************************************************************

function [ii,jj] = mapll(ALAT,ALON)

      SLAT   = 70; %Standard latitude for the SSM/I grids is 70 degrees.
      CellSz = 25.0; NumY = [448,332];
                   %Size of the SSM/I grid ( 12.5 km, 25.0 km)
                   %Number of lines in an SSM/I grid for CellSiz=25.0km; [North, South]
                   %CellSz = 12.5; NumY = [896,664];
      XYDist = [3850.0, 5350.0; 3950.0, 3950.0];
                   %Distance from the origin of the grid in the cartesian plane.
                   %The x-y coordinates for the edge of the lower left pixel
                   %is (3850.0, 5350.0) for the northern grids and 
                   %(3950.0, 3950.0) for the southern grids.
      RE   = 6378.273; %Radius of the earth in kilometers.
      E2   = 0.006693883; %Eccentricity squared
      E    = sqrt(E2); %Eccentricity of the Hughes ellipsoid


      %-- Only defined for lat<90.
      sb=abs(ALAT)>=90;
      ALAT(sb) = NaN;
      ALON(sb) = NaN;

      %-- sign for polar regions (1=North , -1=South)
      SGN = ones(size(ALAT));
      SGN(ALAT<=0) = -1; %lat==0, set to south
      
      %-- Meridian offset for the SSM/I grids (0 degrees for
      %   the South Polar grids; 45 degrees for the North Polar grids. 
      ALON(SGN>0) = ALON(SGN>0) + 45;
      ALON(SGN<0) = ALON(SGN<0) + 0;

      %-- All latitude and longitude need to be positive numbers
      ALAT = abs(ALAT);
      %ALON = abs(ALON); %values are already from 0..360
      
      %-- Degree to radians
      ALAT = ALAT*pi/180;
      ALON = ALON*pi/180;
      

      %---
      T = tan(pi/4-ALAT/2) ./...
          ( (1-E*sin(ALAT)) ./ (1+E*sin(ALAT)) ).^(E/2);
      
      if abs(90-SLAT) < 1e-5
         RHO = 2*RE*T ./ ( (1+E).^(1+E) .* (1-E).^(1-E) ).^0.5;
      else
         SL = SLAT*pi/180;
         TC = tan(pi/4-SL/2) / ( (1-E*sin(SL)) / (1+E*sin(SL)) )^(E/2);
         MC = cos(SL) / sqrt(1-E2*(sin(SL)^2));
         RHO = RE*MC*T/TC;
      end

      Y = -RHO .* SGN .* cos(SGN.*ALON);
      X =  RHO .* SGN .* sin(SGN.*ALON);


      %--Convert x and y distances from origin to I,J pair (ii,jj)
      %
      ii = nan(size(X));
      jj = nan(size(Y));

      for ihem=1:2  %=1 north; =2 south.
         if ihem==1
            sb = SGN>0;
         else
            sb = SGN<0;
         end
         ii(sb) = round( (X(sb)+XYDist(ihem,1)-CellSz/2) / CellSz ) + 1;
         jj(sb) = round( (Y(sb)+XYDist(ihem,2)-CellSz/2) / CellSz ) + 1;
         jj(sb) = NumY(ihem) - (jj(sb)-1);%Flip grid orientation in the 'Y' direction
      end
      
      %ii(isnan(ii) = 0;
      %jj(isnan(jj) = 0;
   
%keyboard
end %---------------------------------------------------
