<< Handling time scales Cookbook Calculation of atmospheric drag >>

celestlab >> - Introduction - > Cookbook > B-plane

B-plane

B-plane calculations

Definition of B-plane and associated frame

The B-plane is a plane that is orthogonal to the hyperbolic trajectory plane and the hyperbolic excess velocity vector. This is used to help define the hyperbolic trajectory in 3D space.

Three orthogonal unit vectors are defined:

The vector T is chosen perpendicular to some reference plane. Here we'll assume this reference plane is the Equator of the planet. More precisely, the B-plane frame will be defined by:

In CelestLab the (T, R, S) frame can be computed by: M = CL_rot_defFrameVec(vinf, [0; 0; -1], 3, 2) where M is the transformation matrix from the implicit frame (where vinf coordinates are defined) to the (T, R, S) frame. Note that the frame in which vinf coordinates are defined is here supposed equatorial because the polar axis has coordinates [0; 0; 1].

The B-point is the intersection between the asymptote and the B-plane and can be defined by its coordinates in the (T, R, S) frame.

Application 1: Keplerian elements to B-point coordinates in (T, R, S) frame

// -----------------------------------------------------
// Application 1: 
// Arrival trajectory at Mars: computation of the B-point coordinates 
// -----------------------------------------------------
// -- Inputs 
// Orbital elements of the incoming trajectory at Mars in equatorial frame
kep = [4600.e3; 1.8; 60 * %pi / 180; 315 * %pi / 180; 0; 0];  
// Gravitational constant 
MU = CL_dataGet("body.Mars.mu"); 

// Conversion to position and velocity (at the pericenter) 
[pos, vel] = CL_oe_kep2car(kep, mu = MU);  

// Conversion to position and velocity vectors at infinity (incoming asymptote)
// NB: pinf is the B-point: intersection between the asymptote and the B-plane. 
[pinf, vinf] = CL_ip_flybyVectors("pv", pos, vel, ["pinfa", "vinfa"], mu = MU);

// Calculation of the B-plane frame 
M_bplane = CL_rot_defFrameVec(vinf, [0;0;-1], 3, 2); 

// B-point spherical coordinates in (T, R, S) frame
// - 1st component: angle with "T" axis in B-plane
// - 3rd component: distance from Mars center 
// NB: 2nd component is 0 as the B-point lies in the B-plane. 
pinf_bplane_sph = CL_co_car2sph(M_bplane * pinf); 

// ------------------------
// Note 
// ------------------------
// Let's call thetaB, the angle between T axis and the B-point in the B-plane. 
// We can check that: 
// -> cos(thetaB) = cos(delta) * cos(inc) + sin(delta) * sin(inc) * sin(alpha - raan)
// -> sin(thetaB) = sin(inc) * cos(alpha - raan)
// where alpha and delta are respectively the right ascension and declination of 
// the excess velocity vector. 
// 
// Excess velocity vector in spherical coordinates
vinf_sph = CL_co_car2sph(vinf); 
alpha = vinf_sph(1); // right ascension 
delta = vinf_sph(2); // declination 
// 
// Some additional variables for simplicity
inc = kep(3); // inclination
raan = kep(5); // right ascension of ascending node 
thetaB = pinf_bplane_sph(1); // thetaB as computed previously
//
// Check formulas (printed results should be 0)  
cos_thetaB = cos(delta) * cos(inc) + sin(delta) * sin(inc) * sin(alpha - raan) 
sin_thetaB = sin(inc) * cos(alpha - raan)
disp(cos_thetaB - cos(thetaB))
disp(sin_thetaB - sin(thetaB))

Application 2: Excess velocity vector and B-point coordinates to Keplerian elements

// -----------------------------------------------------
// Application 2: 
// Arrival trajectory at Mars: computation of the Keplerian elements  
// Note: this is the inverse problem to application 1
// -----------------------------------------------------
// -- Inputs 
// Excess velocity vector in equatorial frame
vinf = [2992.7; 297.7; 515.6]; 
// Spherical coordinates of B-point in (T, R, S) frame 
pinf_bplane_sph = [59.52 * %pi / 180; 0; 6885.e3]; 
// Gravitational constant 
MU = CL_dataGet("body.Mars.mu"); 

// Computation of the (T, R, S) frame
M_bplane = CL_rot_defFrameVec(vinf, [0;0;-1], 3, 2); 

// Coordinates of B-point in initial (implicit) frame 
pinf = M_bplane' * CL_co_sph2car(pinf_bplane_sph); 

// Position and velocity at pericenter
[pos, vel] = CL_ip_flybyVectors("pvinfa", pinf, vinf, ["posp", "velp"], mu = MU);

// Conversion to Keplerian elements 
kep = CL_oe_car2kep(pos, vel, mu = MU);


Report an issue
<< Handling time scales Cookbook Calculation of atmospheric drag >>