FUNCTION PHE1(ID,FREQ,ILINE)
C     ============================
C
C     Absorption profile for four lines of He I, given by
C     Barnard, Cooper, Smith (1974) JQSRT 14, 1025 for the 4471 line;
C     Shamey (1969) PhD thesis, for other lines
C
C     Input: ID    - depth index
C            FREQ  - frequency
C            ILINE - index of the line ( = 1  for 4471,
C                                        = 2  for 4387,
C                                        = 3  for 4026,
C                                        = 4  for 4922)
C
C     Output: PHE1 - profile coefficient in frequency units,
C                    normalized to sqrt(pi) [not unity]
C
      INCLUDE 'PARAMS.FOR'
      INCLUDE 'MODELP.FOR'
      PARAMETER (NT=4)
      COMMON/PROHE1/PRFHE1(50,NT,8,3),DLMHE1(50,8,3),XNEHE1(8),
     *              NWLAM(8,NT)
      COMMON/PRO447/PRF447(80,NT,7),DLM447(80,7),XNE447(7)
      DIMENSION WLAM0(4),XT0(NT),XX(3),WX(3),YY(2),PP(2),ZZ(3),WZ(3)
      DATA WLAM0 / 4471.50, 4387.93, 4026.20, 4921.93/
      DATA XT0/ 3.699, 4.000, 4.301, 4.602/
C
C     temperature is modified in order to account for the
C     effect of turbulent velocity on the Doppler width
C
      T=TEMP(ID)+2.42E-8*VTURB(ID)
      TL=LOG10(T)
      ANE=ELEC(ID)
      ANEL=LOG10(ANE)
      ALAM=2.997925E18/FREQ
      DLAM=ALAM-WLAM0(ILINE)
      DOPL=SQRT(4.125E7*T)*WLAM0(ILINE)/2.997925E10
C
      IF(TL.GT.XT0(NT)+0.1) GO TO 5
      IF(ILINE.EQ.1.AND.ANEL.GE.XNE447(1)) GO TO 10
      IF(ILINE.NE.1.AND.ANEL.GE.XNEHE1(1)) GO TO 10
C
C     isolated line approximation for low electron densities
C
    5 A=WTOT(T,ANE,ID,ILINE)/DOPL
      V=ABS(DLAM)/DOPL
      V1=ABS(ALAM-4471.682)/DOPL
      PHE1=VOIGTK(A,V)
      IF(ILINE.EQ.1) PHE1=(8.*PHE1+VOIGTK(A,V1))/9.
      RETURN
C
C     otherwise, interpolation (or extrapolation) in tables
C
   10 NX=3
      NZ=3
      NY=2
      NE=8
      ILNE=ILINE-1
      IF(ILINE.EQ.1) NE=7
C
C     Interpolation in electron density
C
      DO JZ=1,NE-1
         IPZ=JZ
         IF(ILINE.EQ.1.AND.ANEL.LE.XNE447(JZ+1)) GO TO 30
         IF(ILINE.NE.1.AND.ANEL.LE.XNEHE1(JZ+1)) GO TO 30
      END DO
   30 N0Z=IPZ-NZ/2+1
      IF(N0Z.LT.1) N0Z=1
      IF(N0Z.GT.NE-NZ+1) N0Z=NE-NZ+1
      N1Z=N0Z+NZ-1
      DO 300 JZ=N0Z,N1Z
         I0Z=JZ-N0Z+1
         IF(ILINE.EQ.1) ZZ(I0Z)=XNE447(JZ)
         IF(ILINE.NE.1) ZZ(I0Z)=XNEHE1(JZ)
C
C        Interpolation in temperature
C
         DO IX=1,NT-1
            IPX=IX
            IF(TL.LE.XT0(IX+1)) GO TO 50
         END DO
   50    N0X=IPX-NX/2+1
         IF(N0X.LT.1) N0X=1
         IF(N0X.GT.NT-NX+1) N0X=NT-NX+1
         N1X=N0X+NX-1
         DO 200 IX=N0X,N1X
            I0X=IX-N0X+1
            XX(I0X)=XT0(IX)
C
C           Interpolation in wavelength
C
C           1. For delta lambda beyond tabulated values - special
C              extrapolation (Cooper's suggestion)
C
            NLST=NWLAM(JZ,ILINE)
            IF(ILINE.EQ.1) THEN
               D1=DLM447(1,JZ)
               D2=DLM447(NLST,JZ)
               IF(DLAM.LT.D1) THEN
                  PRF0=EXTPRF(DLAM,IX,ILINE,ZZ(I0Z),D1,PRF447(1,IX,JZ))
                  GO TO 150
                ELSE IF(DLAM.GT.D2) THEN
                  PRF0=EXTPRF(DLAM,IX,ILINE,ZZ(I0Z),D2,
     *                PRF447(NLST,IX,JZ))
                  GO TO 150
               END IF
             ELSE
               D1=DLMHE1(1,JZ,ILNE)
               D2=DLMHE1(NLST,JZ,ILNE)
               IF(DLAM.LT.D1) THEN
                  PRF0=EXTPRF(DLAM,IX,ILINE,ZZ(I0Z),D1,
     *                PRFHE1(1,IX,JZ,ILNE))
                  GO TO 150
                ELSE IF(DLAM.GT.D2) THEN
                  PRF0=EXTPRF(DLAM,IX,ILINE,ZZ(I0Z),D2,
     *                PRFHE1(NLST,IX,JZ,ILNE))
                  GO TO 150
               END IF
            END IF
C
C           normal linear interpolation in wavelength
C           (for 4471, linear interpolation in logarithms)
C
            DO IY=1,NLST-1
               IPY=IY
               IF(ILINE.EQ.1.AND.DLAM.LE.DLM447(IY+1,JZ)) GO TO 70
               IF(ILINE.NE.1.AND.DLAM.LE.DLMHE1(IY+1,JZ,ILNE))
     *            GO TO 70
            END DO
   70       N0Y=IPY-NY/2+1
            IF(N0Y.LT.1) N0Y=1
            IF(N0Y.GT.NLST-NY+1) N0Y=NLST-NY+1
            N1Y=N0Y+NY-1
            DO IY=N0Y,N1Y
               I0=IY-N0Y+1
               IF(ILINE.EQ.1) YY(I0)=DLM447(IY,JZ)
               IF(ILINE.EQ.1) PP(I0)=LOG(PRF447(IY,IX,JZ))
               IF(ILINE.NE.1) YY(I0)=DLMHE1(IY,JZ,ILNE)
               IF(ILINE.NE.1) PP(I0)=PRFHE1(IY,IX,JZ,ILNE)
           END DO
           IF(ILINE.NE.1) THEN
              WX(I0X)=(PP(2)*(DLAM-YY(1))+PP(1)*(YY(2)-DLAM))/
     *                (YY(2)-YY(1))
            ELSE
             WX(I0X)=(PP(2)*(DLAM-YY(1))+PP(1)*(YY(2)-DLAM))/
     *                (YY(2)-YY(1))
             WX(I0X)=EXP(WX(I0X))
           END IF
           GO TO 200
  150      WX(I0X)=PRF0
  200   CONTINUE
        WZ(I0Z)=YINT(XX,WX,TL)
  300 CONTINUE
      W0=YINT(ZZ,WZ,ANEL)
      PHE1=W0*DOPL*1.772454
      RETURN
      END