Matlab Convert Aerospace Coordinates

% Always: convXXX(value, 'from', 'to')
alt_m = convlength(35000, 'ft', 'm');        % 10668.0 m
speed_ms = convvel(250, 'kts', 'm/s');       % 128.61 m/s
angle_rad = convang(45, 'deg', 'rad');       % 0.7854 rad
accel_g = convacc(9.81, 'm/s^2', 'G''s');   % 1.0 G
omega_rpm = convangvel(360, 'deg/s', 'rpm'); % 60 rpm

% Julian Date from components (year, month, day, hour, min, sec)
jd = juliandate(2024, 6, 15, 12, 0, 0);     % 2460477.0

% Modified Julian Date
mjd = mjuliandate(2024, 6, 15, 12, 0, 0);   % 60476.5

% From datetime objects
dt = datetime(2024, 6, 15, 12, 0, 0);
jd = juliandate(dt);

% Decimal year
dy = decyear(2024, 6, 15);                   % 2024.4536

% TDB Julian Date from Terrestrial Time [yr mo day hr min sec]
tt = [2024 6 15 12 0 0];
jdTDB = tdbjuliandate(tt);

% LLA is [latitude_deg, longitude_deg, altitude_m]
lla = [40, -74, 0];
ecef = lla2ecef(lla);           % [1348613.0, -4703172.4, 4077985.6] m

% Round-trip verification
lla_check = ecef2lla(ecef);     % [40.0, -74.0, 0.0]

% ECI transforms require UTC time
lla = [40, -74, 1000];
utc = [2024 6 15 12 0 0];
posECI = lla2eci(lla, utc);

% Back to LLA
lla_check = eci2lla(posECI, utc);

% Convert ECI position/velocity to Keplerian orbital elements
rijk = [-2981784; 5207055; 3161595];   % ECI position (m), 3x1
vijk = [-3384; -4887; 4843];           % ECI velocity (m/s), 3x1
[a, ecc, incl, RAAN, argp, nu] = ijk2keplerian(rijk, vijk);
% a=m, angles in degrees

% Convert Keplerian elements back to ECI position/velocity
[rijk_back, vijk_back] = keplerian2ijk(a, ecc, incl, RAAN, argp, nu);

% For a different central body (e.g., Mars)
[a_mars, ecc_mars, incl_mars, RAAN_mars, argp_mars, nu_mars] = ...
    ijk2keplerian(rijk, vijk, CentralBody="Mars");

% Compute Greenwich Mean Sidereal Time (GMST) and Apparent (GAST)
utcJD = juliandate(2019, 1, 4, 12, 0, 0);  % UTC as Julian date
dUT1 = 0;    % UT1-UTC offset (seconds)
dAT  = 37;   % TAI-UTC offset (leap seconds)
[thGMST, thGAST] = siderealTime(utcJD, dUT1, dAT);
% thGMST in degrees, thGAST in seconds

% Good for short-range simulations (< ~100 km from reference)
llo = [40, -74];    % reference lat/lon (deg)
psio = 0;           % angular direction of flat Earth x-axis (rad, 0=North)
href = 0;           % reference height (m)

lla_point = [40.01, -73.99, 100];
flatPos = lla2flat(lla_point, llo, psio, href);  % [x, y, z] in meters
lla_back = flat2lla(flatPos, llo, psio, href);

% CRITICAL: angle2dcm expects RADIANS, not degrees
yaw = deg2rad(30); pitch = deg2rad(10); roll = deg2rad(5);
dcm = angle2dcm(yaw, pitch, roll, 'ZYX');

% Extract angles back (returns radians)
[y, p, r] = dcm2angle(dcm, 'ZYX');
fprintf('Yaw=%.1f, Pitch=%.1f, Roll=%.1f deg\n', rad2deg(y), rad2deg(p), rad2deg(r));

% Aerospace Toolbox quaternion format: [q0 q1 q2 q3] = [scalar, vector]
q = dcm2quat(dcm);       % 1x4, scalar-first
dcm_back = quat2dcm(q);  % 3x3

