Fortran/Przykłady/Metoda najmniejszych kwadratów

Poniższy program w FORTRAN-77 przybliża punkty pomiarowe prostą o równaniu y = ax + b wyliczając wartości współczynników a i b oraz ich niepewności.

     PROGRAM metoda_najmniejszych_kwadratow
     
     REAL x(20), y(20), a, b, xSredn, D, 
    1      yDelta, aDelta, bDelta
     INTEGER n
 
     PRINT *, 
     PRINT *, 
     PRINT *, 'Program dopasowuje prosta do zbioru punktow'
     PRINT *, 'doswiadczalnych wykorzystujac metode najmniejszych'
     PRINT *, 'kwadratow, ponadto oblicza odchylenia standardowe'
     PRINT *, 'obydwu parametrow prostej.'
     PRINT *, 
 
     PRINT *, 'Wprowadz liczbe punktow dosw. (min=3,maks=20)..'
     READ *, n
     IF ( (n.GT.20) .OR. (n.LT.3) ) THEN
       PRINT *, 'Liczba spoza zakresu !!'
       STOP
     ENDIF
 
     PRINT *, 'Wprowadz wspolrzedne punktow (na przemian: odcieta'
     PRINT *, 'i rzedna oddzielone spacja lub ENTER bez nawiasow'
     PRINT *, 'i przecinkow, wartosci calkowite od dziesietnych'
     PRINT *, 'oddzielone kropka)...'
     READ *, (x(i), y(i), i=1, n)      
 
     xSredn=0
     DO i=1, n
       xSredn = xSredn + x(i)
     ENDDO
     xSredn = xSredn / n
     
     D=0
     DO i=1, n
       D = D + ( x(i) - xSredn )**2  
     ENDDO
     IF (D.EQ.0) THEN
       PRINT *, 'Wprowadzone punkty sugeruja, ze szukana prosta'
       PRINT *, 'bedzie rownolegla do osi Oy, co uniemozliwia'
       PRINT *, 'obliczenie wartosci wspolczynnikow a i b rownania.'
       STOP
     ENDIF
  
     a=0
     DO i=1, n
       a = a + y(i) * ( x(i) - xSredn )
     ENDDO
     a = a / D
 
     b=0
     DO i=1, n
       b = b + ( y(i) - a*x(i) )
     ENDDO
     b = b / n
 
     yDelta = 0
     DO i=1, n
       yDelta = yDelta + ( y(i) - a*x(i) - b )**2
     ENDDO
     yDelta = yDelta / (n-2)
     yDelta = SQRT(yDelta)
 
     aDelta = yDelta / SQRT( D )
     bDelta = yDelta * SQRT( 1/n + (xSredn**2)/D )
 
     PRINT *, 'Wprowadzone punkty najlepiej przybliza prosta:' 
     PRINT *, '     y =', a, ' x +', b, '  ,'
     PRINT *, 'przy czym bledy wspolczynnikow a i b maja wartosc odp.:'
     PRINT *, '     ', aDelta, '  ,  ', bDelta
 
 
     END