This is my first post, and my very first contribution will really be a simple one...
1.Background
Knowledge of the water vapour saturation pressure is paramount in applied meteorology, as it allows converting measurements of relative humidity to absolute water content.
For instance, if E and ES represent respectively the partial pressure of water vapour in air at the reference pressure of 1025.15 hPa, and ES the water vapour saturation pressure, then by definition the relative humidity is
E
Rh = 100 -----
ES
The water saturation pressure depends on temperature, and many approximate formulae have been proposed to fit measurements. One of the most popular, and simple, is the "Magnus formula",
17.1|Ta|
ES = 6.1078 exp ------------- (hPa)
235 + |Ta|
where Ta is air temperature in °C.
This is the formula I've used in the "pbl_met" library ( https://github.com/serv-terr/pbl_met ). Yet there is something better, namely the Goff and Gratch formula, reformulated by Murray:
ES = 7.95357242e10 exp U
where
U = -18.1972839 Q + 5.02808 ln Q - 70242.1852 exp( -26.1205253 / Q) + 58.0691913 exp( -8.03945282 Q )
with
Q = 373.16 / T
In the later expression, T is the absolute temperature in K, and the resulting ES is again in hPa.
One may wonder, why should we use a "better" formula when the popular and the new are both approximations. One reason (at least my one) is, in this current moment, purely aesthetical: in lack of a theory allowing to get an "exact" formula, I don't see any special reason why we should not use a "more accurate version". It will not be perfect, but better, and that's a legitimate motivation...
Sure, the Magnus formula is visually simpler. This is a very nice property, would you have to perform calculations with a slide rule. But to date we can count on more advanced computing devices, like the DM42, and "simplicity" should not be an issue.
To this we might add another point: the Goff and Gratch (Murray) formula is surely more "complicate", yet not really more computationally "complex". Both formulae are algebraic, and their evaluation on a vector of data of size N is anyway O(N), that is, "linear in the problem size".
2.Implementation
The reason I desired an RPN implementation of the Goff and Gratch (Murray) formula is I need to check its numerical stability, in view of inclusion into the pbl_met library. That the formula might be numerically unstable is not completely a theoretical issue: its main exponential factor is the difference between positive quantities, and would their magnitude be similar in some point of the meteorological range ( -40 to +60 °C ) digit cancellation might occur.
I know the DM42, as its HP42s predecessor, allows use of advanced RPN instructions, like variables and user menus. But in view of testing, I've chosen to make the implementation as simple as possible, so to be compatible with my DM41X. The plan is comparing DM42 and DM41X results on same input data, instead of conducing an extensive numerical analysis just right now. You may think to this as the mathematical equivalent of "speditive survey" in geology: get a first idea, and understand whether some need exists to get deeper.
Because of that, my coding is (I admit) a bit "brutal":
Code: Select all
01 LBL "ES"
02 1/X
03 373.16
04 x
05 STO 00
06 -18.1972839
07 RCL 00
08 x
09 5.02808
10 RCL 00
11 LN
12 x
13 +
14 -26.1205253
15 RCL 00
16 ÷
17 E^X
18 70242.1852
19 x
20 -
21 -8.03945282
22 RCL 00
23 x
24 E^X
25 58.0691913
26 x
27 +
28 E^X
29 7.95357242E10
30 x
31 RTN
3.Rationale
The idea of comparing results of the "same" program on the DM42 and DM41X to check for delicate points of ES in the meteorological interval (-40 to +60 °C) may seem a bit crank, and sure it is.
Yet, I feel it is not completely irrational.
Internally, the DM41X and DM42 operate on very different precisions.
You may check this by yourself, by estimating the "machine epsilon" of both calculators. The "machine epsilon" is defined as the smallest floating point number greater than 0, let's say it "EPS", so that 1 + EPS ≠ 1.
You may estimate it with a program like this:
Code: Select all
01 LBL "EPS"
02 1
03 STO 02
04 0.1
05 STO 01
06 LBL 00
07 1
08 STO+ 02
09 10
10 STO÷ 01
11 RCL 01
12 1
13 +
14 1
15 X≠Y?
16 GTO 00
17 RCL 01
18 RCL 02
19 RTN
And here are my results:
DM42: Counter = 34 EPS = 1.000e-34
DM41X: Counter = 10 EPS = 1.000e-10
The smaller the machine epsilon, the more precise a data type (or in our case computer) is. Our DM42 has a precision roughly "three times and half" that of DM41X, and may be used, assuming a calibration mind set, as a sort of "reference instrument". The DM41X, then, can be thought as roughly simulating a conventional computer. Would any difference between the DM42 and DM41X results be found on same input data, and we would have found one numerically delicate point in the formula under test, in this case the one for ES.
In the moment, that's all.
Sorry so sloppy.
Patrizia