| This page uses content fromWikipedia. The original article was atHaversine formula. The list of authors can be seen in thepage history.As with Rosetta Code, the text of Wikipediais available under theGNU FDL. (See links for details on variance) |
Thehaversine formula is an equation important in navigation, giving great-circle distances between two points on a sphere from their longitudes and latitudes.
It is a special case of a more general formula in spherical trigonometry, thelaw of haversines, relating the sides and angles of spherical "triangles".
Implement a great-circle distance function, or use a library function, to show the great-circle distance between:
N 36°7.2',W 86°40.2' (36.12, -86.67) -and-
N 33°56.4',W 118°24.0' (33.94, -118.40)
User Kaimbridge clarified on the Talk page: -- 6371.0 km is the authalic radius based on/extracted from surface area; -- 6372.8 km is an approximation of the radius of the average circumference (i.e., the average great-elliptic or great-circle radius), where the boundaries are the meridian (6367.45 km) and the equator (6378.14 km).Using either of these values results, of course, in differing distances: 6371.0 km -> 2886.44444283798329974715782394574671655 km; 6372.8 km -> 2887.25995060711033944886005029688505340 km; (results extended for accuracy check: Given that the radii are only approximations anyways, .01' ≈ 1.0621333 km and .001" ≈ .00177 km, practical precision required is certainly no greater than about .0000001——i.e., .1 mm!)As distances are segments of great circles/circumferences, it isrecommended that the latter value (r = 6372.8 km) be used (whichmost of the given solutions have already adopted, anyways).
Most of the examples below adopted Kaimbridge's recommended value of6372.8 km for the earth radius. However, the derivation of thisellipsoidal quadratic mean radiusis wrong (the averaging over azimuth is biased). When applying theseexamples in real applications, it is better to use themean earth radius,6371 km. This value is recommended by the International Union ofGeodesy and Geophysics and it minimizes the RMS relative error between thegreat circle and geodesic distance.
F haversine(=lat1, lon1, =lat2, lon2) V r = 6372.8 V dLat = radians(lat2 - lat1) V dLon = radians(lon2 - lon1) lat1 = radians(lat1) lat2 = radians(lat2) V a = sin(dLat / 2) ^ 2 + cos(lat1) * cos(lat2) * sin(dLon / 2) ^ 2 V c = 2 * asin(sqrt(a)) R r * cprint(haversine(36.12, -86.67, 33.94, -118.40))
2887.26
DATA:X1TYPEF,Y1TYPEF,X2TYPEF,Y2TYPEF,YDTYPEF,PITYPEF,PI_180TYPEF,MINUS_1TYPEFVALUE'-1'.PI=ACOS(MINUS_1).PI_180=PI/180.LATITUDE1=36,12.LONGITUDE1=-86,67.LATITUDE2=33,94.LONGITUDE2=-118,4.X1=LATITUDE1*PI_180.Y1=LONGITUDE1*PI_180.X2=LATITUDE2*PI_180.Y2=LONGITUDE2*PI_180.YD=Y2-Y1.DISTANCE=20000/PI*ACOS(SIN(X1)*SIN(X2)+COS(X1)*COS(X2)*COS(YD)).WRITE:'Distance between given points = ',distance,'km .'.
Distance between given points = 2.884,2687 km .
withAda.Text_IO;useAda.Text_IO;withAda.Long_Float_Text_IO;useAda.Long_Float_Text_IO;withAda.Numerics.Generic_Elementary_Functions;procedureHaversine_FormulaispackageMathis newAda.Numerics.Generic_Elementary_Functions(Long_Float);useMath;-- Compute great circle distance, given latitude and longitude of two points, in radiansfunctionGreat_Circle_Distance(lat1,long1,lat2,long2:Long_Float)returnLong_FloatisEarth_Radius:constant:=6371.0;-- in kilometersa:Long_Float:=Sin(0.5*(lat2-lat1));b:Long_Float:=Sin(0.5*(long2-long1));beginreturn2.0*Earth_Radius*ArcSin(Sqrt(a*a+Cos(lat1)*Cos(lat2)*b*b));endGreat_Circle_Distance;-- convert degrees, minutes and seconds to radiansfunctionDMS_To_Radians(Deg,Min,Sec:Long_Float:=0.0)returnLong_FloatisPi_Over_180:constant:=0.017453_292519_943295_769236_907684_886127;beginreturn(Deg+Min/60.0+Sec/3600.0)*Pi_Over_180;endDMS_To_Radians;beginPut_Line("Distance in kilometers between BNA and LAX");Put(Great_Circle_Distance(DMS_To_Radians(36.0,7.2),DMS_To_Radians(86.0,40.2),-- Nashville International Airport (BNA)DMS_To_Radians(33.0,56.4),DMS_To_Radians(118.0,24.0)),-- Los Angeles International Airport (LAX)Aft=>3,Exp=>0);endHaversine_Formula;
Tested with Agena 6.2.2 Win32
scope # compute the distance between places using the Haversine formula local proc distance( th1Deg :: number, ph1Deg :: number, th2Deg :: number, ph2Deg :: number ) :: number local constant ph1 := ( ph1Deg - ph2Deg ) * radians; local constant th1 := th1Deg * radians; local constant th2 := th2Deg * radians; local constant dz := sin( th1 ) - sin( th2 ); local constant dx := cos( ph1 ) * cos( th1 ) - cos( th2 ); local constant dy := sin( ph1 ) * cos( th1 ); return arcsin( sqrt( dx * dx + dy * dy + dz * dz ) / 2 ) * 2 * 6371 end; scope local constant d := distance( 36.12, -86.67, 33.94, -118.4 ); local constant km, constant mi := round( d ), round( d / 1.609344 ); printf( "distance: %4d km (%4d mi.)", km, mi ) epocsepocs
distance: 2886 km (1794 mi.)
File: Haversine_formula.a68
#!/usr/local/bin/a68g --script #REAL r = 20 000/pi + 6.6 # km #, to rad = pi/180;PROC dist = (REAL th1 deg, ph1 deg, th2 deg, ph2 deg)REAL:( REAL ph1 = (ph1 deg - ph2 deg) * to rad, th1 = th1 deg * to rad, th2 = th2 deg * to rad, dz = sin(th1) - sin(th2), dx = cos(ph1) * cos(th1) - cos(th2), dy = sin(ph1) * cos(th1); arc sin(sqrt(dx * dx + dy * dy + dz * dz) / 2) * 2 * r);main:( REAL d = dist(36.12, -86.67, 33.94, -118.4); # Americans don't know kilometers # printf(($"dist: "g(0,1)" km ("g(0,1)" mi.)"l$, d, d / 1.609344)))dist: 2887.3 km (1794.1 mi.)
Using the mean radius value suggested in the task.
begin % compute the distance between places using the Haversine formula % real procedure arcsin( real value x ) ; arctan( x / sqrt( 1 - ( x * x ) ) ); real procedure distance ( real value th1Deg, ph1Deg, th2Deg, ph2Deg ) ; begin real ph1, th1, th2, toRad, dz, dx, dy; toRad := pi / 180; ph1 := ( ph1Deg - ph2Deg ) * toRad; th1 := th1Deg * toRad; th2 := th2Deg * toRad; dz := sin( th1 ) - sin( th2 ); dx := cos( ph1 ) * cos( th1 ) - cos( th2 ); dy := sin( ph1 ) * cos( th1 ); arcsin( sqrt( dx * dx + dy * dy + dz * dz ) / 2 ) * 2 * 6371 end distance ; begin real d; integer mi, km; d := distance( 36.12, -86.67, 33.94, -118.4 ); km := round( d ); mi := round( d / 1.609344 ); writeon( i_w := 4, s_w := 0, "distance: ", km, " km (", mi, " mi.)" ) endend.distance: 2886 km (1794 mi.)
/* fórmula de Haversine para distancias en una superficie esférica */#include<basico.h>#define MIN 60#define SEG 3600#define RADIO 6372.8#define UNAMILLA 1.609344algoritmonúmeros(lat1,lon1,lat2,lon2,dlat,dlon,millas)si'totalargumentoses(9)'// LAT1 M LON1 M LAT2 M LON2 M#basic{ lat1 = narg(2) + narg(3)/MINlon1=narg(4)+narg(5)/MINlat2=narg(6)+narg(7)/MINlon2=narg(8)+narg(9)/MIN}sinosi'totalargumentoses(13)'// LAT1 M LON1 M S LAT2 M LON2 M S#basic{ lat1 = narg(2) + narg(3)/MIN + narg(4)/SEGlon1=narg(5)+narg(6)/MIN+narg(7)/SEGlat2=narg(8)+narg(9)/MIN+narg(10)/SEGlon2=narg(11)+narg(12)/MIN+narg(13)/SEG}sinoimprimir("Modo de uso:\ndist.bas La1 M [S] Lo1 M [S] La2 M [S] Lo2 M [S]\n")términoprematurofinsi#basic{dlat=sin(radian(lat2-lat1)/2)^2dlon=sin(radian(lon2-lon1)/2)^2RADIO*(2*arcsin(sqrt(dlat+cos(radian(lat1))*cos(radian(lat2))*dlon)))}---copiaren'millas'---," km. (",millas,divididopor(UNAMILLA)," mi.)\n"decimales'2',imprimirterminar
$ hopper3 basica/dist.bas -x -o bin/distGenerating binary ‘bin/dist’...Ok! Symbols: 27Total size: 0.75 Kb$ ./bin/dist 36 7.2 86 40.2 33 56.4 118 242887.26 km. (1794.06 mi.)
setlocation;setgeo;paramcoord{iinlocation,jingeo};paramdist{iinlocation,jinlocation};data;setlocation:=BNALAX;setgeo:=LATLON;paramcoord:LATLON:=BNA36.12-86.67LAX33.94-118.4;letdist['BNA','LAX']:=2*6372.8*asin(sqrt(sin(atan(1)/45*(coord['LAX','LAT']-coord['BNA','LAT'])/2)^2+cos(atan(1)/45*coord['BNA','LAT'])*cos(atan(1)/45*coord['LAX','LAT'])*sin(atan(1)/45*(coord['LAX','LON']-coord['BNA','LON'])/2)^2));printf"The distance between the two points is approximately %f km.\n",dist['BNA','LAX'];
The distance between the two points is approximately 2887.259951 km.
r←6371hf←{(pq)←○⍺⍵÷180⋄2×rׯ1○(+/(2*⍨1○(p-q)÷2)×1(×/2○⊃¨pq))*÷2}36.12¯86.67hf33.94¯118.40
2886.44
AppleScript provides no trigonometric functions.
Here we reach through a foreign function interface to a temporary instance of a JavaScript interpreter.
useAppleScriptversion"2.4"-- Yosemite (10.10) or lateruseframework"Foundation"useframework"JavaScriptCore"usescriptingadditionspropertyjs:missing value-- haversine :: (Num, Num) -> (Num, Num) -> Numonhaversine(latLong,latLong2)set{lat,lon}tolatLongset{lat2,lon2}tolatLong2set{rlat1,rlat2,rlon1,rlon2}to¬map(myradians,{lat,lat2,lon,lon2})setdLattorlat2-rlat1setdLontorlon2-rlon1setradiusto6372.8-- kmsetasintomath("asin")setsintomath("sin")setcostomath("cos")|round|((2*radius*¬(asin's|λ|((sqrt(((sin's|λ|(dLat/2))^2)+¬(((sin's|λ|(dLon/2))^2)*¬(cos's|λ|(rlat1))*(cos's|λ|(rlat2))))))))*100)/100endhaversine-- math :: String -> Num -> Numonmath(f)scripton|λ|(x)ifmissing valueisjsthen¬setjstocurrent application'sJSContext'snew()(js'sevaluateScript:("Math."&f&"("&x&")"))'stoDouble()end|λ|endscriptendmath-------------------------- TEST ---------------------------onrunsetdistancetohaversine({36.12,-86.67},{33.94,-118.4})setjstomissing value-- Clearing a c pointer.returndistanceendrun-------------------- GENERIC FUNCTIONS ---------------------- map :: (a -> b) -> [a] -> [b]onmap(f,xs)-- The list obtained by applying f-- to each element of xs.tellmReturn(f)setlngtolengthofxssetlstto{}repeatwithifrom1tolngsetendoflstto|λ|(itemiofxs,i,xs)endrepeatreturnlstendtellendmap-- mReturn :: First-class m => (a -> b) -> m (a -> b)onmReturn(f)-- 2nd class handler function lifted into 1st class script wrapper.ifscriptisclassoffthenfelsescriptproperty|λ|:fendscriptendifendmReturn-- radians :: Float x => Degrees x -> Radians xonradians(x)(pi/180)*xendradians-- round :: a -> Inton|round|(n)roundnend|round|-- sqrt :: Num -> Numonsqrt(n)ifn≥0thenn^(1/2)elsemissing valueendifendsqrt
2887.26
radians:function[x]->x*pi//180haversine:function[src,tgt][dLat:radianstgt\0-src\0dLon:radianstgt\1-src\1lat1:radianssrc\0lat2:radianstgt\0a:addproduct@[coslat1,coslat2,sindLon/2,sindLon/2](sindLat/2)^2c:2*asinsqrtareturn6372.8*c]printhaversine@[36.12neg86.67]@[33.94,neg118.40]
2887.259950607111
#include"share/atspre_staload.hats"staload "libc/SATS/math.sats"staload _ = "libc/DATS/math.dats"staload "libc/SATS/stdio.sats"staload "libc/SATS/stdlib.sats"#define R 6372.8#define TO_RAD (3.1415926536 / 180)typedef d = doublefundist( th1: d, ph1: d, th2: d, ph2: d) : d = let val ph1 = ph1 - ph2 val ph1 = TO_RAD * ph1 val th1 = TO_RAD * th1 val th2 = TO_RAD * th2 val dz = sin(th1) - sin(th2) val dx = cos(ph1) * cos(th1) - cos(th2) val dy = sin(ph1) * cos(th1)in asin(sqrt(dx*dx + dy*dy + dz*dz)/2)*2*Rend // end of [dist]implementmain0((*void*)) = let val d = dist(36.12, ~86.67, 33.94, ~118.4); /* Americans don't know kilometers */in $extfcall(void, "printf", "dist: %.1f km (%.1f mi.)\n", d, d / 1.609344)end // end of [main0]
dist: 2887.3 km (1794.1 mi.)
MsgBox,%GreatCircleDist(36.12,33.94,-86.67,-118.40,6372.8,"km")GreatCircleDist(La1,La2,Lo1,Lo2,R,U){return,2*R*ASin(Sqrt(Hs(Rad(La2-La1))+Cos(Rad(La1))*Cos(Rad(La2))*Hs(Rad(Lo2-Lo1))))A_SpaceU}Hs(n){return,(1-Cos(n))/2}Rad(Deg){return,Deg*4*ATan(1)/180}
2887.259951 km
# syntax: GAWK -f HAVERSINE_FORMULA.AWK# converted from PythonBEGIN{distance(36.12,-86.67,33.94,-118.40)# BNA to LAXexit(0)}functiondistance(lat1,lon1,lat2,lon2,a,c,dlat,dlon){dlat=radians(lat2-lat1)dlon=radians(lon2-lon1)lat1=radians(lat1)lat2=radians(lat2)a=(sin(dlat/2))^2+cos(lat1)*cos(lat2)*(sin(dlon/2))^2c=2*atan2(sqrt(a),sqrt(1-a))printf("distance: %.4f km\n",6372.8*c)}functionradians(degree){# degrees to radiansreturndegree*(3.1415926/180.)}
distance: 2887.2599 km
100HOME:rem100CLSforGW-BASICandMSXBASIC:DELETEforMinimalBASIC110LETP=ATN(1)*4120LETD=P/180130LETM=36.12140LETK=-86.67150LETN=33.94160LETL=-118.4170LETR=6372.8180PRINT" DISTANCIA DE HAVERSINE ENTRE BNA Y LAX = ";190LETA=SIN((L-K)*D/2)200LETA=A*A210LETB=COS(M*D)*COS(N*D)220LETC=SIN((N-M)*D/2)230LETC=C*C240LETD=SQR(C+B*A)250LETE=D/SQR(1-D*D)260LETF=ATN(E)270PRINT2*R*F;"KM"280END
global radioTierra # radio de la tierra en kmradioTierra = 6372.8function Haversine(lat1, long1, lat2, long2 , radio)d_long = radians(long1 - long2)theta1 = radians(lat1)theta2 = radians(lat2)dx = cos(d_long) * cos(theta1) - cos(theta2)dy = sin(d_long) * cos(theta1)dz = sin(theta1) - sin(theta2)return asin(sqr(dx*dx + dy*dy + dz*dz) / 2) * radio * 2end functionprintprint " Distancia de Haversine entre BNA y LAX = ";print Haversine(36.12, -86.67, 33.94, -118.4, radioTierra); " km"end
Distancia de Haversine entre BNA y LAX = 2887.25994877 km.
100cls110pi=arctan(1)*4:remdefinepi=3.1415...120deg2rad=pi/180:remdefinegradosaradianes0.01745..130lat1=36.12140long1=-86.67150lat2=33.94160long2=-118.4170radio=6372.8180print" Distancia de Haversine entre BNA y LAX = ";190d_long=deg2rad*(long1-long2)200theta1=deg2rad*(lat1)210theta2=deg2rad*(lat2)220dx=cos(d_long)*cos(theta1)-cos(theta2)230dy=sin(d_long)*cos(theta1)240dz=sin(theta1)-sin(theta2)250print(asin(sqr(dx*dx+dy*dy+dz*dz)/2)*radio*2);"km"260end
Publicdeg2radAsFloat=Pi/180' define grados a radianes 0.01745..PublicradioTierraAsFloat=6372.8' radio de la tierra en kmPublicSubMain()Print"\n Distancia de Haversine entre BNA y LAX = ";Haversine(36.12,-86.67,33.94,-118.4,radioTierra);" km"EndFunctionHaversine(lat1AsFloat,long1AsFloat,lat2AsFloat,long2AsFloat,radioAsFloat)AsFloatDimd_longAsFloat=deg2rad*(long1-long2)Dimtheta1AsFloat=deg2rad*lat1Dimtheta2AsFloat=deg2rad*lat2DimdxAsFloat=Cos(d_long)*Cos(theta1)-Cos(theta2)DimdyAsFloat=Sin(d_long)*Cos(theta1)DimdzAsFloat=Sin(theta1)-Sin(theta2)ReturnASin(Sqr(dx*dx+dy*dy+dz*dz)/2)*radio*2EndFunction
Same as FreeBASIC entry.
100CLS:rem100HOMEforApplesoftBASIC:DELETEforMinimalBASIC110LETP=ATN(1)*4120LETD=P/180130LETM=36.12140LETK=-86.67150LETN=33.94160LETL=-118.4170LETR=6372.8180PRINT" DISTANCIA DE HAVERSINE ENTRE BNA Y LAX = ";190LETA=SIN((L-K)*D/2)200LETA=A*A210LETB=COS(M*D)*COS(N*D)220LETC=SIN((N-M)*D/2)230LETC=C*C240LETD=SQR(C+B*A)250LETE=D/SQR(1-D*D)260LETF=ATN(E)270PRINT2*R*F;"KM"280END
110LETP=ATN(1)*4120LETD=P/180130LETM=36.12140LETK=-86.67150LETN=33.94160LETL=-118.4170LETR=6372.8180PRINT" DISTANCIA DE HAVERSINE ENTRE BNA Y LAX = ";190LETA=SIN((L-K)*D/2)200LETA=A*A210LETB=COS(M*D)*COS(N*D)220LETC=SIN((N-M)*D/2)230LETC=C*C240LETD=SQR(C+B*A)250LETE=D/SQR(1-D*D)260LETF=ATN(E)270PRINT2*R*F;"KM"280END
TheGW-BASIC solution works without any changes.
CONSTpi=3.141593' define piCONSTradio=6372.8' radio de la tierra en kmPRINT:PRINT" Distancia de Haversine:";PRINTHaversine!(36.12,-86.67,33.94,-118.4);"km"ENDFUNCTIONHaversine!(lat1!,long1!,lat2!,long2!)deg2rad!=pi/180' define grados a radianes 0.01745..dLong!=deg2rad!*(long1!-long2!)theta1!=deg2rad!*lat1!theta2!=deg2rad!*lat2!dx!=COS(dLong!)*COS(theta1!)-COS(theta2!)dy!=SIN(dLong!)*COS(theta1!)dz!=SIN(theta1!)-SIN(theta2!)Haversine!=(SQR(dx!*dx!+dy!*dy!+dz!*dz!)/2)*radio*2ENDFUNCTION
Distancia de Haversine: 2862.63 km
100CLS110LETp=atan(1)*4120LETd=p/180130LETk=36.12140LETm=-86.67150LETl=33.94160LETn=-118.4170LETr=6372.8180PRINT" Distancia de Haversine entre BNA y LAX = ";190LETg=d*(m-n)200LETt=d*(k)210LETs=d*(l)220LETx=COS(g)*COS(t)-COS(s)230LETy=SIN(g)*COS(t)240LETz=SIN(t)-SIN(s)250PRINT(ASIN(SQR(x*x+y*y+z*z)/2)*r*2);"km"260END
DEFHaversine(lat1,long1,lat2,long2)OPTIONANGLERADIANSLETR=6372.8!radioterrestreenkm.LETdLat=RAD(lat2-lat1)LETdLong=RAD(long2-long1)LETlat1=RAD(lat1)LETlat2=RAD(lat2)LETHaversine=R*2*ASIN(SQR(SIN(dLat/2)^2+SIN(dLong/2)^2*COS(lat1)*COS(lat2)))ENDDEFPRINTPRINT"Distancia de Haversine:";Haversine(36.12,-86.67,33.94,-118.4);"km"END
Distancia de Haversine: 2887.26 km
//pi está predefinido en Yabasicdeg2rad = pi / 180 // define grados a radianes 0.01745..radioTierra = 6372.8 // radio de la tierra en kmsub Haversine(lat1, long1, lat2, long2 , radio) d_long = deg2rad * (long1 - long2) theta1 = deg2rad * lat1 theta2 = deg2rad * lat2 dx = cos(d_long) * cos(theta1) - cos(theta2) dy = sin(d_long) * cos(theta1) dz = sin(theta1) - sin(theta2) return asin(sqr(dx*dx + dy*dy + dz*dz) / 2) * radio * 2end subprint " Distancia de Haversine entre BNA y LAX = ", Haversine(36.12, -86.67, 33.94, -118.4, radioTierra), " km"end
Distancia de Haversine entre BNA y LAX = 259.478 km
Uses BBC BASIC'sMOD(array()) function which calculates the square-root of the sum of the squares of the elements of an array.
PRINT"Distance = ";FNhaversine(36.12,-86.67,33.94,-118.4)" km"ENDDEFFNhaversine(n1,e1,n2,e2)LOCALd():DIMd(2)d()=COSRAD(e1-e2)*COSRAD(n1)-COSRAD(n2),\\SINRAD(e1-e2)*COSRAD(n1),\\SINRAD(n1)-SINRAD(n2)=ASN(MOD(d())/2)*6372.8*2
Distance = 2887.25995 km
… supplied with a small POSIX shell wrapper to feed the arguments tobc. (see also)
#!/bin/sh#-# © 2021 mirabilos Ⓕ CC0 or MirBSD## implementation of Haversine GCD from public sources# output is in metres (rounded to millimetres), error < ¼%## now developed online: (user/pass public)# https://evolvis.org/plugins/scmgit/cgi-bin/gitweb.cgi?p=useful-scripts/mirkarte.git;a=blob;f=geo.sh;hb=HEADiftest"$#"-ne4;thenecho>&2"E: syntax:$0 lat1 lon1 lat2 lon2"exit1fiset-e# make GNU bc use POSIX mode and shut upBC_ENV_ARGS=-qsexportBC_ENV_ARGS# assignment of constants, variables and functions# p: multiply with to convert from degrees to radians (π/180)# r: earth radius in metres# d: distance# h: haversine intermediate# i,j: (lat,lon) point 1# x,y: (lat,lon) point 2# k: delta lat/ΔΦ# l: delta lon/Δλ# m: sin(k/2) (square root of hav(k))# n: sin(l/2) ( partial haversine )# n(x): arcsin(x)# r(x,n): round x to n decimal digits# v(x): sign (Vorzeichen)# w(x): min(1, sqrt(x)) (Wurzel)execbc-l<<-EOFscale=64define n(x) {if (x == -1) return (-2 * a(1))if (x == 1) return (2 * a(1))return (a(x / sqrt(1 - x*x)))}define v(x) {if (x < 0) return (-1)if (x > 0) return (1)return (0)}define r(x, n) {auto oo = scaleif (scale < (n + 1)) scale = (n + 1)x += v(x) * 0.5 * A^-nscale = nx /= 1scale = oreturn (x)}define w(x) {if (x >= 1) return (1)return (sqrt(x))}/* WGS84 reference ellipsoid: große Halbachse (metres), Abplattung *//* (6378137 in WGS84); this radius from Astronomical Almanac 2021 */i = 6378136.600x = 1/298.257223563/* other axis */j = i * (1 - x)/* mean radius resulting */r = (2 * i + j) / 3/* coordinates */p = (4 * a(1) / 180)i = (p * $1)j = (p * $2)x = (p * $3)y = (p * $4)/* calculation */k = (x - i)l = (y - j)m = s(k / 2)n = s(l / 2)h = ((m * m) + (c(i) * c(x) * n * n))d = 2 * r * n(w(h))r(d, 3)EOF
$ sh dist.sh 36.12 -86.67 33.94 -118.42886448.236
Note I used a more precise earth radius; the result matches that of the other implementations when choosing the same earth radius.
#include<stdio.h>#include<stdlib.h>#include<math.h>#define R 6371#define TO_RAD (3.1415926536 / 180)doubledist(doubleth1,doubleph1,doubleth2,doubleph2){doubledx,dy,dz;ph1-=ph2;ph1*=TO_RAD,th1*=TO_RAD,th2*=TO_RAD;dz=sin(th1)-sin(th2);dx=cos(ph1)*cos(th1)-cos(th2);dy=sin(ph1)*cos(th1);returnasin(sqrt(dx*dx+dy*dy+dz*dz)/2)*2*R;}intmain(){doubled=dist(36.12,-86.67,33.94,-118.4);/* Americans don't know kilometers */printf("dist: %.1f km (%.1f mi.)\n",d,d/1.609344);return0;}
publicstaticclassHaversine{publicstaticdoublecalculate(doublelat1,doublelon1,doublelat2,doublelon2){varR=6372.8;// In kilometersvardLat=toRadians(lat2-lat1);vardLon=toRadians(lon2-lon1);lat1=toRadians(lat1);lat2=toRadians(lat2);vara=Math.Sin(dLat/2)*Math.Sin(dLat/2)+Math.Sin(dLon/2)*Math.Sin(dLon/2)*Math.Cos(lat1)*Math.Cos(lat2);varc=2*Math.Asin(Math.Sqrt(a));returnR*2*Math.Asin(Math.Sqrt(a));}publicstaticdoubletoRadians(doubleangle){returnMath.PI*angle/180.0;}}voidMain(){Console.WriteLine(String.Format("The distance between coordinates {0},{1} and {2},{3} is: {4}",36.12,-86.67,33.94,-118.40,Haversine.calculate(36.12,-86.67,33.94,-118.40)));}// Returns: The distance between coordinates 36.12,-86.67 and 33.94,-118.4 is: 2887.25995060711
#define _USE_MATH_DEFINES#include<math.h>#include<iostream>conststaticdoubleEarthRadiusKm=6372.8;inlinedoubleDegreeToRadian(doubleangle){returnM_PI*angle/180.0;}classCoordinate{public:Coordinate(doublelatitude,doublelongitude):myLatitude(latitude),myLongitude(longitude){}doubleLatitude()const{returnmyLatitude;}doubleLongitude()const{returnmyLongitude;}private:doublemyLatitude;doublemyLongitude;};doubleHaversineDistance(constCoordinate&p1,constCoordinate&p2){doublelatRad1=DegreeToRadian(p1.Latitude());doublelatRad2=DegreeToRadian(p2.Latitude());doublelonRad1=DegreeToRadian(p1.Longitude());doublelonRad2=DegreeToRadian(p2.Longitude());doublediffLa=latRad2-latRad1;doubledoffLo=lonRad2-lonRad1;doublecomputation=asin(sqrt(sin(diffLa/2)*sin(diffLa/2)+cos(latRad1)*cos(latRad2)*sin(doffLo/2)*sin(doffLo/2)));return2*EarthRadiusKm*computation;}intmain(){Coordinatec1(36.12,-86.67);Coordinatec2(33.94,-118.4);std::cout<<"Distance = "<<HaversineDistance(c1,c2)<<std::endl;return0;}
(defnhaversine[{lon1:longitudelat1:latitude}{lon2:longitudelat2:latitude}](let[R6372.8; kilometersdlat(Math/toRadians(-lat2lat1))dlon(Math/toRadians(-lon2lon1))lat1(Math/toRadianslat1)lat2(Math/toRadianslat2)a(+(*(Math/sin(/dlat2))(Math/sin(/dlat2)))(*(Math/sin(/dlon2))(Math/sin(/dlon2))(Math/coslat1)(Math/coslat2)))](*R2(Math/asin(Math/sqrta)))))(haversine{:latitude36.12:longitude-86.67}{:latitude33.94:longitude-118.40});=> 2887.2599506071106
haversine=(args...) ->R=6372.8;# kmradians=args.map(deg) ->deg/180.0*Math.PIlat1=radians[0];lon1=radians[1];lat2=radians[2];lon2=radians[3]dLat=lat2-lat1dLon=lon2-lon1a=Math.sin(dLat/2)*Math.sin(dLat/2)+Math.sin(dLon/2)*Math.sin(dLon/2)*Math.cos(lat1)*Math.cos(lat2)R*2*Math.asin(Math.sqrt(a))console.loghaversine(36.12,-86.67,33.94,-118.40)
2887.2599506071124
PETSCII has the pi symbolπ in place of the ASCII tilde~; Commodore BASIC interprets this symbol as the mathematical constant.
10REM================================15REM HAVERSINE FORMULA20REM25REM 2021-09-2430REM EN.WIKIPEDIA.ORG/WIKI/HAVERSINE_FORMULA35REM40REM C64 HAS PI CONSTANT45REM X1 LONGITUDE 150REM Y1 LATITUDE 155REM X2 LONGITUDE 260REM Y2 LATITUDE 265REM70REM V1, 2021-10-02, ALVALONGO75REM ===============================100REM MAIN105DR=π/180:REMDEGREESTORADIANS110PRINTCHR$(147);CHR$(5);"HAVERSINE FORMULA"120PRINT"GREAT-CIRCLE DISTANCE"130R=6372.8:REMAVERAGEEARTHRADIUSINKILOMETERS200REM GET DATA210PRINT220INPUT"LONGITUDE 1=";X1230INPUT"LATITUDE 1=";Y1240PRINT250INPUT"LONGITUDE 2=";X2260INPUT"LATITUDE 2=";Y2270GOSUB500280PRINT290PRINT"DISTANCE=";D;"KM"300GETK$:IFK$=""THEN300310GOTO210490END500REM HAVERSINE FORMULA ------------520A=SIN((X2-X1)*DR/2)530A=A*A540B=COS(Y1*DR)*COS(Y2*DR)550C=SIN((Y2-Y1)*DR/2)560C=C*C570D=SQR(C+B*A)580E=D/SQR(1-D*D)590F=ATN(E)600D=2*R*F610RETURN
HAVERSINE FORMULAGREAT-CIRCLE DISTANCELONGITUDE 1=? -86.67LATITUDE 1=? 36.12LONGITUDE 2=? -118.40LATITUDE 2=? 33.94DISTANCE= 2887.25995 KM
(defparameter*earth-radius*6372.8)(defparameter*rad-conv*(/pi180))(defundeg->rad(x)(*x*rad-conv*))(defunhaversine(x)(expt(sin(/x2))2))(defundist-rad(lat1lng1lat2lng2)(let*((hlat(haversine(-lat2lat1)))(hlng(haversine(-lng2lng1)))(root(sqrt(+hlat(*(coslat1)(coslat2)hlng)))))(*2*earth-radius*(asinroot))))(defundist-deg(lat1lng1lat2lng2)(dist-rad(deg->radlat1)(deg->radlng1)(deg->radlat2)(deg->radlng2)))
CL-USER> (format t "~%The distance between BNA and LAX is about ~$ km.~%" (dist-deg 36.12 -86.67 33.94 -118.40))The distance between BNA and LAX is about 2887.26 km.
includeMathdefhaversine(lat1,lon1,lat2,lon2)r=6372.8# Earth radius in kilometersdeg2rad=PI/180# convert degress to radiansdLat=(lat2-lat1)*deg2raddLon=(lon2-lon1)*deg2radlat1=lat1*deg2radlat2=lat2*deg2rada=sin(dLat/2)**2+cos(lat1)*cos(lat2)*sin(dLon/2)**2c=2*asin(sqrt(a))r*cendputs"distance is#{haversine(36.12,-86.67,33.94,-118.40)} km "
distance is 2887.2599506071106 km
importstd.stdio,std.math;realhaversineDistance(inrealdth1,inrealdph1,inrealdth2,inrealdph2)purenothrow@nogc{enumrealR=6371;enumrealTO_RAD=PI/180;aliasimr=immutablereal;imrph1d=dph1-dph2;imrph1=ph1d*TO_RAD;imrth1=dth1*TO_RAD;imrth2=dth2*TO_RAD;imrdz=th1.sin-th2.sin;imrdx=ph1.cos*th1.cos-th2.cos;imrdy=ph1.sin*th1.cos;returnasin(sqrt(dx^^2+dy^^2+dz^^2)/2)*2*R;}voidmain(){writefln("Haversine distance: %.1f km",haversineDistance(36.12,-86.67,33.94,-118.4));}
Haversine distance: 2887.3 km
An alternate direct implementation of the haversine formula as shown atwikipedia. The same length, but perhaps a little more clear about what is being done.
importstd.stdio,std.math;realtoRad(inrealdegrees)purenothrow@safe@nogc{returndegrees*PI/180;}realhaversin(inrealtheta)purenothrow@safe@nogc{return(1-theta.cos)/2;}realgreatCircleDistance(inreallat1,inreallng1,inreallat2,inreallng2,inrealradius)purenothrow@safe@nogc{immutableh=haversin(lat2.toRad-lat1.toRad)+lat1.toRad.cos*lat2.toRad.cos*haversin(lng2.toRad-lng1.toRad);return2*radius*h.sqrt.asin;}voidmain(){enumrealearthRadius=6372.8L;// Average earth radius.writefln("Great circle distance: %.1f km",greatCircleDistance(36.12,-86.67,33.94,-118.4,earthRadius));}
Great circle distance: 2887.3 km
import'dart:math';classHaversine{staticfinalR=6372.8;// In kilometersstaticdoublehaversine(doublelat1,lon1,lat2,lon2){doubledLat=_toRadians(lat2-lat1);doubledLon=_toRadians(lon2-lon1);lat1=_toRadians(lat1);lat2=_toRadians(lat2);doublea=pow(sin(dLat/2),2)+pow(sin(dLon/2),2)*cos(lat1)*cos(lat2);doublec=2*asin(sqrt(a));returnR*c;}staticdouble_toRadians(doubledegree){returndegree*pi/180;}staticvoidmain(){print(haversine(36.12,-86.67,33.94,-118.40));}}
2887.2599506071106
programHaversineDemo;usesMath;functionHaversineDist(th1,ph1,th2,ph2:double):double;constdiameter=2*6372.8;vardx,dy,dz:double;beginph1:=degtorad(ph1-ph2);th1:=degtorad(th1);th2:=degtorad(th2);dz:=sin(th1)-sin(th2);dx:=cos(ph1)*cos(th1)-cos(th2);dy:=sin(ph1)*cos(th1);Result:=arcsin(sqrt(sqr(dx)+sqr(dy)+sqr(dz))/2)*diameter;end;beginWriteln('Haversine distance: ',HaversineDist(36.12,-86.67,33.94,-118.4):7:2,' km.');end.
Haversine distance: 2887.26 km.
The `spatial` extension includes the ST_Distance_Sphere() function forcomputing the haversine distance in meters between two points based onthe authalic radius (6371.0 km). The input is expected to be in WGS84(EPSG:4326) coordinates, using a [latitude, longitude] axis order.
The function ST_Distance_Spheroid() is similar but uses an ellipsoidal model and is more accurate.
#installspatial;-- if requiredloadspatial;select[lat1,lon1],[lat2,lon2],ST_Distance_sphere(ST_Point(lat1,lon1),ST_Point(lat2,lon2))as"Haversine (m)",ST_Distance_spheroid(ST_Point(lat1,lon1),ST_Point(lat2,lon2))as"ellipsoid (m)"from(select36.12aslat1,-86.67aslon1,33.94atlat2,-118.40aslon2);
┌─────────────────────────────┬─────────────────────────────┬────────────────────┬───────────────────┐│ main.list_value(lat1, lon1) │ main.list_value(lat2, lon2) │ Haversine (m) │ ellipsoid (m) ││ decimal(4,2)[] │ decimal(5,2)[] │ double │ double │├─────────────────────────────┼─────────────────────────────┼────────────────────┼───────────────────┤│ [36.12, -86.67] │ [33.94, -118.40] │ 2886444.4428379815 │ 2892776.957354406 │└─────────────────────────────┴─────────────────────────────┴────────────────────┴───────────────────┘
func dist th1 ph1 th2 ph2 . r = 6371 ph1 -= ph2 dz = sin th1 - sin th2 dx = cos ph1 * cos th1 - cos th2 dy = sin ph1 * cos th1 return 2 * r * pi / 180 * asin (sqrt (dx * dx + dy * dy + dz * dz) / 2).print dist 36.12 -86.67 33.94 -118.4
Here the radius of the Earth is taken as 6371.1 km. This is copied from the MK-61/52 solution in order to compare results for the long hop across Russia. The EDSAC solution is 5 metres out in 7000 kilometres. The other two EDSAC examples are correct to the nearest metre.
[Haversine formula, for Rosetta Code.][EDSAC, initial orders 2.] [Arrange the storage] T49K P56F [L, 35 locations: library s/r R4 to read integer.] [Can overwrite R9 once R9 is no longer needed.] T46K P92F [N, 35 locations: library s/r P7 (print integer)] T52K P128F [A, 33 locations: library subroutine T4 (arccos)] T53K P162F [B, 44 locations: 1ibrary subroutine T1 (cos)] T54K P206F [C, 10 locations: wrapper for library subroutine T1] T51K P216F [G, 24 locations: inverse haversine function] T45K P240F [H, 36 locations: haversine function] T55K P276F [V, 11 locations: constants] T47K P288F [M, 92 locations: main routine][----------------------------------------------------------------------------][Library subroutine R9 to read integers from tape at load time.][Must be placed at location 56; occupies 15 locations.] T56KGKT20FVDL8FA40DUDTFI40FA40FS39FG@S2FG23FA5@T5@E4@[----------------------------------------------------------------------------][Library subroutine M3. Prints header at load time and is then overwritten.][In this program, also sets teleprinter to figures.]PFGKIFAFRDLFUFOFE@A6FG@E8FEZPF *DISTANCES!IN!METRES@&WITH!RADIUS# .. PK [after header, blank tape and PK (WWG, 1951, page 91)][----------------------------------------------------------------------------][Constants. Load at even address.] E25K TV GK E69K [at load time, call library s/r R9 to read 35-bit contants] T#V [tell R9 where to store constants] [0] 1800000000F [180 degrees, in units of 10^-7 degree] [2] 900000000F [90 ditto] [4] 16097821018F [for changing angle unit to radians] [6] 13493037750F [pi/4 as EDSAC integer, i.e. times 2^34] [8] 6371100# [Earth radius in metres, from MK-61/52 solution] T10Z [resume normal loading] [10] RF [17-bit 1/4][----------------------------------------------------------------------------][Main routine. Load at even address.] E25K TM GK [35-bit variables, general workspace] [0] [2] [4] [6] [8] [17-bit variables] [10] [number of examples] [11] [example number] T12Z [variablss are not initialized; skip over them at load time][17-bit constants] [12] PD [1] [13] CF [ teleprinter colon (in figures mode)] [14] @F [15] &F [16] K4096F [CR, LF, null] [Enter with acc = 0] [17] A8#V TD [pass Earth's radius to printer routine] A19@ GN [finish header: print radius followed by CR, LF] O14@ O15@ A23@ GL [call library s/r R4, 0D := number of examples] AF T10@ [store (assume < 2^16)] [Here acc = example number just done (0 first time)] [27] S10@ E98@ [if all done, jump to exit] A10@ A12@ T11@ [inc and store example number] TD A11@ TF [0D := example number, extended to 35 bits] A35@ GN O13@ [print example number followed by colon] [Use library subroutine R4 to read data from tape to 0D.] [Coordinates are lat1, lon1, lat2, lon2.] A38@ GL AD U#@ T4D [lat1 to work{0} and 4D] A43@ GC A4D T2#@ [0.5*cos(lat1) to work{1}] A47@ GL AD T8#@ [lon1 to work{4}] A51@ GL AD U4#@ T4D [lat2 to work{2} and 4D] A56@ GC A4D T6#@ [0.5*cos(lat2) to work{3}] A60@ GL AD S8#@ T4D [lon2 - lon1 to 4D] A65@ GH [0.5*hav(lon2 - lon1) to 4D] H2#@ V4D LD YF T4D [times cos(lat1) to 4D] H6#@ V4D LD YF T8#@ [times cos(lat2) to work{4}] A#@ S4#@ T4D [lat1 - lat2 to 4D] A80@ GH A4D [0.5*hav(lat1 - lat2) to acc] A8#@ T4D [add product from above, sum to 4D] A85@ GG H4D [mult reg := angle/4 in radians] V8#V L1F YF TD [times radius*4, to 0D for printing] A92@ GN O14@ O15@ [print distance followed by CR, LF] A11@ E27@ [acc := example numbr; loop back] [Exit] [98] O16@ [print null to flush teleprinter buffer] ZF [stop][----------------------------------------------------------------------------][Haversine function, for Rosetta Code.][Input: 4D = angle x in range -360..360 deg,] [expressed as integer multiple of 10^-7 deg.][Output: 4D = hav(x)/2 = (1 - cos(x))/4][Load at even address. Requires library subroutines R9 and T1.] E25K TH GK A3F T35@ [plant return link as usual] [2] A4D E6@ [acc := x; jump if x >= 0] T4D [5] S4D [else x := -x] [Here with acc = abs(x). Use symmetry about 180 deg.] [6] S#V G10@ [jump if x < 180] TD SD [else acc := 180 - x] [10] A2#V [combine next 2 orders, which are shown as comments] [A #V] [add 180: either restore x, or set x := 360 - x] [Here 0 <= x <= 180. Use antisymmetry about 90 deg.] [S 2 #V] [subtract 90] E17@ [jump if x >= 90 degrees] TD A5@ T30@ [save acc, plant S4D below] AD G21@ [acc := x - 90, join common code] [17] TD A2@ T30@ [save acc, plant A4D below] SD [acc := 90 - x] [21] A2#V [either restore x, or set x := 180 - x] TD H4#V VD L4F YF [apply angle conversion factor] [Here acc = x/2, now in radians. Pass it to cosine function] T4D A28@ GB [call cosine function; returns 4D := cos(x)/2] [30] XF [(planted) either A4D or S4D] RD YF A10V [divide by 2, round, add 1/4] T4D [return haversine/2 in 4D] [35] ZF [(planted) jump back to caller][----------------------------------------------------------------------------][Inverse haversine function.][Input: 4D = hav(x)/2][Output: 4D = x/4 in radians, where 0 <= x <= pi] E25K TG GK A3F T23@ [plant return link as usual] A10V S4D LD YF [acc := 1/2 - hav(x) = cos(x)/2] G13@ [jumo if cos(x) < 0] T4D [pass cos(x)/2 to library subroutine T4 (arccis)] A8@ GA AD RD [call T4, 0D ;= x/2; then acc := x/4] E21@ [join common code] [13] TD SD T4D [here if cos(x) < 0; set 4D := -cos(x)/2] A16@ GA [call T4 (arccos); 0D := (pi - x)/2] SD RD A6#V [acc := -(pi - x)/4 + pi/4 = x/4] [21] YF T4D [round redult and return it in 4D] [23] ZF [(planted) jump back to caller][----------------------------------------------------------------------------][Wrapper for library subroutine T1 (cosine).][Input: 4D = angle x as integer multiple of 10^-7 degree.][Output: 4D = cos(x)/2] E25K TC GK A3F T9@ [plant return link as usual] H4#V V4D L4F YF T4D [4D := x/2 in radians] A7@ GB [call T1, puts cos(x)/2 into 4D] [9] ZF [jump back to caller][----------------------------------------------------------------------------][Library subroutine T1: calculates cos(x), where abs(x) <= pi/2][Input: 4D = x/2][Output: 4D = cos(x)/2][Requires library subroutine R9.] E25K TBGKT20FVDL8FA40DUDTFI40FA40FS39FG@S2FG23FA5@T5@E4@E13ZT32#@1614F73454F243967F54539267F763549741F5726623061#TZA3FT30@H4DV4DYFT4DH4DN32#@A34#@TDNDA36#@TDNDA38#@TDNDA40#@TDNDA42#@TDNDYFTDNDS4DA31@YFT4DEFIFT44Z[----------------------------------------------------------------------------][Library subroutine T4: calculates arccos][Input: 4D = cos(x)/2, assumed >= 0][Output: 0D = x/2 where 0 <= x <= pi/2] E25K TAGKA3FT28@TDA32@T6DH4DV4DL1FS29@YFE16@T4DSDA6DTDS4DT4DA6DRDYFG4@HDN30#@YFE27@T4DS4DTDEFIFT30#ZD888FT30ZO699DT32ZK4096F[----------------------------------------------------------------------------][Library subroutine P7. Prints 35-bit positive integer in 0D.][Up to 10 digits, right-justified, padded on left with spaces.] E25K TNGKA3FT26@H28#@NDYFLDT4DS27@TFH8@S8@T1FV4DAFG31@SFLDUFOFFFSFL4FT4DA1FA27@G11@T28#ZPFT27ZP1024FP610D@524D!FO30@SFL8FE22@[----------------------------------------------------------------------------][Library subroutine R4. Reads signed integer from tape at run time.][Input : None][Output : 0D = 35-bit signed integer.] E25K TLGKA3FT21@T4DH6@E11@P5DJFT6FVDL4FA4DTDI4FA4FS5@G7@S5@G20@SDTDT6FEF[----------------------------------------------------------------------------] E25K TM GK [M parameter again] E17Z [define entry point] PF [acc = 0 on entry][----------------------------------------------------------------------------][Latitude and longitude of each airport, in units of 10^-7 degree.][Read from tape by library subroutine R4. Note that sign comes after value.]3+ [number of examples]361200000+866700000- [Nashville]339400000+1184000000- [Los Angeles]547166667+205000000+ [Kaliningrad]431166667+1319000000+ [Vladivostok]499133333+62916667- [St Mary's, Scilly Isles]601919444+12436111- [Tingwall, Shetland Isles]DISTANCES IN METRESWITH RADIUS 6371100 1: 2886490 2: 7357457 3: 1186460
ELENA 6.x:
import extensions;import system'math; Haversine(lat1,lon1,lat2,lon2){ var R := 6372.8r; var dLat := (lat2 - lat1).Radian; var dLon := (lon2 - lon1).Radian; var dLat1 := lat1.Radian; var dLat2 := lat2.Radian; var a := (dLat / 2).sin() * (dLat / 2).sin() + (dLon / 2).sin() * (dLon / 2).sin() * dLat1.cos() * dLat2.cos(); ^ R * 2 * a.sqrt().arcsin()} public Program(){ console.printLineFormatted("The distance between coordinates {0},{1} and {2},{3} is: {4}", 36.12r, -86.67r, 33.94r, -118.40r, Haversine(36.12r, -86.67r, 33.94r, -118.40r))}The distance between coordinates 36.12,-86.67 and 33.94,-118.4 is: 2887.259950607
defmoduleHaversinedo@v:math.pi/180@r6372.8# km for the earth radiusdefdistance({lat1,long1},{lat2,long2})dodlat=:math.sin((lat2-lat1)*@v/2)dlong=:math.sin((long2-long1)*@v/2)a=dlat*dlat+dlong*dlong*:math.cos(lat1*@v)*:math.cos(lat2*@v)@r*2*:math.asin(:math.sqrt(a))endendbna={36.12,-86.67}lax={33.94,-118.40}IO.putsHaversine.distance(bna,lax)
2887.2599506071106
haversine:(Float,Float)->(Float,Float)->Floathaversine(lat1,lon1)(lat2,lon2)=letr=6372.8dLat=degrees(lat2-lat1)dLon=degrees(lon2-lon1)a=(sin(dLat/2))^2+(sin(dLon/2))^2*cos(degreeslat1)*cos(degreeslat2)inr*2*asin(sqrta)view=Html.div[][Html.text(toString(haversine(36.12,-86.67)(33.94,-118.4)))]
2887.2599506071106
% Implementer by Arjun Sunel-module(haversine).-export([main/0]).main()->haversine(36.12,-86.67,33.94,-118.40).haversine(Lat1,Long1,Lat2,Long2)->V=math:pi()/180,R=6372.8,% In kilometersDiff_Lat=(Lat2-Lat1)*V,Diff_Long=(Long2-Long1)*V,NLat=Lat1*V,NLong=Lat2*V,A=math:sin(Diff_Lat/2)*math:sin(Diff_Lat/2)+math:sin(Diff_Long/2)*math:sin(Diff_Long/2)*math:cos(NLat)*math:cos(NLong),C=2*math:asin(math:sqrt(A)),R*C.
2887.2599506071106
% Implemented by Claudio LariniPROGRAM HAVERSINE_DEMO!$DOUBLECONST DIAMETER=12745.6FUNCTION DEG2RAD(X) DEG2RAD=X*π/180END FUNCTIONFUNCTION RAD2DEG(X) RAD2DEG=X*180/πEND FUNCTIONPROCEDURE HAVERSINE_DIST(TH1,PH1,TH2,PH2->RES) LOCAL DX,DY,DZ PH1=DEG2RAD(PH1-PH2) TH1=DEG2RAD(TH1) TH2=DEG2RAD(TH2) DZ=SIN(TH1)-SIN(TH2) DX=COS(PH1)*COS(TH1)-COS(TH2) DY=SIN(PH1)*COS(TH1) RES=ASN(SQR(DX^2+DY^2+DZ^2)/2)*DIAMETEREND PROCEDUREBEGIN HAVERSINE_DIST(36.12,-86.67,33.94,-118.4->RES) PRINT("HAVERSINE DISTANCE: ";RES;" KM.")END PROGRAMUsing double-precision variables output is 2887.260209071741 km, while using single-precision variable output is 2887.261 Km.
Euler has a package for spherical geometry, which is used in the following code. The distances are then computed with the average radius between the two positions. Overwriting the rearth function with the given value yields the known result.
>load spherical Spherical functions for Euler. >TNA=[rad(36,7.2),-rad(86,40.2)];>LAX=[rad(33,56.4),-rad(118,24)];>esdist(TNA,LAX)->km 2886.48817482>type esdist function esdist (frompos: vector, topos: vector) r1=rearth(frompos[1]); r2=rearth(topos[1]); xfrom=spoint(frompos)*r1; xto=spoint(topos)*r2; delta=xto-xfrom; return asin(norm(delta)/(r1+r2))*(r1+r2); endfunction>function overwrite rearth (x) := 6372.8*km$>esdist(TNA,LAX)->km 2887.25995061
Binding the name HAVERSINE to the following lambda expression in the Name Manager of the Excel workbook:
(SeeLAMBDA: The ultimate Excel worksheet function)
HAVERSINE=LAMBDA(lla,LAMBDA(llb,LET(REM,"Approximate radius of Earth in km.",earthRadius,6372.8,sinHalfDeltaSquared,LAMBDA(x,SIN(x/2)^2)(RADIANS(llb-lla)),2*earthRadius*ASIN(SQRT(INDEX(sinHalfDeltaSquared,1)+(PRODUCT(COS(RADIANS(CHOOSE({1,2},INDEX(lla,1),INDEX(llb,1))))))*INDEX(sinHalfDeltaSquared,2))))))
Each of the two arguments in the example below is an Excel dynamic array of two adjacent values. The # character yields a reference to the array with the given top-left grid address.
Cell B2 is formatted to display only two decimal places.
| fx | =HAVERSINE(E2#)(H2#) | |||||||||
|---|---|---|---|---|---|---|---|---|---|---|
| A | B | C | D | E | F | G | H | I | ||
| 1 | Distance | BNA | LAX | |||||||
| 2 | 2887.26 | km | 36.12 | -86.67 | 33.94 | -118.4 | ||||
using units of measure
openSystem[<Measure>]typedeg[<Measure>]typerad[<Measure>]typekmlethaversine(θ:float<rad>)=0.5*(1.0-Math.Cos(θ/1.0<rad>))letradPerDeg=(Math.PI/180.0)*1.0<rad/deg>typepos(latitude:float<deg>,longitude:float<deg>)=memberthis.φ=latitude*radPerDegmemberthis.ψ=longitude*radPerDegletrEarth=6372.8<km>lethsDist(p1:pos)(p2:pos)=2.0*rEarth*Math.Asin(Math.Sqrt(haversine(p2.φ-p1.φ)+Math.Cos(p1.φ/1.0<rad>)*Math.Cos(p2.φ/1.0<rad>)*haversine(p2.ψ-p1.ψ)))[<EntryPoint>]letmainargv=printfn"%A"(hsDist(pos(36.12<deg>,-86.67<deg>))(pos(33.94<deg>,-118.40<deg>)))0
2887.259951
USING:arrayskernelmathmath.constantsmath.functionsmath.vectorssequences;:haversin(x--y)cos1swap-2/;:haversininv(y--x)2*1swap-acos;:haversineDist(asbs--d)[[180/pi*]map]bi@[[swap-haversin]2map][[firstcos]bi@*1swap2array]2biv.haversininvR_earth*;
(scratchpad){36.12 -86.67}{33.94 -118.4}haversineDist.2887.259950607113
Based on the Fortran and Groovy versions.
#APPTYPECONSOLEPRINT"Distance = ",Haversine(36.12,-86.67,33.94,-118.4)," km"PAUSEFUNCTIONHaversine(DegLat1ASDOUBLE,DegLon1ASDOUBLE,DegLat2ASDOUBLE,DegLon2ASDOUBLE)ASDOUBLECONSTradius=6372.8DIMdLatASDOUBLE=D2R(DegLat2-DegLat1)DIMdLonASDOUBLE=D2R(DegLon2-DegLon1)DIMlat1ASDOUBLE=D2R(DegLat1)DIMlat2ASDOUBLE=D2R(DegLat2)DIMaASDOUBLE=SIN(dLat/2)*SIN(dLat/2)+SIN(dLon/2)*SIN(dLon/2)*COS(lat1)*COS(lat2)DIMcASDOUBLE=2*ASIN(SQRT(a))RETURNradius*cENDFUNCTION
Distance = 2887.25995060711 kmPress any key to continue...
Fennel/Lua doesn't have a math.round, so we have to define a round function here.
(do;;; compute the distance between places using the Haversine formula(fndistance[th1Degph1Degth2Degph2Deg](localph1(math.rad(-ph1Degph2Deg)))(localth1(math.radth1Deg))(localth2(math.radth2Deg))(localdz(-(math.sinth1)(math.sinth2)))(localdx(-(*(math.cosph1)(math.costh1))(math.costh2)))(localdy(*(math.sinph1)(math.costh1)))(*(math.asin(/(math.sqrt(+(*dxdx)(*dydy)(*dzdz)))2))26371))(do(fnround[n](math.floor(+n0.49999999999)))(locald(distance36.12-86.6733.94-118.4))(local(kmmi)(values(roundd)(round(/d1.609344))))(print(.."distance: "km" km ("mi" mi.)"))))
distance: 2886 km (1794 mi.)
1.01 S BA = 36.12; S LA = -86.671.02 S BB = 33.94; S LB = -118.41.03 S DR = 3.1415926536 / 180; S D = 2 * 6372.81.04 S TA = (LB - LA) * DR1.05 S TB = DR * BA1.06 S TC = DR * BB1.07 S DZ = FSIN(TB) - FSIN(TC)1.08 S DX = FCOS(TA) * FCOS(TB) - FCOS(TC)1.09 S DY = FSIN(TA) * FCOS(TB)1.10 S AS = DX * DX + DY * DY + DZ * DZ1.11 S AS = FSQT(AS) / 21.12 S HDIST = D * FATN(AS / FSQT(1 - AS^2))1.13 T %6.2,"Haversine distance ",HDIST,!
Haversine distance = 2887.26
Note that FOCAL lacks a built-in arcsine function, but appendix D of the FOCAL manual shows how to compute it using arctangent and square root instead.
:s>fs>dd>f;:deg>rad174532925199433e-16f*;:differencef-deg>rad2s>ff/fsinfdupf*;:haversine( lat1 lon1 lat2 lon2 -- haversine)frotdifference( lat1 lat2 dLon^2)frotfrotfoverfover( dLon^2 lat1 lat2 lat1 lat2)fswapdifference( dLon^2 lat1 lat2 dLat^2)fswapdeg>radfcos( dLon^2 lat1 dLat^2 lat2)frotdeg>radfcosf*( dLon^2 dLat2 lat1*lat2)frotf*f+( lat1*lat2*dLon^2+dLat^2)fsqrtfasin127456s>ff*10s>ff/( haversine);36.12e-86.67e33.94e-118.40ehaversinecrf.
2887.25995060711
programexampleimplicit nonereal::dd=haversine(36.12,-86.67,33.94,-118.40)! BNA to LAXprint'(A,F9.4,A)','distance: ',d,' km'! distance: 2887.2600 kmcontains functionto_radian(degree)result(rad)! degrees to radiansreal,intent(in)::degreereal,parameter::deg_to_rad=atan(1.0)/45! exploit intrinsic atan to generate pi/180 runtime constantreal::radrad=degree*deg_to_radend functionto_radianfunctionhaversine(deglat1,deglon1,deglat2,deglon2)result(dist)! great circle distance -- adapted from Matlabreal,intent(in)::deglat1,deglon1,deglat2,deglon2real::a,c,dist,dlat,dlon,lat1,lat2real,parameter::radius=6372.8dlat=to_radian(deglat2-deglat1)dlon=to_radian(deglon2-deglon1)lat1=to_radian(deglat1)lat2=to_radian(deglat2)a=(sin(dlat/2))**2+cos(lat1)*cos(lat2)*(sin(dlon/2))**2c=2*asin(sqrt(a))dist=radius*cend functionhaversineend programexample
Here is a Free Pascal version, works in most Pascal dialects, but also note the Delphi entry that also works in Free Pascal.
programHaversineDemo;usesMath;functionHaversineDistance(constlat1,lon1,lat2,lon2:double):double;inline;constrads=pi/180;dia=2*6372.8;beginHaversineDistance:=dia*arcsin(sqrt(sqr(cos(rads*(lon1-lon2))*cos(rads*lat1)-cos(rads*lat2))+sqr(sin(rads*(lon1-lon2))*cos(rads*lat1))+sqr(sin(rads*lat1)-sin(rads*lat2)))/2);end;beginWriteln('Haversine distance between BNA and LAX: ',HaversineDistance(36.12,-86.67,33.94,-118.4):7:2,' km.');end.
' version 09-10-2016' compile with: fbc -s console' Nashville International Airport (BNA) in Nashville, TN, USA,' N 36°07.2', W 86°40.2' (36.12, -86.67)' Los Angeles International Airport (LAX) in Los Angeles, CA, USA,' N 33°56.4', W 118°24.0' (33.94, -118.40).' 6372.8 km is an approximation of the radius of the average circumference#DefinePiAtn(1)*4' define Pi = 3.1415..#Definedeg2radPi/180' define deg to rad 0.01745..#Defineearth_radius6372.8' earth radius in km.FunctionHaversine(lat1AsDouble,long1AsDouble,lat2AsDouble,_long2AsDouble,radiusAsDouble)AsDoubleDimAsDoubled_long=deg2rad*(long1-long2)DimAsDoubletheta1=deg2rad*lat1DimAsDoubletheta2=deg2rad*lat2DimAsDoubledx=Cos(d_long)*Cos(theta1)-Cos(theta2)DimAsDoubledy=Sin(d_long)*Cos(theta1)DimAsDoubledz=Sin(theta1)-Sin(theta2)ReturnAsin(Sqr(dx*dx+dy*dy+dz*dz)/2)*radius*2EndFunctionPrintPrint" Haversine distance between BNA and LAX = ";_Haversine(36.12,-86.67,33.94,-118.4,earth_radius);" km."' empty keyboard bufferWhileInkey<>"":WendPrint:Print"hit any key to end program"SleepEnd
Haversine distance between BNA and LAX = 2887.259950607111 km.
Frink has built-in constants for the radius of the earth, whether it is the mean radiusearthradius, the equatorial radiusearthradius_equatorial, or the polar radiusearthradius_polar. Below calculates the distance between the points using the haversine formula on a sphere using the mean radius, but we can do better:
haversine[theta] := (1-cos[theta])/2dist[lat1, long1, lat2, long2] := 2 earthradius arcsin[sqrt[haversine[lat2-lat1] + cos[lat1] cos[lat2] haversine[long2-long1]]]d = dist[36.12 deg, -86.67 deg, 33.94 deg, -118.40 deg]println[d-> "km"]
2886.4489734366999158 km
Note that physical constants like degrees, kilometers, and the average radius of the earth (as well as the polar and equatorial radii) are already known to Frink. Also note that units of measure are tracked throughout all calculations, and results can be displayed in a huge number of units of distance (miles, km, furlongs, chains, feet, statutemiles, etc.) by changing the final"km" to something like"miles".
However, Frink's library/sample programnavigation.frink (included in larger distributions) contains a much higher-precision calculation that uses ellipsoidal (not spherical) calculations to determine the distance on earth's geoid with far greater accuracy.
The calculations are due to:
"Direct and Inverse Solutions of Geodesics on the Ellipsoid with Applicationof Nested Equations", T.Vincenty,Survey Review XXII, 176, April 1975.http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf
There is also a slightly higher-accuracy algorithm (closer to nanometers instead of fractions of millimeters LOL):
"Algorithms for geodesics", Charles F. F. Karney,Journal of Geodesy, January 2013, Volume 87, Issue 1, pp 43-55http://link.springer.com/article/10.1007%2Fs00190-012-0578-z
use navigation.frinkd = earthDistance[36.12 deg North, 86.67 deg West, 33.94 deg North, 118.40 deg West]println[d-> "km"]
2892.7769573807044975 km
Which should be treated as the closest-to-right answer for the actual distance on the earth's geoid, based on the WGS84 geoid datum.
import math.*def haversin( theta ) = (1 - cos( theta ))/2def radians( deg ) = deg Pi/180def haversine( (lat1, lon1), (lat2, lon2) ) = R = 6372.8 h = haversin( radians(lat2 - lat1) ) + cos( radians(lat1) ) cos( radians(lat2) ) haversin( radians(lon2 - lon1) ) 2R asin( sqrt(h) )println( haversine((36.12, -86.67), (33.94, -118.40)) )
2887.259950607111
Note: The Haversine function returns an approximate theoretical value of the Great Circle Distance between two points because it does not factor the ellipsoidal shape of Earth -- fat in the middle from centrifugal force, and squashed at the ends. Navigators once relied on trigonometric functions like versine (versed sine) where angle A is 1-cos(A), and haversine (half versine) or ( 1-cos(A) ) / 2.Also, the radius of the Earth varies, at least depending on who you talk to. Here's NASA's take on it:http://nssdc.gsfc.nasa.gov/planetary/factsheet/earthfact.html
Since it was trivial, this functions returns the distance in miles and kilometers.
window 1local fn Haversine( lat1 as double, lon1 as double, lat2 as double, lon2 as double, miles as ^double, kilometers as ^double ) double deg2rad, dLat, dLon, a, c, earth_radius_miles, earth_radius_kilometers earth_radius_miles = 3959.0 // Radius of the Earth in miles earth_radius_kilometers = 6372.8 // Radius of the Earth in kilometers deg2rad = Pi / 180 // Pi is predefined in FutureBasic dLat = deg2rad * ( lat2 - lat1 ) dLon = deg2rad * ( lon2 - lon1 ) a = sin( dLat / 2 ) * sin( dLat / 2 ) + cos( deg2rad * lat1 ) * cos( deg2rad * lat2 ) * sin( dLon / 2 ) * sin( dLon / 2 ) c = 2 * asin( sqr(a) ) miles.nil# = earth_radius_miles * c kilometers.nil# = earth_radius_kilometers * cend fndouble miles, kilometersfn Haversine( 36.12, -86.67, 33.94, -118.4, @miles, @kilometers )print "Distance in miles between BNA and LAX: "; using "####.####"; miles; " miles."print "Distance in kilometers between BNA LAX: "; using "####.####"; kilometers; " km."HandleEvents
Output:
Distance in miles between BNA and LAX: 1793.6640 miles.Distance in kilometers between BNA LAX: 2887.2600 km.
packagemainimport("fmt""math")funchaversine(θfloat64)float64{return.5*(1-math.Cos(θ))}typeposstruct{φfloat64// latitude, radiansψfloat64// longitude, radians}funcdegPos(lat,lonfloat64)pos{returnpos{lat*math.Pi/180,lon*math.Pi/180}}constrEarth=6372.8// kmfunchsDist(p1,p2pos)float64{return2*rEarth*math.Asin(math.Sqrt(haversine(p2.φ-p1.φ)+math.Cos(p1.φ)*math.Cos(p2.φ)*haversine(p2.ψ-p1.ψ)))}funcmain(){fmt.Println(hsDist(degPos(36.12,-86.67),degPos(33.94,-118.40)))}
2887.2599506071097
defhaversine(lat1,lon1,lat2,lon2){defR=6372.8// In kilometersdefdLat=Math.toRadians(lat2-lat1)defdLon=Math.toRadians(lon2-lon1)lat1=Math.toRadians(lat1)lat2=Math.toRadians(lat2)defa=Math.sin(dLat/ 2) * Math.sin(dLat /2)+Math.sin(dLon/ 2) * Math.sin(dLon /2)*Math.cos(lat1)*Math.cos(lat2)defc=2*Math.asin(Math.sqrt(a))R*c}haversine(36.12,-86.67,33.94,-118.40)>2887.25995060711
importControl.Monad(join)importData.Bifunctor(bimap)importText.Printf(printf)-------------------- HAVERSINE FORMULA --------------------- The haversine of an angle.haversine::Float->Floathaversine=(^2).sin.(/2)-- The approximate distance, in kilometers,-- between two points on Earth.-- The latitude and longtitude are assumed to be in degrees.greatCircleDistance::(Float,Float)->(Float,Float)->FloatgreatCircleDistance=distDeg6371wheredistDegradiusp1p2=distRadradius(deg2radp1)(deg2radp2)distRadradius(lat1,lng1)(lat2,lng2)=(2*radius)*asin(min1.0(sqrt$haversine(lat2-lat1)+((coslat1*coslat2)*haversine(lng2-lng1))))deg2rad=joinbimap((/180).(pi*))--------------------------- TEST -------------------------main::IO()main=printf"The distance between BNA and LAX is about %0.f km.\n"(greatCircleDistancebnalax)wherebna=(36.12,-86.67)lax=(33.94,-118.40)
The distance between BNA and LAX is about 2886 km.
linkprintfproceduremain()#: Haversine formulaprintf("BNA to LAX is %d km (%d miles)\n",d:=gcdistance([36.12,-86.67],[33.94,-118.40]),d*3280/5280)# with cute km2mi conversionendproceduregcdistance(a,b)a[2]-:=b[2]every(x:=a|b)[i:=1to2]:=dtor(x[i])dz:=sin(a[1])-sin(b[1])dx:=cos(a[2])*cos(a[1])-cos(b[1])dy:=sin(a[2])*cos(a[1])returnasin(sqrt(dx*dx+dy*dy+dz*dz)/2)*2*6371end
printf.icn provides formatting
BNA to LAX is 2886 km (1793 miles)
moduleMain-- The haversine of an angle.hsin:Double->Doublehsint=letu=sin(t/2)inu*u-- The distance between two points, given by latitude and longtitude, on a-- circle. The points are specified in radians.distRad:Double->(Double,Double)->(Double,Double)->DoubledistRadradius(lat1,lng1)(lat2,lng2)=lethlat=hsin(lat2-lat1)hlng=hsin(lng2-lng1)root=sqrt(hlat+coslat1*coslat2*hlng)in2*radius*asin(min1.0root)-- The distance between two points, given by latitude and longtitude, on a-- circle. The points are specified in degrees.distDeg:Double->(Double,Double)->(Double,Double)->DoubledistDegradiusp1p2=distRadradius(deg2radp1)(deg2radp2)whered2r:Double->Doubled2rt=t*pi/180deg2rad(t,u)=(d2rt,d2ru)-- The approximate distance, in kilometers, between two points on Earth.-- The latitude and longtitude are assumed to be in degrees.earthDist:(Double,Double)->(Double,Double)->DoubleearthDist=distDeg6372.8main:IO()main=putStrLn$"The distance between BNA and LAX is about "++show(floordst)++" km."wherebna:(Double,Double)bna=(36.12,-86.67)lax:(Double,Double)lax=(33.94,-118.40)dst:Doubledst=earthDistbnalax
The distance between BNA and LAX is about 2887 km.
100 PROGRAM "Haversine.bas"110 PRINT "Haversine distance:";HAVERSINE(36.12,-86.67,33.94,-118.4);"km"120 DEF HAVERSINE(LAT1,LON1,LAT2,LON2)130 OPTION ANGLE RADIANS140 LET R=6372.8150 LET DLAT=RAD(LAT2-LAT1):LET DLON=RAD(LON2-LON1)160 LET LAT1=RAD(LAT1):LET LAT2=RAD(LAT2)170 LET HAVERSINE=R*2*ASIN(SQR(SIN(DLAT/2)^2+SIN(DLON/2)^2*COS(LAT1)*COS(LAT2)))190 END DEF
Solution:
require'trig'haversin=:0.5*1-cosRearth=:6372.8haversineDist=:Rearth*haversin^:_1@((1,*&(cos@{.))+/.*[:haversin-)&rfd
Note: J derives the inverse haversin (haversin^:_1 ) from the definition of haversin.
Example Use:
36.12_86.67haversineDist33.94_118.42887.26
publicclassHaversine{publicstaticfinaldoubleR=6372.8;// In kilometerspublicstaticdoublehaversine(doublelat1,doublelon1,doublelat2,doublelon2){lat1=Math.toRadians(lat1);lat2=Math.toRadians(lat2);doubledLat=lat2-lat1;doubledLon=Math.toRadians(lon2-lon1);doublea=Math.pow(Math.sin(dLat/2),2)+Math.pow(Math.sin(dLon/2),2)*Math.cos(lat1)*Math.cos(lat2);doublec=2*Math.asin(Math.sqrt(a));returnR*c;}publicstaticvoidmain(String[]args){System.out.println(haversine(36.12,-86.67,33.94,-118.40));}}
2887.2599506071106
functionhaversine(){varradians=Array.prototype.map.call(arguments,function(deg){returndeg/180.0*Math.PI;});varlat1=radians[0],lon1=radians[1],lat2=radians[2],lon2=radians[3];varR=6372.8;// kmvardLat=lat2-lat1;vardLon=lon2-lon1;vara=Math.sin(dLat/2)*Math.sin(dLat/2)+Math.sin(dLon/2)*Math.sin(dLon/2)*Math.cos(lat1)*Math.cos(lat2);varc=2*Math.asin(Math.sqrt(a));returnR*c;}console.log(haversine(36.12,-86.67,33.94,-118.40));
2887.2599506071124
((x,y)=>{'use strict';// haversine :: (Num, Num) -> (Num, Num) -> Numconsthaversine=([lat1,lon1],[lat2,lon2])=>{// Math lib function namesconst[pi,asin,sin,cos,sqrt,pow,round]=['PI','asin','sin','cos','sqrt','pow','round'].map(k=>Math[k]),// degrees as radians[rlat1,rlat2,rlon1,rlon2]=[lat1,lat2,lon1,lon2].map(x=>x/180*pi),dLat=rlat2-rlat1,dLon=rlon2-rlon1,radius=6372.8;// km// kmreturnround(radius*2*asin(sqrt(pow(sin(dLat/2),2)+pow(sin(dLon/2),2)*cos(rlat1)*cos(rlat2)))*100)/100;};// TESTreturnhaversine(x,y);// --> 2887.26})([36.12,-86.67],[33.94,-118.40]);
2887.26
Also works with gojq and jaq, the Go and Rust implementations of jq
For purposes of comparison, the following uses an unsatisfactory value for the Earth's radius.As noted in the task description, "in real applications, it is better to use the mean earth radius, 6371 km."
def haversine(lat1;lon1; lat2;lon2): def radians: . * (1|atan)/45; def sind: radians|sin; def cosd: radians|cos; def sq: . * .; (((lat2 - lat1)/2) | sind | sq) as $dlat | (((lon2 - lon1)/2) | sind | sq) as $dlon | 2 * 6372.8 * (( $dlat + (lat1|cosd) * (lat2|cosd) * $dlon ) | sqrt | asin) ;
Example:
haversine(36.12; -86.67; 33.94; -118.4)# 2887.2599506071106
From Javascript, ES5, except thearguments value is an Array in jsish, not an Object.
/* Haversine formula, in Jsish */functionhaversine(){varradians=arguments.map(function(deg){returndeg/180.0*Math.PI;});varlat1=radians[0],lon1=radians[1],lat2=radians[2],lon2=radians[3];varR=6372.8;// kmvardLat=lat2-lat1;vardLon=lon2-lon1;vara=Math.sin(dLat/2)*Math.sin(dLat/2)+Math.sin(dLon/2)*Math.sin(dLon/2)*Math.cos(lat1)*Math.cos(lat2);varc=2*Math.asin(Math.sqrt(a));returnR*c;};haversine(36.12,-86.67,33.94,-118.40);/*=!EXPECTSTART!=haversine(36.12, -86.67, 33.94, -118.40) ==> 2887.259950607112=!EXPECTEND!=*/
prompt$ jsish -u haversineFormula.jsi[PASS] haversineFormula.jsi
haversine(lat1,lon1,lat2,lon2)=2*6372.8*asin(sqrt(sind((lat2-lat1)/2)^2+cosd(lat1)*cosd(lat2)*sind((lon2-lon1)/2)^2))@showhaversine(36.12,-86.67,33.94,-118.4)
haversine(36.12, -86.67, 33.94, -118.4) = 2887.2599506071106
Use Unicode characters.
importjava.lang.Math.*constvalR=6372.8// in kilometersfunhaversine(lat1:Double,lon1:Double,lat2:Double,lon2:Double):Double{valλ1=toRadians(lat1)valλ2=toRadians(lat2)valΔλ=toRadians(lat2-lat1)valΔφ=toRadians(lon2-lon1)return2*R*asin(sqrt(pow(sin(Δλ/2),2.0)+pow(sin(Δφ/2),2.0)*cos(λ1)*cos(λ2)))}funmain(args:Array<String>)=println("result: "+haversine(36.12,-86.67,33.94,-118.40))
{defhaversine{defdiameter{*6372.82}}{defradians{lambda{:a}{*{/{PI}180}:a}}}{lambda{:lat1:lon1:lat2:lon2}{let{{:dLat{radians{-:lat2:lat1}}}{:dLon{radians{-:lon2:lon1}}}{:lat1{radians:lat1}}{:lat2{radians:lat2}}}{*{diameter}{asin{sqrt{+{pow{sin{/:dLat2}}2}{*{cos:lat1}{cos:lat2}{pow{sin{/:dLon2}}2}}}}}}}}}->haversine{haversine36.12-86.6733.94-118.40}->2887.2599506071106or,using{defdeg2dec{lambda{:s:w}{let{{:s{if{or{W.equal?:sW}{W.equal?:sS}}then-else+}}{:dm{S.replace°byspacein{S.replace'byin:w}}}}:s{S.get0:dm}.{round{*{/10060}{S.get1:dm}}}}}}->deg2decwecanjustwrite{haversine{deg2decN36°7.2'}{deg2decW86°40.2'}{deg2decN33°56.4'}{deg2decW118°24.0'}}->2887.2599506071106
print "Haversine distance: "; using( "####.###########", havDist( 36.12, -86.67, 33.94, -118.4)); " km."endfunction havDist( th1, ph1, th2, ph2) degtorad = acs(-1)/180 diameter = 2 * 6372.8 LgD = degtorad * (ph1 - ph2) th1 = degtorad * th1 th2 = degtorad * th2 dz = sin( th1) - sin( th2) dx = cos( LgD) * cos( th1) - cos( th2) dy = sin( LgD) * cos( th1) havDist = asn( ( dx^2 +dy^2 +dz^2)^0.5 /2) *diameterend function
Haversine distance: 2887.25995060711 km.
function radians n return n * (3.1415926 / 180)end radiansfunction haversine lat1, lng1, lat2, lng2 local radiusEarth local lat3, lng3 local lat1Rad, lat2Rad, lat3Rad local lngRad1, lngRad2, lngRad3 local haver put 6372.8 into radiusEarth put (lat2 - lat1) into lat3 put (lng2 - lng1) into lng3 put radians(lat1) into lat1Rad put radians(lat2) into lat2Rad put radians(lat3) into lat3Rad put radians(lng1) into lngRad1 put radians(lng2) into lngRad2 put radians(lng3) into lngRad3 put (sin(lat3Rad/2.0)^2) + (cos(lat1Rad)) \ * (cos(lat2Rad)) \ * (sin(lngRad3/2.0)^2) \ into haver return (radiusEarth * (2.0 * asin(sqrt(haver)))) end haversine
Test
haversine(36.12, -86.67, 33.94, -118.40)2887.259923
localfunctionhaversine(x1,y1,x2,y2)r=0.017453292519943295769236907684886127;x1=x1*r;x2=x2*r;y1=y1*r;y2=y2*r;dy=y2-y1;dx=x2-x1;a=math.pow(math.sin(dx/2),2)+math.cos(x1)*math.cos(x2)*math.pow(math.sin(dy/2),2);c=2*math.asin(math.sqrt(a));d=6372.8*c;returnd;end
Usage:
print(haversine(36.12,-86.67,33.94,-118.4));
Output:
2887.2599506071
Inputs assumed to be in radians.
distance:=(theta1,phi1,theta2,phi2)->2*6378.14*arcsin(sqrt((1-cos(theta2-theta1))/2+cos(theta1)*cos(theta2)*(1-cos(phi2-phi1))/2));
If you prefer, you can define a haversine function to clarify the definition:
haversin:=theta->(1-cos(theta))/2;distance:=(theta1,phi1,theta2,phi2)->2*6378.14*arcsin(sqrt(haversin(theta2-theta1)+cos(theta1)*cos(theta2)*haversin(phi2-phi1)));
Usage:
distance(0.6304129261, -1.512676863, 0.5923647483, -2.066469834)
2889.679287
Inputs assumed in degrees. Sin and Haversine expect arguments in radians; the built-in variable 'Degree' converts from degrees to radians.
distance[{theta1_,phi1_},{theta2_,phi2_}]:=2*6378.14ArcSin@Sqrt[Haversine[(theta2-theta1)Degree]+Cos[theta1*Degree]Cos[theta2*Degree]Haversine[(phi2-phi1)Degree]]
Usage:
distance[{36.12, -86.67}, {33.94, -118.4}]2889.68
functionrad=radians(degree)% degrees to radiansrad=degree.*pi/180;end;function[a,c,dlat,dlon]=haversine(lat1,lon1,lat2,lon2)% HAVERSINE_FORMULA.AWK - converted from AWKdlat=radians(lat2-lat1);dlon=radians(lon2-lon1);lat1=radians(lat1);lat2=radians(lat2);a=(sin(dlat./2)).^2+cos(lat1).*cos(lat2).*(sin(dlon./2)).^2;c=2.*asin(sqrt(a));arrayfun(@(x)printf("distance: %.4f km\n",6372.8*x),c);end;[a,c,dlat,dlon]=haversine(36.12,-86.67,33.94,-118.40);% BNA to LAX
distance: 2887.2600 km
dms(d,m,s):=(d+m/60+s/3600)*%pi/180$great_circle_distance(lat1,long1,lat2,long2):=12742*asin(sqrt(sin((lat2-lat1)/2)^2+cos(lat1)*cos(lat2)*sin((long2-long1)/2)^2))$/* Coordinates are found here: http://www.airport-data.com/airport/BNA/ http://www.airport-data.com/airport/LAX/ */great_circle_distance(dms(36,7,28.10),-dms(86,40,41.50),dms(33,56,32.98),-dms(118,24,29.05)),numer;/* 2886.326609413624 */
// compute the distance between places using the Haversine formuladistance=function(th1Deg,ph1Deg,th2Deg,ph2Deg)toRadians=pi/180ph1=toRadians*(ph1Deg-ph2Deg)th1=toRadians*th1Degth2=toRadians*th2Degdz=sin(th1)-sin(th2)dx=cos(ph1)*cos(th1)-cos(th2)dy=sin(ph1)*cos(th1)returnasin(sqrt(dx*dx+dy*dy+dz*dz)/2)*2*6371endfunctiond=distance(36.12,-86.67,33.94,-118.4)km=round(d)mi=round(d/1.609344)print("distance: "+km+" km ("+mi+" mi.)")
distance: 2886 km (1794 mi.)
П3->П2->П1->П0пи180/П4ИП1МГИП3МГ-ИП4*П1ИП0МГИП4*П0ИП2МГИП4*П2ИП0sinИП2sin-П8ИП1cosИП0cos*ИП2cos-П6ИП1sinИП0cos*П7ИП6x^2ИП7x^2ИП8x^2++КвКор2/arcsin2*ИП5*С/П
Input: 6371,1 as a radius of the Earth, taken as the ball, or 6367,554 as an average radius of the Earth, or 6367,562 as an approximation of the radius of the average circumference (by Krasovsky's ellipsoid) to Р5; В/Оlat1 С/Пlong1 С/Пlat2 С/Пlong2 С/П; the coordinates must be entered asdegrees,minutes (example: 46°50' as 46,5).
Test:
DELIMITER$$CREATEFUNCTIONhaversine(lat1FLOAT,lon1FLOAT,lat2FLOAT,lon2FLOAT)RETURNSFLOATNOSQLDETERMINISTICBEGINDECLARErFLOATunsignedDEFAULT6372.8;DECLAREdLatFLOATunsigned;DECLAREdLonFLOATunsigned;DECLAREaFLOATunsigned;DECLAREcFLOATunsigned;SETdLat=ABS(RADIANS(lat2-lat1));SETdLon=ABS(RADIANS(lon2-lon1));SETlat1=RADIANS(lat1);SETlat2=RADIANS(lat2);SETa=POW(SIN(dLat/2),2)+COS(lat1)*COS(lat2)*POW(SIN(dLon/2),2);SETc=2*ASIN(SQRT(a));RETURN(r*c);END$$DELIMITER;
Usage:
SELECT haversine(36.12, -86.67, 33.94, -118.4);
2887.260009765625
On Windows, you may need to usechcp 65001.
/*REXX pgm calculates distance between Nashville & Los Angles airports. */say" Nashville: north 36° 7.2', west 86° 40.2' = 36.12°, -86.67°"say"Los Angles: north 33° 56.4', west 118° 24.0' = 33.94°, -118.40°"saydist=surfaceDistance(36.12,-86.67,33.94,-118.4)saydistbef=Length(dist/1)kdist=format(dist/1,bef,2)mdist=format(dist/1.609344,bef,2)ndist=format(mdist*5280/6076.1,bef,2)say' distance between= 'kdist" kilometers,"say' or 'mdist" statute miles,"say' or 'ndist" nautical or air miles."exit/*stick a fork in it, we're done.*//*----------------------------------SURFACEDISTANCE subroutine----------*/methodradians(x)staticreturnx*Math.PI/180methodsurfaceDistance(lat1,lon1,lat2,lon2)publicstatic/*use haversine formula for dist.*/radius=6372.8/*earth's mean radius in km */dLat=radians(lat2-lat1)dLon=radians(lon2-lon1)lat1=radians(lat1)lat2=radians(lat2)a=sin(dLat/2)**2+cos(lat1)*cos(lat2)*sin(dLon/2)**2c=2*asin(sqrt(a))returnradius*cmethodcos(x)staticreturnRexxMath.cos(x)methodsin(x)staticreturnRexxMath.sin(x)methodasin(x)staticreturnRexxMath.asin(x)methodsqrt(x)staticreturnRexxMath.sqrt(x)
Same as ooRexx.
importstd/mathprochaversine(lat1,lon1,lat2,lon2:float):float=constr=6372.8# Earth radius in kilometersletdLat=degToRad(lat2-lat1)dLon=degToRad(lon2-lon1)lat1=degToRad(lat1)lat2=degToRad(lat2)a=sin(dLat/2)*sin(dLat/2)+cos(lat1)*cos(lat2)*sin(dLon/2)*sin(dLon/2)c=2*arcsin(sqrt(a))result=r*cechohaversine(36.12,-86.67,33.94,-118.40)
2887.259950607111
Works with oo2c version2
MODULEHaversines;IMPORTLRealMath,Out;PROCEDUREDistance(lat1,lon1,lat2,lon2:LONGREAL):LONGREAL;CONSTr=6372.8D0;(* Earth radius as LONGREAL *)to_radians=LRealMath.pi/180.0D0;VARd,ph1,th1,th2:LONGREAL;dz,dx,dy:LONGREAL;BEGINd:=lon1-lon2;ph1:=d*to_radians;th1:=lat1*to_radians;th2:=lat2*to_radians;dz:=LRealMath.sin(th1)-LRealMath.sin(th2);dx:=LRealMath.cos(ph1)*LRealMath.cos(th1)-LRealMath.cos(th2);dy:=LRealMath.sin(ph1)*LRealMath.cos(th1);RETURNLRealMath.arcsin(LRealMath.sqrt(LRealMath.power(dx,2.0)+LRealMath.power(dy,2.0)+LRealMath.power(dz,2.0))/2.0)*2.0*r;ENDDistance;BEGINOut.LongRealFix(Distance(36.12,-86.67,33.94,-118.4),6,10);Out.LnENDHaversines.
Output:
2887.2602975600
bundle Default { class Haversine { function : Dist(th1 : Float, ph1 : Float, th2 : Float, ph2 : Float) ~ Float { ph1 -= ph2; ph1 := ph1->ToRadians(); th1 := th1->ToRadians(); th2 := th2->ToRadians(); dz := th1->Sin()- th2->Sin(); dx := ph1->Cos() * th1->Cos() - th2->Cos(); dy := ph1->Sin() * th1->Cos(); return ((dx * dx + dy * dy + dz * dz)->SquareRoot() / 2.0)->ArcSin() * 2 * 6371.0; } function : Main(args : String[]) ~ Nil { IO.Console->Print("distance: ")->PrintLine(Dist(36.12, -86.67, 33.94, -118.4)); } }}distance: 2886.44
+(double)distanceBetweenLat1:(double)lat1lon1:(double)lon1lat2:(double)lat2lon2:(double)lon2{//degrees to radiansdoublelat1rad=lat1*M_PI/180;doublelon1rad=lon1*M_PI/180;doublelat2rad=lat2*M_PI/180;doublelon2rad=lon2*M_PI/180;//deltasdoubledLat=lat2rad-lat1rad;doubledLon=lon2rad-lon1rad;doublea=sin(dLat/2)*sin(dLat/2)+sin(dLon/2)*sin(dLon/2)*cos(lat1rad)*cos(lat2rad);doublec=2*asin(sqrt(a));doubleR=6372.8;returnR*c;}
The core calculation is fairly straightforward, but with an eye toward generality and reuse, this is how I might start:
(* Preamble -- some math, and an "angle" type which might be part of a common library. *)letpi=4.*.atan1.letradians_of_degrees=(*.)(pi/.180.)lethaversintheta=0.5*.(1.-.costheta)(* The angle type can track radians or degrees, which I'll use for automatic conversion. *)typeangle=Degoffloat|Radoffloatletas_radians=function|Degd->radians_of_degreesd|Radr->r(* Demonstrating use of a module, and record type. *)moduleLatLong=structtypet={lat:float;lng:float}letof_angleslatlng={lat=as_radianslat;lng=as_radianslng}letsubab={lat=a.lat-.b.lat;lng=a.lng-.b.lng}letdistradiusab=letd=subbainleth=haversind.lat+.haversind.lng*.cosa.lat*.cosb.latin2.*.radius*.asin(sqrth)end(* Now we can use the LatLong module to construct coordinates and calculate * great-circle distances. * NOTE radius and resulting distance are in the same measure, and units could * be tracked for this too... but who uses miles? ;) *)letearth_dist=LatLong.dist6372.8andbna=LatLong.of_angles(Deg36.12)(Deg(-86.67))andlax=LatLong.of_angles(Deg33.94)(Deg(-118.4))inearth_distbnalax;;
If the above is fed to the REPL, the last line will produce this:
# earth_dist bna lax;;- : float = 2887.25995060711102
import: math: haversine(lat1, lon1, lat2, lon2)| lat lon | lat2 lat1 - asRadian ->lat lon2 lon1 - asRadian ->lon lon 2 / sin sq lat1 asRadian cos * lat2 asRadian cos * lat 2 / sin sq + sqrt asin 2 * 6372.8 * ;haversine(36.12, -86.67, 33.94, -118.40) println
2887.25995060711
The rxmath library provides the required functions.
/*REXX pgm calculates distance between Nashville & Los Angles airports. */say" Nashville: north 36º 7.2', west 86º 40.2' = 36.12º, -86.67º"say"Los Angles: north 33º 56.4', west 118º 24.0' = 33.94º, -118.40º"saydist=surfaceDistance(36.12,-86.67,33.94,-118.4)kdist=format(dist/1,,2)/*show 2 digs past decimal point.*/mdist=format(dist/1.609344,,2)/* " " " " " " */ndist=format(mdist*5280/6076.1,,2)/* " " " " " " */say' distance between= 'kdist" kilometers,"say' or 'mdist" statute miles,"say' or 'ndist" nautical or air miles."exit/*stick a fork in it, we're done.*//*----------------------------------SURFACEDISTANCE subroutine----------*/surfaceDistance:argth1,ph1,th2,ph2/*use haversine formula for dist.*/radius=6372.8/*earth's mean radius in km */ph1=ph1-ph2x=cos(ph1)*cos(th1)-cos(th2)y=sin(ph1)*cos(th1)z=sin(th1)-sin(th2)returnradius*2*aSin(sqrt(x**2+y**2+z**2)/2)cos:ReturnRxCalcCos(arg(1))sin:ReturnRxCalcSin(arg(1))asin:ReturnRxCalcArcSin(arg(1),,'R')sqrt:ReturnRxCalcSqrt(arg(1))::requiresrxMathlibrary
Nashville: north 36º 7.2', west 86º 40.2' = 36.12º, -86.67ºLos Angles: north 33º 56.4', west 118º 24.0' = 33.94º, -118.40º distance between= 2887.26 kilometers, or 1794.06 statute miles, or 1559.00 nautical or air miles.
dist(th1, th2, ph)={ my(v=[cos(ph)*cos(th1)-cos(th2),sin(ph)*cos(th1),sin(th1)-sin(th2)]); asin(sqrt(norml2(v))/2)};distEarth(th1, ph1, th2, ph2)={ my(d=12742, deg=Pi/180); \\ Authalic diameter of the Earth d*dist(th1*deg, th2*deg, (ph1-ph2)*deg)};distEarth(36.12, -86.67, 33.94, -118.4)%1 = 2886.44444
ProgramHaversineDemo(output);usesMath;functionhaversineDist(th1,ph1,th2,ph2:double):double;constdiameter=2*6372.8;vardx,dy,dz:double;beginph1:=degtorad(ph1-ph2);th1:=degtorad(th1);th2:=degtorad(th2);dz:=sin(th1)-sin(th2);dx:=cos(ph1)*cos(th1)-cos(th2);dy:=sin(ph1)*cos(th1);haversineDist:=arcsin(sqrt(dx**2+dy**2+dz**2)/2)*diameter;end;beginwriteln('Haversine distance: ',haversineDist(36.12,-86.67,33.94,-118.4):7:2,' km.');end.
Haversine distance: 2887.26 km.
constr=6372.8;functionhaversine(lat1,lon1,lat2,lon2:real):real;beginvardLat:=degToRad(lat2-lat1);vardLon:=degToRad(lon2-lon1);lat1:=degToRad(lat1);lat2:=degToRad(lat2);vara:=sin(dLat/2)*sin(dLat/2)+cos(lat1)*cos(lat2)*sin(dLon/2)*sin(dLon/2);varc:=2*arcsin(sqrt(a));result:=r*c;end;beginhaversine(36.12,-86.67,33.94,-118.40).Println;end.
2887.25995060711
usentheoryqw/Pi/;subasin{my$x=shift;atan2($x,sqrt(1-$x*$x));}subsurfacedist{my($lat1,$lon1,$lat2,$lon2)=@_;my$radius=6372.8;my$radians=Pi()/180;;my$dlat=($lat2-$lat1)*$radians;my$dlon=($lon2-$lon1)*$radians;$lat1*=$radians;$lat2*=$radians;my$a=sin($dlat/2)**2 + cos($lat1) * cos($lat2) * sin($dlon/2)**2;my$c=2*asin(sqrt($a));return$radius*$c;}my@BNA=(36.12,-86.67);my@LAX=(33.94,-118.4);printf"Distance: %.3f km\n",surfacedist(@BNA,@LAX);
Distance: 2887.260 km
Contrary to ntheory, Math::Trig is part of the Perl core distribution.It comes with a great circle distance built-in.
useMath::Trigqw(great_circle_distance deg2rad);# Notice the 90 - latitude: phi zero is at the North Pole.# Parameter order is: LON, LATmy@BNA=(deg2rad(-86.67),deg2rad(90-36.12));my@LAX=(deg2rad(-118.4),deg2rad(90-33.94));print"Distance: ",great_circle_distance(@BNA,@LAX,6372.8)," km\n";
Distance: 2887.25995060711 km
withjavascript_semanticsconstantMER=6371,-- mean earth radius, authalic(km)--constant MER = 6372.8, -- mean earth radius, average(km)DEG_TO_RAD=PI/180functionhaversine(atomlat1,long1,lat2,long2)lat1*=DEG_TO_RADlat2*=DEG_TO_RADlong1*=DEG_TO_RADlong2*=DEG_TO_RADreturnMER*arccos(sin(lat1)*sin(lat2)+cos(lat1)*cos(lat2)*cos(long2-long1))endfunctionatomd=haversine(36.12,-86.67,33.94,-118.4)printf(1,"Distance is %.10f km (%.10f miles)\n",{d,d/1.609344})
Distance is 2886.4444428380 km (1793.5534247731 miles)
Or using the average radius of 6372.8:
Distance is 2887.2599506071 km (1794.0601578078 miles)
classPOI{private$latitude;private$longitude;publicfunction__construct($latitude,$longitude){$this->latitude=deg2rad($latitude);$this->longitude=deg2rad($longitude);}publicfunctiongetLatitude(){return$this->latitude;}publicfunctiongetLongitude(){return$this->longitude;}publicfunctiongetDistanceInMetersTo(POI$other){$radiusOfEarth=6371;// Earth's radius in kilometers.$diffLatitude=$other->getLatitude()-$this->latitude;$diffLongitude=$other->getLongitude()-$this->longitude;$a=sin($diffLatitude/2)**2+cos($this->latitude)*cos($other->getLatitude())*sin($diffLongitude/2)**2;$c=2*asin(sqrt($a));$distance=$radiusOfEarth*$c;return$distance;}}
Test:
$bna=newPOI(36.12,-86.67);// Nashville International Airport$lax=newPOI(33.94,-118.40);// Los Angeles International Airportprintf('%.2f km',$bna->getDistanceInMetersTo($lax));
2886.44 km
(scl 12)(load "@lib/math.l")(de haversine (Th1 Ph1 Th2 Ph2) (setq Ph1 (*/ (- Ph1 Ph2) pi 180.0) Th1 (*/ Th1 pi 180.0) Th2 (*/ Th2 pi 180.0) ) (let (DX (- (*/ (cos Ph1) (cos Th1) 1.0) (cos Th2)) DY (*/ (sin Ph1) (cos Th1) 1.0) DZ (- (sin Th1) (sin Th2)) ) (* `(* 2 6371) (asin (/ (sqrt (+ (* DX DX) (* DY DY) (* DZ DZ))) 2 ) ) ) ) )
Test:
(prinl "Haversine distance: " (round (haversine 36.12 -86.67 33.94 -118.4)) " km" )
Haversine distance: 2,886.444 km
test: procedure options (main); /* 12 January 2014. Derived from Fortran version */ declare d float; d = haversine(36.12, -86.67, 33.94, -118.40); /* BNA to LAX */ put edit ( 'distance: ', d, ' km') (A, F(10,3)); /* distance: 2887.2600 km */degrees_to_radians: procedure (degree) returns (float); declare degree float nonassignable; declare pi float (15) initial ( (4*atan(1.0d0)) ); return ( degree*pi/180 );end degrees_to_radians; haversine: procedure (deglat1, deglon1, deglat2, deglon2) returns (float); declare (deglat1, deglon1, deglat2, deglon2) float nonassignable; declare (a, c, dlat, dlon, lat1, lat2) float; declare radius float value (6372.8); dlat = degrees_to_radians(deglat2-deglat1); dlon = degrees_to_radians(deglon2-deglon1); lat1 = degrees_to_radians(deglat1); lat2 = degrees_to_radians(deglat2); a = (sin(dlat/2))**2 + cos(lat1)*cos(lat2)*(sin(dlon/2))**2; c = 2*asin(sqrt(a)); return ( radius*c );end haversine;end test;
distance: 2887.260 km
do-- compute the distance between places using the Haversine formulalocalfunctiondistance(th1Deg:number,ph1Deg:number,th2Deg:number,ph2Deg:number):numberlocalph1<const>=math.rad(ph1Deg-ph2Deg)localth1<const>=math.rad(th1Deg)localth2<const>=math.rad(th2Deg)localdz<const>=math.sin(th1)-math.sin(th2)localdx<const>=math.cos(ph1)*math.cos(th1)-math.cos(th2)localdy<const>=math.sin(ph1)*math.cos(th1)returnmath.asin(math.sqrt(dx*dx+dy*dy+dz*dz)/2)*2*6371enddolocald<const>=distance(36.12,-86.67,33.94,-118.4)localkm<const>,mi<const>=math.round(d),math.round(d/1.609344)print($"distance: {km} km ({mi} mi.)")endend
distance: 2886 km (1794 mi.)
Add-Type-AssemblyNameSystem.Device$BNA=New-ObjectSystem.Device.Location.GeoCoordinate36.12,-86.67$LAX=New-ObjectSystem.Device.Location.GeoCoordinate33.94,-118.40$BNA.GetDistanceTo($LAX)/1000
2888.93627213254
functionGet-GreatCircleDistance($Coord1,$Coord2){# Convert decimal degrees to radians$Lat1=$Coord1[0]/180*[math]::Pi$Long1=$Coord1[1]/180*[math]::Pi$Lat2=$Coord2[0]/180*[math]::Pi$Long2=$Coord2[1]/180*[math]::Pi# Mean Earth radius (km)$R=6371# Haversine formula$ArcLength=2*$R*[math]::Asin([math]::Sqrt([math]::Sin(($Lat1-$Lat2)/2)*[math]::Sin(($Lat1-$Lat2)/2)+[math]::Cos($Lat1)*[math]::Cos($Lat2)*[math]::Sin(($Long1-$Long2)/2)*[math]::Sin(($Long1-$Long2)/2)))return$ArcLength}$BNA=36.12,-86.67$LAX=33.94,-118.40Get-GreatCircleDistance$BNA$LAX
2886.44444283799
Up until now there is no 64bit float in Pure Data, so the result of the calculation might not be completely accurate.
#N canvas 527 1078 450 686 10;#X obj 28 427 atan2;#X obj 28 406 sqrt;#X obj 62 405 sqrt;#X obj 28 447 * 2;#X obj 62 384 -;#X msg 62 362 1 \$1;#X obj 28 339 t f f;#X obj 28 210 sin;#X obj 83 207 sin;#X obj 138 206 cos;#X obj 193 206 cos;#X obj 28 179 / 2;#X obj 83 182 / 2;#X obj 28 74 unpack f f;#X obj 28 98 t f f;#X obj 28 301 expr $f1 + ($f2 * $f3 * $f4);#X obj 28 148 deg2rad;#X obj 83 149 deg2rad;#X obj 138 148 deg2rad;#X obj 193 149 deg2rad;#X obj 28 232 t f f;#X obj 28 257 *;#X obj 83 232 t f f;#X obj 83 257 *;#X obj 83 98 t f b;#X obj 28 542 * 6372.8;#X obj 193 120 f 33.94;#X obj 28 125 - 33.94;#X msg 28 45 36.12 -86.67;#X obj 83 123 - -118.4;#X floatatom 28 577 8 0 0 0 - - -, f 8;#X connect 0 0 3 0;#X connect 1 0 0 0;#X connect 2 0 0 1;#X connect 3 0 25 0;#X connect 4 0 2 0;#X connect 5 0 4 0;#X connect 6 0 1 0;#X connect 6 1 5 0;#X connect 7 0 20 0;#X connect 8 0 22 0;#X connect 9 0 15 2;#X connect 10 0 15 3;#X connect 11 0 7 0;#X connect 12 0 8 0;#X connect 13 0 14 0;#X connect 13 1 24 0;#X connect 14 0 27 0;#X connect 14 1 18 0;#X connect 15 0 6 0;#X connect 16 0 11 0;#X connect 17 0 12 0;#X connect 18 0 9 0;#X connect 19 0 10 0;#X connect 20 0 21 0;#X connect 20 1 21 1;#X connect 21 0 15 0;#X connect 22 0 23 0;#X connect 22 1 23 1;#X connect 23 0 15 1;#X connect 24 0 29 0;#X connect 24 1 26 0;#X connect 25 0 30 0;#X connect 26 0 19 0;#X connect 27 0 16 0;#X connect 28 0 13 0;#X connect 29 0 17 0;
#DIA=2*6372.8Procedure.dHaversine(th1.d,ph1.d,th2.d,ph2.d)Definedx.d,dy.d,dz.dph1=Radian(ph1-ph2)th1=Radian(th1)th2=Radian(th2)dz=Sin(th1)-Sin(th2)dx=Cos(ph1)*Cos(th1)-Cos(th2)dy=Sin(ph1)*Cos(th1)ProcedureReturnASin(Sqr(Pow(dx,2)+Pow(dy,2)+Pow(dz,2))/2)*#DIAEndProcedureOpenConsole("Haversine distance")Print("Haversine distance: ")Print(StrD(Haversine(36.12,-86.67,33.94,-118.4),7)+" km.")Input()
Haversine distance: 2887.2599506 km.
frommathimportradians,sin,cos,sqrt,asindefhaversine(lat1,lon1,lat2,lon2):R=6372.8# Earth radius in kilometersdLat=radians(lat2-lat1)dLon=radians(lon2-lon1)lat1=radians(lat1)lat2=radians(lat2)a=sin(dLat/2)**2+cos(lat1)*cos(lat2)*sin(dLon/2)**2c=2*asin(sqrt(a))returnR*c>>>haversine(36.12,-86.67,33.94,-118.40)2887.2599506071106>>>
SCREEN _NEWIMAGE(800, 100, 32)'*** Units: K=kilometers M=miles N=nautical milesDIM UNIT AS STRINGDIM Distance AS STRINGDIM Result AS DOUBLEDIM ANSWER AS DOUBLE'*** Change the To/From Latittude/Logitudes for your run'*** LAT/LON for Nashville International Airport (BNA)lat1 = 36.12Lon1 = -86.67'*** LAT/LONG for Los Angeles International Airport (LAX)Lat2 = 33.94Lon2 = -118.40'*** Initialize ValuesUNIT = "K"Distance = ""'Radius = 6378.137Radius = 6372.8'*** Calculate distance using Haversine Functionlat1 = (lat1 * _PI / 180)Lon1 = (Lon1 * _PI / 180)Lat2 = (Lat2 * _PI / 180)Lon2 = (Lon2 * _PI / 180)DLon = Lon1 - Lon2ANSWER = _ACOS(SIN(lat1) * SIN(Lat2) + COS(lat1) * COS(Lat2) * COS(DLon)) * Radius'*** Adjust Answer based on Distance Unit (kilometers, miles, nautical miles)SELECT CASE UNIT CASE "M" Result = ANSWER * 0.621371192 Distance = "miles" CASE "N" Result = ANSWER * 0.539956803 Distance = "nautical miles" CASE ELSE Result = ANSWER Distance = "kilometers"END SELECT'*** Change PRINT statement with your labels for FROM/TO locationsPRINT "The distance from Nashville International to Los Angeles International in "; Distance;PRINT USING " is: ##,###.##"; Result;PRINT "."END
dms_to_rad<-function(d,m,s)(d+m/60+s/3600)*pi/180# Volumetric mean radius is 6371 km, see http://nssdc.gsfc.nasa.gov/planetary/factsheet/earthfact.html# The diameter is thus 12742 kmgreat_circle_distance<-function(lat1,long1,lat2,long2){a<-sin(0.5*(lat2-lat1))b<-sin(0.5*(long2-long1))12742*asin(sqrt(a*a+cos(lat1)*cos(lat2)*b*b))}# Coordinates are found here:# http://www.airport-data.com/airport/BNA/# http://www.airport-data.com/airport/LAX/great_circle_distance(dms_to_rad(36,7,28.10),dms_to_rad(86,40,41.50),# Nashville International Airport (BNA)dms_to_rad(33,56,32.98),dms_to_rad(118,24,29.05))# Los Angeles International Airport (LAX)# Output: 2886.327
Almost the same as the Scheme version.
#langracket(requiremath)(defineearth-radius6371)(define(distancelat1long1lat2long2)(define(hab)(sqr(sin(/(-ba)2))))(*2earth-radius(asin(sqrt(+(hlat1lat2)(*(coslat1)(coslat2)(hlong1long2)))))))(define(deg-to-raddms)(*(/pi180)(+d(/m60)(/s3600))))(distance(deg-to-rad367.20)(deg-to-rad8640.20)(deg-to-rad3356.40)(deg-to-rad11824.00))
2886.444442837984
(formerly Perl 6)
classEarthPoint {has$.lat;# latitudehas$.lon;# longitudehas$earth_radius =6371;# mean earth radiushas$radian_ratio =pi /180;# accessors for radiansmethodlatR {$.lat *$radian_ratio }methodlonR {$.lon *$radian_ratio }methodhaversine-dist(EarthPoint$p) {myEarthPoint$arc .=new(lat =>$!lat -$p.lat,lon =>$!lon -$p.lon );my$a =sin($arc.latR/2) **2 +sin($arc.lonR/2) **2 *cos($.latR) *cos($p.latR);my$c =2 *asin(sqrt($a) );return$earth_radius *$c; }}myEarthPoint$BNA .=new(lat =>36.12,lon => -86.67);myEarthPoint$LAX .=new(lat =>33.94,lon => -118.4);say$BNA.haversine-dist($LAX);# 2886.44444099822
define PI -1 acosdefine toRadians use $degree $degree PI * 180 /define haversine use $lat1, $lon1, $lat2, $lon2 6372.8 as $R # In kilometers $lat2 $lat1 - toRadians as $dLat $lon2 $lon1 - toRadians as $dLon $lat1 toRadians as $lat1 $lat2 toRadians as $lat2 $dLat 2 / sin $dLat 2 / sin * $dLon 2 / sin $dLon 2 / sin * $lat1 cos * $lat2 cos * + as $a $a sqrt asin 2 * as $c $R $c *} -118.40 33.94 -86.67 36.12 haversine "haversine: %.15g\n" print
haversine: 2887.25995060711
The use of normalization for angles isn't required for the Haversine formula, but those normalization functions were included
herein anyway (to support normalization of input arguments to the trigonometric functions for the general case).
/*REXX program calculates the distance between Nashville and Los Angles airports.*/callpi;numericdigitslength(pi)%2/*use half of PI dec. digits for output*/say" Nashville: north 36º 7.2', west 86º 40.2' = 36.12º, -86.67º"say" Los Angles: north 33º 56.4', west 118º 24.0' = 33.94º, -118.40º"@using_radius='using the mean radius of the earth as '/*a literal for SAY.*/radii.=.;radii.1=6372.8;radii.2=6371/*mean radii of the earth in kilometers*/say;m=1/0.621371192237/*M: one statute mile in " */doradius=1whileradii.radius\==./*calc. distance using specific radii. */d=surfaceDist(36.12,-86.67,33.94,-118.4,radii.radius);saysaycenter(@using_radiusradii.radius' kilometers',75,'─')say' Distance between: 'format(d/1,,2)" kilometers,"say' or 'format(d/m,,2)" statute miles,"say' or 'format(d/m*5280/6076.1,,2)" nautical (or air miles)."end/*radius*//*show──┘ 2 dec. digs past dec. point*/exit/*stick a fork in it, we're all done. *//*──────────────────────────────────────────────────────────────────────────────────────*/surfaceDist:parseargth1,ph1,th2,ph2,r/*use haversine formula for distance.*/numericdigitsdigits()*2/*double number of decimal digits used.*/ph1=d2r(ph1-ph2)/*convert degrees ──► radians & reduce.*/th1=d2r(th1)/* " " " " " */th2=d2r(th2)/* " " " " " */cosTH1=cos(th1)/*compute a shortcut (it's used twice).*/x=cos(ph1)*cosTH1-cos(th2)/* " X coordinate. */y=sin(ph1)*cosTH1/* " Y " */z=sin(th1)-sin(th2)/* " Z " */returnAsin(sqrt(x*x+y*y+z*z)*.5)*r*2/*compute the arcsin and return value. *//*──────────────────────────────────────────────────────────────────────────────────────*/Acos:returnpi()*.5-aSin(arg(1))/*calculate the ArcCos of an argument. */d2d:returnarg(1)//360/*normalize degrees to a unit circle. */d2r:returnr2r(arg(1)*pi()/180)/*normalize and convert deg ──► radians*/r2d:returnd2d((arg(1)*180/pi()))/*normalize and convert rad ──► degrees*/r2r:returnarg(1)//(pi()*2)/*normalize radians to a unit circle. */pi:pi=3.141592653589793238462643383279502884197169399375105820975;returnpi/*──────────────────────────────────────────────────────────────────────────────────────*/Asin:procedure;parseargx1z1o1p;a=abs(x);aa=a*aifa>=sqrt(2)*.5thenreturnsign(x)*Acos(sqrt(1-aa))doj=2by2untilp=z;p=z;o=o*aa*(j-1)/j;z=z+o/(j+1)end/*j*/;returnz/* [↑] compute until no more noise. *//*──────────────────────────────────────────────────────────────────────────────────────*/cos:procedure;parseargx;x=r2r(x);a=abs(x);Hpi=pi*.5numericfuzzmin(6,digits()-3);ifa=pithenreturn-1ifa=Hpi|a=Hpi*3thenreturn0;ifa=pi/3thenreturn.5ifa=pi*2/3thenreturn-.5;q=x*x;p=1;z=1;_=1dok=2by2;_=-_*q/(k*(k-1));z=z+_;ifz=pthenleave;p=z;end;returnz/*──────────────────────────────────────────────────────────────────────────────────────*/sin:procedure;parseargx;x=r2r(x);numericfuzzmin(5,digits()-3)ifabs(x)=pithenreturn0;q=x*x;p=x;z=x;_=xdok=2by2;_=-_*q/(k*(k+1));z=z+_;ifz=pthenleave;p=z;end;returnz/*──────────────────────────────────────────────────────────────────────────────────────*/sqrt:procedure;parseargx;ifx=0thenreturn0;d=digits();m.=9;numericform;h=d+6numericdigits;parsevalueformat(x,2,1,,0)'E0'withg"E"_.;g=g*.5'e'_%2doj=0whileh>9;m.j=h;h=h%2+1;end/*j*/dok=j+5to0by-1;numericdigitsm.k;g=(g+x/g)*.5;end/*k*/;returng
REXX doesn't have most of the higher math functions, so they are included here (above) as subroutines (functions).
╔════════════════════════════════════════════════════════════════════════╗ ║ A note on built─in functions: REXX doesn't have a lot of mathematical ║ ║ or (particularly) trigonometric functions, so REXX programmers have ║ ║ to write their own. Usually, this is done once, or most likely, one ║ ║ is borrowed from another program. Knowing this, the one that is used ║ ║ has a lot of boilerplate in it. ║ ║ ║ ║ Programming note: the "general 1─liner" subroutines are taken from ║ ║ other programs that I wrote, but I broke up their one line of source ║ ║ so it can be viewed without shifting the viewing window. ║ ║ ║ ║ The pi constant (as used here) is actually a much more robust ║ ║ function and will return up to one million digits in the real version. ║ ║ ║ ║ One bad side effect is that, like a automobile without a hood, you see ║ ║ all the dirty stuff going on. Also, don't visit a sausage factory. ║ ╚════════════════════════════════════════════════════════════════════════╝
Nashville: north 36º 7.2', west 86º 40.2' = 36.12º, -86.67º Los Angles: north 33º 56.4', west 118º 24.0' = 33.94º, -118.40º─────────using the mean radius of the earth as 6372.8 kilometers───────── Distance between: 2887.26 kilometers, or 1794.06 statute miles, or 1559.00 nautical (or air miles).──────────using the mean radius of the earth as 6371 kilometers────────── Distance between: 2886.44 kilometers, or 1793.55 statute miles, or 1558.56 nautical (or air miles).
decimals(8)see haversine(36.12, -86.67, 33.94, -118.4) + nlfunc haversine x1, y1, x2, y2 r=0.01745 x1= x1*r x2= x2*r y1= y1*r y2= y2*r dy = y2-y1 dx = x2-x1 a = pow(sin(dx/2),2) + cos(x1) * cos(x2) * pow(sin(dy/2),2) c = 2 * asin(sqrt(a)) d = 6372.8 * c return d
| Code | Comments |
|---|---|
≪ ROT - 2 / DEG SIN SQ OVER COS * 3 PICK COS * ROT ROT - 2 / SIN SQ + √ RAD ASIN 6372.8 * 2 *≫ 'AHAV' STO | ( lat1 lon1 lat2 lon2 -- distance ) Start by the end of the formula, in degree mode Switch to radian mode to compute Arcsin |
The following line of command delivers what is required:
36.12 -86.67 33.94 -118.4 AHAV
Due to the uncertainty in values of Earth radius and airports coordinates, the result shall be announced as 2887 ± 1 km even if the calculation provides many digits after the decimal point
1: 2887.25995061
includeMathRadius=6372.8# rough radius of the Earth, in kilometersdefspherical_distance(start_coords,end_coords)lat1,long1=deg2rad*start_coordslat2,long2=deg2rad*end_coords2*Radius*asin(sqrt(sin((lat2-lat1)/2)**2+cos(lat1)*cos(lat2)*sin((long2-long1)/2)**2))enddefdeg2rad(lat,long)[lat*PI/180,long*PI/180]endbna=[36.12,-86.67]lax=[33.94,-118.4]puts"%.1f"%spherical_distance(bna,lax)
2887.3
Alternatively:
includeMathdefhaversine(lat1,lon1,lat2,lon2)r=6372.8# Earth radius in kilometersdeg2rad=PI/180# convert degress to radiansdLat=(lat2-lat1)*deg2raddLon=(lon2-lon1)*deg2radlat1=lat1*deg2radlat2=lat2*deg2rada=sin(dLat/2)**2+cos(lat1)*cos(lat2)*sin(dLon/2)**2c=2*asin(sqrt(a))r*cendputs"distance is#{haversine(36.12,-86.67,33.94,-118.40)} km "
distance is 2887.2599506071106 km
D2R = atn(1)/45 diam = 2 * 6372.8Lg1m2 = ((-86.67)-(-118.4)) * D2RLt1 = 36.12 * D2R ' degrees to radLt2 = 33.94 * D2R dz = sin(Lt1) - sin(Lt2) dx = cos(Lg1m2) * cos(Lt1) - cos(Lt2) dy = sin(Lg1m2) * cos(Lt1) hDist = asn((dx^2 + dy^2 + dz^2)^0.5 /2) * diamprint "Haversine distance: ";using("####.#############",hDist);" km." 'Tips: ( 36 deg 7 min 12 sec ) = print 36+(7/60)+(12/3600). Produces: 36.12 deg. ' ' http://maps.google.com ' Search 36.12,-86.67 ' Earth. ' Center the pin, zoom airport. ' Directions (destination). ' 36.12.-86.66999 ' Distance is 35.37 inches.Output
Haversine distance: 2887.2599506071104 km.
structPoint{lat:f64,lon:f64,}fnhaversine(origin:Point,destination:Point)->f64{constR:f64=6372.8;letlat1=origin.lat.to_radians();letlat2=destination.lat.to_radians();letd_lat=lat2-lat1;letd_lon=(destination.lon-origin.lon).to_radians();leta=(d_lat/2.0).sin().powi(2)+(d_lon/2.0).sin().powi(2)*lat1.cos()*lat2.cos();letc=2.0*a.sqrt().asin();R*c}#[cfg(test)]modtest{usesuper::{Point,haversine};#[test]fntest_haversine(){letorigin:Point=Point{lat:36.12,lon:-86.67};letdestination:Point=Point{lat:33.94,lon:-118.4};letd:f64=haversine(origin,destination);println!("Distance: {} km ({} mi)",d,d/1.609344);assert_eq!(d,2887.2599506071106);}}
Output
Distance: 2887.2599506071106 km (1794.060157807846 mi)
options minoperator;%macro haver(lat1, long1, lat2, long2, type=D, dist=K);%if%upcase(&type)in (D DEG DEGREE DEGREES)%then%do;%let convert = constant('PI')/180;%end;%else%if%upcase(&type)in (R RAD RADIAN RADIANS)%then%do;%let convert =1;%end;%else%do;%put ERROR - Enter RADIANSor DEGREES for type.;%goto exit;%end;%if%upcase(&dist)in (M MILE MILES)%then%do;%let distrat =1.609344;%end;%else%if%upcase(&dist)in (K KM KILOMETER KILOMETERS)%then%do;%let distrat =1;%end;%else%do;%put ERROR - Enter Mon KM for dist;%goto exit;%end;data_null_;convert =&convert;lat1 =&lat1 * convert;lat2 =&lat2 * convert;long1 =&long1 * convert;long2 =&long2 * convert;diff1 = lat2 - lat1;diff2 = long2 - long1;part1 =sin(diff1/2)**2;part2 =cos(lat1)*cos(lat2);part3 =sin(diff2/2)**2;root =sqrt(part1 + part2*part3);dist =2 *6372.8 /&distrat *arsin(root);put"Distance is " dist"%upcase(&dist)";run;%exit:%mend;%haver(36.12, -86.67,33.94, -118.40);
Distance is 2887.2599506 K
importmath._objectHaversine{valR=6372.8//radius in kmdefhaversine(lat1:Double,lon1:Double,lat2:Double,lon2:Double)={valdLat=(lat2-lat1).toRadiansvaldLon=(lon2-lon1).toRadiansvala=pow(sin(dLat/2),2)+pow(sin(dLon/2),2)*cos(lat1.toRadians)*cos(lat2.toRadians)valc=2*asin(sqrt(a))R*c}defmain(args:Array[String]):Unit={println(haversine(36.12,-86.67,33.94,-118.40))}}
2887.2599506071106
(defineearth-radius6371)(definepi(acos-1))(define(distancelat1long1lat2long2)(define(hab)(expt(sin(/(-ba)2))2))(*2earth-radius(asin(sqrt(+(hlat1lat2)(*(coslat1)(coslat2)(hlong1long2)))))))(define(deg-to-raddms)(*(/pi180)(+d(/m60)(/s3600))))(distance(deg-to-rad367.20)(deg-to-rad8640.20)(deg-to-rad3356.40)(deg-to-rad11824.00)); 2886.444442837984
$ include "seed7_05.s7i"; include "float.s7i"; include "math.s7i";const func float: greatCircleDistance (in float: latitude1, in float: longitude1, in float: latitude2, in float: longitude2) is func result var float: distance is 0.0; local const float: EarthRadius is 6372.8; # Average great-elliptic or great-circle radius in kilometers begin distance := 2.0 * EarthRadius * asin(sqrt(sin(0.5 * (latitude2 - latitude1)) ** 2 + cos(latitude1) * cos(latitude2) * sin(0.5 * (longitude2 - longitude1)) ** 2)); end func;const func float: degToRad (in float: degrees) is return degrees * 0.017453292519943295769236907684886127;const proc: main is func begin writeln("Distance in kilometers between BNA and LAX"); writeln(greatCircleDistance(degToRad(36.12), degToRad(-86.67), # Nashville International Airport (BNA) degToRad(33.94), degToRad(-118.4)) # Los Angeles International Airport (LAX) digits 2); end func;2887.26
classEarthPoint(lat,lon){constearth_radius=6371# mean earth radiusconstradian_ratio=Num.pi/180# accessors for radiansmethodlatR{self.lat*radian_ratio}methodlonR{self.lon*radian_ratio}methodhaversine_dist(EarthPointp){vararc=EarthPoint(self.lat-p.lat,self.lon-p.lon,)vara=Math.sum((arc.latR/2).sin**2,(arc.lonR/2).sin**2*self.latR.cos*p.latR.cos)earth_radius*a.sqrt.asin*2}}varBNA=EarthPoint.new(lat:36.12,lon:-86.67)varLAX=EarthPoint.new(lat:33.94,lon:-118.4)sayBNA.haversine_dist(LAX)#=> 2886.444442837983299747157823945746716...
'*** LAT/LONG for Nashville International Airport (BNA)lat1=36.12Lon1=-86.67'*** LAT/LONG for Los Angeles International Airport (LAX)Lat2=33.94Lon2=-118.40'*** Units: K=kilometers M=miles N=nautical milesUnit$ = "K"Result=HAVERSINE(Lat1,Lon1,Lat2,Lon2,Unit$)R$=STR$(Result,"#,###.##")PRINT "The distance between Nashville International Airport and Los Angeles International Airport in kilometers is: "&R$STOPDEF HAVERSINE(Lat1,Lon1,Lat2,Lon2,Unit$)'---------------------------------------------------------------'*** Haversine Formula - Calculate distances by LAT/LONG''*** Pass to it the LAT/LONG of the two locations, and then unit of measure'*** Usage: X=HAVERSINE(Lat1,Lon1,Lat2,Lon2,Unit$) PI=3.14159265358979323846 Radius=6372.8 Lat1=(Lat1*PI/180) Lon1=(Lon1*PI/180) Lat2=(Lat2*PI/180) Lon2=(Lon2*PI/180) DLon=Lon1-Lon2 Answer=ACOS(SIN(Lat1)*SIN(Lat2)+COS(Lat1)*COS(Lat2)*COS(DLon))*Radius IF UNIT$="M" THEN Answer=Answer*0.621371192 IF UNIT$="N" THEN Answer=Answer*0.539956803 RETURN AnswerENDDEF
The distance between Nashville International Airport and Los Angeles International Airport in kilometers is: 2,887.26
First, a program to add a distance variable to a dataset, given variables for LAT/LON of two points.
program spheredistversion 15.0syntax varlist(min=4 max=4 numeric), GENerate(namelist max=1) ///[Radius(real 6371) ALTitude(real 0) LABel(string)]confirm new variable `generate'local lat1 : word 1 of `varlist'local lon1 : word 2 of `varlist'local lat2 : word 3 of `varlist'local lon2 : word 4 of `varlist'local r=2*(`radius'+`altitude'/1000)local k=_pi/180gen `generate'=`r'*asin(sqrt(sin((`lat2'-`lat1')*`k'/2)^2+ ///cos(`lat1'*`k')*cos(`lat2'*`k')*sin((`lon2'-`lon1')*`k'/2)^2))if `"`label'"' != "" {label variable `generate' `"`label'"'}endIllustration with a sample dataset.
import delimited airports.csv, clearformat %9.4f l*list +----------------------------------------------------------------------------------------------------+ | iata airport city country lat lon | |----------------------------------------------------------------------------------------------------| 1. | AMS Amsterdam Airport Schiphol Amsterdam Netherlands 52.3086 4.7639 | 2. | BNA Nashville International Airport Nashville United States 36.1245 -86.6782 | 3. | CDG Charles de Gaulle International Airport Paris France 49.0128 2.5500 | 4. | CGN Cologne Bonn Airport Cologne Germany 50.8659 7.1427 | 5. | LAX Los Angeles International Airport Los Angeles United States 33.9425 -118.4080 | |----------------------------------------------------------------------------------------------------| 6. | MEM Memphis International Airport Memphis United States 35.0424 -89.9767 | +----------------------------------------------------------------------------------------------------+
MEM/CGN joins two Fedex Express hubs. The line AMS/LAX is operated by KLM Royal Dutch Airlines.We will compute the distance between each pair of airports, both at sea level and at typical cruising flight level (35000 ft).
Bear in mind that the actual route of an airliner is usually not a piece of great circle, so this will only give an idea. For instance, according toFlightAware, the route of a Fedex flight from Memphis to Paris is 7852 km long, at FL300 altitude (9150 m). The program given here would yield 7328.33 km instead.
keep iata lat lonrename (iata lat lon) =2gen k=0tempfile tmpsave "`tmp'"rename *2 *1joinby k using `tmp'drop if iata1>=iata2drop klist +-----------------------------------------------------------+ | iata1 lat1 lon1 iata2 lat2 lon2 | |-----------------------------------------------------------| 1. | AMS 52.3086 4.7639 BNA 36.1245 -86.6782 | 2. | AMS 52.3086 4.7639 CGN 50.8659 7.1427 | 3. | AMS 52.3086 4.7639 LAX 33.9425 -118.4080 | 4. | AMS 52.3086 4.7639 CDG 49.0128 2.5500 | 5. | AMS 52.3086 4.7639 MEM 35.0424 -89.9767 | |-----------------------------------------------------------| 6. | BNA 36.1245 -86.6782 CGN 50.8659 7.1427 | 7. | BNA 36.1245 -86.6782 CDG 49.0128 2.5500 | 8. | BNA 36.1245 -86.6782 LAX 33.9425 -118.4080 | 9. | BNA 36.1245 -86.6782 MEM 35.0424 -89.9767 | 10. | CDG 49.0128 2.5500 LAX 33.9425 -118.4080 | |-----------------------------------------------------------| 11. | CDG 49.0128 2.5500 MEM 35.0424 -89.9767 | 12. | CDG 49.0128 2.5500 CGN 50.8659 7.1427 | 13. | CGN 50.8659 7.1427 LAX 33.9425 -118.4080 | 14. | CGN 50.8659 7.1427 MEM 35.0424 -89.9767 | 15. | LAX 33.9425 -118.4080 MEM 35.0424 -89.9767 | +-----------------------------------------------------------+
Now compute the distances and print the result.
spheredist lat1 lon1 lat2 lon2, gen(dist) lab(Distance at sea level)spheredist lat1 lon1 lat2 lon2, gen(fl350) alt(10680) lab(Distance at FL350 altitude)format %9.2f dist fl350list iata* dist fl350 +-----------------------------------+ | iata1 iata2 dist fl350 | |-----------------------------------| 1. | AMS CGN 229.64 230.03 | 2. | AMS CDG 398.27 398.94 | 3. | AMS MEM 7295.19 7307.56 | 4. | AMS BNA 7004.61 7016.48 | 5. | AMS LAX 8955.95 8971.13 | |-----------------------------------| 6. | BNA LAX 2886.32 2891.21 | 7. | BNA CGN 7222.75 7234.99 | 8. | BNA CDG 7018.39 7030.29 | 9. | BNA MEM 321.62 322.16 | 10. | CDG LAX 9102.51 9117.94 | |-----------------------------------| 11. | CDG CGN 387.82 388.48 | 12. | CDG MEM 7317.82 7330.23 | 13. | CGN LAX 9185.47 9201.04 | 14. | CGN MEM 7514.96 7527.70 | 15. | LAX MEM 2599.71 2604.12 | +-----------------------------------+
Notice that the distance from Nashville to Los Angeles is given as 2886.32 km, which is slightly different from the task description. The coordinates come fromOpenFlights and are supposably more accurate. Using the data in the task description, one gets 2886.44 as expected.
import Foundationfunc haversine(lat1:Double, lon1:Double, lat2:Double, lon2:Double) -> Double { let lat1rad = lat1 * Double.pi/180 let lon1rad = lon1 * Double.pi/180 let lat2rad = lat2 * Double.pi/180 let lon2rad = lon2 * Double.pi/180 let dLat = lat2rad - lat1rad let dLon = lon2rad - lon1rad let a = sin(dLat/2) * sin(dLat/2) + sin(dLon/2) * sin(dLon/2) * cos(lat1rad) * cos(lat2rad) let c = 2 * asin(sqrt(a)) let R = 6372.8 return R * c}print(haversine(lat1:36.12, lon1:-86.67, lat2:33.94, lon2:-118.40))2887.25995060711
lat1 : 36.12lon1 : -86.67lat2 : 33.94lon2 : -118.4dx : 0.dy : 0.dz : 0.kms : 0. {degtorad(lon2 - lon1)} lon1 {degtorad lat1} lat1 {degtorad lat2} lat2 {sin lat1 - sin lat2} dz {cos lon1 * cos lat1 - cos lat2} dx {sin lon1 * cos lat1} dy {arcsin(sqrt(dx^2 + dy^2 + dz^2)/2) * 12745.6} kms "'Haversine distance: ' kms ' kms'" []Haversine distance: 2887.259951 kms
option angle radians ' the defaultsub haversine(lat1, lon1, lat2, lon2)dim EarthRadiusKm = 6372.8 ' Earth radius in kilometersdim latRad1 = RAD(lat1)dim latRad2 = RAD(lat2)dim lonRad1 = RAD(lon1)dim lonRad2 = RAD(lon2)dim _diffLa = latRad2 - latRad1dim _doffLo = lonRad2 - lonRad1dim sinLaSqrd = sin(_diffLa / 2) ^ 2dim sinLoSqrd = sin(_doffLo / 2) ^ 2dim computation = asin(sqrt(sinLaSqrd + cos(latRad1) * cos(latRad2) * sinLoSqrd))return 2 * EarthRadiusKm * computationend subprint using "Nashville International Airport to Los Angeles International Airport ####.########### km", haversine(36.12, -86.67, 33.94, -118.40)print using "Perth, WA Australia to Baja California, Mexico #####.########### km", haversine(-31.95, 115.86, 31.95, -115.86)
Nashville International Airport to Los Angeles International Airport 2887.25995060712 kmPerth, WA Australia to Baja California, Mexico 15188.70229560390 km
package require Tcl 8.5proc haversineFormula {lat1 lon1 lat2 lon2} { set rads [expr atan2(0,-1)/180] set R 6372.8 ;# In kilometers set dLat [expr {($lat2-$lat1) * $rads}] set dLon [expr {($lon2-$lon1) * $rads}] set lat1 [expr {$lat1 * $rads}] set lat2 [expr {$lat2 * $rads}] set a [expr {sin($dLat/2)**2 + sin($dLon/2)**2*cos($lat1)*cos($lat2)}] set c [expr {2*asin(sqrt($a))}] return [expr {$R * $c}]}# Don't bother with too much inappropriate accuracy!puts [format "distance=%.1f km" [haversineFormula 36.12 -86.67 33.94 -118.40]]distance=2887.3 km
{{trans|BASIC}}FUNCTION HAVERSINE!---------------------------------------------------------------!*** Haversine Formula - Calculate distances by LAT/LONG!!*** LAT/LON of the two locations and Unit of measure are GLOBAL!*** as they are defined in the main logic of the program, so they!*** available for use in the Function.!*** Usage: X=HAVERSINE Radius=6378.137 Lat1=(Lat1*MATH.PI/180) Lon1=(Lon1*MATH.PI/180) Lat2=(Lat2*MATH.PI/180) Lon2=(Lon2*MATH.PI/180) DLon=Lon1-Lon2 ANSWER=ACOS(SIN(Lat1)*SIN(Lat2)+COS(Lat1)*COS(Lat2)*COS(DLon))*Radius DISTANCE="kilometers" SELECT CASE UNIT CASE "M" HAVERSINE=ANSWER*0.621371192 Distance="miles" CASE "N" HAVERSINE=ANSWER*0.539956803 Distance="nautical miles" END SELECT END FUNCTION
The following is the main code that invokes the function. It takes your location and determines how far away you are from Tampa, Florida. You can change UNIT to either M=Miles, N=Nautical Miles, or K (or leave blank) as default is in Kilometers:
!*** In techBASIC, all variables defined in the main program act as GLOBAL!*** variables and are available to all SUBROUTINES and FUNCTIONS. So in the!*** HAVERSINE Function being used, no paramaters need to be passed to it, so!*** it acts as a variable when I use it as Result=HAVERSINE. The way that!*** the Function is setup, it returns its value back as HAVERSINE.BASE 1!*** Get the GPS LAT/LONG of current locationlocation = sensors.location(30)Lat1=location(1) Lon1=location(2) !*** LAT/LONG For Tampa, FLLat2=27.9506Lon2=-82.4572!*** Units: K=kilometers M=miles N=nautical milesDIM UNIT AS STRING DIM Distance AS STRINGDIM Result AS SINGLEUNIT = "M"!*** Calculate distance using Haversine FunctionResult=HAVERSINEPRINT "The distance from your current location to Tampa, FL in ";Distance;" is: ";PRINT USING "#,###.##";Result;"."STOP
OUTPUT: *** NOTE: When I run this, I am in my house in Venice, Florida, and that distance is correct (as the crow flies). ***
The distance from your current location to Tampa, FL in miles is: 57.94
# syntax: call SP_HAVERSINE(36.12,33.94,-86.67,-118.40,x);CREATE PROCEDURE SP_HAVERSINE(IN lat1 FLOAT,IN lat2 FLOAT,IN lon1 FLOAT,IN lon2 FLOAT,OUT distance FLOAT) BEGIN DECLARE dLat FLOAT; DECLARE dLon FLOAT; DECLARE c FLOAT; DECLARE a FLOAT; DECLARE km FLOAT; SET dLat = RADIANS(lat2-lat1); SET dLon = RADIANS(lon2-lon1); SET a = SIN(dLat / 2) * SIN(dLat / 2) + SIN(dLon / 2) * SIN(dLon / 2) * COS(RADIANS(lat1)) * COS(RADIANS(lat2)); SET c = 2 * ASIN(SQRT(a)); SET km = 6372.8 * c; select km into distance;END;
distance: 2887.2599 km
CREATE FUNCTION [dbo].[Haversine](@Lat1 AS DECIMAL(9,7), @Lon1 AS DECIMAL(10,7), @Lat2 AS DECIMAL(9,7), @Lon2 AS DECIMAL(10,7))RETURNS DECIMAL(12,7)ASBEGINDECLARE @RDECIMAL(11,7);DECLARE @dLatDECIMAL(9,7);DECLARE @dLonDECIMAL(10,7);DECLARE @aDECIMAL(10,7);DECLARE @cDECIMAL(10,7);SET @R= 6372.8;SET @dLat= RADIANS(@Lat2 - @Lat1);SET @dLon= RADIANS(@Lon2 - @Lon1);SET @Lat1= RADIANS(@Lat1);SET @Lat2= RADIANS(@Lat2);SET @a= SIN(@dLat / 2) * SIN(@dLat / 2) + SIN(@dLon / 2) * SIN(@dLon / 2) * COS(@Lat1) * COS(@Lat2);SET @c= 2 * ASIN(SQRT(@a));RETURN @R * @c;ENDGOSELECT dbo.Haversine(36.12,-86.67,33.94,-118.4)
2887.2594934
let radians = function (degree: number) { // degrees to radians let rad: number = degree * Math.PI / 180; return rad;}export const haversine = (lat1: number, lon1: number, lat2: number, lon2: number) => { // var dlat: number, dlon: number, a: number, c: number, R: number; let dlat, dlon, a, c, R: number; R = 6372.8; // km dlat = radians(lat2 - lat1); dlon = radians(lon2 - lon1); lat1 = radians(lat1); lat2 = radians(lat2); a = Math.sin(dlat / 2) * Math.sin(dlat / 2) + Math.sin(dlon / 2) * Math.sin(dlon / 2) * Math.cos(lat1) * Math.cos(lat2) c = 2 * Math.asin(Math.sqrt(a)); return R * c;}console.log("Distance:" + haversine(36.12, -86.67, 33.94, -118.40));Distance: 2887.2599506071106
10 Point 7 'Sets decimal display to 32 places (0+.1^56) 20 Rf=#pi/180 'Degree -> Radian Conversion 100 ?Using(,7),.DxH(36+7.2/60,-(86+40.2/60),33+56.4/60,-(118+24/60));" km" 999 End 1000 '*** Haversine Distance Function *** 1010 .DxH(Lat_s,Long_s,Lat_f,Long_f) 1020 L_s=Lat_s*rf:L_f=Lat_f*rf:LD=L_f-L_s:MD=(Long_f-Long_s)*rf 1030 Return(12745.6*asin( (sin(.5*LD)^2+cos(L_s)*cos(L_f)*sin(.5*MD)^2)^.5)) '' '' Run 2887.2599506 km OK
Const MER = 6371 '-- mean earth radius(km)Public DEG_TO_RAD As Double Function haversine(lat1 As Double, long1 As Double, lat2 As Double, long2 As Double) As Double lat1 = lat1 * DEG_TO_RAD lat2 = lat2 * DEG_TO_RAD long1 = long1 * DEG_TO_RAD long2 = long2 * DEG_TO_RAD haversine = MER * WorksheetFunction.Acos(Sin(lat1) * Sin(lat2) + Cos(lat1) * Cos(lat2) * Cos(long2 - long1))End Function Public Sub main() DEG_TO_RAD = WorksheetFunction.Pi / 180 d = haversine(36.12, -86.67, 33.94, -118.4) Debug.Print "Distance is "; Format(d, "#.######"); " km ("; Format(d / 1.609344, "#.######"); " miles)."End SubDistance is 2886,444443 km (1793,553425 miles).
If you read the fine print in the Wikipedia article, you will find that the Haversine method of finding distances may have an error of up to 0.5%. This could lead one to believe that discussion about whether to use 6371.0 km or 6372.8 km for an approximation of the Earth's radius is moot.
Imports System.MathModule Module1 Const deg2rad As Double = PI / 180 Structure AP_Loc Public IATA_Code As String, Lat As Double, Lon As Double Public Sub New(ByVal iata_code As String, ByVal lat As Double, ByVal lon As Double) Me.IATA_Code = iata_code : Me.Lat = lat * deg2rad : Me.Lon = lon * deg2rad End Sub Public Overrides Function ToString() As String Return String.Format("{0}: ({1}, {2})", IATA_Code, Lat / deg2rad, Lon / deg2rad) End Function End Structure Function Sin2(ByVal x As Double) As Double Return Pow(Sin(x / 2), 2) End Function Function calculate(ByVal one As AP_Loc, ByVal two As AP_Loc) As Double Dim R As Double = 6371, ' In kilometers, (as recommended by the International Union of Geodesy and Geophysics) a As Double = Sin2(two.Lat - one.Lat) + Sin2(two.Lon - one.Lon) * Cos(one.Lat) * Cos(two.Lat) Return R * 2 * Asin(Sqrt(a)) End Function Sub ShowOne(pntA As AP_Loc, pntB as AP_Loc) Dim adst As Double = calculate(pntA, pntB), sfx As String = "km" If adst < 1000 Then adst *= 1000 : sfx = "m" Console.WriteLine("The approximate distance between airports {0} and {1} is {2:n2} {3}.", pntA, pntB, adst, sfx) Console.WriteLine("The uncertainty is under 0.5%, or {0:n1} {1}." & vbLf, adst / 200, sfx) End Sub' Airport coordinate data excerpted from the data base at http://www.partow.net/miscellaneous/airportdatabase/' The four additional airports are the furthest and closest pairs, according to the "Fun Facts..." section.' KBNA, BNA, NASHVILLE INTERNATIONAL, NASHVILLE, USA, 036, 007, 028, N, 086, 040, 041, W, 00183, 36.124, -86.678' KLAX, LAX, LOS ANGELES INTERNATIONAL, LOS ANGELES, USA, 033, 056, 033, N, 118, 024, 029, W, 00039, 33.942, -118.408' SKNV, NVA, BENITO SALAS, NEIVA, COLOMBIA, 002, 057, 000, N, 075, 017, 038, W, 00439, 2.950, -75.294' WIPP, PLM, SULTAN MAHMUD BADARUDDIN II, PALEMBANG, INDONESIA, 002, 053, 052, S, 104, 042, 004, E, 00012, -2.898, 104.701 ' LOWL, LNZ, HORSCHING INTERNATIONAL AIRPORT (AUS - AFB), LINZ, AUSTRIA, 048, 014, 000, N, 014, 011, 000, E, 00096, 48.233, 14.183' LOXL, N/A, LINZ, LINZ, AUSTRIA, 048, 013, 059, N, 014, 011, 015, E, 00299, 48.233, 14.188 Sub Main() ShowOne(New AP_Loc("BNA", 36.124, -86.678), New AP_Loc("LAX", 33.942, -118.408)) ShowOne(New AP_Loc("NVA", 2.95, -75.294), New AP_Loc("PLM", -2.898, 104.701)) ShowOne(New AP_Loc("LNZ", 48.233, 14.183), New AP_Loc("N/A", 48.233, 14.188)) End SubEnd ModuleThe approximate distance between airports BNA: (36.124, -86.678) and LAX: (33.942, -118.408) is 2,886.36 km.The uncertainty is under 0.5%, or 14.4 km.The approximate distance between airports NVA: (2.95, -75.294) and PLM: (-2.898, 104.701) is 20,009.28 km.The uncertainty is under 0.5%, or 100.0 km.The approximate distance between airports LNZ: (48.233, 14.183) and N/A: (48.233, 14.188) is 370.34 m.The uncertainty is under 0.5%, or 1.9 m.
Looking at the altitude difference between the last two airports, (299 - 96 = 203), the reported distance of 370 meters ought to be around 422 meters if you actually went there and saw it for yourself.
import mathfn haversine(h f64) f64 { return .5 * (1 - math.cos(h))} struct Pos { lat f64 // latitude, radians long f64 // longitude, radians} fn deg_pos(lat f64, lon f64) Pos { return Pos{lat * math.pi / 180, lon * math.pi / 180}} const r_earth = 6372.8 // km fn hs_dist(p1 Pos, p2 Pos) f64 { return 2 * r_earth * math.asin(math.sqrt(haversine(p2.lat-p1.lat)+ math.cos(p1.lat)*math.cos(p2.lat)*haversine(p2.long-p1.long)))} fn main() { println(hs_dist(deg_pos(36.12, -86.67), deg_pos(33.94, -118.40)))}2887.2599506071
var R = 6372.8 // Earth's approximate radius in kilometers./* Class containing trig methods which work with degrees rather than radians. */class D { static deg2Rad(deg) { (deg*Num.pi/180 + 2*Num.pi) % (2*Num.pi) } static sin(d) { deg2Rad(d).sin } static cos(d) { deg2Rad(d).cos }}var haversine = Fn.new { |lat1, lon1, lat2, lon2| var dlat = lat2 - lat1 var dlon = lon2 - lon1 return 2 * R * (D.sin(dlat/2).pow(2) + D.cos(lat1) * D.cos(lat2) * D.sin(dlon/2).pow(2)).sqrt.asin}System.print(haversine.call(36.12, -86.67, 33.94, -118.4))2887.2599506071
Assemble with tasm /m /l; tlink /t
0000 .model tiny0000 .code .486 org 100h ;.com files start here0100 9B DB E3 start: finit ;initialize floating-point unit (FPU) ;Great circle distance = ; 2.0*Radius * ASin( sqrt( Haversine(Lat2-Lat1) + ; Haversine(Lon2-Lon1)*Cos(Lat1)*Cos(Lat2) ) )0103 D9 06 0191r fld Lat2 ;push real onto FPU stack0107 D8 26 018Dr fsub Lat1 ;subtract real from top of stack (st(0) = st)010B E8 0070 call Haversine ;(1.0-cos(st)) / 2.0010E D9 06 0199r fld Lon2 ;repeat for longitudes0112 D8 26 0195r fsub Lon10116 E8 0065 call Haversine ;st(1)=Lats; st=Lons0119 D9 06 018Dr fld Lat1011D D9 FF fcos ;replace st with its cosine011F D9 06 0191r fld Lat20123 D9 FF fcos ;st=cos(Lat2); st(1)=cos(Lat1); st(2)=Lats; st(3)=Lons0125 DE C9 fmul ;st=cos(Lat2)*cos(Lat1); st(1)=Lats; st(2)=Lons0127 DE C9 fmul ;st=cos(Lat2)*cos(Lat1)*Lats; st(1)=Lons0129 DE C1 fadd ;st=cos(Lat2)*cos(Lat1)*Lats + Lons012B D9 FA fsqrt ;replace st with its square root ;asin(x) = atan(x/sqrt(1-x^2))012D D9 C0 fld st ;duplicate tos012F D8 C8 fmul st, st ;x^20131 D9 E8 fld1 ;get 1.00133 DE E1 fsubr ;1 - x^20135 D9 FA fsqrt ;sqrt(1-x^2)0137 D9 F3 fpatan ;take atan(st(1)/st)0139 D8 0E 019Dr fmul Radius2 ;*2.0*Radius ;Display value in FPU's top of stack (st) =0004 before equ 4 ;places before =0002 after equ 2 ; and after decimal point =0001 scaler = 1 ;"=" allows scaler to be redefined, unlike equ rept after ;repeat block "after" times scaler = scaler*10 endm ;scaler now = 10^after013D 66| 6A 64 push dword ptr scaler;use stack for convenient memory location0140 67| DA 0C 24 fimul dword ptr [esp] ;st:= st*scaler0144 67| DB 1C 24 fistp dword ptr [esp] ;round st to nearest integer0148 66| 58 pop eax ; and put it into eax014A 66| BB 0000000A mov ebx, 10 ;set up for idiv instruction0150 B9 0006 mov cx, before+after;set up loop counter0153 66| 99 ro10: cdq ;convert double to quad; i.e: edx:= 00155 66| F7 FB idiv ebx ;eax:= edx:eax/ebx; remainder in edx0158 52 push dx ;save least significant digit on stack0159 E2 F8 loop ro10 ;cx--; loop back if not zero015B B1 06 mov cl, before+after;(ch=0)015D B3 00 mov bl, 0 ;used to suppress leading zeros015F 58 ro20: pop ax ;get digit0160 0A D8 or bl, al ;turn off suppression if not a zero0162 80 F9 03 cmp cl, after+1 ;is digit immediately to left of decimal point?0165 75 01 jne ro30 ;skip if not0167 43 inc bx ;turn off leading zero suppression0168 04 30 ro30: add al, '0' ;if leading zero then ' ' else add 0016A 84 DB test bl, bl016C 75 02 jne ro40016E B0 20 mov al, ' '0170 CD 29 ro40: int 29h ;display character in al register0172 80 F9 03 cmp cl, after+1 ;is digit immediately to left of decimal point?0175 75 04 jne ro50 ;skip if not0177 B0 2E mov al, '.' ;display decimal point0179 CD 29 int 29h017B E2 E2 ro50: loop ro20 ;loop until all digits displayed017D C3 ret ;return to OS017E Haversine: ;return (1.0-Cos(Ang)) / 2.0 in st017E D9 FF fcos0180 D9 E8 fld10182 DE E1 fsubr0184 D8 36 0189r fdiv N20188 C3 ret0189 40000000 N2 dd 2.0018D 3F21628D Lat1 dd 0.63041 ;36.12*pi/1800191 3F17A4E8 Lat2 dd 0.59236 ;33.94*pi/1800195 BFC19F80 Lon1 dd -1.51268 ;-86.67*pi/1800199 C004410B Lon2 dd -2.06647 ;-118.40*pi/180019D 46472666 Radius2 dd 12745.6 ;6372.8 average radius of Earth (km) times 2 ;(TASM isn't smart enough to do floating point constant calculations) end start
2887.25
include c:\cxpl\codes; \intrinsic 'code' declarationsfunc real Haversine(Ang);real Ang;return (1.0-Cos(Ang)) / 2.0;func real Dist(Lat1, Lat2, Lon1, Lon2); \Great circle distancereal Lat1, Lat2, Lon1, Lon2;def R = 6372.8; \average radius of Earth (km)return 2.0*R * ASin( sqrt( Haversine(Lat2-Lat1) + Cos(Lat1)*Cos(Lat2)*Haversine(Lon2-Lon1) ));def D2R = 3.141592654/180.0; \degrees to radiansRlOut(0, Dist(36.12*D2R, 33.94*D2R, -86.67*D2R, -118.40*D2R ));
2887.25995
declare namespace xsd = "http://www.w3.org/2001/XMLSchema";declare namespace math = "http://www.w3.org/2005/xpath-functions/math";declare function local:haversine($lat1 as xsd:float, $lon1 as xsd:float, $lat2 as xsd:float, $lon2 as xsd:float) as xsd:float{ let $dlat := ($lat2 - $lat1) * math:pi() div 180 let $dlon := ($lon2 - $lon1) * math:pi() div 180 let $rlat1 := $lat1 * math:pi() div 180 let $rlat2 := $lat2 * math:pi() div 180 let $a := math:sin($dlat div 2) * math:sin($dlat div 2) + math:sin($dlon div 2) * math:sin($dlon div 2) * math:cos($rlat1) * math:cos($rlat2) let $c := 2 * math:atan2(math:sqrt($a), math:sqrt(1-$a)) return xsd:float($c * 6371.0)};local:haversine(36.12, -86.67, 33.94, -118.4)2886.444
When a Zigstruct type can be inferred then anonymous structs .{} can be used for initialisation. This can be seen on the lines where the constantsbna andlax are instantiated.
A Zigstruct can have methods, the same as anenum and or aunion.They are only namespaced functions that can be called with dot syntax.
const std = @import("std");const math = std.math; // Save some typing, reduce clutter. Otherwise math.sin() would be std.math.sin() etc.pub fn main() !void { // Coordinates are found here: // http://www.airport-data.com/airport/BNA/ // http://www.airport-data.com/airport/LAX/ const bna = LatLong{ .lat = .{ .d = 36, .m = 7, .s = 28.10 }, .long = .{ .d = 86, .m = 40, .s = 41.50 }, }; const lax = LatLong{ .lat = .{ .d = 33, .m = 56, .s = 32.98 }, .long = .{ .d = 118, .m = 24, .s = 29.05 }, }; const distance = calcGreatCircleDistance(bna, lax); std.debug.print("Output: {d:.6} km\n", .{distance}); // Output: 2886.326609 km}const LatLong = struct { lat: DMS, long: DMS };/// degrees, minutes, decimal secondsconst DMS = struct { d: f64, m: f64, s: f64, fn toRadians(self: DMS) f64 { return (self.d + self.m / 60 + self.s / 3600) * math.pi / 180; }};// Volumetric mean radius is 6371 km, see http://nssdc.gsfc.nasa.gov/planetary/factsheet/earthfact.html// The diameter is thus 12742 kmfn calcGreatCircleDistance(lat_long1: LatLong, lat_long2: LatLong) f64 { const lat1 = lat_long1.lat.toRadians(); const lat2 = lat_long2.lat.toRadians(); const long1 = lat_long1.long.toRadians(); const long2 = lat_long2.long.toRadians(); const a = math.sin(0.5 * (lat2 - lat1)); const b = math.sin(0.5 * (long2 - long1)); return 12742 * math.asin(math.sqrt(a * a + math.cos(lat1) * math.cos(lat2) * b * b));}haversine(36.12, -86.67, 33.94, -118.40).println(); fcn haversine(Lat1, Long1, Lat2, Long2){ const R = 6372.8; // In kilometers; Diff_Lat := (Lat2 - Lat1) .toRad(); Diff_Long := (Long2 - Long1).toRad(); NLat := Lat1.toRad(); NLong := Lat2.toRad(); A := (Diff_Lat/2) .sin().pow(2) + (Diff_Long/2).sin().pow(2) * NLat.cos() * NLong.cos(); C := 2.0 * A.sqrt().asin(); R*C;}2887.26
10 LET diam=2*6372.820 LET Lg1m2=FN r((-86.67)-(-118.4))30 LET Lt1=FN r(36.12)40 LET Lt2=FN r(33.94)50 LET dz=SIN (Lt1)-SIN (Lt2)60 LET dx=COS (Lg1m2)*COS (Lt1)-COS (Lt2)70 LET dy=SIN (Lg1m2)*COS (Lt1)80 LET hDist=ASN ((dx*dx+dy*dy+dz*dz)^0.5/2)*diam90 PRINT "Haversine distance: ";hDist;" km."100 STOP 1000 DEF FN r(a)=a*0.017453293: REM convert degree to radians