function getlevpop, nbar, kT, beta, no_beta=no_beta

; Compute the level populations for the given molecule for the given
; mean density, temperature, and matrix of escape probabilities beta.
; If no_beta is set, then beta is not used and the escape probabilities
; are all set to unity.

; Block for molecule data
COMMON LINELUM_MOLDATA, $
	molfile, H2opr, nlev, levenergy, levwgt, radupper, radlower, $
	radfreq, radTupper, temptable, colupper, collower, coltable, $
	einstein_a, wgt

On_Error, 2

kb = 1.38d-16

; Initialize matrix to hold collision rates
q = dblarr(nlev, nlev)

; Get interpolated collision rates based on temperature
t=kT/kb
tidx=min(where(t gt temptable))
if tidx lt 0 then tidx = 0
if tidx gt n_elements(temptable-1) then tidx=n_elements(temptable-1)
wgt=alog(t/temptable[tidx]) / alog(temptable[tidx+1]/temptable[tidx])
colrates=reform( coltable[tidx,*] * $
	         (coltable[tidx+1,*]/coltable[tidx,*])^wgt )

; Create matrix of collision rates, both upper to lower and lower
; to upper. Get collision rates from upper to lower out of table, and
; from lower to upper using detailed balance
q[colupper,collower]=colrates
gj=levwgt[colupper]
gi=levwgt[collower]
ej=levenergy[colupper]
ei=levenergy[collower]
q[collower,colupper]=colrates * gj/gi * exp(-(ej-ei)/kT)

; Construct the matrix M such that 
; M_ij = (n_H2 q_ji + beta_ji * A_ji) / [Sum_k (n_H2 q_ik + A_ik)]
if not keyword_set(no_beta) then m=nbar*q+einstein_a*beta $
else m=nbar*q+einstein_a
m1=m
for i=0,nlev-1 do m[*,i]=m[*,i]/total(m1[i,*])

;; Solve matrix equation
;eval=la_eigenproblem(m, eigenvectors=evec)
;; We want the solution with eigenvalue 1.0, normalized to give
;; a total population of 1.0
;dummy=min(abs(eval-1.0), idx)
;levpop=double(evec[*,idx])
;levpop/=total(levpop)

; about 10x faster than the above!
d=indgen(nlev)
m[d,d] -= 1.d
levpop = la_least_square_equality(m, DBLARR(nlev)+1, DBLARR(nlev), DBLARR(1)+1)


; Return level population
return, levpop

end