FUNCTION STARKA(BETA,FAC)
C     =========================
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.761
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   - Multiplicative factor (2. for H I; 1. for He II)
C
      INCLUDE 'IMPLIC.FOR'
      INCLUDE 'BASICS.FOR'
      INCLUDE 'MODELQ.FOR'
      PARAMETER (F0=-0.5758228,F1=0.4796232,F2=0.07209481,AL=1.26)
      PARAMETER (SD=0.5641895,SLO=-2.5,THRA=1.5,BL1=1.14,BL2=11.4)
      PARAMETER (SAC=0.08, PISQ1=UN/1.77245385090551D0)
C
C     for a > 1 Doppler core + asymptotic Holtzmark wing with division
C               point DIV
C
      BETAD1=UN/BETAD
      IF(ADH.GT.AL) THEN
         XD=BETA*BETAD1
         IF(XD.LE.DIVH) THEN
            STARKA=EXP(-XD*XD)*BETAD1*PISQ1
          ELSE
            STARKA=THRA*FAC*EXP(SLO*LOG(BETA))
         END IF
       ELSE
C
C     empirical formula for a < 1
C
         IF(BETA.LE.BL1) THEN
            STARKA=SAC
          ELSE IF(BETA.LT.BL2) THEN
            XL=LOG(BETA)
            FL=(F0*XL+F1)*XL
            STARKA=F2*EXP(FL)
          ELSE
            STARKA=THRA*FAC*EXP(SLO*LOG(BETA))
         END IF
      END IF
      RETURN
      END