PROGRAM series3
!
! Author:  Br. David Carlson
!
! Date:  December 31, 1999
!
! Revised:  May 2, 2004
!
! This program prints out partial sums for a series, stopping when no further
! change is found by adding on another term (or when a set number of terms
! have been added -- to prevent an infinite loop).  It also computes the sum
! in a second way, by adding the terms in reverse order to see if adding the
! smaller ones first is an improvement.  The particular series used is the
! well-known series for sin(x).

IMPLICIT NONE

INTEGER, PARAMETER::max = 100
REAL::x, xsquared, term, sum, sumold
REAL, DIMENSION(max)::hold
INTEGER::n, count, k

WRITE (*, 20)
20 FORMAT (1X, 'Program to investigate a series')

WRITE (*, *) 'Enter an angle in degrees: '
READ (*, *) x

! Convert x to radians:
x = x * ASIN(1.0) / 90.0

sum = 0.0
count = 0
n = 1
term = x
xsquared = x * x

WRITE (*, 50)
50 FORMAT (' ', T5, 'count', T15, 'term', T35, 'partial sum')

DO 
   count = count + 1
   sumold = sum
   hold(count) = term
   sum = sum + term
   WRITE (*, 100) count, term, sum
   IF ((sum == sumold) .OR. (count == max)) THEN
      EXIT
   END IF
   term = -term * xsquared / ((n + 1) * (n + 2))
   n = n + 2
END DO

100 FORMAT (' ', T5, I3, T15, E16.8, T35, E16.8)

WRITE (*, *) ' '
WRITE (*, 200) SIN(x)
WRITE (*, *) ' '
200 FORMAT (1X, 'The built-in sine function gives: ', E16.8)

! Add in reverse order:
sum = 0.0
DO k = count, 1, -1
   sum = sum + hold(k)
   WRITE (*, 300) sum
END DO

300 FORMAT (1X, 'Reverse order partial sum: ', E16.8)

END PROGRAM