Home > IT Answers > AS/400 > Haversine formula in RPG or SQL
Brenger
170 pts.
 Haversine formula in RPG or SQL
AS/400 development, Great-circle distance, Haversine formula, RPG, SQL
I need to put the haversine formula in an sql or rpg function on the as400 to calc the distance between to zip codes. I do not have much experience with sql. can anyone give me some help?

Software/Hardware used:
ASKED: October 31, 2008  10:01 PM
UPDATED: November 6, 2008  5:42 PM
RELATED QUESTIONS
Answer Wiki:
Here's the formula in java ---------------------------------------------------------------------- var R = 6371; // km var dLat = (lat2-lat1).toRad(); var dLon = (lon2-lon1).toRad(); var a = Math.sin(dLat/2) * Math.sin(dLat/2) + Math.cos(lat1.toRad()) * Math.cos(lat2.toRad()) * Math.sin(dLon/2) * Math.sin(dLon/2); var c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a)); var d = R * c; ------------------------------------------------------------------- Not RPG or SQL (I think) need to use C math stuff in RPG This is one for a C math - RPG dud/dudette well, based on the JAVA code, then this should work <pre> <b>SQL</b> h option(*nodebugio:*srcstmt) dftactgrp(*no) actgrp(*caller) indent('|') d test s 8f d Haversine pr 8f d lat1 8f const d lon1 8f const d lat2 8f const d lon2 8f const c eval test = haversine(1:1:2:2) c dump(a) c eval *inlr = *on p Haversine b d Haversine pi 8f d lat1 8f const d lon1 8f const d lat2 8f const d lon2 8f const d dLat s 8f d dLon s 8f d a s 8f . d c s 8f c/exec sql c+ set :dLat = radians(:lat2 - :lat1) c/end-exec c/exec sql c+ set :dLon = radians(:lon2 - :lon1) c/end-exec c/exec sql c+ set :a = sin(:dLat/2) * sin(:dLat/2) + c+ cos(radians(:lat1)) * cos(radians(:lat2)) * c+ sin(:dLon/2) * sin(:dLon/2) c/end-exec c/exec sql c+ set :c = 2 * atan2(sqrt(:a), sqrt(1 - :a)) c/end-exec c return 6371 * c p Haversine e <b>C functions</b> h option(*nodebugio:*srcstmt) dftactgrp(*no) actgrp(*caller) indent('|') h bnddir('QC2LE') expropts(*maxdigits) fltdiv d test s 8f d sin pr 8f extproc('sin') d parm1 8f const d cos pr 8f extproc('cos') d parm1 8f const d atan2 pr 8f extproc('atan2') d parm1 8f const d parm2 8f const d sqrt pr 8f extproc('sqrt') d parm1 8f const d Haversine pr 8f d lat1 8f const d lon1 8f const d lat2 8f const d lon2 8f const c eval test = haversine(1:1:2:2) c dump(a) c eval *inlr = *on p Haversine b d Haversine pi 8f d lat1 8f const d lon1 8f const d lat2 8f const d lon2 8f const d dlat s 8f d dlon s 8f d a s 8f d c s 8f d pi c const(3.1415926) /free dlat = (lat2 - lat1) * pi / 180; dlon = (lon2 - lon1) * pi / 180; a = sin(dlat / 2)**2 + cos(lat1 * pi / 180) * cos(lat2 * pi / 180) * sin(dlon / 2)**2; c = 2 * atan2(sqrt(a):sqrt(1 - a)); return 6371 * c; /end-free p Haversine e </pre> Kevin C. Ketzler - <a href="http://www.aresgrp.com">Affiliated</a>
Last Wiki Answer Submitted:  November 6, 2008  3:51 pm  by  BigKat   7,185 pts.
All Answer Wiki Contributors:  BigKat   7,185 pts. , philpl1jb   44,180 pts.
To see all answers submitted to the Answer Wiki: View Answer History.


RSS - Discussions Help
Discuss This Question:
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _

 

There was a similar discussion that may be relevant on devshed.com:

C program for Haversine formula

Leah Rosin
 490 pts.
 

Cool but that doesn’t do the RPG bit.

philpl1jb
 44,180 pts.
 

C math in RPG
http://www.rpgiv.com/rpgnews/Feb99a/highmath.html
C library
http://publib.boulder.ibm.com/iseries/v5r2/ic2924/books/c4156070.pdf

Anyone got a sample RPGLE -C math
need to radians, sin, cos, atan, sqrt

philpl1jb
 44,180 pts.
 

to convert from degrees to radians I used:

degrees * pi/180 = radians

where pi = 3.1415926 (approx.) since I wasn’t sure if there was a c function that did the conversion.

BigKat
 7,185 pts.
 

also return 6371 * c returns kilometers
use 3956 for miles

BigKat
 7,185 pts.
 

Original request was Zip to distiance. Anyone know a table of zip to Lat/Long.
Computations are in km, you probably want them in miles

philpl1jb
 44,180 pts.
 

got the zip table file from the census…thanks. you guys are great!

Brenger
 170 pts.
 

OK guys, one more question.

C+ set :c = 2 * atan2(sqrt(:a), sqrt(1 – :a))

C+ set :c = 2 * asin(min(sqrt(:a)))

It doesn’t seem to calc correctly with stmt one. another version i’ve seen out there is stmt 2. but RPG doesn’t like the min part. any ideas on that?

Brenger
 170 pts.
 

The min function returns the smaller of two values, looks like the formula you found
should have been min(1,sqrt(a)) ??

heres a site discussing it
http://www.cs.nyu.edu/visual/home/proj/tiger/gisfaq.html
but the original source seems to be gone.

dlon = lon2 – lon1
dlat = lat2 – lat1
a = (sin(dlat/2))^2 + cos(lat1) * cos(lat2) * (sin(dlon/2))^2
c = 2 * arcsin(min(1,sqrt(a)))
d = R * c

where R is the radius of the Earth.

philpl1jb
 44,180 pts.
 

if the JAVA program is not correct, what I posted will not be correct as I just converted what was there

copy the sin prototype and change “sin” to “asin” both times, 

d asin            pr             8f   extproc('asin')      
d  parm1                         8f   const

then 

if a < 1;
c = 2 * asin(sqrt(a));
else;
c = 2 * asin(1);
endif;
BigKat
 7,185 pts.
 

that was it! thanks so much!
c = 2 * arcsin(min(1,sqrt(a)))

Brenger
 170 pts.