Re: New ringmap - help needed.
Posted by Roger Hunter on April 01, 2002 at 05:24:50:

.02

No problem, here's the subroutine which creates the circles.

What happens is a gap once a circle gets big enough. Happens in all bigger rings and part is reflected, suggesting a sign error.

This is the same routine the other ringmaps use which really puzzles me. The only difference is that this routine does not convert to a map projection since I'm plotting lat-lon directly.


SUB create_circle (rlat,rlon,radius)
LOCAL i
!
! ----------------------------------------------------
! COMPUTE A CIRCLE OF RADIUS radius (km) AT CENTER rlat, rlon
LET radlat = rlat*dtr ! CENTER LAT IN RADIANS
LET radlon = rlon*dtr ! CENTER LON IN RADIANS
LET radterra = 6371.0 ! RADIUS OF THE EARTH (km)

LET theta = radius/radterra
LET st = sin(theta)
LET ct = cos(theta)

FOR i= 0 to 360
LET angl = i*dtr

LET sl = sin(radlon)
LET cl = cos(radlon)

LET cp = cos(radlat)
LET sp = sin(radlat)

LET x = st*cos(angl)
LET y = st*sin(angl)

LET xp = x*sp + ct*cp
LET zp = -x*cp + ct*sp
LET cpf = (1. - zp*zp)^0.5

LET vlat = asin(zp)*rtd ! LATITUDE

LET sinang = y/cpf
LET ang = asin(sinang)

! DETERMINE QUADRANTS

! FIRST OR THIRD QUADRANT
IF xp < 0 then
IF y < 0. then LET ang = abs(ang) + pi
IF y >= 0. then LET ang = pi - ang
END IF
!
LET ang = ang + radlon
IF ang > 2.*pi then LET ang = ang - 2.*pi

LET ulon = ang*rtd ! LONGITUDE
IF ulon > 180.0 then LET ulon = ulon - 360.0

LET cords(i,1) = ulon
LET cords(i,2) = vlat
NEXT i

END SUB
!

Roger