/*
 * 
 *  RINGS OF SATURN 
 *        converted from BASIC to C 
 *              BASIC Code from Sky & Telescope, May 1995 
 *              Corrected lines 144 and 184 on November 26, 1996
 */
 
#include <stdio.h>
#include <math.h>
 
#define atn(d) atan(d)
 
float rd=M_PI/(double)180.;
float dl=9.;
float year;
 
float sr, sb, sl, er, eb, el;  /* saturn and earth positions */
float t;
 
main (int argc, char *argv[])
{
    int i;

    if (argc != 2) {
      printf("\n   This program will compute Saturn's ring inclination as");
      printf("\n seen from Earth (b) and as seen from the Sun (b') for");
      printf("\n each day of a given year.");
      printf("\n\n                  usage: satring year\n\n");
      exit(1); }
 
    year = atoi (argv[1]);
    do_mainloop();
}
 
do_mainloop()
{
    float a, b, jd, je, i, om, x, z, la, be;
    float s, sp, bp, xb=0., xp=0, al;
    float aj, bj, cj, dj, ej, d, m, ci, id, mv;
 
    a = floor ((year-1.)/100.);
    b = 2. - a + floor(a/4.);
    if (year<1583.) b=0.;
    jd = floor(365.25*(year+4715))+floor(30.6001*(14));
    jd += (b-1523.5);
    je = jd+365.;
 
    do {
	/* loop over days in the year */
 
	t = (jd-2451545.)/365250.;
	i = (28.04922-.13*t+.0004*t*t)*rd;
	om = (169.53+13.826*t+.04*t*t)*rd;
 
	get_planet_pos();
 
	x = sr*cos(sb)*cos(sl)-er*cos(el);
	year = sr*cos(sb)*sin(sl)-er*sin(el);
	z = sr*sin(sb)-er*sin(eb);
 
	dl = sqrt(x*x+year*year+z*z);
	la = atn(year/x);
	if (x<0) la+=M_PI;
 
	be = atn(z/sqrt(x*x+year*year));
	s = sin(i)*cos(be)*sin(la-om)-cos(i)*sin(be);
	b = atn(s/sqrt(1.-s*s));
	sp = sin(i)*cos(sb)*sin(sl-om)-cos(i)*sin(sb);
	bp = atn(sp/sqrt(1.-sp*sp));
 
	if (b*xb<0) {
	    printf ("< Earth crosses ring plane (%s) >\n",
		    ((b>0.) ? "s to n" : "n to s"));
	}
	if (bp*xp<0) {
	    printf ("< Sun crosses ring plane (%s) >\n",
		    ((bp>0.) ? "s to n" : "n to s"));
	}
 
	z = floor(jd+1);
	al = floor((z-1867216.25)/36524.25);
	aj = z+1.+al-floor(al/4.);
	if (z<2299161.) aj=z;
	bj = aj+1524;
	cj = floor((bj-122.1)/365.25);
	dj = floor(365.25*cj);
	ej = floor((bj-dj)/30.6001);
	d = bj-dj-floor(30.6001*ej);
	m = ej-1;
	if (ej>13.5) m=m-12.;
	year = cj-4715;
	if (m>2.) year=year-1.;
 
	printf ("%2d %2d %5d ", (int)m, (int)d, (int)year);
	printf ("  b =%5.2f deg", b/rd);
	printf ("  b'=%5.2f deg", bp/rd);
 
	ci = (sr*sr+dl*dl-er*er)/(2*sr*dl);
	id = atn(sqrt(1-ci*ci)/ci)/rd;
	mv = -8.88+5.*log(sr*dl)/log(10);
	mv = mv+.044*id-2.6*fabs(s)+1.25*s*s;
	printf ("  Mv =%3.1f\n", mv);
 
	jd = jd+1.;
	xb = b;
	xp = bp;
 
    } while (jd <= je);
}
 
