PROGRAM Weighted_Quintic ! Fit a quintic polynomial using exponential weighting. ! Full weight = 1.0 is given to the latest case, weight = lambda ! is given to the previous case, lambda^2 to the previous case, etc. ! Also calculates first & 2nd derivatives of the polynomial. ! Uses the author's least-squares package in file lsq.f90 ! Latest revision - 11 November 2001 ! Alan Miller (amiller @ bigpond.net.au) USE lsq IMPLICIT NONE CHARACTER (LEN=50) :: infile, outfile CHARACTER (LEN=10) :: text INTEGER :: ier, iostatus, pos REAL (dp) :: lambda, xrow(0:5), t, t0, y, b(0:5), d1, d2, fitted REAL (dp), PARAMETER :: one = 1.0_dp DO WRITE(*, '(a)', ADVANCE='NO') ' Enter file name: ' READ(*, *) infile OPEN(UNIT=8, FILE=infile, STATUS='OLD', IOSTAT=iostatus) IF (iostatus == 0) EXIT WRITE(*, '(a, a, a)') ' *** File: ', TRIM(infile), ' not found -- TRY AGAIN!' END DO ! Name of output file = name of input file with extension .OUT pos = INDEX(infile, '.') IF (pos > 0) THEN outfile = infile(:pos) // 'out' ELSE outfile = TRIM(infile) // '.out' END IF OPEN(UNIT=9, FILE=outfile) ! Enter the exponential weight DO WRITE(*, '(a)', ADVANCE='NO') ' Enter exponential weight (default = 0.95): ' READ(*, '(a)') text IF (LEN_TRIM(text) == 0) THEN lambda = 0.95_dp EXIT ELSE READ(text, '(f9.3)', IOSTAT=iostatus) lambda IF (iostatus /= 0) CYCLE IF (lambda > 0.0_dp .AND. lambda < one) EXIT WRITE(*, *) ' Exponential weight must be between 0 and 1' END IF END DO WRITE(UNIT=9, '(a, f9.4)') 'Exponential weight = ', lambda ! Initialize the least-squares calculations CALL startup(6, .FALSE.) ! Read in the data, one line at a time, and update the quintic WRITE(UNIT=9, '(a)') ' Time 1st deriv. 2nd deriv. Fitted value Actual' DO READ(UNIT=8, *, IOSTAT=iostatus) t, y IF (iostatus < 0) EXIT IF (iostatus > 0) CYCLE IF (nobs == 0) t0 = t xrow(0) = one xrow(1) = t - t0 xrow(2) = (t - t0) * xrow(1) xrow(3) = (t - t0) * xrow(2) xrow(4) = (t - t0) * xrow(3) xrow(5) = (t - t0) * xrow(4) CALL includ(one, xrow, y) IF (nobs > 5) THEN CALL regcf(b, 6, ier) d2 = ((20*b(5)*(t-t0) + 12*b(4))*(t-t0) + 6*b(3))*(t-t0) + 2*b(2) d1 = (((5*b(5)*(t-t0) + 4*b(4))*(t-t0) + 3*b(3))*(t-t0) + 2*b(2))*(t-t0) + b(1) fitted = ((((b(5)*(t-t0) + b(4))*(t-t0) + b(3))*(t-t0) + b(2))*(t-t0) + b(1))*(t-t0) + b(0) WRITE(UNIT=9, '(f8.3, 4g13.4)') t, d1, d2, fitted, y END IF ! Now downweight the calculations up to this point d = lambda * d END DO STOP END PROGRAM Weighted_Quintic