!     Last change:  VGC  11 Mar 2003    9:48 pm


MODULE SUBS
  CONTAINS


  !-----------------------------------------!
  ! Computes yhat from M, K, a and x vector !
  !-----------------------------------------!
  SUBROUTINE MODEL (M,K,A,X,YHAT)
    IMPLICIT NONE
    REAL, INTENT(IN):: M, K, A
    REAL, DIMENSION(1:19), INTENT(IN):: X
    REAL, DIMENSION(1:19), INTENT(OUT):: YHAT
    INTEGER:: I
    DO I = 1,19
       YHAT(I) = M/(1 + K*EXP(-A*X(I)))   ! Logistic model !
    END DO
  RETURN
  END SUBROUTINE MODEL


  !--------------------------------------!
  ! Computes square of distance between  !
  ! vectors y and yhat in 19-dim space.  !
  ! Square of distance called SSE.       !
  !--------------------------------------!
  FUNCTION SSE(Y, YHAT)
    IMPLICIT NONE
    INTEGER:: I
    REAL,INTENT(IN), DIMENSION(1:19):: Y, YHAT
    REAL:: SSE
    SSE = 0.0
    DO I = 1,19
      SSE = SSE + (Y(I) - YHAT(I))**2
    END DO
  RETURN
  END FUNCTION SSE


  !-------------------------------------------------------------!
  ! Makes error table for best logistic curve fit to empirical  !
  ! US population data.                                         !
  !-------------------------------------------------------------!
  SUBROUTINE TABLE (M, K, A, X, Y)
    IMPLICIT NONE
    REAL, INTENT(IN):: M, K, A
    REAL, DIMENSION(1:19), INTENT(IN):: X, Y
    REAL, DIMENSION(1:19):: YHAT
    REAL:: ERROR, RELERROR
    INTEGER:: I
    OPEN (UNIT = 2, FILE = 'MODTABLE.DAT', STATUS = 'UNKNOWN')
    CALL MODEL (M,K,A,X,YHAT)
    WRITE (2, '(A4,5X,5A15)') 'Year','            t  ','            y  ','         yhat  ', &
    '        Error ','      Rel Error'
    DO I = 1,19
    ERROR = YHAT(I) - Y(I)
    RELERROR = ERROR/Y(I)
    WRITE (2, '(I4,5X,5F15.2)')INT(X(I)+1790.5),X(I),Y(I),YHAT(I),ERROR,RELERROR
    END DO
  RETURN
  END SUBROUTINE TABLE

END MODULE SUBS


! --------------------------------------------------------!
! Finds parameters M, K and a for the best logistic fit   !
! to the US Population data by complete coverage method.  !
!                                                         !
! by J R Nanney                                           !
!                                                         !
!---------------------------------------------------------!
PROGRAM USPOPMODEL
  USE SUBS
  IMPLICIT NONE
  REAL:: M,K,A,Z,MMIN,KMIN,AMIN,SSEMIN
  REAL, DIMENSION(1:19):: X, Y, YHAT
  CHARACTER(LEN=80):: F1
  INTEGER:: MCOUNTER, KCOUNTER, ACOUNTER
  OPEN (UNIT = 1, FILE = 'USPOPMOD.DAT', STATUS = 'UNKNOWN')

  X(1:10)  = (/  0,     10,    20,    30,    40,    50,    60,    70,    80,    90  /)
  Y(1:10)  = (/ 3.929, 5.308, 7.240, 9.638,12.866,17.069,23.192,31.443,38.558,50.156/)

  X(11:19) = (/  100,    110,    120,    130,    140,    150,    160,    170,    180  /)
  Y(11:19) = (/ 62.948, 75.995, 91.972,105.711,122.775,131.669,150.697,179.323,203.185/)

  SSEMIN = 1.0E20
  DO MCOUNTER = 100,500
  WRITE (*,*) MCOUNTER           ! Complete coverage of the rectangular solid !
  DO KCOUNTER = 0,800            !         100 <= M <= 500                    !
  DO ACOUNTER = 100,500          !           0 <= K <= 80                     !
    M = MCOUNTER                 !         .01 <= a <= .05                    !
    K = KCOUNTER/10.0            ! From the values determined by 3-point      !
    A = ACOUNTER/10000.0         ! fit, it seems that SSEmin must occur       !
    CALL MODEL(M,K,A,X,YHAT)     ! in this region. We are attempting to       !
    Z = SSE(Y,YHAT)              ! compute each of these parameters to        !
    IF (Z <= SSEMIN) THEN        ! 3 significant digits.                      !
       SSEMIN = Z
       MMIN = M
       KMIN = K
       AMIN = A
    END IF
  END DO
  END DO
  END DO

  F1 = '(A5,F10.3,A6,F10.2,A11,F6.4,A11,F12.0)'
  WRITE (1,F1) 'M =  ',MMIN, 'K = ',KMIN,'a = ',AMIN,'SSEmin = ',SSEMIN
  WRITE (*,*)
  CALL TABLE (MMIN, KMIN, AMIN, X, Y)
STOP
END PROGRAM USPOPMODEL  ! Amen!