get_planet_pos ()
{
    float l0, l1, l2, l3, l4;
    float r0, r1, r2, r3;
    float b0, b1, b2;
 
    l0 = 175347046.;
    l0 = l0+3341656.*cos(4.66926+6283.07585*t);
    l0 = l0+34894.*cos(4.6261+12566.1517*t);
    l1 = 628331966747.;
    l1 = l1+206059.*cos(2.67824+6283.07585*t);
    l2 = 52919.;
    el = (l0+l1*t+l2*t*t)/(1E+08);
    eb = 0;
 
    r0 = 100013989.;
    r0 = r0+1670700.*cos(3.09846+6283.07585*t);
    r0 = r0+13956.*cos(3.05525+12566.1517*t);
    r1 = 103019*cos(1.10749+6283.07585*t);
    r2 = 4359.*cos(5.7846+6283.0758*t);
    er = (r0+r1*t+r2*t*t)/(1e+08);
 
    t = t-(.00578*dl)/365250.;
    l0 = 87401354.;
    l0 = l0+11107660.*cos(3.96205+213.2991*t);
    l0 = l0+1414151.*cos(4.58582+7.11355*t);
    l0 = l0+398379.*cos(.52112+206.18555*t);
    l0 = l0+350769.*cos(3.3033+426.598191*t);
    l0 = l0+206816.*cos(.24658+103.09277*t);
    l0 = l0+79271.*cos(3.84007+220.41264*t);
    l0 = l0+23990*cos(4.66977+110.20632*t);
    l0 = l0+16574*cos(.43719+419.48464*t);
    l0 = l0+15820*cos(.93809+632.78374*t);
    l0 = l0+15054*cos(2.7167+639.89729*t);
    l0 = l0+14907*cos(5.76903+316.39187*t);
    l0 = l0+14610*cos(1.56519+3.93215*t);
    l0 = l0+13160*cos(4.44891+14.22709*t);
    l0 = l0+13005*cos(5.98119+11.0457*t);
    l0 = l0+10725*cos(3.1294+202.2534*t);
 
    l1 = 21354295596.;
    l1 = l1+1296855.*cos(1.82821+213.2991*t);
    l1 = l1+564348.*cos(2.885+7.11355*t);
    l1 = l1+107679.*cos(2.2777+206.18555*t);
    l1 = l1+98323.*cos(1.0807+426.59819*t);
    l1 = l1+40255.*cos(2.04128+220.41264*t);
    l2 = 116441.*cos(1.17988+7.11355*t);
    l2 = l2+91921.*cos(.07435+213.2991*t);
    l2 = l2+90592.;
    l2 = l2+15277*cos(4.06492+206.18555*t);
    l3 = 16039*cos(5.73945+7.11355*t);
    l4 = 1662*cos(3.9983+7.1135*t);
    sl = (l0+l1*t+l2*t*t+l3*t*t*t+l4*t*t*t*t)/(1E+08);
 
    b0 = 4330678.*cos(3.60284+213.2991*t);
    b0 = b0+240348.*cos(2.85238+426.59819*t);
    b0 = b0+84746.;
    b0 = b0+34116.*cos(.57297+206.18555*t);
    b0 = b0+30863*cos(3.48442+220.41264*t);
    b0 = b0+14734*cos(2.11847+639.89729*t);
    b0 = b0+9917*cos(5.79+419.4846*t);
    b0 = b0+6994*cos(4.736+7.1135*t);
    b1 = 397555.*cos(5.3329+213.2991*t);
    b1 = b1+49479.*cos(3.14159);
    b1 = b1+18572*cos(6.09919+426.59819*t);
    b1 = b1+14801*cos(2.30586+206.18555*t);
    b1 = b1+9644*cos(1.6967+220.4126*t);
    b2 = 20630*cos(.50482+213.2991*t);
    sb = (b0+b1*t+b2*t*t)/(1E08);
 
    r0 = 955758136.;
    r0 = r0+52921382.*cos(2.39226+213.2991*t);
    r0 = r0+1873680.*cos(5.2355+206.18555*t);
    r0 = r0+1464664.*cos(1.64763+426.59819*t);
    r0 = r0+821891.*cos(5.9352+316.39187*t);
    r0 = r0+547507.*cos(5.01533+103.09277*t);
    r0 = r0+371684.*cos(2.27115+220.41264*t);
    r0 = r0+361778.*cos(3.13904+7.11355*t);
    r1 = 6182981.*cos(.25844+213.2991*t);
    r1 = r1+506578.*cos(.71115+206.18555*t);
    r1 = r1+341394.*cos(5.79636+426.95819*t);
    r2 = 436902.*cos(4.78672+213.2991*t);
    r3 = 20315*cos(3.02187+213.2991*t);
    sr = (r0+r1*t+r2*t*t+r3*t*t*t)/(1e+08);
 
}