PROGRAM quad_surf ! Fit a quadratic surface to a set of data USE lsq IMPLICIT NONE INTEGER :: iostatus, n, ier, row, col INTEGER, PARAMETER :: nmax = 100 REAL (dp), PARAMETER :: wt = 1.0_dp REAL (dp) :: rh, temp, yval, x(nmax,5), y(nmax), xrow(0:5), b(0:5), & yfit(3) OPEN(UNIT=8, FILE='dry_bulb.dat', STATUS='OLD') n = 0 CALL startup(5, .TRUE.) DO READ(8, *, IOSTAT=iostatus) rh, temp, yval IF (iostatus < 0) EXIT IF (iostatus > 0) CYCLE n = n + 1 rh = rh - 70._dp temp = temp - 95._dp y(n) = yval x(n,1) = rh x(n,2) = temp x(n,3) = rh**2 x(n,4) = rh*temp x(n,5) = temp**2 xrow(1:5) = x(n,:) xrow(0) = 1.0_dp CALL includ(wt, xrow, y(n)) END DO CALL regcf(b, 6, ier) WRITE(*, *) 'Fitted quadratic surface:' WRITE(*, '(a, 2f8.3, a, f8.4, a, f8.4, a, f8.4, a, f8.4, a)') & ' Y = ', b(0), b(1), '*X1 ', b(2), '*X2 ', b(3), '*X1^2 ', & b(4), '*X1*X2 ', b(5), '*X2^2' WRITE(*, *) WRITE(*, *) 'Fitted values:' WRITE(*, *) ' %RH 90 95 100' WRITE(*, *) '-------------------------------------' DO row = 1, 4 xrow(0) = 1.0_dp xrow(1) = 20.0 - 10.0*row xrow(3) = xrow(1)**2 DO col = 1, 3 xrow(2) = -10.0 + 5.0*col xrow(4) = xrow(1)*xrow(2) xrow(5) = xrow(2)**2 yfit(col) = DOT_PRODUCT( b, xrow ) END DO WRITE(*, '(f8.0, " ", 3f9.1)') xrow(1) + 70., yfit END DO STOP END PROGRAM quad_surf