FUNCTION STARKA(BETA,A,DIV,FAC)
C
C     Approximate expressions for the hydrogen Stark profile
C
C     Input: BETA  - delta lambda in beta units,
C            BETAD - Doppler width in beta units
C            A     - auxiliary parameter
C                    A=1.5*LOG(BETAD)-1.671
C            DIV   - only for A > 1; division point between Doppler
C                    and asymptotic Stark wing, expressed in units
C                    of betad.
C                    DIV = solution of equation
C                    exp(-(beta/betad)**2)/betad/sqrt(pi)=
C                     = 1.5*FAC*beta**-5/2
C                    (ie. the point where Doppler profile is equal to
C                     the asymptotic Holtsmark)
C                    In order to save computer time, the division point
C                    DIV is calculated in advance by routine DIVSTR.
C            FAC   - factor by which the Holtsmark profile is to be
C                    multiplied to get total Stark Profile
C                    FAC should be taken to 2 for hydrogen, (and =1
C                    for He II)
C
      INCLUDE 'PARAMS.FOR'
      PARAMETER (F0=-0.5758228,F1=0.4796232,F2=0.07209481/2.,AL=1.26)
      PARAMETER (SD=0.5641895,SLO=-2.5,TRHA=1.5,BL1=1.52,BL2=8.325)
      PARAMETER (SAC=0.07966/2.)
      XD=BETA/BETAD
C
C     for a > 1 Doppler core + asymptotic Holtzmark wing with division
C               point DIV
C
      IF(A.GT.AL) THEN
         IF(XD.LE.DIV) THEN
            STARKA=SD*EXP(-XD*XD)/BETAD
          ELSE
            STARKA=TRHA*FAC*EXP(SLO*LOG(BETA))
         END IF
       ELSE
C
C     empirical formula for a < 1
C
         IF(BETA.LE.BL1) THEN
            STARKA=SAC*FAC
          ELSE IF(BETA.LT.BL2) THEN
            XL=LOG(BETA)
            FL=(F0*XL+F1)*XL
            STARKA=F2*FAC*EXP(FL)
          ELSE
            STARKA=TRHA*FAC*EXP(SLO*LOG(BETA))
         END IF
      END IF
      RETURN
      END