C***********************************************************************
      SUBROUTINE  KMAG(TIME,EPOCH,R,THETA,PHI,RLT,BY_IMF,BZ_IMF,
     ~  BRM,BTM,BPM,Dp)
C     R,  THETA AND PHI ARE RADIAL DISTANCE, COLATITUDE (IN RADIANS) AND LONGITUDE (IN RADIANS)
C     OF THE SPACECRAFT. THE PROGRAM RETURNS, LOCAL TIME (LT, IN HOURS) AND THE MODEL 
C     MAGNETIC FIELD IN SYSTEM 3 SPHERICAL COORDINATES.

      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION PVECIN(3), PVECOU(3)
      DIMENSION BVECIN(3), BVECOU(3)
      REAL*8, INTENT(IN) :: TIME
      REAL*8, INTENT(IN) :: R
      REAL*8, INTENT(IN) :: THETA
      REAL*8, INTENT(IN) :: PHI
      REAL*8, INTENT(OUT) :: RLT
      REAL*8, INTENT(IN) :: BY_IMF
      REAL*8, INTENT(IN) :: BZ_IMF
      REAL*8, INTENT(OUT) :: BRM
      REAL*8, INTENT(OUT) :: BTM
      REAL*8, INTENT(OUT) :: BPM
      REAL*8, INTENT(IN) :: Dp
      CHARACTER*5, INTENT(IN) :: EPOCH
      INTEGER is_inside_mp
      INTEGER ics_on
      REAL*8 Dp0, a1, a2, dK, dhD
      REAL*8 D, r0, r0or
      DIMENSION C(5)
C     NL is the number of current sheet modes
      INTEGER, PARAMETER :: NL = 5
      REAL*8 Brscsarr(NL), Btscsarr(NL), Bpscsarr(NL)
C     DIMENSION Brscsarr(5), Btscsarr(5), Bpscsarr(5)
      DATA dql /1.41/

C     Strength of tail modes
      ics_on = 1
C     NL = 5 !number of current sheet modes
      DATA C /
     ~   2.023453,  
     ~   1.500525, 
     ~   0.441146,  
     ~   0.286629,  
     ~   0.152069/ 
      
C     DATA C/
C     ~  -0.553754,
C     ~   0.714655,
C     ~   0.119748,
C     ~   1.237026/

      DATA Dp0/0.017/ !nominal dynamic pressure

      DATA a1/9.696/
      DATA a2/0.2351/
      DATA dK/0.71516/ !for selfsimilar mp

      DATA dhD/2.0/ !nominal half thickness of the cs

      r0 = a1 * Dp0**(-a2)

      r0or = r0 / (a1 * Dp**(-a2))

      D = dhD / r0or !Half thickness of the CS

C FIRST CALCULATE THE INTERNAL FIELD HIGHER THAN DIPOLE
      CALL  KRONIAN_HIGHER(3,R,THETA,PHI,BRI,BTI,BPI)

C NOW CALCULATE THE AXISYMMETRIC DIPOLE AND ITS SHIELD FIELD
      XS3 = R * DSIN(THETA) * DCOS(PHI)
      YS3 = R * DSIN(THETA) * DSIN(PHI)
      ZS3 = R * DCOS(THETA)

C ROTATE TRAJECTORY INTO DIS COORDINATES
      pvecin(1) = xs3
      pvecin(2) = ys3
      pvecin(3) = zs3 
      CALL KROT('S3C','DIS',pvecin,pvecou,time,Epoch)
      XDIS = pvecou(1)
      YDIS = pvecou(2)
      ZDIS = pvecou(3)

C CALCULATE THE MAPPED LOCATION
      CALL mapit(time,XDIS,YDIS,ZDIS,XMAP,YMAP,ZMAP,r0or,Epoch)

      RMAP = dsqrt(XMAP**2 + YMAP**2 + ZMAP**2)
      
      CALL checkIfInsideMpSATURN(RMAP,XMAP,is_inside_mp,Dp,a1,a2,dK)

      CALL dipole_shielded(XMAP, YMAP, ZMAP, BXD, BYD, BZD, r0or)
C     CALL dipole_shielded(XMAP, YMAP, ZMAP, BXD, BYD, BZD, r0or)
C     WARNING : XMAPP to XMAP and added YMAP to arguments

C NOW MAP IT AND ROTATE INTO SYSTEM III
      CALL Mapped_field(time,XMAP,YMAP,ZMAP,BXD,BYD,BZD,BXDMAP,BYDMAP,
     ~BZDMAP,r0or,Epoch)
      bvecin(1) = BXDMAP
      bvecin(2) = BYDMAP
      bvecin(3) = BZDMAP
      CALL KROT('DIS','S3C',bvecin,bvecou,time,Epoch)
      BXDS3 = bvecou(1)
      BYDS3 = bvecou(2)
      BZDS3 = bvecou(3)
      BRD =  BXDS3*DSIN(THETA)*DCOS(PHI) + BYDS3*DSIN(THETA)*
     .DSIN(PHI)+BZDS3*DCOS(THETA)
      BTD =  BXDS3*DCOS(THETA)*DCOS(PHI) + BYDS3*DCOS(THETA)*DSIN(PHI)-
     .BZDS3*DSIN(THETA)
      BPD = -BXDS3*DSIN(PHI)             + BYDS3*DCOS(PHI)
    
C NOW CALCULATE THE FIELD FROM CURRENT SHEET AND ITS SHIELDING
      CALL shielded_csheetfield(time,R,THETA,PHI,RLT,Brscs,Btscs,Bpscs,
     ~  NL,C,D,r0or,dql,Epoch,Brscsarr,Btscsarr,Bpscsarr)

C NOW GET THE INFLUENCE OF THE IMF
  
      CALL getIMF_penetration(BY_IMF,BZ_IMF,BY_p,BZ_p)
      bvecin(1) = 0.0
      bvecin(2) = BY_p
      bvecin(3) = BZ_p
      CALL KROT('KSM','S3C',bvecin,bvecou,time,Epoch)
      BX_ps3 = bvecou(1)
      BY_ps3 = bvecou(2)
      BZ_ps3 = bvecou(3)
  
      CALL CAR2SPH_MAG(BX_ps3,BY_ps3,BZ_ps3,BR_p,BT_p,BP_p,THETA,PHI)

C NOW ADD ALL THE FIELDS

      IF (is_inside_mp .eq. 1) THEN

          ! higher order internal field + shieded dipole + shielded current sheet field + IMF penetration
          IF (ics_on .eq. 1) THEN
              Brm = BRI + BRD + Brscs + Br_p
              Btm = BTI + BTD + Btscs + Bt_p
              Bpm = BPI + BPD + Bpscs + Bp_p
          ELSE
              Brm = BRI + BRD + Br_p
              Btm = BTI + BTD + Bt_p
              Bpm = BPI + BPD + Bp_p
          END IF

      !Brm=BRD + Brscs + Br_p
      !Btm=BTD + Btscs + Bt_p
      !Bpm=BPD + Bpscs + Bp_p

      !Brm=BRI
      !Btm=BTI
      !Bpm=BPI

      !Brm=Brscs 
      !Btm=Btscs 
      !Bpm=Bpscs 
      ELSE
          Brm = Br_p
          Btm = Bt_p
          Bpm = Bp_p
          
          !set current sheet field to zero if outside mp
          Brscs = 0.0
          Btscs = 0.0
          Bpscs = 0.0
          
          DO i = 1, NL
              Brscsarr(i) = 0.0
              Btscsarr(i) = 0.0
              Bpscsarr(i) = 0.0
          END DO

      END IF

      RETURN
      END SUBROUTINE KMAG
C***********************************************************************

      SUBROUTINE KRONIAN_HIGHER (NM,R,T,F,BR,BT,BF)

C     BASED ON THE SUBROUTINE IGRF WRITTEN BY N.A. TSYGANENKO (1979)
C     MODIFIED BY KRISHAN KHURANA, JULY, 1996. AND NOV. 2004.
C
C     CALCULATES COMPONENTS OF MAIN KRONIAN FIELD IN RIGHT HANDED SPHERICAL
C     COORD SYSTEM. BASED ON THE  SPHERICAL HARMONIC COEFFICIENTS GIVEN BY 
C     ACUNA ET AL. [1983] (Z3 MODEL,  V1+V2 DATA) 
C     MAXIMUM ORDER OF HARMONICS TAKEN INTO ACCOUNT (NOT MORE THAN 0RDER 3)
C
C
C     IT IS ASSUMED THAT THE TRAJECTORY IS IS IN RIGHT HANDED S III COORDINATES.
C     THE OUTPUT IS ALSO IN RTP (RH) COORDINATES.
C
C            INPUT:  NM (INTEGER)- MAXIMUM ORDER OF HARMONICS TAKEN
C                                  INTO ACCOUNT (NM.LE.12)
C
C                    R,T,F (REAL)- POSITION OF DESIRED FIELD VALUE IN
C                                  SPHERICAL JOVIGRAPHIC COORDINATE SYSTEM
C                                  (R IN PLANET RADII, COLATITUDE T AND
C                                  RH LONGITUDE F IN RADIANS)
C
C           OUTPUT: BR,BT,BF (REAL)- COMPONENTS OF THE INTERNAL PORTION
C                                    OF THE MAIN MAGNETIC FIELD IN
C                                    SPHERICAL S III COORD SYSTEM
C                                    (VALUES GIVEN IN GAMMA)
C
C***********************************************************************
      IMPLICIT REAL*8 (A-H,O-Z)
      INTEGER, INTENT(IN)  :: NM
      REAL*8,  INTENT(IN)  :: R,T,F
      REAL*8,  INTENT(OUT) :: BR, BT, BF
      LOGICAL BK, BM
      DIMENSION A(13), B(13), G(91), H(91), REC(91)
      DATA REC/ 91 * 1.0/

C     DATA G/0.,0.,24.,1449.,-179.,-35.,2012.,
C         ~       53.,92.,-40.,
C         ~       81*0.0/

C     DATA H/0.,0.,-0.9110,0.,43.5680,-44.7762,0.,
C         ~       140.9273,-92.8876,68.2759,
C         ~       81*0.0/

      ! newer coefficients 2012
      DATA G/0.,0.,0.,1586.,0.,0.,2374.,
     ~       0.,0.,0.,65.,0.,0.,0.,0.,185.,
     ~       75*0.0/

      DATA H/91*0.0/
C      1.Cao, H., C.T. Russell, U.R. Christensen, M.K. Dougherty, and M.E. Burton(2011), Saturn's very axisymmetric magnetic field: No detectable secular variation or tilt, Earth and Planet. Sci. Lett., 304, 22-28 
C      2.Cao, H., Russell, C.T., Christensen, U.R., Wicht, J., Dougherty, M.K., 2012. Saturn's High Degree Magnetic Moments: Evidence for a Unique Planetary Dynamo. Submitted to Icarus)

      DATA FIRSTI/0.0/
C     WRITE(1,'(5F15.3)')(G(I),I=1,10)
C     WRITE(1,'(5F14.3)')(H(I),I=1,10)
      IF (FIRSTI.EQ.0.0) GO TO 1
      GO TO 6
  1   FIRSTI=1.0
      G(1)=0.
      H(1)=0.
      KNM=15
      DO 2 N=1,13
      N2=2*N-1
      N2=N2*(N2-2)
      DO 2 M=1,N
      MN=N*(N-1)/2+M
  2   REC(MN)=DBLE((N-M)*(N+M-2))/DBLE(N2)
      S=1.
      DO 5 N=2,13
      MN=N*(N-1)/2+1
      S=S*DBLE(2*N-3)/DBLE(N-1)
      G(MN)=G(MN)*S
      H(MN)=H(MN)*S
      P=S
      DO 5 M=2,N
      AA=1.
      IF (M.EQ.2) AA=2.
      P=P*DSQRT(AA*DBLE(N-M+1)/DBLE(N+M-2))
      MNN=MN+M-1
      G(MNN)=G(MNN)*P
  5   H(MNN)=H(MNN)*P
  6   IF(KNM.EQ.NM) GOTO 61
      KNM=NM
      K=KNM+1
 61   PP=1./R
      P=PP
      DO 7 N=1,K
      P=P*PP
      A(N)=P
  7   B(N)=P*N
      P=1.
      D=0.
      BBR=0.
      BBT=0.
      BBF=0.
      U=T
      CF=DCOS(F)
      SF=DSIN(F)
      C=DCOS(U)
      S=DSIN(U)
      BK=(S.LT.1.D-5)
      DO 12 M=1,K
      BM=(M.EQ.1)
      IF(BM) GOTO 8
      MM=M-1
      W=X
      X=W*CF+Y*SF
      Y=Y*CF-W*SF
      GOTO 9
  8   X=0.
      Y=1.
  9   Q=P
      Z=D
      BI=0.
      P2=0.
      D2=0.
      DO 11 N=M,K
      AN=A(N)
      MN=N*(N-1)/2+M
      E=G(MN)
      HH=H(MN)
      W=E*Y+HH*X
      IF (DABS(P2).LT.1.D-38) P2=0.0
      IF (DABS(Q).LT.1.D-38) Q=0.0
      BBR=BBR+B(N)*W*Q
      BBT=BBT-AN*W*Z
      IF(BM) GOTO 10
      QQ=Q
      IF(BK) QQ=Z
      BI=BI+AN*(E*X-HH*Y)*QQ
  10  XK=REC(MN)
      DP=C*Z-S*Q-XK*D2
      PM=C*Q-XK*P2
      D2=Z
      P2=Q
      Z=DP
      Q=PM
  11  CONTINUE
      D=S*D+C*P
      P=S*P
      IF(BM) GOTO 12
      BI=BI*MM
      BBF=BBF+BI
  12  CONTINUE
      BR=BBR
      BT=BBT
      IF(BK) GOTO 13
      BF=BBF/S
      GOTO 14
  13  IF(C.LT.0.) BBF=-BBF
      BF=BBF
  14  CONTINUE
      RETURN
      END SUBROUTINE KRONIAN_HIGHER

C***********************************************************************

      SUBROUTINE dipole_shielded(x,y,z,Bx,By,Bz,r0or)
C           WRITTEN BY K. K. KHURANA     11/2004

C           INPUT: x,y,z input position
C           OUTPUT: mag .filed Bx,By,Bz at x,y,z
      IMPLICIT REAL*8(A-H,O-Z)
      REAL*8, INTENT(IN) :: X,Y,Z
      REAL*8, INTENT(IN) :: r0or
      REAL*8, INTENT(OUT) :: Bx, By, Bz
      DIMENSION pvecin(3),pvecou(3)

      DATA B0x/0.0/
      DATA B0y/0.0/
      !DATA B0z/21160.0/
      !new dipole strength 2012
C     DATA B0z/26191.0/  !for references see KRONIAN_HIGHER
      DATA B0z/21191.0/  !for references see KRONIAN_HIGHER

C     We will  first calculate the field of the dipole
      CALL DIPOLE(B0x,B0y,B0z,x,y,z,Bxd,Byd,Bzd)

C     Now calculate the perpendicular dipole shielding field
      IF( (z .eq. 0) .and. (y .eq. 0)) THEN
        phi=0.0
        GO TO 1
      END IF
      phi = datan2(z,y)
  1   rho=y*dcos(phi)+z*dsin(phi)

      Call B_mp_perp(rho,phi,x,Brho1,Bphi1,Bx1,8,r0or)

      By1 = Brho1*dcos(phi) - Bphi1*dsin(phi)
      Bz1 = Brho1*dsin(phi) + Bphi1*dcos(phi)

      Bx = Bxd - Bx1
      By = Byd - By1
      Bz = Bzd - Bz1

      RETURN
      END SUBROUTINE dipole_shielded

C***********************************************************************

      SUBROUTINE DIPOLE(B0x,B0y,B0z,x,y,z,Bx,By,Bz)
      IMPLICIT REAL*8(A-H,O-Z)
      REAL*8, INTENT(IN) :: B0x,B0y,B0z
      REAL*8, INTENT(IN) :: x,y,z
      REAL*8, INTENT(OUT) :: Bx,By,Bz
C     Calculates the field of Jupiter's dipole for shielding in the magnetopause
C     INPUT: (B0x, B0y, B0z) is the dipole moment
C             x, y, z is the position vector
C     OUTPUT: Bx, By, Bz is the output field vector

      r=dsqrt(x**2+y**2+z**2)
      a11=(3*x**2-r**2)/r**5
      a12=(3*x*y)/r**5
      a13=(3*x*z)/r**5
      a21=(3*x*y)/r**5
      a22=(3*y**2-r**2)/r**5
      a23=(3*y*z)/r**5
      a31=(3*x*z)/r**5
      a32=(3*y*z)/r**5
      a33=(3*z**2-r**2)/r**5
      Bx=a11*B0x+a12*B0y+a13*B0z
      By=a21*B0x+a22*B0y+a23*B0z
      Bz=a31*B0x+a32*B0y+a33*B0z
      RETURN
      END SUBROUTINE DIPOLE

