Ostermädchen

Calculating Easter in MATLAB

Every year, Easter Sunday lies in a different week. But which one?

The answer is as simple as profound: Determine the Sunday following the paschal full moon, the full moon. (Wikipedia, License CC-BY-SA)

This can be calculated thanks to Carl Friedrich Gauß. (Wikipedia).

My means of calculation is MATLAB, therefore here the calculation of Easter Sunday according to the refined Gauß method as MATLAB Code:

function eastersunday = easter(year)   
    K = fix(year/100);                  % 1.  secular number
    M = 15 + fix((3*K + 3)/4) - fix((8*K + 13)/25); % 2. secular moon shift
    S = 2- fix((3*K + 3)/4);            % 3.  secular sun shift
    A = mod(year,19);                   % 4.  moon parameter
    D = mod((19*A + M),30);             % 5.  seed for 1st full moon in spring
    R = fix((D + fix(A/11))/29);        % 6.  calendarian correction quantity
    OG = 21 + D - R;                    % 7.  easter limit
    SZ = 7 - mod((year + fix(year/4) + S),7); % 8.  first sunday in March
    OE = 7 - mod((OG - SZ),7);          % 9.  distance easter sunday from easter limit
    OS = OG + OE;                       % 10. easter sunday as March date
    eastersunday = datetime(year,3,1) + OS - 1; % 11. easter sunday as date
    eastersunday.Format = 'dd.MM.yyyy';
end

Easter Sunday 2015 I determine with

eastersun = easter(2015)
eastersun = 

   05.04.2015

If I want to know the dates of all Easter Sundays from years 1000 to 2500, I call this function as follows:

eastersun = ostern(1000:2500);

To determine which days occur most often, we just need month and date in the result. The datetime class allows easy calculation and manipulation with dates.

The year we set to 0 in order for the sorting to work:

eastersun.Format='dd.MM.';
eastersun.Year=0;

The histogram reflects the count of Easter Sundays for a given date:

hist(categorical(eastersun));

Some beautifications of this histogram:

distri = gca;
distri.XTickLabelRotation = 60;
distri.YGrid = 'on';
distri.YMinorGrid = 'on';
distri.Parent.Position(3)=800;

Now it’s your turn. Play around with the algorithm, and with dates in general.

Photo: Joachim Schlosser

Sharing is caring

Kommentare

Leave a Reply

Your email address will not be published. Required fields are marked *