subroutine rhonen(id,t,rho,an,ane)
c     ==================================
c
c     iterative determination of N and Ne from given T and RHO
C
C     Input:  T   - temperature
C             RHO - mass density
C     Output: AN  - total particle density
C             ANE - elctron density
C             ENRG- internal energy
C             ENTT- entropy
c
      INCLUDE 'IMPLIC.FOR'
      INCLUDE 'BASICS.FOR'
      INCLUDE 'MODELQ.FOR'
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
         if(t.lt.5000.) anerel=1.e-5
         if(t.lt.4000.) anerel=1.e-6
      end if
      wm=wmm(id)*(un-anerel)/hmass
      wm0=wm
   10 continue
      it=it+1
      an=rho/wm/hmass
      ane0=anerel*an
      wm0=wm
      call eldens(id,t,an,ane,enrgi,entt,wm,0)
      anerel=ane/an
      if(abs((ane-ane0)/ane0).lt.1.e-5.and.
     *   abs((wm-wm0)/wm0).lt.1.e-5) go to 20
      if(it.lt.30) go to 10
   20 continue
      return
      end