subroutine rhonen(id,t,rho,an,ane)
c     ==================================
c
c     iterative determination of N and Ne from given T and RHO
c
C
C     Input:  T   - temperature
C             RHO - mass density
C     Output: AN  - total particle density
C             ANE - elctron density
C
      INCLUDE 'PARAMS.FOR'
      common/nerela/anerel
c
      it=0
      if(id.eq.1.and.anerel.eq.0.) then
         anerel=0.5
         if(t.lt.9000.) anerel=0.4
         if(t.lt.8000.) anerel=0.1
         if(t.lt.7000.) anerel=0.01
         if(t.lt.6000.) anerel=0.001
         if(t.lt.5500.) anerel=0.0001
c         if(t.lt.5000.) anerel=1.e-5
c         if(t.lt.4000.) anerel=1.e-6
      end if
   10 continue
      it=it+1
      an=rho/wmm(id)/(1.d0-anerel)
      ane0=anerel*an
      call eldens(id,t,an,ane)
      anerel=ane/an
      write(6,602) it,id,t,rho,an,ane,wmm(id),anerel
  602 format(/' **** rhonen it,id,t,r,N,Ne,wmm,ner',2i4,f7.0,1p5e11.4)
      if(abs((ane-ane0)/ane0).lt.1.e-5) go to 20
      if(it.lt.50) go to 10
c      write(6,601) an,ane,ane0
c 601 format(/' slow convergence of RHONEN - N,Ne,Nep=',1p3e11.3)
   20 continue
c
      return
      end