✅作者简介：热爱科研的Matlab仿真开发者，修心和技术同步精进，matlab项目合作可私信。
🍎个人主页：Matlab科研工作室
🍊个人信条：格物致知。
更多Matlab仿真内容点击👇
⛄ 内容介绍
基于甚长基线干涉测量技术(Very Long Base Interferometer,VLBI)的时延,时延率以及USB(United S-band)/UXB(United X-band)的测距特点,采用联合统计定位及月面高程约束策略,实现了高精度的探测器的实时单点定位和准实时联合统计定位,且实时单点定位不受力学约束能够快速准确地给出三维位置信息
⛄ 部分代码
%
% function vlbiCoverageMap([{stations}, 'sourcename'])
%
% Plot 2D maps of baseline-elevation against (RA x Dec) for every baseline
% formed by {stations}, and adds the track of source 'sourcename'.
% Allowed station names are those returned by [~,names]=get_vlbi_stations_xyz()
% A NED database lookup is performed to get the source RA,Dec coodinates.
% The "baseline-elevation" is the minimum of the source elevations at
% either antenna on a given baseline.
%
% Requires: Aerospace Toolbox and Mapping Toolbox
%
% Example plots:
% vlbiCoverageMap({'KT','KU','KY'}, 'Sgr A*')
% The plots are exported to files that are
% vlbi-cover-<plotNr>-<station1>-<station2>.pdf
%
function [figfiles] = vlbiSkyCoverage(varargin)
min_elev = 10.0; % minimum elevation (degrees)
freq = 86e9; % observing frequency (Hz) used to estimate baseline length
%% Helpers
fnt = struct();
fnt.FontSize = 10;
fnt.FontWeight = 'bold';
fnt.FontName = 'Times New Roman';
figfiles = {};
%% Args
if ~(length(varargin) == 0 || (length(varargin) == 2))
fprintf(1, 'vlbiCoverageMap() : specify either 2 or zero arguments!\n');
return
end
% Defaults
if (length(varargin) == 0),
stations = {'KU', 'KT', 'Eb', 'Pv', 'CARMA', 'SMA', 'APEX', 'JCMT'};
target = 'NRAO530';
else
stations = varargin{1};
target = varargin{2};
end
%% Geometry
[xyz,~] = get_vlbi_stations_xyz(stations);
lla = ecef2lla(xyz, 'WGS84');
fprintf(1, 'Looking up %s on NED...\n',target);
[~,src] = m_getNED(target);
%% Make grid of sky RA & Dec positions
Nst = size(lla,1);
if (Nst <= 1),
fprintf(1, 'Error: need 2 or more stations, got %d (%s)\n', Nst, char(stations));
return;
end
Npt = 100;
decs = deg2rad(linspace(-90,+90,Npt)); %decs(end+1) = decs(1);
ras = deg2rad(linspace(0,360,Npt)); %ras(end+1) = ras(1);
lats = deg2rad(lla(:,1));
lons = deg2rad(lla(:,2));
[decsMN,rasMN] = meshgrid(decs,ras);
%% Evaluate Az,El across sky RA&Dec grid at site of every antenna/baseline
baselines = nchoosek(1:Nst, 2);
Nbl = size(baselines,1);
for bb=1:Nbl,
s1 = baselines(bb,1);
s2 = baselines(bb,2);
% station 1: alt, az
HA = rasMN + lons(s1);
lat = lats(s1);
[alt1,az1] = m_pos2altaz(lat,HA,decsMN);
% station 2: alt, az
HA = rasMN + lons(s2);
lat = lats(s2);
[alt2,az2] = m_pos2altaz(lat,HA,decsMN);
% baseline: alt clip
if 0,
altbl = (alt1 + alt2) / 2;
else
altbl = min(alt1,alt2);
end
% altbl(altbl < min_elev) = NaN;
% sky map of baseline
figure(1), close(gcf);
figure(1), clf;
axesm('MapProjection','hammer', 'AngleUnits','degrees', 'grid','on', ...
'frame','on', 'meridianlabel','on', 'labelformat','signed', ...
'parallellabel','on', 'mlabelparallel',0);
geoshow(rad2deg(decsMN), rad2deg(-rasMN), altbl, 'DisplayType', 'contour');
hgeo = gca();
box off; hold on;
% overlay the station names
x = rad2deg(lats);
y = rad2deg(lons);
plotm(x(s1), y(s1), 'bo', 'MarkerFaceColor', 'r');
plotm(x(s2), y(s2), 'bo', 'MarkerFaceColor', 'r');
h1 = textm(x(s1), y(s1)+180/50, stations{s1});
h2 = textm(x(s2), y(s2)+180/50, stations{s2});
set(h1,fnt);
set(h2,fnt);
% overlay source track
x = ones(size(ras)) * src.dec;
y = rad2deg(ras);
x(end+1) = x(1);
y(end+1) = y(1);
plotm(x,y, 'k-.');
plotm(src.dec, src.RA, 'bo', 'MarkerFaceColor', 'g');
h = textm(src.dec+5*sign(src.dec), src.RA, target);
set(h,fnt);
% add color scale of elevation
axh = colorbar('Peer',gca());
set(get(axh,'YLabel'),'String','Elevation on baseline (deg)','FontSize',12);
% add title, add axis labels
dlonH = (12/pi)*abs(lons(s1)-lons(s2));
dlatD = (180/pi)*abs(lats(s1)-lats(s2));
lambda = 299792458.0/freq;
dlambda = norm(xyz(s1,:)-xyz(s2,:), 2) / lambda;
[maxEl,maxElIdx] = max(altbl(:));
[M1,M2] = ind2sub(size(altbl),maxElIdx);
s={};
s{end+1} = ['\makebox[4in][c]{', sprintf('%s -- %s', stations{s1}, stations{s2}), '}'];
s{end+1} = ['\makebox[4in][c]{', sprintf('dHA=%.2f\\,h dDec=%.0f\\,deg B=%.0f\\,M$\\lambda$ @ %.0f\\,GHz', dlonH,dlatD,dlambda*1e-6,freq*1e-9), '}'];
s{end+1} = ['\makebox[4in][c]{', sprintf('Maximum El=%.0f\\,deg for Dec=%.0f\\,deg', maxEl,rad2deg(decsMN(M1,M2))), '}'];
title(s, 'Interpreter','latex','FontSize',12,'HorizontalAlignment','center');
xlabel('Celestial RA (deg)');
ylabel('Celestial Dec (deg)');
% fix the declination color scale
set(hgeo,'clim',[-90 90]);
% export to graphics file (PNG: low quality, SVG: not supported in Matlab R2014)
figfiles{end+1} = sprintf('vlbi-cover-%02d-%s-%s', numel(figfiles)+1, stations{s1}, stations{s2});
saveas(gcf(),[figfiles{end} '.pdf'],'pdf');
W = 29.7; H = 21.0; Z = 3.0;
set(gcf,'PaperPositionMode','auto', 'PaperOrientation','portrait', 'PaperUnits','centimeters', 'Units','centimeters');
set(gcf,'PaperSize',[W H], 'PaperPosition',[-Z, -Z, W+Z, H+Z], 'InvertHardcopy','off', 'Resize','on');
print(gcf(),[figfiles{end} '.png'],'-dpng','-r300');
if (bb ~= Nbl), input('next'); end
end
end
function [alt1,az1] = m_pos2altaz(lat,HA,dec)
% http://www.stargazing.net/kepler/altaz.html
% note: input args are in radians, and can be matrices
alt = asin(sin(dec).*sin(lat) + cos(dec).*cos(lat).*cos(HA));
az = (sin(dec) - sin(alt).*sin(lat)) ./ (cos(alt).*cos(lat));
% remove roundoff errors
[i,j] = find(abs(az)>1.0);
az(i,j) = sign(az(i,j));
% change hemisphere
az = acos(az);
az(sin(HA)>0) = deg2rad(360) - az(sin(HA)>0);
alt1 = rad2deg(alt);
az1 = rad2deg(az);
end
function [entry,sentry] = m_getNED(srcname)
% Query that returns ASCII:
nedURL = 'http://ned.ipac.caltech.edu/cgi-bin/objsearch';
nedopts = { 'objname', srcname', ...
'extend','no', ...
'out_csys','Equatorial', ...
'out_equinox','J2000.0', ...
'obj_sort','RA or Longitude', ...
'of','ascii_bar', ...
'img_stamp','NO' };
[str,status] = urlread(nedURL, 'get',nedopts);
% fprintf(1, 'NED returned: %d - %s\n', status, str);
str = strtrim(str);
% Note:
% When a source is found in NED, status==1 and NED's reply looks like:
% No.|Object Name|RA(deg)|DEC(deg)|Type|Velocity|Redshift|Redshift Flag|Magnitude and Filter|Distance (arcmin)|References|Notes|Photometry Points|Positions|Redshift Points|Diameter Points|Associations
% 1|MESSIER 081|148.88822 | 69.06529 |G| -34|-0.000113 | | 7.89||1879|32|186|9|18|8|1
entry = {'', 0, 0, 0, 'xxxx-xx', {}, {}};
sentry = struct();
sentry.name = '';
sentry.RA = 0.0;
sentry.dec = 0.0;
if (status==1),
% TODO: the NED reply *can* contain several lines. For example a
% look-up of NGC 4649 returns 536 source names. How to treat this...
obj = regexp(str,'\n','split');
obj = obj{end};
data = regexp(char(obj),'\|','split');
if (numel(data)<7),
fprintf(1, 'surfDb_getNED: unexpected reply from NED:\n');
data,
return;
end
name = data{2};
name = upper(name);
name = strrep(name,'GROUP','');
vsys = floor(str2double(data{6}));
sra = str2double(data{3});
sdec = str2double(data{4});
sz = str2double(data{7});
entry = {name, vsys, sz, 0.0, 'xxxx-xx', {}, {}, sra, sdec};
sentry.name = name;
sentry.RA = sra;
sentry.dec = sdec;
else
entry = {'', 0, 0.0, 0.0, 'xxxx-xx', {}, {}, 0.0, 0.0};
end
end
⛄ 运行结果
⛄ 参考文献
[1]王霄, 雷辉, 杨旭海,等. 基于VLBI观测的基线改正方法[J]. 测绘科学, 2019, 44(10):7.
