Skip to content

Commit

Permalink
ENH: API break for g704 kh and kd. T is now in Kelvin.
Browse files Browse the repository at this point in the history
  • Loading branch information
@MilanSkocic
MilanSkocic committed Jul 1, 2024
1 parent 66ad77a commit 72b1cf9
Show file tree
Hide file tree
Showing 10 changed files with 33 additions and 33 deletions.
6 changes: 3 additions & 3 deletions example/example.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@

int main(void){

double T = 25.0; /* in C*/
double T = 25.0 + 273.15; /* in C*/
char *gas = "O2";
double kh, kd;
char **gases_list;
Expand All @@ -31,10 +31,10 @@ int main(void){
printf("%s\n", "########################## IAPWS G7-04 ##########################");
/* Compute kh and kd in H2O*/
iapws_g704_kh(&T, gas, heavywater, &kh, strlen(gas), 1);
printf("Gas=%s\tT=%fC\tkh=%+10.4f\n", gas, T, kh);
printf("Gas=%s\tT=%fK\tkh=%+10.4f\n", gas, T, kh);

iapws_g704_kd(&T, gas, heavywater, &kd, strlen(gas), 1);
printf("Gas=%s\tT=%fC\tkd=%+15.4f\n", gas, T, kd);
printf("Gas=%s\tT=%fK\tkd=%+15.4f\n", gas, T, kd);

/* Get and print the available gases */
ngas = iapws_g704_ngases(heavywater);
Expand Down
6 changes: 3 additions & 3 deletions example/example.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,12 +25,12 @@ program example_in_f

print *, '########################## IAPWS G7-04 ##########################'
! Compute kh and kd in H2O
T(1) = 25.0d0
T(1) = 25.0_dp + 273.15_dp
call kh(T, gas, heavywater, kh_res)
print "(A10, 1X, A10, 1X, A2, F10.1, A, 4X, A3, SP, F10.4)", "Gas=", gas, "T=", T, "C", "kh=", kh_res
print "(A10, 1X, A10, 1X, A2, F10.1, A, 4X, A3, SP, F10.4)", "Gas=", gas, "T=", T, "K", "kh=", kh_res

call kd(T, gas, heavywater, kd_res)
print "(A10, 1X, A10, 1X, A2, F10.1, A, 4X, A3, SP, F15.4)", "Gas=", gas, "T=", T, "C", "kh=", kd_res
print "(A10, 1X, A10, 1X, A2, F10.1, A, 4X, A3, SP, F15.4)", "Gas=", gas, "T=", T, "K", "kh=", kd_res

! Get and print available gases
heavywater = 0
Expand Down
8 changes: 4 additions & 4 deletions example/example.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,15 +22,15 @@

print("########################## IAPWS G7-04 ##########################")
gas = "O2"
T = array.array("d", (25.0,))
T = array.array("d", (25.0+273.15,))

# Compute kh and kd in H2O
heavywater = False
k = pyiapws.g704.kh(T, "O2", heavywater)
print(f"Gas={gas}\tT={T[0]}C\tkh={k[0]:+10.4f}\n")
print(f"Gas={gas}\tT={T[0]}K\tkh={k[0]:+10.4f}\n")

k = pyiapws.g704.kd(T, "O2", heavywater)
print(f"Gas={gas}\tT={T[0]}C\tkh={k[0]:+10.4f}\n")
print(f"Gas={gas}\tT={T[0]}K\tkd={k[0]:+10.4f}\n")

# Get and print the available gases
heavywater = False
Expand All @@ -53,7 +53,7 @@

style = {"marker":".", "ls":"", "ms":2}
T_KELVIN = 273.15
T = np.linspace(0.0, 360.0, 1000)
T = np.linspace(0.0, 360.0, 1000) + 273.15

solvent = {True: "D2O", False: "H2O"}

Expand Down
Binary file modified media/g704-kd_D2O.png
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file modified media/g704-kd_H2O.png
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file modified media/g704-kh_D2O.png
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file modified media/g704-kh_H2O.png
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
4 changes: 2 additions & 2 deletions py/src/pyiapws/g704.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@ def kh(T: np.ndarray, gas: str, heavywater: bool=False)->Union[np.ndarray, float
Parameters
----------
T: int, float or 1d-array.
Temperature in °C.
Temperature in K.
gas: str
Gas.
heavywater: bool
Expand Down Expand Up @@ -45,7 +45,7 @@ def kd(T: np.ndarray, gas: str, heavywater: bool=False)->Union[np.ndarray, float
Parameters
----------
T: int, float, or 1d-array.
Temperature in °C.
Temperature in K.
gas: str
Gas.
heavywater: bool
Expand Down
34 changes: 17 additions & 17 deletions src/iapws_g704.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -173,14 +173,14 @@ pure elemental function f_p1star_H2O(T)result(value)
!! Compute p1* in H2O.
implicit none
real(dp), intent(in) :: T
!! Temperature in °C.
!! Temperature in K.
real(dp) :: value
!! p1* in MPa.

real(dp) :: Tr
real(dp) :: tau

Tr = (T+T_KELVIN)/Tc1_H2O
Tr = T/Tc1_H2O
tau = 1 - Tr
value = exp(1/(Tr) * sum(aibi_H2O(:,1)*tau**(aibi_H2O(:,2)))) * pc1_H2O
end function
Expand All @@ -189,14 +189,14 @@ pure elemental function f_p1star_D2O(T)result(value)
!! Compute p1* in D2O.
implicit none
real(dp), intent(in) :: T
!! Temperature in °C.
!! Temperature in K.
real(dp) :: value
!! p1* in MPa.

real(dp) :: Tr
real(dp) :: tau

Tr = (T+T_KELVIN)/Tc1_D2O
Tr = T/Tc1_D2O
tau = 1 - Tr
value = exp(1/(Tr) * sum(aibi_D2O(:,1)*tau**(aibi_D2O(:,2)))) * pc1_D2O
end function
Expand All @@ -205,7 +205,7 @@ pure elemental function f_kh_p1star_H2O(T, abc)result(value)
!! Compute kh/p1* in H2O.
implicit none
real(dp), intent(in) :: T
!! Temperature in °C.
!! Temperature in K.
type(abc_t), intent(in) :: abc
!! ABC coefficients.
real(dp) :: value
Expand All @@ -214,7 +214,7 @@ pure elemental function f_kh_p1star_H2O(T, abc)result(value)
real(dp) :: Tr
real(dp) :: tau

Tr = (T+T_KELVIN)/Tc1_H2O
Tr = T/Tc1_H2O
tau = 1 - Tr
value = exp(abc%A/Tr + abc%B*(tau**0.355_dp)/Tr + abc%C*exp(tau)*Tr**(-0.41_dp))
end function
Expand All @@ -223,7 +223,7 @@ pure elemental function f_kh_p1star_D2O(T, abc)result(value)
!! Compute kh/p1* in D2O.
implicit none
real(dp), intent(in) :: T
!! Temperature in °C.
!! Temperature in K.
type(abc_t), intent(in) :: abc
!! ABC coefficients.
real(dp) :: value
Expand All @@ -232,7 +232,7 @@ pure elemental function f_kh_p1star_D2O(T, abc)result(value)
real(dp) :: Tr
real(dp) :: tau

Tr = (T+T_KELVIN)/Tc1_D2O
Tr = T/Tc1_D2O
tau = 1 - Tr
value = exp(abc%A/Tr + abc%B*(tau**0.355_dp)/Tr + abc%C*exp(tau)*Tr**(-0.41_dp))
end function
Expand Down Expand Up @@ -261,7 +261,7 @@ pure elemental function f_kh_H2O(T, abc)result(value)
!! Compute kH in H2O.
implicit none
real(dp), intent(in) :: T
!! Temperature in °C.
!! Temperature in K.
type(abc_t), intent(in) :: abc
!! ABC coefficients.
real(dp) :: value
Expand All @@ -273,7 +273,7 @@ pure elemental function f_kh_D2O(T, abc)result(value)
!! Compute kH in D2O.
implicit none
real(dp), intent(in) :: T
!! Temperature in °C.
!! Temperature in K.
type(abc_t), intent(in) :: abc
!! ABC coefficients.
real(dp) :: value
Expand All @@ -285,7 +285,7 @@ pure elemental function f_kd_H2O(T, efgh) result(value)
!! Compute kd in H2O.
implicit none
real(dp), intent(in) :: T
!! Temperature in °C.
!! Temperature in K.
type(efgh_t), intent(in) :: efgh
!! EFGH coefficients.
real(dp) :: value
Expand All @@ -302,7 +302,7 @@ pure elemental function f_kd_H2O(T, efgh) result(value)
tau = 1-Tr

p1 = q_H2O*efgh%F
p2 = efgh%E/(T+T_KELVIN)*ft_H2O(tau)
p2 = efgh%E/T*ft_H2O(tau)
p3 = (efgh%F + efgh%G*tau**(2.0_dp/3.0_dp) + efgh%H*tau)
p4 = exp(-T/100.0_dp)

Expand All @@ -314,7 +314,7 @@ pure elemental function f_kd_D2O(T, efgh) result(value)
!! Compute kd in D2O.
implicit none
real(dp), intent(in) :: T
!! Temperature in °C.
!! Temperature in K.
type(efgh_t), intent(in) :: efgh
!! EFGH coefficients.
real(dp) :: value
Expand All @@ -327,11 +327,11 @@ pure elemental function f_kd_D2O(T, efgh) result(value)
real(dp) :: p3
real(dp) :: p4

Tr = (T+T_KELVIN)/Tc1_D2O
Tr = T/Tc1_D2O
tau = 1-Tr

p1 = q_D2O*efgh%F
p2 = efgh%E/(T+T_KELVIN)*ft_D2O(tau)
p2 = efgh%E/T*ft_D2O(tau)
p3 = (efgh%F + efgh%G*tau**(2.0_dp/3.0_dp) + efgh%H*tau)
p4 = exp(-T/100.0_dp)

Expand All @@ -344,7 +344,7 @@ pure subroutine kh(T, gas, heavywater, k)
implicit none

! arguments
real(dp), intent(in) :: T(:) !! Temperature in °C.
real(dp), intent(in) :: T(:) !! Temperature in K.
character(len=*), intent(in) :: gas !! Gas.
integer(int32), intent(in) :: heavywater !! Flag if D2O (1) is used or H2O(0).
real(dp), intent(out) :: k(:) !! Henry constant. Filled with NaNs if gas not found.
Expand Down Expand Up @@ -375,7 +375,7 @@ pure subroutine kd(T, gas, heavywater, k)
implicit none

! arguments
real(dp), intent(in) :: T(:) !! Temperature in °C.
real(dp), intent(in) :: T(:) !! Temperature in K.
character(len=*), intent(in) :: gas !! Gas.
integer(int32), intent(in) :: heavywater !! Flag if D2O (1) is used or H2O(0).
real(dp), intent(out) :: k(:) !! Vapor-liquid constant. Filled with NaNs if gas not found.
Expand Down
8 changes: 4 additions & 4 deletions test/testsuite_g704.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -152,7 +152,7 @@ subroutine test_kh_H2O(error)
T_C = T_K - 273.15d0

do i=1, ngas
call kh(T_C, gases(i), heavywater, k)
call kh(T_K, gases(i), heavywater, k)
do j=1, nT
value = log(k(j)/1000d0)
expected = expected_khs(i, j)
Expand Down Expand Up @@ -189,7 +189,7 @@ subroutine test_kh_D2O(error)
T_C = T_K - 273.15d0

do i=1, ngas
call kh(T_C, gases(i), heavywater, k)
call kh(T_K, gases(i), heavywater, k)
do j=1, nT
value = log(k(j)/1000d0)
expected = expected_khs(i, j)
Expand Down Expand Up @@ -233,7 +233,7 @@ subroutine test_kd_H2O(error)
T_C = T_K - 273.15d0

do i=1, ngas
call kd(T_C, gases(i), heavywater, k)
call kd(T_K, gases(i), heavywater, k)
do j=1, nT
value = log(k(j)/1000d0)
expected = expected_khs(i, j)
Expand Down Expand Up @@ -270,7 +270,7 @@ subroutine test_kd_D2O(error)
T_C = T_K - 273.15d0

do i=1, ngas
call kd(T_C, gases(i), heavywater, k)
call kd(T_K, gases(i), heavywater, k)
do j=1, nT
value = log(k(j)/1000d0)
expected = expected_khs(i, j)
Expand Down

0 comments on commit 72b1cf9

Please sign in to comment.