Flugrouten

Auf der duchsichtig dargestellten Erdkugel werden die Flugrouten gezeigt, welche von Frankfurt aus entlang von Grosskreisen zu verschiedenen Zielen führen. Man beachte die Tatsache, dass für Verbindungen auf der Nordhalbkugel die Grosskreise sehr weit nach Norden Richtung Pol ausholen.

function voidbk = airportfilms(varargin)
  % uses:  grosskreisfct, breitkreisfct, grosskreis, flugfct
  %  plot3ers
  figure(1); clf
     voidbk = 0;
     rad = 8; flypaus = 0.04; townpaus = 1.5; townnxpaus = 0.4;
% showflug.m Beipiel-Plots f\"ur Flugrouten
hdtx = text(-9.3,8.5,'From Airport Frankfurt to', ...
    'FontUnits','points','FontSize',20 , ...
   'Fontweight','demi');
axis([-10 10 -10 10]); axis square;
hold on
axis off
wipl = 8*pi/180; tepl = -18*pi/180;
% meridiannetz.m  Zeichnen eines Meridian-Netzes
for klg = 1:4
  long = (klg-1)*30;
  pl = rad*grosskreisfct(long, 360);
  if long == 0
% Muster  hdkr1 = plot3ers(x,y,zc1, 0,tet,0,0,'g',2) ; hold on; 
    hdzm = plot3ers(pl(1,:),pl(2,:),pl(3,:), wipl,tepl, 0,-1,'k',2) ;
  else 
    hdgke(klg-1) =   plot3ers(pl(1,:),pl(2,:),pl(3,:), wipl,tepl, 0,-1,'b',1) ;
    hdgkw(klg-1) =  plot3ers(-pl(1,:),pl(2,:),pl(3,:), wipl,tepl, 0,-1,'b',1) ;
  end
end
for klat = -4:4
  lat = klat*15;
 pb = rad*breitkreisfct(lat);
  if lat == 0
     hdbr(klat+5) = plot3ers(pb(1,:),pb(2,:),pb(3,:), wipl,tepl, 0,-1,'k',2) ;
  else 
    hdbr(klat+5) =  plot3ers(pb(1,:),pb(2,:),pb(3,:), wipl,tepl, 0,-1,'b',1) ;
  end
end

%Start = Frankfurt
longs = 8.6833; lats = 50.1167; 
cities{1} = 'Neapel';    cities{2} = 'Athen';
cities{3} = 'Bejing';    cities{4} = 'Madrid';
cities{5} = 'Chicago';   cities{6} = 'Mumbay';
cities{7} = 'Kapstadt';  cities{8} = 'Mexico';
cities{9} = 'Moskau';    cities{10}= 'Miami';
cities{11} = 'Rio';      cities{12} = 'Sydney';
cities{13} = 'Tokyo';    cities{14} = 'Frisco';

%      Neapel,   Athen,  Bejing,  Madrid,  Chicago  Bombay,   Kapstadt,
%      London
%       Mexico, Moskau, NYC , Rio, Sydney , Tokyo,  Frisco,  
lat = [40.833  37.9667   39.9167   40.433   41.7833  19.0    -35.9167 ...
         19.4332  55.75  26.76  -22.95  -34.0  35.6667   37.6167    ];
long = [14.25   23.7167  116.25    -3.7    -95.35  72.8     18.3667  ...
        -99.1167  37.6  -80.66   -43.2  151.0  139.75  -122.3833   ];
%
 for k = 1:length(long)
   if k > 1
      delete(hdcity);
   end
 % Stadt anzeigen
   hdcity = text(-9.6, -9.6,char(cities{k}), ...
    'FontUnits','points','FontSize',20 , ...
   'Fontweight','demi');
  [dist,azi,maxbreit,norvec] = grosskreis(longs,lats,long(k),lat(k));
  pl = rad*flugfct(longs,lats,azi,dist/6336.2*180/pi);

  ll = length(pl);
  for dkl = 1:20
    dpart = 5*dkl;
    if dkl > 1 
      delete(hdflac);
    end
    hdflac =   plot3ers(pl(1,1:dpart), pl(2,1:dpart), ...
        pl(3,1:dpart), wipl,tepl, 0,-1,'r',2) ;
    pause(flypaus)
  end
  hdflg(k) = hdflac;
  pause(townpaus);
