implicit none

    integer::i
    integer, parameter::n = 13
    real*8::sxin, sxin2, sy1n, sy1n2, sxiny1n, DA, a0, a1, DC, c0, c1, seps,seps2
    real*8::samplemean
    real, dimension(n)::inputxi, rawoutputy1, x1, eps
    data inputxi/80.000, 82.500,  85.000,  87.500,  90.000,  92.500, 95.000,&
    97.500, 100.000, 102.500, 105.000, 107.500, 110.000/
    data rawoutputy1/3.0033, 3.4949, 3.9882, 4.4937, 5.0049, 5.5004, 6.0055,&
    6.4839, 6.9801, 7.5008, 7.9971, 8.5061, 9.0029/


     !########## STEP1  ##########!

    sxin = 0
    sxin2 = 0
    sy1n = 0
    sy1n2 = 0
    sxiny1n = 0

    do i = 1,n
        sxin = sxin + inputxi(i)
        sxin2 = sxin2 + inputxi(i) * inputxi(i)
        sy1n = sy1n + rawoutputy1(i)
        sy1n2 = sy1n2 + rawoutputy1(i) * rawoutputy1(i)
        sxiny1n = sxiny1n + inputxi(i) * rawoutputy1(i)
    end do

    print*, 'sxin = ',   sxin
    print*, 'sxin2 = ',  sxin2
    print*, 'sy1n = ',   sy1n
    print*, 'sy1n2 = ',  sy1n2
    print*, 'sxiny1n = ',sxiny1n



    !########## STEP2 ###########!
    DA = n * sxin2 - sxin * sxin
    print*, 'DA = ', DA

    a0 = (sy1n * sxin2 - sxin * sxiny1n) / DA
    print*, 'a0 = ', a0

    a1 = (n * sxiny1n - sxin * sy1n) / DA
    print*, 'a1 = ', a1

    print*
    print*, "@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@"
    print*

    !########### STEP3 #########!

    DC = n * sy1n2 - sy1n * sy1n
    print*, 'DC = ', DC

    c0 = (sxin * sy1n2 - sy1n * sxiny1n)/DC
    print*, 'c0 = ', c0

    c1 = (n * sxiny1n - sxin * sy1n)/DC
    print*, 'c1 = ', c1


    !########## STEP4 ##########!


    open(1, file = "table3-1.dat")

    do i = 1,n
        x1(i) =  c0 + c1 * rawoutputy1(i)

        eps(i) = (x1(i) - inputxi(i)) * 1000
        !mnozim sa 1000 jer je eps u tabeli u paskalima

        write(1, * ) inputxi(i), rawoutputy1(i), x1(i), eps(i)
    end do

    !########### STEP5 #########!
    seps = 0
    seps2 = 0

    do i = 1, n
        seps = seps + eps(i) / 1000

        seps2 = seps2 + (eps(i) * eps(i)) / 1000000
        print*, 'seps2 = ', seps2, 'eps = ', eps(i) / 1000000

    end do

    !ponovo je ovo u kilopaskalima

    print*, 'seps = ', seps
    print*, 'seps2 = ', seps2


    !########### STEP6 ###########!

    samplemean = seps/n
    print*, 'samplemean = ', samplemean
    print*



end program