c     implicit none
c     double precision tc,scatt
c     write(6,*)'tc?'
c     read(5,*)tc
c     call scattf(tc,scatt)
c     write(6,*)scatt
c     end

c This subroutine calculates tc in the homogenoeus
c scattering model. See (ajop6,11)
c input:
c   scatt= xi_0/mean free path
c output:
c   tc=T_c/T_c0
      subroutine tcimp(scatt,tc)
      implicit none
      double precision scatt,tc,t,f,fprime,xh,eps,change,vak
      integer nsplit,iter,ie
      parameter(nsplit=10,vak=3.141592654**2/4.)
      iter=0
      if(scatt.gt.0.2807299454)then
        tc=0.
        goto 200
      end if
      t=1.-vak*scatt
      xh=0.5*scatt
10    iter=iter+1
      f=log(t+xh/float(nsplit))
     x -(1./float(nsplit)**2-t**2/(t*float(nsplit)+xh)**2)/24.
      fprime=float(nsplit)/(t*float(nsplit)+xh)
     x +t*xh/(12.*(t*float(nsplit)+xh)**3)
      do 100 ie=1,nsplit
      eps=float(ie)-0.5
      f=f+1./eps-t/(t*eps+xh)
100   fprime=fprime-xh/(t*eps+xh)**2
      change=-f/fprime
c     write(6,*)'t, change=',t,change
      t=t+change
      if(abs(change).gt.0.000000001)goto 10
      tc=t
200   continue
c     write(6,1)tc,iter,scatt
c1    format(' tc=',f8.6,' iter=',i3,' for scatt=',f8.6)
      return
      end

c This subroutine calculates scatt for a given tc in the homogenoeus
c scattering model. (ajop6.11)
c input:
c   tc=T_c/T_c0
c output:
c   scatt= xi_0/mean free path
      subroutine scattf(tc,scatt)
      implicit none
      double precision scatt,tc,t,f,fprime,xh,eps,change,vak
      integer nsplit,iter,ie
      parameter(nsplit=10,vak=3.141592654**2/4.)
      iter=0
      t=tc
      scatt=(1.-t)/vak
10    iter=iter+1
      xh=0.5*scatt
      f=log(t+xh/float(nsplit))
     x -(1./float(nsplit)**2-t**2/(t*float(nsplit)+xh)**2)/24.
      fprime=0.5*(1./(t*float(nsplit)+xh)
     x    -t**2/(12.*(t*float(nsplit)+xh)**3))
      do 100 ie=1,nsplit
      eps=float(ie)-0.5
      f=f+1./eps-t/(t*eps+xh)
100   fprime=fprime+0.5*t/(t*eps+xh)**2
      change=-f/fprime
c     write(6,*)'scatt, change=',scatt,change
      scatt=scatt+change
      if(abs(change).gt.0.000000001)goto 10
      write(6,1)scatt,iter,tc
 1    format(' scatt=',f8.6,' iter=',i3,' for tc=',f8.6)
      return
      end