The below program uses the concept of the Gall Equal-Area rectangular map projection (coopted by someone named Peters, as the Peters projection). Each choice of point takes a longitude at random from 0 to 360, and the sine of the latitude as uniformly distributed between that of the arctic circle and 1 (the sine of 90°, the North Pole). The arcsine of that number is taken as the latitude, as the equal area projection plots according to the sine of the latitude, relating to the derivative of the cosine being -sine.
The distance is found using the spherical law of cosines, and multiplying by the number of miles in a degree.
clearvars
degree=3958.8*pi/180;
dists=0; ct=0;
for trial=1:100000000
long=360*rand;
sinlat=sind(66.55)+(1-sind(66.55))*rand;
lat=asind(sinlat);
pt1=[lat,long];
long=360*rand;
sinlat=sind(66.55)+(1-sind(66.55))*rand;
lat=asind(sinlat);
pt2=[lat,long];
deltaLong=abs(pt2(2)-pt1(2));
dist=acosd(sind(pt1(1))*sind(pt2(1))+cosd(pt1(1))*cosd(pt2(1))*cosd(deltaLong))*degree;
dists=dists+dist; ct=ct+1;
end
disp(dists/ct)
The average comes out consistently as slightly under 1458 miles.
Five runs of 10 million each:
>> averageArcticDistance
1457.73576591568
>> averageArcticDistance
1457.66899209651
>> averageArcticDistance
1457.73790271195
>> averageArcticDistance
1457.88318029836
>> averageArcticDistance
1457.75403639454
Part 2:
degree=3958.8*pi/180;
dists=0; ct=0;
for trial=1:100000000
long=360*rand;
sinlat=rand*2-1;
lat=asind(sinlat);
pt1=[lat,long];
long=360*rand;
sinlat=sind(66.55)+(1-sind(66.55))*rand;
lat=asind(sinlat);
pt2=[lat,long];
deltaLong=abs(pt2(2)-pt1(2));
dist=acosd(sind(pt1(1))*sind(pt2(1))+cosd(pt1(1))*cosd(pt2(1))*cosd(deltaLong))*degree;
dists=dists+dist; ct=ct+1;
end
disp(dists/ct)
finds about 6218 or 6219 miles. This also makes sense as 6218.468 is the assumed quarter of the circumference of the earth, and, from a chosen point, it's uniformly possible to be a given distance from the original point or from its antipode. |