{"id":211,"date":"2017-09-13T11:14:49","date_gmt":"2017-09-13T14:14:49","guid":{"rendered":"http:\/\/www.professores.uff.br\/diomarcesarlobao\/?page_id=211"},"modified":"2017-09-13T11:14:49","modified_gmt":"2017-09-13T14:14:49","slug":"lsq2-f","status":"publish","type":"page","link":"https:\/\/www.professores.uff.br\/diomarcesarlobao\/lsq2-f\/","title":{"rendered":"lsq2.f"},"content":{"rendered":"<pre>      program leastsq\r\n      implicit none\r\nc\r\nc   Program designed to perform a least squares fit of a quadratic equation\r\nc   to some data.  In this specific example the data is the result of\r\nc   measuring the location of a falling object at various times.  The result\r\nc   of interest is gravitational acceleration obtained from the coefficient\r\nc   off t**2.\r\nc\r\n      real t(50),z(50),coef(0:2),g\r\n      integer np,npoly\r\nc\r\nc   Variables\r\nc\r\nc   t   -  array containing times (seconds) at which positions are measured\r\nc   z   -  measured distance of the falling object from the nominal point of\r\nc          release (fit may let you estimate any offset in the actual release\r\nc          point\r\nc   np  -  number of data points\r\nc   coef - Coefficients in the second order polynomial approximating the data\r\nc   g   -  approximation to the gravitational acceleration deduced from data\r\nc\r\nc\r\n      npoly=2\r\n      call getdata(t,z,np)\r\n      call fit(t,z,np,npoly,coef(0))\r\n      g=2*coef(2)\r\n      write(6,2000) g\r\n 2000 format (' Predicted value of g = ', f6.3,'m\/s**2')\r\n      write(6,2001) coef(0),coef(1)\r\n 2001 format(' Predicted initial offset = ',f7.4,'m\/s'\/\r\n     #         ' Predicted initial velocity = ',f7.4,'m')\r\n      stop\r\n      end\r\n      subroutine getdata(t,s,n )\r\n      implicit none\r\nc\r\nc   Get Experimental Data\r\nc\r\nc   Input Arguments: NONE\r\nc   Output Arguments:\r\nc     t     -    time (sec) of measurement\r\nc     s     -    distance of fall measured (meters)\r\nc     n     -    number of data points\r\nc\r\n      integer n\r\n      real t(*),s(*)\r\n      open(11,file='fall.data')\r\n      n=0\r\n   10 read (11,*,end=20 ) t(n+1),s(n+1)\r\n      n=n+1\r\n      go to 10\r\n   20 return\r\n      end\r\n      subroutine fit (xd,yd,ndata,npoly,c)\r\n      implicit none\r\nc\r\nc       Take corresponding data points from the arrays xd and yd and fit\r\nc       them with the following equation:\r\nc\r\nc       y = c(1) + c(2) * x + c(3) * x**2  + ... + c(npoly+1)*x**npoly\r\nc\r\n      real xd(ndata),yd(ndata),aa(npoly+1,npoly+1),c(*),fs,\r\n     &amp; sqrt,rsid,xdsum(0:(npoly+npoly)),xdp\r\n      integer ipvt(npoly+1),ir,is,ij,ndata, info, npoly,i,j,npow, neq\r\nc\r\nc    NOTE:  The use of aa(npoly+1,npoly+1), and ipvt(npoly+1) above is\r\nc           an example of the FORTRAN 90 AUTOMATIC ARRAY feature.  It is\r\nc           a short-hand way of using the ALLOCATE statement.  Space for\r\nc           aa, c, and ipvt is assigned when SUBROUTINE FIT is entered and\r\nc           is DEALLOCATED automatically when a RETURN is executed.\r\nc           Note the flexibility available in assigning the array sizes.\r\nc           A compiler without FORTRAN 90 will not permit using npoly to\r\nc           specify the sizes of these arrays unless they (aa, c, ipvt)\r\nc           also appear in the argument list.  Beware of use of this\r\nc           convenient FORTRAN 90 feature in frequently used subroutines.\r\nc           A fair amount of machine time is associated with each ALLOCATE\r\nc           and DEALLOCATE.\r\nc\r\nc    Input Arguments:\r\nc      xd   -   x values for data points\r\nc      yd   -   y valuess for data point pairs\r\nc      ndata -  number of data points\r\nc      npoly -  order of polynomial fit (highest power of x in polynomial)\r\nc\r\nc    Output Arguements:\r\nc      c    -   array containing the coefficients in the chosen order\r\nc               polynomial that provide the \"best\" fit to the data from a\r\nc               least squares method.  Note that it also temporarily holds\r\nc               the values of the right hand sides of each Least Squares\r\nc               equation.\r\nc    Other Key Variables:\r\nc      xdsum -  array containing the sums of all needed powers of the elements\r\nc               in xd\r\nc      aa   -   matrix containing coefficients of the system of equations\r\nc               generated by the least squares method\r\nc      ipvt -   an array used by Linpack routines to keep track of \"pivoting\"\r\nc               operations during Gauss Elimination.\r\nc      neq   -  number of linear equations to be solved\r\nc\r\nc      Zero Sum counters\r\nc\r\n      npow=npoly+npoly\r\n      neq=npoly+1\r\n      xdsum(0)=ndata\r\n      do 20 j=1,npow\r\n   20   xdsum(j)=0.\r\n      do 21 j=1,neq\r\n   21   c(j)=0.\r\nc\r\nc   Sum the Powers of xd, note that I duck having to use exponentiation (**)\r\nc\r\n      do 40 i=1,ndata\r\n         xdp=1.\r\n         do 40 j=1,npow\r\n      \t    xdp=xdp*xd(i)\r\n   40       xdsum(j)=xdsum(j)+xdp\r\nc\r\nc       DO loop 45 generates the right hand sides for all equations\r\nc\r\n      do 45 i=1,ndata\r\n         xdp=1.\r\n         c(1)=c(1)+yd(i)\r\n         do 45 j=2,neq\r\n            xdp=xdp*xd(i)\r\n  45        c(j)=c(j)+yd(i)*xdp\r\nc\r\nc       DO loop 50 generates the coefficients a given equation\r\nc\r\n      do 50 i=1,neq\r\n         do 50 j=1,neq\r\n  50        aa(i,j)=xdsum(i+j-2)\r\nc\r\nc    Solve the Least Squares Equations\r\nc\r\n      call sgefa(aa,neq,neq,ipvt,info)\r\n      call sgesl(aa,neq,neq,ipvt,c ,0)\r\nc\r\nc    Calculate a measure of mean error between the data and the curve\r\nc\r\n            rsid=0.\r\n            do 65 ir=1,ndata\r\n               fs=0.\r\n               do 60 is=1,neq\r\n   60             fs=fs+c(is)*xd(ir)**(is-1)\r\n   65          rsid=rsid+(yd(ir)-fs)**2\r\n            rsid=sqrt(rsid\/float(ndata-1))\r\n            write(6,2001) rsid\r\n 2001       format(' Fit to 2nd order polynomial has a mean error of',\r\n     $             1p,e12.5)\r\n      return\r\n      end<\/pre>\n","protected":false},"excerpt":{"rendered":"<p>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 [&hellip;]<\/p>\n","protected":false},"author":22,"featured_media":0,"parent":0,"menu_order":0,"comment_status":"closed","ping_status":"closed","template":"","meta":{"_exactmetrics_skip_tracking":false,"_exactmetrics_sitenote_active":false,"_exactmetrics_sitenote_note":"","_exactmetrics_sitenote_category":0,"footnotes":""},"categories":[],"tags":[],"class_list":["post-211","page","type-page","status-publish","hentry"],"_links":{"self":[{"href":"https:\/\/www.professores.uff.br\/diomarcesarlobao\/wp-json\/wp\/v2\/pages\/211","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/www.professores.uff.br\/diomarcesarlobao\/wp-json\/wp\/v2\/pages"}],"about":[{"href":"https:\/\/www.professores.uff.br\/diomarcesarlobao\/wp-json\/wp\/v2\/types\/page"}],"author":[{"embeddable":true,"href":"https:\/\/www.professores.uff.br\/diomarcesarlobao\/wp-json\/wp\/v2\/users\/22"}],"replies":[{"embeddable":true,"href":"https:\/\/www.professores.uff.br\/diomarcesarlobao\/wp-json\/wp\/v2\/comments?post=211"}],"version-history":[{"count":1,"href":"https:\/\/www.professores.uff.br\/diomarcesarlobao\/wp-json\/wp\/v2\/pages\/211\/revisions"}],"predecessor-version":[{"id":212,"href":"https:\/\/www.professores.uff.br\/diomarcesarlobao\/wp-json\/wp\/v2\/pages\/211\/revisions\/212"}],"wp:attachment":[{"href":"https:\/\/www.professores.uff.br\/diomarcesarlobao\/wp-json\/wp\/v2\/media?parent=211"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/www.professores.uff.br\/diomarcesarlobao\/wp-json\/wp\/v2\/categories?post=211"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/www.professores.uff.br\/diomarcesarlobao\/wp-json\/wp\/v2\/tags?post=211"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}