From: lew@cbnewsd.cb.att.com (lewis.h.mammel..jr)
Newsgroups: sci.astro
Subject: Re: Help.Analemma question???
Summary: C program and output
Date: 30 Sep 92 22:45:37 GMT
Organization: AT&T


/*
Here is a program which prints out the equation of time.
It separately prints the contribution due to the obliquity
of the eclipitic ( a geometric effect ) and that due to
the eccentricity of the earth's orbit ( a dynamical effect ).

Note that the iteration is over evenly spaced positions of the
sun, rather than over evenly spaced times. This allows for
solving Kepler's equation directly, rather than inverting it.
This means that the values in column one of the output ( clock time )
"float". Zero clock time occurs when the dynamic mean sun
crosses the First Point of Aries. The real vernal equinox
is at a negative clock time.

The above is just by way of explaining how the output
can stay exact while retaining some sort of simplicity.

I think the last three columns of the output give the idea
of how the analemma attains the shape that it has. The second
and last columns could be used to plot an analemma.

ref: Marion's CLASSICAL DYNAMICS ( for Kepler's equation )

*/
#include <math.h>
#define PI 3.1415926535
main(argc,argv)
int argc;
char **argv;
{
        double psi,m,alpha,theta,delta;
        double Psi,beta,gamma;
        int i;

        /* obliquity of the ecliptic is 23.5 degrees */
        alpha = cos(23.5/180*PI);

        /* heliocentric longitude of perihelion is 102 degrees -
                convert to geocentric longitude */

        delta = (102.0-180.0)/180*PI;
        beta = sqrt((1.0 - 0.0167)/(1.0 + 0.0167));
        gamma = sin(23.5/180*PI);
        for( i= -24 ; i<24 ; i++ ) {

                /* psi is the ecliptic longitude of the sun */

                psi = PI/24*i;

                /* theta is the Right Ascension of the sun */

                theta = atan2(alpha*sin(psi),cos(psi));

                /* Psi is the eccentric anomaly for Kepler's equation */
                Psi = 2*atan2(beta*sin((psi-delta)/2),cos((psi-delta)/2));

                /* m is mean anomaly, or clock time.
                        This is Kepler's equation */

                m = Psi - 0.0167*sin(Psi);

                /* Displace zero back to equinox from perigee */
                m += delta;

                printf("%12.4f %12.4f %12.4f %12.4f %12.4f\n",

        m*365.25/PI/2, /* clock time in days from MEAN vernal equinox. */
        asin(gamma*sin(theta))*180.0/PI, /* declination of sun, in degrees. */
        (psi-theta)*24*30/PI, /* geometrical term in minutes */
        (m-psi)*24*30/PI, /* dynamic term in minutes */
        (m-theta)*24*30/PI); /* EQUATION OF TIME in minutes */

        }
}

   -180.7210       0.0000       0.0000       7.5067       7.5067
   -173.0782      -2.7395       2.4633       7.6383      10.1015
   -165.4687      -5.4601       4.7780       7.6389      12.4168
   -157.8922      -8.1403       6.8011       7.5091      14.3102
   -150.3481     -10.7535       8.4014       7.2517      15.6531
   -142.8351     -13.2661       9.4661       6.8717      16.3378
   -135.3515     -15.6350       9.9092       6.3760      16.2853
   -127.8949     -17.8073       9.6802       5.7735      15.4538
   -120.4628     -19.7212       8.7726       5.0746      13.8472
   -113.0522     -21.3091       7.2301       4.2913      11.5214
   -105.6595     -22.5039       5.1498       3.4367       8.5866
    -98.2813     -23.2474       2.6781       2.5253       5.2034
    -90.9137     -23.5000       0.0000       1.5722       1.5722
    -83.5527     -23.2474      -2.6781       0.5932      -2.0849
    -76.1941     -22.5039      -5.1498      -0.3957      -5.5455
    -68.8339     -21.3091      -7.2301      -1.3780      -8.6081
    -61.4679     -19.7212      -8.7726      -2.3375     -11.1101
    -54.0921     -17.8073      -9.6802      -3.2585     -12.9388
    -46.7027     -15.6350      -9.9092      -4.1256     -14.0348
    -39.2959     -13.2661      -9.4661      -4.9243     -14.3904
    -31.8683     -10.7535      -8.4014      -5.6411     -14.0425
    -24.4169      -8.1403      -6.8011      -6.2638     -13.0649
    -16.9389      -5.4601      -4.7780      -6.7816     -11.5595
     -9.4319      -2.7395      -2.4633      -7.1853      -9.6486
     -1.8941       0.0000       0.0000      -7.4677      -7.4677
      5.6758       2.7395       2.4633      -7.6233      -5.1600
     13.2786       5.4601       4.7780      -7.6489      -2.8709
     20.9148       8.1403       6.8011      -7.5434      -0.7423
     28.5838      10.7535       8.4014      -7.3080       1.0934
     36.2850      13.2661       9.4661      -6.9462       2.5199
     44.0168      15.6350       9.9092      -6.4636       3.4456
     51.7772      17.8073       9.6802      -5.8682       3.8120
     59.5637      19.7212       8.7726      -5.1700       3.6026
     67.3732      21.3091       7.2301      -4.3808       2.8493
     75.2024      22.5039       5.1498      -3.5143       1.6355
     83.0473      23.2474       2.6781      -2.5857       0.0924
     90.9038      23.5000       0.0000      -1.6112      -1.6112
     98.7676      23.2474      -2.6781      -0.6082      -3.2863
    106.6342      22.5039      -5.1498       0.4057      -4.7441
    114.4989      21.3091      -7.2301       1.4123      -5.8178
    122.3572      19.7212      -8.7726       2.3939      -6.3787
    130.2048      17.8073      -9.6802       3.3330      -6.3472
    138.0374      15.6350      -9.9092       4.2132      -5.6960
    145.8512      13.2661      -9.4661       5.0190      -4.4471
    153.6425      10.7535      -8.4014       5.7364      -2.6650
    161.4084       8.1403      -6.8011       6.3533      -0.4478
    169.1460       5.4601      -4.7780       6.8591       2.0812
    176.8534       2.7395      -2.4633       7.2456       4.7824