!********o*********o*********o*********o*********o*********o*********o**
! This program reads a specific file containing rocket thrust data in
! a known format, and approximates the integral using three different
! integration methods.
!---------------------
! Inputs
!   file rocket.dat, containing two lines of gibberish followed by
!   ... some number of lines containing time and thrust
!
! Outputs
!   inputCount, outputCount:        count of lines in each file
!   rectangle, trapezoid, simpsons: integration results
!---------------------
! NOTE: The inputs such as file names and integration limits have been
!       hard-wired into this program;  it is not menu-driven, either.
!       The point is just to showan example of how to do the
!       integration tehcniques in project 2.
!---------------------
! Author:   Bob Harris
! Class:    CMPSC 201F
! Section:  1234
! Date Due: 5 Nov 2003
!********o*********o*********o*********o*********o*********o*********o**

PROGRAM RocketSurgery

character(32) :: inputFilename
character(32) :: outputFilename
integer       :: inputStatus, outputStatus   ! file i/o status
real          :: startTime, endTime          ! integration interval
integer       :: startNum,  endNum           ! integration sample range
integer       :: inputCount, outputCount
real          :: time, thrust                ! values from file
real          :: thrustStart, thrustEnd, thrustEven, thrustOdd
real          :: sumEven, sumOdd             ! sums of thrust values
real          :: dx, rectangle, trapezoid, simpsons

!----------
! open the input file
!----------

! open it (file name is hardwired into program)

inputFilename = "rocket.dat"
open (unit=10, file=inputFilename, &
      iostat=inputStatus, &
      position="rewind", action="read", status="old")
if (inputStatus /= 0) stop "problem opening input file"

! skip the first two lines, which contain text

read (10,*)
read (10,*)

!----------
! open the output file
!----------

! open it (file name is hardwired into program)

outputFilename = "interval.dat"
open (unit=11, file=outputFilename, &
      iostat=outputStatus, &
      action="write", status="replace")
if (outputStatus /= 0) stop "problem opening output file"

!----------
! determine the integration interval
!----------

! hardwire start and end points, in terms of time

startTime = 2.423
endTime   = 3.387

! compute corresponding sample numbers in the input file;  we know that
! the file beings at time 0 and steps time by .001

startNum  = 1 + (startTime * 1000)
endNum    = 1 + (endTime   * 1000)

!----------
! make one pass through the input file, creating the output file and
! computing the integral
!----------

! reset sample counters;  note that every time we read a valid sample,
! we will immediatly increase inputCount;  similarly, every time we
! write a sample, we will immediatly increase outputCount

inputCount  = 0
outputCount = 0

! skip input samples prior to the integration interval

do
    read (10,*) time, thrust
    inputCount = inputCount + 1
    if (inputCount >= startNum) exit
end do

! read the first sample in the integration interval

read (10,*) time, thrustStart
inputCount = inputCount + 1
write (11,*) time, thrustStart
outputCount = outputCount + 1

! read the samples in the interior of the integration interval, and
! sum the odd and even samples separately (we must sum them separately
! to have the right information for simpson's rule)

sumOdd  = 0
sumEven = 0

do
    read (10,*) time, thrustOdd
    inputCount = inputCount + 1
    write (11,*) time, thrustOdd
    outputCount = outputCount + 1

    read (10,*) time, thrustEven
    inputCount = inputCount + 1
    write (11,*) time, thrustEven
    outputCount = outputCount + 1

    sumOdd  = sumOdd  + thrustOdd
    sumEven = sumEven + thrustEven

    if (inputCount > endNum) exit
end do

! read the final sample in the integration interval

read (10,*) time, thrustEnd
inputCount = inputCount + 1
write (11,*) time, thrustEnd
outputCount = outputCount + 1

! compute the integral, with three different formulas

dx = .001
rectangle = (thrustStart   +   sumOdd +   sumEven              ) * dx
trapezoid = (thrustStart/2 +   sumOdd +   sumEven + thrustEnd/2) * dx
simpsons  = (thrustStart   + 4*sumOdd + 2*sumEven + thrustEnd  ) * dx/3

! skip input samples after the integration interval

do
    read (10,*, iostat=inputStatus) time, thrust
    if (inputStatus /= 0) exit
    inputCount = inputCount + 1
end do

! close the files

close (10)
close (11)

!----------
! tell the user what we think of her rocket data
!----------

print *, "there were ", inputCount,  " input data points"
print *, "there were ", outputCount, " output data points"
print *, "rectangle rule says integral is about ", rectangle
print *, "trapezoid rule says integral is about ", trapezoid
print *, "homer simpson  says integral is about ", simpsons

end program RocketSurgery
