module arrays real, allocatable :: t(:), z(:) integer np 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 end module c program leastsq use arrays 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 c John Mahaffy 2/6/95 c real coef(0:2),g 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 z = coef(0) + coef(1)*t + coef(2)*t**2 c g - approximation to the gravitational acceleration deduced from data c c call getdata call quadfit(t,z,np,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'/ # ' Predicted initial velocity = ',f7.4,' m/s') stop end subroutine getdata use arrays 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 np - number of data points c integer i, iend open(11,file='fall.data') c c Count Data Pairs in fall.data c np=-1 do while (iend.eq.0) read(11,*,iostat=iend) np = np + 1 enddo if (iend.gt.0.or. np.lt.1) then print *, 'Empty or Bad File, check fall.data' stop endif c c Rewind, Allocate Arrays, and Read data c rewind (11) print *, np, ' data points read.' allocate (t(1:np),z(1:np)) do i = 1,np read (11,*) t(i),z(i) enddo c end c subroutine quadfit (xd,yd,ndata,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 real xd(ndata),yd(ndata),aa(3,3),c(3),fs,suma,sumb,sqrt,rsid integer ipvt(3),ir,is,ij,ndata, info 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 c Output Arguements: c c - array containing the three coefficients in the 2nd 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 aa - matrix containing coefficients of the system of equations c generated by the least squares method c suma - a variable that tallies the sum that is needed to generate c each element of aa. c sumb - a variable that tallies the sum that is needed for the right c hand side of each Least Squares equation. c ir - an index that is used to keep track of the equation number c within the system of equations . c is - an index used to track the coefficient number within a given c Least Squares equation. c c c DO loop 55 generates terms for each of the 3 Least Squares Equations c do 55 ir=1,3 c c DO loop 45 generates the right hand side for a given equation c sumb=0. do 45 ij=1,ndata 45 sumb=sumb+yd(ij)*xd(ij)**(ir-1) c(ir)=sumb c c DO loop 50 generates the coefficients a given equation c do 50 is=1,3 suma=0. do 52 ij=1,ndata 52 suma=suma+xd(ij)**(is-1)*xd(ij)**(ir-1) 50 aa(ir,is)=suma 55 continue c c Solve the Least Squares Equations c call sgefa(aa,3,3,ipvt,info) call sgesl(aa,3,3,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,3 60 fs=fs+c(is)*xd(ir)**(is-1) 65 rsid=rsid+(yd(ir)-fs)**2 c rsid=sqrt(rsid/float(ndata-1)) c write(6,2001) rsid 2001 format(' Fit to 2nd order polynomial has a mean error of', $ 1p,e12.5) return end cc c