subroutine istat1(p2,t2,iw,iswch,converg)
*
* compute the thermodynamic variables for a specific P, T
* point.  This version searches either the C or O eos table, depending
* on the composition and returns the thermodynamic properties.
*
      implicit double precision(a-h,o-z)

      common/temp/s1,r1,s2,r2,b2,px2,tx2,e2,xc2,xo2,fca2,fx2,q2,w2,c
      common/tablec/tablec(6,39,87)
      common/tableo/tableo(6,36,87)
      common/sizec/tminc,tmaxc,deltc,pminc,pmaxc,delpc,pmc(50),ipc(50),
     1   iphasec(50),ioverc,itc
      common/sizeo/tmino,tmaxo,delto,pmino,pmaxo,delpo,pmo(50),ipo(50),
     1   iphaseo(50),iovero,ito
      common/thermo/d2,dp,dt,e,ep,et,psi,pg,o2,op,ot,fp,ft, 
     1   en2(5),fn2,fcc2,ce,ci,cif,w,nu 
      common/neut/en(5),fn
      common/iii/ iii
      common/ixswch/ixswch
      common/je/je

      logical converg
      dimension porder(4),torder(4),iorder(4),jorder(4),
     1          xorder(4), yorder(4)
      dimension thermo(24)
      equivalence (d2,thermo(1))
      data delta/2.e-3/
      data fnat/4.342945e-1/
*
* initialize some variables
*
      xc = xc2
      xo = xo2
      je=3
      iii=0
      p=p2
      t=t2
      if ( iswch .eq. 12 ) then
         ci=1./(xc/12.+xo/16.) 
      else if ( iswch .eq. 16 ) then
         ci=1./(xo/16.+(1-xo)/21.45)
      else
         stop
      endif
*
* find eos box and interpolate.
* We search carbon if ISWCH is 12 and oxygen if ISWCH is 16.
*
      if( iswch .eq. 12 ) then
        call wdselect(t,tminc,deltc,torder,iorder,1,itc,converg)
        kmax=6
        if( iw .ge. 1 )  kmax=1
          do 5 k=1,kmax
            do 4 i=1,4
              indexc=iorder(i)
              pfrzc = phasec(t)
              if( p .lt. pfrzc ) then 
                jminc=iphasec(indexc)+1
                jmaxc=ipc(indexc)
                pstartc=pmaxc-ioverc*delpc
               else 
                jminc=1
                jmaxc=iphasec(indexc) 
                pstartc=pmaxc
              endif
 
              call wdselect(p,pstartc,delpc,porder,jorder,jminc,jmaxc, 
     1          converg)
 
              do 3 j=1,4
                xorder(j)=tablec(k,indexc,jorder(j))
3             continue
              call dali(p,porder,xorder,yorder(i))
4           continue
            l=k
            if(k.gt.4) l=k+1
            call dali(t,torder,yorder,thermo(l))
5        continue
*
*  or if oxygen, then come here
*
      else if ( iswch .eq. 16 ) then
         call wdselect(t,tmino,delto,torder,iorder,1,ito,converg)
         kmax=6
         if ( iw .ge. 1 )  kmax=1
         do 35 k=1,kmax
            do 34 i=1,4
               indexo=iorder(i)
               pfrzo = phaseo(t)
               pfrzc = pfrzo
               if ( p .lt. pfrzo ) then 
                  jmino=iphaseo(indexo)+1
                  jmaxo=ipo(indexo)
                  pstarto=pmaxo-iovero*delpo
               else 
                  jmino=1
                  jmaxo=iphaseo(indexo) 
                  pstarto=pmaxo
               endif
 
               call wdselect(p,pstarto,delpo,porder,jorder,jmino,jmaxo, 
     1          converg)
 
               do 33 j=1,4
                  xorder(j)=tableo(k,indexo,jorder(j))
33             continue
               call dali(p,porder,xorder,yorder(i))
34          continue
            l=k
            if(k.gt.4) l=k+1
            call dali(t,torder,yorder,thermo(l))
35       continue
      endif
      d  = d2
*
* Now that we've interpolated the main quantities, we need to
* calculate the opacity, carbon burning and neutrino luminosities
* and their derivitives wrt rho and t
* Note that the carbon burning luminosity is reduced by a factor
* of 5 arbitrarily (otherwise, runaway).
*
      n  = -1
20    x1=49.4+d-2*t/3-36550*(1+7e-1*10.**(t-10))**(1./3)*10.**(-t/3)
 
      if ( x1 .le. -30. ) then
         fcc=0
      else
         xtmp = .999999
         fcc=4.368*xtmp*xtmp*10.**x1
         fcc=fcc/5. 
      endif
 
      if (iw .ge. 1) then
         fcc2=fcc
         return
      endif
 
      if (nu .le. 0) then
         fn=0
      elseif(nu .eq. 1)then
*
* call venerable BPS neutrino rates
*
         call neu(d,t)
      elseif(nu .eq. 2)then
*
* call Munkata etal Plasmon, Photo, and Pair rates
*
         call neutrino(d,t)
      elseif(nu .eq. 3)then
*
* call Itoh etal neutrino rates
*
         call gnutrino(d,t)
      endif

      f=fn-fcc*w
      o=10.**opac(d,t)
      if(n) 14,15,16

14    o2=o
      f2=f
      fn2=fn
      fcc2=fcc
      do 17 i=1,5
17    en2(i)=en(i)
      e  = 10.**e
      et = e * et
      ep = (dt/fnat)*10.**(p-d-t)
      pg=dlog10(10.**p-10.**(4*t-14.59836))
      d=d2+delta
      n=0 
      go to 20

15    od=(o-o2)/delta
      fd=(f-f2)/delta
      d=d2
      t=t2+delta
      n=1 
      go to 20

16    ot=(o-o2)/delta
      ft=(f-f2)/delta
      op=od*dp
      ot=ot+od*dt
      fp=fd*dp
      ft=ft+fd*dt

      return
      end 

************************************************************************