!      program rungekutta

	common /param/ pi
	external fcn
      open(1,file='inpaula.txt')
	open(2,file='outnum.txt')
	open(3,file='outteo.txt')
	open(4,file='outexemplo.txt')
!
!	integrador de EDO's com passo fixo - RK 4a ordem
!
	read(1,*)n,tau0
	write(*,*)'n=',n,'tau0=',tau0
	read(1,*) tend,dtau,teta0
	write(*,*)'tend=',tend,'dtau=',dtau,'teta0=',teta0
!
!	calculo de pi
	pi=4*atan(1.d0)
	write(*,*)'pi=',pi

!
!	initial values
!
	erro=0.
	teta=teta0
	tetateo=teta0
	time=tau0
	k=0 
	write(*,*) '     Tabela de resultados'
	write(*,*)'---------------------------------------------------------------------'
	write(*,*)'  Passo Nr     t          x(t)          xteo(t)          Erro (ETL)'
	write(*,*)'---------------------------------------------------------------------'
	write(*,*)k,time,teta-273.15 ! ,tetateo,erro

	write(4,*) '     Tabela de resultados'
	write(4,*)'---------------------------------------------------------------------'
	write(4,*)'  Passo Nr     t          x(t)          xteo(t)          Erro (ETL)'
	write(4,*)'---------------------------------------------------------------------'
	write(4,*)k,time,teta-273.15 ! ,tetateo,erro
	write(2,*)time,teta-273.15
!	write(3,*)time,teta
 50	k=k+1 
	th=time+0.5*dtau 
	f1=dtau*fcn(time,teta)
	f2=dtau*fcn(th,teta+0.5*f1)
	f3=dtau*fcn(th,teta+0.5*f2)
	f4=dtau*fcn(time+dtau,teta+f3)
	teta=teta+(f1+2.*f2+2.*f3+f4)/6.
	time=time+dtau
	if (time.le.tend) then
			write(*,*)k,time,teta-273.15
			write(2,*)time,teta-273.15
		!	tetateo=exp(-time**2/2)  
		!	tetateo=time/(0.5+log(time))
		!	erro=abs(teta-tetateo)
		!	write(3,*)time,tetateo
			write(4,*)k,time,teta-273.15  ! ,tetateo,erro
			goto 50 
		else
			goto 300
		endif

300 continue
	close(2)
	close(3)
	close(4)


  call system('notepad outnum.txt') 

 ! call system('notepad outteo.txt') 

  call system('notepad outexemplo.txt') 

  call system('wgnuplot dados.gnu') ! gráfico

	stop
	end
!**************************************
	function fcn(t,f)
	common /param/ pi

	! avalia o valor da derivada em cada passo
	! aux=-t*f  
	! aux=(t*f-f**2)/t/t
	!	aux=exp(2*t)+2*f

	fi=0.019
	d=0.1
	rl=0.1
	qr=5400.
	hint=10.
	hext=1.
	eiso=0.01
	rk_iso=0.1
	tinf=293.15
	rmdot_ag=0.08
	c_ag=4180.
	rho_ag=1000.

!	calculo de U
	u=1./(1./hint+eiso/rk_iso+1./hext)

!	calculo da area lateral
	area_lat=pi*d*rl+pi*(d**2-fi**2)/2

!	calculo da perda de calor
	qinf=u*area_lat*(tinf-f)

! calculo de qconv - fluxo de agua
	qconv=rmdot_ag*c_ag*(tinf-f)

!	calculo da massa de agua
	vol=pi*d**2/4*rl
	rm_ag=rho_ag*vol

	aux=(qr+qinf+qconv)/(rm_ag*c_ag)

	fcn=aux
	return
	end
!***************************************
