C======================================================================================== C SUBROUTINE TAG14_EQUAT_SHEET (XGSM,YGSM,ZGSM,PSI,PDYN,ByIMF,BzIMF, * XGSM_S,YGSM_S,ZGSM_S) C C INPUT: XGSM,YGSM,ZGSM: position of a point inside the magnetosphere C PSI: geodipole tilt angle in radians C PDYN: solar wind ram pressure in nanoPascals c ByIMF, BzIMF: Y- and Z- GSM components of the IMF in nT c (averaged over previous 30 minutes) C c OUTPUT: XGSM_S,YGSM_S,ZGSM_S: GSM coordinates of a point of the TAG14 equatorial sheet, c located at the same geocentric distance R=sqrt(XGSM^2+YGSM^2+ZGSM^2) c and lying in the same GSM meridian plane as the original point {XGSM,YGSM,ZGSM} c c Author: N. A. Tsyganenko, Dec.18, 2014 c Reference: Tsyganenko, N. A., V. A. Andreeva, and E. I. Gordeev, Internally and externally induced deformations c of the magnetospheric equatorial current as inferred from spacecraft data, c Ann.Geophys., v.33, pp.1-11, doi:10.5194/angeo-33-1-2015, 2015. c IMPLICIT REAL*8 (A-H,O-Z) DATA PI/3.1415926539D0/,ERR/1.D-8/ C DATA RH0,RH1,RH2,RH3,RH4,RH5,T0,T1,A00,A01,A02,A10,A11,A12,ALPHA0, * DALPHA1,DALPHA2,DALPHA3,XAPPA,BETA0,BETA1 * /11.01595255D0,6.054799445D0,0.8376136340D0,-2.283371357D0, * -0.2505011095D0,-0.9648465280D0,0.2865647076D0,0.1847283161D0, * 2.909423481D0,-0.1585372900D0,0.5648121327D0,1.889466989D0, * 0.6086006001D-01,0.4853873874D0,7.128509839D0,4.867577218D0, * -0.2248219212D0,-0.1441988549D0,-0.2868279445D0,2.179888855D0, * 0.4035940369D0/ C XGSM1=XGSM YGSM1=YGSM ZGSM1=ZGSM R=DSQRT(XGSM1**2+YGSM1**2+ZGSM1**2) THETA=DACOS(ZGSM1/R) IF (XGSM1.EQ.0.D0.AND.YGSM1.EQ.0.D0) THEN PHI=0.D0 ELSE PHI=DATAN2(YGSM1,XGSM1) ENDIF BYFACT=BYIMF/5.D0 BZFACT=BZIMF/5.d0 C PDYN_0=2.D0 PFACT=(PDYN/PDYN_0)**XAPPA-1.D0 CPS=DCOS(PSI) SPS=DSIN(PSI) CP=DCOS(PHI) SP=DSIN(PHI) TH1=0.D0 ! initial bounding values of THETA_S_GSM TH2=PI C 1 TH=0.5D0*(TH1+TH2) XGSM1=R*DSIN(TH)*CP YGSM1=R*DSIN(TH)*SP ZGSM1=R*DCOS(TH) XSM=XGSM1*CPS-ZGSM1*SPS YSM=YGSM1 ZSM=XGSM1*SPS+ZGSM1*CPS C RHO=DSQRT(XSM**2+YSM**2) C IF (DABS(RHO).LT.1.D-6) THEN COSPHI=1.D0 SINPHI=0.D0 ELSE COSPHI=XSM/RHO SINPHI=YSM/RHO ENDIF ALPHA=ALPHA0+DALPHA1*COSPHI+DALPHA2*PFACT+DALPHA3*BZFACT BETA=BETA0+BETA1*BZFACT C RH=RH0+RH1*PFACT+RH2*BZFACT+(RH3+RH4*PFACT+RH5*BZFACT)*COSPHI T=T0+T1*PFACT A0=A00+A01*PFACT+A02*BZFACT A1=A10+A11*PFACT+A12*BZFACT C F=1.D0-(1.D0+(RHO/RH)**ALPHA)**(1.D0/ALPHA) ! AXISYMMETRIC FACTOR (BOWL) G=A0+A1*COSPHI ! LONGITUDINAL ASYMMETRY FACTOR F1=T*BYFACT*(RHO/10.D0)**BETA*SINPHI ! TWISTING TERM ZSM_SHEET=RH*DTAN(PSI)*F*G +F1 FF=ZSM-ZSM_SHEET C PRINT *, TH1,TH2 C IF (FF.GT.0.D0) THEN TH1=TH ELSE TH2=TH ENDIF IF (DABS(TH1-TH2).GT.ERR) GOTO 1 C THETA_S_GSM=0.5D0*(TH1+TH2) C XGSM_S=R*DSIN(THETA_S_GSM)*CP YGSM_S=R*DSIN(THETA_S_GSM)*SP ZGSM_S=R*DCOS(THETA_S_GSM) C RETURN END C C=========================================================================================