1 #include <math.h>
2
3 #include "astro.h"
4
5 #undef sqr
6 #define sqr(x) ((x)*(x))
7
8 /* given a planet, the sun, the planet's eq pole position and a
9 * position of a satellite (as eq x=+e y=+s z=front in planet radii) find x,y
10 * position of shadow.
11 * return 0 if ok else -1 if shadow not on planet
12 */
13 int
plshadow(Obj * op,Obj * sop,double polera,double poledec,double x,double y,double z,float * sxp,float * syp)14 plshadow (Obj *op, Obj *sop, double polera, double poledec, double x,
15 double y, double z, float *sxp, float *syp)
16 {
17 /* equatorial to ecliptic sky-plane rotation */
18 double sa = cos(op->s_dec) * cos(poledec) *
19 (cos(op->s_ra)*sin(polera) - sin(op->s_ra)*cos(polera));
20 double ca = sqrt (1.0 - sa*sa);
21
22 /* rotate moon from equatorial to ecliptic */
23 double ex = x*ca + y*sa;
24 double ey = -x*sa + y*ca;
25
26 /* find angle subtended by earth-sun from planet */
27 double a = asin (sin(op->s_hlong - sop->s_hlong)/op->s_edist);
28 double b = asin (-sin(op->s_hlat)/op->s_edist);
29
30 /* find displacement in sky plane */
31 double x0 = ex - z*tan(a);
32 double y0 = ey - z*tan(b);
33
34 /* projection onto unit sphere */
35 double x1 = x0 + (ex-x0)/sqrt(sqr(ex-x0)+sqr(z));
36 double y1 = y0 + (ey-y0)/sqrt(sqr(ey-y0)+sqr(z));
37
38 /* check behind or off edge */
39 if (z < 0 || sqr(x1) + sqr(y1) > 1)
40 return (-1);
41
42 /* rotate back to equatorial */
43 *sxp = x1*ca - y1*sa;
44 *syp = x1*sa + y1*ca;
45
46 return (0);
47 }
48
49