MODULE My_Data
IMPLICIT NONE
INTEGER, PARAMETER  :: dp = SELECTED_REAL_KIND(12, 60)

REAL (dp), SAVE  :: x(100), y(100)
END MODULE My_Data



PROGRAM Test_AS290
USE Nonlin_Conf_Regions
USE My_Data
IMPLICIT NONE

CHARACTER (LEN=20)  :: fname
CHARACTER (LEN=80)  :: text
INTEGER             :: dim1, ifault, fstatus, i, length, line1, ncols, nconf, &
                       ngrid1, ngrid2, nobs, np, np1, np2, row
REAL (dp)           :: conf(6), f(6), grid(31,31), p1lo, p1hi, p2lo, p2hi,  &
                       param(4)
LOGICAL             :: lindep(4)

INTERFACE
  SUBROUTINE Logistic4(param, np, ncols, iobs, grad, resid, wt)
    USE My_Data
    IMPLICIT NONE
    REAL (dp), INTENT(IN)   :: param(:)
    INTEGER, INTENT(IN)     :: np, ncols, iobs
    REAL (dp), INTENT(OUT)  :: grad(:), resid, wt
  END SUBROUTINE Logistic4
END INTERFACE

DO
  WRITE(*, *) 'Enter name of data file: '
  READ(*, *) fname
  OPEN(10, file=fname, IOSTAT=fstatus)
  IF (fstatus == 0) EXIT
  WRITE(*, *) ' !! File not found -- try again !!'
END DO

WRITE(*, *) ' This is the start of the file:'
DO row = 1, 4
  READ(10, '(a)') text
  length = LEN_TRIM(text)
  WRITE(*, '(" ", a75)') text(1:length)
END DO
REWIND(10)
WRITE(*, *)
WRITE(*, *) ' What line does data start on? '
READ(*, *) line1

IF (line1 > 1) THEN
  DO row = 1, line1-1
    READ(10, *)
  END DO
END IF

nobs = 1
DO
  READ(10, *, iostat=fstatus) x(nobs), y(nobs)
  IF (fstatus < 0) EXIT
  nobs = nobs + 1
END DO
nobs = nobs - 1
WRITE(*, *) ' No. of cases read = ', nobs

! Set up the 2_D grid of variance ratios.
! The values below are for the case in the Applied Statistics paper.
! Change for your needs.

param(1) = 0.88_dp
param(2) = 21.3_dp
np = 4
ncols = 4
np1 = 3
np2 = 4
p1lo = 0.4_dp
p1hi = 1.5_dp
p2lo = 5.8_dp
p2hi = 7.2_dp
conf(1:6) = (/ 10._dp, 50._dp, 90._dp, 95._dp, 98._dp, 99._dp /)
nconf = 6
dim1 = 31
ngrid1 = 31
ngrid2 = 31
CALL halprn(param, np, ncols, nobs, np1, np2, p1lo, p1hi, p2lo, p2hi,  &
                  conf, nconf, f, grid, dim1, ngrid1, ngrid2, Logistic4,  &
                  lindep, ifault)
WRITE(*, *) ' Percentage points of the F-distribution'
WRITE(*, *) '   %    F-value'
DO i = 1, nconf
  WRITE(*, '(f6.1, f8.2)') conf(i), f(i)
END DO

! Write the grid to an output file for later use with a contouring program.
OPEN(11, file='confgrid.dat')
DO row = 1, ngrid1
  WRITE(11, '(10g13.5)') grid(row,1:ngrid2)
END DO

STOP
END PROGRAM Test_AS290



SUBROUTINE Logistic4(param, np, ncols, iobs, grad, resid, wt)

! Calculate residual and first derivatives for the 4-parameter logistic model:
!      Y = a + (b-a) / [1 + exp{-c(X-d)}]
! for one case IOBS.

USE My_Data
IMPLICIT NONE

REAL (dp), INTENT(IN)   :: param(:)
INTEGER, INTENT(IN)     :: np, ncols, iobs
REAL (dp), INTENT(OUT)  :: grad(:), resid, wt

! Local variables
REAL (dp)  :: denom, diff, expntl
REAL (dp), PARAMETER  :: one = 1.0_dp

wt = one
diff = x(iobs) - param(4)
expntl = EXP(-param(3)*diff)
denom = one + expntl
resid = y(iobs) - param(1) - (param(2) - param(1))/denom
grad(1) = -one + one/denom
grad(2) = -one/denom
grad(3) = -(param(2) - param(1))*diff*expntl/denom**2
grad(4) = (param(2) - param(1))*param(3)*expntl/denom**2

RETURN
END SUBROUTINE Logistic4