SUBROUTINE SETRAY
C     =================
C
C     Setup impact rays and angles
C      (assumes one impact ray tangent to every depth layer)
C
      INCLUDE 'PARAMS.FOR'
      INCLUDE 'MODELP.FOR'
      INCLUDE 'WINCOM.FOR'
      PARAMETER (PI4=4.*3.141592654)
      PARAMETER (UN=1., TWO=2., HALF=0.5)
      DIMENSION RS(MDEPF ),RDX(MDEPF )
      DIMENSION ZIU(MDEPTH),VIU(MDEPTH),ZIUF(MDEPF ),VIUF(MDEPF )
C
C     Fine radial grid
C
      if(ndf.eq.0.or.ndf.eq.nd) then
         ndf=nd
         DO ID=1,NDF
            DENSF(ID)=DENS(ID)
         END DO
      else
         XR1=LOG(DENS(1))
         XR2=LOG(DENS(ND))
         DXR=(XR2-XR1)/FLOAT(NDF-1)
         DO ID=1,NDF
            DENSF(ID)=EXP(XR1+FLOAT(ID-1)*DXR)
         END DO
      end if
C
C
C     Impact rays
C
      NREXT=ND
      DO ID=1,NREXT
        PIM(ID)=RD(ID)
        NUD(ID)=ID
      END DO
      DO IU=1,NRCORE
        PIM(NREXT+IU)=FLOAT(NRCORE-IU)/FLOAT(NRCORE)*RCORE
        NUD(NREXT+IU)=ND
      END DO
      KMU=NREXT+NRCORE
C
C     Angles
C
      DO ID=1,ND
        RD1=UN/RD(ID)
        DO IU=ID,KMU
          PRR=PIM(IU)*RD1
          BMU(IU,ID)=SQRT(UN-PRR*PRR)
        END DO
      END DO
C
C     Depth increments along each ray
C
      DELZ(1,1)=0.
      DFRQ(1,1)=0.
      DO IU=2,KMU
        NUDF(IU)=NUD(IU)
        IU1=IU
        IF(IU.GT.ND) IU1=ND
        DO ID=1,IU1-1
          DELZ(IU,ID)=BMU(IU,ID)*RD(ID)-BMU(IU,ID+1)*RD(ID+1)
          DFRQ(IU,ID)=BMU(IU,ID)*VEL(ID)/CL
          JD=2*NUD(IU)-ID
          DFRQ(IU,JD)=-DFRQ(IU,ID)
        END DO
        DELZ(IU,IU1)=DELZ(IU,IU1-1)
        DFRQ(IU,IU1)=0.
        IF(IU.GT.NREXT) DFRQ(IU,ND)=BMU(IU,ND)*VEL(ND)/CL
      END DO
C
C Finer grid along the NFIRY most external rays
C   velocity steps DVD(ID)
C
      XMD4=XMDOT/PI4
      CLV=UN/CL
      DO ID=1,ND
        DVD(ID)=SQRT(1.6D7*TEMP(ID)+VTURB(ID)) * 0.3
