PROGRAM Quintic_Spline ! Fit a quintic spline with user control of knot positions. ! If the knots are at tk1, tk2, ..., then the fitted spline is ! b0 + b1.t + b2.t^2 + b3.t^3 + b4.t^4 + b5.t^5 for t <= tk1 ! b0 + ... + b5.t^5 + b6.(t-tk1)^5 for tk1 < t <= tk2 ! b0 + ... + b5.t^5 + b6.(t-tk1)^5 + b7.(t-tk2)^5 for tk2 < t <= tk3 ! b0 + ... + b5.t^5 + b6.(t-tk1)^5 + b7.(t-tk2)^5 + b8.(t-tk3)^5 ! for tk3 < t <= tk4, etc. ! In this version, the knots are evenly spaced. ! Also calculates first & 2nd derivatives of the spline. ! Uses the author's least-squares package in file lsq.f90 ! Latest revision - 2 November 2003 ! Alan Miller (amiller @ bigpond.net.au) USE lsq IMPLICIT NONE CHARACTER (LEN=50) :: infile, outfile INTEGER :: i, ier, iostatus, j, n, next_knot, nk, pos REAL (dp) :: t, t1, t2, y, d1, d2, dist, fitted REAL (dp), PARAMETER :: one = 1.0_dp REAL (dp), ALLOCATABLE :: knot(:), xrow(:), b(:) 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) ! Display first & last cases in data file n = 0 DO READ(UNIT=8, *, IOSTAT=iostatus) t, y IF (iostatus < 0) EXIT IF (iostatus > 0) CYCLE IF (n == 0) THEN WRITE(*, '(a, 2g13.4)') ' First t, y = ', t, y t1 = t END IF n = n + 1 END DO t2 = t WRITE(*, '(a, 2g13.4)') ' Last t, y = ', t, y WRITE(*, '(a, i6)') ' Number of cases read = ', n REWIND (UNIT=8) DO WRITE(*, '(a)', ADVANCE='NO') ' How many knots do you want? ' READ(*, *) nk IF (nk > n/5) THEN WRITE(*, *) ' ** Too many knots requested - TRY AGAIN' CYCLE ELSE IF (nk < 0) THEN WRITE(*, *) ' ** Must be >= 0 - TRY AGAIN' CYCLE ELSE EXIT END IF END DO ! Calculate knot positions, evenly spaced. ALLOCATE ( knot(nk) ) dist = (t2 - t1) / (nk + 1) DO i = 1, nk knot(i) = t1 + dist * i END DO WRITE(UNIT=9, '(a, i4)') 'Number of knots = ', nk ! Read in the data & update the least-squares calculations. ALLOCATE ( xrow(0:5+nk), b(0:5+nk) ) next_knot = 1 ! Initialize the least-squares calculations CALL startup(6+nk, .FALSE.) DO READ(UNIT=8, *, IOSTAT=iostatus) t, y IF (iostatus < 0) EXIT IF (iostatus > 0) CYCLE xrow(0) = one xrow(1) = (t - t1) xrow(2) = (t - t1) * xrow(1) xrow(3) = (t - t1) * xrow(2) xrow(4) = (t - t1) * xrow(3) xrow(5) = (t - t1) * xrow(4) IF (t > knot(next_knot)) THEN next_knot = MIN(nk, next_knot + 1) END IF DO j = 1, next_knot-1 xrow(5+j) = (t - knot(j))**5 END DO xrow(5+next_knot:5+nk) = 0.0_dp CALL includ(one, xrow, y) END DO CALL regcf(b, 6+nk, ier) WRITE(*, *) ' Coefficient Value' WRITE(*, '(a, g13.5)') ' Constant ', b(0) WRITE(*, '(a, g13.5)') ' Linear ', b(1) WRITE(*, '(a, g13.5)') ' Quadratic ', b(2) WRITE(*, '(a, g13.5)') ' Cubic ', b(3) WRITE(*, '(a, g13.5)') ' Quartic ', b(4) WRITE(*, '(a, g13.5)') ' Quintic ', b(5) WRITE(*, *) ' Knot position Quintic Coefficient' DO j = 1, nk WRITE(*, '(g13.5, t17, g13.5)') knot(j), b(5+j) END DO ! Calculate fitted values and derivatives next_knot = 1 REWIND (UNIT=8) WRITE(UNIT=9, '(a)') ' Time 1st deriv. 2nd deriv. Fitted value Actual' DO i = 1, n READ(UNIT=8, *) t, y xrow(0) = one xrow(1) = (t - t1) xrow(2) = (t - t1) * xrow(1) xrow(3) = (t - t1) * xrow(2) xrow(4) = (t - t1) * xrow(3) xrow(5) = (t - t1) * xrow(4) IF (t > knot(next_knot)) THEN next_knot = next_knot + 1 IF (next_knot <= nk) THEN WRITE(UNIT=9, '(a, g13.5)') 'New knot at t = ', knot(next_knot-1) ELSE next_knot = nk + 1 END IF END IF DO j = 1, next_knot-1 xrow(5+j) = (t - knot(j))**5 END DO fitted = DOT_PRODUCT( b(0:5+next_knot-1), xrow(0:5+next_knot-1) ) d2 = ((20*b(5)*(t-t1) + 12*b(4))*(t-t1) + 6*b(3))*(t-t1) + 2*b(2) d1 = (((5*b(5)*(t-t1) + 4*b(4))*(t-t1) + 3*b(3))*(t-t1) + 2*b(2))*(t-t1) + b(1) DO j = 1, next_knot-1 d1 = d1 + 5*b(j+5)*(t - knot(j))**4 d2 = d2 + 20*b(j+5)*(t - knot(j))**3 END DO WRITE(UNIT=9, '(f8.3, 4g13.4)') t, d1, d2, fitted, y END DO STOP END PROGRAM Quintic_Spline