root / synthbench / euroben-dm / mod1k / lsq.f @ 0:839f52ef7657
History | View | Annotate | Download (2.9 kB)
1 |
Subroutine lsq ( k, x, y, slope, incpt, perc ) |
---|---|
2 |
! ---------------------------------------------------------------------- |
3 |
! --- Forms Least Squares fit of straight line to data received so |
4 |
! far. |
5 |
! --- Input data in pairs (x,y). Best straight line fitted is described |
6 |
! by the output parameters slope and incpt. |
7 |
! |
8 |
! --- k - If k = 0 initialize, else to fit line. |
9 |
! x - Abscissa. |
10 |
! y - Ordinate. |
11 |
! slope - Inverse slope of best line of all data up to now. |
12 |
! incpt - Negative intercept of best line on X-axis. |
13 |
! perc - Percentage of error. |
14 |
! ---------------------------------------------------------------------- |
15 |
Implicit None |
16 |
|
17 |
Integer :: k |
18 |
Real(8) :: x, y, slope ,incpt, perc |
19 |
|
20 |
Integer :: ndata |
21 |
Real(8), Parameter :: zero = 0.0d0, one = 1.0d0 |
22 |
Real(8) :: permax, permin |
23 |
Real(8) :: dx, dy, rmse, denom, arg |
24 |
Real(8) :: xndata, sumx, sumx2, sumy, sumy2, sumxy |
25 |
|
26 |
Save xndata, sumx, sumx2, sumy, sumy2, sumxy |
27 |
! ---------------------------------------------------------------------- |
28 |
! --- Initialise. |
29 |
|
30 |
If ( k == 0 ) Then |
31 |
xndata = zero |
32 |
sumx = zero |
33 |
sumx2 = zero |
34 |
sumy = zero |
35 |
sumy2 = zero |
36 |
sumxy = zero |
37 |
permin = Huge( one ) |
38 |
permax = zero |
39 |
perc = zero |
40 |
Return |
41 |
End If |
42 |
! ---------------------------------------------------------------------- |
43 |
! --- Update sums with new data x and y. |
44 |
|
45 |
xndata = xndata + one |
46 |
ndata = Int( xndata ) |
47 |
sumx = sumx + x |
48 |
sumx2 = sumx2 + x*x |
49 |
sumy = sumy + y |
50 |
sumy2 = sumy2 + y*y |
51 |
sumxy = sumxy + x*y |
52 |
If ( ndata < 3 ) Then |
53 |
slope = zero |
54 |
incpt = zero |
55 |
rmse = zero |
56 |
Else ! --- Calculate new slope and incpt. |
57 |
denom = xndata*sumx2 - sumx*sumx |
58 |
If ( Abs( denom ) < Epsilon( one ) ) Then |
59 |
slope = zero |
60 |
incpt = zero |
61 |
rmse = zero |
62 |
Else |
63 |
slope = ( xndata*sumxy - sumx*sumy )/denom |
64 |
incpt = ( sumy*sumx2 - sumxy*sumx )/denom |
65 |
arg = ( sumy2 + xndata*incpt**2 + sumx2*slope**2 - |
66 |
& 2.0d0*( incpt*sumy - incpt*slope*sumx + |
67 |
& slope*sumxy ) )/( xndata - 2.0d0 ) |
68 |
If ( arg >= 0 ) Then |
69 |
rmse = Sqrt( arg*sumx2/denom ) |
70 |
Else |
71 |
rmse = -Sqrt(-arg*sumx2/denom ) |
72 |
End If |
73 |
End If |
74 |
End If |
75 |
! ---------------------------------------------------------------------- |
76 |
! --- Update errors of percentages. |
77 |
|
78 |
If ( ndata > 2 ) Then |
79 |
perc = 100.0d0*rmse/y |
80 |
If ( ndata == 3 ) permin = perc |
81 |
permax = Max( permax, perc ) |
82 |
End If |
83 |
! --------------------------------------------------------------------- |
84 |
End Subroutine lsq |