function wn(xn,a,id,z)
C     ======================
c
c     evaluation of the occupation probablities for a hydrogenic ion
c     using eqs (4.26), and (4.39) of Hummer,Mihalas Ap.J. 331, 794, 1988.
c     approximate evaluation of Q(beta) - Hummer
c
c     Input: xn  - real number corresponding to quantum number n
c            a   - correlation parameter
c            id  - depth index
c            z   - ionic charge
c
      INCLUDE 'PARAMS.FOR'
      INCLUDE 'MODELP.FOR'
      parameter (p1=0.1402,p2=0.1285,p3=1.,p4=3.15,p5=4.,un=1.)
      parameter (tkn=3.01,ckn=5.33333333,cb=8.59e14)
      parameter (f23=-2./3.)
      parameter (a0=0.529177e-8,wa0=-3.1415926538/6.*a0*a0*a0)
c
c     evaluation of k(n)
c
      if(xn.le.tkn) then
         xkn=un
       else
         xkn=ckn*xn/(xn+un)/(xn+un)
      end if
c
c     evaluation of beta
c
c     beta=cb*bergfc*z*z*z*xkn/(xn*xn*xn*xn)*exp(f23*log(elec(id)))
      beta=cb*z*z*z*xkn/(xn*xn*xn*xn)*exp(f23*log(elec(id)))
c
c     approximate expression for Q(beta)
c
      x=exp(p4*log(un+p3*a))
c     c1=p1*(x+p5*z*a*a*a)    ! previous expression -ERROR !!!!!!
      c1=p1*(x+p5*(z-un)*a*a*a)
      c2=p2*x
      f=(c1*beta*beta*beta)/(un+c2*beta*sqrt(beta))
      wp=f/(un+f)
c
c     contribution from neutral particles
c
      xn2=xn*xn+un
      xnh=0.
      xnhe1=0.
      if(ielh.gt.0) xnh=popul(nfirst(ielh),id)
      if(ielhe1.gt.0) xnhe1=popul(nfirst(ielhe1),id)
      w0=exp(wa0*xn2*xn2*xn2*(xnh+xnhe1))
      W0=1.
      wn=wp*w0
      return
      end