! questa e' la finestra dell'editor
program main_mmq
implicit none
! dichiarazioni
integer, parameter		:: n_max = 100;
real, dimension(n_max)	:: x,y
integer					:: i,n
real					:: a11,a12,a21,a22,b1,b2
real                    :: d, c1, c2
real, parameter			:: g = 9.80665E+02 ! acc. gravita' cgs [cm/s^2]
real                    :: s
!...

! inizio del programma
!a) lettura file di dati
!
write(*,*) 'Lettura dati'
open(11,file='DATI.DAT',status='OLD', err=100)
! istruzioni di lettura
i=1
do while(.not. EOF(11))
	read(11,*, err=110) x(i),y(i)
	write(*,*) 'riga=',i,'        ',x(i), y(i)
	i = i + 1
end do
n = i - 1
close(11)

!b) calcoli di mmq
! Costruzione della matrice di coefficienti e del vettore noto
! del sistema di eq. normale 
a11=0.
a12=0.
a21=0.
a22=0.
b1 = 0.
b2 = 0.
do i=1,n
	a11 = a11 + 1.
	a12 = a12 + x(i)
	a22 = a22 + x(i)*x(i)
    b1 = b1 + y(i)
	b2 = b2 + y(i) * x(i)
end do
a21 = a12  ! per simmetria
write(*,*)
write(*,*) ' Sistema da risolvere'
write(*,*) a11, a12, '   c1    ',b1
write(*,*) a21, a22, '   c2  = ',b2
! Il sistema si risolve direttamente, come nelle dispense
d = a11*a22-a21*a12
write(*,*) 'Determinante = ',d
c1 = (b1*a22-a12*b2)/d
c2 = (a11*b2-b1*a21)/d
! 
write(*,*)
write(*,*) 'Soluzione'
write(*,*) 'c1=',c1, ' [cm]'
write(*,*) 'c2=',c2, ' [cm/g]'
write(*,*)
write(*,*) 'Costante della molla=',1.e-3*g/c2, '[N/m]'
write(*,*)
write(*,'(1x,t5,a,t20,a,t35,a)') 'm [g]','L_exp','L_calc'
s=0.
do i=1,n
	write(*,*) x(i), y(i), c1+c2*x(i)
	s = s + (c1+c2*x(i)-y(i))**2.
enddo
write(*,*)
write(*,*) 'S =',s
write(*,*) "S'=",s/(n-2)
stop
! Errori
100 write(*,*) 'Manca file DATI.DAT'
stop
110 write(*,*) "Errore di lettura in file 'DATI.DAT', riga ",i
stop

end program
