In the May 95 Sky & Telescope, there was a program listing in BASIC which computes the angle from Saturn's ring plane to the Earth, to the Sun, and the apparent magnitude of Saturn as viewed from Earth.

I entered the program on my PC, and then converted it to C. Here it is in case anyone wants to use it.

The original BASIC code is commented out using the cpp construct "#ifdef BASIC". If you delete all the C code, you should have the original basic code left, if that's what you want. The BASIC code prompts for a year, then displays one line of output for each day of the year, pausing when it detects a crossing of the ring plane.

The C code takes the year off the command line and uses the year 0 if no year appears on the command line. It also prints one line of data for each day of the year but doesn't bother pausing when it detects a ring plane crossing.

The output for each day is of the form: month day year b= b'= Mv=

where "b" is the angle from Saturn's ring plane to the earth, "b'" the angle from the ring plane to the Sun, and "Mv" the apparent Magnitude of Saturn as viewed from the Earth.

I have taken this code and wrapped it with OpenGL calls which create a 3D model of Saturn and display it at the correct angle. If you have a system with an OpenGL library, mail me and I'll send you the source for this program.

Enjoy...

 
   -paul	pmartz@dsd.es.com    Evans & Sutherland
		"Be bright! Feel right! Take Eno!"
 
-------------------------------------------------------------
 
#ifdef BASIC
10 ' from Sky & Telescope, May 1995 */
#endif /* BASIC */
 
 
#include 
#include 
 
 
#define atn(d) atan(d)
 
#ifdef BASIC
20 defdbl a-z
30 cls: pi=4#*atn(1#): rd=pi/180#: dl=9#
#endif
 
float rd=M_PI/(double)180.;
float dl=9.;
float year;
 
float sr, sb, sl, er, eb, el;  /* saturn and earth positions */
float t;
 
 
#ifdef BASIC
40 input "Year of Interest "; y
#endif
 
