An algorithm from the Jean Meeus book-CodePudding

I am writing a simple program to calculate the phases of the moon; I use the algorithms presented in Jean Meeus' book "Astronomical algorithms" 2nd edition.

On p. 350 (chapter 49) he writes:

Calculate [...] the following angles, which are expressed in degrees and may be reduced to the interval 0-360 degrees and, if necessary, to radians before going further on

So the conversion should look like this:

reduction of the angle to the range 0-360
conversion of the angle from degrees to radians

Next, on page 353, he gives a series of examples:

M = -8234.2625 = 45.7375

M1 = -108984.6278 = 95.3722

But I don't understand how the conversions were made. I have tried the following calculations (in pascal-like) for the value M

// 1. reduction of the angle to the range 0-360
deg := -8234.2625;
reducedDeg := Abs(deg mod 360) // 314.2625

// 2. conversion of the angle from degrees to radians
convertedDeg := reducedDeg * Pi / 180 // 5.4849

So the result of the conversion is 5.4849 while the expected result should be 45.7375. Same thing for the second example M1:

// 1. reduction of the angle to the range 0-360
deg: = -108984.6278;
reducedDeg := Abs(deg mod 360); // 264.6278;

// 2. conversion of the angle from degrees to radians
convertedDeg := reducedDeg * Pi / 180; // 4,6186 should be 95.3722

What could be the problem?

For clarity:

Abs(val mod x) in those examples is the positive (absolute) remainder of dividing val by x; is an abbreviation for the following sequence:

// reduction of the angle to the range 0-360 => Abs(-8234.2625 mod 360)
val := -8234.2625;
val := val / 360; // −22.872951389
If (val < 0) // in the range 0-360 there are only positive numbers
begin
  val := val * -1; // 22.872951389
end
val := val - Int(val); // 0.872951389
val := val * 360; // 314.2625

CodePudding user response：

Your “mod” algorithm (reduce angle to [0.0°, 360°) range) is wrong, in particular:

If (val < 0) // in the range 0-360 there are only positive numbers
begin
  val := val * -1; // 22.872951389
end

This will mirror, “flip” the angle along the x-axis. You do not want to do that. What you meant to do is

val := val   ord(val < 0) * 360; { conditionally add a complete turn }

You can avoid all of this troubles if you have a compiler complying to the ISO standard 10206, “Extended Pascal”. You can then solve this sub‑task using the built‑in data type complex:

program degreeToRadianDemo(output);
const
    { This is meant to familiarize you with the functions. }
                pi = 2 * arg(cmplx(0.0, maxReal));
    { 1° ≘ π / 180 rad }
    degreeInRadian = pi / 180;
var
    degree, radian: real;
begin
    degree := -108984.6278;
    
    { `arg` returns principal argument, i. e. `real` in (−π,  π] }
    radian := arg(polar(maxReal, degree * degreeInRadian));
    { Add one turn (2π) if we’re in the (−π, 0.0) range. }
    radian := radian   ord(radian < 0.0) * 2 * pi;
    
    writeLn(                 radian:8:4);
    { Convert back to degrees. }
    writeLn(radian / degreeInRadian:8:4);
end.

This works with the GPC (GNU Pascal Compiler). The FPC (FreePascal Compiler) furnishes a ucomplex unit, but it only supports a rectangular complex (i. e. there is no arg and polar). [FPC version 3.2.0]

Furthermore – I know this is meant to be an exercise for you – but don’t program what’s already been programmed for you. Delphi and the FPC deliver a math unit. This unit has the DegToRad and RadToDeg functions. There is also a FMod function – mod operator for real numbers – and mod operator overload, because using the DegToRad/RadToDeg functions you will still need to eliminate complete turns. Have a look at FMod’s implementation.

Since there seems to have been some confusion regarding the mod operator, a quote from the ISO standard 7185 (“Standard Pascal”), page 48:

A term of the form i mod j shall be an error if j is zero or negative; otherwise, the value of i mod j shall be that value of (i − (k × j)) for integral k such that 0 ≤ i mod j < j.

Thus the result of the mod operator is guaranteed to be non-negative. Unfortunately, though, not all compilers adhere to the ISO standards. For instance, in the FPC (FreePascal Compiler) will only return the proper result if {$modeSwitch isoMod } is set. Rest assured, however, Delphi and the GPC (GNU Pascal Compiler) do work correctly.