!**************************************************** !* Program to demonstrate the general integration * !* subroutine. Example is the integral of SIN(X) * !* from X1 to X2 in real*8 precision. * !* ------------------------------------------------ * !* Reference: BASIC Scientific Subroutines, Vol. II * !* By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1].* !* * !* F90 version by J-P Moreau, Paris * !* (www.jpmoreau.fr) * !* ------------------------------------------------ * !* SAMPLE RUN: * !* * !* Integration of SIN(X) * !* * !* Input start point end point: 0 0.5 * !* Integral from 0.000 to 0.500 equals 0.1224231 * !* Error check: 1 * !* * !**************************************************** PROGRAM Integra parameter(SIZE=14) real*8 X(0:SIZE),Y(0:SIZE) real*8 x1,x2,zz integer z1 iv=14; ! Number of pooints in table print *,' ' print *,'Integration of SIN(X)' print *,' ' ! Input sine table !----------------------------------------------------------------- ! Sine table values from Handbook of mathematical functions ! by M. Abramowitz and I.A. Stegun, NBS, june 1964 !----------------------------------------------------------------- X(1)=0.000; Y(1)=0.00000000; X(2)=0.125; Y(2)=0.12467473; X(3)=0.217; Y(3)=0.21530095; X(4)=0.299; Y(4)=0.29456472; X(5)=0.376; Y(5)=0.36720285; X(6)=0.450; Y(6)=0.43496553; X(7)=0.520; Y(7)=0.49688014; X(8)=0.589; Y(8)=0.55552980; X(9)=0.656; Y(9)=0.60995199; X(10)=0.721; Y(10)=0.66013615; X(11)=0.7853981634; Y(11)=0.7071067812; X(12)=0.849; Y(12)=0.75062005; X(13)=0.911; Y(13)=0.79011709; X(14)=0.972; Y(14)=0.82601466; !---------------------------------------------------------------- write(*,"(' Input start point end point: ')",advance='no') read *, x1,x2 call Integrale(iv,x1,x2,X,Y,zz,z1) ! Call integral function write(*,50) x1,x2,zz write(*,60) z1 stop 50 format(/' Integral from ',f6.3,' to ',f6.3,' equals ',f10.7) 60 format(' Error check: ',i1/) end !******************************************************** !* Akima spline fitting subroutine * !* ---------------------------------------------------- * !* The input table is X(i), Y(i), where Y(i) is the * !* dependant variable. The interpolation point is xx, * !* which is assumed to be in the interval of the table * !* with at least one table value to the left, and three * !* to the right. The interpolated returned value is yy. * !* n is returned as an error check (n=0 implies error). * !* It is also assumed that the X(i) are in ascending * !* order. * !******************************************************** Subroutine Interpol_Akima(iv,n,xx,yy,X,Y) !Labels: 100,200,300 parameter(SIZE=14) integer i,iv,n real*8 a,b,xx,yy real*8 X (0:SIZE), Y (0:SIZE) real*8 XM(0:SIZE+3) real*8 Z (0:SIZE) n=1 !special case xx=0 if (xx.eq.0.0) then yy=0.d0; return end if !Check to see if interpolation point is correct if (xx=X(iv-3)) then n=0 ; return end if X(0)=2.d0*X(1)-X(2) !Calculate Akima coefficients, a and b do i = 1, iv-1 !Shift i to i+2 XM(i+2)=(Y(i+1)-Y(i))/(X(i+1)-X(i)) end do XM(iv+2)=2.d0*XM(iv+1)-XM(iv) XM(iv+3)=2.d0*XM(iv+2)-XM(iv+1) XM(2)=2.d0*XM(3)-XM(4) XM(1)=2.d0*XM(2)-XM(3) do i = 1, iv a=dabs(XM(i+3)-XM(i+2)) b=dabs(XM(i+1)-XM(i)) if (a+b.ne.0.d0) goto 100 Z(i)=(XM(i+2)+XM(i+1))/2.d0 goto 200 100 Z(i)=(a*XM(i+1)+b*XM(i+2))/(a+b) 200 end do !Find relevant table interval i=0 300 i=i+1 if (xx>X(i)) goto 300 i=i-1 !Begin interpolation b=X(i+1)-X(i) a=xx-X(i) yy=Y(i)+Z(i)*a+(3.d0*XM(i+2)-2.d0*Z(i)-Z(i+1))*a*a/b yy=yy+(Z(i)+Z(i+1)-2.d0*XM(i+2))*a*a*a/(b*b) return end !****************************************************** !* Routine for the first trapezoidal integration, xi1 * !* n1 keeps track of the number of intervals. * !****************************************************** Subroutine Trapez1(x1,x2,i,iv,X,Y,xi1) ! Labels: 100,150,200 parameter(SIZE=14) real*8 X(0:SIZE),Y(0:SIZE) real*8 d,x1,x2,xi1,xx,yy integer i,iv,n,j1 !n not used here xi1=0 ; n1 = 0 ; xx=x1 ! Call interpolation subroutine call Interpol_Akima(iv,n,xx,yy,X,Y) ! Is there at least one table interval? if (x2>X(i+1)) goto 100 !If not, integral is simple n1=n1+1 ; d=yy ; xx=x2 ! Find end point yy value call Interpol_Akima(iv,n,xx,yy,X,Y) xi1=(yy+d)*(x2-x1)/2.d0 return ! At least one table interval must be summed over 100 j1=i xi1=xi1+(yy+Y(i+1))*(X(i+1)-xx)/2.d0 ! Any more intervals? If not, finish integral with end point 150 if (x2X(i+1)) goto 600 xx=x1+(x2-x1)/2.d0 call Interpol_Akima(iv,n,xx,yy,X,Y) xi2=xi2+(d+yy)*(x2-x1)/4.d0 return 600 xx=x1+(X(i+1)-x1)/2.d0 j1=i call Interpol_Akima(iv,n,xx,yy,X,Y) xi2=xi2+(yy+d)*(xx-x1)/2.d0 xi2=xi2+(yy+Y(j1+1))*(X(j1+1)-xx)/2.d0 650 if (x2error). * !******************************************************** Subroutine Integrale(iv,x1,x2,X,Y,zz,z1) ! Labels: 100, 200 parameter(SIZE=14) real*8 X(0:SIZE),Y(0:SIZE) real*8 xi1,xi2,x1,x2,x3,zz integer i,iv,z1 zz=0 ; z1=0 ! Check to see if end points are correct if (x1X(iv-3)) return ! If x1>x2 then switch and set flag if (x1