SUBROUTINE LOG2PT(KFILDI,KFILDO,X1,Y1,X2,Y2,GL)
C
C JUNE 2018 GLAHN MOS-2000
C JANUARY 2020 GLAHN CHANGED y (LC) TO Y (UC) IN
C SEVERAL PLACES
C
C
C PURPOSE
C TO FIT A LOGIT CURVE TO TWO POINTS. THIS CAN BE USED
C TO BLEND TWO GRIDS TOGETHER, POINT BY POINT, AFTER
C OVERLAPPING THE GRIDS. THE FUNCTION DETERMINED IS
C EVALLUATED FOR EACH INTEGER BETWEEN AND INCLUDING
C X1 AND X2 INNCLUSIVELY AND THE RESULT RETURNED
C IN GL( ).
C
C A LINEAR LINE IS FIRST FIT TO THE TWO POINTS TO FURNISH
C A FIRST GUESS FOR THE ITERATIVE LOGISTIC SOLUTION.
C
C DIAGNOSTIC INFORMATION AND THE TIME STAMPS CAN BE
C ELIMINATED, IF DESIRED.
C
C DATA SET USE
C KFILDI - UNIT NUMBER OF INPUT FILE. (INPUT) (NOT USED)
C KFILDO - UNIT NUMBER OF OUTPUT (PRINT) FILE. (OUTPUT)
C
C VARIABLES
C KFILDI = UNIT NUMBER OF INPUT FILE. (INPUT) (NOT USED)
C KFILDO = UNIT NUMBER OF OUTPUT (PRINT) FILE. (OUTPUT)
C X1 = FIRST X POSITION. (INPUT)
C Y1 = FIRST Y POSITION. (INPUT)
C X2 = SECOND X POSITION. (INPUT)
C Y2 = SECOND Y POSITION. (INPUT)
C GL(L) = LOGIT VALUES FOR S = X1 THROUGH X2. THE ARRAY
C SIZE OF GL( ) IS NINT(X2). (OUTPUT)
C 1 2 3 4 5 6 7 X
C
C NONSYSTEM SUBROUTINES USED
C NONE
C
DIMENSION BSTAR(2),B(2),C(2,2),D(2),X(2),Y(2),P(2),
1 CINV(2,2),CINVD(2),GL(NINT(X2)),HL(NINT(X2))
C
C FIRST FIT A LINEAR LINE 'Y = A + BX' TO THE TWO POINTS TO
C FURNISH THE FIRST GUESS FOR THE LOGIT.
C
CALL TIMPR(KFILDO,KFILDO,'START LOG2PT ')
C
CCCC WRITE(KFILDO,105)X1,Y1,X2,Y2
CCCC 105 FORMAT(/' AT 105 IN LOG2PT--X1,Y1,X2,Y2',4F10.4)
C
B1=(Y2-Y1)/(X2-X1)
B0=Y1-X1*B1
CCCC WRITE(KFILDO,110)B0,B1
CCCC 110 FORMAT(/' LINE IS B0 + B1*X WITH B0 = ',F10.4,
CCCC 1 ' AND B1 =',F10.4)
C
C VERIFY.
C
Y11=B0+B1*X1
Y22=B0+B1*X2
CCCC WRITE(KFILDO,115)Y11,Y22
CCCC 115 FORMAT(/' VERIFY, Y11, Y22 =',2F10.3)
C
C COMPUTE INITIAL LOGIT FIT FOR Y11 AND Y22
C
X(1)=X1
X(2)=X2
Y(1)=Y1
Y(2)=Y2
B(1)=B0
B(2)=B1
C
Y11=EXP(B(1)+B(2)*X(1))
Y22=EXP(B(1)+B(2)*X(2))
Y11=Y11/(1.+Y11)
Y22=Y22/(1.+Y22)
C
CCCC WRITE(KFILDO,118)Y11,Y22
CCCC 118 FORMAT(/' Y11,Y22',2F10.4)
C
ICOUNT=0
C
120 ICOUNT=ICOUNT+1
C(1,1)=0.
C(2,1)=0.
C(1,2)=0.
C(2,2)=0.
D(1)=0.
D(2)=0.
C
DO 140 L=1,2
P(L)=EXP(B(1)+B(2)*X(L))
P(L)=P(L)/(1.+P(L))
C(1,1)=C(1,1)+P(L)**2-P(L)
C(2,1)=C(2,1)+X(L)*(P(L)**2-P(L))
C(1,2)=C(1,2)+X(L)*(P(L)**2-P(L))
C(2,2)=C(2,2)+X(L)**2*(P(L)**2-P(L))
D(1)=D(1)+Y(L)-P(L)
D(2)=D(2)+X(L)*(Y(L)-P(L))
C
140 CONTINUE
C
DETC=C(1,1)*C(2,2)-C(1,2)*C(2,1)
CINV(1,1)=C(2,2)/DETC
CINV(2,1)=-C(2,1)/DETC
CINV(1,2)=-C(1,2)/DETC
CINV(2,2)=C(1,1)/DETC
C
CINVD(1)=CINV(1,1)*D(1)+CINV(1,2)*D(2)
CINVD(2)=CINV(2,1)*D(1)+CINV(2,2)*D(2)
C
BSTAR(1)=B(1)-CINVD(1)
BSTAR(2)=B(2)-CINVD(2)
C
Y11=EXP(BSTAR(1)+BSTAR(2)*X(1))
Y22=EXP(BSTAR(1)+BSTAR(2)*X(2))
Y11=Y11/(1.+Y11)
Y22=Y22/(1.+Y22)
C
CCCC WRITE(KFILDO,150)BSTAR(1),BSTAR(2),Y11,Y22
CCCC 150 FORMAT(/' BSTAR(1),BSTAR(2),Y11,Y22',4F10.4)
C
B(1)=BSTAR(1)
B(2)=BSTAR(2)
IF(ICOUNT.LT.10)GO TO 120
C
CCCC WRITE(KFILDO,164)
CCCC 164 FORMAT(/' X LOGIT LINEAR')
C
DO 160 L=1,NINT(X2)
GL(L)=EXP(B(1)+B(2)*L)
GL(L)=GL(L)/(1.+GL(L))
HL(L)=L/(X2+1)
CCCC WRITE(KFILDO,165)L,GL(L),HL(L)
CCCCC THE PRINT FOR EACH VALUE OF X IS (1) THE LOGIT VALUE,
CCCCC AND (2) A LINEAR FIT.
CCCC 165 FORMAT(/' L,GL(L),HL(L)',I4,2F10.4)
160 CONTINUE
C
CCCC CALL TIMPR(KFILDO,KFILDO,'END LOG2PT ')
C
RETURN
END