program TM2003_test
c
c   Coded:      N. A. Tsyganenko
c   Date:       Nov.28, 2017

 1    print *, 
     *' Enter distance RHO from Z-axis (Re) and angle PHI from midnight 
     *(radians)'
      read *, RHO,PHI
      print*,' Enter solar wind speed (km/s) and proton density (cc^-1)'
      read *, Vsw,Densw
      print *, ' Enter IMF By and Bz (nT) '
      read *, ByIMF,BzIMF

      T  =Plasmasheet_T (RHO,PHI,Densw,Vsw,BzIMF)
      Den=Plasmasheet_Den (RHO,PHI,Densw,Vsw,BzIMF)
      P  =Plasmasheet_P (RHO,PHI,Densw,Vsw,ByIMF,BzIMF)

      print *,''
      print *, '  Temperature = ',T,' (keV)'
      print *, '  Density = ',Den,' (cc^-1)'
      print *, '  Pressure = ',P,' (nPa)'
      print *,''

      goto 1
      end
c
c==========================================================================================
c
      FUNCTION Plasmasheet_T (RHO,PHI,Densw,Vsw,BzIMF)
c
c   Input: RHO = SQRT(XGSM**2+YGSM**2) (all coordinates in Earth's radii)
c          PHI = ATAN2 (YGSM,-ZGSM) (azimuthal angle measured from midnight meridian toward dusk, radians)
c          Densw:  solar wind proton density in cc^-1
c          Vsw:     solar wind speed in km/s
c          BzIMF:  IMF Bz (GSM, in nT)
c
c   Output: plasma sheet proton temperature in keV
c
c   Reference:  Tsyganenko, N. A., and T. Mukai, Tail plasma sheet models derived from Geotail particle data, 
c               J. Geophys. Res., 108(A3), 1136, doi:10.1029/2002JA009707, 2003.  
c
c   Coded:      N. A. Tsyganenko
c   Date:       Nov.28, 2017
c
      DIMENSION A(16)
c
      DATA A/1.67870,-0.16057,1.66874,4.81969,2.85451,-0.60195,-0.83629,
     * -2.49105,0.25679,0.22486,0.18870,-0.44579,-0.03306,-0.02412,
     * -2.68930,1.22152/
c
      X=-RHO*COS(PHI)
      Y= RHO*SIN(PHI)
      Z=0.
      Psw=1.94E-6*Densw*Vsw**2    ! ram pressure; coefficient 1.94e-6 corresponds to 4% of He++ ions in the solar wind

      CALL SHUEETAL (X,Y,Z,Psw,BzIMF,ID)   !  Shue et al (1998) magnetopause model

       IF (RHO.LT.10..OR.ABS(PHI).GT.1.5708.OR.ID.LT.0) THEN
         PRINT *, '  Coordinates outside the model validity region'
         PRINT *,''
         Plasmasheet_T=0.
         RETURN
       ENDIF

       VL=Vsw/500.
       Bz=BzIMF/5.
       BS=-Bz
       IF (BS.LT.0.) BS=0.
       BN= Bz
       IF (BN.LT.0.) BN=0.

       RO=RHO/10.-1.
c
       S=SIN(PHI)
       S2=S**2

       ARG=-(A(9)*VL**A(15)+A(10)*BN+A(11)*BS)*RO
       E1=EXP(ARG)

       ARG=-(A(12)*VL**A(16)+A(13)*BN+A(14)*BS)*RO
       E2=EXP(ARG)

       Plasmasheet_T=A(1)*VL+A(2)*BN+A(3)*BS+A(4)*E1+(A(5)*VL+A(6)*BN
     *  +A(7)*BS+A(8)*E2)*S2
c
      RETURN
      END
c
c==========================================================================================
c
      FUNCTION Plasmasheet_Den (RHO,PHI,Densw,Vsw,BzIMF)
c
c   Input: RHO = SQRT(XGSM**2+YGSM**2) (all coordinates in Earth's radii)
c          PHI = ATAN2 (YGSM,-ZGSM) (azimuthal angle measured from midnight meridian toward dusk, radians)
c          Densw:  solar wind proton density in cc^-1
c          Vsw:    solar wind speed in km/s
c          BzIMF:  IMF Bz (GSM, in nT)
c
c   Output: plasma sheet proton density in cc^-1
c
c   Reference:  Tsyganenko, N. A., and T. Mukai, Tail plasma sheet models derived from Geotail particle data, 
c               J. Geophys. Res., 108(A3), 1136, doi:10.1029/2002JA009707, 2003.  
c
c   Coded:      N. A. Tsyganenko
c   Date:       Nov.28, 2017
c
      DIMENSION A(11)
      DATA A/-0.15930,0.60805,0.50555,0.07959,0.27462,0.03611,-0.03419,
     *  -0.79347,1.16222,0.47559,0.71166/
c
      X=-RHO*COS(PHI)
      Y= RHO*SIN(PHI)
      Z=0.
      Psw=1.94E-6*Densw*Vsw**2   ! ram pressure; coefficient 1.94e-6 corresponds to 4% of He++ ions in the solar wind

      CALL SHUEETAL (X,Y,Z,Psw,BzIMF,ID)    !  Shue et al (1998) magnetopause model

       IF (RHO.LT.10..OR.ABS(PHI).GT.1.5708.OR.ID.LT.0) THEN
         PRINT *, '  Coordinates outside the model validity region'
         PRINT *,''
         Plasmasheet_Den=0.
         RETURN
       ENDIF

      RHS=RHO/10.           
      DNS=Densw/10.
      BNS=BzIMF/5.
      IF (BNS.LT.0.) BNS=0.
      BSS=-BzIMF/5. *(Vsw/500.)
      IF (BSS.LT.0.) BSS=0.
      Plasmasheet_Den=(A(1)+A(2)*DNS**A(10)+A(3)*BNS+A(4)*BSS)*RHS**A(8)
     * +(A(5)*DNS**A(11)+A(6)*BNS+A(7)*BSS)*RHS**A(9)*SIN(PHI)**2
c
      RETURN
      END
c
c============================================================================
c
      FUNCTION Plasmasheet_P (RHO,PHI,Densw,Vsw,ByIMF,BzIMF)
c
c   Input: RHO = SQRT(XGSM**2+YGSM**2) (all coordinates in Earth's radii)
c          PHI = ATAN2 (YGSM,-ZGSM) (azimuthal angle measured from midnight meridian toward dusk, radians)
c          Densw:  solar wind proton density in cc^-1
c          Vsw:    solar wind speed in km/s
c          ByIMF:  IMF By (GSM, in nT)
c          BzIMF:  IMF Bz (GSM, in nT)
c
c   Output: plasma sheet proton pressure in nPa
c
c   Reference:  Tsyganenko, N. A., and T. Mukai, Tail plasma sheet models derived from Geotail particle data, 
c               J. Geophys. Res., 108(A3), 1136, doi:10.1029/2002JA009707, 2003.  
c
c   Coded:      N. A. Tsyganenko
c   Date:       Nov.28, 2017
c
      DIMENSION A(14)
      DATA A/0.05701,0.52401,0.09075,0.52720,0.07814,-4.42226,-1.53325,
     *  -1.21666,2.53964,0.31988,0.75433,1.04862,-0.07390,1.01540/

      X=-RHO*COS(PHI)
      Y= RHO*SIN(PHI)
      Z=0.
      Psw=1.94E-6*Densw*Vsw**2   ! ram pressure; coefficient 1.94e-6 corresponds to 4% of He++ ions in the solar wind

      CALL SHUEETAL (X,Y,Z,Psw,BzIMF,ID)   !  Shue et al (1998) magnetopause model

       IF (RHO.LT.10..OR.ABS(PHI).GT.1.5708.OR.ID.LT.0) THEN
         PRINT *, '  Coordinates outside the model validity region'
         PRINT *,''
         Plasmasheet_p=0.
         RETURN
       ENDIF

      RHO10=RHO/10.D0

      S=SIN(PHI)
      S2=S*S

      Clock=ATAN2(ByIMF,BzIMF)
      IF (Clock.LT.0.) Clock=Clock+6.283185307

      FAC=SQRT((ByIMF**2+BzIMF**2)*SIN(Clock/2.))/5.
      Psw=1.94e-6*Densw*Vsw**2/3.

      ARG1=-A(9)*RHO10
      ARG2=-A(10)*RHO10

       Plasmasheet_P =
     *    A(1)*RHO10**A(6)+A(2)*Psw**A(11)*RHO10**A(7)
     *  + A(3)*FAC**A(12)*RHO10**A(8)
     *  +(A(4)*Psw**A(13)*EXP(ARG1)
     *  +A(5)*FAC**A(14)*EXP(ARG2))*S2

      RETURN
      END
C
C=========================================================================
c
        SUBROUTINE SHUEETAL (XGSM,YGSM,ZGSM,PDYN,BZIMF,ID)

        ID=-1
        R=SQRT(XGSM**2+YGSM**2+ZGSM**2)
        THETA=ACOS(XGSM/R)
        CT=COS(THETA)
        IF (CT.LT.-.99999) CT=-0.99999
        R0=(10.22+1.29*TANH(0.184*(BZIMF+8.14)))*PDYN**(-1./6.6)
        ALPHA=(0.58-0.007*BZIMF)*(1.+0.024*ALOG(PDYN))
        RM=R0*(2./(1.+CT))**ALPHA
        IF (R.LT.RM) ID=+1
        RETURN
        END
C
c--------------------------------------------------------------------------------