! 4.13.f90
! created on 2024-10-20
! $Id: 4.13.f90 1.1 2024/11/02 09:22:40 s Exp $

implicit none

! declaration
integer i, j, n
parameter (n=8)
integer t(n)
real(8) a(n), aa(2,2),aa1(2,2), a2(8), bb(2),bb1(2), df(n),df2(n),dfx(n), s(n), &
        s_percent(n), txdf2(n), x(n), xx1(2)
real(8) p1, d1, duration, p2, d2, pv, sumtxdf2, x1, x2
character*40 fmt
character*10 keisen
keisen='----------'

! input
s_percent = [7.67, 8.27, 8.81, 9.31, 9.75, 10.16, 10.52, 10.85]
x =         [ 500,  900,  600,  500,  100,   100,   100,    50]
p1 =  65.95
d1 =   7.07
p2 = 101.66
d2 =   3.80

! calculating
s   = s_percent / 100
!n   = 8
t   = [1, 2, 3, 4, 5, 6, 7, 8]
a   = (1 + s) ** t
df  = 1 / a
a2  = (1 + s) ** (t + 1)
df2 = 1 / a2

dfx      = df * x
pv       = sum(dfx)
txdf2    = t * x * df2
sumtxdf2 = sum(txdf2)
duration = sumtxdf2 / pv

aa(1,1) = p1
aa(1,2) = p2
aa(2,1) = p1 * d1
aa(2,2) = p2 * d2
bb = [pv, pv * duration]



aa1=aa
bb1=bb
call GAUSS(aa1, bb1, xx1, 2, 2)
x1 = xx1(1)
x2 = xx1(2)

! output
print *,'4.13.f90'
print '(6a10)', "t", "s", "df", "x", "dfx", "txdf2"
print '(a)', repeat('-',60)
fmt='(i10,f10.4,f10.3,f10.2,f10.2,f10.2)'
print fmt, (t(i), s(i), df(i), x(i), dfx(i), txdf2(i),i=1,n)
print '(a)', repeat('-',60)

print '(a40,a3,f7.2,f10.2)', "","PV=", pv, sumtxdf2
print '(a50a5f5.2)', "", "D=", duration

print '(a10,f10.2)', "p1= ", p1
print '(a10,f10.2)', "p2= ", p2
print '(a10,f10.2)', "d1= ", d1
print '(a10,f10.2)', "d2= ", d2
print *, 'A'
print '(2f10.2)', ((aa(i, j), j=1,2), i=1,2)
print '(a)', repeat('-',20)
print *, 'b'
print '(2f10.2)', (bb(i), i=1,2)
print '(a)', repeat('-',20)
print '(a10,f10.2)', "x1= ", x1
print '(a10,f10.2)', "x2= ", x2
stop
end

SUBROUTINE GAUSS(A,B,X,N,NN) 
    ! 出典  戸川隼人『数値計算』岩波書店, 1991, p.80, プログラム4.1
    implicit none
    real(8) A(NN,NN),B(NN),X(NN),Q,S
    INTEGER I,J,K,N,NN
    DO 10 K=1,N-1
        DO 20 I=K+1,N
            Q=A(I,K)/A(K,K)
            DO 30 J=K,N
                A(I,J)=A(I,J)-Q*A(K,J)
30          CONTINUE
            B(I)=B(I)-Q*B(K)
20      CONTINUE
10  CONTINUE
    X(N)=B(N)/A(N,N)
    DO 40 K=N-1,1,-1
        S=B(K)
        DO 50 J=K+1,N
            S=S-A(K,J)*X(J)
50      CONTINUE
        X(K)=S/A(K,K)
40  CONTINUE
    RETURN
END 
! eof