c ctrapezoid.f c c
program trapezoid
c
c    Test Trapezoidal Rule of Integration
c
c      John Mahaffy 4/1/95
c
implicit none
c
c   PARAMETERS
c    ntrys - Number of integrations, each with half the x step of the last
c    nint1 - Number of steps between lowest and highest x values
c    npoints - Largest number of points evaluated for use in the integration
c
c   Variables
c    x   -   x values at which the function is evaluated
c    y   -   values of the function for corresponding values of x
c    np  -   number of points currently used in x and y
c    ansreal - Actual Value of the integral
c    ansnum  - Value of the integral obtained from a Trapezoidal Rule
c
integer ntrys,nint1,npoints
parameter (ntrys=8,nint1=10,npoints=nint1*2**(ntrys-1)+1 )
c
common x(npoints),y(npoints)
common/misc/np
c
real x,y,ansreal,ansnum
integer np
integer i
do 10 i=1,ntrys
np=nint1*2**(i-1)+1
call setcurve
call integrate
call compare
10  continue
stop
c
end
c
subroutine setcurve
c
c    Set up a Curve and range for integration
c    Here we will integrate sin(x) from x=0 to x=3.0
c
c      John Mahaffy 4/1/95
c
implicit none
integer ntrys,nint1,npoints
parameter (ntrys=8,nint1=10,npoints=nint1*2**(ntrys-1)+1 )
common x(npoints),y(npoints)
common/misc/np
real x,y,ansreal,ansnum
integer np
c
integer i
real dx,xlow,xhigh
data xlow,xhigh/0.0,3.0/
c
dx=(xhigh-xlow)/(np-1)
x(1)=xlow
y(1)=sin(x(1))
do 10 i=2,np
x(i)=x(i-1)+dx
y(i)=sin(x(i))
10     continue
c
ansreal=cos(xlow)-cos(xhigh)
c
return
end
subroutine integrate
c
c    Use Trapezoidal Rule to Integrate the area under the curve
c    specified by the data points in x and y
c
c      John Mahaffy 4/1/95
c
implicit none
integer ntrys,nint1,npoints
parameter (ntrys=8,nint1=10,npoints=nint1*2**(ntrys-1)+1 )
common x(npoints),y(npoints)
common/misc/np
real x,y,ansreal,ansnum
integer np
c
c   Local Variables
c
c    esterr   -   estimated error for Trapizoidal integration
c    sum2     -   Trapizoidal integration using every other available point
c
integer i
real sum2, esterr
ansnum=0.0
do 10 i=1,np-1
10    ansnum=ansnum+.5*(y(i)+y(i+1))*(x(i+1)-x(i))
c
c   The following only works if np is an odd integer
c
sum2=0.0
c
do 20 i=1,np-1,2
c
20    sum2=sum2+.5*(y(i)+y(i+2))*(x(i+2)-x(i))
esterr=(sum2-ansnum)*.333333
write(6,2000) np,ansnum,esterr
2000 format (' For ',i4,' points the integral is ',1p,e12.5,
& ', Estimated Error=',e10.2)
return
end
subroutine compare
c
c  Compare results of the numerical integration with the actual result
c
c      John Mahaffy 4/1/95
c
implicit none