0
%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.

Kara
  • 6,115
  • 16
  • 50
  • 57
Pierson Sargent
  • 175
  • 2
  • 18

2 Answers2

1

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.

stachyra
  • 4,423
  • 4
  • 20
  • 34
  • I used the distance function and got an answer of 66.2987, that was without dictating the ellipsoidal Earth model though. Is there a way to syntaxually include the Earth model within the distance function? In addition, in what format is the 66.2987? Is that the delta sigma value that represents the central angle between them? If so taking that multiplied by the mean radius of the earth gives me 422390km which is incorrect. – Pierson Sargent Jan 03 '14 at 04:10
  • Output units are degrees, so to calculate distance in km in a spherical Earth model, you would take 6371 * 66.2987 * 3.14159 / 180 = 7372.07. On the website that you gave, I get a distance of 7367.11 km if I tell it to assume a spherical Earth, or I get 7381.57 (the value that you originally quoted) if I tell it to use the WGS84 ellipsoidal geoid. The spherical Earth figure (7367.11 km) is the one that should be expected to match the MATLAB result of 7372.07 km. I am not sure why you are still off by 5 km but at least now you should be hitting a lot closer to the mark than you were before. – stachyra Jan 03 '14 at 04:34
  • How do you dictate the WGS84 ellipsoidal geoid? And I guess my mistake was that I didn't divide by 180 to get radians, doh!! – Pierson Sargent Jan 03 '14 at 04:37
  • In the MATLAB distance() function, I believe you just specify 'wgs84' for the ellipsoid parameter. See [here](http://www.mathworks.com/help/map/ref/referenceellipsoidclass.html) for a complete list of allowed options. BTW, it appears that distance() only returns arc length, or "sigma" in [this notation](http://en.wikipedia.org/wiki/Vincenty%27s_formulae#Notation). If you want the true geodetic distance ("s" in same notation) I think you may have to instead download and install [this package](http://www.mathworks.com/matlabcentral/fileexchange/5379-geodetic-distance-on-wgs84-earth-ellipsoid). – stachyra Jan 03 '14 at 05:01
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.

cffk
  • 1,839
  • 1
  • 12
  • 17