8
\$\begingroup\$

Write a program that gets coordinates of two objects on Earth, and calculates how far they are from each other directly in space (a straight line through Earth) and on the surface (through the shortest route on the Earth surface). Don't take into account that polar and equatorial radii differ, consider Earth's radius 6371.0 km.

Write a program that gets 4 numbers:

  • Obj1 latitude ( [ -90.0° ; 90.0° ] )
  • Obj1 longitude ( ( -180.0° ; 180.0° ] )
  • Obj2 latitude ( [ -90.0° ; 90.0° ] )
  • Obj2 longitude ( ( -180.0° ; 180.0° ] )

and gives 2 numbers: the physical distance in space between two objects and the distance of the shortest route to go through on the Earth surface. In the context of the task Earth is a perfect sphere with no highlands/mountains or lowlands, etc..

Positive and negative values are for different hemispheres, degrees with decimal part representing minutes and seconds (e.g. -38.20 means -38°12'00", so no need to overthink). No need to check the input values, they're guaranteed to be in the range.

Hardcore. Calculate two values (in degrees) that represent:

  1. Direction of the shortest route from Obj1 to Obj2 (0° is the north, 180° is the south (-180° is represented like 180° and so on), positive degrees are on the east side (90° is the east), negative - on the west (-90° is the west)), of course before any movement, only in the starting point;
  2. From the perspective of the one of the objects, angle between its horizon and the physical position of the second object below the horizon (inclination, here always below the horizon, e.g. 90° means directly on the opposite side of Earth).

Overall output:

physical_distance shortest_route shortest_route_direction inclination_below_horizon

\$\endgroup\$
2

4 Answers 4

2
\$\begingroup\$

Jelly, 25 bytes

I;°ḊÆS,ÆẠƲP€SÆAHÆSḤ,Ʋ×⁽ßx

A dyadic Link that accepts the longitudes on the left and the latitudes, in the same order, on the right and yields the two distances, [direct, surface] in kilometres.

Try it online!

How?

Labelling the two longitudes as \$\lambda_i\$ and the two latitudes as \$\phi_i\$, the angle between the lines from each of the two points and the centre of the sphere, the central angle, is[1]:

$$\Delta\sigma = \arccos{(\cos{|\lambda_1-\lambda_2|}\cos{\phi_1}\cos{\phi_2}+\sin{\phi_1}\sin{\phi_2})}$$

Note: Since \$\cos\$ is symmetric around the y-axis, we can use \$\lambda_2-\lambda_1\$ in place of \$|\lambda_1-\lambda_2|\$.

The shortest distance between two points in Euclidean space is a straight line, so we have an isosceles triangle with the equal side lengths being the sphere's radius, \$r=6371\$ kilometres, so the direct path in space is (using \$\text{opposite}=\text{hypotenuse}\sin\alpha\$):

$$d = 2r \sin \frac{\Delta\sigma}{2}$$

The shortest distance between two points on a sphere is the minor arc of the great circle on which they lie. This is the arc projected by the angle \$\Delta\sigma\$ and thus has length:

$$s = r \Delta\sigma$$

The code takes \$(\lambda_1, \lambda_2)\$ on the left and \$(\phi_1, \phi_2)\$ on the right and produces \$(d, s)\$:

I;°ḊÆS,ÆẠƲP€SÆAHÆSḤ,Ʋ×⁽ßx - Link: Longitudes; Latitudes
I                         - forward difference {Longitudes} -> [l2 - l1]
 ;                        - concatenate {Latitudes} -> [l2 - l1, p1, p2]
  °                       - convert from degrees to radians
         Ʋ                - last four links as a monad f(X=that):
   Ḋ                      -   dequeue {X} -> [p1, p2]
    ÆS                    -   sine -> [sin(p1), sin(p2)]
       ÆẠ                 -   cosine {X} -> [cos(l2 - l1), cos(p1), cos(p2)]
      ,                   -   pair -> [[sin(p1), sin(p2)],
                                       [cos(l2 - l1), cos(p1), cos(p2)]]
          P€              - product of each
            S             - sum
             ÆA           - arccosine -> CentralAngle
                    Ʋ     - last four links as a monad - f(CentralAngle):
               H          -   halve
                ÆS        -   sine
                  Ḥ       -   double
                   ,      -   pair -> [2 × sin(CentralAngle / 2),
                                       CentralAngle]
                      ⁽ßx - 6371
                     ×    - multiply -> [6371 × 2 × sin(CentralAngle / 2),
                                         6371 × CentralAngle]
                                      = [direct, surface]
\$\endgroup\$
2
  • 1
    \$\begingroup\$ Interesting, so the direct distance through the interior doesn't actually require an \$ \arccos \$ operation, since you can use \$ d = r \sqrt { 2 - 2 cos Δσ } \$. \$\endgroup\$
    – Neil
    May 10 at 20:54
  • \$\begingroup\$ Ugh, I meant \$ d = r \sqrt { 2 - 2 \cos Δσ } \$ of course. \$\endgroup\$
    – Neil
    May 11 at 11:03
2
\$\begingroup\$

Ruby, 105 bytes

->z,p,q{a,b,c,d,e=(z<<p-q).flat_map{|j|(1i**j/=90r).rect}
p Math.acos(m=a*c*e+b*d)*r=6371,(2-m*2)**0.5*r}

4 bytes saved thanks to Dingus.

Try it online!

Ruby, 109 bytes

A function taking 3 arguments: an array of latitudes, followed by the longitudes as separate arguments.

->z,p,q{a,b,c,d,e=(z<<p-q).map{|j|(1i**j/=90.0).rect}.flatten
p Math.acos(m=a*c*e+b*d)*r=6371,(2-m*2)**0.5*r}

Try it online!

It uses essentially the same formula as Jonathan Allen's answer, which can be found at Haversine formula . Ruby syntax for trig functions requires either include Math at the beginning of the program or Math. before every instance of a trig function. Therefore extensive use is made of complex numbers, which are available by default. The formula works as follows.

The cosine of the angle between the two points is calculated as a dot product (for a pair of vectors of magnitude 1, the dot product is equal to cosine of the angle between them):

u1*u2 + v1*v2 + z1*z2

In order to do this the latitude information for each point must first be converted from an angle into vectors r1+z1i and r2+z2i.

The longitudes may be converted to u+vi coordinates r1*(x1+y1i) and r2*(x2+y2i) but by rotating the coordinate system (subtracting longitude 1 from both longitudes 1 and 2) longitude 1 becomes effectively zero, hence there is no y contribution and the x contribution becomes 1. The dot product is therefore

cos theta = (r1*1)*(r2*x2) + 0 + z1*z2 or equivalently a*1*c*e + b*d

To avoid using the longwinded Math.sin function, the sine of half the angle can be found using the identity sin(theta/2) = (1/2-(cos theta)/2)**0.5. Because we want to know double this value, the code considers coefficients 4 times larger (2-m*2) to give the required value directly.

Commented code

->z,p,q{                        #Take input as an array of latitudes, plus longitudes as separate arguments
a,b,c,d,e=                      #Set up 5 variables for the complex numbers (the 6th will be discarded)
(z<<p-q).map{|j|                #Calculate the difference between longitudes and append to array. Iterate through 3 elements
(1i**j/=90.0).rect              #For each element convert angle into rectangular coordinates by raising i to a power
}.flatten                       #Flatten the array of 3 pairs and assign to a,b,c,d,e
p Math.acos(m=a*c*e+b*d)*r=6371,#Calculate the cosine of the angle and the surface distance
(2-m*2)**0.5*r}                 #Convert the cosine into the straight line distance.
\$\endgroup\$
1
  • \$\begingroup\$ map{...}.flatten can be flat_map{...} and 90.0 can be 90r to save 4 bytes. \$\endgroup\$
    – Dingus
    May 11 at 12:02
1
\$\begingroup\$

Charcoal, 66 bytes

⊞θ⁻⊟θ⊟θ≔ΣE⪪.imag.real⁵ΠE…θ⁺²κUV⁺IXI1j∕λ⁹⁰ιθI×⁶³⁷¹⟦▷math.acosθ₂⁻²⊗θ

Attempt This Online! Link is to verbose version of code. Takes input as an array of two latitudes and two longitudes. Explanation: Port of @LevelRiverSt's Ruby answer.

⊞θ⁻⊟θ⊟θ

Subtract one longitude from the other.

≔ΣE⪪.imag.real⁵ΠE…θ⁺²κUV⁺IXI1j∕λ⁹⁰ιθ

Obtain the sine and cosine of each angle by dividing by 90, raising i to the resulting power, and extracting the real or imaginary part as appropriate. Multiply all three cosines and the sines of the latitudes together, then take the sum. This results in the cosine of the angle between the two points.

I×⁶³⁷¹⟦▷math.acosθ₂⁻²⊗θ

Take the arccosine to obtain the angle and also subtract the the doubled cosine from 2 and take the square root of the result. Multiply both values by the Earth's radius to obtain the desired output.

Would be 56 bytes in a future version of Charcoal with a function that obtains the real and/or imaginary part of a complex number:

⊞θ⁻⊟θ⊟θ≔ΣE²ΠE…θ⁺²ι§⪪XI1j∕λ⁹⁰I1jιθI×⁶³⁷¹⟦▷math.acosθ₂⁻²⊗θ

Don't Attempt This Online! because ATO doesn't have a new enough version of Charcoal yet. Link is to verbose version of code.

Would be only 52 bytes if AtIndex could split a complex number into real and imaginary parts:

⊞θ⁻⊟θ⊟θ≔ΣE²ΠE…θ⁻³ι§XI1j∕λ⁹⁰ιθI×⁶³⁷¹⟦▷math.acosθ₂⁻²⊗θ
\$\endgroup\$
0
\$\begingroup\$

JavaScript (ES6), 106 bytes

Same formulas (and same test case) as all other answers so far.

Expects (longitude1, longitude2, latitude1, latitude2).

with(Math)f=(p,q,r,s,R=12742)=>[R*sin(t=acos(cos(r*=x=PI/180)*cos(s*=x)*cos(x*=q-p)+sin(r)*sin(s))/2),R*t]

Try it online!

Formulas

Given the longitudes \$\lambda_1\$, \$\lambda_2\$ and the latitudes \$\phi_1\$, \$\phi_2\$:

$$t = \frac{1}{2}\arccos{(\cos{\phi_1}\cos{\phi_2}\cos{(\lambda_2-\lambda_1)}+\sin{\phi_1}\sin{\phi_2})}\\ d = R\sin t\\ s = Rt$$

where \$R=2\times 6371=12742\$

Degrees to radians conversions

In order to avoid extra parentheses, we use the *= operator for all conversions with carefully ordered operations:

cos(r *= x = PI / 180) * // define the conversion ratio x
                         // and apply the conversion to r
cos(s *= x) *            // apply the conversion to s
cos(x *= q - p) +        // x can now be trashed:
                         // this time we multiply x by q - p
sin(r) *                 // r has already been converted
sin(s)                   // s has already been converted
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.