function [LANaE, LANaW, LANdE, LANdW] = ...
                nscat50Q(CO, minlat, maxlat, minlon, maxlon)
%---------------------------------------------------------------------------
%
% Copyright 1997 University of Rhode Island.
%
% Compute the ascending and descending eastern and western boundries of the
% equatorial crossing for the given minimum and maximum latitude values. Use
% a `cannonical orbit' to perform these calculations. 
%
% Based on an algorithm by Peter Cornillon. This code is from an
% implementation by John Ryan. jhrg 11/19/97.
%  
% $Id: nsc50EQ.m,v 1.1 1999/05/13 02:55:20 dbyrne Exp $
% 
% a completely new version based on the older versions.
% begun 12 december 1998 by paul hemenway
%---------------------------------------------------------------------------

% The answer is:  
% We want the ascending equator crossing longitudes 
% Longitude of the ascending node==LAN) as follows:
%
% LANaW   LAN of the ascending pass whose swath just clips the SW corner.
% LANaE - LAN of the ascending pass whose swath just clips the NE corner.
% LANdW - LAN of the descending pass whose swath just clips the NW corner.
% LANdE - LAN of the descending pass whose swath just clips the SE corner.
%
% So, we want the ascending orbits with longitudes: LANaW<=long(asc)<=LANaE
%                     and
%     we want the descending orbits with longitudes LANdW<=long(des)<=LANdE.


% start with the latitude and longitude ranges of the requested dataset.
% (minlat, maxlat, minlon, maxlon)

% by convention, we want ALL longitudes between -360 and +720 degrees.
% by convention, we want the minimum longitude between 0 and 359.99 degrees.

if minlon<0; minlon=minlon + 360; end;
if maxlon<minlon; maxlon=maxlon + 360*(fix((minlon-maxlon)/360) + 1);
end;
%fprintf('     ********************************************************************\n')
%fprintf('     *                         maxlat: %f                        *\n',maxlat)
%disp('     *                                                                  *')
%disp('     *                                                                  *')
%disp('     *                                                                  *')
%disp('     *                                                                  *')
%disp('     *                                                                  *')
%disp('     *                                                                  *')
%disp('     *                                                                  *')
%disp('     *                                                                  *')
%fprintf('     *minlon: %f                              maxlon: %f*\n',minlon,maxlon)
%disp('     *                                                                  *')
%disp('     *                                                                  *')
%disp('     *                                                                  *')
%disp('     *                                                                  *')
%disp('     *                                                                  *')
%disp('     *                                                                  *')
%disp('     *                                                                  *')
%fprintf('     *                          minlat: %f                        *\n',minlat)
%%%%%%fprintf('     ********************************************************************\n')


% Load canonical orbit. This matrix contains the canonic orbit variables: 
%           
% CO_eqX1, CO_eqX2: The longitude of the ascending and descending equatorial
%           crossings.
% LD_AD: The longitude difference between the ascending and descending
%           equatorial crossings of the canonical orbit. 
%
% LON_1,LAT_1         are the coordinates of the SOUTHERN SWATH EDGE.
% LON_NADIR,LAT_NADIR are the coordinates of the SUB-SATELLITE TRACK.
% LON_24,LAT_24       are the coordinates of the NORTHERN SWATH EDGE.

%  Offset these values to find the equatorial crossing of an actual orbit.

load(CO);

% Find [ascending and descending] first-equatorial-crossing bounds for 
%     data to exist within user specified latitude,longitude limits.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Ascending Pass %%%%%%%%%%%%%%%%%%%%%%%%%%%

% assume the Canonical orbit given by "CO" as of 12 December 1998.
%fprintf('start Southwest Corner\n\n');

ascnd = [1:410];

% SOUTHERN-WESTERN CORNER (remember that the satellite is retrograde:
%
%           True Orbit    minlon                 Canonic Orbit
%             1  n  24       *                      n  24
%               1  n  24     *                        n  24
%                 1  n  24   *                          n  24
%                  1  n  24  *                           n  24
%            ^      1  n  24 *                            n  24
%           /|\     1  N  24 *   -------------------------N--24------lat=0
%            |      1  n  24 *                            n  24
%            |       1  n  24***************** minlat*****|n  24 ***********
%                     1  n  24                            | n  24
%                      1  n  24                           |  n  24
%                       1  n  24                          |   n  24
%                         1  n  24                        |   | n  24
%                          I                              |   |
%                          I                              |   |
%                          I                              |   |
%                          I                              |   LON_24(idx(1))
% LANaW =  {minlon -  [LON_24(idx(1)) - CO_eqX1]}          CO_eqX1
%
%  so we want ASCENDING orbits with longitudes GREATER than LANaW.

% What is the index of the scan of the northern swath edge at the lower
%          latitude limit?
[SL,idx] = sort(abs(LAT_24(ascnd) - minlat));

% Calculate longitude difference between [lon_nadir at equator] and
% [lon_24 at minlat]

LD1 =  LON_24(idx(1)) - CO_eqX1;

%fprintf('lon(northern edge @ minlat) - lon(ascend nadir-track @ equator):\n ')
%fprintf(' LD1  =  LLON_24(idx(1)) - CO_eqX1\n')
%fprintf(' %f = %f  -  %f\n\n',LD1,LON_24(idx(1)),CO_eqX1) 
 
%  Add this to minlon to get the south-western boundary for the ascending orbit.
%     LANaW is "unnormalized" at this point.

LANaW = minlon - LD1;

%fprintf('Lon. Ascend node of swath at south-west corner: \n')
%fprintf('        LANaW = minlon -LD1\n')
%fprintf('        %f = %f - %f\n\n',LANaW,minlon,LD1)

%disp(' end of southwestern corner (ascending-west)')
%disp(' ')
%disp('***********************************************')
%---------------------------------------------------------------------------
% NORTH-EASTERN CORNER
%disp(' start North-east corner')
%disp(' ')
% 
%
%          Canonic Orbit    minlon                True Orbit
%             1  n  24                              1  n 
%               1  n  24                              1  n 
%   **************1**n**24 ************ maxlat *********1  n 
%                 |1  n  24                            * 1  n 
%            ^    | 1  n  24                           *  1  n 
%           /|\   | 1  N  24     ----------------------*--1--n ------lat=0
%            |    | 1  n  24                           *  1  n 
%            |    |  1 |n  24                          *   1 |n             
%                 |   1| n  24                         *    1| n 
%                 |    1  n  24                        *     1  n 
%                 |    |1  n  24                       *     |1  n 
%                 |    |  1  n  24                     *     |  1  n 
%                 |    |                               *     | 
%                 |    |                               *     | 
%                 |    CO_eqX1                               |
%                 LON_1(idx(1))                      maxlon  |
%                                                            |
%                               LANaE =  {maxlon + [CO_eqX1 - LON_1(idx(1)) ] } 
%
%  so we want ASCENDING orbits with longitudes LESS than LANaE.

% What is the index of the scan of the southern swath edge at the northern
%          latitude limit?
[SL,idx] = sort(abs(LAT_1(ascnd) - maxlat));

% Calculate longitude difference between [lon_nadir at equator] and
% [lon_1 at maxlat]

LD2 = CO_eqX1 - LON_1(idx(1));

%fprintf('lon(southern edge @ maxlat) - lon(ascend nadir-track @ equator):\n ')
%fprintf(' LD2 = CO_eqX1 - LON_1(idx(1))')
%fprintf(' %f = %f  -  %f\n\n',LD2,LON_1(idx(1)),CO_eqX1) 
 
 
%  Add this to maxlon to get the north-eastern boundary for the ascending orbit.
%     LANaE is "unnormalized" at this point.

LANaE = maxlon + LD2 ;

%fprintf('Lon. Ascend node of swath at North-east corner: \n')
%fprintf('        LANaE = maxlon + LD1\n')
%fprintf('        %f = %f + %f\n\n',LANaE,maxlon,LD2)

%disp(' end of north-east corner (ascending orbit)' )
%disp(' ')
%disp('***********************************************')
%---------------------------------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%% Descending Pass %%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
descnd = [411:820];

% NORTH-WESTERN CORNER
%disp(' start North-west Corner')

[SL,idx] = sort(abs(LAT_1(descnd) - maxlat));

% Calculate longitude difference between [lon_nadir at equator] and
% [lon_1 at maxlat]
%
%  note that the index is offset from zero for the descending orbit.

LD3 = LON_1( idx(1)+descnd(1)-1 )  -  CO_eqX2;

%fprintf('lon(southren edge @ maxlat) - lon(descend nadir-track @ equator):\n ')
%fprintf(' LD3  =  LON_1(idx(1)) - CO_eqX2\n')
%fprintf(' %f = %f  -  %f\n\n',LD3,LON_1(idx(1)),CO_eqX2) 
 

%  Subtract from minlon to get western longitude boundary for the 
%                descending orbit, and correct to get the ascending node.
%  LANdW is unnormalized at this point.

LANdW = minlon - LD3 + (CO_eqX1 - CO_eqX2);

%fprintf('Lon. Ascend node of swath at North-west corner: \n')
%fprintf('        LANdW = minlon - LD3 + (CO_eqX1 - CO_eqX2)\n')
%fprintf('        %f = %f - %f + %f )\n\n',LANdW,minlon,LD3,...
%               (CO_eqX1 - CO_eqX2))


%disp(' end of North-west corner (descending orbit)' )
%disp(' ')
%disp('***********************************************')

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

% SOUTH_EASTERN CORNER

[SL,idx] = sort(abs(LAT_24(descnd) - minlat));

% Calculate longitude difference between [lon_nadir at equator] and
% [lon_24 at maxlat]
%
%  note that the index is offset from zero for the descending orbit.

LD4 = CO_eqX2 - LON_24(idx(1) + descnd(1) - 1) ;

%  Add this to maxlon to get eastern longitude boundary for the
%             descending orbit, and correct to get the ascending node.

LANdE = maxlon + LD4 + (CO_eqX1 - CO_eqX2);

%disp(' end of south-east corner (descending orbit)' )
%disp(' ')
%disp('***********************************************')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% rectify the eastern boundaries to be greater than the western boundaries.
% since the western boundaries are bound to be between 0-360,
%  then the eastern boundaries will be between 0-720.

% rectify LANaE and LANdE to the maximum (eastern) longitude:
  LANaE = xrectang(maxlon,LANaE,180);
  LANdE = xrectang(maxlon+LD_AD,LANdE,180);

% rectify LANaW and LANdW to the minimum (western) longitude:
  LANaW = xrectang(minlon,LANaW,180);
  LANdW = xrectang(minlon+LD_AD,LANdW,180);

%disp(' end of rectification')
%disp(' next return with return to the calling program')
%disp(' ')NaE = xrectang(maxlon,LANaE,180);


%NOW REALLY CONFUSE THE ISSUE:
%if the Western edge (aW or dW) is more than 360, subtract 360 from BOTH.
if (LANaW>360)
   LANaW=LANaW-360;
   LANaE=LANaE-360;
end
if (LANdW>360)
   LANdW=LANdW-360;
   LANdE=LANdE-360;
end


%if the Eastern edge (aE or dE) is less than zero, add 360 to BOTH quantities.

if LANaE<0
    LANaE=LANaE + 360.
    LANaW=LANaW + 360.
end
if LANdE<0
    LANdE=LANdE + 360.
    LANdW=LANdW + 360.
end

return;


