RE: calculation of time-to-nadir
From: "Olofsen, E. (ANST)"
Subject: RE:[NMusers] calculation of time-to-nadir
Date: Fri, July 30, 2004 6:52 am
Dear Anthe,
Perhaps the idea implemented in the control file below can be of help.
It considers the case where the analytical solution exists: TMAXE is
the time where the concentration in compartment 2 attains a maximum.
The derivative will never be exactly zero, but it will become
negative. The higher the value of TOL, the more often $DES will be
called at different time values, and the more accurate the estimate,
available as A(4) (at the end of TIME), will be. So the trick is to
integrate 1 until Tmax. Cmax is also available as A(3).
Best regards,
Erik Olofsen
$PROBLEM CMAX
$INPUT TIME AMT RATE MDV X1 X2 DV
$DATA cmax.d
$SUBROUTINES ADVAN6 TOL=6
$MODEL
COMP = CENTRAL
COMP = EFFECT
COMP = CMAX
COMP = TMAX
$PK
K = THETA(1)
KE0 = THETA(2)
TMAXE = (LOG(K)-LOG(KE0))/(K-KE0)
$DES
DADT(1) = -K*A(1)
D2 = KE0*(A(1)-A(2))
DADT(2) = D2
IF (D2.GT.0) THEN
DADT(3) = D2
DADT(4) = 1
ELSE
DADT(3) = 0
DADT(4) = 0
ENDIF
$ERROR
Y1 = A(1)
Y2 = A(2)
Y3 = A(3)
Y4 = A(4)
Y = Y2*(1+ERR(1))
$THETA
0.01
0.03
$OMEGA
0.001
$EST
$TABLE FILE=stream NOPRINT NOHEADER TIME AMT RATE MDV Y1 Y2 Y3 Y4 TMAXE
A data file to test this with:
0.0000E+00 1.0000E+00 0.0000E+00 1.0000E+00 1.0000E+00 0.0000E+00 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00 1.0000E+00 1.0000E+00 0.0000E+00 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00
1.0000E+01 0.0000E+00 0.0000E+00 0.0000E+00 9.0484E-01 2.4603E-01 2.3131E-01
2.4603E-01 -1.4716E-02 0.0000E+00
2.0000E+01 0.0000E+00 0.0000E+00 0.0000E+00 8.1873E-01 4.0488E-01 4.1299E-01
4.0488E-01 8.1097E-03 0.0000E+00
3.0000E+01 0.0000E+00 0.0000E+00 0.0000E+00 7.4082E-01 5.0137E-01 4.9237E-01
5.0137E-01 -9.0010E-03 0.0000E+00
4.0000E+01 0.0000E+00 0.0000E+00 0.0000E+00 6.7032E-01 5.5369E-01 5.3405E-01
5.5369E-01 -1.9643E-02 0.0000E+00
5.0000E+01 0.0000E+00 0.0000E+00 0.0000E+00 6.0653E-01 5.7510E-01 5.8321E-01
5.7510E-01 8.1137E-03 0.0000E+00
6.0000E+01 0.0000E+00 0.0000E+00 0.0000E+00 5.4881E-01 5.7527E-01 5.8127E-01
5.7527E-01 6.0053E-03 0.0000E+00
7.0000E+01 0.0000E+00 0.0000E+00 0.0000E+00 4.9659E-01 5.6119E-01 5.4722E-01
5.6119E-01 -1.3975E-02 0.0000E+00
8.0000E+01 0.0000E+00 0.0000E+00 0.0000E+00 4.4933E-01 5.3792E-01 5.4652E-01
5.3792E-01 8.6028E-03 0.0000E+00
9.0000E+01 0.0000E+00 0.0000E+00 0.0000E+00 4.0657E-01 5.0905E-01 4.9859E-01
5.0905E-01 -1.0452E-02 0.0000E+00
1.0000E+02 0.0000E+00 0.0000E+00 0.0000E+00 3.6788E-01 4.7714E-01 4.9331E-01
4.7714E-01 1.6172E-02 0.0000E+00