c------------------------------------------------------------------------------
c ... code to prepare tables of transport coefficients partial
c ... differentials with respect to temperature (e.g. amu_t = partial mu wrt T
c ... flag: in_n2 has the following purpose
c ... in_n2 = .false. ==> transport properties are evaluated based on
c ...                     the 1D flame mixture composition 
c ...                     (of course the binary diffusion coefficients
c ...                     are D_i_N2, but given the rest of the mixture)
c ...                     so: amu_t(T) = amu_t_1D-flame-mixture (T)
c ...                         htcond_t(T) = htcond_t_1D-flame-mixture (T)
c ...                         cp_t(T) = cp_t_1D-flame-mixture (T)
c ...
c ... in_n2 = .true.  ==> transport properties are evaluated based on 
c ...                     an N2 gas only.
c ...                     so: amu_t(T)    = amu_t_N2 (T)
c ...                         htcond_t(T) = htcond_t_N2 (T)
c ...                         cp_t(T)     = cp_t_N2 (T)
c ...                         
c------------------------------------------------------------------------------
c ... unaffected by in_n2: are:
c ... cpsp_k_t (T), hsp_k_t (T), gsp_k_t (T)
c ... capd_k_t(T) = capd_k_N2_t (T) binary diffusion coeff. indep of mixture
c ... anyway. At any rate, mechanism of in_n2 coded in gcapd also.
c ... gcapd gives identical answers for in_n2 false or true
c------------------------------------------------------------------------------
c ... the only meaningful result of in_n2=.true. is to change
c ... the tabulated/used transport coefficients for momentum and heat
c ... namely: mu (viscosity) and lambda (thermal conductivity) of the
c ... mixture
c------------------------------------------------------------------------------

	subroutine dfitder(temp_min,temp_max)

        include 'parameter.par' 
#include "precision.com"
	include 'common.com'

	character*1024 fmt

c-------------------------------------------------------------------------------
	epstemp = 1.e-2		! epsilon Temp in K
	epso2   = 0.5*epstemp
c-------------------------------------------------------------------------------c ... read chemistry input file and initialize chemical data bases

	call cheminit
	call setcon

c-------------------------------------------------------------------------------
c ... read input data

	call init		! initialize 

c------------------------------------------------------------------------------
c ... evaluate detailed transport and thermodynamic properties
c ... these are dimensional quantities in SI units

	call getmu		! get amu_t(j), j=1,2,...,n
	call gthtcond		! get htcond_t(j), j=1,2,...,n
	call getcpsp            ! get cpsp_t(j,k), j=1,2,...,n, k=1,2,...,nsp
	call gcapd              ! get D_t(j,k)_nsp, j=1,2,...,n, k=1,2,...,nsp
	call gethsp             ! get hsp_t(j,k), j=1,2,...,n, k=1,2,...,nsp
	call getgsp             ! get gsp_t(j,k), j=1,2,...,n, k=1,2,...,nsp

c------------------------------------------------------------------------------
c ... writeout

	write(*,*)' writing out data'

c	call chk

c------------------------------------------------------------------------------
c ... fit each property against temperature
c------------------------------------------------------------------------------
c------------------------------------------------------------------------------
c ... set the span of the tabulation of temperature

	t1 = temp_min
	t2 = temp_max

c ... set the number of points in the table

	dt = 0.2      
c	dt = 2.0              !Use this when making tables for the demo

	ng = nint((t2-t1)/dt)+1
	if (ng.gt.ntabmx) then
	   write(*,*) 'dfit_t: Error you need to increase ntabmx or dicrease '
	   write(*,*) 'the temperature range' 
	endif
	nc = ng - 1

c ... tabulate temperature 

	do k=1,ng
	   tv      = t1 + dble(k-1)*dt
	   ttab(k) = tv
	end do

c------------------------------------------------------------------------------
c ... fit cp_k(T)

	open(unit=50,file='cpsp_T.tab',status='unknown')
	write(50,'(''# 1:T 2:cp_1_T 3:cp_2_T ... nsp+1:cp_nsp_T'')')
	write(50,'(''# '',e22.14)')dt
	write(50,'(''# '',2i9)')ng,nsp+1
	do k=1,ng
	   vtemp = ttab(k)
           call getcpor(vtemp+epso2,nsp,cpspj_p)
           call getcpor(vtemp-epso2,nsp,cpspj_m)
	   do kk=1,nsp
	      cpspj(kk) = rg(kk)*(cpspj_p(kk)-cpspj_m(kk))/epstemp
	   end do
	   write(fmt,'(a,i3,a)') '(', nsp+1, 'e22.14)'
	   write(50,fmt) ttab(k),(cpspj(kk),kk=1,nsp)
	end do
	close(50)

c------------------------------------------------------------------------------
c ... fit h_k (T)

	open(unit=50,file='hk_T.tab',status='unknown')
	write(50,'(''# 1:T 2:h_1_T 3:h_2_T ... nsp+1:h_nsp_T'')')
	write(50,'(''# '',e22.14)')dt
	write(50,'(''# '',2i9)')ng,nsp+1
	do k=1,ng
	   vtemp = ttab(k)
           call gethort(vtemp+epso2,nsp,hspj_p)
           call gethort(vtemp-epso2,nsp,hspj_m)
	   do kk=1,nsp
	      hspj_p(kk) = hspj_p(kk)*rg(kk)*(vtemp+epso2)
	      hspj_m(kk) = hspj_m(kk)*rg(kk)*(vtemp-epso2)
	      hspj(kk)   = (hspj_p(kk)-hspj_m(kk))/epstemp
	   end do
	   write(fmt,'(a,i3,a)') '(', nsp+1, 'e22.14)'
	   write(50,fmt)ttab(k),(hspj(kk),kk=1,nsp)
	end do
	close(50)

c------------------------------------------------------------------------------
c ... fit g_k (T)

	open(unit=50,file='gk_T.tab',status='unknown')
	write(50,'(''# 1:T 2:g_1_T 3:g_2_T ... nsp+1:g_nsp_T'')')
	write(50,'(''# '',e22.14)')dt
	write(50,'(''# '',2i9)')ng,nsp+1
	do k=1,ng
	   vtemp = ttab(k)
           call getgibbs(vtemp+epso2,nsp,gspj_p)
           call getgibbs(vtemp-epso2,nsp,gspj_m)
	   do kk=1,nsp
	      gspj(kk) = ( gspj_p(kk) - gspj_m(kk) ) / epstemp 
	   end do
	   write(fmt,'(a,i3,a)') '(', nsp+1, 'e22.14)'
	   write(50,fmt)ttab(k),(gspj(kk),kk=1,nsp)
	end do
	close(50)

c------------------------------------------------------------------------------

	return
	end

c*******************************************************************************