Método de integração de Ridder em F90 (código)

Esta é uma pré-visualização de arquivo. Entre para ver o arquivo original

! Rootfinder using Ridders method
PROGRAM ridder
 IMPLICIT NONE
 REAL zr
 zr=zriddr(sen,3.,3.2,1.0e-6)
 WRITE(*,*) zr
CONTAINS
REAL FUNCTION zriddr(func, x1, x2, xacc)
! Using Ridders’ method, return the root of a function func known to lie 
! between x1 and x2. The root, returned as zriddr, will be refined to an
! approximate accuracy xacc.
 IMPLICIT NONE
 INTEGER, PARAMETER :: MAXIT=50 ! Maximum allowed number of iterations.
 REAL, PARAMETER :: UNUSED=-1.11e30
 REAL, EXTERNAL :: func !Função para qualquer variável independente da subrotina (EXTERNAL).
 REAL :: x1,x2,xacc
 INTEGER j
 REAL :: fh,fl,fm,fnew,s,xh,xl,xm,xnew
 fl=func(x1)
 fh=func(x2)
 IF(fl*fh < 0.0) THEN !Teste para ver se um ponto é negativo e o outro positivo.
 xl=x1
 xh=x2
 zriddr=UNUSED; ! Any highly unlikely value, to simplify logic
 DO j=1,MAXIT 
 xm=0.5*(xl+xh)
 fm=func(xm); ! First of two function evaluations per iteration 
 s=sqrt(fm*fm-fl*fh)
 IF (s == 0.0) RETURN
 xnew=xm+(xm-xl)*SIGN(fl,1.0)*fm/s ! Updating formula. 
 IF (ABS(xnew-zriddr) <= xacc) RETURN
 zriddr=xnew
 fnew=func(zriddr)
 IF (fnew == 0.0) RETURN
 IF (fnew*fm<0.0) THEN
 xl=xm
 fl=fm
 xh=zriddr
 fh=fnew
 ELSE IF (fnew*fl<0.0) THEN 
 xh=zriddr
 fh=fnew
 ELSE
 xl=zriddr
 fl=fnew
 END IF
 IF (ABS(xh-xl) <= xacc) RETURN
 END DO
 WRITE(*,*) 'Error! zriddr exceeded maximum iterations'
 ELSE
 IF (fl==0.0) THEN 
 zriddr=x1
 RETURN
 ELSE IF (fh==0.0) THEN
 zriddr=x2
 RETURN
 END IF
 WRITE(*,*) 'Error! Root must be bracketed in zriddr.'
 ENDIF
END FUNCTION zriddr
REAL FUNCTION sen(x)
 IMPLICIT NONE
 REAL :: x
 sen=sin(x)
END FUNCTION sen
END PROGRAM ridder