end
pause(townnxpaus)
%
% grosskreisfct.m Funktion zum Erzeugen der Koordinatentripel
%    eines Laengenmeridian-Sektors zur L\"ange long 
%    und der Winkeldistanz dist ab dem Nordpol
function coordmat = grosskreisfct(long,dist)
t = linspace(0,dist*pi/180);
coordmat = [sin(t)*cos(long*pi/180);  sin(t)*sin(long*pi/180); cos(t)];
end
%
% breitkreisfct.m Funktion zum Erzeugen der Koordinatentripel
%    eines Breitnmeridians der Breite lati 
function coordmat = bpeitkreisfct(lati)
t = linspace(0,2*pi);
coordmat=[cos(t)*cos(lati*pi/180); sin(t)*cos(lati*pi/180);...
 sin(lati*pi/180)*ones(1,length(t)) ];
end
%
function [dist,azi,maxbreit,norvec] = grosskreis(long1,lat1,long2,lat2)
%[dist,azi,maxbreit,norvec] = grosskreis(long1,lati1,long2,lati2)
% Berechnung von Distanz, Startazimut, maximaler Breite und Normalenvektor
% zur Grosskreis-Flugroute zwischen zwei Orten auf der Welt, deren
% Laengenposition (Ost=+, West= -) und Breitenposition (N=+) bekannt ist.
dgr = pi/180;
v1=[cos(lat1*dgr)*cos(long1*dgr) cos(lat1*dgr)*sin(long1*dgr) ...
     sin(lat1*dgr)]';
v1n=[-sin(lat1*dgr)*cos(long1*dgr) -sin(lat1*dgr)*sin(long1*dgr) ...
     cos(lat1*dgr)]';
v2=[cos(lat2*dgr)*cos(long2*dgr) cos(lat2*dgr)*sin(long2*dgr) ...
      sin(lat2*dgr)]';
dist=abs(acos(v1'*v2))*6366.2;  % Winkel in rad = acos(Skalarprodukt)
if dist > 0   % Ausschliessen identische Punkte
  norvec = cross(v1,v2);
  if norm(norvec) > 0;   % Ausschliessen Antipoden
     norvec = norvec/norm(norvec);
     maxbreit =  acos(norvec(3)/norm(norvec))/dgr;
     takoff = cross(norvec,v1);
     norcomp =  (takoff'*v1n/norm(v1n))*v1n;
     westcomp = takoff-norcomp;
     azi = acos(takoff'*v1n/norm(takoff)/norm(v1n))/dgr;
     if long2 > long1 ;
        azi = -azi;
     else
        maxbreit = 180-maxbreit;
     end;
   else; disp('Antipoden!'); end;
else;  disp('identische Punkte!');  end
end
%
% plk=flugfct(long,lat,azi,dist) Flugroute ab Startort long lat
%   mit Abflug-Azimuth azi und Flugdistanz dist
function plk=flugfct(long,lat,azi,dist)
pm = grosskreisfct(180+azi, dist);
wy = (-90+lat)*pi/180;
Rym = [ cos(wy) 0 -sin(wy); 0 1 0; sin(wy) 0 cos(wy)];
wz = long*pi/180;
Rzm = [ cos(wz)  -sin(wz) 0 ;  sin(wz)  cos(wz) 0 ; 0 0 1];
plk = Rzm*Rym*pm;
end
%
function plhdb = plot3ers(x,y,z,w,t,xof,yof,col,thik)
% function plhdb = plot3ers(x,y,z,w,t,xof,yof,col,thik)
% Ersatz fuer plot3
%  mit Zusatz-Parametern   w , t horizontaler und vertikaler Drehwinkel
%  xof, yof  x- und y- Offset (Zusatz) nach Projektion in 2D
%  col  Farbe und thik Liniendicke
%  Rueckgabe plhdb  plot-handle des 2D plot-Aufrufs
zrm = [cos(w) -sin(w) 0;  sin(w) cos(w) 0 ; 0 0 1];
yrm = [  cos(t) 0 -sin(t) ; 0 1 0;  sin(t) 0  cos(t) ];
trko = yrm*zrm*[x;y;z];
plhdb = plot(trko(2,:)+xof, trko(3,:)+yof,col,'LineWidth',thik); 
end
end