OS/2 e-Zine! - The Rexx File

The Rexx File

- by Dr. Dirk Terrell

Dr. Dirk Terrell is an astronomer at the University of Florida specializing in interacting binary stars. His hobbies include cave diving, martial arts, painting and writing OS/2 software such as HTML Wizard.

Related Articles
Recent Rexx Files:
Vol 3, Nr. 12
Vol 3, Nr. 9
Vol 3, Nr. 6
Vol 3, Nr. 4

Blast Back! Send a private message directly to Dr. Dirk Terrell with your thoughts:

Contact OS/2 e-Zine!

Go to a Printer Friendly version of this page

Summary: A look at advanced math functions in Rexx

Being an interpreted language, REXX is usually not a good choice when it comes to complex numerical calculations. There are situations, however, where REXX is a good choice. Fortran is normally the language of choice for scientists like myself because of its optimizations for numerical work. But Fortran is not very useful when it comes to reading and writing text files. REXX, of course, excels at handling text and piecing things together. And with a good mathematical library, it can do a decent job for numerical work.

This month a reader at NASA's Marshall Space Flight Center asked me if I had any OS/2 REXX code to calculate the position of the Sun at a given time and location. They had a program that would do such calculations, but it wasn't very convenient to use. With a REXX program, users could more easily customize it for their needs. I hadn't written such a program in REXX, but I thought it would give me a chance to talk about the numerical capabilities of REXX.

REXX comes with the basic mathematical operators:

Operator Symbol Example
Addition + 4 + 2 = 6
Subtraction - 4 - 2 = 2
Multiplication * 4 * 2 = 8
Division / 4 / 2 = 2
Exponentiation ** 4 ** 2 = 16

You also have a couple of other useful operators at you disposal. The // operator returns the remainder of a division. For example, nine divided by two doesn't yield a whole number result, but two plus a remainder of one. The // operator would return one in that case. The other useful operator, which is complementary to the // operator is the % operator which returns only the whole number part of a division. So, 9 % 2 would return four.

And that's where the numerical functions of REXX end unfortunately. For anything else, you will have to either write a REXX subroutine to perform the function or use an external library. For example, numerical work will involve the trigonometric functions like sine and cosine. These functions can be written as the sums of infinite series. The cosine function, for example, can be written as:

cosine(x) = 1 - (x**2)/2! + (x**4)/4! - (x**6)/6! - ...

Where the ! means the factorial function, e.g. 6! = 6 * 5 * 4 * 3 * 2 * 1. So, a REXX routine that implements this algorithm looks like this:


/* Routine to calculate the cosine of an angle */
Say "Enter the angle in radians:"
Parse Pull Angle
Say "How many digits of precision?"
Parse Pull NDigits
Numeric Digits NDigits+1
Limit=10**-(NDigits+1)
C=cos(Angle)
Say "Cos("||Angle||") =" C

Exit

Cos:
/* Compute the cosine of an angle in radians */
Procedure Expose Limit
Parse Arg Angle
Sum=1
Sign=1
Epsilon=1
Old=0
i=0
Do While Epsilon > Limit
i=i+2
Sign=-1*Sign
Sum=Sum+Sign*(Angle**i/Factorial(i))
Say "Sum:" Sum
Epsilon=Old-Sum
/* Make epsilon always positive for comparison */
If Epsilon < 0  then
Epsilon=-1*Epsilon
Old=Sum
end /* do */

Return Sum

and the factorial routine looks like this:

Factorial:
Procedure
Parse Arg N
Prod=1
Do j=N to 2 By -1
Prod=Prod*j
end /* do */
Return Prod

This program asks the user for the angle in radians and the number of digits to use in the calculations. That is one advantage to using REXX for numerical calculations: the ability to easily change the precision of the calculation. Now, above I said that the series we are using is infinite, but obviously we are not doing an infinite number of calculations. How many do we need to do? The answer lies in the program's use of the number of digits of precision. We only need to do the calculations until the computed number doesn't change from one iteration to the next. So, we ask the user how many digits of precision are necessary and then compute the cosine series until the difference between one iteration and the next is less than the number of significant digits dictates.

The above routine can be used to calculate the cosine of an angle, but it may not be the optimal solution if you have a large number of cosines to calculate because of the interpreted nature of REXX. If you need a faster method or don't want to write the routines yourself, then you need to use an external library for math functions. There are several of them around and I use REXXMATH by Zhitao Zeng.

REXXMATH has many mathematical functions. Some of the common ones are:

Function Description
acos( x ) Inverse cosine
asin( x ) Inverse sine
atan( x ) Inverse tangent (-Pi/2 to +Pi/2 version)
atan2( x, y ) Inverse tangent (0 to 2*Pi version)
ceil( x ) Next integer larger than x
cos( x ) Cosine
cosh( x ) Hyperbolic cosine
exp( x ) Exponential function, e.g. e**x)
fabs( x ) Absolute value
floor( x ) Next integer less than x
fmod( x, y ) Modulo function (same as // operator in REXX)
log( x ) Natural (base e) logarithms
log10( x ) Common (base 10) logarithms
sin( x ) Sine
sinh( x ) Hyperbolic sine
sqrt( x ) Square root
tan( x ) Tangent
tanh( x ) Hyperbolic tangent

To use these functions, simply load the library in the usual way:

call RxFuncAdd 'MathLoadFuncs', 'REXXMATH', 'MathLoadFuncs'
call MathLoadFuncs

and afterwards they will be available to you. Remember that the argument to a function like the cosine is an angle and that the angle must be measured in radians, not degrees. Most people are used to thinking of angles in degrees, so you must convert angles measured in degrees to radians before passing them to these functions. The conversion is pretty simple though. There are 2*Pi radians in a full circle, 360 degrees. So, to convert an angle X in degrees to radians, you multiply it by 360/(2*Pi) or 180/Pi. To convert from radians to degrees, just invert the conversion factor. In REXX, two routines that do these conversions look like this:

Rad2Deg:
Procedure Expose Pi
/* Convert an input angle in radians to degrees */
Parse Arg Radians
Radians = Radians * (180/Pi)
Return Radians

and

Deg2Rad:
Procedure Expose Pi
/* Convert an input angle in degrees to radians */
Parse Arg Degrees
Radians = Degrees * (Pi/180)
Return Radians

Obviously, a library like REXXMATH saves you a lot of programming time as opposed to writing everything yourself in REXX.

Finally, let me make a few comments on the code to calculate the position of the Sun at a given location and time. I won't bore you with the details of how such a calculation is done since it is, as you might imagine, somewhat complicated and not really germane to this column. But I will make it available for the curious. Note that this code is only accurate to a few minutes and doesn't include corrections for various effects that would be necessary for higher accuracy. It computes the altitude and azimuth of the Sun given the date and time and location of the observer. The date and time must be in Universal Time, so make sure to apply the correction for your time zone before entering it. The altitude of an object is the angle that it makes to the horizon, 0 degrees meaning that the object is on the horizon and 90 degrees meaning that it is directly overhead. The azimuth is an indication of the compass direction the observer is facing, with 0 degrees being north, 90 degrees east, 180 degrees south, and 270 degrees west. Have fun with it.

Download the source (.ZIP, 5K) for the examples in this week's Rexx File.

ISSN 1203-5696