Distance and Bearing: Difference between revisions

=={{header|PARI/GP}}==

<syntaxhighlight lang="parigp">/* copy CSV aiports.dat from URL */

csv=externstr("curl 'https://raw.githubusercontent.com/jpatokal/openflights/master/data/airports.dat'");

/* We want an easy convert from CSV to matrix. */

/* Get rid of all those backslash in airports.dat */

{

for(i = 1 , #csv, tmp = Vecsmall(csv[i]);

for(j = 1 , #tmp, if(tmp[j]==92, tmp[j]=78));

csv[i] = strchr(tmp)

);

};

{

convert(csv,mat)=

write1(mat,"[");

for(i=1, #csv-1,

write1(mat,csv[i]); write1(mat,";")

);

write1(mat,csv[#csv]); write(mat,"]")

};

convert(csv,airports);

M = read(airports);

{

distance(lat1, lon1, lat2, lon2, unit = "mi") =

my(dist = 0, rlat1, rlat2, rtheta);

if(lat1!=lat2 || lon1!=lon2,

rlat1 = lat1*Pi/180; rlat2 = lat2*Pi/180;

rtheta = (lon1-lon2)*Pi/180;

dist = sin(rlat1)*sin(rlat2) + cos(rlat1)*cos(rlat2)*cos(rtheta);

dist = acos(min(dist,1))*180/Pi*60*1.1515576;

if(unit=="Km", dist*= 1.609344, unit=="NM", dist*= 0.868976)

);dist

};

{

bearing(lat1, lon1, lat2, lon2)=

my(bear, rlat1, rlat2, rdlon, x, y);

if (lat1!=lat2 || lon1!=lon2,

rlat1 = lat1*Pi/180; rlat2 = lat2*Pi/180;

rdlon = (lon2-lon1)*Pi/180;

y = sin(rdlon)*cos(rlat2);

x = cos(rlat1)*sin(rlat2) - sin(rlat1)*cos(rlat2)*cos(rdlon);

/* Pari/GP is developed for 30 years without built-in atan2() function, very unfortunate. */

if( x>0, bear=atan(y/x), x<0 && y>0, bear=atan(y/x)+Pi,

x<0 && y==0, bear=Pi, x<0 && y<0, bear=atan(y/x)-Pi,

x==0 && y>0, bear=Pi/2, x==0 && y<0, bear=-Pi/2);

bear = (bear*180/Pi+360)%360;

);bear

};

{

request(lat, lon, r = 20 , unit = "NM")=

my(D = vector(#M[,1]), R = matrix(r+1, 5));

R[1,1] = "NAME"; R[1,2] = "COUNTRY"; R[1,3] = "ICAO"; R[1,4] = "DISTANCE"; R[1,5] = "BEARING";

printf("%-36s%-18s\t%-5s\t%s\t%s\n", R[1,1], R[1,2], R[1,3], R[1,4], R[1,5]);

for(i = 1, #M[,1],

D[i] = [distance(lat, lon, M[i,7], M[i,8], unit) , i];

);

D=vecsort(D);

for(i = 1, r,

R[i+1,1] = M[D[i][2],2]; R[i+1,2] = M[D[i][2],4]; R[i+1,3] = M[D[i][2],6];

R[i+1,4] = D[i][1]; R[i+1,5] = bearing(lat, lon, M[D[i][2],7], M[D[i][2],8]);

printf("%-36s%-18s\t%-5s\t%.1f\t\t%.0f\n", R[i+1,1], R[i+1,2], R[i+1,3], R[i+1,4], R[i+1,5]);

);

}

request(51.514669, 2.198581)</syntaxhighlight>

{{out}}

<pre>NAME COUNTRY ICAO DISTANCE BEARING

Koksijde Air Base Belgium EBFN 30.7 146

Ostend-Bruges International Airport Belgium EBOS 31.3 127

Kent International Airport United Kingdom EGMH 33.5 252

Calais-Dunkerque Airport France LFAC 34.4 196

Westkapelle heliport Belgium EBKW 42.6 105

Lympne Airport United Kingdom EGMK 51.6 240

Ursel Air Base Belgium EBUL 52.8 114

Southend Airport United Kingdom EGMC 56.2 274

Merville-Calonne Airport France LFQT 56.4 163

Wevelgem Airport Belgium EBKT 56.5 137

Midden-Zeeland Airport Netherlands EHMZ 57.3 90

Lydd Airport United Kingdom EGMD 58.0 235

RAF Wattisham United Kingdom EGUW 59.0 309

Beccles Airport United Kingdom EGSM 59.3 339

Lille/Marcq-en-Baroeul Airport France LFQO 59.7 146

Lashenden (Headcorn) Airfield United Kingdom EGKH 62.2 250

Le Touquet-Côte d'Opale Airport France LFAT 63.7 200

Rochester Airport United Kingdom EGTO 64.2 262

Lille-Lesquin Airport France LFQQ 66.2 149

Thurrock Airfield United Kingdom EGMT 68.4 272

cpu time = 183 ms, real time = 183 ms.</pre>