      program integra
	implicit real *8 (a-h,o-z)
	external f1,f2,f3,f4,f5,f6
      open(unit=1,file='inpint.txt',status='old')
	open(unit=2,file='outint.txt',status='unknown')
	open(unit=3,file='distnorm.txt',status='unknown')
	open(unit=4,file='curvanorm.txt',status='unknown')
!
!
!
	read(1,*)ndiv,a,b
	write(*,*)'ndiv=',ndiv,'a=',a,'b=',b
!
!
!
	call simp(f1,10.42d0,3.d0,100,c1)
	write(*,*)'o resultado da integral 1-simp eh',c1
	write(2,*)'o resultado da integral 1-simp eh',c1
	call trap(f1,10.42d0,3.d0,100,c1)
	write(*,*)'o resultado da integral 1-trap eh',c1
	write(2,*)'o resultado da integral 1-trap eh',c1
	call simp(f2,3.d0,0.61d0,100,c2)
	write(*,*)'o resultado da integral 2-simp eh',c2
	write(2,*)'o resultado da integral 2-simp eh',c2
	call trap(f2,3.d0,0.61d0,100,c2)
	write(*,*)'o resultado da integral 2-trap eh',c2
	write(2,*)'o resultado da integral 2-trap eh',c2
	call simp(f3,0.61d0,0.81d0,100,comb)
	write(*,*)'o resultado da integral 3-simp eh',comb
	write(2,*)'o resultado da integral 3-simp eh',comb
	call trap(f3,0.61d0,0.81d0,100,comb)
	write(*,*)'o resultado da integral 3-trap eh',comb
	write(2,*)'o resultado da integral 3-trap eh',comb
	call simp(f4,-5.d0,5.d0,100,res1)
	write(*,*)'o resultado da integral 4-simp eh',res1
	write(2,*)'o resultado da integral 4-simp eh',res1
	call trap(f4,-5.d0,5.d0,100,res1)
	write(*,*)'o resultado da integral 4-trap eh',res1
	write(2,*)'o resultado da integral 4-trap eh',res1

	call simp(f5,0.81d0,10.42d0,100,exp)
	write(*,*)'o resultado da integral 5-simp eh',exp
	write(2,*)'o resultado da integral 5-simp eh',exp
	call trap(f5,0.81d0,10.42d0,100,exp)
	write(*,*)'o resultado da integral 5-trap eh',exp
	write(2,*)'o resultado da integral 5-trap eh',exp

	wliq=c1+c2+comb+exp
	write(*,*)'o trabalho liquido do motor ICO eh',wliq
	write(2,*)'o trabalho liquido do motor ICO eh',wliq


	call simp(f6,-1.d0,6.d0,6,res1)
	write(*,*)'o resultado da integral 6-simp eh',res1
	write(2,*)'o resultado da integral 6-simp eh',res1
	call trap(f6,-1.d0,6.d0,6,res1)
	write(*,*)'o resultado da integral 6-trap eh',res1
	write(2,*)'o resultado da integral 6-trap eh',res1

!	graficos da curva e da distribuicao normal

	pi=4*atan(1.)

	dx=(b-a)/ndiv
	xb=a+dx
	temp=0.
	write(3,*)a,temp
	temp=f4(a)
	write(4,*)a,temp
	do i=1,ndiv
		call simp(f4,a,xb,50,res)
		temp=res/sqrt(2*pi)
		write(3,*)xb,temp
		temp=f4(xb)/sqrt(2*pi)
		write(4,*)xb,temp
		xb=xb+dx	
	enddo

	close(2)
	close(3)
	close(4)

	call system("notepad curvanorm.txt")

	call system("notepad distnorm.txt")

	call system("wgnuplot dados-curva.gnu")


	stop
	end
!	
!***********************************
!
	subroutine simp(func,a,b,n,sum)
	implicit real *8 (a-h,o-z)
	external func

	sum=0.d0
	h=(b-a)/n
	nend=n/2
	do i=1,nend
		a1=(2*i-2)*h+a
		a2=(2*i-1)*h+a
		a3=2*i*h+a
		fa=func(a1)
		fb=func(a2)
		fc=func(a3)
		sum=sum+(fa+4*fb+fc)
	enddo
	sum=h*sum/3
	return
	end
!
!*************************
!
	subroutine trap(func,a,b,n,sum)
	implicit real *8 (a-h,o-z)
	external func

	sum=0.d0
	h=(b-a)/n
	x=a
	do i=1,n
		sum=sum+(func(x)+func(x+h))
		x=x+h
	enddo
	sum=h*sum/2
	return
	end
!
!**************************
	function f1(x)
	implicit real *8 (a-h,o-z)

	! f1=cos(2.d0*x)

!	f1=cos(x)/exp(x)+asin(x)

!	f1=cos(x)/exp(x)+2.d0**(-x*x)

	f1=1.31-0.128*(x-10.42)		 ! MCI - compressao 1

	return
	end

!*************************
	function f2(x)
	implicit real *8 (a-h,o-z)

!	f2=x**2

!	f2=log(x)*sin(x)/x**2+exp(-x*x/2)

	f2=2.26-1.14*(x-3.)+3.21*(x-3.)*(x-2.5)	 ! MCI - compressao 2

	return
	end
!************************
	function f3(x)
	implicit real *8 (a-h,o-z)

!	f3=sin(x)/x+atan(x)

!	f3=sin(x)/x+atan(x*x)-log(x)

	f3=19.5+60.*(x-0.61)-340.*(x-0.61)*(x-0.71)  ! MCI - combustao

	return
	end
!***********************
	function f4(x)
	implicit real *8 (a-h,o-z)

	f4=exp(-x*x/2)

!	f4=24.7-2.78*(x-0.81)+0.255*(x-0.81)*(x-9.07)  ! MCI - expansao

	return
	end
!***********************
	function f5(x)
	implicit real *8 (a-h,o-z)

!	f5=sin(x)*sin(x)*log(x)

	f5=24.7-2.78*(x-0.81)+0.255*(x-0.81)*(x-9.07)  ! MCI - expansao

	return
	end
!***********************
	function f6(x)
	implicit real *8 (a-h,o-z)

	f6=sin(x)/x*exp(-x*x)
	return
	end
!***********************
