-
Notifications
You must be signed in to change notification settings - Fork 2
Expand file tree
/
Copy pathEX214.FOR
More file actions
145 lines (145 loc) · 4.11 KB
/
Copy pathEX214.FOR
File metadata and controls
145 lines (145 loc) · 4.11 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
C
C***************************** ABSTRACT *******************************
C
C THIS PROGRAM APPLIES THE SECANT METHOD TO SOLVE FOR THE BUBBLE
C POINT TEMPERATUR OF AN IDEAL SYSTEM.
C
C**************************** NOMENCLATURE ****************************
C
C A(I)- ANTOINE CONSTANT FOR THE ITH COMPONENT
C B(I)- ANTOINE CONSTANT FOR THE ITH COMPONENT
C C(I)- ANTOINE CONSTANT FOR THE ITH COMPONENT
C DXMAX- THE MAXIMUM ALLOWED CHANGE IN X BETWEEN ITERATIONS
C (USER SPECIFIED)
C ERLIM- CONVERGENCE CRITERIA
C ICMAX- THE MAXIMUM NUMBER OF ITERATIONS ALLOWED
C IPRINT- PRINT FLAG =1 PRINT INTERMEDIATE VALUES; =0 NO PRINT
C P- SYSTEM PRESSURE
C T- THE INITIAL GUESS FOR THE ROOT (INPUT); THE SOLUTION VALUE OF
C T (OUTPUT)
C X(I)- MOLE FRACTION OF THE ITH COMPONENT
C
C************************************************************************
C
IMPLICIT REAL*8(A-H,O-Z)
COMMON /DATA/ A(4),B(4),C(4),X(4),P
C SPECIFY NUMERICAL PARAMETERS
ERLIM=1.E-6
ICMAX=100
IPRINT=1
DXMAX=.1
C SPECIFY ANTOINE CONSTANTS FOR EACH COMPONENT
A(1)=15.8366
A(2)=15.8737
A(3)=15.9426
A(4)=15.9671
B(1)=2697.5
B(2)=2911.32
B(3)=3120.29
B(4)=3291.45
C(1)=-48.78
C(2)=-56.51
C(3)=-63.63
C(4)=-71.33
C SET LIQUID COMPOSITION
X(1)=.15
X(2)=.25
X(3)=.45
X(4)=.15
C SET PRESSURE
P=1.0
C SPECIFY STARTING POINT
TI=1./(273.15)
C
C CALL SECANT
C
CALL SECANT(ERLIM,ICMAX,DXMAX,IPRINT,TI)
C
C WRITE BUBBLE POINT TEMPERATURE
C
TEMP=1./TI -272.15
WRITE(6,11)TEMP
11 FORMAT( 5X,' BUBBLE POINT TEMPERATURE=',F7.2,' DEG C')
C
C END
C
STOP
END
C
C***************************** ABSTRACT *******************************
C
C THIS SUBROUTINE APPLIES THE SECANT METHOD FOR THE SOLUTION OF A
C SINGLE NONLINEAR EQUATION.
C
C**************************** NOMENCLATURE ****************************
C
C DX- THE CALCULATED CHANGE IN THE VALUE OF X
C DXMAX- THE MAXIMUM ALLOWED CHANGE IN X BETWEEN ITERATIONS
C (USER SPECIFIED)
C ERLIM- CONVERGENCE CRITERIA
C FO- THE NEXT TO THE LATEST FUNCTION VALUE
C F1- THE LATEST FUNCTION VALUE
C FV- THE FUNCTION VALUE
C IC- THE ITERATION COUNTER
C ICMAX- THE MAXIMUM NUMBER OF ITERATIONS ALLOWED
C IPRINT- PRINT FLAG =1 PRINT INTERMEDIATE VALUES; =0 NO PRINT
C X1- THE LATEST VALUE THE ROOT OF THE NONLINEAR EQUATION
C XO- THE NEXT TO THE LATEST VALUE THE ROOT OF THE NONLINEAR EQUATION
C
C************************************************************************
C
SUBROUTINE SECANT(ERLIM,ICMAX,DXMAX,IPRINT,X1)
IMPLICIT REAL*8(A-H,O-Z)
IC=0
C
C BEGIN ITERATIVE LOOP
C
10 IC=IC+1
C CALL SUBROUTINE F
CALL F(X1,F1)
C INITIALIZE SECANT METHOD
IF(IC.EQ.1)XO=X1*.95
IF(IC.EQ.1)CALL F(XO,FO)
C APPLY SECANT APPROXIATION
DX=-F1*(X1-XO)/(F1-FO)
IF(DABS(DX).GT.DXMAX)DX=DXMAX*DABS(DX)/DX
XC=X1+DX
C TRANSFER X1,F1 TO XO,FO
XO=X1
FO=F1
X1=XC
C PRINT OUT INTERMEDIATE RESULTS
IF(IPRINT.EQ.1)WRITE(6,11)IC,XC,F1
11 FORMAT( 5X,' IC=',I3,5X,' X=',F7.5,5X,' F(X)=',D10.3)
C CHECK FOR EXCESSIVE ITERATIONS
IF(IC.GT.ICMAX)WRITE(6,22)
22 FORMAT( ' SECANT METHOD DID NOT CONVERGE')
IF(IC.GT.ICMAX)STOP
C CHECK FOR CONVERGENCE
IF(DABS(DX).GT.ERLIM)GO TO 10
C
C CONVERGED SOLUTION OBTAINED; RETURNING TO CALLING PROGRAM
C
RETURN
END
C
C***************************** ABSTRACT *******************************
C
C THIS SUBROUTINE CALCULATES THE CALCULATES THE VALUE OF F(X)
C GIVEN THE VALUE OF X.
C
C************************************************************************
C
SUBROUTINE F(TI,FV)
IMPLICIT REAL*8(A-H,O-Z)
COMMON /DATA/ A(4),B(4),C(4),X(4),P
C CONVERT TO ABSOLUTE TEMPERATURE
T=1./TI
C CALCULATE SUM OF Y(I)
SUM=0.0
DO 10 I=1,4
10 SUM=SUM+X(I)*DEXP(A(I)-B(I)/(C(I)+T))/P/760.
C LOG OF SUM SHOULD EQUAL TO 1 SO LOG OF SUM IS FUNCTION VALUE
FV=DLOG(SUM)
RETURN
END