| angles=satang4(strSys1,latA,longA,longsat,strSys2,latB,longB)
|
function angles=satang4(strSys1,latA,longA,longsat,strSys2,latB,longB)
%function angles=satang4(strSys1,latA,longA,longsat,strSys2,latB,longB)
%calculation of satellite angles: elevation, azimuth, and also relative azimuth and baseline orientation if input data
%belongs to two sites
%Sources:
%Digital satellite communications, 2nd Ed. Author Tri Ha, for satellite azimuth and elevation calculation
%http://williams.best.vwh.net/avform.htm#Rhumb for distance and relative azimuth between two coordinates
%tested using various satellites and earth stations, and using http://www.sfwx.com/softworks/geoBearing.html and
%http://www.geocities.com/senol_gulgonul/linkbudget/azel.html
%differences are in part due to the selection of the earth radius and orbit heigth values, and also to possible mistakes in the
%code...please point them out
%INPUTS
%strSys1: string: 'E' for longitude in degrees east, 'W' for longitude in degrees West. no negative values
%strSys2: string: 'GC' for GREAT CIRCLE , 'RH' for RHUMB (CONSTANT AZIMUTH: recommended for site diversity purposes)
%latXX: latitude of site A or B. + north, - south
%longXX: longitude of site A and B: following strSys1.
%longsat: longitude of subsatellite point, E or W according to srtSys
%OUTPUTS
%angles: a structure containing:
% .elevation elevation angle to satelite,
% .azimuth azimuth to satelite
% .poltilt & poltilt_c feed rotation to compensate geographic position
%NOTE standing behind the antenna,
%rotation is CLOCKWISE for a station
%located east of the satellite meridian,
%and ANTICLOCKWISE for a station west of
%the satellite meridian
%poltilt_c is more reliable than poltilt.
% .Az_AB true course, azimuth in the direction of site B (rhumb line, constant azimuth)
% .distance distance between the two sites
% .baselineorientation for site diversity, BOA acording to ITU-R P618-8 definition.
%EXAMPLE OF USE
%angles=satang4('E',-8.83,13.25,0) | angles=satang('E',-33.93,18.47,0)
%angles = |
% elevation: 71.2774 | elevation: 45.7078
% azimuth: 303.1005 | azimuth: 329.1038
%Meteosat satellite, cities of cape town and luanda.
%
%angles=satang4('E',6.4,360-76.4,304.5)
%angles =
% elevation: 64.3472
% azimuth: 106.2730 Intelsat 8A, city of medellin
%geographical constants
% earthradius=6378.155; %according to TRI Ha book, page 41
earthradius=6371; %http://www.movable-type.co.uk/scripts/latlong.html
satheigth=41264.2;
I=find(longA<0);
if isempty(I)==0
longA(I)=360+longA(I);
end
clear I
[Az El poltilt poltilt_c]=pointangle(strSys1,latA,longA,longsat);
angles.Az_tosatA=Az;
angles.elevationA=El;
angles.poltiltA=poltilt;
angles.poltiltA_c=poltilt_c; %poltilt using complex method
if nargin>5 && strcmp('GC',strSys2)
IB=find(longB<0);
if isempty(IB)==0
latB(IB)=360+longB(IB);
end
clear IB
[Az El poltilt poltilt_c]=pointangle(strSys1,latB,longB,longsat);
%calculate distance, orientation
angles.Az_tosatB=Az;
angles.elevationB=El;
angles.poltiltB=poltilt;
angles.poltiltB_c=poltilt_c;
dlatn=latB-latA;
dlongn=longB-longA;
an=sin(deg2rad(dlatn/2))^2+cos(deg2rad(latA))*cos(deg2rad(latB))*sin(deg2rad(dlongn/2))^2;
cn=2*atan2(sqrt(an),sqrt(1-an));
baselineL=earthradius*cn;
Az_AB=atan2(sin(deg2rad(dlongn)).*cos(deg2rad(dlatn)),cos(deg2rad(latA)).*sin(deg2rad(latB))-sin(deg2rad(latA)).*cos(deg2rad(latB)).*cos(deg2rad(dlongn)));
Az_AB=mod(Az_AB*(180/pi)+360,360);
angles.Az_AB=Az_AB;
angles.distance=baselineL;
if abs(angles.Az_tosatA-Az_AB)>90
baselineO=abs(180-abs(angles.Az_tosatA-Az_AB)); %reference site A
else
baselineO=abs(angles.Az_tosatA-Az_AB);
end
angles.baselineorientation=baselineO;
elseif nargin>5 && strcmp('RH',strSys2)
% dfin=log(tan(deg2rad((latB/2)+(pi/4)))/tan(deg2rad((latA/2)+(pi/4))));
[Az El poltilt poltilt_c]=pointangle(strSys1,latB,longB,longsat);
%calculate distance, orientation
angles.Az_tosatB=Az;
angles.elevationB=El;
angles.poltiltB=poltilt;
angles.poltiltB_c=poltilt_c;
latBr=deg2rad(latB);latAr=deg2rad(latA);
%calculate distance between sites, rhumb, assuming SPHERICAL earth
dlon_W=mod(deg2rad(longB)-deg2rad(longA),2*pi);
dlon_E=mod(deg2rad(longA)-deg2rad(longB),2*pi);
% (*)
dphi=log(tan(latBr./2+pi/4)./tan(latAr./2+pi/4));
% (*)
if abs(latBr-latAr) < sqrt(1E-15)
q=cos(latAr);
else
q=(latBr-latAr)./dphi;
%(*)
end
if (dlon_W < dlon_E)%// Westerly rhumb line is the shortest
Az_AB=rad2deg(mod(atan2(-dlon_W,dphi),2*pi));
% DISTANCIA(*)
baselineL= (sqrt(q.^2.*dlon_W.^2 + (latBr-latAr).^2)).*earthradius;
%(*)
else
%Az_AB=rad2deg(mod(atan2(dlon_E,dphi),2*pi));
Az_AB=rad2deg(atan2(deg2rad(longB-longA),dphi));
if Az_AB<0
Az_AB=360+Az_AB;
end
baselineL= earthradius.*(sqrt(q.^2*dlon_E.^2 + (latBr-latAr).^2)); %distance is arc length:radius times angle subtended
end
if abs(angles.Az_tosatA-Az_AB)>90
baselineO=180-abs(angles.Az_tosatA-Az_AB); %reference site A
else
baselineO=abs(angles.Az_tosatA-Az_AB);
end
angles.Az_AB=Az_AB;
angles.distance=baselineL;
angles.baselineorientation=baselineO;
end
% Suppose point 1 is LAX: (33deg 57min N, 118deg 24min W)
% Suppose point 2 is JFK: (40deg 38min N, 73deg 47min W)
%
return
%subroutine for pointing angle calculation
function [Az El Poltilt Poltilt_c]=pointangle(strSys1,lat,long,longsat)
earthradius=6378.155; %according to TRI Ha book, page 41
satheight=41264.2; %distance from satellite to center of earth in kilometers
satB=longsat-long; %relative longitude of ES wrt to satellite
c=90-lat;
satcosb=cosd(satB).*cosd(lat);
satb=acosd(satcosb); %intelsat's alfa angle
d_tosat=sqrt(satheight^2+earthradius^2-2.*satheight.*earthradius.*satcosb);
coselev_tosat=satheight.*sind(satb)./d_tosat;
elev_tosat=acosd(coselev_tosat);
%calculation of azimuth and elevation angle to a satellite from site A
%azimuth according to tri ha
aprime=atand(tand(abs(satB))./sind(abs(lat)));
%determine position of es relative to sat
for rc=1:size(aprime)%aca toca meter un for....
if (strcmp('E',strSys1)&&lat(rc)>0&&satB(rc)>0)||(strcmp('W',strSys1)&&lat(rc)>0&&satB(rc)<0) %if satB>0 (deg E) or satB<0 (deg W), station is at the west of the satellite, both long values are in E degrees
Az_tosat(rc)=180-aprime(rc);
elseif (strcmp('E',strSys1)&&lat(rc)>0&&satB(rc)<0)||(strcmp('W',strSys1)&&lat(rc)>0&&satB(rc)>0)%NH, east of satellite
Az_tosat(rc)=180+aprime(rc);
elseif (strcmp('E',strSys1)&&lat(rc)<0&&satB(rc)>0)||(strcmp('W',strSys1)&&lat(rc)<0&&satB(rc)<0) %southern hemisphere, west of satellite
Az_tosat(rc)=aprime(rc);
else
Az_tosat(rc)=360-aprime(rc);
end
end
El=elev_tosat;
Az=Az_tosat';
Poltilt=90-atand( sind(abs(satB))./tand(lat) ); %simplified method Maral, Bousquet: Satellite communication systems
Re_factor=earthradius /( earthradius+(satheight-earthradius) ) ;
Poltiltc_num=sind(lat).*( 1-Re_factor.*cosd(abs(satB)).*cosd(lat) );
Poltiltc_den_1=( ( cosd(lat).*sind(abs(satB)) ).^2 + sind(lat).^2 ) ;
Poltiltc_den_2=Re_factor.^2.*cosd(lat).^2 + 1 - 2.*( Re_factor.*cosd(abs(satB)).*cosd(lat) );
Poltiltc = Poltiltc_num ./sqrt( Poltiltc_den_1.*Poltiltc_den_2 )
Poltilt_c = 90 - acosd(Poltiltc);
%complex method: Maral, Bousquet: Satellite communication systems
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
return |
|
