From fa080de7afc95aa1c19a6e6fc0e0708ced2eadc4 Mon Sep 17 00:00:00 2001 From: Joseph Hunkeler Date: Wed, 8 Jul 2015 20:46:52 -0400 Subject: Initial commit --- math/bevington/integ.f | 58 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 58 insertions(+) create mode 100644 math/bevington/integ.f (limited to 'math/bevington/integ.f') diff --git a/math/bevington/integ.f b/math/bevington/integ.f new file mode 100644 index 00000000..153bafe1 --- /dev/null +++ b/math/bevington/integ.f @@ -0,0 +1,58 @@ +c subroutine integ.f +c +c source +c Bevington, page 274. +c +c purpose +c integrate the area beneath two data points +c +c usage +c call integ (x, y, nterms, i1, x1, x2, sum) +c +c description of parameters +c x - array of data points for independent variable +c y - array of data points for dependent variable +c nterms - number of terms in fitting polymonial +c i1 - first data point for fitting polynomial +c x1 - first value of x for integration +c x2 - final value of x for integration +c +c subroutines and function subprograms required +c none +c +c comments +c dimension statement valid for nterms up to 10 +c + subroutine integ (x,y,nterms,i1,x1,x2,sum) + double precision xjk,array,a,denom,deltax,sum + dimension x(1),y(1) + dimension array(10,10) +c +c construct square matrix and invert +c +11 do 17 j=1,nterms + i=j+i1-1 + deltax=x(i)-x(i1) + xjk=1. + do 17 k=1,nterms + array(j,k)=xjk +17 xjk=xjk*deltax +21 call matinv (array,nterms,det) + if (det) 31,23,31 +23 imid=i1+nterms/2 + sum=sum+y(imid)*(x2-x1) + goto 40 +c +c evaluate coefficients and integrate +c +31 dx1=x1-x(i1) + dx2=x2-x(i1) +33 do 39 j=1,nterms + i=j+i1-1 + a=0. + do 37 k=1,nterms +37 a=a+y(i)*array(j,k) + denom=j +39 sum=sum+(a/denom)*(dx2**j-dx1**j) +40 return + end -- cgit