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

      PARAMETER (NSPEC=4)
      real*8 path(8),dpath(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,
     &               BB    ,
     &               DBBDRR,DBBDTH,DBBDPH
      COMMON /QDENS /XNS   (NSPEC),
     &               DNSDRR(NSPEC),
     &               DNSDTH(NSPEC),
     &               DNSDPH(NSPEC)

      real*8 DMUDRR,DMUDTH,DMUDPH,
     &       DK2DF ,DA3DF ,DENOM 

      CALL REF(PATH)

      IF(IMDFY.EQ.1) THEN

      PATH(4) = PATH(4)*PMU/PMU0
      PATH(5) = PATH(5)*PMU/PMU0
      PATH(6) = PATH(6)*PMU/PMU0

      DEL = dATAN2(PATH(5),PATH(4))
      EPS = dATAN (PATH(6)/dSQRT(PATH(4)**2+PATH(5)**2))
      DEL = DEL*DEG
      EPS = EPS*DEG

      END IF

      DMUDRR = 0.0D0
      DMUDTH = 0.0D0
      DMUDPH = 0.0D0
      DK2DF  = 0.0D0
      DA3DF  = 0.0D0
      DENOM  = (4.0D0*PA*PMUS+2.0D0*PB)*PMU

      DO 20 I=1,NSPEC
      DADX   = -                       COPSIS
      DBDX   = -(PK1/PZ(I)-PK2)*(1.0D0+COPSIS)
      DCDX   = - PK5
      PK2A3Y = 2.0D0*(PK2-PA3*PY(I))
      DADX   = DADX+           SIPSIS/PZ(I)
      DBDX   = DBDX-    PK2A3Y*SIPSIS/PZ(I)
      DCDX   = DCDX+PK1*PK2A3Y       /PZ(I)
      DMUDX  = -(DADX*PMUS**2+DBDX*PMUS+DCDX)/DENOM

      DK2DY  = -2.0D0    *PX(I)* PY(I)       /PZ(I)**2
      DK5DY  =  2.0D0*PA3*PX(I)*(PZ(I)+2.0D0)/PZ(I)**2
      DK5DY  = DK5DY+2.0D0*DK2DY*PK2
      DADY   =                           DK2DY*SIPSIS
      DBDY   = -PK1*DK2DY*(1.0D0+COPSIS)-DK5DY*SIPSIS
      DCDY   =  PK1*DK5DY
      DMUDY  = -(DADY*PMUS**2+DBDY*PMUS+DCDY)/DENOM

      DXDRR  = PX(I)*DNSDRR(I)
      DXDTH  = PX(I)*DNSDTH(I)
      DXDPH  = PX(I)*DNSDPH(I)
      DYDRR  = PY(I)*DBBDRR
      DYDTH  = PY(I)*DBBDTH
      DYDPH  = PY(I)*DBBDPH
      DMUDRR = DMUDRR+DMUDX*DXDRR+DMUDY*DYDRR
      DMUDTH = DMUDTH+DMUDX*DXDTH+DMUDY*DYDTH
      DMUDPH = DMUDPH+DMUDX*DXDPH+DMUDY*DYDPH

      DK2DFI =  2.0D0       *PX(I)      /PZ(I)**2
      DA3DFI = (2.0D0-PZ(I))*PX(I)*PY(I)/PZ(I)**2
      DK2DF  = DK2DF+DK2DFI/FKC
      DA3DF  = DA3DF+DA3DFI/FKC
   20 CONTINUE

      DADPS  = PA1+PA2
      DBDPS  = PK4-PK5
      DMUDPS = (DADPS*PMUS+DBDPS)*PMUS
      DMUDPS = -2.0D0*COPSI*DMUDPS/DENOM

      DPSDRR = DYRDRR*PATH(4)
     &         +DYTDRR*PATH(5)
     &         +DYPDRR*PATH(6)
      DPSDTH = DYRDTH*PATH(4)
     &         +DYTDTH*PATH(5)
     &         +DYPDTH*PATH(6)
      DPSDPH = DYRDPH*PATH(4)
     &         +DYTDPH*PATH(5)
     &         +DYPDPH*PATH(6)
      DPSDRR = -DPSDRR/PMU
      DPSDTH = -DPSDTH/PMU
      DPSDPH = -DPSDPH/PMU
      
      if (XNS(1).NE.0.0d0) then
         DMUDRR = (DMUDRR+DMUDPS*DPSDRR)/PMU
         DMUDTH = (DMUDTH+DMUDPS*DPSDTH)/PMU
         DMUDPH = (DMUDPH+DMUDPS*DPSDPH)/PMU
         DMUDVR = DMUDPS*(PATH(4)*COPSI/PMU-YR)
         DMUDVT = DMUDPS*(PATH(5)*COPSI/PMU-YT)
         DMUDVP = DMUDPS*(PATH(6)*COPSI/PMU-YP)
      else
         DMUDRR = 0.0d0
         DMUDTH = 0.0d0
         DMUDPH = 0.0d0
         DMUDVR = 0.0d0
         DMUDVT = 0.0d0
         DMUDVP = 0.0d0
      endif

      DK1DF     = 2.0D0* PA1/FKC
      DK4DF     =        PK1*DK2DF+PK2*DK1DF
      DK5DF     = 2.0D0*(PK2*DK2DF-PA3*DA3DF)
      DCDF      =        PK1*DK5DF    +PK5*DK1DF
      DADF      =  DK1DF*       COPSIS +DK2DF*SIPSIS
      DBDF      = -DK4DF*(1.0D0+COPSIS)-DK5DF    *SIPSIS
      DMUDF     = -((DADF*PMUS+DBDF)*PMUS+DCDF)/DENOM
      if (XNS(1).EQ.0.0d0) dmudf=0.0d0
      GMU       = PMU+FKC*DMUDF

      IF (PATH(1).LT.6460.0) THEN
       DMUDVR = 0.0d0
       DMUDVR = 0.0d0
       DMUDVP = 0.0d0
       DMUDRR = 0.0d0 
       DMUDTH = 0.0d0 
       DMUDPH = 0.0d0 
      END IF 

      DPATH(1) = PATH(4)-DMUDVR
      DPATH(2) = PATH(5)-DMUDVT
      DPATH(3) = PATH(6)-DMUDVP
      DPATH(1) = DPATH(1)/PMUS
      DPATH(2) = DPATH(2)/PMUS/PATH(1)
      DPATH(3) = DPATH(3)/PMUS/PATH(1)/dSIN(PATH(2)*RAD)
      DPATH(4) = DMUDRR
     &           +PATH(5)*DPATH(2)
     &           +PATH(6)*DPATH(3)*dSIN  (PATH(2)*RAD)
      DPATH(5) = DMUDTH/PATH(1)
     &           -PATH(5)*DPATH(1)/PATH(1)
     &           +PATH(6)*DPATH(3)*dCOS  (PATH(2)*RAD)
      DPATH(6) = DMUDPH/PATH(1)/dSIN(PATH(2)*RAD)
     &           -PATH(6)*DPATH(1)/PATH(1)
     &           -PATH(6)*DPATH(2)/dTAN(PATH(2)*RAD)
      DPATH(2) = DPATH(2)*DEG
      DPATH(3) = DPATH(3)*DEG
      DPATH(7) = 1.0D0/PMU/CVLCTY
      DPATH(8) = 1.0D0/CVLCTY
     
      if (XNS(1).eq.0.0d0) alf=0.0d0
      ALF = -ATAN(DMUDPS*SIPSI/PMU)*DEG
      RAY = PSI+ALF

      RETURN
      END