00001
00002 c
00003 c Copyright (c) 1986,1987,1988,1989,1990,1991,1992,1993,
00004 c by Steve McMillan, Drexel University, Philadelphia, PA.
00005 c
00006 c All rights reserved.
00007 c
00008 c Redistribution and use in source and binary forms are permitted
00009 c provided that the above copyright notice and this paragraph are
00010 c duplicated in all such forms and that any documentation,
00011 c advertising materials, and other materials related to such
00012 c distribution and use acknowledge that the software was developed
00013 c by the author named above.
00014 c
00015 c THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
00016 c IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
00017 c WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE.
00018 c
00019
00020 c------------------------------------------------------------------------
00021 c
00022 c Integration and differentiation of the input arrays.
00023 c
00024 c------------------------------------------------------------------------
00025
00026 subroutine integrate(x,y,z,n,iprompt,*)
00027 c
00028 c Integrate y with respect to x, using the trapezoid rule,
00029 c assuming that there are no pathologies in the data, placing
00030 c the result in z.
00031 c
00032 dimension x(n),y(n),z(n)
00033 c
00034 if (n.le.0) return
00035 c
00036 if (n.eq.1) then
00037 c
00038 c Result = 0 if only one point is supplied.
00039 c
00040 z(1) = 0.
00041 c
00042 else
00043 sum = 0.
00044 do i=2,n
00045 sum1 = sum + .5*(y(i)+y(i-1))*(x(i)-x(i-1))
00046 z(i-1) = sum
00047 sum = sum1
00048 end do
00049 z(n) = sum
00050 end if
00051 c
00052 end
00053
00054
00055 subroutine differentiate(x,y,z,n,iprompt,*)
00056 c
00057 c Differentiate y with respect to x, assuming that there are no
00058 c pathologies in the data, placing the result in z. Note that
00059 c the array "z" may really be x or y.
00060 c
00061 dimension x(n),y(n),z(n)
00062 c
00063 if (n.le.0) return
00064 c
00065 nerr = 0
00066 if (n.eq.1) then
00067 c
00068 c Result = 0 if only one point is supplied.
00069 c
00070 z(1) = 0.
00071 c
00072 else
00073 c
00074 c Don't (can't) center the first and last derivatives.
00075 c The rest are centered and hence second-order accurate.
00076 c
00077 zz = (y(2) - y(1))/(x(2) - x(1))
00078 c
00079 do i=2,n-1
00080 zz0 = zz
00081 c
00082 dx = x(i+1) - x(i-1)
00083 dy = y(i+1) - y(i-1)
00084 c
00085 if (dx.eq.0.) then
00086 zz = 0.
00087 nerr = nerr + 1
00088 else
00089 zz = dy/dx
00090 end if
00091 c
00092 c Set z(i-1) at the end of the i loop, since location
00093 c i-1 is not going to be used again.
00094 c
00095 z(i-1) = zz0
00096 end do
00097 c
00098 z(n) = (y(n) - y(n-1))/(x(n) - x(n-1))
00099 z(n-1) = zz
00100 c
00101 end if
00102 c
00103 if (nerr.gt.0) then
00104 if (iprompt.ne.0)then
00105 call devoff
00106 write(6,'(i6,a)'),nerr,' errors'
00107 end if
00108 return 1
00109 end if
00110 c
00111 end