C***********************************************************************

      SUBROUTINE B_mp_perp(rho,phi,x,Bperpr,Bperpf,Bperpx,M,r0or)
      IMPLICIT REAL*8(A-H,O-Z)
      REAL*8, INTENT(IN) :: rho,phi,x
      INTEGER, INTENT(IN) :: M
      REAL*8, INTENT(IN) :: r0or
      REAL*8, INTENT(OUT) :: Bperpr,Bperpf,Bperpx

      Dimension a(8),b(8)

      Data a/3.32125816385087D-004, -0.784083155434928,     
     ~  -0.416208789540178, -3.60127324820496D-002, -0.273188407067209,     
     ~   0.736497845500708, -0.939713175408542,0.433616725727916/     
      Data b/9.24095494266823, 30.4656008995161,     
     ~     18.1008047927882, 60.6014301282221,  125.347578447821,     
     ~     253.707136094045, 509.877831023220, 1021.96079781869/ 

      Bperpr=0.
      Bperpf=0.
      Bperpx=0.

      Do 4 K=1,M
          Term1=Dsin(phi)*Dexp(x/(b(K)/r0or))
          Term2=bessj1(rho/(b(K)/r0or))/(rho/(b(K)/r0or))
     ~         -bessj0(rho/(b(K)/r0or))
          Term3=-Dcos(phi)*Dexp(x/(b(K)/r0or))
          Term4=bessj1(rho/(b(K)/r0or))/(rho/(b(K)/r0or))
          Term5=-Dsin(phi)*Dexp(x/(b(K)/r0or))
          Term6=bessj1(rho/(b(K)/r0or))
          Bperpr=Bperpr + a(K)*r0or**3*Term1*Term2
          Bperpf=Bperpf + a(K)*r0or**3*Term3*Term4
          Bperpx=Bperpx + a(K)*r0or**3*Term5*Term6
  4   CONTINUE

      RETURN
      END SUBROUTINE B_mp_perp

C***********************************************************************

      REAL*8 FUNCTION DBSJ2(x)
      IMPLICIT REAL*8(A-H,O-Z)
      REAL*8, INTENT(IN) :: x

      DBSJ2 = (2.0/x) * bessj1(x) - bessj0(x)
      RETURN
      END FUNCTION DBSJ2

C***********************************************************************

      REAL*8 FUNCTION DBSJ3(x)
      IMPLICIT REAL*8(A-H,O-Z)
      REAL*8, INTENT(IN) :: x

      DBSJ3 = (4.0/x) * DBSJ2(x) - bessj1(x)
      RETURN
      END FUNCTION DBSJ3

C***********************************************************************
C     Hannes needs his own Bessl functions
C***********************************************************************

