| circles on the sky [message #65879] |
Fri, 27 March 2009 10:53  |
Christopher Thom
Messages: 66
Registered: October 2006
|
Member |
|
|
Hi all,
I'm drawing some points on a sky map. The map is <0.5deg across, so I
thought a flat approximation would be ok...this may not be true. Around my
central point, I draw a circle using a flat geometry relation [x =
x0+r*cos(theta); y = y0 + r*sin(theta)], but i see a point outside the
circle that I expect to be inside.
I expect this point to be inside the circle, because the radius of the
circle (in arcsec) is *greater* than the great-circle angular distance
from the centre of the circle to the point.
So...I'm thinking that my flat-geometry assumption is false. My question:
can anyone point me towards forumlae/code that will calculate this circle
on the sky (i.e. all points which have a fixed great-circle distance from
the centre)? I'm using the astro library gcirc.pro to calculte my
great-circle angular distances...I kind of want the "inverse" of that
routine, I guess.
Or...is there a better way to do it? [Or maybe my bug is elsewhere?]
cheers
chris
|
|
|
|
| Re: circles on the sky [message #65933 is a reply to message #65879] |
Tue, 31 March 2009 15:45  |
pgrigis
Messages: 436
Registered: September 2007
|
Senior Member |
|
|
Christopher Thom wrote:
> Quoth Paolo:
>
>> Christopher Thom wrote:
>>> Quoth Kenneth P. Bowman:
>>>
>>>> In article <alpine.OSX.1.10.0903311335490.8491@kanangra.uchicago.edu>,
>>>> Christopher Thom <cthom@oddjob.uchicago.edu> wrote:
>>>>
>>>> > Given a co-ordinate position (ra/dec or lat/long), a direction (e.g an
>>>> > angle east of north, for instance), and a great circle angular distance,
>>>> > how do I compute the coordinate of the final position?
>>>>
>>>> LL_ARC_DISTANCE.
>>>>
>>>> What! That wasn't obvious? :-)
>>>>
>>>> (This function should be referenced in the manual page for MAP_2POINTS,
>>>> and vice versa.)
>>>
>>> AHA!!! Missed this one. Now, by just passing all azimuths 0 -> 360deg, i
>>> have the coordinates of the "circles" i'm trying to draw (where, by
>>> "circle", i mean "the set of all points that are r distance from my
>>> lon/lat").
>>
>> Is that significantly different than a circle with radius r drawn
>> in the projected map, if r is about 0.5 degree as yous said
>> in the original post?
>
> Well...I think so. Map projections continually confuse me, and getting
> them right in IDL confuses me even more! what I can say for sure is this:
>
> If i just calculate a cartesian circle, using the following code:
>
> theta = findgen(361)/!DRADEG
> xx = x0 + r*cos(theta)
> yy = y0 + r*sin(theta)
> plot, x0, y0
> oplot, xx, yy
>
> I get a very circular object in my plots, both on an equirectangular plot
> of points, as well as a projected map, made using map_set.
>
> BUT...if i now calculate the great circle distance to each of the 361
> points in my cartesian circle from the centre of the circle, the distance
> is NOT constant, as I expect. Rather, it is sinusoidal, approaching r at
> the maximum of the curve.
>
> OTOH, using ll_arc_distance gives me a rather egg-like "circle", but at
> least the distance from the centre to all the points on my "circle" is
> constant (i.e. r), as expected.
Well, I tried:
;go to the equator at central meridian
x=0
y=0
;0.5 degrees distance
arc_dist=0.5/360*!Pi*2
;define azimuts from 0 to 2Pi
az=findgen(100)/99*2*!Pi
;result coordinates
resx=az*0
resy=az*0
;ll_arc_distance seems not to be vectorized?
;or maybe I misread the docs
.run
FOR i=0,n_elements(az)-1 DO BEGIN
res=ll_arc_distance([x,y],arc_dist,az[i])
resx[i]=res[0]
resy[i]=res[1]
ENDFOR
end
;plot result coordinates
plot,resx,resy,/iso
;looks very circular to me
I did not do map projections because I'd like to keep
what remains of my sanity :) but I don't expect them
to distort that circle too much...
Ciao,
Paolo
>
> I must have spent 2 or 3 days digging through my code, convinced that I
> must have screwed up the object locations, rather than just the "drawing a
> circle" part.
>
> cheers
> chris
|
|
|
|
| Re: circles on the sky [message #65934 is a reply to message #65879] |
Tue, 31 March 2009 15:27 |
Christopher Thom
Messages: 66
Registered: October 2006
|
Member |
|
|
Quoth Paolo:
> Christopher Thom wrote:
>> Quoth Kenneth P. Bowman:
>>
>>> In article <alpine.OSX.1.10.0903311335490.8491@kanangra.uchicago.edu>,
>>> Christopher Thom <cthom@oddjob.uchicago.edu> wrote:
>>>
>>>> Given a co-ordinate position (ra/dec or lat/long), a direction (e.g an
>>>> angle east of north, for instance), and a great circle angular distance,
>>>> how do I compute the coordinate of the final position?
>>>
>>> LL_ARC_DISTANCE.
>>>
>>> What! That wasn't obvious? :-)
>>>
>>> (This function should be referenced in the manual page for MAP_2POINTS,
>>> and vice versa.)
>>
>> AHA!!! Missed this one. Now, by just passing all azimuths 0 -> 360deg, i
>> have the coordinates of the "circles" i'm trying to draw (where, by
>> "circle", i mean "the set of all points that are r distance from my
>> lon/lat").
>
> Is that significantly different than a circle with radius r drawn
> in the projected map, if r is about 0.5 degree as yous said
> in the original post?
Well...I think so. Map projections continually confuse me, and getting
them right in IDL confuses me even more! what I can say for sure is this:
If i just calculate a cartesian circle, using the following code:
theta = findgen(361)/!DRADEG
xx = x0 + r*cos(theta)
yy = y0 + r*sin(theta)
plot, x0, y0
oplot, xx, yy
I get a very circular object in my plots, both on an equirectangular plot
of points, as well as a projected map, made using map_set.
BUT...if i now calculate the great circle distance to each of the 361
points in my cartesian circle from the centre of the circle, the distance
is NOT constant, as I expect. Rather, it is sinusoidal, approaching r at
the maximum of the curve.
OTOH, using ll_arc_distance gives me a rather egg-like "circle", but at
least the distance from the centre to all the points on my "circle" is
constant (i.e. r), as expected.
I must have spent 2 or 3 days digging through my code, convinced that I
must have screwed up the object locations, rather than just the "drawing a
circle" part.
cheers
chris
|
|
|
|
| Re: circles on the sky [message #65935 is a reply to message #65879] |
Tue, 31 March 2009 14:20 |
pgrigis
Messages: 436
Registered: September 2007
|
Senior Member |
|
|
Christopher Thom wrote:
> Quoth Kenneth P. Bowman:
>
>> In article <alpine.OSX.1.10.0903311335490.8491@kanangra.uchicago.edu>,
>> Christopher Thom <cthom@oddjob.uchicago.edu> wrote:
>>
>>> Given a co-ordinate position (ra/dec or lat/long), a direction (e.g an
>>> angle east of north, for instance), and a great circle angular distance,
>>> how do I compute the coordinate of the final position?
>>
>> LL_ARC_DISTANCE.
>>
>> What! That wasn't obvious? :-)
>>
>> (This function should be referenced in the manual page for MAP_2POINTS,
>> and vice versa.)
>
> AHA!!! Missed this one. Now, by just passing all azimuths 0 -> 360deg, i
> have the coordinates of the "circles" i'm trying to draw (where, by
> "circle", i mean "the set of all points that are r distance from my
> lon/lat").
Is that significantly different than a circle with radius r drawn
in the projected map, if r is about 0.5 degree as yous said
in the original post?
Ciao,
Paolo
>
> thanks all for the help
> chris
|
|
|
|
| Re: circles on the sky [message #65936 is a reply to message #65879] |
Tue, 31 March 2009 14:07 |
Christopher Thom
Messages: 66
Registered: October 2006
|
Member |
|
|
Quoth Kenneth P. Bowman:
> In article <alpine.OSX.1.10.0903311335490.8491@kanangra.uchicago.edu>,
> Christopher Thom <cthom@oddjob.uchicago.edu> wrote:
>
>> Given a co-ordinate position (ra/dec or lat/long), a direction (e.g an
>> angle east of north, for instance), and a great circle angular distance,
>> how do I compute the coordinate of the final position?
>
> LL_ARC_DISTANCE.
>
> What! That wasn't obvious? :-)
>
> (This function should be referenced in the manual page for MAP_2POINTS,
> and vice versa.)
AHA!!! Missed this one. Now, by just passing all azimuths 0 -> 360deg, i
have the coordinates of the "circles" i'm trying to draw (where, by
"circle", i mean "the set of all points that are r distance from my
lon/lat").
thanks all for the help
chris
|
|
|
|
|
|
|
|
| Re: circles on the sky [message #66018 is a reply to message #65934] |
Tue, 07 April 2009 07:17 |
JDS
Messages: 94
Registered: March 2009
|
Member |
|
|
On Mar 31, 6:27 pm, Christopher Thom <ct...@oddjob.uchicago.edu>
wrote:
> Quoth Paolo:
>
>
>
>> Christopher Thom wrote:
>>> Quoth Kenneth P. Bowman:
>
>>>> In article <alpine.OSX.1.10.0903311335490.8...@kanangra.uchicago.edu>,
>>>> Christopher Thom <ct...@oddjob.uchicago.edu> wrote:
>
>>>> > Given a co-ordinate position (ra/dec or lat/long), a direction (e.g an
>>>> > angle east of north, for instance), and a great circle angular distance,
>>>> > how do I compute the coordinate of the final position?
>
>>>> LL_ARC_DISTANCE.
>
>>>> What! That wasn't obvious? :-)
>
>>>> (This function should be referenced in the manual page for MAP_2POINTS,
>>>> and vice versa.)
>
>>> AHA!!! Missed this one. Now, by just passing all azimuths 0 -> 360deg, i
>>> have the coordinates of the "circles" i'm trying to draw (where, by
>>> "circle", i mean "the set of all points that are r distance from my
>>> lon/lat").
>
>> Is that significantly different than a circle with radius r drawn
>> in the projected map, if r is about 0.5 degree as yous said
>> in the original post?
>
> Well...I think so. Map projections continually confuse me, and getting
> them right in IDL confuses me even more! what I can say for sure is this:
>
> If i just calculate a cartesian circle, using the following code:
>
> theta = findgen(361)/!DRADEG
> xx = x0 + r*cos(theta)
> yy = y0 + r*sin(theta)
> plot, x0, y0
> oplot, xx, yy
>
> I get a very circular object in my plots, both on an equirectangular plot
> of points, as well as a projected map, made using map_set.
>
> BUT...if i now calculate the great circle distance to each of the 361
> points in my cartesian circle from the centre of the circle, the distance
> is NOT constant, as I expect. Rather, it is sinusoidal, approaching r at
> the maximum of the curve.
That's not surprising. Only certain projections preserve shape in the
2D plane. Conformal projections like Mercator are among them (but in
that case relative areas are not preserved: the Greenland effect).
JD
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
