Re: longitude / latitude distance calculations

[Scott Klement <rpg400-l@xxxxxxxxxxxxxxxx>]

Hi Karen,

> Working on a project and need to know the distance in miles between two
> locations.  The information that I have in the file to use to do this is the
> longitude and latitude.  I have seen the code in PC jargon, but have not
> been able to find anything about how to accomplish this in RPG.
[SNIP]

I took the visual basic code you posted, and converted (literally, 
line-by-line) to RPG code.  I haven't tested this code, so you'll want to 
compare the results you get from your VB code to this to make sure I 
didn't goof something up.

Though, one thing didn't make sense to me.  In the VB code you've got
"Select Case ucase(unit)" but, the values you're comparing that you 
compare it to are lowercase?  Maybe I'm misremembering, but doesn't 
ucase() convert to uppercase? (I haven't done any serious VB code in about 
8 years)  So I did that part slightly differently...  If you find that the 
VB code always returns it's result in miles, you know why :)

Anyway, here's the code.  You should be able to compile it with
CRTBNDRPG MYPGM SRCFILE(mylib/QRPGLESRC), and then run it...

     H DFTACTGRP(*NO) BNDDIR('QC2LE')

     D distance        PR             8F
     D   lat1                         8F   value
     D   lon1                         8F   value
     D   lat2                         8F   value
     D   lon2                         8F   value
     D   unit                         1A   const
     D deg2rad         PR             8F
     D   deg                          8F   value
     D rad2deg         PR             8F
     D   rad                          8F   value

     D pi              s              8F   inz(3.14159265358979323846)
     D result          s             20P10
     D wait            s              1A

      /free
          result = distance(32.9697: -96.80322: 29.46786: -98.53506: 'm');
          dsply (%char(result) + ' miles');

          result = distance(32.9697: -96.80322: 29.46786: -98.53506: 'k');
          dsply (%char(result) + ' kilometers');

          result = distance(32.9697: -96.80322: 29.46786: -98.53506: 'n');
          dsply (%char(result) + ' nautical miles');

          dsply ' ' ' ' wait;
          *inlr = *on;
      /end-free

      *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      * distance(): Calculate distance between two points
      *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     P distance        B                   export
     D distance        PI             8F
     D   lat1                         8F   value
     D   lon1                         8F   value
     D   lat2                         8F   value
     D   lon2                         8F   value
     D   unit                         1A   const

     D sin             PR             8F   extproc('sin')
     D  x                             8F   value
     D cos             PR             8F   extproc('cos')
     D  x                             8F   value
     D acos            PR             8F   extproc('acos')
     D  x                             8F   value

     D distance        s              8F
     D dist            s              8F
     D theta           s              8F

      /free
          theta = lon1 - lon2;
          dist = sin(deg2rad(lat1)) * sin(deg2rad(lat2))
               + cos(deg2rad(lat1)) * cos(deg2rad(lat2))
               * cos(deg2rad(theta));

      /if defined(DEBUGING)
          dsply ('dist = ' + %editflt(dist));
      /endif

          dist = acos(dist);
          dist = rad2deg(dist);

      /if defined(DEBUGING)
          dsply ('dist = ' + %editflt(dist));
      /endif

          distance = dist * 60 * 1.1515;

          select;
          when (unit = 'k' or unit = 'K');
             distance = distance * 1.609344;
          when (unit = 'n' or unit = 'N');
             distance = distance * 0.8684;
          endsl;

          return distance;
      /end-free
     P                 E


      *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      * deg2rad(): Convert degrees to radians
      *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     P deg2rad         B
     D deg2rad         PI             8F
     D   deg                          8F   value
      /free
         return (deg * pi / 180);
      /end-free
     P                 E


      *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      * rad2deg(): Convert radians to degrees
      *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     P rad2deg         B
     D rad2deg         PI             8F
     D   rad                          8F   value
      /free
         return (rad * 180 / pi);
      /end-free
     P                 E