The inputs are:
	The geocentric position of the sat.
	The velocity of the sat.
	The topocentric position of the sat.
	The RA and Dec of the object to be reflected.  Usually the sun
		but if you want to use the moon be sure to use the
		position as seen by the satellite.
	The two mirror rotation angles.

FUNCTION SAT_REFLECTION(
		EQ_X, EQ_Y, EQ_Z,
		V_X, V_Y, V_Z,
                SAT_X_TOPO_EQ, SAT_Y_TOPO_EQ, SAT_Z_TOPO_EQ,
                SUN_ALPHA, SUN_DELTA,
		ROT_1, ROT_2 : REAL) : REAL;
VAR
    SIN_ROT_1, COS_ROT_1, SIN_ROT_2, COS_ROT_2 : REAL;
    RR: ARRAY [1..3, 1..3] OF REAL;
    XX, YY, ZZ, TT, NTT : ARRAY [1..3] OF REAL;
    SUN_REF_X, SUN_REF_Y, SUN_REF_Z, SUN_REF_ALPHA, SUN_REF_DELTA : REAL;
    SUN_REF_ANG : REAL;

BEGIN
     { Calculate where the reflection of the observers gaze would point to. }
     { Assume that a mirror is attached to the sat and is pointing }
     { into the direction of motion.  It is now swiveled in the }
     { orbit plane (like an altitude axis, towards earth is negative) and }
     { then again about the (almost) earth-sat direction (like an azimuth }
     { axis, a positive value would point the mirror West if the sat was }
     { moving North).  By "almost" I really mean about the in-plane vector }
     { normal to motion). }

     { Step 1 - Create a new coordinate system oriented to the satellite's }
     {		local situation.  The vectors for x, y and z in the }
     {		equatorial system are : }
     {             x is the velocity vector }
     {		   y is the position-vector cross x }
     {		   z is x cross y }
     {		Normalize vectors. }
     {		Create the transform matrix RR using direction cosines. }
     {		(see Fund. of Astrod., p. 82). }
     XX[1] := V_X;
     XX[2] := V_Y;
     XX[3] := V_Z;
     TEMP := SQRT(SQR(XX[1]) + SQR(XX[2]) + SQR(XX[3]));
     XX[1] := V_X / TEMP;
     XX[2] := V_Y / TEMP;
     XX[3] := V_Z / TEMP;

     YY[1] := EQ_Y * V_Z - EQ_Z * V_Y;
     YY[2] := EQ_Z * V_X - EQ_X * V_Z;
     YY[3] := EQ_X * V_Y - EQ_Y * V_X;
     TEMP := SQRT(SQR(YY[1]) + SQR(YY[2]) + SQR(YY[3]));
     YY[1] := YY[1] / TEMP;
     YY[2] := YY[2] / TEMP;
     YY[3] := YY[3] / TEMP;

     ZZ[1] := XX[2] * YY[3] - XX[3] * YY[2];
     ZZ[2] := XX[3] * YY[1] - XX[1] * YY[3];
     ZZ[3] := XX[1] * YY[2] - XX[2] * YY[1];
     TEMP := SQRT(SQR(ZZ[1]) + SQR(ZZ[2]) + SQR(ZZ[3]));
     ZZ[1] := ZZ[1] / TEMP;
     ZZ[2] := ZZ[2] / TEMP;
     ZZ[3] := ZZ[3] / TEMP;

{DEFINE FLARE }  { This stuff is for testing only. }
{$IFDEF FLARE }
	{ Let's fake the orientation of the sat. }
	{ The sat is moving towards RA=90 degrees. }
	{ The sat's left side is pointing to DEC=90 degrees. }
	{ The sat's vertical (up) axis is pointing to RA=0 degrees. }
	{ Later I'll place the sun at RA=0 too so the sun will shine on the sat's top. }
	{ The mirrors our on the sat's sides.  We will have to rotate them. }
XX[1] := 0;
XX[2] := 1;
XX[3] := 0;
YY[1] := 0;
YY[2] := 0;
YY[3] := 1;
ZZ[1] := 1;
ZZ[2] := 0;
ZZ[3] := 0;
{$ENDIF }

	{ The xx, yy, zz vectors are mutually perpendicular and define the coordinate }
	{ system of the sat. }
	{   xx is forward (velocity vector) }
	{   yy is left (normal to the orbit plane) }
	{   zz is up (well, not up, but perpendicular to the other two; almost 'up') }
	{ Everything that we know about the MMA's is given with respect to this }
	{ coordinate system. Since the vectors are equatorial, we now know the }
	{ orientation of the sat in equatorial coordinates. The rr matrix (below) can }
	{ transform a vector in sat coordinates into equatorial coordinates. It }
	{ should be able to transform (1, 0, 0) back into your original velocity }
	{ vector (but normalized). I hope I got all that right. }

     RR[1, 1] := XX[1];
     RR[1, 2] := YY[1];
     RR[1, 3] := ZZ[1];
     RR[2, 1] := XX[2];
     RR[2, 2] := YY[2];
     RR[2, 3] := ZZ[2];
     RR[3, 1] := XX[3];
     RR[3, 2] := YY[3];
     RR[3, 3] := ZZ[3];

     { Step 2 - Initialize an x unit vector to (1, 0, 0). }
     {          Rotate it about the y-axis the desired amount. }
     {          Rotate it about the z-axis the desired amount. }
     {		Transform it to equatorial using RR and call it XX. }
     {		Repeat for y (0, 1, 0) and z (0, 0, 1). }
     {		These vectors define the coordinate system for the mirror. }
     {		The XX vector is normal to the surface of the mirror. }
     {		Recreate the transform matrix RR. }

     SIN_ROT_1 := SIN_D(ROT_1);	{ Rotate the y-axis. }
     COS_ROT_1 := COS_D(ROT_1);
     SIN_ROT_2 := SIN_D(ROT_2);	{ Rotate the z-axis. }
     COS_ROT_2 := COS_D(ROT_2);

{$IFDEF FLARE }
	{ Let's deploy a mirror 90 degrees.  That's backwards from what an Iridium }
	{ would do (-40).  So sunlight coming from above the sat would be reflected back up. }
	{ Since the first rotation is 90 the second doesn't change the orientation of }
	{ the plane of the mirror. }
     SIN_ROT_1 := SIN_D(90);
     COS_ROT_1 := COS_D(90);
     SIN_ROT_2 := SIN_D(45);
     COS_ROT_2 := COS_D(45);
{$ENDIF }

     TT[1] := 1;
     TT[2] := 0;
     TT[3] := 0;
     NTT[1] := TT[1] * COS_ROT_1 -  TT[3]  * SIN_ROT_1;
     NTT[2] := TT[2];
     NTT[3] := TT[1] * SIN_ROT_1 +  TT[3]  * COS_ROT_1;
     TT[1] := NTT[1] * COS_ROT_2 +  NTT[2] * SIN_ROT_2;
     TT[2] :=-NTT[1] * SIN_ROT_2 +  NTT[2] * COS_ROT_2;
     TT[3] := NTT[3];
     XX[1] := TT[1] * RR[1, 1] + TT[2] * RR[1, 2] + TT[3] * RR[1, 3];
     XX[2] := TT[1] * RR[2, 1] + TT[2] * RR[2, 2] + TT[3] * RR[2, 3];
     XX[3] := TT[1] * RR[3, 1] + TT[2] * RR[3, 2] + TT[3] * RR[3, 3];

     TT[1] := 0;
     TT[2] := 1;
     TT[3] := 0;
     NTT[1] := TT[1] * COS_ROT_1 -  TT[3]  * SIN_ROT_1;
     NTT[2] := TT[2];
     NTT[3] := TT[1] * SIN_ROT_1 +  TT[3]  * COS_ROT_1;
     TT[1] := NTT[1] * COS_ROT_2 +  NTT[2] * SIN_ROT_2;
     TT[2] :=-NTT[1] * SIN_ROT_2 +  NTT[2] * COS_ROT_2;
     TT[3] := NTT[3];
     YY[1] := TT[1] * RR[1, 1] + TT[2] * RR[1, 2] + TT[3] * RR[1, 3];
     YY[2] := TT[1] * RR[2, 1] + TT[2] * RR[2, 2] + TT[3] * RR[2, 3];
     YY[3] := TT[1] * RR[3, 1] + TT[2] * RR[3, 2] + TT[3] * RR[3, 3];

     TT[1] := 0;
     TT[2] := 0;
     TT[3] := 1;
     NTT[1] := TT[1] * COS_ROT_1 -  TT[3]  * SIN_ROT_1;
     NTT[2] := TT[2];
     NTT[3] := TT[1] * SIN_ROT_1 +  TT[3]  * COS_ROT_1;
     TT[1] := NTT[1] * COS_ROT_2 +  NTT[2] * SIN_ROT_2;
     TT[2] :=-NTT[1] * SIN_ROT_2 +  NTT[2] * COS_ROT_2;
     TT[3] := NTT[3];
     ZZ[1] := TT[1] * RR[1, 1] + TT[2] * RR[1, 2] + TT[3] * RR[1, 3];
     ZZ[2] := TT[1] * RR[2, 1] + TT[2] * RR[2, 2] + TT[3] * RR[2, 3];
     ZZ[3] := TT[1] * RR[3, 1] + TT[2] * RR[3, 2] + TT[3] * RR[3, 3];

     RR[1, 1] := XX[1];
     RR[1, 2] := YY[1];
     RR[1, 3] := ZZ[1];
     RR[2, 1] := XX[2];
     RR[2, 2] := YY[2];
     RR[2, 3] := ZZ[2];
     RR[3, 1] := XX[3];
     RR[3, 2] := YY[3];
     RR[3, 3] := ZZ[3];

     { Step 3 - Convert the observer-sat vector to the mirror coord system. }
     {		Negate the x coordinate (the reflection). }
     {		Convert it back to the equatorial system. }
     {		Convert to RA and Dec and find the angle between that }
     {		position and the sun's. }
     {		Make sure that the sun is in front of the mirror. }
{$IFDEF FLARE }
	{ The observer sees the sat at RA=180 degrees at a distance of 10 units (it's the vector that matters). }
SAT_X_TOPO_EQ := -10;
SAT_Y_TOPO_EQ := 0;
SAT_Z_TOPO_EQ := 0;
{$ENDIF }
     TT[1] := SAT_X_TOPO_EQ * RR[1, 1] +
              SAT_Y_TOPO_EQ * RR[2, 1] +
	      SAT_Z_TOPO_EQ * RR[3, 1];
     TT[2] := SAT_X_TOPO_EQ * RR[1, 2] +
              SAT_Y_TOPO_EQ * RR[2, 2] +
	      SAT_Z_TOPO_EQ * RR[3, 2];
     TT[3] := SAT_X_TOPO_EQ * RR[1, 3] +
              SAT_Y_TOPO_EQ * RR[2, 3] +
	      SAT_Z_TOPO_EQ * RR[3, 3];

     TT[1] := -TT[1];

     IF TT[1] >=0 THEN
     	BEGIN
     	SUN_REF_X := TT[1] * RR[1, 1] + TT[2] * RR[1, 2] + TT[3] * RR[1, 3];
     	SUN_REF_Y := TT[1] * RR[2, 1] + TT[2] * RR[2, 2] + TT[3] * RR[2, 3];
     	SUN_REF_Z := TT[1] * RR[3, 1] + TT[2] * RR[3, 2] + TT[3] * RR[3, 3];
     	SUN_REF_ALPHA := ARCTAN_D2(SUN_REF_Y, SUN_REF_X);
     	SUN_REF_DELTA := ASIN_D(SUN_REF_Z /
     				SQRT(SQR(SUN_REF_X) +
				     SQR(SUN_REF_Y) +
				     SQR(SUN_REF_Z)));
     	SUN_REF_ANG := ANG_DIS(SUN_ALPHA, SUN_DELTA,
     			       SUN_REF_ALPHA, SUN_REF_DELTA);
        END
     ELSE
     	SUN_REF_ANG := 180;
{$IFDEF FLARE }
	{ Fake the sun's position at RA=0 degrees. }
SUN_REF_ANG := ANG_DIS(0, 0, SUN_REF_ALPHA, SUN_REF_DELTA);
	{ Let's sum up: }
	{ The observer is directly between the sun and the sat. }
	{ The sun is at RA=0 and the sat at RA=180. }
	{ The sat's up points at us and the mirror is pointing straight up. }
	{ This had better reflect our gaze right back at the sun.  SUN_REF_ANG = 0. }
	{ Now try rotations of (45, 0) or (45, 90).  SUN_REF_ANG should be 90. }
{$ENDIF }
     SAT_REFLECTION := SUN_REF_ANG;
END;
=====================
This is how it is called.
     SAT_REF_MIRROR := 1;
     SAT_REF_ANG :=
     	SAT_REFLECTION(EQ_X, EQ_Y, EQ_Z,
     		       V_X, V_Y, V_Z,
        	       SAT_X_TOPO_EQ, SAT_Y_TOPO_EQ, SAT_Z_TOPO_EQ,
        	       ALPHA, DELTA,
		       -40, 0);

     TEMP :=
     	SAT_REFLECTION(EQ_X, EQ_Y, EQ_Z,
     		       V_X, V_Y, V_Z,
        	       SAT_X_TOPO_EQ, SAT_Y_TOPO_EQ, SAT_Z_TOPO_EQ,
        	       ALPHA, DELTA,
		       -40, 120);
     IF TEMP < SAT_REF_ANG THEN
     	BEGIN
        SAT_REF_MIRROR := 2;
     	SAT_REF_ANG := TEMP;
        END;

     TEMP :=
     	SAT_REFLECTION(EQ_X, EQ_Y, EQ_Z,
     	 	       V_X, V_Y, V_Z,
        	       SAT_X_TOPO_EQ, SAT_Y_TOPO_EQ, SAT_Z_TOPO_EQ,
        	       ALPHA, DELTA,
		       -40, -120);
     IF TEMP < SAT_REF_ANG THEN
        BEGIN
        SAT_REF_MIRROR := 3;
     	SAT_REF_ANG := TEMP;
        END;