Great Circle Distance Problems (Spherical Vincenty Case and a simplified version)

Question

%Great Circle Distnace -- Simplified
%% 12.18993,133.45898 %% point 1 (lat/long) 
%% 14.34243,65.12750 %% point 2 (lat/long)

%%VARIABLES%%
phi_1=12.18993; %lat_1
phi_2=14.34243; %lat_2
gam_1=133.45898; %long_1
gam_2=65.12750;  %long_2
delt_gam = abs(gam_1 - gam_2); %absoulte difference in longitudes
R_Earth = 6371000; %mean radius of the earth in meters, change to FT to get distance accordingly

%%Unsimplified Great-Circle Equation -- Breaking it up into numerator and
%%denominator sections to avoid more problems -- Spherical Case of the
%%Vincenty Formula
Numer_sec1= ((cos(phi_2))*(sin(delt_gam))^2);
Numer_sec2=(((((cos(phi_1))*(sin(phi_2)))+((sin(phi_1))*(cos(phi_2))*(delt_gam))))^2);
Denom_1= (((sin(phi_1))*(sin(phi_2)))+((cos(phi_1))*(cos(phi_2))*delt_gam));

delt_sig2=atan((sqrt(Numer_sec1+Numer_sec2))/(Denom_1));

delt_GC2=R_Earth*delt_sig2;

disp(delt_GC2)

Hey guys, so currently I'm trying to get my distance between two Lat/Long Points hammered out with the Spherical Case of the Vincenty formula in MatLab. I've been referencing http://en.wikipedia.org/wiki/Great-circle_distance

And from that have created the above MatLab code. I tried the first equation that was given (a more simplified version, but to no avail either), so I'm going with the Vincenty case. Given the two lat/long points (decimal format) that are listed at the beginning of the code I've yet to calculate the correct distance in between the two points with my program. I can't seem to find out what's going on, so I'm asking you guys if there's any way you could help me figure this out.

Thank you very much in advance, and I'll be looking at the post frequently so that I can help you help me by answering any questions you may have about my code thus far.

Based on this website: http://williams.best.vwh.net/gccalc.htm the distance should be 7381.56km.

The first answer below has reminded me that I have the mapping toolbox, yet I'm not sure how to interpret the results that I'm getting, so please check the commment that I posted below. [ARCLEN, AZ] = distance(LAT1,LON1,LAT2,LON2)

this does in-fact work, but I'm not sure what I do with the arc-length or the azimuth that's produced.

Thank you and Happy New Year to all.

Answer 1

If you just want an answer for the WGS84 without programming up the algorithm and without paying for the Mapping Toolbox, download the Matlab package Geodesics on an ellipsoid of revolution. This includes an improvement in the Mapping Toolbox function, called geoddistance. To solve your problem

format long;
geoddistance(12.18993,133.45898,14.34243,65.12750)
->
7381566.23351761

The arguments to geoddistance are in degrees and the result is in meters. This does the calculation for the WGS84 ellipsoid. If you want to use a difference ellipsoid specify a 5th argument [a,e] (equatorial radius, eccentricity). (For a sphere, set e = 0; if you want to specify a prolate ellipsoid, set e to a pure imaginary. Accurate answers are returned for |e| < 0.2.)

Incidentally many of the pictures of geodesics shown in the Wikipedia article on ellipsoidal geodesics are drawn with this package.

Answer 2

The default units for trigonometric functions in MATLAB are radians. You appear to be specifying your latitudes and longitudes in degrees. Either translate to radians or else use the sind() and cosd() functions.

Or, if you happen to have the mapping toolbox installed (Mathworks does charge extra for it, however), you can just use the distance() function. The distance() function should in principle actually be the superior way, if it is available to you, because it can accept an ellipsoidal Earth model.