## Haversine formula in RPG or SQL

170 pts.
Tags:
AS/400 development
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?

Thanks. We'll let you know when a new response is added.

Here’s the formula in java:
———————————————————————-

```var R = 6371; // km
var a = Math.sin(dLat/2) * Math.sin(dLat/2) +
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

```SQL
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+ 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

C functions
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```

Kevin C. Ketzler –

## Discuss This Question: 11 Replies

Thanks. We'll let you know when a new response is added.
• There was a similar discussion that may be relevant on devshed.com: C program for Haversine formula
report
• Cool but that doesn't do the RPG bit.
report
• 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
report
• 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.
report
• also return 6371 * c returns kilometers use 3956 for miles
report
• 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
report
• got the zip table file from the census...thanks. you guys are great!
report
• 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?
report
• 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.
report
• 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;
```