Sometimes calculation of the reflections in astronomy is more simple than in the crystallography. One of the example are the flares from Iridium satellites. This is because all calculations are done directly in Euclidean space. The code presented below helps to calculate all kind of reflections from satellites if the main satellite axis and the panel position in respect to the main satellite axis is known.
This code is written in C# and works perfectly. Additional functions are added below.
The code can be freely used.
The MathEx.Vector4 can be defined as:
public struct Vector4
{
public double x;
public double y;
public double z;
public double w;
};
/**
Module: Satellite flares
Description: Calculate satellite flares for given mirror angles. This code can calculate reflections not only for Iridium
Reference:
Created by: Daniel
Modyfied by:
Data: 2008
Additional : 2013 small acceleration
Input: sun_vector -ECI vector of sun
pos_vector, vel_vector, - satellite vector in space in ECI
obs_vector - observer vector in ECI
alpha - antena angle in xz satellite plane [rad]
beta - antena angle in xy satellite plane [rad]
For Iridium it was: alpha=-40.0deg /angle to the satellite z axis (pointing to the Earth)
3 panels with angles around z axis beta: 0, 120, 240deg
Output: sun_sat_obs angle in [rad] (this is the most important number for 0 it is max of intensity. It is difference between sun-sat-user and real reflection
***************************************************************************************/
double Iridium_Flares(MathEx.Vector4 sun_vector, /*sun vector in ECI*/
MathEx.Vector4 pos_vector, /*satellite vector in ECI*/
MathEx.Vector4 vel_vector, /*satellite velocity vector in ECI*/
MathEx.Vector4 obs_vector, /*observer vector in ECI*/
double alpha, /*antena angle in xz plane*/
double beta, /*antena angle in xy plane*/
out double sun_sat_obs /*out of angle between sun sat and observer*/
)
{
MathEx.Vector4 n1, n2, n3;
MathEx.Vector4 reflected, reflected_back;
double[,] PP = new double[3, 3];
double angle;
//normalize vector v satelite x
MaEx.Scalar_Multiply(1.0, vel_vector, out MathEx.Vector4 xx); MaEx.Normalize(ref xx);
//normalize vector h satelite zz
MaEx.Scalar_Multiply(1.0, pos_vector, out MathEx.Vector4 zz); MaEx.Normalize(ref zz);
//cross xx and yy and normalize
MaEx.Cross(xx, zz, out MathEx.Vector4 yy);
double temp1, temp2;
double nx, ny, nz;
double sin_alpha = Math.Sin(alpha);
double cos_alpha = Math.Cos(alpha);
double sin_beta = Math.Sin(beta);
double cos_beta = Math.Cos(beta);
//x of the mirror in sat unit frame
nx = 1.0; ny = 0.0; nz = 0.0;
//rotate alpha around y (around h vector)
n1.x = nx * cos_alpha - nz * sin_alpha;
n1.y = ny;
n1.z = nx * sin_alpha + nz * cos_alpha;
//rotate around beta (ortogonal axe)
temp1 = n1.x; temp2 = n1.y;
n1.x = temp1 * cos_beta + temp2 * sin_beta;
n1.y = temp1 * -sin_beta + temp2 * cos_beta;
PP[0, 0] = n1.x * xx.x + n1.y * yy.x + n1.z * zz.x;
PP[1, 0] = n1.x * xx.y + n1.y * yy.y + n1.z * zz.y;
PP[2, 0] = n1.x * xx.z + n1.y * yy.z + n1.z * zz.z;
//y of the mirror in sat unit frame
nx = 0.0; ny = 1.0; nz = 0.0;
//rotate alpha around y (around h vector)
n2.x = nx * cos_alpha - nz * sin_alpha;
n2.y = ny;
n2.z = nx * sin_alpha + nz * cos_alpha;
//rotate around beta (ortogonal axe)
temp1 = n2.x; temp2 = n2.y;
n2.x = temp1 * cos_beta + temp2 * sin_beta;
n2.y = temp1 * -sin_beta + temp2 * cos_beta;
PP[0, 1] = n2.x * xx.x + n2.y * yy.x + n2.z * zz.x;
PP[1, 1] = n2.x * xx.y + n2.y * yy.y + n2.z * zz.y;
PP[2, 1] = n2.x * xx.z + n2.y * yy.z + n2.z * zz.z;
//z of the mirror in sat unit frame
nx = 0.0; ny = 0.0; nz = 1.0;
//rotate alpha around y (around h vector)
n3.x = nx * cos_alpha - nz * sin_alpha;
n3.y = ny;
n3.z = nx * sin_alpha + nz * cos_alpha;
//rotate around beta (ortogonal z axe)
temp1 = n3.x; temp2 = n3.y;
n3.x = temp1 * cos_beta + temp2 * sin_beta;
n3.y = temp1 * -sin_beta + temp2 * cos_beta;
PP[0, 2] = n3.x * xx.x + n3.y * yy.x + n3.z * zz.x;
PP[1, 2] = n3.x * xx.y + n3.y * yy.y + n3.z * zz.y;
PP[2, 2] = n3.x * xx.z + n3.y * yy.z + n3.z * zz.z;
//calculate difference between sat vposition vector and observer vector
MaEx.Vec_Sub(pos_vector, obs_vector, out MathEx.Vector4 sat_obs);
MaEx.Normalize(ref sat_obs);
// Convert the observer-sat vector to the mirror coordination system.
reflected.x = sat_obs.x * PP[0, 0] + sat_obs.y * PP[1, 0] + sat_obs.z * PP[2, 0];
reflected.y = sat_obs.x * PP[0, 1] + sat_obs.y * PP[1, 1] + sat_obs.z * PP[2, 1];
reflected.z = sat_obs.x * PP[0, 2] + sat_obs.y * PP[1, 2] + sat_obs.z * PP[2, 2];
reflected.x *= -1.0; //mirror observer location
if (reflected.x >= 0)
{ //reflection in x must be in positive direction
reflected_back.x = reflected.x * PP[0, 0] + reflected.y * PP[0, 1] + reflected.z * PP[0, 2];
reflected_back.y = reflected.x * PP[1, 0] + reflected.y * PP[1, 1] + reflected.z * PP[1, 2];
reflected_back.z = reflected.x * PP[2, 0] + reflected.y * PP[2, 1] + reflected.z * PP[2, 2];
reflected_back.w = 0.0;
MaEx.Vec_Sub(sun_vector, pos_vector, out MathEx.Vector4 sun_sat);//magnitude is calculated inside Vec_Sub
MaEx.Normalize(ref sun_sat);
angle = MaEx.Angle(sun_sat, reflected_back);
sun_sat_obs = pi - MaEx.Angle(sat_obs, sun_sat);
return angle; //angle between sun vector and reflected observer position
}
sun_sat_obs = 0.0;
return 100.0;//nothing calculated
} public void Scalar_Multiply(double k, Vector4 v1, out Vector4 v2)
{
/* Multiplies the vector v1 by the scalar k to produce the vector v2 */
v2.x = k * v1.x;
v2.y = k * v1.y;
v2.z = k * v1.z;
v2.w = Math.Abs(k) * v1.w;
}
public void Vec_Sub(Vector4 v1, Vector4 v2, out Vector4 v3)
{
/* Subtracts vector v2 from v1 to produce v3 */
v3.x = v1.x - v2.x;
v3.y = v1.y - v2.y;
v3.z = v1.z - v2.z;
v3.w = 0.0;
Magnitude(ref v3);
}
public void Normalize(ref Vector4 v)
{
/* Normalize a vector */
v.x /= v.w;
v.y /= v.w;
v.z /= v.w;
}
public double Angle(Vector4 v1, Vector4 v2)
{
/* Calculates the angle between vectors v1 and v2 */
Magnitude(ref v1);
Magnitude(ref v2);
return (Math.Acos(Dot(v1, v2) / (v1.w * v2.w)));
}