Path: news.mathworks.com!not-for-mail
From: "Rob Comer" <rob.comer.nospam@mathworks.com>
Newsgroups: comp.soft-sys.matlab
Subject: Re: do reckon and scircle give different results?
Date: Sun, 15 May 2011 01:21:05 +0000 (UTC)
Organization: The MathWorks Inc
Lines: 55
Message-ID: <iqn9q1$4p8$1@newscl01ah.mathworks.com>
References: <iqj6e8$498$1@newscl01ah.mathworks.com>
Reply-To: "Rob Comer" <rob.comer.nospam@mathworks.com>
NNTP-Posting-Host: www-00-blr.mathworks.com
Content-Type: text/plain; charset=UTF-8; format=flowed
Content-Transfer-Encoding: 8bit
X-Trace: newscl01ah.mathworks.com 1305422465 4904 172.30.248.45 (15 May 2011 01:21:05 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Sun, 15 May 2011 01:21:05 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 1184234
Xref: news.mathworks.com comp.soft-sys.matlab:726837

"Fatih" wrote in message <iqj6e8$498$1@newscl01ah.mathworks.com>...
> I am trying to create a circular area on earth. I have reference latitude and longitude values to be at the center. This circular area will have many circles in each other at the end. The way it works is from the reference point every 2.5 meters for 0 degree azimuth, I get the latitude and longitude values. I repeat this for 9000 points. Then for 0.25 degree azimuth I repeat this whole process again, until I complete 360 degree azimuth. This is what I am trying to do.
>  My initial approach was using reckon function from Mapping Toolbox:
> geoid = almanac('earth', 'wgs84');
> geoid(1) = geoid(1)*1000; % to make radius use meter unit
> for azimuth = 1:1440
>     for distance = 1:9000
>         [latitude, longitude] = reckon('rh', ref_latitude, ref_longitude, 2.5*distance, 0.25*azimuth, geoid);
>      end
> end
> 
> Later I thought I can do it by using scircle1 instead in a lot shorter time(18 sec comparing to 1 hour):
> 
> for radius = 1:9000
>         [latitude, longitude] = scircle1('rh', ref_latitude, ref_longitude, 2.5*radius, [], geoid, 'degrees', 1440);
> end
> 
> Then I decided to compare the results to see if they produce same results. Here is my question:
> 
> The results are NOT same. I do not understand the reason. Is it what is epected or Should I do something to make them work exactly same?

Side comment: to use meters, just call almanac like this:

  ellipsoid = almanac('earth','wgs84','meters')  % It's really an ellipsoid, not a geoid

You should be able either to use scircle1 (which calls reckon), or to use reckon directly.

In what way do the result differ? Do the circles look different if you plot them, or are the numerical values just in a different order?

I'd actually expect your first approach to produce radial lines, because you iterate over distance from the center point in your inner loop. Is that what's happening?

With scircle1, you can't control the way the azimuth angles are sorted, but when you call reckon you are looping over azimuth with a specific starting value, ending value, and order. Perhaps that's the only difference.

It also looks as if you are unaware that reckon is "vectorized." The reason why your computation with reckon is taking so long is that you are calling it over 10 million times, with scalar inputs. That's hugely inefficient.

Instead of this:

for azimuth = 1:1440
  for distance = 1:9000
     [latitude, longitude] = reckon('rh', ref_latitude, ref_longitude, 2.5*distance, 0.25*azimuth, geoid);
    end
end

try something like this:

az = 0.25*(1:1440);
for distance = 1:9000
  [latitude, longitude] = reckon('rh', ref_latitude, ref_longitude, 2.5*distance, az, ellipsoid); % Avoiding use of 'geoid'
end

You could, in fact, do the whole computation with just one call to reckon (eliminating both loops), but you'd use more memory because you'd have to expand distance and az into 1440-by-9000 matrices.

Rob Comer
Mapping Toolbox Development
MathWorks