!*************************************
!* Solve HW4, MAE 185. Spring 2006
!* UCI.
!* all blames to Nasser Abbasi
!*
!* Compiler: g95 on windows XP
!* (GCC 4.0.2 (g95) Mar  3 2006)
!************************************
PROGRAM hw4
IMPLICIT NONE

REAL*8 result,solve_problem
INTEGER I

     OPEN(UNIT=999,FILE='hw4_output.txt')
     DO I=1,6
        result=solve_problem(I)
        WRITE(UNIT=999,FMT='(A,I1,A,F14.10)') '\n result problem ',I,'=',result
        WRITE(UNIT=999,FMT='(A)') '============================='
     END DO

END PROGRAM

!*************
!*
!*************
REAL*8 function solve_problem(I)
IMPLICIT NONE

real*8 :: romberg,TOL,xa,xb
real*8,external:: f_1,f_2,f_3,f_4,f_5,f_6
PARAMETER (TOL=1D5*EPSILON(1.0D0))  ! pick some tolerance value
integer :: RMAX,I
parameter (RMAX=12 )

     SELECT CASE(I)
         CASE(1)
            xa=0.0d0
            xb=1.d0              
            solve_problem = romberg(xa,xb,f_1,TOL,RMAX)
         CASE(2)
            xa=0.0d0
            xb=1.d0              
            solve_problem = romberg(xa,xb,f_2,TOL,RMAX)
         CASE(3)
            xa=0.0d0
            xb=1.d0                
            solve_problem = romberg(xa,xb,f_3,TOL,RMAX)
         CASE(4)
            xa=0.0d0
            xb=10.d0              
            solve_problem = xb*LOG(xb) - romberg(xa,xb,f_4,TOL,RMAX)
         CASE(5)
            xa=0.0d0
            xb=1.d0        
            solve_problem = 2*SQRT(xb)*EXP(xb**2) - 4*romberg(xa,xb,f_5,TOL,RMAX)
         CASE(6)
            xa=1.d0
            xb=5.d0        
            solve_problem = romberg(xa,xb,f_6,TOL,RMAX)

     END SELECT      
END  

!************************
!* ROMBERG 
!*
!************************
REAL*8 function romberg(xa,xb,f,EPSILON,RMAX)
IMPLICIT NONE

real*8, external:: f
real*8  :: xa,xb,dx,epsilon,x,s,dI
integer :: RMAX,r,k,m,c

real*8  :: I(0:RMAX,0:RMAX)  ! the romberg table

  r  = 0
  dx = xb-xa
  I(0,0) = dx* ( f(xa)+f(xb) )/2

  WRITE(UNIT=999,FMT='(F14.10)',advance='no') I(0,0)

  DO WHILE( .TRUE. )

    r  = r+1
    dx = dx/2
    k  = 2**(r-1)
    s  = 0

    DO m = 1,k
        x = xa + (2*m-1)*dx
        s = s  + f(x)
    ENDDO
 
    dI = dx*s - I(r-1,0)/2
    I(r,0) = I(r-1,0) + dI
   
    WRITE(UNIT=999,FMT='(F14.10,A)',advance='no') I(r,0),' '

    DO c = 1,r
       I(r,c) = I(r,c-1)+(  I(r,c-1) - I(r-1,c-1) )/(4**c - 1)
       WRITE(UNIT=999,FMT='(F14.10,A)',advance='no') I(r,c),' '
    ENDDO
    IF( r>1 ) WRITE(UNIT=999,FMT='(A)') '' 

    romberg = I(r,r)
    IF( r==RMAX .OR. ABS(romberg - I(r-1,r-1))<EPSILON ) EXIT

  ENDDO

end function

!************************
!* problem 1
!*
!************************
real*8 function f_1(x)
implicit none

real*8 x

f_1=exp(x**2)

end function
!************************
!* problem 2
!*
!************************
real*8 function f_2(x)
implicit none

real*8 x

f_2=x**x

end function
!************************
!* problem 3
!*
!************************
real*8 function f_3(x)
implicit none

real*8 x

f_3=sin(50*x)

end function
!************************
!* problem 4
!*
!************************
real*8 function f_4(x)
implicit none

real*8 x

!f_4=log(x) ===> integrate by parts to change to:

f_4=1

end function
!************************
!* problem 5
!*
!************************
real*8 function f_5(x)
implicit none

real*8 x

f_5 =x**(3./2.)*exp(x**2)

end function
!************************
!* problem 6
!*
!************************
real*8 function f_6(x)
implicit none

real*8 x

f_6 =(x-1)**0.2d0 /(x**2+1)

end function