function [XMGNP,YMGNP,ZMGNP,DIST,ID] = GEOPACK_SHUETAL_MGNP(XN_PD,VEL,BZIMF,XGSM,YGSM,ZGSM);
% function [XMGNP,YMGNP,ZMGNP,DIST,ID] = GEOPACK_SHUETAL_MGNP(XN_PD,VEL,BZIMF,XGSM,YGSM,ZGSM);
%       SUBROUTINE SHUETAL_MGNP(XN_PD,VEL,BZIMF,XGSM,YGSM,ZGSM,
%      *  XMGNP,YMGNP,ZMGNP,DIST,ID)
% C
% C  FOR ANY POINT OF SPACE WITH COORDINATES (XGSM,YGSM,ZGSM) AND SPECIFIED CONDITIONS
% C  IN THE INCOMING SOLAR WIND, THIS SUBROUTINE:
% C
% C (1) DETERMINES IF THE POINT (XGSM,YGSM,ZGSM) LIES INSIDE OR OUTSIDE THE
% C      MODEL MAGNETOPAUSE OF SHUE ET AL. (JGR-A, V.103, P. 17691, 1998).
% C
% C (2) CALCULATES THE GSM POSITION OF A POINT {XMGNP,YMGNP,ZMGNP}, LYING AT THE MODEL
% C      MAGNETOPAUSE AND ASYMPTOTICALLY TENDING TO THE NEAREST BOUNDARY POINT WITH
% C      RESPECT TO THE OBSERVATION POINT {XGSM,YGSM,ZGSM}, AS IT APPROACHES THE MAGNETO-
% C      PAUSE.
% C
% C  INPUT: XN_PD - EITHER SOLAR WIND PROTON NUMBER DENSITY (PER C.C.) (IF VEL>0)
% C                    OR THE SOLAR WIND RAM PRESSURE IN NANOPASCALS   (IF VEL<0)
% C         BZIMF - IMF BZ IN NANOTESLAS
% C
% C         VEL - EITHER SOLAR WIND VELOCITY (KM/SEC)
% C                  OR ANY NEGATIVE NUMBER, WHICH INDICATES THAT XN_PD STANDS
% C                     FOR THE SOLAR WIND PRESSURE, RATHER THAN FOR THE DENSITY
% C
% C         XGSM,YGSM,ZGSM - GSM POSITION OF THE OBSERVATION POINT IN EARTH RADII
% C
% C  OUTPUT: XMGNP,YMGNP,ZMGNP - GSM POSITION OF THE BOUNDARY POINT
% C          DIST - DISTANCE (IN RE) BETWEEN THE OBSERVATION POINT (XGSM,YGSM,ZGSM)
% C                 AND THE MODEL NAGNETOPAUSE
% C          ID -  POSITION FLAG:  ID=+1 (-1) MEANS THAT THE OBSERVATION POINT
% C          LIES INSIDE (OUTSIDE) OF THE MODEL MAGNETOPAUSE, RESPECTIVELY.
% C
% C  OTHER SUBROUTINES USED: T96_MGNP
% C
% c          AUTHOR:  N.A. TSYGANENKO,
% C          DATE:    APRIL 4, 2003.
% C
if (VEL < 0.) ,
    PD=XN_PD;
else
    PD=1.94E-6*XN_PD*VEL^2;  %! PD IS THE SOLAR WIND DYNAMIC PRESSURE (IN nPa)
end

% c
% c  DEFINE THE ANGLE PHI, MEASURED DUSKWARD FROM THE NOON-MIDNIGHT MERIDIAN PLANE;
% C  IF THE OBSERVATION POINT LIES ON THE X AXIS, THE ANGLE PHI CANNOT BE UNIQUELY
% C  DEFINED, AND WE SET IT AT ZERO:
% c
if (YGSM ~= 0.) | (ZGSM ~= 0.) ,
    PHI=atan2(YGSM,ZGSM);
else
    PHI=0.;
end
% C
% C  FIRST, FIND OUT IF THE OBSERVATION POINT LIES INSIDE THE SHUE ET AL BDRY
% C  AND SET THE VALUE OF THE ID FLAG:
% C
ID=-1;
R0=(10.22+1.29*tanh(0.184*(BZIMF+8.14)))*PD^(-.15151515);
ALPHA=(0.58-0.007*BZIMF)*(1.+0.024*log(PD));
R=sqrt(XGSM^2+YGSM^2+ZGSM^2);
RM=R0*(2./(1.+XGSM/R))^ALPHA;
if (R <= RM), ID=+1; end
% C
% C  NOW, FIND THE CORRESPONDING T96 MAGNETOPAUSE POSITION, TO BE USED AS
% C  A STARTING APPROXIMATION IN THE SEARCH OF A CORRESPONDING SHUE ET AL.
% C  BOUNDARY POINT:
% C
[XMT96,YMT96,ZMT96,DIST,ID96] = GEOPACK_T96_MGNP(PD,-1.,XGSM,YGSM,ZGSM);
% C
RHO2=YMT96^2+ZMT96^2;
R=sqrt(RHO2+XMT96^2);
ST=sqrt(RHO2)/R;
CT=XMT96/R;
% C
% C  NOW, USE NEWTON'S ITERATIVE METHOD TO FIND THE NEAREST POINT AT THE
% C   SHUE ET AL.'S BOUNDARY:
% C
NIT=0;

while 1,
    T=atan2(ST,CT); % 1
    RM=R0*(2./(1.+CT))^ALPHA;
    
    F=R-RM;
    GRADF_R=1.;
    GRADF_T=-ALPHA/R*RM*ST/(1.+CT);
    GRADF=sqrt(GRADF_R^2+GRADF_T^2);
    
    DR=-F/GRADF^2;
    DT= DR/R*GRADF_T;
    
    R=R+DR;
    T=T+DT;
    ST=sin(T);
    CT=cos(T);
    
    DS=sqrt(DR^2+(R*DT)^2);
    
    NIT=NIT+1;
    
    if (NIT > 1000) ,
        disp('BOUNDARY POINT COULD NOT BE FOUND; ITERATIONS DO NOT CONVERGE');
    end
    
    %      IF (DS > 1.E-4) GOTO 1
    if ~ (DS > 1.E-4), break; end;
end % while 1

XMGNP=R*cos(T);
RHO=  R*sin(T);

YMGNP=RHO*sin(PHI);
ZMGNP=RHO*cos(PHI);

DIST=sqrt((XGSM-XMGNP)^2+(YGSM-YMGNP)^2+(ZGSM-ZMGNP)^2);

%       RETURN
%       END

% end of function SHUEETAL_MGNP
% C
% C=======================================================================================
% C