main (int argc, char *argv[])
{
    int i;
 
    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;
 
 
#ifdef BASIC
50 a=int((y-1#)/100#): b=2#-a+int(a/4#)
60 if y<1583 then b=0
70 jd=int(365.25*(y+4715))+int(30.6001*(14))
80 jd=jd+b-1523.5
90 je=jd+365#
#endif
 
    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 */
 
#ifdef BASIC
100 t=(jd-2451545#)/365250#
110 i=(28.04922#-.13#*t+.0004#*t*t)*rd
120 om=(169.53#+13.826#*t+.04#*t*t)*rd
130 gosub 550: ' get planet positions
#endif
 
	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();
 
#ifdef BASIC
140 x=sr*cos(sb)*cos(sl)-er*cos(el)
150 y=sr*cos(sb)*sin(sl)-er*sin(el)
160 z=sr*sin(sb)-er*sin(eb)
#endif
 
	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);
 
#ifdef BASIC
170 dl=sqr(x*x+y*y+z*z)
180 la=atn(y/x): if x<0 then la=la+pi
#endif
 
	dl = sqrt(x*x+year*year+z*z);
	la = atn(year/x);
	if (x<0) la+=M_PI;
 
#ifdef BASIC
190 be=atn(z/sqr(x*x+y*y))
200 s=sin(i)*cos(be)*sin(la-om)-cos(i)*sin(be)
210 b=atn(s/sqr(1#-s*s))
220 sp=sin(i)*cos(sb)*sin(sl-om)-cos(i)*sin(sb)
230 bp=atn(sp/sqr(1#-sp*sp))
#endif
 
	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));
 
#ifdef BASIC
240 ' check for crossing of ring plane
250 if b*xb>= 0 then goto 280
260 n$="(n to s)": if b>0 then n$="(s to n)"
270 print "< Earth crosses ring plane "; n$; " >"
280 if bp*xp>=0 then goto 310
290 n$="(n to s)": if bp>0 then n$="(s to n)"
300 print "< Sun crosses ring plane "; n$; " >"
#endif
 
	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"));
	}
 
#ifdef BASIC
310 z=int(jd+1): al=int((z-1867216.25#)/36524.25#)
320 aj=z+1#+al-int(al/4#): if z<2299161# then aj=z
330 bj=aj+1524: cj=int((bj-122.1)/365.25)
340 dj=int(365.25#*cj): ej=int((bj-dj)/30.6001)
350 d=bj-dj-int(30.6001*ej)
360 m=ej-1: if ej>13.5 then m=m-12
370 y=cj-4715: if m>2 then y=y-1
#endif
 
	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.;
 
#ifdef BASIC
380 print using "## ## #####"; m; d; y;
390 print using "  b =###.## deg"; b/rd;
400 print using "  b'=###.## deg"; bp/rd;
#endif
 
	printf ("%2d %2d %5d ", (int)m, (int)d, (int)year);
	printf ("  b =%5.2f deg", b/rd);
	printf ("  b'=%5.2f deg", bp/rd);
 
#ifdef BASIC
440 '
450 ci=(sr^2+dl^2-er^2)/(2*sr*dl)
460 id=atn(sqr(1-ci*ci)/ci)/rd
470 mv=-8.88+5*log(sr*dl)/log(10)
480 mv=mv+.044*id-2.6*abs(s)+1.25*s*s
490 print using "  Mv =##.#"; mv
#endif
 
	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);
 
#ifdef BASIC
500 if b*xb<0 then while inkey$="": wend
510 if bp*xp<0 then while inkey$="": wend
 
520 jd=jd+1#: xb=b: xp=bp
#endif
 
	jd = jd+1.;
	xb = b;
	xp = bp;
 
#ifdef BASIC
530 if jd<=je then goto 100
540 end
#endif
 
    } while (jd <= je);
}
 
 
 
 
 
get_planet_pos ()
{
    float l0, l1, l2, l3, l4;
    float r0, r1, r2, r3;
    float b0, b1, b2;
 
#ifdef BASIC
550 ' calculate position of earth
560 l0=175347046#
570 l0=l0+3341656#*cos(4.66926+6283.07585#*t)
580 l0=l0+34894#*cos(4.6261+12566.1517#*t)
590 l1=628331966747#
600 l1=l1+206059#*cos(2.67824+6283.07585#*t)
610 l2=52919#
620 el=(l0+l1*t+l2*t*t)/(1E+08): eb=0
#endif
 
    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;
 
#ifdef BASIC
630 r0=100013989#
640 r0=r0+1670700#*cos(3.09846+6283.07585#*t)
650 r0=r0+13956#*cos(3.05525+12566.1517#*t)
660 r1=103019*cos(1.10749+6283.07585#*t)
670 r2=4359#*cos(5.7846+6283.0758#*t)
680 er=(r0+r1*t+r2*t*t)/(1e+08)
#endif
 
    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);
 
#ifdef BASIC
690 ' calculate position of saturn
700 t=t-(.00478#*dl)/365250#
710 l0=87401354#
720 l0=l0+11107660#*cos(3.96205+213.2991#*t)
730 l0=l0+1414151#*cos(4.58582+7.11355*t)
740 l0=l0+398379#*cos(.52112+206.18555#*t)
750 l0=l0+350769#*cos(3.3033+426.598191#*t)
760 l0=l0+206816#*cos(.24658+103.09277#*t)
770 l0=l0+79271#*cos(3.84007+220.41264#*t)
780 l0=l0+23990*cos(4.66977+110.20632#*t)
790 l0=l0+16574*cos(.43719+419.48464#*t)
800 l0=l0+15820*cos(.93809*632.78374#*t)
810 l0=l0+15054*cos(2.7167+639.89729#*t)
820 l0=l0+14907*cos(5.76903+316.39187#*t)
830 l0=l0+14610*cos(1.56519+3.93215*t)
840 l0=l0+13160*cos(4.44891+14.22709*t)
850 l0=l0+13005*cos(5.98119+11.0457*t)
860 l0=l0+10725*cos(3.1294+202.2534*t)
#endif
 
    t = t-(.00478*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);
 
#ifdef BASIC
870 l1=21354295596#
880 l1=l1+1296855#*cos(1.82821+213.2991#*t)
890 l1=l1+564348#*cos(2.885+7.11355*t)
900 l1=l1+107679#*cos(2.2777+206.18555#*t)
910 l1=l1+98323#*cos(1.0807+426.59819#*t)
920 l1=l1+40255#*cos(2.04128+220.41264#*t)
930 l2=116441#*cos(1.17988+7.11355*t)
940 l2=l2+91921#*cos(.07435+213.2991*t)
950 l2=l2+90592#
960 l2=l2+15277*cos(4.06492+206.18555#*t)
970 l3=16039*cos(5.73945+7.11355*t)
980 l4=1662*cos(3.9983+7.1135*t)
990 sl=(l0+l1*t+l2*t*t+l3*t*t*t+l4*t^4)/(1E+08)
#endif
 
    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);
 
#ifdef BASIC
1000 b0=4330678#*cos(3.60284+213.2991#*t)
1010 b0=b0+240348#*cos(2.85238+426.59819#*t)
1020 b0=b0+84746#
1030 b0=b0+34116#*cos(.57297+206.18555#*t)
1040 b0=b0+30863*cos(3.48442+220.41264#*t)
1050 b0=b0+14734*cos(2.11847+639.89729#*t)
1060 b0=b0+9917*cos(5.79+419.4846*t)
1070 b0=b0+6994*cos(4.736+7.1135*t)
1080 b1=397555#*cos(5.3329+213.2991#*t)
1090 b1=b1+49479#*cos(3.14159)
1100 b1=b1+18572*cos(6.09919+426.59819#*t)
1110 b1=b1+14801*cos(2.30586+206.18555#*t)
1120 b1=b1+9644*cos(1.6967+220.4126*t)
1130 b2=20630*cos(.50482+213.2991*t)
1140 sb=(b0+b1*t+b2*t*t)/(1E08)
#endif
 
    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);
 
#ifdef BASIC
1150 r0=955758136#
1160 r0=r0+52921382#*cos(2.39226+213.2991#*t)
1170 r0=r0+1873680#*cos(5.2355*206.18555#*t)
1180 r0=r0+1464664#*cos(1.64763+426.59819#*t)
1190 r0=r0+821891#*cos(5.9352+316.39187#*t)
1200 r0=r0+547507#*cos(5.01533+103.09277#*t)
1210 r0=r0+371684#*cos(2.27115+220.41264#*t)
1220 r0=r0+361778#*cos(3.13904+7.11355*t)
1230 r1=6182981#*cos(.25844+213.2991#*t)
1240 r1=r1+506578#*cos(.71115+206.18555#*t)
1250 r1=r1+341394#*cos(5.79636+426.95819#*t)
1260 r2=436902#*cos(4.78672+213.2991#*t)
1270 r3=20315*cos(3.02187+213.2991#*t)
1280 sr=(r0+r1*t+r2*t*t+r3*t*t*t)/(1e+08)
#endif
 
    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);
 
/* 1290 return */
 
}