The function has worked great till I added the xm and ym output arguments (and the "if" structure). Why does this occur now?
function [glongitude, date, longitude, latitude, gpt, xm, ym] = GroundTrack(VettXs, ... VettYs, VettZs, omegarot, dt, Nt, RT)[longitude, latitude, ~] = cart2sph(VettXs, VettYs, VettZs);longitude = longitude(:); % expressed in radians
latitude = latitude(:);tom = evalin('base', 'tom'); % identifies the type of map
alphadate = evalin('base', 'alphadate');d = evalin('base', 'd'); ds = string(d);date = strings(Nt); gpt = nan(Nt, 3); date(1) = ds; date(2) = date(1); % Till the "MOLLWEIDE PROJECTION" section a simple conversion from cartesian
% coordinates to polar coordinates is implemented. This conversion is
% then used to obtain longitude (glongitude) and latitude of the projection on Earth of
% the satellite's position. From longitude and latitude is then possible to
% get all the coordinates of the ground track on other maps obtained with
% other cartographical projections.
modulus = sqrt(VettXs(2)^2 + VettYs(2)^2 + VettZs(2)^2);gpt(1, :) = 0;gpt(2, 1) = VettXs(2)*RT/modulus; gpt(2, 2) = VettYs(2)*RT/modulus; gpt(2, 3) = VettZs(2)*RT/modulus; it = seconds(dt);glongitude = nan(Nt, 1); glongitude(1) = 0;glongitude(2) = longitude(2) - alphadate;n = glongitude(2)/(2*pi); % number of complete turns
if (n - floor(n)) <= 0.5 glongitude(2) = 2*pi*(n - floor(n)); glongitude(2) = glongitude(2)*180/pi;else glongitude(2) = -2*pi*(1 - (n - floor(n))); glongitude(2) = glongitude(2)*180/pi;endfor i = 3 : Nt beta = omegarot*dt*(i-2); glongitude(i) = longitude(i) - beta - alphadate; n = glongitude(i)/(2*pi); d = d + it; ds = string(d); date(i) = ds;if (n - floor(n)) <= 0.5 glongitude(i) = 2*pi*(n - floor(n)); glongitude(i) = glongitude(i)*180/pi;else glongitude(i) = -2*pi*(1 - (n - floor(n))); glongitude(i) = glongitude(i)*180/pi;end modulus = sqrt(VettXs(i)^2 + VettYs(i)^2 + VettZs(i)^2); gpt(i, 1) = VettXs(i)*RT/modulus; gpt(i, 2) = VettYs(i)*RT/modulus; gpt(i, 3) = VettZs(i)*RT/modulus;end% MOLLWEIDE PROJECTION
if (tom == 2) teta = nan(length(latitude), 1); xm = nan(length(longitude), 1); ym = nan(length(latitude), 1); % Iterative cycle
for i = 1 : Nt % First step
tetaf = latitude(i); dtetaf = -(tetaf + sin(tetaf) - pi*sin(latitude(2)))/(1 + cos(tetaf)); % Other steps
while abs(dtetaf) >= 1e-6 dtetaf = -(tetaf+sin(tetaf)-pi*sin(latitude(i)))/(1+cos(tetaf)); tetaf = tetaf + dtetaf; end teta(i) = tetaf/2; xm(i) = (RT*2*sqrt(2)/pi)*(glongitude(i)-glongitude(2))*cos(teta(i)); ym(i) = RT*sqrt(2)*sin(teta(i)); endendend
Best Answer