//  NYD.C
    //  Author: Peter Meyer
    //  This program calculates the statistics for NYD occurrences
    //  (dates in March/April) over a 4000-year period.
    //  Last mod.: 1999-03-29
    
    #include <stdio.h>
    
    #include <mplc.h>
    
    #define N 15
    
    long freq[64][3];
    int n[3] = { 201, 1001, 4001 };
    double sum_pc[3] = { 0.0, 0.0, 0.0 };
    
    //  Common Era    ------ Meyer-Palmen Solilunar Calendar ------     Julian day
    //    0-03-22    069-06-01-01 Wednesday, Aristarchus  1, 069-06      1,721,141
    
    /*-----------*/
    void main(void)
    {
    int i, j, diy, month, day;
    long year, ce_yr;
    double pc;
    MPDate mpdt;
    
    for ( i=0; i<64; i++ )
        for ( j=0; j<3; j++ )
            freq[i][j] = 0;
    
    mpdt.jdn = 1721141L;
    //  This is the JD of the MP NYD in Year 0 CE.
    jdn_to_mpdate(&mpdt);
    //  mpdt now holds the MP era date correspondng to the MP NYD in Year 0 CE.
    
    //  For each year in the range 0 through 4000 CE
    //  get the date CE of the MP NYD.
    for ( ce_yr=0; ce_yr<=4000; ce_yr++ )
        {
        //  Convert MP era date to CE date.
        mpdate_to_ce_date(&mpdt,&year,&month,&day);
    
        diy = day + 31*(month==4);
        //  diy = days since last day of February.
    
        freq[diy][2]++;
        if ( ce_yr >=1500 && ce_yr <= 2500 )
            freq[diy][1]++;
        if ( ce_yr >= 1900 && ce_yr <= 2100 )
            freq[diy][0]++;
    
        //  Advance mpdt to next NYD.
        if ( ++mpdt.year > 6840 )
            {
            mpdt.year = 1;
            mpdt.cycle++;
            }
        //  mpdt.month and mpdt.day both remain 1.
        }
    
    //  Now print the results.
    
    printf("\n%*s%*s%*s%*s",
        N,"",N,"2000CE+/-100",N,"2000CE+/-500",N,"2000CE+/-2000");
    printf("\n%*s%*s%*s%*s\n\n",
        N,"",N,"201 NYDs",N,"1001 NYDs",N,"4001 NYDs");
    
    for ( diy=1; diy<=61; diy++ )
        {
        if ( freq[diy][2] )
            {
            printf("%*s",N-3,( diy <= 31 ? "March" : "April" ));
            printf("%3d", ( diy <= 31 ? diy : diy-31 ));
            for ( i=0; i<3; i++ )
                {
                pc = (((double)freq[diy][i])*100)/n[i];
                sum_pc[i] += pc;
                printf("%*ld, %4.*f%%",N-7,freq[diy][i],(i?2:1),pc);
                }
            printf("\n");
            }
        }
    
    printf("%*s%*s%*s%*s\n",
        N,"",N,"-----",N,"-----",N,"-----");
    
    printf("%*s",N,"");
    for ( i=0; i<3; i++ )
        printf("%*.2f%%",N-1,sum_pc[i]);
    printf("\n");
    }