program leastsq implicit none c c Program designed to perform a least squares fit of a quadratic equation c to some data. In this specific example the data is the result of c measuring the location of a falling object at various times. The result c of interest is gravitational acceleration obtained from the coefficient c off t**2. c real t(50),z(50),coef(0:2),g integer np,npoly c c Variables c c t - array containing times (seconds) at which positions are measured c z - measured distance of the falling object from the nominal point of c release (fit may let you estimate any offset in the actual release c point c np - number of data points c coef - Coefficients in the second order polynomial approximating the data c g - approximation to the gravitational acceleration deduced from data c c npoly=2 call getdata(t,z,np) call fit(t,z,np,npoly,coef(0)) g=2*coef(2) write(6,2000) g 2000 format (' Predicted value of g = ', f6.3,'m/s**2') write(6,2001) coef(0),coef(1) 2001 format(' Predicted initial offset = ',f7.4,'m/s'/ # ' Predicted initial velocity = ',f7.4,'m') stop end subroutine getdata(t,s,n ) implicit none c c Get Experimental Data c c Input Arguments: NONE c Output Arguments: c t - time (sec) of measurement c s - distance of fall measured (meters) c n - number of data points c integer n real t(*),s(*) open(11,file='fall.data') n=0 10 read (11,*,end=20 ) t(n+1),s(n+1) n=n+1 go to 10 20 return end subroutine fit (xd,yd,ndata,npoly,c) implicit none c c Take corresponding data points from the arrays xd and yd and fit c them with the following equation: c c y = c(1) + c(2) * x + c(3) * x**2 + ... + c(npoly+1)*x**npoly c real xd(ndata),yd(ndata),aa(npoly+1,npoly+1),c(*),fs, & sqrt,rsid,xdsum(0:(npoly+npoly)),xdp integer ipvt(npoly+1),ir,is,ij,ndata, info, npoly,i,j,npow, neq c c NOTE: The use of aa(npoly+1,npoly+1), and ipvt(npoly+1) above is c an example of the FORTRAN 90 AUTOMATIC ARRAY feature. It is c a short-hand way of using the ALLOCATE statement. Space for c aa, c, and ipvt is assigned when SUBROUTINE FIT is entered and c is DEALLOCATED automatically when a RETURN is executed. c Note the flexibility available in assigning the array sizes. c A compiler without FORTRAN 90 will not permit using npoly to c specify the sizes of these arrays unless they (aa, c, ipvt) c also appear in the argument list. Beware of use of this c convenient FORTRAN 90 feature in frequently used subroutines. c A fair amount of machine time is associated with each ALLOCATE c and DEALLOCATE. c c Input Arguments: c xd - x values for data points c yd - y valuess for data point pairs c ndata - number of data points c npoly - order of polynomial fit (highest power of x in polynomial) c c Output Arguements: c c - array containing the coefficients in the chosen order c polynomial that provide the "best" fit to the data from a c least squares method. Note that it also temporarily holds c the values of the right hand sides of each Least Squares c equation. c Other Key Variables: c xdsum - array containing the sums of all needed powers of the elements c in xd c aa - matrix containing coefficients of the system of equations c generated by the least squares method c ipvt - an array used by Linpack routines to keep track of "pivoting" c operations during Gauss Elimination. c neq - number of linear equations to be solved c c Zero Sum counters c npow=npoly+npoly neq=npoly+1 xdsum(0)=ndata do 20 j=1,npow 20 xdsum(j)=0. do 21 j=1,neq 21 c(j)=0. c c Sum the Powers of xd, note that I duck having to use exponentiation (**) c do 40 i=1,ndata xdp=1. do 40 j=1,npow xdp=xdp*xd(i) 40 xdsum(j)=xdsum(j)+xdp c c DO loop 45 generates the right hand sides for all equations c do 45 i=1,ndata xdp=1. c(1)=c(1)+yd(i) do 45 j=2,neq xdp=xdp*xd(i) 45 c(j)=c(j)+yd(i)*xdp c c DO loop 50 generates the coefficients a given equation c do 50 i=1,neq do 50 j=1,neq 50 aa(i,j)=xdsum(i+j-2) c c Solve the Least Squares Equations c call sgefa(aa,neq,neq,ipvt,info) call sgesl(aa,neq,neq,ipvt,c ,0) c c Calculate a measure of mean error between the data and the curve c rsid=0. do 65 ir=1,ndata fs=0. do 60 is=1,neq 60 fs=fs+c(is)*xd(ir)**(is-1) 65 rsid=rsid+(yd(ir)-fs)**2 rsid=sqrt(rsid/float(ndata-1)) write(6,2001) rsid 2001 format(' Fit to 2nd order polynomial has a mean error of', $ 1p,e12.5) return end