FUNCTION nint,x
  ix=long(x+.5)
  j=where(x lt 0,n)
  if n gt 0 then ix(j)=ix(j)-1
  return,ix
END

FUNCTION lookup,lnrho,ee,dvdecr,dvdrce,iv=iv,mz=mz,status=status

      common cdat,x,y,z,nx,ny,nz,mp,ntmax,date0,time0
      common ctable,ttmean,rhm,drhm,eemin,eemax,dee,mtab,itab,tab

      if n_params() eq 0 then begin
	print,'lookup,lnrho,ee[,dvdecr,dvdrce],mz=mz,iv=iv'
	print,'  dvdecr is the partial der of the answ wrt ee at const lnr'
	print,'  dvdrce is the partial der of the answ wrt lnr at const ee'
	return,0
      endif

      if n_elements(iv) eq 0 then iv=0
      if n_elements(mz) eq 0 then mz=n_elements(rhm)

      rhm1=alog(rhm(0))
      rhm2=alog(rhm(mz-1))
      drhm=(rhm2-rhm1)/(mz-1)
      s=size(lnrho)
      lnrho1=lnrho
      ee1=ee
      rhm11=lnrho1
      drhm1=rhm11
      if s(0) eq 3 then begin
        rhm11(*,*,*)=rhm1
        drhm1(*,*,*)=drhm
      endif else if s(0) eq 2 then begin
        rhm11(*,*)=rhm1
        drhm1(*,*)=drhm
      endif else begin
        rhm11=rhm1
        drhm1=drhm
      endelse
;
;  density index
;
      algrk=(lnrho1-rhm11)/drhm1
      np=((long(algrk) > 0) < (mz-2))
      px=algrk-np
      qx=1.-px
      pxqx=px*qx

      pxmax = max(px)
      qxmax = max(qx)
      status=0
      if pxmax gt 1 or qxmax gt 1 then begin
	status=status+1
	print,'warning: extrapolation in lnrho, pxmax =',pxmax,' qxmax =',qxmax
        px=px < 1
        qx=qx < 1
      endif
;
;  Interpolate the vars and dvar/dlnrho at smaller lnrho
;
      ntab=mtab(np)
      ltab=3*ntab
      eek=(ee1-eemin(np))/dee(np)
      kee=((long(eek) > 0) < (ntab-2))
      py=eek-kee
      qy=1.-py
      pyqy=py*qy

      pymax = max(py)
      qymax = max(qy)
      if pymax gt 1 or qymax gt 1 then begin
	status=status+1
	print,'warning: extrapolation in ee, pymax =',pymax,' qymax =',qymax
        print,'ee1' & stat,ee1
        print,'eemin' & stat,eemin(np)
;       print,'eek' & stat,eek
;       print,'kee' & stat,kee
        py=py < 1
        qy=qy < 1
      endif

      ik=itab(np)+kee-1
      ntab1=mtab(np+1)
      eek1=(ee1-eemin(np+1))/dee(np+1)
      ;kee1=((long(eek1) > 0) < (ntab1-2))
      kee1=kee-nint((eemin(np+1)-eemin(np))/dee(np))
      py1=eek1-kee1
      if max(abs(py1-py)) gt 0.1 then begin
;	print,'address problem in lookup'
;	stop
        print,'warning: extrapolating in ee, py1max =',max(py1)
      end
      kee1=(kee1 > 0) < (ntab1-2)
      ik1=itab(np+1)+kee1-1
;
;  Get ik for variable desired
;
      ik=ik+iv*3*ntab
      ik1=ik1+iv*3*ntab1
;
      dfx1=tab(ik1)-tab(ik)
      dfx2=tab(ik1+1)-tab(ik+1)
      dfy1=tab(ik+1)-tab(ik)
      dfy2=tab(ik1+1)-tab(ik1)
;
      tabvar = $
          qy * (qx * (tab(ik           ) + $
               pxqx* (tab(ik +2*ntab   ) - dfx1) + $
               pyqy* (tab(ik +  ntab   ) - dfy1) ) + $
                px * (tab(ik1          ) - $
               pxqx* (tab(ik1+2*ntab1  ) - dfx1) + $
               pyqy* (tab(ik1+  ntab1  ) - dfy2) ) ) + $
          py * (px * (tab(ik1        +1) - $
               pxqx* (tab(ik1+2*ntab1+1) - dfx2) - $
               pyqy* (tab(ik1+  ntab1+1) - dfy2) ) + $
                qx * (tab(ik         +1) + $
               pxqx* (tab(ik +2*ntab +1) - dfx2) - $
               pyqy* (tab(ik +  ntab +1) - dfy1) ) )

;
;  Energy derivative
;
      if n_params() lt 3 then return,tabvar
      dvdecr = $
               qx * (   - (tab(ik           )       ) + $
     (1.-4.*py+3.*py^2) * (tab(ik +  ntab   ) - dfy1) + $
                (-pxqx) * (tab(ik +2*ntab   ) - dfx1) + $
                          (tab(ik         +1)       ) - $
        (2.*py-3.*py^2) * (tab(ik +  ntab +1) - dfy1) + $
                   pxqx * (tab(ik +2*ntab +1) - dfx2) ) + $
               px * (     (tab(ik1        +1)       ) - $
        (2.*py-3.*py^2) * (tab(ik1+  ntab1+1) - dfy2) - $
                   pxqx * (tab(ik1+2*ntab1+1) - dfx2) - $
                          (tab(ik1          )       ) + $
     (1.-4.*py+3.*py^2) * (tab(ik1+  ntab1  ) - dfy2) - $
                (-pxqx) * (tab(ik1+2*ntab1  ) - dfx1) )
;
      dvdecr = dvdecr/dee(np)
      ; stat,dvdecr
      if n_params() lt 4 then return,tabvar
;  Density derivative
;
      dvdrce = $
               qy * (   - (tab(ik           )       ) + $
     (1.-4.*px+3.*px^2) * (tab(ik +2*ntab   ) - dfx1) + $
                (-pyqy) * (tab(ik +  ntab   ) - dfy1) + $
                          (tab(ik1          )       ) - $
        (2.*px-3.*px^2) * (tab(ik1+2*ntab1  ) - dfx1) + $
                   pyqy * (tab(ik1+  ntab1  ) - dfy2) ) + $
               py * (     (tab(ik1        +1)       ) - $
        (2.*px-3.*px^2) * (tab(ik1+2*ntab1+1) - dfx2) - $
                   pyqy * (tab(ik1+  ntab1+1) - dfy2) - $
                          (tab(ik         +1)       ) + $
     (1.-4.*px+3.*px^2) * (tab(ik +2*ntab +1) - dfx2) - $
                (-pyqy) * (tab(ik +  ntab +1) - dfy1) )
;
;  Dimensional factors
;
      dlnrho = drhm
      dvdrce = dvdrce/dlnrho
      ; stat,dvdecr
;
      return,tabvar
end