diff options
author | Joseph Hunkeler <jhunkeler@gmail.com> | 2015-07-08 20:46:52 -0400 |
---|---|---|
committer | Joseph Hunkeler <jhunkeler@gmail.com> | 2015-07-08 20:46:52 -0400 |
commit | fa080de7afc95aa1c19a6e6fc0e0708ced2eadc4 (patch) | |
tree | bdda434976bc09c864f2e4fa6f16ba1952b1e555 /math/bevington/area.f | |
download | iraf-linux-fa080de7afc95aa1c19a6e6fc0e0708ced2eadc4.tar.gz |
Initial commit
Diffstat (limited to 'math/bevington/area.f')
-rw-r--r-- | math/bevington/area.f | 79 |
1 files changed, 79 insertions, 0 deletions
diff --git a/math/bevington/area.f b/math/bevington/area.f new file mode 100644 index 00000000..ae89d35a --- /dev/null +++ b/math/bevington/area.f @@ -0,0 +1,79 @@ +c function area.f +c +c source +c Bevington, pages 272-273. +c +c purpose +c integrate the area beneath a set of data points +c +c usage +c result = area (x, y, npts, nterms) +c +c description of parameters +c x - array of data points for independent variable +c y - array of data points for dependent variable +c npts - number of pairs of data points +c nterms - number of terms in fitting polynomial +c +c subroutines and function subprograms required +c integ (x, y, nterms, i1, x1, x2 sum) +c fits a polynomial with nterms starting at i1 +c and integrates area from x1 to x2 +c + function area (x,y,npts,nterms) + double precision sum + dimension x(1),y(1) +11 sum=0. + if (npts-nterms) 21,21,13 +13 neven=2*(nterms/2) + idelta=nterms/2-1 + if (nterms-neven) 31,31,51 +c +c fit all points with one curve +c +21 x1=x(1) + x2=x(npts) +23 call integ (x,y,npts,1,x1,x2,sum) + goto 71 +c +c even number of terms +c +31 x1=x(1) + j=nterms-idelta + x2=x(j) + call integ (x,y,nterms,1,x1,x2,sum) + i1=npts-nterms+1 + j=i1+idelta + x1=x(j) + x2=x(npts) +39 call integ (x,y,nterms,i1,x1,x2,sum) + if (i1-2) 71,71,41 +41 imax=i1-1 + do 46 i=2,imax + j=i+idelta + x1=x(j) + x2=x(j+1) +46 call integ (x,y,nterms,i,x1,x2,sum) + goto 71 +c +c odd number of terms +c +51 x1=x(1) + j=nterms-idelta + x2=(x(j)+x(j-1))/2. + call integ (x,y,nterms,1,x1,x2,sum) + i1=npts-nterms+1 + j=i1+idelta + x1=(x(j)+x(j+1))/2. + x2=x(npts) +59 call integ (x,y,nterms,i1,x1,x2,sum) + if (i1-2) 71,71,61 +61 imax=i1-1 + do 66 i=2,imax + j=i+idelta + x1=(x(j+1)+x(j))/2. + x2=(x(j+2)+x(j+1))/2. +66 call integ (x,y,nterms,i,x1,x2,sum) +71 area=sum + return + end |