c http://Evgenii.Rudnyi.Ru/ c PROGRAM DKCl_ML3 IMPLICIT DOUBLE PRECISION (A-H, O-Z) EXTERNAL SS24, GatSolve, GbtSolve PARAMETER (IXJAC = 600) DIMENSION X(2), PARM(4), F(600), XJAC(IXJAC,4), XJTJ(10), X1(2), * WORK(1230), XtXinv(4), SrOld(40), GaOld(40), GbOld(40), B(4) CHARACTER Name*25, CurTime*8 LOGICAL Eq, Q, Qa, Qb COMMON /Dat/ NexT, NT, NexTr, NTr, NexEf, NEf, NexK5, NK5, * Name(40), NP(40), NFP(40), T(600), P(600), P1(600), Eq(40) COMMON /Work/ AinvT(600), HS1(600), HS2(600), R COMMON /Variance/ Sr(40), Ga(40), Gb(40), Sum(40), Suma(40), * Sumb(40), AvrInvT(40), PT(40), * Q(40), Qa(40), Qb(40), Srt COMMON /XConst/ X2, X4 R = 8.31441 CALL SetKCl1(.True.) DO i = 1, NTr AinvT(i) = 1e3/T(i) CALL HSCOR (T(i), DH1, DS1, DH2, DS2) HS1(i) = (-DH1/T(i) + DS1)/R HS2(i) = EXP((-DH2/T(i) + DS2)/R) END DO DO i = 1, NexTr AvrInvT(i) = 0. DO j = NFP(i), NFP(i) + NP(i) - 1 AvrInvT(i) = AvrInvT(i) + AinvT(j) END DO AvrInvT(i) = AvrInvT(i)/NP(i) PT(i) = 0. DO j = NFP(i), NFP(i) + NP(i) - 1 PT(i) = PT(i) + (AinvT(j) - AvrInvT(i))**2 END DO END DO OPEN (1, NAME='dkcl_ml1.q') DO i = 1, NexTr READ (1, '(3L5)') Q(i), Qa(i), Qb(i) END DO CLOSE (1) OPEN (3, NAME='kcl_ml3.lst') CALL UGETIO (3, 5, 3) N = 2 NSIG = 4 EPS = 0. DELTA = 0. MAXFN = 100 IOPT = 1 X(1) = 206.818744 X2 = 133.030d0 X(2) = 25.739014 X4 = 14.483d0 50 DO i = 1, NexTr Ga(i) = 0. Gb(i) = 0. Sr(i) = 1. END DO WRITE (3, *) 'Initial estimates' WRITE (3, '(4F11.3)') X WRITE (3, *) CALL SS24(X, NTr, N, F) AL = 0. SSQ = 0. 70 Srt = 0. NPt = 0 ALOld = AL + SSQ DO i = 1, NexTr Sr(i) = (Sum(i) - Ga(i)/(1. + NP(i)*Ga(i))*Suma(i) * - Gb(i)/(1. + PT(i)*Gb(i))*Sumb(i))/NP(i) IF (Q(i)) THEN Srt = Srt + Sr(i)*NP(i) NPt = NPt + NP(i) END IF END DO Srt = Srt/NPt CALL BISECT(GatSolve, 0d0, 100d0, Gat) CALL BISECT(GbtSolve, 0d0, 5000d0, Gbt) AL = 0. SSQ = 0. DO i = 1, NexTr IF (Q(i)) Sr(i) = Srt IF (PT(i) .GT. 0) THEN IF (Qb(i)) THEN Gb(i) = Gbt ELSE Gb(i) = (Sumb(i)/PT(i)/Sr(i) - 1.)/PT(i) END IF ELSE Gb(i) = 0. END IF IF (Qa(i)) THEN Ga(i) = Gat ELSE Ga(i) = (Suma(i)/NP(i)/Sr(i) - 1.)/NP(i) END IF IF (Ga(i) .LT. 0) Ga(i) = 0. IF (Gb(i) .LT. 0) Gb(i) = 0. AL = AL + NP(i)*LOG(Sr(i)) + * LOG(1. + NP(i)*Ga(i)) + LOG(1. + PT(i)*Gb(i)) SSQ = SSQ + (Sum(i) - Ga(i)/(1. + NP(i)*Ga(i))*Suma(i) * - Gb(i)/(1. + PT(i)*Gb(i))*Sumb(i))/Sr(i) END DO PRINT*, 'Srt ', Srt, ' Gat ', Gat, ' Gbt ', Gbt WRITE (3, '(3(A,F10.5))') * 'Srt ', Srt, ' Gat ', Gat, ' Gbt ', Gbt PRINT*, 'Likelihood ', SSQ + AL, ' SSQ ', SSQ, ' AL ', AL WRITE (3, '(3(A,F10.3))') * 'Likelihood ', AL + SSQ, ' SSQ ', SSQ, ' AL ', AL IF (AL+SSQ .LT. ALOld - 0.00001) GO TO 70 WRITE (3, *) WRITE (3, *) 'Data in minimization' WRITE (3, 3) WRITE (3, *) 'Total pressure' DO i = 1, NexTr Sr(i) = SQRT(Sr(i)) WRITE (3, 2) Name(i), i, NP(i), Eq(i), Q(i), Qa(i), * Qb(i), 1e3/AvrInvT(i), SQRT(PT(i)), * Sr(i), Sr(i)*SQRT(Ga(i)), Sr(i)*SQRT(Gb(i)) IF (i .EQ. NexT) THEN WRITE (3, *) WRITE (3, *) 'Knudsen effusion' ELSE IF (i .EQ. NexEf) THEN WRITE (3, *) WRITE (3, *) 'Transpiration' END IF END DO WRITE (3, *) Iter = 1 100 CALL TIME(CurTime) PRINT*, CurTime, ' ', X2, X4 PRINT*, 'Iteration ', Iter WRITE (3, '(2A,F10.3,A,F10.3)') * CurTime, ' DelSm ', X2, ' DelSdm ', X4 WRITE (3, *) 'Iteration ', Iter CALL ZXSSQ(SS24,NTr,N,NSIG,EPS,DELTA,MAXFN,IOPT, * PARM,X,SSQ,F,XJAC,IXJAC,XJTJ,WORK,INFER,IER) PRINT*, 'The norm of the gradient ', WORK(1) PRINT*, 'The number of function evaluation ', WORK(2) PRINT*, 'The number of significant digits in parameters ' * , WORK(3) PRINT*, 'The MARQUARDT scaling parameter', WORK(4) PRINT*, 'The number of iteration ', WORK(5) PRINT*, 'INFER = ', INFER, ' IER = ', IER WRITE (3, '(A,F11.7,F5.0,F5.1,F11.7,F5.0)') * 'WORK ', (WORK(j), j = 1, 5) WRITE (3, *) 'INFER = ', INFER, ' IER = ', IER PRINT*, X WRITE (3, 1) X 1 FORMAT (1X, 4F11.3) DO i = 1, NexTr IF (Eq(i)) SSQ = SSQ + 0.03**2*(NP(i)-2)/Sr(i)**2 END DO AL = 0. DO i = 1, NexTr AL = AL + NP(i)*LOG(Sr(i)**2) + * LOG(1. + NP(i)*Ga(i)) + LOG(1. + PT(i)*Gb(i)) END DO PRINT*, 'Likelihood ', AL + SSQ, ' SSQ ', SSQ, ' AL ', AL WRITE (3, '(3(A,F10.3))') * 'Likelihood ', AL + SSQ, ' SSQ ', SSQ, ' AL ', AL WRITE (3, *) DO i = 1, NexTr SrOld(i) = Sr(i) GaOld(i) = Ga(i) GbOld(i) = Gb(i) END DO 80 Srt = 0. NPt = 0 ALOld = AL + SSQ DO i = 1, NexTr Sr(i) = (Sum(i) - Ga(i)/(1. + NP(i)*Ga(i))*Suma(i) * - Gb(i)/(1. + PT(i)*Gb(i))*Sumb(i))/NP(i) IF (Q(i)) THEN Srt = Srt + Sr(i)*NP(i) NPt = NPt + NP(i) END IF END DO Srt = Srt/NPt CALL BISECT(GatSolve, 0d0, 100d0, Gat) CALL BISECT(GbtSolve, 0d0, 5000d0, Gbt) AL = 0. SSQ = 0. DO i = 1, NexTr IF (Q(i)) Sr(i) = Srt IF (PT(i) .GT. 0) THEN IF (Qb(i)) THEN Gb(i) = Gbt ELSE Gb(i) = (Sumb(i)/PT(i)/Sr(i) - 1.)/PT(i) END IF ELSE Gb(i) = 0. END IF IF (Qa(i)) THEN Ga(i) = Gat ELSE Ga(i) = (Suma(i)/NP(i)/Sr(i) - 1.)/NP(i) END IF IF (Ga(i) .LT. 0) Ga(i) = 0. IF (Gb(i) .LT. 0) Gb(i) = 0. AL = AL + NP(i)*LOG(Sr(i)) + * LOG(1. + NP(i)*Ga(i)) + LOG(1. + PT(i)*Gb(i)) SSQ = SSQ + (Sum(i) - Ga(i)/(1. + NP(i)*Ga(i))*Suma(i) * - Gb(i)/(1. + PT(i)*Gb(i))*Sumb(i))/Sr(i) END DO PRINT*, 'Srt ', Srt, ' Gat ', Gat, ' Gbt ', Gbt WRITE (3, '(3(A,F10.5))') * 'Srt ', Srt, ' Gat ', Gat, ' Gbt ', Gbt PRINT*, 'Likelihood ', SSQ + AL, ' SSQ ', SSQ, ' AL ', AL WRITE (3, '(3(A,F10.3))') * 'Likelihood ', AL + SSQ, ' SSQ ', SSQ, ' AL ', AL IF (AL+SSQ .LT. ALOld - 0.0001) GO TO 80 Conv = 0. DO i = 1, NexTr Sr(i) = SQRT(Sr(i)) Conv = Conv + ABS(Sr(i)-SrOld(i)) + * ABS(Sr(i)*SQRT(Ga(i))-SrOld(i)*SQRT(GaOld(i))) + * ABS(Sr(i)*SQRT(Gb(i))-SrOld(i)*SQRT(GbOld(i)))*0.001 END DO PRINT*, 'Convergence ', Conv WRITE (3, *) 'Convergence ', Conv Iter = Iter + 1 IF (Conv. GT. 0.001) GO TO 100 WRITE (3, '(A,F11.7)') 'The norm of the gradient ', WORK(1) WRITE (3, '(A,F5.0)') * 'The number of function evaluation ', WORK(2) WRITE (3, '(A,F5.1)') * 'The number of significant digits in parameters ', WORK(3) WRITE (3, '(A,F11.7)') 'The MARQUARDT scaling parameter', WORK(4) WRITE (3, '(A,F5.0)') 'The number of iteration ', WORK(5) WRITE (3, *) 'INFER = ', INFER, ' IER = ', IER PRINT*, X WRITE (3, 4) X, SSQ, NTr, SQRT(SSQ/(NTr-4)) 4 FORMAT (/,1X, ' DelH1 DelH2 ',/, * 1X, 2F11.3,//, 'Sum of squares', F10.4, ' Number of points', * I5,/, 'Average standard deviation ', F10.4) IJOB = 1 CALL LINV3P (XJTJ, B, IJOB, N, IER) WRITE (3, *) 'LINV3P - IER ', IER WRITE (3, *) WRITE (3, *) ' Parameter sqrt(Cii) ' DO i = 1, N XtXinv(i) = SQRT(XJTJ(i*(i+1)/2)) WRITE (3, '(4F11.2)') X(i), XtXinv(i) END DO WRITE (3, *) WRITE (3, *) 'Correlation matrix' WRITE (3, '(F10.3)') 1 DO i = 2, N DO j = 1, i-1 Loc = i*(i-1)/2+j XJTJ(Loc) = XJTJ(Loc)/SQRT(XJTJ(i*(i+1)/2)*XJTJ(j*(j+1)/2)) END DO WRITE (3, '(4F10.3)') (XJTJ(j), j = i*(i-1)/2+1, i*(i+1)/2-1), 1 END DO WRITE (3, *) WRITE (3, *) 'Data used in minimization' WRITE (3, 3) WRITE (3, *) 'Total pressure' DO i = 1, NexTr WRITE (3, 2) Name(i), i, NP(i), Eq(i), Q(i), Qa(i), * Qb(i), 1e3/AvrInvT(i), SQRT(PT(i)), * Sr(i), Sr(i)*SQRT(Ga(i)), Sr(i)*SQRT(Gb(i)) IF (i .EQ. NexT) THEN WRITE (3, *) WRITE (3, *) 'Knudsen effusion' ELSE IF (i .EQ. NexEf) THEN WRITE (3, *) WRITE (3, *) 'Transpiration' END IF END DO WRITE (3, *) WRITE (3, *) '**********************************************' WRITE (3, *) '**********************************************' 2 FORMAT (A, I2, I5, 4L3, F6.0, F6.3, 3F6.3) 3 FORMAT (' i NP Eq Q ', * 'Qa Qb T X2 ', * 'Sr Sa Sb') 99 X1(1) = X(1) X1(2) = X(2) X2 = X2 + 0.5d0 X4 = X4 CALL TIME(CurTime) PRINT*, CurTime, ' ', X2, X4 WRITE (3, '(2A,F10.3,A,F10.3)') * CurTime, ' DelSm ', X2, ' DelSdm ', X4 CALL ZXSSQ(SS24,NTr,N,NSIG,EPS,DELTA,MAXFN,IOPT, * PARM,X,SSQ,F,XJAC,IXJAC,XJTJ,WORK,INFER,IER) DO i = 1, NexTr IF (Eq(i)) SSQ = SSQ + 0.03**2*(NP(i)-2)/Sr(i)**2 END DO WRITE (3, '(A,F11.7)') 'The norm of the gradient ', WORK(1) WRITE (3, '(A,F5.0)') * 'The number of function evaluation ', WORK(2) WRITE (3, '(A,F5.1)') * 'The number of significant digits in parameters ', WORK(3) WRITE (3, '(A,F11.7)') 'The MARQUARDT scaling parameter', WORK(4) WRITE (3, '(A,F5.0)') 'The number of iteration ', WORK(5) WRITE (3, *) 'INFER = ', INFER, ' IER = ', IER PRINT*, X WRITE (3, 4) X, SSQ, NTr, SQRT(SSQ/(NTr-4)) X(1) = X1(1) X(2) = X1(2) X2 = X2 - 2*0.5d0 X4 = X4 CALL TIME(CurTime) PRINT*, CurTime, ' ', X2, X4 WRITE (3, '(2A,F10.3,A,F10.3)') * CurTime, ' DelSm ', X2, ' DelSdm ', X4 CALL ZXSSQ(SS24,NTr,N,NSIG,EPS,DELTA,MAXFN,IOPT, * PARM,X,SSQ,F,XJAC,IXJAC,XJTJ,WORK,INFER,IER) DO i = 1, NexTr IF (Eq(i)) SSQ = SSQ + 0.03**2*(NP(i)-2)/Sr(i)**2 END DO WRITE (3, '(A,F11.7)') 'The norm of the gradient ', WORK(1) WRITE (3, '(A,F5.0)') * 'The number of function evaluation ', WORK(2) WRITE (3, '(A,F5.1)') * 'The number of significant digits in parameters ', WORK(3) WRITE (3, '(A,F11.7)') 'The MARQUARDT scaling parameter', WORK(4) WRITE (3, '(A,F5.0)') 'The number of iteration ', WORK(5) WRITE (3, *) 'INFER = ', INFER, ' IER = ', IER PRINT*, X WRITE (3, 4) X, SSQ, NTr, SQRT(SSQ/(NTr-4)) X(1) = X1(1) X(2) = X1(2) X2 = X2 + 0.5d0 X4 = X4 + 6.32d0 CALL TIME(CurTime) PRINT*, CurTime, ' ', X2, X4 WRITE (3, '(2A,F10.3,A,F10.3)') * CurTime, ' DelSm ', X2, ' DelSdm ', X4 CALL ZXSSQ(SS24,NTr,N,NSIG,EPS,DELTA,MAXFN,IOPT, * PARM,X,SSQ,F,XJAC,IXJAC,XJTJ,WORK,INFER,IER) DO i = 1, NexTr IF (Eq(i)) SSQ = SSQ + 0.03**2*(NP(i)-2)/Sr(i)**2 END DO WRITE (3, '(A,F11.7)') 'The norm of the gradient ', WORK(1) WRITE (3, '(A,F5.0)') * 'The number of function evaluation ', WORK(2) WRITE (3, '(A,F5.1)') * 'The number of significant digits in parameters ', WORK(3) WRITE (3, '(A,F11.7)') 'The MARQUARDT scaling parameter', WORK(4) WRITE (3, '(A,F5.0)') 'The number of iteration ', WORK(5) WRITE (3, *) 'INFER = ', INFER, ' IER = ', IER PRINT*, X WRITE (3, 4) X, SSQ, NTr, SQRT(SSQ/(NTr-4)) X(1) = X1(1) X(2) = X1(2) X2 = X2 X4 = X4 - 2*6.32d0 CALL TIME(CurTime) PRINT*, CurTime, ' ', X2, X4 WRITE (3, '(2A,F10.3,A,F10.3)') * CurTime, ' DelSm ', X2, ' DelSdm ', X4 CALL ZXSSQ(SS24,NTr,N,NSIG,EPS,DELTA,MAXFN,IOPT, * PARM,X,SSQ,F,XJAC,IXJAC,XJTJ,WORK,INFER,IER) DO i = 1, NexTr IF (Eq(i)) SSQ = SSQ + 0.03**2*(NP(i)-2)/Sr(i)**2 END DO WRITE (3, '(A,F11.7)') 'The norm of the gradient ', WORK(1) WRITE (3, '(A,F5.0)') * 'The number of function evaluation ', WORK(2) WRITE (3, '(A,F5.1)') * 'The number of significant digits in parameters ', WORK(3) WRITE (3, '(A,F11.7)') 'The MARQUARDT scaling parameter', WORK(4) WRITE (3, '(A,F5.0)') 'The number of iteration ', WORK(5) WRITE (3, *) 'INFER = ', INFER, ' IER = ', IER PRINT*, X WRITE (3, 4) X, SSQ, NTr, SQRT(SSQ/(NTr-4)) CLOSE (3) END SUBROUTINE BISECT(F, A, B, RES) IMPLICIT DOUBLE PRECISION (A-H, O-Z) EXTERNAL F DATA EPS/2E-8/ XL = A XR = B FL = F(XL) IF(FL.EQ.0D0) GO TO 4 7 FR = F(XR) IF(FR.EQ.0D0) GO TO 5 IF(SIGN(1.,FL)*SIGN(1.,FR).LT.0.) GO TO 1 WRITE (3, *) 'BISECT FAILED' RES = 0. RETURN 1 XM = (XR+XL) * 0.5 FM = F(XM) IF(FM.EQ.0D0.OR.ABS(XR-XL).LE.EPS*ABS(XM)) GO TO 3 IF(SIGN(1.,FM)*SIGN(1.,FR).GT.0.) GO TO 2 XL = XM FL = FM GO TO 1 2 XR = XM FR = FM GO TO 1 3 RES = XM RETURN 4 RES = XL RETURN 5 RES = XR RETURN END DOUBLE PRECISION FUNCTION GatSolve(X) IMPLICIT DOUBLE PRECISION (A-H, O-Z) CHARACTER Name*25 LOGICAL Eq, Q, Qa, Qb COMMON /Dat/ NexT, NT, NexTr, NTr, NexEf, NEf, NexK5, NK5, * Name(40), NP(40), NFP(40), T(600), P(600), P1(600), Eq(40) COMMON /Variance/ Sr(40), Ga(40), Gb(40), Sum(40), Suma(40), * Sumb(40), AvrInvT(40), PT(40), * Q(40), Qa(40), Qb(40), Srt S = 0. DO i = 1, NexTr IF (Qa(i)) THEN S = S + Suma(i)/Srt/(1.+NP(i)*X)**2 - NP(i)/(1.+NP(i)*X) END IF END DO GatSolve = S RETURN END DOUBLE PRECISION FUNCTION GbtSolve(X) IMPLICIT DOUBLE PRECISION (A-H, O-Z) CHARACTER Name*25 LOGICAL Eq, Q, Qa, Qb COMMON /Dat/ NexT, NT, NexTr, NTr, NexEf, NEf, NexK5, NK5, * Name(40), NP(40), NFP(40), T(600), P(600), P1(600), Eq(40) COMMON /Variance/ Sr(40), Ga(40), Gb(40), Sum(40), Suma(40), * Sumb(40), AvrInvT(40), PT(40), * Q(40), Qa(40), Qb(40), Srt S = 0. DO i = 1, NexTr IF (Qb(i)) THEN S = S + Sumb(i)/Srt/(1.+PT(i)*X)**2 - PT(i)/(1.+PT(i)*X) END IF END DO GbtSolve = S RETURN END SUBROUTINE SS24(X, M, N, F) IMPLICIT DOUBLE PRECISION (A-H, O-Z) DIMENSION X(N), F(M) CHARACTER Name*25 LOGICAL Eq, Q, Qa, Qb COMMON /Dat/ NexT, NT, NexTr, NTr, NexEf, NEf, NexK5, NK5, * Name(40), NP(40), NFP(40), T(600), P(600), P1(600), Eq(40) COMMON /Work/ AinvT(600), HS1(600), HS2(600), R COMMON /Variance/ Sr(40), Ga(40), Gb(40), Sum(40), Suma(40), * Sumb(40), AvrInvT(40), PT(40), * Q(40), Qa(40), Qb(40), Srt COMMON /XConst/ X2, X4 DO i = 1, NT PM = (-X(1)*AinvT(i) + X2)/R + HS1(i) PDM = EXP((-X(2)*AinvT(i) + X4)/R)*HS2(i) F(i) = PM + LOG(1. + PDM) - P(i) END DO DO i = NT + 1, NEf PM = (-X(1)*AinvT(i) + X2)/R + HS1(i) PDM = EXP((-X(2)*AinvT(i) + X4)/R)*HS2(i) F(i) = PM + LOG(1. + SQRT(2.d0)*PDM) - P(i) END DO DO i = NEf + 1, NTr PM = (-X(1)*AinvT(i) + X2)/R + HS1(i) PDM = EXP((-X(2)*AinvT(i) + X4)/R)*HS2(i) F(i) = PM + LOG(1. + 2.*PDM) - * LOG(1. + EXP(PM)*PDM/P1(i)) - P(i) END DO DO i = 1, NexTr Sum(i) = 0. Suma(i) = 0. Sumb(i) = 0. DO j = NFP(i), NFP(i) + NP(i) - 1 Sum(i) = Sum(i) + F(j)**2 Suma(i) = Suma(i) + F(j) Sumb(i) = Sumb(i) + F(j)*(AinvT(j) - AvrInvT(i)) END DO IF (Eq(i)) Sum(i) = Sum(i) + 0.03**2*(NP(i)-2) Da = (1. - 1./SQRT(1. + NP(i)*Ga(i)))/NP(i) IF (PT(I) .GT. 0) THEN Db = (1. - 1./SQRT(1. + PT(i)*Gb(i)))/PT(i) ELSE Db = 0. END IF DO j = NFP(i), NFP(i) + NP(i) - 1 F(j) = (F(j) - Da*Suma(i) - * Db*(AinvT(j) - AvrInvT(i))*Sumb(i))/Sr(i) END DO Suma(i) = Suma(i)**2 Sumb(i) = Sumb(i)**2 END DO RETURN END