C     Taken from: http://audiolab.uwaterloo.ca/~jeffb/thesis/node50.html
      REAL*8 FUNCTION bessj1(x)
      IMPLICIT REAL*8(A-H,O-Z)
      REAL*8, INTENT(IN) :: x

      data r1,r2,r3,r4,r5,r6/72362614232.d0,-7895059235.d0,
     #242396853.1d0,
     #-2972611.439d0,15704.48260d0,-30.160366606d0/,
     #s1,s2,s3,s4,s5,s6/144725228442.d0,2300535178.d0,
     #18583304.74d0,99447.43394d0,376.9991397d0,1.d0/
      data p1,p2,p3,p4,p5/1.d0,.183105d-2,-.3516396496d-4,
     #.2457520174d-5,
     #-.240337019d-6/, q1,q2,q3,q4,q5
     #/.04687499995d0,-.2002690873d-3,
     #.8449199096d-5,-.88228987d-6,.105787412d-6/

      IF (dabs(x).lt.8) THEN 
          y=x**2
          bessj1=x*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))))
     #            /(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6)))))
      ELSE
          ax=dabs(x)
          z=8./ax
          y=z**2
          xx=ax-2.356194491
          bessj1=dsqrt(.636619772/ax)*(dcos(xx)*(p1+y*(p2+y*(p3+y*
     #        (p4+y*p5))))-z*dsin(xx)*(q1+y*(q2+y*(q3+y*(q4+y*q5)))))
     #        *dsign(dble(1.),x)
      END IF

      RETURN
      END FUNCTION bessj1

C***********************************************************************

C     Taken from: http://audiolab.uwaterloo.ca/~jeffb/thesis/node50.html
      REAL*8 FUNCTION bessj0(x)
      IMPLICIT REAL*8(A-H,O-Z)
      REAL*8, INTENT(IN) :: x

      data p1,p2,p3,p4,p5/1.d0,-.1098628627d-2,.2734510407d-4,
     #-.2073370639d-5,.2093887211d-6/
     #,q1,q2,q3,q4,q5/-.1562499995d-1,
     #.1430488765d-3,-.6911147651d-5,.7621095161d-6,-.934945152d-7/
      data r1,r2,r3,r4,r5,r6/57568490574.d0,-13362590354.d0,
     #651619640.7d0,
     #-11214424.18d0,77392.33017d0,-184.9052456d0/,
     #s1,s2,s3,s4,s5,s6/57568490411.d0,1029532985.d0,
     #9494680.718d0,59272.64853d0,267.8532712d0,1.d0/

      IF(abs(x).lt.8) THEN 
        y=x**2
        bessj0=(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))))
     #        /(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6)))))
      else
        ax=dabs(x)
        z=8./ax
        y=z**2
        xx=ax-.785398164
        bessj0=dsqrt(.636619772/ax)*(dcos(xx)*(p1+y*(p2+y*(p3+
     #        y*(p4+y*p5))))-z*dsin(xx)*(q1+y*(q2+y*(q3+y*(q4+y*q5)))))
      END IF

      RETURN
      END FUNCTION bessj0

C***********************************************************************
      SUBROUTINE KROT(From,To,Vecin,Vecout,time,epoch)
C***********************************************************************
C INPUTS: From, a chracter*3 variable denoting the incoming coordinate system
C     To:   a chracter*3 variable denoting the outgoing coordinate system
C     Vecin a variable of dimension 3 containing the incoming vector
C     Vecout: a variable of dimension 3 containing the outgoing vector
C     time  A double precision variable denoting number of seconds from an epoch 
C     epoch   A five letter character which can be either 'ctime', 'etime' or 'J2000'/upper or lower case
C
C
C Rotates the vector vecin from the coordinate system "From" to the
C system "To". The rotated vector is placed in the vector Vecout.
C
C We always first rotate the Vecin in "from" system to system III. Then we rotate 
C it to the "To"coordinate system
C
C The supported coordinate systems are:
C     S3C System III Cartesian (right-handed)
C     KSO Kronian  -Sun-Orbital
C     KSM Kronian-Sun-Magnetic (So far same as KSO)
C         DIP Dipole (cartesian)
C***********************************************************************
      IMPLICIT Real*8(A-H,O-Z)
      REAL*8, INTENT(IN):: time
      CHARACTER*3, INTENT(IN) :: From, To
      REAL*8, INTENT(IN) :: Vecin(3)
      CHARACTER*5, INTENT(IN) :: epoch
      REAL*8, INTENT(OUT) :: Vecout(3)
      REAL*8 :: vector(3)
      REAL*8 :: dipole(3,3), first(3,3), second(3,3), dummy(3,3)
C     REAL*4 stheta,sphi
      REAL*8 :: jtime,eetime
      PARAMETER (PI=3.1415927,twopi=2.*PI,radian=PI/180.,degree=180./PI)
C     PARAMETER (TH_DIP=0.0,PH_DIP=0.0) !The dipole parameters
      DATA dipole(3,1),dipole(3,2)/0.,0./
      DATA dipole(3,3)/1.0/                              
      DATA dipole(2,1),dipole(2,2),dipole(2,3)/0.,1.0,0/
      DATA dipole(1,1),dipole(1,2)/1.0,0./
      DATA dipole(1,3)/0.0/

C     write(6,*)"dipole: ",dipole

C     First ascertain the time epoch

      IF ((epoch(1:5).eq."j2000").or.(epoch(1:5).eq."J2000")) then
          jtime=time
          GO TO 101
      END IF

      IF ((epoch(1:5).eq."i2000").or.(epoch(1:5).eq."I2000")) then
          jtime=time+32.184D0
          GO TO 101
      END IF

      IF (((epoch(1:1) .eq. "c") .or. epoch(1:1) .eq. "C") .and.
     ~     ((epoch(2:2) .eq. "t") .or. epoch(2:2) .eq. "T") .and.
     ~     ((epoch(3:3) .eq. "i") .or. epoch(3:3) .eq. "I") .and.
     ~     ((epoch(4:4) .eq. "m") .or. epoch(4:4) .eq. "M") .and.
     ~     ((epoch(5:5) .eq. "e") .or. epoch(5:5) .eq. "E")) then
          jtime=eetime(time)
          GO TO 101
      END IF
C
      IF (((epoch(1:1) .eq. "e") .or. epoch(1:1) .eq. "E") .and.
     ~     ((epoch(2:2) .eq. "t") .or. epoch(2:2) .eq. "T") .and.
     ~     ((epoch(3:3) .eq. "i") .or. epoch(3:3) .eq. "I") .and.
     ~     ((epoch(4:4) .eq. "m") .or. epoch(4:4) .eq. "M") .and.
     ~     ((epoch(5:5) .eq. "e") .or. epoch(5:5) .eq. "E")) then
          jtime=time
          GO TO 101
      END IF

      Write(*,*)" I do not understand the time epoch."
      Go to 100
  
C     Next ascertain the incoming coordinate system
C     Is it system 3?
  101 IF(((From(1:1) .eq. "s") .or. from(1:1) .eq. "S") .and.
     ~   ((From(2:2) .eq. "3") .or. from(2:2) .eq. "3") .and.
     ~   ((From(3:3) .eq. "c") .or. from(3:3) .eq. "C")) Go to 1
C     Is it JSO?
      IF(((From(1:1) .eq. "k") .or. from(1:1) .eq. "K") .and.
     ~   ((From(2:2) .eq. "s") .or. from(2:2) .eq. "S") .and.
     ~   ((From(3:3) .eq. "o") .or. from(3:3) .eq. "O")) Go to 2
C     Is it JSM?
      IF(((From(1:1) .eq. "k") .or. from(1:1) .eq. "K") .and.
     ~   ((From(2:2) .eq. "s") .or. from(2:2) .eq. "S") .and.
     ~   ((From(3:3) .eq. "m") .or. from(3:3) .eq. "M")) Go to 3
C     Is it Dipole?
      IF(((From(1:1) .eq. "d") .or. from(1:1) .eq. "D") .and.
     ~   ((From(2:2) .eq. "i") .or. from(2:2) .eq. "I") .and.
     ~   ((From(3:3) .eq. "p") .or. from(3:3) .eq. "P")) Go to 4
C     Is it Dipole Sun?
      IF(((From(1:1) .eq. "d") .or. from(1:1) .eq. "D") .and.
     ~   ((From(2:2) .eq. "i") .or. from(2:2) .eq. "I") .and.
     ~   ((From(3:3) .eq. "s") .or. from(3:3) .eq. "S")) Go to 5

      Write(*,*)" I do not understand the incoming coordinate system"
      GO TO 100
   1  CONTINUE
C     Write(*,*) " Incoming coordinate system is S3 Cartesian"
      First(1,1)=1.0
      First(1,2)=0.0
      First(1,3)=0.0
      First(2,1)=0.0
      First(2,2)=1.0
      First(2,3)=0.0
      First(3,1)=0.0
      First(3,2)=0.0
      First(3,3)=1.0
      GO TO 60
   2  CONTINUE

C     Write(*,*) " Incoming coordinate system is KSO"

C     Calculate the rotation matrix to go from KSO to S3R
      call KSun(jtime,stheta,sphi,Ztheta,Zphi)

C     Define X component in system III of XKSO unit vector etc.
      dummy(1,1)=cos(stheta)*cos(sphi)
      dummy(1,2)=cos(stheta)*sin(sphi)
      dummy(1,3)=sin(stheta)
C     Calculate the Z axis of the KSO coordinate system.

C     Define X component in system III of ZKSO unit vector etc.
      dummy(3,1)=cos(Ztheta)*cos(Zphi)
      dummy(3,2)=cos(Ztheta)*sin(Zphi)
      dummy(3,3)=sin(Ztheta)

C     Define X component in system III of YKSO unit vector etc.
      dummy(2,1)=dummy(3,2)*dummy(1,3)-dummy(3,3)*dummy(1,2)
      dummy(2,2)=dummy(3,3)*dummy(1,1)-dummy(3,1)*dummy(1,3)
      dummy(2,3)=dummy(3,1)*dummy(1,2)-dummy(3,2)*dummy(1,1)

C     The transpose of the dummy matrix takes KSO to System 3.
      DO 21 i=1,3
      DO 21 j=1,3
  21    first(i,j)=dummy(j,i)
      GO TO 60

  3   CONTINUE
C     Write(*,*) " Incoming coordinate system is KSM"

      CALL KSun(jtime,stheta,sphi,Ztheta,Zphi)

C     Now define the KSM transpose vector
C     Define X component in system III of XJSM unit vector etc.

      dummy(1,1)=cos(stheta)*cos(sphi)
      dummy(1,2)=cos(stheta)*sin(sphi)
      dummy(1,3)=sin(stheta)
C     Now define the Y vector so that it is perpendicular to the dipole vector and X 
      dummy(2,1)=dummy(1,3)*dipole(3,2)-dummy(1,2)*dipole(3,3)
      dummy(2,2)=dummy(1,1)*dipole(3,3)-dummy(1,3)*dipole(3,1)
      dummy(2,3)=dummy(1,2)*dipole(3,1)-dummy(1,1)*dipole(3,2)
      denom=sqrt(dummy(2,1)**2+dummy(2,2)**2+dummy(2,3)**2)
      dummy(2,1)=dummy(2,1)/denom
      dummy(2,2)=dummy(2,2)/denom
      dummy(2,3)=dummy(2,3)/denom
C     Now define the z vector
      dummy(3,1)=dummy(1,2)*dummy(2,3)-dummy(1,3)*dummy(2,2)
      dummy(3,2)=dummy(1,3)*dummy(2,1)-dummy(1,1)*dummy(2,3)
      dummy(3,3)=dummy(1,1)*dummy(2,2)-dummy(1,2)*dummy(2,1)
      DO 30 i=1,3
      DO 30 j=1,3
   30   first(i,j)=dummy(j,i)
      GO TO 60
C
    4 CONTINUE
C     Write(*,*) " Incoming coordinate system is DIP"
      DO 40 i=1,3
      DO 40 j=1,3
   40   first(i,j)=dipole(j,i)
      GO TO 60
C
    5 CONTINUE
C     Write(*,*) " Incoming coordinate system is DIS"

C     Define  ZDIS unit vector 
      dummy(3,1)= dipole(3,1)
      dummy(3,2)= dipole(3,2)
      dummy(3,3)= dipole(3,3)

      CALL KSun(jtime,stheta,sphi,Ztheta,Zphi)


C     Define  XKSO unit vector 
      dummy(1,1)=cos(stheta)*cos(sphi)
      dummy(1,2)=cos(stheta)*sin(sphi)
      dummy(1,3)=sin(stheta)

C     Define YDIS unit vector (as normal to ZDIS and XKSO)
      dummy(2,1)=dummy(1,3)*dummy(3,2)-dummy(1,2)*dummy(3,3)
      dummy(2,2)=dummy(1,1)*dummy(3,3)-dummy(1,3)*dummy(3,1)
      dummy(2,3)=dummy(1,2)*dummy(3,1)-dummy(1,1)*dummy(3,2)
      denom=sqrt(dummy(2,1)**2+dummy(2,2)**2+dummy(2,3)**2)
      dummy(2,1)=dummy(2,1)/denom
      dummy(2,2)=dummy(2,2)/denom
      dummy(2,3)=dummy(2,3)/denom

C     Define XDIS unit vector as normal to YDIS and ZDIS
C
      dummy(1,1)=dummy(2,2)*dummy(3,3)-dummy(2,3)*dummy(3,2)
      dummy(1,2)=dummy(2,3)*dummy(3,1)-dummy(2,1)*dummy(3,3)
      dummy(1,3)=dummy(2,1)*dummy(3,2)-dummy(2,2)*dummy(3,1)


C     The transpose of the dummy matrix takes DIS to System 3.
      DO 22 i=1,3
      DO 22 j=1,3
   22   first(i,j)=dummy(j,i)
C
   60 CONTINUE
C
C     Now Figure out the outgoing coordinate system

      IF(((To(1:1) .eq. "s") .or. To(1:1) .eq. "S") .and.
     ~   ((To(2:2) .eq. "3") .or. To(2:2) .eq. "3") .and.
     ~   ((To(3:3) .eq. "c") .or. To(3:3) .eq. "C")) GO TO 6
C     Is it JSO?
      IF(((To(1:1) .eq. "k") .or. To(1:1) .eq. "K") .and.
     ~   ((To(2:2) .eq. "s") .or. To(2:2) .eq. "S") .and.
     ~   ((To(3:3) .eq. "o") .or. To(3:3) .eq. "O")) GO TO 7
C     Is it JSM?
      IF(((To(1:1) .eq. "k") .or. To(1:1) .eq. "K") .and.
     ~   ((To(2:2) .eq. "s") .or. To(2:2) .eq. "S") .and.
     ~   ((To(3:3) .eq. "m") .or. To(3:3) .eq. "M")) GO TO 8
C     Is it Dipole?
      IF(((To(1:1) .eq. "d") .or. To(1:1) .eq. "D") .and.
     ~   ((To(2:2) .eq. "i") .or. To(2:2) .eq. "I") .and.
     ~   ((To(3:3) .eq. "p") .or. To(3:3) .eq. "P")) GO TO 9
C     Is it DIS?
      IF(((To(1:1) .eq. "d") .or. To(1:1) .eq. "D") .and.
     ~   ((To(2:2) .eq. "i") .or. To(2:2) .eq. "I") .and.
     ~   ((To(3:3) .eq. "s") .or. To(3:3) .eq. "S")) GO TO 10
C
      Write(*,*)" I do not understand the outgoing coordinate system"
      GO TO 100

    6 CONTINUE
C     Write(*,*) " The outgoing coordinate system is S3 cartesian"
      Second(1,1)=1.0
      Second(1,2)=0.0
      Second(1,3)=0.0
      Second(2,1)=0.0
      Second(2,2)=1.0
      Second(2,3)=0.0
      Second(3,1)=0.0
      Second(3,2)=0.0
      Second(3,3)=1.0
      GO TO 70
C
    7 CONTINUE
C     Write(*,*) " The outgoing coordinate system is KSO"

C     Get the matrix that Rotates from System III cartesian into KSO. 
      CALL KSun(jtime,stheta,sphi,Ztheta,Zphi)

C     Define X component in system III of XKSO unit vector etc.
      Second(1,1)=cos(stheta)*cos(sphi)
      Second(1,2)=cos(stheta)*sin(sphi)
      Second(1,3)=sin(stheta)

C     Define X component in system III of ZKSO unit vector etc.
      Second(3,1)=cos(Ztheta)*cos(Zphi)
      Second(3,2)=cos(Ztheta)*sin(Zphi)
      Second(3,3)=sin(Ztheta)

C     Check orthogonality of the new coordinate system
      dotprod=second(3,1)*second(1,1)+second(3,2)*second(1,2)+
     .second(3,3)*second(1,3)
C     Write(*,*) dotprod
C     Define X component in system III of YKSO unit vector etc.

      Second(2,1)=Second(3,2)*Second(1,3)-Second(3,3)*Second(1,2)
      Second(2,2)=Second(3,3)*Second(1,1)-Second(3,1)*Second(1,3)
      Second(2,3)=Second(3,1)*Second(1,2)-Second(3,2)*Second(1,1)
C     seconm= sqrt(second(2,1)**2+second(2,2)**2+second(2,3)**2)
C     Write(*,*)seconm
      GO TO 70
   8  CONTINUE
C     Write(*,*) " The outgoing coordinate system is KSM"

C     Get the matrix that Rotates from System III cartesian into KSM. 
      CALL KSun(jtime,stheta,sphi,Ztheta,Zphi)
C     Define X component in system III of XKSM unit vector etc.
      Second(1,1)=cos(stheta)*cos(sphi)
      Second(1,2)=cos(stheta)*sin(sphi)
      Second(1,3)=sin(stheta)
C     Now define the Y vector so that it is perpendicular to  the dipole vector and X 
      Second(2,1)=second(1,3)*dipole(3,2)-second(1,2)*dipole(3,3)
      Second(2,2)=second(1,1)*dipole(3,3)-second(1,3)*dipole(3,1)
      Second(2,3)=second(1,2)*dipole(3,1)-second(1,1)*dipole(3,2)
      denom=sqrt(second(2,1)**2+second(2,2)**2+second(2,3)**2)
      Second(2,1)=second(2,1)/denom
      Second(2,2)=second(2,2)/denom
      Second(2,3)=second(2,3)/denom
C     Now define the z vector
      Second(3,1)=second(1,2)*second(2,3)-second(1,3)*second(2,2)
      Second(3,2)=second(1,3)*second(2,1)-second(1,1)*second(2,3)
      Second(3,3)=second(1,1)*second(2,2)-second(1,2)*second(2,1)
      GO TO 70
C
    9 CONTINUE
C     Write(*,*) " The outgoing coordinate system is Dipole"
      DO 80 i=1,3
      DO 80 j=1,3
   80   Second(i,j)=dipole(i,j)  
      GO TO 70
   10 CONTINUE
C     Write(*,*) " Outgoing coordinate system is DIS"

C     Define  ZDIS unit vector 
      dummy(3,1)= dipole(3,1)
      dummy(3,2)= dipole(3,2)
      dummy(3,3)= dipole(3,3)

      ! write(6,*)"DUMMY1: ",dummy

      CALL KSun(jtime,stheta,sphi,Ztheta,Zphi)
C     Define  XKSO unit vector 
      dummy(1,1)=cos(stheta)*cos(sphi)
      dummy(1,2)=cos(stheta)*sin(sphi)
      dummy(1,3)=sin(stheta)
C
      !write(6,*)"stheta: ",stheta
      !write(6,*)"sphi: ",sphi
      !write(6,*)"DUMMY2: ",dummy

C     Define YDIS unit vector (as normal to ZDIS and XKSO)
      dummy(2,1)=dummy(1,3)*dummy(3,2)-dummy(1,2)*dummy(3,3)
      dummy(2,2)=dummy(1,1)*dummy(3,3)-dummy(1,3)*dummy(3,1)
      dummy(2,3)=dummy(1,2)*dummy(3,1)-dummy(1,1)*dummy(3,2)
      denom=sqrt(dummy(2,1)**2+dummy(2,2)**2+dummy(2,3)**2)
      dummy(2,1)=dummy(2,1)/denom
      dummy(2,2)=dummy(2,2)/denom
      dummy(2,3)=dummy(2,3)/denom

      !write(6,*)"DUMMY3: ",dummy

C     Define XDIS unit vector as normal to YDIS and ZDIS
      dummy(1,1)=dummy(2,2)*dummy(3,3)-dummy(2,3)*dummy(3,2)
      dummy(1,2)=dummy(2,3)*dummy(3,1)-dummy(2,1)*dummy(3,3)
      dummy(1,3)=dummy(2,1)*dummy(3,2)-dummy(2,2)*dummy(3,1)

C     The  dummy matrix takes  System 3 to DIS.
      DO 23 i=1,3
      DO 23 j=1,3
   23 second(i,j)=dummy(i,j)
C
   70 CONTINUE
C     Now multiply vecin with first and second matrices to get the vecout
C     Write(*,*)" calling matmult 1"
      CALL matmult(first,vecin,vector,3,3,1)
C     Write(*,*)" calling matmult 2"
      CALL matmult(second,vector,vecout,3,3,1)

      !write(6,*)"dipole: ",dipole
      !write(6,*)"DUMMY: ",dummy
      !write(6,*)"FIRST: ",first
      !write(6,*)"sECOND: ",second

  100 CONTINUE
      RETURN
      END SUBROUTINE KROT
C***********************************************************************

      SUBROUTINE KSUN(time,stheta,sphi,Ztheta,Zphi)
C     INPUT:  J2000 time of the data point 
C     OUTPUTS: stheta, sphi, latitude and longitude  (in radians) of the Sun in system III (RH).  
C     OUTPUTS: Ztheta, Zphi, latitude and longitude  (in radians) of the Zkso in system III (RH).  
C     The equations are written in J2000 convention followed by PDS. 
C
C     Last updated August 26, 2004.
C
C     theta=a1*cos(omegay*t)+a2*sin(omegay*t)+
C     a3*cos(2.*omegay*t)+a4*sin(2.*omegay*t)+
C     a5*cos(3.*omegay*t)+a6*sin(3.*omegay*t)+
C     a7*cos(4.*omegay*t)+a8*sin(4.*omegay*t)+
C         a9*t**2+a10*t+a11
C     fphi=b1*cos(omegay*t)+b2*sin(omegay*t)+
C         b3*cos(2.*omegay*t)+b4*sin(2.*omegay*t)+
C         b5*cos(3.*omegay*t)+b6*sin(3.*omegay*t)+
C         b7*cos(4.*omegay*t)+b8*sin(4.*omegay*t)+
C     b9*t**2+b10*t+b11 (I assume I have despun Saturn first.)
C
C     Then we rotate into the System III coordinates
C
C     Time is a double precision variable and is assumed to be in J2000.
C     Theta and phi are single precision variables.
C   
C     omega is Saturn's rotation rate
C     omegaz is saturn's rotation rate for use of Zkso lat and lon
C     omegay is Saturn's yearly orbital rate
      Implicit Real*8(A-H,O-Z)
      REAL*8, INTENT(IN) :: time
      REAL*8, INTENT(OUT) :: stheta,sphi,Ztheta,Zphi
      REAL*8 :: aa(11),bb(11), cc(11),x(11)
      REAL*8 :: fphi,t
C     REAL*4 stheta,sphi,ztheta,Zphi
      REAL*8, PARAMETER :: PI=3.141592654d0
      REAL*8, PARAMETER :: twopi=2.*PI, radian=PI/180., degree=180./PI
      REAL*8, PARAMETER :: yrsat=10759.22d0*86400.0d0
      REAL*8, PARAMETER :: omega=360.D0/(10.65622222D0*3600.D0)
C     REAL*8, PARAMETER :: omega=360.D0/(10.53333333D0*3600.D0) ![Helled,2015]
      REAL*8, PARAMETER :: omegay=2.*PI/yrsat
      REAL*8, PARAMETER :: omegaz=360.D0/(10.65622222D0*3600.D0)+3.880688963D-7
C     REAL*8, PARAMETER :: omegaz=360.D0/(10.53333333D0*3600.D0)+3.880688963D-7
      REAL*8, PARAMETER :: year=86400D0*.36525D3, D360=360.D0
      REAL*8, PARAMETER :: jtime1=-.825767955816D+09
      REAL*8, PARAMETER :: Zthetd=63.271d0
C     REAL*8, PARAMETER :: Zthetd=-26.729
      DATA aa  /-.26237934D+02,.30525104D+01,-.13686733D+01,
     +.10182726D+00,-.30897805D+00,.89277949D-01,-.39638771D-01,
     +.97499653D-02,.40974159D-04,.14075684D-02,.13505775D+01/
C
      DATA bb /-.32108124D+00,.58569866D+01,.71266272D+00,.33244539D+01,
     +.47529267D-01,.32770362D+00,.46935622D-01,.10720536D+00,
     + -.50594764D-03,.29439612D-01,.26423581D+03/
C
      DATA cc/-.13200649D-02,-.18358429D-02,.64658927D-03,.12800487D-02,
     +.17618936D-04,-.58790898D-03,.49804081D-04,.42372137D-03,
     +.14744891D-04,.25369744D-01,.77943328D+02/

C     First calculate the latitude and longitude in non-rotating Saturn coordinates.

C     Calculate the best fit theta and fphi
      t=time-jtime1
      x(1)=cos(omegay*t)
      x(2)=sin(omegay*t)
      x(3)=cos(2.*omegay*t)
      x(4)=sin(2.*omegay*t)
      x(5)=cos(3.*omegay*t)
      x(6)=sin(3.*omegay*t)
      x(7)=cos(4.*omegay*t)
      x(8)=sin(4.*omegay*t)
      x(9)=(t/year)**2
      x(10)=t/year
      x(11)=1.0
      stheta=0.
C     fphi is phi in Saturn fixed (non-rotating) coordinate
      fphi=0.0
      Zfphi=0.0
      DO 1 j=1,11
          fphi=fphi+bb(j)*x(j)
          Zfphi=Zfphi+cc(j)*x(j)
    1 stheta=stheta+aa(j)*x(j)

C     Now rotate the longitude to Saturn System III
C     First Add the rotation of Saturn around the Sun.
C     fphi is the phi of the Sun as unspinning Saturn goes around the sun
      fphi=DMod(fphi+t/yrsat*360.d0, D360)
C     Next add the rotation of Saturn around its axis.
      sphi=DMod(fphi-t*omega, D360)
      IF (sphi .lt. 0.0) sphi = sphi+360.0
      sphi=sphi*radian
      stheta=stheta*radian

C     Next rotate the longitude of Zfphi to Saturn System III
C     First Add the rotation of Saturn around the Sun.
C     Zfphi is the phi of Zkso as unspinning Saturn goes around the sun
      Zfphi=DMod(Zfphi+t/yrsat*360.d0+180.D0, D360)
C     Next add the rotation of Saturn around its axis.
      Zphi=DMod(Zfphi-t*omegaz, D360)
      IF (Zphi .lt. 0.0) Zphi = Zphi+360.0
      Zphi=Zphi*radian
      Ztheta=Zthetd*radian

      RETURN
      END SUBROUTINE KSUN

C***********************************************************************

      REAL*8 FUNCTION eetime(ctime)
      IMPLICIT real*8(A-H,O-z)
      REAL*8, INTENT(IN) :: ctime
      Tcor=0.
      if(ctime .ge. 189302400.000) Tcor = Tcor+10 
      if(ctime .ge. 205027200.000) Tcor = Tcor+1  
      if(ctime .ge. 220924800.000) Tcor = Tcor+1
      if(ctime .ge. 252460800.000) Tcor = Tcor+1
      if(ctime .ge. 283996800.000) Tcor = Tcor+1
      if(ctime .ge. 315532800.000) Tcor = Tcor+1
      if(ctime .ge. 347155200.000) Tcor = Tcor+1
      if(ctime .ge. 378691200.000) Tcor = Tcor+1
      if(ctime .ge. 410227200.000) Tcor = Tcor+1
      if(ctime .ge. 441763200.000) Tcor = Tcor+1
      if(ctime .ge. 489024000.000) Tcor = Tcor+1
      if(ctime .ge. 520560000.000) Tcor = Tcor+1
      if(ctime .ge. 552096000.000) Tcor = Tcor+1
      if(ctime .ge. 615254400.000) Tcor = Tcor+1
      if(ctime .ge. 694224000.000) Tcor = Tcor+1
      if(ctime .ge. 757382400.000) Tcor = Tcor+1
      if(ctime .ge. 788918400.000) Tcor = Tcor+1
      if(ctime .ge. 836179200.000) Tcor = Tcor+1
      if(ctime .ge. 867715200.000) Tcor = Tcor+1
      if(ctime .ge. 899251200.000) Tcor = Tcor+1
      if(ctime .ge. 946684800.000) Tcor = Tcor+1
      if(ctime .ge. 993945600.000) Tcor = Tcor+1
      if(ctime .ge. 1041379200.000) Tcor = Tcor+1
      eetime=ctime+dble(Tcor)-.1072958367816D10
      RETURN
      END FUNCTION eetime

C***********************************************************************

      SUBROUTINE matmult(xx,yy,zz,nxrow,nxcol,nycol)
      IMPLICIT REAL*8(A-H,O-Z)
      INTEGER, INTENT(IN)  :: nxrow,nxcol,nycol
      REAL*8,  INTENT(IN)  :: xx(nxrow,nxcol)
      REAL*8,  INTENT(IN)  :: yy(nxcol,nycol)
      REAL*8,  INTENT(OUT) :: zz(nxrow,nycol)
      ! DIMENSION(:), INTENT(IN) :: xx, yy
      ! DIMENSION xx(nxrow,nxcol),yy(nxcol,nycol),zz(nxrow,nycol)
      Do 1 i=1, nxrow
        Do 1 j=1, nycol
          zz(i,j) = 0.
          Do 1 k=1, nxcol
            zz(i,j) = zz(i,j) + xx(i,k) * yy(k,j)
    1 CONTINUE

      RETURN
      END SUBROUTINE matmult

C******************************************************************
      REAL*8 FUNCTION ctimer(iyr,imon,iday,ihr,imin,sec)
      INTEGER, INTENT(IN) :: iyr,imon,iday,ihr,imin
      INTEGER*4           :: ndays,doy
      REAL*8              :: sec

C     First calculate the number of days from Jan 1, 1966
      ndays = 0
      IF  (iyr .ge. 1966) THEN
        DO 1 i=1966, iyr-1
        ndays=ndays+365
        IF(((mod(i,4) .eq. 0) .and. (mod(i,100) .ne. 0)) .or. 
     +      (mod(i,400)) .eq. 0) ndays=ndays+1
   1    CONTINUE

C       Now add the number of days of the current year
        ndays = ndays+doy(iyr,imon,iday)-1  
        ctimer = dble(ndays)*86400D0+dble(ihr)*3600D0+dble(imin)*60D0+
     +           dble(sec)
        GO TO 4
      END IF

C     Calculate the seconds for the negative years
      IF  (iyr .lt. 1966) THEN
        DO 2 i=iyr, 1965
        ndays=ndays-365
        IF (((mod(i,4) .eq. 0) .and. (mod(i,100) .ne. 0)) .or. 
     +       (mod(i,400)) .eq. 0) ndays=ndays-1
    2   CONTINUE
C       Now subtract the number of days of the current year 
        ndays=ndays+doy(iyr,imon,iday)
      END IF

      ctimer = dble(ndays-1)*86400D0+dble(ihr)*3600D0+dble(imin)*60D0+
     +dble(sec)

    4 CONTINUE
      idyr=doy(iyr,imon,iday)
      Write(*,*) " doy is ",idyr
      Write(*,*)" ndays is ",ndays

      RETURN 
      END FUNCTION ctimer

C***********************************************************************

      INTEGER*4 FUNCTION doy(iyr,imon,iday)
      INTEGER, INTENT(IN)  :: iyr,imon,iday
      INTEGER mon(12)
      data mon /31,28,31,30,31,30,31,31,30,31,30,31/
      doy = 0
      Do 1 i = 2,imon
        doy = doy + mon(i-1)
C       Add an extra day for February
        IF(i .eq. 3) THEN
          IF(((mod(iyr,4) .eq. 0) .and. (mod(iyr,100) .ne. 0)) .or. 
     +        (mod(iyr,400)) .eq. 0) doy=doy+1
        END IF
    1 CONTINUE
      doy = doy + iday
      RETURN 
      END FUNCTION doy

C***********************************************************************

      SUBROUTINE mapit(time,XDIS,YDIS,ZDIS,XMAP,YMAP,ZMAP,r0or,Epoch)
      IMPLICIT REAL*8(A-H,O-Z)
      REAL*8, INTENT(IN) :: time
      REAL*8, INTENT(IN) :: XDIS,YDIS,ZDIS
      REAL*8, INTENT(OUT) :: XMAP,YMAP,ZMAP
      REAL*8, INTENT(IN) :: r0or
      CHARACTER*5, INTENT(IN) :: Epoch
      data RH0/23.0/
      Dimension pvecin(3),pvecou(3)

C     r=sqrt(XKSM**2+YKSM**2+ZKSM**2)

C     r0or = 1.0
      !write(*,*)r0or

      RH = RH0 / r0or !hinging as function of pressure

C     First calculate the dipole tilt (si) in the KSM coordinates
      pvecin(1) = 0.
      pvecin(2) = 0.
      pvecin(3) = 1.
      CALL KROT('S3C','KSM',pvecin,pvecou,time,epoch)

      si = dacos(pvecou(3))
      if(pvecou(1) .lt. 0) si = -si

C     Now calculate the height of the current sheet
      rhoDIS = dsqrt(XDIS**2 + YDIS**2)
      phiLT  = datan2(YDIS,XDIS)
      theta  = datan2(ZDIS,rhoDIS)
      Zcs    = (rhoDIS-RH*dtanh(rhoDIS/RH))*dtan(-si) !sign of si?
      Thcs   = datan2(Zcs,rhoDIS) !sign of Thcs?
  
      if (XDIS .gt. 0.0) Thcs = -Thcs ! change 11/09/2012
  
      rhomap =  rhoDIS*dcos(thcs) + ZDIS*dsin(thcs) 
      Zmap   = -rhoDIS*dsin(thcs) + ZDIS*dcos(thcs)
      Xmap   =  rhomap*dcos(phiLT)
      Ymap   =  rhomap*dsin(phiLT)

      RETURN
      END SUBROUTINE mapit

C***********************************************************************
      SUBROUTINE mapitMP(time,XDIS,YDIS,ZDIS,XMAP,YMAP,ZMAP,r0or,Epoch)
      IMPLICIT REAL*8(A-H,O-Z)
      REAL*8, INTENT(IN) :: time
      REAL*8, INTENT(IN) :: XDIS,YDIS,ZDIS
      REAL*8, INTENT(OUT) :: XMAP,YMAP,ZMAP
      REAL*8, INTENT(IN) :: r0or
      CHARACTER*5, INTENT(IN) :: Epoch
      data RH0/23.0/
      Dimension pvecin(3),pvecou(3)

C     r=sqrt(XKSM**2+YKSM**2+ZKSM**2)

C     r0or=1.0
      !write(*,*)r0or

      RH=RH0/r0or !hinging as function of pressure
      !RH=0.0
  
C     First calculate the dipole tilt (si) in the KSM coordinates
      pvecin(1) = 0.
      pvecin(2) = 0.
      pvecin(3) = 1.
      CALL KROT('S3C','KSM',pvecin,pvecou,time,epoch)
      si = dacos(pvecou(3))
      IF(pvecou(1) .lt. 0) si = -si

C     Now calculate the height of the current sheet
      rhoDIS = dsqrt(XDIS**2+YDIS**2)
      phiLT  = datan2(YDIS,XDIS)
      theta  = datan2(ZDIS,rhoDIS)
      Zcs    = (rhoDIS-RH*dtanh(rhoDIS/RH))*dtan(+si) !sign of si?
      Thcs   = datan2(Zcs,rhoDIS) !sign of Thcs?
  
      IF (XDIS .gt. 0.0) Thcs = -Thcs ! change 11/09/2012
  
      rhomap =  rhoDIS*dcos(thcs)+ZDIS*dsin(thcs) 
      Zmap   = -rhoDIS*dsin(thcs)+ZDIS*dcos(thcs)
  
      Xmap=rhomap*dcos(phiLT)
      Ymap=rhomap*dsin(phiLT)

      RETURN
      END SUBROUTINE mapitMP

C***********************************************************************

      SUBROUTINE Mapped_field(
     ~  time,X,Y,Z,BX,BY,BZ,BXMAP,BYMAP,BZMAP,r0or,Epoch)
      IMPLICIT REAL*8(A-H,O-Z)
      REAL*8, INTENT(IN) :: time
      REAL*8, INTENT(IN) :: X,Y,Z,BX,BY,BZ
      REAL*8, INTENT(OUT) :: BXMAP,BYMAP,BZMAP
      REAL*8, INTENT(IN) :: r0or
      CHARACTER*5, INTENT(IN) :: Epoch

C     DIMENSION pvecin(3),pvecou(3)
C     r0or=1.0
C
C     Now define the nine derivatives
C     These are calculated at the original location
      dx=0.01
      dy=0.01
      dz=0.01
      xp=x+dx
      xm=x-dx
      yp=y+dy
      ym=y-dy
      zp=z+dz
      zm=z-dz
C     We begin with the x derivatives
      CALL mapit(time,xp,y,z,xpp,ypp,zpp,r0or,Epoch)
      CALL mapit(time,xm,y,z,xpm,ypm,zpm,r0or,Epoch)

      dxpdx=(xpp-xpm)/(2.*dx)
      dypdx=(ypp-ypm)/(2.*dx)
      dzpdx=(zpp-zpm)/(2.*dx)

C     Next the y derivatives
C     dxpdy=0.
C     dypdy=1.
C     dzpdy=0.
      CALL mapit(time,x,yp,z,xpp,ypp,zpp,r0or,Epoch)
      CALL mapit(time,x,ym,z,xpm,ypm,zpm,r0or,Epoch)

      dxpdy=(xpp-xpm)/(2.*dy)
      dypdy=(ypp-ypm)/(2.*dy)
      dzpdy=(zpp-zpm)/(2.*dy)

C     Next the z gradients 
      CALL mapit(time,x,y,zp,xpp,ypp,zpp,r0or,Epoch)
      CALL mapit(time,x,y,zm,xpm,ypm,zpm,r0or,Epoch)

      dxpdz=(xpp-xpm)/(2.*dz)
      dypdz=(ypp-ypm)/(2.*dz)
      dzpdz=(zpp-zpm)/(2.*dz)

C     Now calculate the T matrix
C     Calculate the mapped location
      Txx=dypdy*dzpdz-dypdz*dzpdy
      Txy=dxpdz*dzpdy-dxpdy*dzpdz
      Txz=dxpdy*dypdz-dxpdz*dypdy
      Tyx=dypdz*dzpdx-dypdx*dzpdz
      Tyy=dxpdx*dzpdz-dxpdz*dzpdx
      Tyz=dxpdz*dypdx-dxpdx*dypdz
      Tzx=dypdx*dzpdy-dypdy*dzpdx
      Tzy=dxpdy*dzpdx-dxpdx*dzpdy
      Tzz=dxpdx*dypdy-dxpdy*dypdx
    
C     Now calculate the field at the mapped location 
      Bxmap=Txx*Bx+Txy*By+Txz*Bz  
      Bymap=Tyx*Bx+Tyy*By+Tyz*Bz  
      Bzmap=Tzx*Bx+Tzy*By+Tzz*Bz  
      RETURN 
      END SUBROUTINE Mapped_field

C***********************************************************************

      SUBROUTINE shielded_csheetfield(time,R1,T1,P1,RLT,Brm1,Btm1,
     ~Bpm1,NL,C,D,r0or,dql,epoch,Brm1arr,Btm1arr,Bpm1arr)
C     This subroutine calculates the best fit current sheet by modeling only
C     Brho and Bz components. Once the coefficients for the 7 modes are known, 
C     the other subroutines MAPIT and Mapped_field introduce the stretch transformations 
C     and calculate the new values of Brho,Bphi and Bz. 
C     Real*4 stheta,sphi,ztheta,Zphi

      IMPLICIT REAL*8(A-H,O-Z)
      REAL*8, INTENT(IN) :: time
      CHARACTER*5, INTENT(IN) :: Epoch
      REAL*8, INTENT(IN) :: R1, T1, P1
      REAL*8, INTENT(OUT) :: RLT, Brm1, Btm1, Bpm1
      INTEGER, INTENT(IN) :: NL
      REAL*8, INTENT(IN) :: C(NL)
      REAL*8, INTENT(IN) :: D
      REAL*8, INTENT(IN) :: r0or
      REAL*8, INTENT(IN) :: dql
      REAL*8, INTENT(OUT) :: Brm1arr(NL), Btm1arr(NL), Bpm1arr(NL)
      REAL*8 pvecin(3), pvecou(3)
      REAL*8 bvecin(3), bvecou(3)
      DIMENSION f(6,6), beta(6,6)
      DIMENSION Xrho(NL), Xz(NL)

C     RING CURRENT
C     ---------------------------
C     Peak AT 12
C     f(1,1)=    -68706.3686478971
C     f(1,2)=       -2.9523400227
C     f(1,3)=   100047.0971381391
C     f(1,4)=   -30783.7506416514
C     f(1,5)=       -0.0000642075
C     f(1,6)=       -0.0031683487
C     beta(1,1)=     9.2005569531
C     beta(1,2)=     2.1784062411
C     beta(1,3)=     8.8767954666
C     beta(1,4)=     8.2680804765
C     beta(1,5)=     0.1100281701
C     beta(1,6)=     0.4863609871

C     PEAK AT 8
C     f(1,1)=    99323.9271557975
C     f(1,2)=       22.6189468473
C     f(1,3)=      720.7929798854
C     f(1,4)=     -485.4538206945
C     f(1,5)=   -37918.4097869931
C     f(1,6)=   -61520.1671330672
C     beta(1,1)=   5.4851103313
C     beta(1,2)=   1.9438277867
C     beta(1,3)=   2.9025159249
C     beta(1,4)=   2.6449016349
C     beta(1,5)=   5.1921884427
C     beta(1,6)=   5.6729005214

C     PEAK AT 6
C     f(1,1)=     -96480.9970060387
C     f(1,2)=     -733.2997300981
C     f(1,3)=     -17474.3949850542
C     f(1,4)=    13878.3552501633
C     f(1,5)=    99937.8587320000
C     f(1,6)=      988.4579332559
C     beta(1,1)=   4.5216256036
C     beta(1,2)=   2.5132995869
C     beta(1,3)=   3.2135379771
C     beta(1,4)=   3.0509806322
C     beta(1,5)=   4.4629794036
C     beta(1,6)=   6.6822375235

C     PEAK AT 4
C     f(1,1)=     69605.0534869644
C     f(1,2)=     149.0772338679
C     f(1,3)=     78909.3792557924
C     f(1,4)=  -54703.2221035385
C     f(1,5)=  -93857.8599750019
C     f(1,6)=     -59.8779082335
C     beta(1,1)=    3.7051342090
C     beta(1,2)=    2.3672348371
C     beta(1,3)=    3.2236662867
C     beta(1,4)=    3.1391528825
C     beta(1,5)=    3.6275112423
C     beta(1,6)=    7.8003148590

C ------------------------------
C     fi: new mode Aug 2012 
C ------------------------------
C     f(1,1)=   828.7022879010
C     f(1,2)= -2953.8174466019
C     f(1,3)=   587.3538467185
C     f(1,4)=  -237.6901233573
C     f(1,5)=  -303.3940846623
C     f(1,6)=  2500.5187349814
C     beta(1,1)=  59.2734758066
C     beta(1,2)=  34.4268323961
C     beta(1,3)=   4.0855094297
C     beta(1,4)=   3.1278017056
C     beta(1,5)=   8.7602255391
C     beta(1,6)=  30.9667525321
C ------------------------------
C     fi: new mode Aug 2012 start at 1.2 peak at 5
C     f(1,1)=     559.6197083702
C     f(1,2)=   -1968.8800437331
C     f(1,3)=     584.4840512599
C     f(1,4)=    -364.3594622078
C     f(1,5)=    -166.6634660654
C     f(1,6)=    1659.0247186932
C     beta(1,1)=  59.2711630145
C     beta(1,2)=  34.4177680800
C     beta(1,3)=   3.1524952826
C     beta(1,4)=   2.7593883598
C     beta(1,5)=   8.4957542421
C     beta(1,6)=  30.9550941820
C ------------------------------
C     fi: new mode Aug 2012 start at 2 peak at 5
C     f(1,1)=       9.7672928003
C     f(1,2)=      63.5144165850
C     f(1,3)=   -2975.8264925754
C     f(1,4)=    3960.8018682005
C     f(1,5)=     -65.1828634143
C     f(1,6)=    -863.4772958890
C     beta(1,1)=    1.5894275795
C     beta(1,2)=   34.3759030779
C     beta(1,3)=    3.6428139222
C     beta(1,4)=    3.4574753261
C     beta(1,5)=    8.2962286103
C     beta(1,6)=    2.8921770952
C ------------------------------
C     fi: new mode Aug 2012 start at 5 peak at 8
C     fi
C     f(1,1)=   58371.7584883787
C     f(1,2)=      67.2696628374
C     f(1,3)=    4829.5415810648
C     f(1,4)=   -3247.0372511400
C     f(1,5)=  -59240.7786762627
C     f(1,6)=    -748.9828504765
C     beta(1,1)=   4.3882059995
C     beta(1,2)=   2.5209029894
C     beta(1,3)=   3.2052095595
C     beta(1,4)=   3.0487351464
C     beta(1,5)=   4.3379025266
C     beta(1,6)=   8.2846806078

      f(1,1)=     -100.5052977296
      f(1,2)=     5366.2712701130
      f(1,3)=   -15885.9490677644
      f(1,4)=    40732.9959830497
      f(1,5)=   -32129.0816517276
      f(1,6)=     2215.9395715640
      beta(1,1)=     2.7117925416
      beta(1,2)=     4.0941642385
      beta(1,3)=     4.6702786290
      beta(1,4)=     5.7582130292
      beta(1,5)=     6.2062614171
      beta(1,6)=     8.4355483099

      f(2,1)=   -78718.3907963260
      f(2,2)=      -54.8132316673
      f(2,3)=    -2414.9302620727
      f(2,4)=     1616.7325480546
      f(2,5)=    74549.4133726164
      f(2,6)=     5883.6406782211
      beta(2,1)=    5.0017041440
      beta(2,2)=    2.6863215815
      beta(2,3)=    3.5010631982
      beta(2,4)=    3.3102786530
      beta(2,5)=    4.9394705994
      beta(2,6)=    6.1620526042

      f(3,1)=     828.2831671297
      f(3,2)=  -61442.3855659502
      f(3,3)=  -22900.6452859245
      f(3,4)=   40528.1409725691
      f(3,5)=  -42635.0448098891
      f(3,6)=   87296.6633326159
      beta(3,1)=    72.3455716311
      beta(3,2)=    10.1486476486
      beta(3,3)=     8.0002370990
      beta(3,4)=     8.4224589108
      beta(3,5)=    12.8834880849
      beta(3,6)=    11.6138431830

      f(4,1)=        -52.9340216641
      f(4,2)=          4.9006893928
      f(4,3)=        103.5467051654
      f(4,4)=     -16743.1838794522
      f(4,5)=      26132.5659128414
      f(4,6)=      -7214.7617456977
      beta(4,1)=      4.7199163645
      beta(4,2)=      2.5733685428
      beta(4,3)=      9.5232221496
      beta(4,4)=     18.2512781842
      beta(4,5)=     20.4896316537
      beta(4,6)=     59.5697679599

C     f(5,1)=             9.6382346663
C     f(5,2)=            -1.9403694819
C     f(5,3)=          -110.7011194690
C     f(5,4)=          -212.0255553191
C     f(5,5)=        -15859.9202979478
C     f(5,6)=         55821.0012327883
C     beta(5,1)=          4.5856791348
C     beta(5,2)=          2.5276198276
C     beta(5,3)=          8.3627067607
C     beta(5,4)=         17.9703567386
C     beta(5,5)=         36.2969029093
C     beta(5,6)=         46.4652858681

      f(5,1)=     -23.3052813129
      f(5,2)=       2.2671216463
      f(5,3)=      -4.0457878393
      f(5,4)=   -2345.2488224382
      f(5,5)=   89556.6128080047
      f(5,6)=  -70560.3264679953
      beta(5,1)=    4.6355942963
      beta(5,2)=    2.5383712259
      beta(5,3)=    8.4531654549
      beta(5,4)=   20.0059847755
      beta(5,5)=   42.8693646158
      beta(5,6)=   47.3861513994
       
C     L=4
C     taken: use old mode again 
C     f(2,1)=     -60.6568564967
C     f(2,2)=  -73849.9534670401
C     f(2,3)=   97442.8437444251
C     f(2,4)=  -23178.2043425704
C     f(2,5)=      66.3747390965
C     f(2,6)=      -1.4848752731
C     beta(2,1)=   48.1885954914
C     beta(2,2)=    5.7107159217
C     beta(2,3)=    5.6345236716
C     beta(2,4)=    5.3879744819
C     beta(2,5)=   89.0583961619
C     beta(2,6)=   22.6415533854
      
C     L=4
C     new mode Aug 2012 start at 2 peak at 5
C     f(2,1)=      1.0753834757
C     f(2,2)=    -87.8886394188
C     f(2,3)=    213.0301443067
C     f(2,4)=    180.5413030125
C     f(2,5)=    114.7177473372
C     f(2,6)=   -196.5847020514
C     beta(2,1)=     1.5383102306
C     beta(2,2)=    34.3790120138
C     beta(2,3)=     3.4830213559
C     beta(2,4)=     3.2635807863
C     beta(2,5)=     8.3908399544
C     beta(2,6)=     2.6750562404

C     L=8            
C     taken: use old mode again
C     f(3,1)=      310.2640136545
C     f(3,2)=    24690.5295574541
C     f(3,3)=      786.0977568405
C     f(3,4)=   -24122.7960738661
C     f(3,5)=     -259.2663294173
C     f(3,6)=     -598.5252855461
C     beta(3,1)=      48.9404306370
C     beta(3,2)=       5.2832855659
C     beta(3,3)=      16.8568900868
C     beta(3,4)=       5.2434883472
C     beta(3,5)=      89.2258188514
C     beta(3,6)=      25.7902108924
      
   
C     new mode Aug 2012 start at 2 peak at 5
C     f(3,1)=    34039.2023345358
C     f(3,2)=       60.7894247475
C     f(3,3)=    36273.8783023426
C     f(3,4)=   -25286.2224420778
C     f(3,5)=   -44812.4245264180
C     f(3,6)=      460.6339790425
C     beta(3,1)=       3.8750997980
C     beta(3,2)=       2.4309215256
C     beta(3,3)=       3.3669595839
C     beta(3,4)=       3.2782539511
C     beta(3,5)=       3.7923429792
C     beta(3,6)=      10.7125633476      

C     L=15
C     f(4,1)=           -410.3260290054
C     f(4,2)=           1168.8579712825
C     f(4,3)=         -98495.0510431024
C     f(4,4)=          98918.3033927164
C     f(4,5)=             18.2392467328
C     f(4,6)=           1105.3849889265
C     beta(4,1) =          48.1078229991
C     beta(4,2)  =         13.7459627225
C     beta(4,3)   =         4.9997421200
C     beta(4,4)    =        5.0077621877
C     beta(4,5)     =      89.3921544487
C     beta(4,6)      =     27.4857443437

C     L = 16
C     new mode Aug 2012 start at 2 peak at 5      
C     f(4,1)=           812.9918681787
C     f(4,2)=           441.2821316145
C     f(4,3)=          1301.6739693211
C     f(4,4)=          -143.4003781212
C     f(4,5)=            69.1750175418
C     f(4,6)=           175.9294532933
C     beta(4,1) =           48.1662482305
C     beta(4,2)  =           7.4414931962
C     beta(4,3)   =         15.8508455321
C     beta(4,4)    =         1.9564847414
C     beta(4,5)     =        1.7260534657
C     beta(4,6)      =       3.2147164196      

C     L=14
C     f(4,1)=             -93.7744100883
C     f(4,2)=            1229.3279885934
C     f(4,3)=          -22458.0945711871
C     f(4,4)=           91059.8513053717
C     f(4,5)=          -70480.8918194039
C     f(4,6)=            2864.4839104274
C     beta(4,1) =          45.0729200914
C     beta(4,2)  =         20.3745381016
C     beta(4,3)   =         6.4705187514
C     beta(4,4)    =        5.9188525245
C     beta(4,5)     =       5.8018896209
C     beta(4,6)      =      8.5158398479

C     L=16
C     taken: use old mode again
C     f(4,1)=     -359.4857268465
C     f(4,2)=  -718318.9196872473
C     f(4,3)=      952.9729049448
C     f(4,4)=   718848.3348118553
C     f(4,5)=       29.9449422503
C     f(4,6)=     1422.1605650097
C     beta(4,1)=    48.9404373087
C     beta(4,2)=     5.4911324770
C     beta(4,3)=    13.9191935499
C     beta(4,4)=     5.4924655928
C     beta(4,5)=    89.2258188514
C     beta(4,6)=    25.7560973761

C     L=32
C     taken: use old mode again
C     f(5,1)=     4285.5052331319
C     f(5,2)=  -600859.2054973751
C     f(5,3)=     1283.3234583067
C     f(5,4)=   601269.3489620079
C     f(5,5)=     -421.4303718400
C     f(5,6)=     3031.9399845112
C     beta(5,1)=  48.9404306370  
C     beta(5,2)=   5.1544511385  
C     beta(5,3)=  13.9220169737  
C     beta(5,4)=   5.1557936291  
C     beta(5,5)=  89.2258188514  
C     beta(5,6)=  25.7556569433  
      
C     L=32 
C     new mode Aug 2012 start at 2 peak at 5 
C     f(5,1)=   5336.9276606427
C     f(5,2)=    239.4345954380
C     f(5,3)=    696.6336721901
C     f(5,4)=    -55.6358364446
C     f(5,5)=  -2581.4422206295
C     f(5,6)=   2848.1479327744
C     beta(5,1)=   48.2614188891 
C     beta(5,2)=    4.1173201057 
C     beta(5,3)=    9.4612155568 
C     beta(5,4)=    2.4503125863 
C     beta(5,5)=   89.1104591707 
C     beta(5,6)=   21.1589243453       

C     L=64
C     taken: use old mode again
C     f(6,1)=      9938.2744323159
C     f(6,2)=   -500764.5439833679
C     f(6,3)=      1287.6877590235
C     f(6,4)=    501147.4432669706
C     f(6,5)=     14912.2836873664
C     f(6,6)=      3963.5336376826
C     beta(6,1)=   48.9368813776
C     beta(6,2)=    5.1508684107
c     beta(6,3)=   13.9644917853
C     beta(6,4)=    5.1524290713
C     beta(6,5)=   89.2282691400
C     beta(6,6)=   25.7556037045
      
C     L=64
C     new mode Aug 2012 start at 2 peak at 5 
C     f(6,1)=     9145.1487257103
C     f(6,2)=   -99402.7338129873
C     f(6,3)=     1370.4350614189
C     f(6,4)=    99735.9196457199
C     f(6,5)=    16595.1966646073
C     f(6,6)=     3784.6648254325
C     beta(6,1)=     46.9580300712
C     beta(6,2)=      4.9433214823
C     beta(6,3)=     13.6979549612
C     beta(6,4)=      4.9502924999
C     beta(6,5)=     88.8751266506
C     beta(6,6)=     25.9235159381     

C     L=12
C     f(6,1)=   -171.2655001298
C     f(6,2)=    856.1629846247
C     f(6,3)= -41925.1558800406
C     f(6,4)=  94908.0078132388
C     f(6,5)= -54688.5441954559
C     f(6,6)=   2646.2528663652
C     beta(6,1)=  45.1551416237
C     beta(6,2)=  19.3395079889
C     beta(6,3)=   6.3257807905
C     beta(6,4)=   5.9930084444
C     beta(6,5)=   5.8046579879
C     beta(6,6)=   8.4762240472

      DATA NI/6/
C     NL number of modes
C     NI mumber of coefficients per mode
C     We work in the KSM coordinate system

C     In this subroutine, index I=1-6 corresponds to each component of a mode
C     Index L = 1,7 denotes the seven L modes
C     Index j sums over observations (Br and Btheta) and is twice the number of
C     observations
C     Index M is used to obtain least squares equations (same dimension as L)
      PI     = 4.0 * datan(1.D0)
      twopi  = 2. * PI
      PI2    = PI  / 2.
      DEGREE = 180./ PI
      RADIAN = PI  / 180.
      drho   = 0.05
      dz     = 0.05
C     D is the half thickness of current sheet
      M  = 5
      M1 = 25
      M2 = 2 * M
C
C     END OF DEFINITIONS    
C ----------------------------------------------------------------------
      R_S3  = R1
      THETA = T1
      PHI   = P1
      rho   = R_S3*DSIN(THETA)
C     time  = time
C     RLT=(posLT(j)-12.0)*15.0
C     If (RLT .LT. 0.0) RLT=RLT+360.0
C     RLT=RLT*Radian

      XS3 = R_S3 * DSIN(THETA) * DCOS(PHI)
      YS3 = R_S3 * DSIN(THETA) * DSIN(PHI)
      ZS3 = R_S3 * DCOS(THETA)

C     NOW ROTATE the trajectory INTO KSM COORDINATES
C     write(6,*)"Xs3: ",Xs3
C     write(6,*)"Ys3: ",Ys3
C     write(6,*)"Zs3: ",Zs3
C     write(6,*)"time: ",time,epoch

      pvecin(1) = xs3
      pvecin(2) = ys3
      pvecin(3) = zs3 
      CALL KROT('S3C','DIS',pvecin,pvecou,time,epoch)
      XDIS = pvecou(1)
      YDIS = pvecou(2)
      ZDIS = pvecou(3)

      RLT = datan2(YDIS,XDIS)

C     Start calculating the field for each unit mode  at the mapped location
C     First map the trajectory

C     write(6,*)"XDIS: ",XDIS
C     write(6,*)"YDIS: ",YDIS
C     write(6,*)"ZDIS: ",ZDIS

      CALL mapit(time,XDIS,YDIS,ZDIS,XMAP,YMAP,ZMAP,r0or,Epoch)

      Z = ZMAP
      RHOMAG = sqrt(XMAP**2 + YMAP**2)

      DO 2 L = 1, NL
      ZM = dabs(Z-dz)
      IF (ZM < D) ZM= 0.5*(ZM**2/D+D)
      ZP = dabs(Z+dz)
      IF (ZP < D) ZP= 0.5*(ZP**2/D+D)
      xlpp = 0.
      xlpm = 0.
      DO 3 i=1, NI
          S1p=dsqrt((beta(L,i)/r0or+ZP)**2+(RHOMAG+Beta(L,i)/r0or)**2)
          S2p=dsqrt((beta(L,i)/r0or+ZP)**2+(RHOMAG-Beta(L,i)/r0or)**2)
          S1m=dsqrt((beta(L,i)/r0or+ZM)**2+(RHOMAG+Beta(L,i)/r0or)**2)
          S2m=dsqrt((beta(L,i)/r0or+ZM)**2+(RHOMAG-Beta(L,i)/r0or)**2)
          tp=2*(Beta(L,i)/r0or)/(S1p+S2p)
          tm=2*(Beta(L,i)/r0or)/(S1m+S2m)
          AAp=tp*dsqrt(1.-tp**2)/(S1p*S2p)
          AAm=tm*dsqrt(1.-tm**2)/(S1m*S2m)
          xlpp=xlpp+f(L,i)*AAp*rhomag
          xlpm=xlpm+f(L,i)*AAm*rhomag
    3 CONTINUE
      dxpldz=(xlpp-xlpm)/(2.*dz)
C     write(6,*)"dxpldz: ",dxpldz
      Xrho(L)=-dxpldz
    2 CONTINUE

C     NOW CALCULATE THE BZ COMPONENT
      DO 20 L = 1, NL
      rhom = RHOMAG-drho
      rhop = RHOMAG+drho

      xi = dabs(Z)
      IF (dabs(Z) .le. D) xi = 0.5*(Z**2/D+D)

      xlpp=0.
      xlpm=0.
      DO 30 i=1,NI
          S1p=dsqrt((beta(L,i)/r0or+xi)**2+(rhop+Beta(L,i)/r0or)**2)
          S2p=dsqrt((beta(L,i)/r0or+xi)**2+(rhop-Beta(L,i)/r0or)**2)
          S1m=dsqrt((beta(L,i)/r0or+xi)**2+(rhom+Beta(L,i)/r0or)**2)
          S2m=dsqrt((beta(L,i)/r0or+xi)**2+(rhom-Beta(L,i)/r0or)**2)
          tp=2*(Beta(L,i)/r0or)/(S1p+S2p)
          tm=2*(Beta(L,i)/r0or)/(S1m+S2m)
          AAp=tp*dsqrt(1.-tp**2)/(S1p*S2p)
          AAm=tm*dsqrt(1.-tm**2)/(S1m*S2m)
          xlpp=xlpp+f(L,i)*AAp*rhop
          xlpm=xlpm+f(L,i)*AAm*rhom
   30 CONTINUE
      dxpldr=(rhop*xlpp-rhom*xlpm)/(2.*drho)
C     write(6,*)"dxpldr: ",dxpldr
      Xz(L)=dxpldr/RHOMAG
C     write(6,*)"RHOMAG: ",RHOMAG
   20 CONTINUE

C     Now calculate the new mapped field for each  UNIT mode
      Bx=0.0
      By=0.0
      Bz=0.0

C     DO 60 L=5, 5
      DO 60 L=1, NL
        IF(( ymap .eq. 0.) .and. (xmap .eq. 0.)) THEN
          phimap=0.
          GO TO 69
        END IF
        phimap=atan2(YMAP,XMAP)
   69 Bx1 = Xrho(L) * Dcos(phimap)
      By1 = Xrho(L) * Dsin(phimap)
      Bz1 = Xz(L)

C     write(6,*)"Bx1: ",Bx1
C     write(6,*)"By1: ",By1
C     write(6,*)"Bz1: ",Bz1

C     Now Add the ACTUAL shield field to the current sheet field
      XMAP1 = XMAP 
      YMAP1 = YMAP 
      ZMAP1 = ZMAP 
      CALL shieldfield_csheet(XMAP1,YMAP1,ZMAP1,L,M,M1,M2,
     ~BX2,BY2,BZ2,r0or)

C     BX2=0.0
C     BY2=0.0
C     BZ2=0.0

      Bx=Bx+C(L)*r0or**dql*(Bx1*(1/(r0or))**2.0-Bx2/r0or)
      By=By+C(L)*r0or**dql*(By1*(1/(r0or))**2.0-By2/r0or)
      Bz=Bz+C(L)*r0or**dql*(Bz1*(1/(r0or))**2.0-Bz2/r0or)

C     -------------------------------------
C     Writing out single shielded tail modes
      Bxi = C(L)*r0or**dql*(Bx1*(1/(r0or))**2.0-Bx2/r0or)
      Byi = C(L)*r0or**dql*(By1*(1/(r0or))**2.0-By2/r0or)
      Bzi = C(L)*r0or**dql*(Bz1*(1/(r0or))**2.0-Bz2/r0or)
  
      CALL Mapped_field(
     ~  time,XMAP,YMAP,ZMAP,Bxi,Byi,Bzi,Bxdisi,Bydisi,Bzdisi,r0or,Epoch)
      bvecin(1) = Bxdisi
      bvecin(2) = Bydisi
      bvecin(3) = Bzdisi

C     write(6,*)"Bxdis: ",Bxdis
C     write(6,*)"Bydis: ",Bydis
C     write(6,*)"Bzdis: ",Bzdis

C     Now rotate back into spherical coordinate system
      CALL KROT('DIS','S3C',bvecin,bvecou,time,epoch)
      Bxs3i = bvecou(1)
      Bys3i = bvecou(2)
      Bzs3i = bvecou(3)

C     write(6,*)"Bxs3: ",Bxs3
C     write(6,*)"Bys3: ",Bys3
C     write(6,*)"Bzs3: ",Bzs3

      Brm1arr(L) =  Bxs3i*sin(theta)*cos(phi)+Bys3i*sin(theta)*sin(phi)+
     ~              Bzs3i*cos(theta)
      Btm1arr(L) =  Bxs3i*cos(theta)*cos(phi)+Bys3i*cos(theta)*sin(phi)-
     ~              Bzs3i*sin(theta)
      Bpm1arr(L) = -Bxs3i*sin(phi)+Bys3i*cos(phi)
C     End for writing out single tail modes
C     -------------------------------------

C     write(6,*)"Bx: ",Bx
C     write(6,*)"By: ",By
C     write(6,*)"Bz: ",Bz

   60 CONTINUE

C     write(6,*)"XMAP: ",XMAP
C     write(6,*)"YMAP: ",YMAP
C     write(6,*)"ZMAP: ",ZMAP


      CALL Mapped_field(
     ~  time,XMAP,YMAP,ZMAP,Bx,By,Bz,Bxdis,Bydis,Bzdis,r0or,Epoch)
C
      bvecin(1) = Bxdis
      bvecin(2) = Bydis
      bvecin(3) = Bzdis

C     write(6,*)"Bxdis: ",Bxdis
C     write(6,*)"Bydis: ",Bydis
C     write(6,*)"Bzdis: ",Bzdis

C     Now rotate back into spherical coordinate system
      CALL KROT('DIS','S3C',bvecin,bvecou,time,epoch)

      Bxs3 = bvecou(1)
      Bys3 = bvecou(2)
      Bzs3 = bvecou(3)

C     write(6,*)"Bxs3: ",Bxs3
C     write(6,*)"Bys3: ",Bys3
C     write(6,*)"Bzs3: ",Bzs3

      Brm1 =  Bxs3*sin(theta)*cos(phi) + Bys3*sin(theta)*sin(phi) +
     ~        Bzs3*cos(theta)
      Btm1 =  Bxs3*cos(theta)*cos(phi) + Bys3*cos(theta)*sin(phi) -
     ~        Bzs3*sin(theta)
      Bpm1 = -Bxs3*sin(phi)            + Bys3*cos(phi)

C     write(6,*)"Brm1: ",Brm1
C     write(6,*)"Btm1: ",Btm1
C     write(6,*)"Bpm1: ",Bpm1

      RETURN
      END SUBROUTINE shielded_csheetfield

C***********************************************************************

      SUBROUTINE shieldfield_csheet(X,Y,Z,L,M,M1,M2,
     ~Bxm,Bym,Bzm,r0or)
      IMPLICIT REAL*8(A-H,O-Z)
      REAL*8,  INTENT(IN)  :: X, Y, Z
      INTEGER, INTENT(IN)  :: L, M, M1, M2
      REAL*8,  INTENT(IN)  :: r0or
      REAL*8,  INTENT(OUT) :: Bxm,Bym,Bzm
      DIMENSION ABCD(6,25), PQ(6,10)
      DIMENSION a(M,M), p(M2)
      DIMENSION b(M,M), c(M,M), d(M,M)
      DIMENSION q(M2), r(M2), s(M2)
      DIMENSION XXP(M1), XXM(M1)

C     MODE 1 Aug 2012 peak at 5
C     Data (abcd (1, j), j = 1, 25) /
C     ~    -0.35153006639006954571D-05,
C     ~     0.23575800697128829436D-03,
C     ~    -0.68523032683502002604D-02,
C     ~     0.50584160772160019492D-01,
C     ~    -0.75694100518610099582D-01,
C     ~     0.22456816943470985577D-03,
C     ~    -0.10445488138626792506D-02,
C     ~     0.80112778867057421461D-01,
C     ~    -0.64183468172852997657D+00,
C     ~     0.98809223559946512960D+00,
C     ~    -0.52762855732372444806D-02,
C     ~    -0.30779140045062824971D-01,
C     ~    -0.16265439828762690011D-01,
C     ~     0.14316815615149347529D+00,
C     ~    -0.43849234995754932109D+00,
C     ~     0.22216049621179477513D-01,
C     ~     0.13161835247679370964D+00,
C     ~     0.99049377626076032244D-01,
C     ~     0.26013968527285968868D+00,
C     ~     0.32873575249230019280D+00,
C     ~    -0.17098828184387373369D-01,
C     ~    -0.10079006634277780374D+00,
C     ~    -0.11356834871457452407D+00,
C     ~    -0.29911921101585359217D+00,
C     ~    -0.48562354745533009570D+00/
C       Data (pq (1, j), j = 1, 10) /
C     ~        0.40000081062316894531D+01,
C     ~        0.80000104904174804688D+01,
C     ~        0.16000005722045898438D+02,
C     ~        0.32000000000000000000D+02,
C     ~        0.64000000000000000000D+02,
C     ~        0.39998939037322998047D+01,
C     ~        0.79999985694885253906D+01,
C     ~        0.16000000000000000000D+02,
C     ~        0.32000000000000000000D+02,
C     ~        0.64000000000000000000D+02/

C     MODE 2 Aug 2012 peak at 5
C     Data (abcd (2, j), j = 1, 25) /
C     ~     0.10678961640307714712D-04,
C     ~    -0.79277805987504465154D-03,
C     ~     0.20251500952258305599D-01,
C     ~    -0.14263622699416847461D+00,
C     ~     0.21053465945964125461D+00,
C     ~    -0.18532209610370927204D-03,
C     ~     0.53288622109895090323D-02,
C     ~    -0.23662576366437076958D+00,
C     ~     0.18459082306133383078D+01,
C     ~    -0.28062869665645946426D+01,
C     ~     0.36776835832946536305D-02,
C     ~     0.34357790313751211675D-01,
C     ~    -0.23388070299514979178D-01,
C     ~     0.33422041691935078234D-02,
C     ~     0.33369653661920106735D+00,
C     ~    -0.14215363077573411194D-01,
C     ~    -0.11038974961262414121D+00,
C     ~    -0.98945745397397172338D-01,
C     ~    -0.52453160708230328613D+00,
C     ~    -0.27580417633032239255D+00,
C     ~     0.10387043558141385857D-01,
C     ~     0.68529221559159858579D-01,
C     ~     0.26274248425618644542D-01,
C     ~    -0.35336925055235512971D+00,
C     ~     0.22982151697208269558D+00/
C     Data (pq (2, j), j = 1, 10) /
C     ~        0.39999945163726806641D+01,
C     ~        0.79999995231628417969D+01,
C     ~        0.16000001907348632812D+02,
C     ~        0.32000000000000000000D+02,
C     ~        0.64000000000000000000D+02,
C     ~        0.39993464946746826172D+01,
C     ~        0.80000047683715820313D+01,
C     ~        0.15999997138977050781D+02,
C     ~        0.32000000000000000000D+02,
C     ~        0.64000000000000000000D+02/

C     MODE 3 Aug 2012 peak at 5
C     Data (abcd (3, j), j = 1, 25) /
C     ~     0.56805659376579659158D-04,
C     ~    -0.49826416171172595304D-02,
C     ~     0.12526839464659417223D+00,
C     ~    -0.87409583160880188757D+00,
C     ~     0.12856485167171063377D+01,
C     ~     0.11216145207766179542D-02,
C     ~     0.35698754430595128206D-01,
C     ~    -0.11826555719375608788D+01,
C     ~     0.91620587898341590005D+01,
C     ~    -0.13838867662866388031D+02,
C     ~    -0.33466353340866311639D-01,
C     ~    -0.94500618941237457538D-01,
C     ~    -0.48684028075674673497D+00,
C     ~     0.14607883360870395073D+01,
C     ~    -0.13643921002377688012D+01,
C     ~     0.15327105053855760719D+00,
C     ~     0.77253223257713465877D+00,
C     ~     0.66066721067945421986D+00,
C     ~    -0.86601218815076608237D+00,
C     ~     0.10860893867610756924D+01,
C     ~    -0.12297786782076687573D+00,
C     ~    -0.72241771816889999158D+00,
C     ~    -0.77546966054031540860D+00,
C     ~    -0.22509195030543880378D+01,
C     ~    -0.19436733171165148093D+01/
C     Data (pq (3, j), j = 1, 10) /
C     ~        0.39999673366546630859D+01,
C     ~        0.79999799728393554688D+01,
C     ~        0.16000000000000000000D+02,
C     ~        0.32000000000000000000D+02,
C     ~        0.64000000000000000000D+02,
C     ~        0.39958155155181884766D+01,
C     ~        0.80000371932983398438D+01,
C     ~        0.15999984741210937500D+02,
C     ~        0.32000000000000000000D+02,
C     ~        0.64000000000000000000D+02/

C     MODE 4 Aug 2012 peak at 5
C     Data (abcd (4, j), j = 1, 25) /
C     ~     0.19638341783294487237D-03,
C     ~    -0.19807955436661750515D-01,
C     ~     0.50533167084825092452D+00,
C     ~    -0.35380207109579977498D+01,
C     ~     0.52069369271108874031D+01,
C     ~     0.83297371091017383127D-02,
C     ~     0.12861211658879426190D+00,
C     ~    -0.37415022082172026963D+01,
C     ~     0.28872340578355697005D+02,
C     ~    -0.43434963349345771633D+02,
C     ~    -0.22717934357869418172D+00,
C     ~    -0.81333065020493922059D+00,
C     ~    -0.26455667944533032276D+01,
C     ~     0.11153823484896403073D+02,
C     ~    -0.16209356038027308955D+02,
C     ~     0.10139287180009441958D+01,
C     ~     0.51585299887162561205D+01,
C     ~     0.42923466779422323469D+01,
C     ~    -0.21809683187273272154D+01,
C     ~     0.14371299951134901107D+02,
C     ~    -0.80264576829420197779D+00,
C     ~    -0.44983967922243719428D+01,
C     ~    -0.42669104398308794757D+01,
C     ~    -0.83595341151041111516D+01,
C     ~    -0.16784350456480702007D+02/
C     Data (pq (4, j), j = 1, 10) /
C     ~        0.39999535083770751953D+01,
C     ~        0.79999637603759765625D+01,
C     ~        0.16000024795532226562D+02,
C     ~        0.32000011444091796875D+02,
C     ~        0.64000007629394531250D+02,
C     ~        0.39839067459106445312D+01,
C     ~        0.80001392364501953125D+01,
C     ~        0.15999942779541015625D+02,
C     ~        0.31999998092651367188D+02,
C     ~        0.64000000000000000000D+02/

C     MODE 5 Aug 2012 peak at 5
C     Data (abcd (5, j), j = 1, 25) /
C     ~     0.57602181502077927514D-08,
C     ~    -0.27198178058757190846D-06,
C     ~     0.14413677102796104145D-04,
C     ~    -0.11014990050355253606D-03,
C     ~     0.16067283476392421300D-03,
C     ~     0.13193945273706952556D-05,
C     ~    -0.89828918286282105447D-04,
C     ~     0.21515136058656212509D-01,
C     ~    -0.19682427702020968696D+00,
C     ~     0.32279192414489094043D+00,
C     ~    -0.67184429268400664261D-04,
C     ~    -0.14284490361830121802D-01,
C     ~    -0.35629192133510734219D+00,
C     ~     0.30705765682217287527D+01,
C     ~    -0.61419772176518296902D+01,
C     ~     0.38732153569629476468D-03,
C     ~     0.11458722673551952065D+00,
C     ~     0.17929207140451239866D+01,
C     ~    -0.89188944523992432778D+01,
C     ~     0.32228805685562107897D+02,
C     ~    -0.32942021863423845375D-03,
C     ~    -0.11155626404876511915D+00,
C     ~    -0.19432724644718772922D+01,
C     ~     0.46048240963900362033D+01,
C     ~    -0.52986016965987758454D+02/
C     Data (pq (5, j), j = 1, 10) /
C     ~        0.20007674694061279297D+01,
C     ~        0.40067391395568847656D+01,
C     ~        0.80072870254516601563D+01,
C     ~        0.16005931854248046875D+02,
C     ~        0.32003528594970703125D+02,
C     ~        0.19872297048568725586D+01,
C     ~        0.39964907169342041016D+01,
C     ~        0.80006351470947265625D+01,
C     ~        0.16000080108642578125D+02,
C     ~        0.32000152587890625000D+02/

C     MODE 6 Aug 2012 peak at 5
C     Data (abcd (6, j), j = 1, 25) /
C     ~     0.12600037039950878022D-03,
C     ~    -0.10286886633016696646D-01,
C     ~    -0.12772309288678795891D+00,
C     ~     0.19989257632239989526D+01,
C     ~    -0.34321939046075433666D+01,
C     ~     0.55499013547815743974D-01,
C     ~     0.31173553004160692304D-01,
C     ~     0.61195985306330076270D+01,
C     ~    -0.52035160741602140888D+02,
C     ~     0.82976563376006623685D+02,
C     ~    -0.13749201131347199567D+01,
C     ~    -0.53997441327526773946D+01,
C     ~    -0.33885918241222711345D+02,
C     ~     0.22880910222910569018D+03,
C     ~    -0.38340140380148346821D+03,
C     ~     0.59693886907215576443D+01,
C     ~     0.31333698356997615520D+02,
C     ~     0.59768185557434726718D+02,
C     ~    -0.17342429700241308410D+03,
C     ~     0.45716000890957502634D+03,
C     ~    -0.46705148472034592189D+01,
C     ~    -0.26291656938929079956D+02,
C     ~    -0.47428705936755726214D+02,
C     ~     0.50799236585669959254D+02,
C     ~    -0.36579301953188536345D+03/
C     Data (pq (6, j), j = 1, 10) /
C     ~        0.40019817352294921875D+01,
C     ~        0.80021457672119140625D+01,
C     ~        0.16002088546752929688D+02,
C     ~        0.32001159667968750000D+02,
C     ~        0.64000549316406250000D+02,
C     ~        0.39563248157501220703D+01,
C     ~        0.80006628036499023438D+01,
C     ~        0.15999974250793457031D+02,
C     ~        0.32000030517578125000D+02,
C     ~        0.64000007629394531250D+02/

C     MODE 1 Peak at 6
C     Data (abcd (1, j), j = 1, 25) / 
C     ~    0.57900389912988492469D-05,
C     ~   -0.33867095260193265104D-02,
C     ~   -0.17885402929184579079D+00,
C     ~   -0.76105861802226257850D-01,
C     ~   -0.21707320837598329532D-01,
C     ~   -0.22806696336541141256D-04,
C     ~    0.13605855599644953723D-01,
C     ~    0.72655745785709626716D+00,
C     ~    0.33979077925435658968D+00,
C     ~    0.84015518189067677212D-01,
C     ~    0.38251019250246102387D-04,
C     ~   -0.23331745620323913747D-01,
C     ~   -0.12632616230609083896D+01,
C     ~   -0.66501059469106635901D+00,
C     ~   -0.11557383129450049530D+00,
C     ~   -0.46761660834837321942D-04,
C     ~    0.29188863483653651798D-01,
C     ~    0.16062192067320260058D+01,
C     ~    0.10008527972927729887D+01,
C     ~    0.96874076440719001368D-01,
C     ~    0.25527293306380216542D-04,
C     ~   -0.16076374426793411664D-01,
C     ~   -0.89077972230583615242D+00,
C     ~   -0.61419050551211071820D+00,
C     ~   -0.48351231840305022302D+00/

C     MODE 1 Peak at 4 (valid before Aug 2012)
C     Data (abcd (1, j), j = 1, 25) /
C     ~   -0.14004329253036342350D-07,
C     ~    0.37688393599679170797D-05,
C     ~    0.49241640240943258532D-04,
C     ~    0.27365751661153279172D-03,
C     ~   -0.29637572764194883845D-03,
C     ~    0.23608493437412536586D-06,
C     ~   -0.99590524882565070186D-04,
C     ~   -0.13716499714981711388D-02,
C     ~   -0.68651650806999411358D-02,
C     ~    0.13086363094019834996D-01,
C     ~   -0.17741658850262338020D-05,
C     ~    0.95936069671529971003D-03,
C     ~    0.14635178751353254966D-01,
C     ~    0.13851244022798143618D-01,
C     ~   -0.50732363273520704183D-01,
C     ~    0.59728889093050563374D-05,
C     ~   -0.35847147114331114892D-02,
C     ~   -0.60417634422173911445D-01,
C     ~   -0.24014471471244385192D-02,
C     ~   -0.87454103956755879778D-02,
C     ~   -0.44206961724389151058D-05,
C     ~    0.27227339013033198256D-02,
C     ~    0.47812648728699302935D-01,
C     ~   -0.49208189730297080544D-02,
C     ~    0.84257347502468000755D-02/

C     PEAK at 8 start at 5 (new mode Aug 2012)     
C      Data (abcd (1, j), j = 1, 25) /
C     ~   -0.22806650633874154589D-04,
C     ~   0.19981441572494572541D-02 ,
C     ~  -0.50917717624770160445D-01 ,
C     ~   0.35714413408705841801D+00 ,
C     ~  -0.52608610022227830783D+00 ,
C     ~  -0.39669626460951413474D-03 ,
C     ~  -0.14042301137871120859D-01 ,
C     ~   0.47028862177184993509D+00 ,
C     ~  -0.36578470766660089453D+01 ,
C     ~   0.55325325060575014646D+01 ,
C     ~   0.12202297433719602737D-01 ,
C     ~   0.35921039006745036914D-01 ,
C     ~   0.19367172629944376383D+00 ,
C     ~  -0.47038727079262165898D+00 ,
C     ~   0.31050626588489294022D+00 ,
C     ~  -0.56170941159746065607D-01 ,
C     ~  -0.29912370171159619048D+00 ,
C     ~  -0.30322589282814016443D+00 ,
C     ~   0.27108065541804371446D+00 ,
C     ~  -0.34160605247767883430D+00 ,
C     ~   0.45136745871513728190D-01 ,
C     ~   0.27576166498336385224D+00 ,
C     ~   0.23869604996261831187D+00 ,
C     ~   0.18925998312564620840D+00 ,
C     ~   0.44434492288249283209D+00 /
C       
      Data (abcd (1, j), j = 1, 25) /
     ~   0.52034884821428676664D-08,
     ~  -0.68983393140626007927D-06 ,
     ~   0.26013060818085627384D-04 ,
     ~  -0.20220185555762154508D-03 ,
     ~   0.30564180727029997440D-03 ,
     ~   0.19930189328513339478D-06 ,
     ~  -0.12444044081689495742D-04 ,
     ~   0.14540867652184175443D-02 ,
     ~  -0.12059368770663559875D-01 ,
     ~   0.18972457961433673856D-01 ,
     ~  -0.12712315070435666379D-04 ,
     ~  -0.10948713314089635051D-02 ,
     ~  -0.11325779706928780499D-01 ,
     ~   0.10236958766108690544D+00 ,
     ~  -0.19081140235653648585D+00 ,
     ~   0.78157561366167476293D-04 ,
     ~   0.81610134021544745997D-02 ,
     ~   0.34775577517361599966D-01 ,
     ~  -0.13775781931913483325D+00 ,
     ~   0.59210130685600637079D+00 ,
     ~  -0.67653598302448326595D-04 ,
     ~  -0.77757671746549184011D-02 ,
     ~  -0.43300045770593376115D-01 ,
     ~  -0.72121864531289603395D-01 ,
     ~  -0.10966408811380061739D+01 /                   

      Data (abcd (2, j), j = 1, 25) /
     ~   0.12545399310341492571D-08,
     ~  -0.11363771161960078800D-06,
     ~   0.55343644774000541920D-05,
     ~  -0.44032796493657441067D-04,
     ~   0.66582238656644651925D-04,
     ~   0.55943482168823719799D-07,
     ~  -0.97499376833970518089D-05,
     ~   0.18348839506297323111D-02,
     ~  -0.17055331333206070338D-01,
     ~   0.27862314450439854691D-01,
     ~  -0.49597135329461795178D-05,
     ~  -0.11392557185046436278D-02,
     ~  -0.24260911596492629183D-01,
     ~   0.24245633528367632747D+00,
     ~  -0.49494927634042140685D+00,
     ~   0.33779783670884271877D-04,
     ~   0.10006602582479727503D-01,
     ~   0.12919577290627834443D+00,
     ~  -0.61024568956482516668D+00,
     ~   0.27989901087956345904D+01,
     ~  -0.30128015057845978483D-04,
     ~  -0.10625414457624882125D-01,
     ~  -0.18652399199560046328D+00,
     ~  -0.47424636188494684319D-01,
     ~  -0.61415416903532733528D+01/
     
C     -----------------------------------
C     New mode Aug 2012
C     MODE 1 Start at 2; Peak at 6
C     -----------------------------------
C     Data (abcd (1, j), j = 1, 25) /
C     ~  -0.13774757805365866521D-05,
C     ~  -0.70971199764301553040D-03,
C     ~   0.27634465364814016219D-01,
C     ~  -0.22098384462359346725D+00,
C     ~   0.33708836353742760439D+00,
C     ~   0.10832025970333686131D-02,
C     ~  -0.53530485360866290889D-02,
C     ~   0.18944184087282156526D+00,
C     ~  -0.12716464964495621803D+01,
C     ~   0.18215778704469625104D+01,
C     ~  -0.26286979174857352803D-01,
C     ~  -0.48444880496650646884D-01,
C     ~  -0.52660660932882952245D+00,
C     ~   0.45082237972458321096D+01,
C     ~  -0.65842130925306410205D+01,
C     ~   0.11411754551073770292D+00,
C     ~   0.31979110223247919631D+00,
C     ~   0.19856752256512044674D+00,
C     ~  -0.49506264985306094317D+01,
C     ~   0.77412277628375987248D+01,
C     ~  -0.89074049080843556436D-01,
C     ~  -0.26584905850947632144D+00,
C     ~   0.29268786738959606808D-02,
C     ~   0.22086426652376278668D+01,
C     ~  -0.47584288669132392968D+01/

C     -----------------------------------
C     New mode Aug 2012
C     MODE 1 Start at 1.2; Peak at 5
C     -----------------------------------
C     Data (abcd (1, j), j = 1, 25) /
C     ~    0.27520534657227329261D-06,
C     ~   -0.56479711896642841723D-03,
C     ~    0.20710029128770408857D-01,
C     ~   -0.16323345825271243226D+00,
C     ~    0.24811642320274354678D+00,
C     ~    0.71681996468643278259D-03,
C     ~   -0.30471710664889337106D-02,
C     ~    0.10516844777484216089D+00,
C     ~   -0.68502745701355149954D+00,
C     ~    0.96848183576033330944D+00,
C     ~   -0.17509138573501854125D-01,
C     ~   -0.30505093121171573262D-01,
C     ~   -0.36239860918944366563D+00,
C     ~    0.30638637104303150238D+01,
C     ~   -0.44606235835225174213D+01,
C     ~    0.76215353970331842226D-01,
C     ~    0.21157611817542446020D+00,
C     ~    0.14114391054107056078D+00,
C     ~   -0.33790814687815076134D+01,
C     ~    0.52659027560436868498D+01,
C     ~   -0.59563477929664765398D-01,
C     ~   -0.17737804165846532412D+00,
C     ~   -0.45599955713709894312D-03,
C     ~    0.15056277265595638948D+01,
C     ~   -0.32357335890645630982D+01/

C     MODE 2
C     Data (abcd (2, j), j = 1, 25) /
C     ~    0.77668182761180606377D-08,
C     ~   -0.20196291280988356575D-05,
C     ~   -0.49273322968357806672D-04,
C     ~    0.21539440444312103473D-03,
C     ~   -0.83705446882959897436D-04,
C     ~    0.33536606594698534777D-04,
C     ~   -0.12866462057380196881D-01,
C     ~   -0.31621223228494135248D+00,
C     ~    0.61757501135689674143D+00,
C     ~   -0.11476830385717948779D+00,
C     ~   -0.17431745682254278229D-04,
C     ~    0.65663659425452749474D-02,
C     ~    0.15888921816088288352D+00,
C     ~   -0.33062421999133104755D+00,
C     ~    0.13245199895207815377D+00,
C     ~   -0.26812881293819263994D-04,
C     ~    0.10672183150391700845D-01,
C     ~    0.27287667133660016283D+00,
C     ~   -0.42942834594719494489D+00,
C     ~   -0.32044936014889375819D+00,
C     ~    0.10700204255233096706D-04,
C     ~   -0.43707712411099670646D-02,
C     ~   -0.11613927454426622443D+00,
C     ~    0.71015669683247804044D-01,
C     ~   -0.13462690359270563789D+01/
 
C     MODE 2
      Data (abcd (3, j), j = 1, 25) /
     ~   0.24824112592626616873D-08,
     ~  -0.22636768136616835943D-06,
     ~   0.11627747929157943881D-04,
     ~  -0.93185828541469004570D-04,
     ~   0.14182203142413209369D-03,
     ~  -0.22688551303297704850D-06,
     ~  -0.46061637003510315412D-04,
     ~   0.22179188550876874353D-02,
     ~  -0.21306193978350578683D-01,
     ~   0.34153314431143159957D-01,
     ~   0.57496428715017037080D-05,
     ~   0.97441205991077970046D-03,
     ~  -0.18779683259496859904D-01,
     ~   0.32105889582564206552D+00,
     ~  -0.60933710405131935595D+00,
     ~  -0.15042258847903980102D-04,
     ~  -0.20826093187111503757D-02,
     ~   0.73347474035392712377D-01,
     ~  -0.10606675688097344512D+01,
     ~   0.37511030714926558716D+01,
     ~   0.76517215572531172669D-05,
     ~  -0.16417219343249312910D-02,
     ~  -0.18174979954668427085D+00,
     ~   0.16023634510000572817D+00,
     ~  -0.87233263461601922018D+01/

C     MODE 3
C     Data (abcd (3, j), j = 1, 25) /
C     ~    0.92353191216470111868D-02,
C     ~   -0.47065512098254735917D+01,
C     ~   -0.22554133829293983026D+03,
C     ~   -0.23457137432873151894D+02,
C     ~    0.46111664280590520803D+02,
C     ~   -0.12239552063973810902D-01,
C     ~    0.62483152611269074938D+01,
C     ~    0.29997204397613383974D+03,
C     ~    0.31604341484015141539D+02,
C     ~   -0.57823703323961463951D+02,
C     ~    0.41170674967848803760D-02,
C     ~   -0.21219401709444536408D+01,
C     ~   -0.10293144491195860279D+03,
C     ~   -0.11947843043877295343D+02,
C     ~    0.14488934852695538602D+02,
C     ~   -0.18939743624590780868D-02,
C     ~    0.99324763584517850034D+00,
C     ~    0.49141758650519626883D+02,
C     ~    0.77481048808290280405D+01,
C     ~   -0.39887477569330695104D+01,
C     ~    0.78113976972219205663D-03,
C     ~   -0.41307234546100435323D+00,
C     ~   -0.20642069925114081563D+02,
C     ~   -0.40978583394053940125D+01,
C     ~   -0.35575024161355539575D+01/

C     MODE 3
      Data (abcd (4, j), j = 1, 25) /
     ~   0.27771830099459217676D-06,
     ~  -0.81341899879401334040D-04,
     ~   0.29282572672455508579D-02,
     ~  -0.22583395250099128071D-01,
     ~   0.34010285131856290985D-01,
     ~   0.99291569751668582507D-04,
     ~   0.26052436207605949789D-02,
     ~   0.10471966853985351098D+00,
     ~  -0.83263064611646853130D+00,
     ~   0.13428602164692611254D+01,
     ~  -0.36787553054767226091D-02,
     ~  -0.26316046069385012895D+00,
     ~  -0.93764343790019077929D+00,
     ~   0.50555823544919160639D+01,
     ~  -0.11196287302734999614D+02,
     ~   0.18441461144297600916D-01,
     ~   0.16391989138877585130D+01,
     ~   0.46270302869078419050D+01,
     ~   0.18359289360326336293D+01,
     ~   0.32155291790796020734D+02,
     ~  -0.15031144448049583301D-01,
     ~  -0.14580925618602209592D+01,
     ~  -0.54630082835657294282D+01,
     ~  -0.12099960332237163385D+02,
     ~  -0.53861339558258123361D+02/
   
C     MODE 4
C     Data (abcd (4, j), j = 1, 25) /
C     ~    0.21992617916451395743D+02,
C     ~   -0.11399307544820429516D+05,
C     ~   -0.76667036088618019107D+06,
C     ~   -0.39185539354281475610D+06,
C     ~   -0.59239538019716295735D+03,
C     ~    0.47524457523937337910D+00,
C     ~   -0.24323272755420672908D+03,
C     ~   -0.16297746944208286734D+05,
C     ~   -0.78974933554931041612D+04,
C     ~   -0.34576773059694034806D+02,
C     ~   -0.67509882718360936237D+00,
C     ~    0.34607166925607515395D+03,
C     ~    0.23198896555695553978D+05,
C     ~    0.11312410520290614446D+05,
C     ~    0.45910310222559074233D+02,
C     ~   -0.21888994250814901576D+02,
C     ~    0.11346764342961157545D+05,
C     ~    0.76316098499373987706D+06,
C     ~    0.39024560870819433588D+06,
C     ~    0.58817978798588868016D+03,
C     ~    0.96230585996413289251D-01,
C     ~   -0.50295745279677515071D+02,
C     ~   -0.33917785686196171290D+04,
C     ~   -0.18056934409611232084D+04,
C     ~   -0.22829873613259561437D+02/

C     MODE 4
C     Data (abcd (5, j), j = 1, 25) /
C     ~  -0.45915862230592986108D-03,
C     ~   0.93854011964554895364D-01,
C     ~  -0.38314667487779807686D+01,
C     ~   0.30828311901267987594D+02,
C     ~  -0.47112777682193481610D+02,
C     ~   0.53309485417775998028D-01,
C     ~  -0.18620821471067905417D+00,
C     ~   0.24001863168454320885D+02,
C     ~  -0.21574676433475329418D+03,
C     ~   0.34315573083493438844D+03,
C     ~  -0.13148027970186495494D+01,
C     ~  -0.76415371208712992157D+01,
C     ~  -0.21148444432406684257D+02,
C     ~   0.23239199306603961759D+03,
C     ~  -0.48969517910810122885D+03,
C     ~   0.56621555612493432719D+01,
C     ~   0.35897685020618332885D+02,
C     ~   0.40435286626670226440D+02,
C     ~  -0.91416647341691259498D+01,
C     ~   0.50909741444276653510D+03,
C     ~  -0.44092277778652038123D+01,
C     ~  -0.28613448638516562283D+02,
C     ~  -0.44306202263644621553D+02,
C     ~  -0.53377514196322152884D+02,
C     ~  -0.47645787621277645485D+03/
 
C     Data (abcd (5, j), j = 1, 25) /    
C     ~   0.29184501829407924565D-03,
C     ~  -0.35174081943346383050D-01,
C     ~   0.82480246693769743160D+00,
C     ~  -0.55640837422291227199D+01,
C     ~   0.80920395915216154492D+01,
C     ~   0.33524297003114426896D-01,
C     ~   0.27659122634700405996D+00,
C     ~  -0.51499706442952035701D+01,
C     ~   0.37955245686153382678D+02,
C     ~  -0.55685624735124136464D+02,
C     ~  -0.86509761148266395292D+00,
C     ~  -0.37778516452277250970D+01,
C     ~  -0.74686459122406816391D+01,
C     ~   0.37523091820409682384D+02,
C     ~  -0.67750312112703170442D+02,
C     ~   0.37807149835021629514D+01,
C     ~   0.20453142129327371634D+02,
C     ~   0.15788458773079129926D+02,
C     ~   0.88840198407937549518D+01,
C     ~   0.59253145460556318369D+02,
C     ~  -0.29649080045840339181D+01,
C     ~  -0.17103608884133656431D+02,
C     ~  -0.15987670755539895140D+02,
C     ~  -0.31371844952042430066D+02,
C     ~  -0.71049048051729997155D+02/
 
      Data (abcd (5, j), j = 1, 25) /        
     ~  -0.99911155916475476429D-04,
     ~  -0.16497770784682154899D-02,
     ~  -0.43294948559606633021D+00,
     ~   0.44153048101428264971D+01,
     ~  -0.71275189851081943715D+01,
     ~   0.79586008001811162083D-01,
     ~   0.41014729592350235299D+00,
     ~   0.16633384478053871991D+01,
     ~  -0.22369799721501067324D+02,
     ~   0.41291100537732880582D+02,
     ~  -0.19782143403119247793D+01,
     ~  -0.11556993496696394530D+02,
     ~  -0.81831519476911562805D+01,
     ~   0.11663048191732711345D+02,
     ~  -0.74936198497502246596D+02,
     ~   0.84827948435163005314D+01,
     ~   0.55050326212560698025D+02,
     ~   0.47358311675180331690D+02,
     ~   0.14227065859524753932D+03,
     ~   0.27681529848104631242D+02,
     ~  -0.65980337848883907625D+01,
     ~  -0.44385714456835309250D+02,
     ~  -0.50397503517239158555D+02,
     ~  -0.12514882115042821908D+03,
     ~  -0.12233187994479561667D+03/  

C     L=32
C     Data (abcd (5, j), j = 1, 25)  / 
C     ~       0.38370706169389450224D+02,
C     ~      -0.32479962500340833209D+05,
C     ~      -0.94660043127081774372D+06,
C     ~      -0.56864021527375161468D+06,
C     ~      -0.27814993022411260703D+04,
C     ~       0.20600370116089345984D+00,
C     ~      -0.17025940328754780139D+03,
C     ~      -0.49395843669987780089D+04,
C     ~      -0.30000800342043403290D+04,
C     ~      -0.97079695288155960497D+01,
C     ~      -0.14579821620379142643D+01,
C     ~       0.12204100420071475330D+04,
C     ~       0.35488135357577208495D+05,
C     ~       0.21388387530527799285D+05,
C     ~       0.76407552757476420168D+02,
C     ~      -0.37682332024326319341D+02,
C     ~       0.31912372589853306159D+05,
C     ~       0.93015110639191078689D+06,
C     ~       0.55873046694301766734D+06,
C     ~       0.28055544715537447331D+04,
C     ~       0.56360430778218297831D+00,
C     ~      -0.48256080155038478807D+03,
C     ~      -0.14099258801109653127D+05,
C     ~      -0.84807776710395010866D+04,
C     ~      -0.13800109969955343203D+03/

C     L=64
      Data (abcd (6, j), j = 1, 25)  / 
     ~   -0.88135632848555722773D+01,
     ~    0.11514098078690224724D+05,
     ~    0.21560988124946391941D+06,
     ~    0.45832673059045454522D+05,
     ~    0.21572565347949170799D+03,
     ~    0.54760001008211061090D+00,
     ~   -0.70240562027308603987D+03,
     ~   -0.13144093980819422373D+05,
     ~   -0.29206319879254337479D+04,
     ~    0.16213317181975170910D+02,
     ~   -0.19639794571317899851D+03,
     ~    0.25945417987723824460D+06,
     ~    0.48620565446258918740D+07,
     ~    0.10126003837642234550D+07,
     ~    0.84424842477639039373D+04,
     ~    0.18558502462981778080D+03,
     ~   -0.24493795838444620827D+06,
     ~   -0.45896954969516352562D+07,
     ~   -0.95728945603260022778D+06,
     ~   -0.75585692609229706562D+04,
     ~    0.19078884327836625800D+02,
     ~   -0.25327914385527452445D+05,
     ~   -0.47482700183867816434D+06,
     ~   -0.98233311870422994616D+05,
     ~   -0.12493968682442764261D+04/

C     Peak at 6
C     Data (pq (1, j), j = 1, 10) /
C     ~       0.14671022415161132812D+02,
C     ~       0.17176361083984375000D+02,
C     ~       0.22174074172973634588D+02,
C     ~       0.38154026031494140625D+02,
C     ~       0.73457305908203123223D+02,
C     ~       0.11755837202072143554D+01,
C     ~       0.19926652908325195312D+01,
C     ~       0.38097207546234130859D+01,
C     ~       0.89784660339355468750D+01,
C     ~       0.32120651245117186611D+02/

C     MODE1 PEAK at 4 (valid before Aug 2012)
C     Data (pq (1, j), j = 1, 10) /
C     ~    0.67008705139160156250D+01,
C     ~    0.10257417678833007368D+02,
C     ~    0.17875274658203124111D+02,
C     ~    0.36698616027832029473D+02,
C     ~    0.73508399963378909802D+02,
C     ~    0.14384461641311645507D+01,
C     ~    0.27737910747528076171D+01,
C     ~    0.63922395706176757812D+01,
C     ~    0.35446662902832031250D+02,
C     ~    0.73347290039062498223D+02/

C     MODE1 PEAK at 8 strart at 5
C     Data (pq (1, j), j = 1, 10) /
C     ~   0.39999859333038330078D+01,
C     ~   0.79999923706054687500D+01,
C     ~   0.16000000000000000000D+02,
C     ~   0.32000000000000000000D+02,
C     ~   0.64000000000000000000D+02,
C     ~   0.39984464645385742188D+01,
C     ~   0.80000057220458984375D+01,
C     ~   0.15999995231628417969D+02,
C     ~   0.32000000000000000000D+02,
C     ~   0.64000000000000000000D+02/

      Data (pq (1, j), j = 1, 10) /
     ~   0.23783802986145019531D+01,
     ~   0.47568554878234863281D+01,
     ~   0.95136852264404296875D+01,
     ~   0.19027336120605468750D+02,
     ~   0.38054645538330078125D+02,
     ~   0.23773388862609863281D+01,
     ~   0.47568039894104003906D+01,
     ~   0.95136604309082031250D+01,
     ~   0.19027313232421875000D+02,
     ~   0.38054626464843750000D+02/

C     -----------------------------------
C     Aug 2012
C     MODE1 start at 2; PEAK at 6
C     -----------------------------------
C     Data (pq (1, j), j = 1, 10) /
C     ~   0.40000061988830566406D+01,
C     ~   0.80000143051147460938D+01,
C     ~   0.16000020980834960938D+02,
C     ~   0.32000015258789062500D+02,
C     ~   0.64000007629394531250D+02,
C     ~   0.39996228218078613281D+01,
C     ~   0.79999980926513671875D+01,
C     ~   0.16000001907348632812D+02,
C     ~   0.32000000000000000000D+02,
C     ~   0.64000000000000000000D+02/
     
C     -----------------------------------
C     Aug 2012
C     MODE1 start at 1.2; PEAK at 5
C     -----------------------------------
C     Data (pq (1, j), j = 1, 10) /
C     ~   0.40000061988830566406D+01,
C     ~   0.80000143051147460938D+01,
C     ~   0.16000020980834960938D+02,
C     ~   0.32000015258789062500D+02,
C     ~   0.64000007629394531250D+02,
C     ~   0.39996228218078613281D+01,
C     ~   0.79999980926513671875D+01,
C     ~   0.16000001907348632812D+02,
C     ~   0.32000000000000000000D+02,
C     ~   0.64000000000000000000D+02/     

C     MODE 2
C     Data (pq (2, j), j = 1, 10) / 
C     ~      0.62806210517883300781D+01,
C     ~      0.19055501937866210937D+02,
C     ~      0.16761466979980468750D+02,
C     ~      0.28816335678100584161D+02,
C     ~      0.77478836059570310723D+02,
C     ~      0.12525197267532348632D+01,
C     ~      0.21577339172363281250D+01,
C     ~      0.42458128929138183593D+01,
C     ~      0.10150924682617188388D+02,
C     ~      0.32478603363037108486D+02/
     
C     MODE 2
      Data (pq (2, j), j = 1, 10) /     
     ~   0.19998002052307128906D+01,
     ~   0.40004315376281738281D+01,
     ~   0.80004920959472656250D+01,
     ~   0.16000474929809570312D+02,
     ~   0.32000373840332031250D+02,
     ~   0.19961978197097778320D+01,
     ~   0.39991734027862548828D+01,
     ~   0.80000810623168945313D+01,
     ~   0.16000000000000000000D+02,
     ~   0.32000003814697265625D+02/     
     
C     MODE 2
      Data (pq (3, j), j = 1, 10) /         
     ~   0.19997159242630004883D+01,
     ~   0.40001826286315917969D+01,
     ~   0.80002622604370117188D+01,
     ~   0.16000364303588867188D+02,
     ~   0.32000404357910156250D+02,
     ~   0.19948312044143676758D+01,
     ~   0.39987049102783203125D+01,
     ~   0.80001277923583984375D+01,
     ~   0.16000005722045898438D+02,
     ~   0.32000026702880859375D+02/
      
C     MODE 3
C     Data (pq (3, j), j = 1, 10) / 
C     ~     -0.20120578765869141513D+02,
C     ~      0.20586240768432615411D+02,
C     ~      0.23876194000244139736D+02,
C     ~      0.38100570678710937500D+02,
C     ~      0.71865959167480468750D+02,
C     ~      0.11369531154632568359D+01,
C     ~      0.18996990919113159179D+01,
C     ~      0.35767178535461425781D+01,
C     ~      0.82503862380981445312D+01,
C     ~      0.30094644546508790838D+02/

C     MODE 3
      Data (pq (4, j), j = 1, 10) /
     ~  0.24642193317413330078D+01,
     ~  0.49311742782592773438D+01,
     ~  0.98530044555664062500D+01,
     ~  0.19699954986572265625D+02,
     ~  0.39397197723388671875D+02,
     ~  0.23889598846435546875D+01,
     ~  0.49206199645996093750D+01,
     ~  0.98495092391967773438D+01,
     ~  0.19698287963867187500D+02,
     ~  0.39396602630615234375D+02/
       
C     MODE 4
C     Data (pq (4, j), j = 1, 10) / 
C     ~      0.38784885406494140625D+02,
C     ~      0.26744506835937498223D+02,
C     ~      0.27683023452758788174D+02,
C     ~      0.38958595275878904473D+02,
C     ~      0.72836189270019531250D+02,
C     ~      0.11952087879180908203D+01,
C     ~      0.20218055248260498046D+01,
C     ~      0.38115730285644531250D+01,
C     ~      0.88158550262451171875D+01,
C     ~      0.33938060760498047763D+02/

C     MODE 4
C     Data (pq (5, j), j = 1, 10) /
C     ~   0.40072717666625976562D+01,
C     ~   0.80101032257080078125D+01,
C     ~   0.16006877899169921875D+02,
C     ~   0.32003021240234375000D+02,
C     ~   0.64001159667968750000D+02,
C     ~   0.39908657073974609375D+01,
C     ~   0.80005998611450195313D+01,
C     ~   0.16000095367431640625D+02,
C     ~   0.32000057220458984375D+02,
C     ~   0.63999942779541015625D+02/
     
C     Data (pq (5, j), j = 1, 10) /
C     ~   0.40006184577941894531D+01,
C     ~   0.80003767013549804688D+01,
C     ~   0.16000297546386718750D+02,
C     ~   0.32000118255615234375D+02,
C     ~   0.64000045776367187500D+02,
C     ~   0.39639279842376708984D+01,
C     ~   0.80005216598510742188D+01,
C     ~   0.15999864578247070312D+02,
C     ~   0.32000000000000000000D+02,
C     ~   0.63999992370605468750D+02/

      Data (pq (5, j), j = 1, 10) /                           
     ~   0.40074467658996582031D+01,
     ~   0.80055150985717773438D+01,
     ~   0.16002401351928710938D+02,
     ~   0.32000724792480468750D+02,
     ~   0.64000205993652343750D+02,
     ~   0.39702200889587402344D+01,
     ~   0.80015344619750976563D+01,
     ~   0.15999733924865722656D+02,
     ~   0.32000061035156250000D+02,
     ~   0.63999969482421875000D+02/

C     L=32
C     Data (pq (5, j), j = 1, 10) /
C     ~      0.38836956024169921875D+02,
C     ~      0.24311889648437499111D+02,
C     ~      0.29589153289794921875D+02,
C     ~      0.39441062927246095526D+02,
C     ~      0.74735847473144527697D+02,
C     ~      0.13389585018157958984D+01,
C     ~      0.23527803421020507812D+01,
C     ~      0.45696606636047363281D+01,
C     ~      0.10911170959472655894D+02,
C     ~      0.44013973236083980822D+02/

C     L=64
      Data (pq (6, j), j = 1, 10) /
     ~     0.36148849487304688388D+02,
     ~     0.26422071456909179687D+02,
     ~     0.52445343017578123223D+02,
     ~     0.50159130096435546875D+02,
     ~     0.72792282104492187500D+02,
     ~     0.14043723344802856445D+01,
     ~     0.26591603755950927734D+01,
     ~     0.55725831985473632812D+01,
     ~     0.15859733581542969638D+02,
     ~     0.70801498413085939276D+02/


C     Separate out the a. b, c, d and p, q.
      KA = 1
      DO 10 I = 1, M
        DO 10 J = 1, M
          a(i,j) = ABCD(L, KA)
          KA = KA + 1
  10  CONTINUE 
C     Now the PQs
      DO 14 I = 1, M
        J = I + M
        p(i) = pq(L, I) / r0or
        p(j) = pq(L, J) / r0or
  14  CONTINUE

C     Now compute the magnetic field from differentiating the potential U.
      CALL gradU(Bxm,Bym,Bzm,x,y,z,a,p,M,M1,M2) 
    
   3  CONTINUE ! NEVER USED (WARNING)
      RETURN
      END SUBROUTINE shieldfield_csheet

C***********************************************************************
      SUBROUTINE gradU(Bxm1,Bym1,Bzm1,x1,y1,z1,a,p,M,M1,M2)
      IMPLICIT REAL*8(A-H,O-Z)
      REAL*8,   INTENT(OUT) :: Bxm1, Bym1, Bzm1
      REAL*8,   INTENT(IN)  :: x1, y1, z1
      INTEGER,  INTENT(IN)  :: M, M1, M2
      REAL*8,   INTENT(IN)  :: a(M,M), p(M2)
C     REAL*8,  INTENT(IN)  :: a,p
C     Dimension a(M,M),p(M2)
C     Calculates the gradient of the potential.

      Bxm1 = 0.
      Bym1 = 0.
      Bzm1 = 0.
      DO 111 i = 1, M
      DO 222 k = M+1, M2
      J = K-M
C     Write(*,*)M,M1,M2,p(i),p(k)
      Term1 = dsqrt(1.D0/p(i)**2 + 1.D0/p(k)**2)
      Term2 = dexp(Term1*x1)
      Bxm1 = Bxm1 + Term1 * a(i,j) * term2 * dcos(y1/p(i))*dsin(z1/p(k))
      Bym1 = Bym1 -    a(i,j)/p(i) * term2 * dsin(y1/p(i))*dsin(z1/p(k))
      Bzm1 = Bzm1 +    a(i,j)/p(k) * term2 * dcos(y1/p(i))*dcos(z1/p(k))
  222 CONTINUE
  111 CONTINUE
      Bxm1 = -Bxm1
      Bym1 = -Bym1
      Bzm1 = -Bzm1
      RETURN
      END SUBROUTINE gradU

C***********************************************************************

      SUBROUTINE CAR2SPH_MAG(BX,BY,BZ,BR,BTH,BPHI,TH,PHI)
      IMPLICIT REAL*8(A-H,O-Z)
      REAL*8,  INTENT(IN)  :: BX,BY,BZ
      REAL*8,  INTENT(IN)  :: TH,PHI
      REAL*8,  INTENT(OUT) :: BR,BTH,BPHI
C     (ARFKEN)
      BR  = BX*dsin(TH)*dcos(PHI)+BY*dsin(TH)*dsin(PHI)+BZ*dcos(TH)
      BTH = BX*dcos(TH)*dcos(PHI)+BY*dcos(TH)*dsin(PHI)-BZ*dsin(TH)
      BPHI=-BX*dsin(PHI)+BY*dcos(PHI) 
      RETURN
      END SUBROUTINE CAR2SPH_MAG

C***********************************************************************

      SUBROUTINE SPH2CAR_MAG(BR,BTH,BPHI, BX,BY,BZ,TH,PHI)
      IMPLICIT REAL*8(A-H,O-Z)
      REAL*8,  INTENT(IN)  :: BR,BTH,BPHI
      REAL*8,  INTENT(IN)  :: TH,PHI
      REAL*8,  INTENT(OUT) :: BX,BY,BZ
C     (LIUTAURAS RUSAITIS 02.07.2018)
C     (ARFKEN)
      BX  = BR*dsin(TH)*dcos(PHI)+BTH*dcos(TH)*dcos(PHI)-BPHI*dsin(PHI)
      BY  = BR*dsin(TH)*dsin(PHI)+BTH*dcos(TH)*dsin(PHI)+BPHI*dcos(PHI)
      BZ  = BR*dcos(TH)-BTH*dsin(TH) 
      RETURN
      END SUBROUTINE SPH2CAR_MAG

C***********************************************************************

      SUBROUTINE getIMF_penetration(BY_IMF,BZ_IMF,BY_p,BZ_p)
      IMPLICIT REAL*8(A-H,O-Z)
      REAL*8,  INTENT(IN)  :: BY_IMF,BZ_IMF
      REAL*8,  INTENT(OUT) :: BY_p,BZ_p
      DATA e1/0.06816/
      DATA e2/0.55417/

      pi=3.14159265358979

      IF ((BY_IMF .eq. 0.0) .AND. (BZ_IMF .eq. 0.0)) THEN
        THETA=0.0
      ELSE
        THETA=DATAN2(BY_IMF,BZ_IMF)
        IF (THETA .lt. 0.0 ) THEN
          THETA = THETA + 2*pi
        ENDIF
      ENDIF

      BY_p=(e1+e2*(dcos(THETA/2.0))**2)*BY_IMF
      BZ_p=(e1+e2*(dcos(THETA/2.0))**2)*BZ_IMF

      RETURN
      END SUBROUTINE getIMF_penetration

C***********************************************************************

      SUBROUTINE checkIfInsideMpSATURN(RR, x, NOUT, Dp, a1, a2, dK)
      IMPLICIT REAL*8(A-H,O-Z)
      REAL*8,  INTENT(IN)   :: RR, x
      REAL*8,  INTENT(IN)   :: Dp
      REAL*8,  INTENT(IN)   :: a1, a2, dK
      INTEGER, INTENT(OUT)  :: NOUT
      r0 = a1 * Dp**(-a2)

      cosTheta = x / RR
      Rmp = r0 * (2 / (1 + cosTheta))**dK

      IF (Rmp .lt. RR) THEN
          NOUT = 0
      ELSE
          NOUT = 1
      ENDIF

      RETURN
      END SUBROUTINE checkIfInsideMpSATURN

C***********************************************************************

      SUBROUTINE KRONIAN(NM,R,T,F,BR,BT,BF)

C     BASED ON THE SUBROUTINE IGRF WRITTEN BY N.A. TSYGANENKO (1979)
C     MODIFIED BY KRISHAN KHURANA, JULY, 1996. AND   NOV. 2004.
C
C     CALCULATES COMPONENTS OF MAIN KRONIAN FIELD IN RIGHT HANDED SPHERICAL
C     COORD SYSTEM. BASED ON THE  SPHERICAL HARMONIC COEFFICIENTS GIVEN BY 
C     ACUNA ET AL. [1983] (Z3 MODEL,  V1+V2 DATA) 
C     MAXIMUM ORDER OF HARMONICS TAKEN INTO ACCOUNT (NOT MORE THAN 0RDER 3)
C
C
C     IT IS ASSUMED THAT THE TRAJECTORY IS IS IN RIGHT HANDED S III COORDINATES.
C     THE OUTPUT IS ALSO IN RTP (RH) COORDINATES.
C
C            INPUT:  NM (INTEGER)- MAXIMUM ORDER OF HARMONICS TAKEN
C                                  INTO ACCOUNT (NM.LE.12)
C
C                    R,T,F (REAL)- POSITION OF DESIRED FIELD VALUE IN
C                                  SPHERICAL JOVIGRAPHIC COORDINATE SYSTEM
C                                  (R IN PLANET RADII, COLATITUDE T AND
C                                   RH LONGITUDE F IN RADIANS)
C
C           OUTPUT: BR,BT,BF (REAL)- COMPONENTS OF THE INTERNAL PORTION
C                                    OF THE MAIN MAGNETIC FIELD IN
C                                    SPHERICAL S III COORD SYSTEM
C                                    (VALUES GIVEN IN GAMMA)
C
C***********************************************************************
      IMPLICIT REAL*8 (A-H,O-Z)
      INTEGER, INTENT(IN)  :: NM
      REAL*8,  INTENT(IN)  :: R, T, F
      REAL*8,  INTENT(OUT) :: BR, BT, BF
      LOGICAL BK, BM
      DIMENSION A(13), B(13), G(91), H(91), REC(91)

      DATA G/0.,21144.,24.,1449.,-179.,-35.,2012.,
     ~       53.,92.,-40.,
     ~       81*0.0/

      DATA H/0.,0.,-0.9110,0.,43.5680,-44.7762,0.,
     ~       140.9273,-92.8876,68.2759,
     ~       81*0.0/

      DATA FIRSTI/0.0/
C     WRITE(1,'(5F15.3)')(G(I),I=1,10)
C     WRITE(1,'(5F14.3)')(H(I),I=1,10)
      IF (FIRSTI.EQ.0.0) GO TO 1
      GO TO 6
  1   FIRSTI=1.0
      G(1)=0.
      H(1)=0.
      KNM=15
      DO 2 N=1,13
      N2=2*N-1
      N2=N2*(N2-2)
      DO 2 M=1,N
      MN=N*(N-1)/2+M
  2   REC(MN)=DBLE((N-M)*(N+M-2))/DBLE(N2)
      S=1.
      DO 5 N=2,13
      MN=N*(N-1)/2+1
      S=S*DBLE(2*N-3)/DBLE(N-1)
      G(MN)=G(MN)*S
      H(MN)=H(MN)*S
      P=S
      DO 5 M=2,N
      AA=1.
      IF (M.EQ.2) AA=2.
      P=P*DSQRT(AA*DBLE(N-M+1)/DBLE(N+M-2))
      MNN=MN+M-1
      G(MNN)=G(MNN)*P
  5   H(MNN)=H(MNN)*P
  6   IF(KNM.EQ.NM) GOTO 61
      KNM=NM
      K=KNM+1
 61   PP=1./R
      P=PP
      DO 7 N=1,K
      P=P*PP
      A(N)=P
  7   B(N)=P*N
      P=1.
      D=0.
      BBR=0.
      BBT=0.
      BBF=0.
      U=T
      CF=DCOS(F)
      SF=DSIN(F)
      C=DCOS(U)
      S=DSIN(U)
      BK=(S.LT.1.D-5)
      DO 12 M=1,K
      BM=(M.EQ.1)
      IF(BM) GOTO 8
      MM=M-1
      W=X
      X=W*CF+Y*SF
      Y=Y*CF-W*SF
      GOTO 9
  8   X=0.
      Y=1.
  9   Q=P
      Z=D
      BI=0.
      P2=0.
      D2=0.
      DO 11 N=M,K
      AN=A(N)
      MN=N*(N-1)/2+M
      E=G(MN)
      HH=H(MN)
      W=E*Y+HH*X
      IF (DABS(P2).LT.1.D-38) P2=0.0
      IF (DABS(Q).LT.1.D-38) Q=0.0
      BBR=BBR+B(N)*W*Q
      BBT=BBT-AN*W*Z
      IF(BM) GOTO 10
      QQ=Q
      IF(BK) QQ=Z
      BI=BI+AN*(E*X-HH*Y)*QQ
  10  XK=REC(MN)
      DP=C*Z-S*Q-XK*D2
      PM=C*Q-XK*P2
      D2=Z
      P2=Q
      Z=DP
      Q=PM
   11 CONTINUE
      D=S*D+C*P
      P=S*P
      IF(BM) GOTO 12
      BI=BI*MM
      BBF=BBF+BI
  12  CONTINUE
      BR=BBR
      BT=BBT
      IF(BK) GOTO 13
      BF=BBF/S
      GOTO 14
  13  IF(C.LT.0.) BBF=-BBF
      BF=BBF
  14  CONTINUE
      RETURN
      END SUBROUTINE KRONIAN
  
C***********************************************************************