% Verify orthogonality
err = max(abs(dcm' * dcm - eye(3)), [], 'all');
assert(err < 1e-14, 'DCM is not orthogonal');

% Create from Euler angles — specify sequence AND frame type
q = quaternion([30 10 5], 'eulerd', 'ZYX', 'frame');

% Rotate a point (active rotation)
pt = [1 0 0];
pt_rotated = rotatepoint(q, pt);   % [0.8529, 0.4924, -0.1736]

% Compose rotations by multiplication
q_total = q2 * q1;  % applies q1 first, then q2

% Always normalize after arithmetic accumulation
q = normalize(q);

% Smooth interpolation between two orientations
q1 = quaternion([0 0 0], 'eulerd', 'ZYX', 'frame');
q2 = quaternion([90 0 0], 'eulerd', 'ZYX', 'frame');
t = linspace(0, 1, 5)';
q_interp = slerp(q1, q2, t);

% Verify: yaw progresses linearly for pure yaw rotation
e = eulerd(q_interp, 'ZYX', 'frame');
% e(:,1) = [0, 22.5, 45, 67.5, 90]

% Compose two rotations: q2 applied after q1
q1 = angle2quat(deg2rad(30), 0, 0, 'ZYX');
q2 = angle2quat(0, deg2rad(10), 0, 'ZYX');
q_total = quatmultiply(q2, q1);  % q2 * q1 (apply q1 first)

% Inverse rotation
q_inv = quatinv(q1);
q_identity = quatmultiply(q1, q_inv);  % [1 0 0 0]

% Rotate a vector
v = [1 0 0];
v_rot = quatrotate(q1, v);  % rotate v by q1

% Normalize after accumulation
q_accumulated = quatnormalize(q_total);

% Interpolation (slerp, lerp, nlerp)
q_mid = quatinterp(q1, q2, 0.5, 'slerp');

% Relative rotation: what rotation takes q1 to q2?
q_rel = quatdivide(q2, q1);  % q_rel such that q2 = q_rel * q1

% Exponential/logarithm (useful for angular velocity integration)
q_log = quatlog(q1);         % maps to tangent space
q_back = quatexp(q_log);     % back to quaternion

% Fractional rotation (half the rotation of q1)
q_half = quatpower(q1, 0.5);

% Compact 3-element representation (singular at 180 deg)
rod = dcm2rod(dcm);           % 1x3 vector
dcm_back = rod2dcm(rod);

% Convert between all representations
rod = quat2rod(q_array);      % q is M×4 (scalar-first)
q_back = rod2quat(rod);

% Body-to-wind DCM from angle of attack and sideslip
alpha = deg2rad(5);
beta = deg2rad(2);
dcm_bw = dcmbody2wind(alpha, beta);

% Extract alpha/beta from a DCM
[alpha_out, beta_out] = dcm2alphabeta(dcm_bw);

% Body-to-stability (alpha only, no sideslip)
dcm_bs = dcmbody2stability(alpha);

% Get the DCM to rotate vectors from ECEF to local NED
lat = deg2rad(40);
lon = deg2rad(-74);
dcm_ecef2ned = dcmecef2ned(lat, lon);

% Transform an ECEF velocity to NED
v_ecef = [10; 20; 30];
v_ned = dcm_ecef2ned * v_ecef;

% Compute look angles from a ground station to a satellite in ECI
posECI = [-2981784, 5207055, 3161595];  % satellite ECI position (m)
utc = [2019 1 4 12 0 0];               % observation time
lla0 = [28.5, -80.5, 0];               % ground station [lat, lon, alt] (deg, deg, m)

aer = eci2aer(posECI, utc, lla0);
% aer = [azimuth_deg, elevation_deg, slant_range_m]

% Detect gimbal lock risk
[~, pitch, ~] = dcm2angle(dcm, 'ZYX');
if abs(pitch) > deg2rad(85)
    warning('Near gimbal lock — use quaternion representation');
end