c        DVD(ID)=SQRT(1.6D7*TEMP(ID))
      END DO
      NUDX=ND
      DO IU=2,NFIRY
        IF(PIM(IU).GT.0.) THEN
          DO ID=1,NUD(IU)
            IID=NUD(IU)-ID+1
            ZIU(ID)=VEL(IID)
            VIU(ID)=DFRQ(IU,IID)*CL
          ENDDO
        ELSE
          DO ID=1,NUD(IU)
            IID=NUD(IU)-ID+1
            ZIU(ID)=RD(IID)
            VIU(ID)=DFRQ(IU,IID)*CL
          ENDDO
        ENDIF
        NUDF(IU)=1
        VIUF(1)=DFRQ(IU,1)*CL
        DO ID=1,NUD(IU)-1
          VZ1=DFRQ(IU,ID)*CL
          VZ2=DFRQ(IU,ID+1)*CL
          NFG=int((VZ1-VZ2)/DVD(ID))+1
          XFG=(VZ1-VZ2)/DFLOAT(NFG)
          IV0=NUDF(IU)
          DO IV=1,NFG
            VIUF(IV0+IV)=VZ1-DFLOAT(IV)*XFG
          ENDDO
          NUDF(IU)=NUDF(IU)+NFG
          IF(NUDF(IU).GT.MDEPF )
     +      CALL quit('Too many points in fine grid - SETRAY')
        END DO
        IF(NUDF(IU).GT.NUDX) NUDX=NUDF(IU)
        INRP=2
        IF(IU.GT.8) INRP=4
        CALL INTERP(VIU,ZIU,VIUF,ZIUF,NUD(IU),NUDF(IU),INRP,0,0)
        IF(PIM(IU).GT.0.) THEN
          DO ID=1,NUDF(IU)
            DMU=VIUF(ID)/ZIUF(ID)
            RS(ID)=PIM(IU)/SQRT(UN-DMU*DMU)
            DFRQF(IU,ID)=VIUF(ID)*CLV
            VELF(IU,ID)=ZIUF(ID)
            RDX(ID)=XMD4/(RS(ID)*RS(ID)*VELF(IU,ID))
            ZIUF(ID)=DMU*RS(ID)
          END DO
        ELSE
          DO ID=1,NUDF(IU)
            RS(ID)=ZIUF(ID)
            DFRQF(IU,ID)=VIUF(ID)*CLV
            VELF(IU,ID)=VIUF(ID)
            RDX(ID)=XMD4/(RS(ID)*RS(ID)*VELF(IU,ID))
          END DO
        END IF
        IF(IU.LE.NREXT) THEN
          DO ID=1,NUDF(IU)
            JD=2*NUDF(IU)-ID
            DFRQF(IU,JD)=-DFRQF(IU,ID)
          END DO
        END IF
        DO ID=1,NUDF(IU)-1
          DELZF(IU,ID)=ZIUF(ID)-ZIUF(ID+1)
        END DO
        DELZF(IU,NUDF(IU))=DELZF(IU,NUDF(IU)-1)
C
C   Assign depth index
C
        KRAY(IU,1)=2
        DRAY(IU,1)=0.
        IDK=1
        DO ID=2,NUDF(IU)
          DO WHILE (RDX(ID).GE.DENSF(IDK).and.idk.le.ndf)
            IDK=IDK+1
          END DO
c          IDK=IDK+1
          IF(IDK.GT.NDF) IDK=NDF
          KRAY(IU,ID)=IDK
          DRAY(IU,ID)=(RDX(ID)-DENSF(IDK-1))/(DENSF(IDK)-DENSF(IDK-1))
        END DO
        IF(IU.LE.NREXT) THEN
          DO ID=1,NUDF(IU)
            JD=2*NUDF(IU)-ID
            KRAY(IU,JD)=KRAY(IU,ID)
            DRAY(IU,JD)=DRAY(IU,ID)
          END DO
        END IF
      END DO
C
C    remaining rays (without finer grid)
C      
      IF(NFIRY.LT.KMU) THEN
        IU=KMU
        KRAY(IU,1)=2
        DRAY(IU,1)=0.
        IDK=1
        DO ID=2,NUDF(IU)
          DO WHILE (DENS(ID).GE.DENSF(IDK).and.idk.le.ndf)
            IDK=IDK+1
          END DO
c          IDK=IDK+1
          IF(IDK.GT.NDF) IDK=NDF
          KRAY(IU,ID)=IDK
          DRAY(IU,ID)=(DENS(ID)-DENSF(IDK-1))/(DENSF(IDK)-DENSF(IDK-1))
        END DO
        DO IU=NFIRY+1,KMU
          DO ID=1,NUDF(IU)
            KRAY(IU,ID)=KRAY(KMU,ID)
            DRAY(IU,ID)=DRAY(KMU,ID)
            DFRQF(IU,ID)=DFRQ(IU,ID)
            DELZF(IU,ID)=DELZ(IU,ID)
          ENDDO
          IF(IU.LE.NREXT) THEN
            DO ID=1,NUDF(IU)
              JD=2*NUDF(IU)-ID
              KRAY(IU,JD)=KRAY(IU,ID)
              DRAY(IU,JD)=DRAY(IU,ID)
              DFRQF(IU,JD)=-DFRQF(IU,ID)
            END DO
          END IF
        END DO
      END IF
C
      NFTOT=0
      DO IU=2,KMU
        IUD=NUDF(IU)
        IF(IU.LE.NREXT) IUD=2*NUDF(IU)-1
        NFTOT=NFTOT+IUD
      ENDDO
      write(10,*) 'NFTOT=',NFTOT
C
      RETURN
      END