************************************************************************
      SUBROUTINE REF(PATH)
************************************************************************
      IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N)

      PARAMETER (NSPEC=4)
      real*8 path(8)
      COMMON /QCONST/PAI,RAD,DEG,RE,CVLCTY
      COMMON /QPATH /FKC,DEL,EPS,GMU,ALF,RAY
      COMMON /QPARAM/RREF,RBTM,RREFP,
     &               PF0(NSPEC),GF0(NSPEC),
     &               XN0,THERM,ETA(NSPEC),SHINV(NSPEC),SHINVP,
     &               RRLI,HWLI,
     &               XLPP,HWPP,ERPP,DRR1,DTH1,DPH1,DRR2,DTH2,DPH2,
     &               RRMIN,RRMAX,ERRR,ERRA,HMIN,DISOUT,
     &               NPATH,NLPMAX,INDOUT,INDMSH,IEOF,
     &               IYEAR,IMNTH,IDATE
      COMMON /QCOND /L,LL,NN,NLOOP,a
      COMMON /QREF  /PMU,PMUS,PMU0,
     &               PSI   ,
     &               COPSI ,SIPSI ,
     &               COPSIS,SIPSIS,
     &               PX(NSPEC),PY(NSPEC),PZ(NSPEC),
     &               PA1,PA2,PA3,
     &               PK1,PK2,PK4,PK5,
     &               PA ,PB ,PC ,PD 
      COMMON /QIGRF /YR    ,YT    ,YP    ,
     &               DYRDRR,DYTDRR,DYPDRR,
     &               DYRDTH,DYTDTH,DYPDTH,
     &               DYRDPH,DYTDPH,DYPDPH,
     &               FH    ,
     &               DBBDRR,DBBDTH,DBBDPH
      COMMON /QDENS /XNS   (NSPEC),
     &               DNSDRR(NSPEC),
     &               DNSDTH(NSPEC),
     &               DNSDPH(NSPEC)

      CALL IGRF(PATH)
c      CALL DIPOLE(PATH)
      CALL DENS(PATH)
      
      DO 10 I=1,NSPEC      
      PX(I) = PF0(I)*XNS(I)/FKC**2
      PY(I) = GF0(I)*FH  /FKC
      PZ(I) = PY(I)**2-1.0D0
   10 CONTINUE
          
      PMU0   = dSQRT(PATH(4)**2+PATH(5)**2+PATH(6)**2)
      COPSI  = (YR*PATH(4)
     &          +YT*PATH(5)
     &          +YP*PATH(6))/PMU0
      COPSIS = COPSI**2 
      SIPSIS = dABS(1.0D0-COPSIS)
      SIPSI  = dSQRT(SIPSIS)
      if (dabs(copsi).gt.1.0) copsi=1.0d0
      PSI    = dACOS(COPSI)*DEG

      PA1 = 0.0D0
      PA2 = 0.0D0
      PA3 = 0.0D0
 
      DO 20 I=1,NSPEC     
      PA1 = PA1+PX(I)
      PA2 = PA2+PX(I)      /PZ(I)
      PA3 = PA3+PX(I)*PY(I)/PZ(I)
   20 CONTINUE

      PK1 = 1.0D0-PA1
      PK2 = 1.0D0+PA2
      PK4 = PK1*PK2
      PK5 = PK2**2-PA3**2
      PA  = PK1*COPSIS+PK2*SIPSIS
      PB  =-PK4*(COPSIS+1.0D0)-PK5*SIPSIS
      PC  = PK1*PK5
      PD  = DABS(PB**2-4.0D0*PA*PC)
      if (pd.lt.0.0d0) go to 1040
      PMUS= (-PB-SQRT(PD))/2.0/PA
c       Above is Extraordinary mode(whistler mode)
      if (pmus.lt.0.0d0) go to 1050
      PMU = dSQRT(PMUS)
 1040 CONTINUE
      IF (PATH(1).LT.6370.0d0+90.0d0) THEN
      PMU=1.0d0
      gmu=1.0d0
      END IF
 1050 CONTINUE

      RETURN
      END