SUBROUTINE ELDENS(ID,T,AN,ANE)
C     ==============================
C
C     Evaluation of the electron density and the total hydrogen
C     number density for a given total particle number density
C     and temperature;
C     by solving the set of Saha equations, charge conservation and
C     particle conservation equations (by a Newton-Raphson method)
C
C     Input parameters:
C     T    - temperature
C     AN   - total particle number density
C
C     Output:
C     ANE   - electron density
C     ANP   - proton number density
C     AHTOT - total hydrogen number density
C     AHMOL - relativer number of hydrogen molecules with respect to the
C             total number of hydrogens
C     ENERG - part of the internal energy: excitation and ionization
C
      INCLUDE 'PARAMS.FOR'
      INCLUDE 'MODELP.FOR'
      common/hydmol/anhmi,ahmol
      common/hydato/ah,anh,anp
      common/nerela/anerel
      parameter (un=1.d0,two=2.d0,half=0.5d0)
      DIMENSION R(3,3),S(3),P(3)
C
      TK=BOLK*T
      if(ifmol.gt.0.and.t.lt.tmolim) then
         aein=an*anerel
         call moleq(id,t,an,aein,ane,0)
         anerel=ane/an
         return
      end if
c
      QM=0.
      Q2=0.
      QP=0.
      Q=0.
      DQN=0.
      TK=BOLK*T
      THET=5.0404D3/T
C
C     Coefficients entering ionization (dissociation) balance of:
C     atomic hydrogen          - QH;
C     negative hydrogen ion    - QM
C     hydrogen molecule        - Q2
C     ion of hydrogen molecule - QP
C
      IF(IATREF.EQ.IATH) THEN
         QM=1.0353D-16/T/SQRT(T)*EXP(8762.9/T)
         QH0=EXP((15.38287+1.5*LOG10(T)-13.595*THET)*2.30258509299405)
c
         if(t.gt.16000.) then
            ih2=0
          else
            ih2=1
            QP=TK*EXP((-11.206998+THET*(2.7942767+THET*
     *      (0.079196803-0.024790744*THET)))*2.30258509299405)
            Q2=TK*EXP((-12.533505+THET*(4.9251644+THET*
     *      (-0.056191273+0.0032687661*THET)))*2.30258509299405)
         end if
      END IF
C
C     Initial estimate of the electron density
C
      if(anerel.le.0.) then
         if(t.gt.1.e4) then
            anerel=0.5
          else
            if(elec(id).gt.0..and.dens(id).gt.0.) then
               anerel=elec(id)/(elec(id)+dens(id)/wmm(id))
             else
               anerel=0.1
            end if
         end if
      end if
c
      ANE=AN*ANEREL
      IT=0
C
C     Basic Newton-Raphson loop - solution of the non-linear set
C     for the unknown vector P, consistiong of AH, ANH (neutral
C     hydrogen number density) and ANE.
C
   10 IT=IT+1
C
C     procedure STATE determines Q (and DQN) - the total charge (and its
C     derivative wrt temperature) due to ionization of all atoms which
C     are considered (both explicit and non-explicit), by solving the set
C     of Saha equations for the current values of T and ANE
C
      CALL STATE(ID,T,ANE,Q)
      QH=QH0*2./PFSTD(1,1)
C
C     Auxiliary parameters for evaluating the elements of matrix of
C     linearized equations.
C     Note that complexity of the matrix depends on whether the hydrogen
C     molecule is taken into account
C     Treatment of hydrogen ionization-dissociation is based on
C     Mihalas, in Methods in Comput. Phys. 7, p.10 (1967)
C
      IF(IATREF.EQ.IATH) THEN
      G2=QH/ANE
      G3=0.
      G4=0.
      G5=0.
      D=0.
      E=0.
      G3=QM*ANE
      A=UN+G2+G3
      D=G2-G3
      IF(IT.LE.1) THEN
         IF(IH2.EQ.0) THEN
            F1=UN/A
            FE=D/A+Q
          ELSE
            E=G2*QP/Q2
            B=TWO*(UN+E)
            GG=ANE*Q2
            C1=B*(GG*B+A*D)-E*A*A
            C2=A*(TWO*E+B*Q)-D*B
            C3=-E-B*Q
            F1=(SQRT(C2*C2-4.*C1*C3)-C2)*HALF/C1
            FE=F1*D+E*(UN-A*F1)/B+Q
         END IF
         AH=ANE/FE
         ANH=AH*F1
      END IF
      AE=ANH/ANE
      GG=AE*QP
      E=ANH*Q2
      B=ANH*QM
C
C     Matrix of the linearized system R, and the rhs vector S
C
      R(1,1)=YTOT(ID)
c     R(1,2)=0.
      r(1,2)=-two*(anh*q2+gg)
      R(1,3)=UN
      R(2,1)=-Q
      R(2,2)=-D-TWO*GG
      R(2,3)=UN+B+AE*(G2+GG)-DQN*AH
      R(3,1)=-UN
      R(3,2)=A+4.*(anh*q2+GG)
      R(3,3)=B-AE*(G2+TWO*GG)
      S(1)=AN-ANE-YTOT(ID)*AH+anh*(anh*q2+gg)
      S(2)=ANH*(D+GG)+Q*AH-ANE
      S(3)=AH-ANH*(A+TWO*(anh*q2+GG))
C
C     Solution of the linearized equations for the correction vector P
C
      CALL LINEQS(R,S,P,3,3)
C
C     New values of AH, ANH, and ANE
C
      AH=AH+P(1)
      ANH=ANH+P(2)
      DELNE=P(3)
      ANE=ANE+DELNE
C
C     hydrogen is not the reference atom
C
      ELSE
C
C     Matrix of the linearized system R, and the rhs vector S
C
      IF(IT.EQ.1) THEN
         ANE=AN*HALF
         AH=ANE/YTOT(ID)
      END IF
      R(1,1)=YTOT(ID)
      R(1,2)=UN
      R(2,1)=-Q-QREF
      R(2,2)=UN-(DQN+DQNR)*AH
      S(1)=AN-ANE-YTOT(ID)*AH
      S(2)=(Q+QREF)*AH-ANE
C
C     Solution of the linearized equations for the correction vector P
C
      CALL LINEQS(R,S,P,2,3)
      AH=AH+P(1)
      DELNE=P(2)
      ANE=ANE+DELNE
      END IF
C
C     Convergence criterion
C
      IF(ANE.LE.0.) ANE=1.D-7*AN
      IF(ABS(DELNE/ANE).GT.1.D-6.AND.IT.LE.20) GO TO 10
C
C     ANEREL is the exact ratio betwen electron density and total
C     particle density, which is going to be used in the subseguent
C     call of ELDENS
C
      ANEREL=ANE/AN
      AHTOT=AH
      IF(IATREF.EQ.IATH) THEN
c     AHMOL=TWO*ANH*(ANH*Q2+ANH/ANE*QP)/AH
      AHMOL=ANH*ANH*Q2
      ANP=ANH/ANE*QH
      ANHMI=ANH*ANE*QM
      anhn=anh+anp+anhmi+2.*ahmol
      wmm(id)=wmy(id)/(ytot(id)-ahmol/anhn)*hmass
      END IF
C
      RETURN
      END