pro readmol, name, lev, rad, temp, col, wgt, npart=npart, verbose=verbose

; Read a molecular data file in the LAMDA (Leiden atomic and molecular 
; database) format.
;
; name = file name
; nlev = number of levels in data file
; lev = info on each level; column 1 = level energy (in cm^-1, i.e. h c / E),
;	column 2 = statistical weight, column 3 = J, column 4 = K, if 
;	applicable. Note that transitions are sorted by spectroscopic
;	classification, not by energy.
; rad = info on each radiative transition; column 1 = upper level, column
;	2 = lower level, column 3 = Einstein A (in s^-1), column 4 =
; 	transition frequency (in GHz), column 5 = energy of upper level
;	(in K)
; temp = table of temperatures at which collision strength is tabulated;
;	if npart > 1, an array of npart pointers to arrays is returned
; col = info on each collision; column 1 = upper level, column 2 = lower
;	level, columns 3 to 3+ntemp-1 = collision rate (in cm^-3 s^-1)
;	at the given temperatures; if npart > 1, an array of npart pointers
;	to arrays of this form is returned

; return on error
On_Error, 2

; open file
openr, fp, name, /get_lun

; skip first 3 lines
dummy=''
for n=0,2 do readf, fp, dummy

; read molecular weight, then skip
wgt=0.0
readf, fp, wgt
readf, fp, dummy

; read number of levels, then skip a line
nlev=0
readf, fp, nlev
readf, fp, dummy

; see if we have just J, or J, K levels
pos=strpos(dummy, "J,")

; read level data
if pos eq -1 then begin
	lev=fltarr(4,nlev)
	readf, fp, lev
	lev=lev[1:3,*]
endif else begin
	lev=fltarr(5,nlev)
	for n=0,nlev-1 do begin
		readf, fp, dummy
		parts=strsplit(dummy, /extract)
		lev[0,n]=fix(parts[0])
		lev[1,n]=float(parts[1])
		lev[2,n]=float(parts[2])
		parts1=strsplit(parts[3], '_', /extract)
		lev[3,n]=fix(parts1[0])
		lev[4,n]=fix(parts1[1])
	endfor
	lev=lev[1:4,*]
end

; skip a line, read number of radiative transitions, skip a line
nrad=0
readf, fp, dummy
readf, fp, nrad
readf, fp, dummy

; read transition data
rad=fltarr(6,nrad)
readf, fp, rad
rad=rad[1:5,*]

; subtract 1 from all state numbers to convert them to 0 offset,
; instead of 1 offset
rad[0:1,*]=rad[0:1,*]-1

; read number of collision partners and create array of collisions
npart=0
readf, fp, dummy
readf, fp, npart
if npart gt 1 then begin
	ntemp=intarr(npart)
	temp=ptrarr(npart)
	col=ptrarr(npart)
endif

; loop over partners
for n=0,npart-1 do begin

	; read number of collisions and temperatures
	ncol=0
	ntemp1=0
	readf, fp, dummy
	readf, fp, dummy
	readf, fp, dummy
	readf, fp, ncol
	readf, fp, dummy
	readf, fp, ntemp1
	readf, fp, dummy

	; read temperatures
	temp1=fltarr(ntemp1)
	readf, fp, temp1
	readf, fp, dummy

	; read collision table
	col1=fltarr(3+ntemp1,ncol)
	readf, fp, col1
	col1=col1[1:3+ntemp1-1,*]

	; adjust state numbers to 0 offset
	col1[0:1,*]=col1[0:1,*]-1

	; store data
	if npart eq 1 then begin
		ntemp=ntemp1
		temp=temp1
		col=col1
	endif else begin
		ntemp[n]=ntemp1
		temp[n]=ptr_new(temp1)
		col[n]=ptr_new(col1)
	endelse
endfor

free_lun, fp


end