############################################################################
#  Patrick Marsh
#  23 February 2008
#
#  NAME:  circlepy.py
#
#  DESCRIPTION:  This script generates a circle placefile for use with GRx
#                Radar viewers.  The user specifies the center of the
#                circle, the radius of the circle in kilometers, and the 
#                number of points along the circle to generate.  The output
#                is a simple text file.  Please remember, the more points
#                on the circle, the longer you will need to let the script
#                run to compute the points.  If you open the text file and
#                the last line is not "End:" you need to re-run the script
#                and wait a little longer before opening the text file.
#
#                The default THRESHOLD is 999
#
#                The default COLOR is 255 255 255
#
#                Unless you are familiar with PYTHON, please do not attempt
#                to change anything below the first block of code.
#
############################################################################

import numpy
from pylab import *


# DECLRATIONS THAT CAN BE CHANGED
latcent    = 34.62389        # LATITUDE OF CIRCLE CENTER
loncent    = -98.270833      # LONGITUDE OF CIRCLE CENTER
radius     =   40.           # RADIUS IN KILOMETERS (END WITH A .)
numpoints  = 1000.           # NUMBER OF POINTS TO USE (END WITH A .)
f = open("circle.txt","w")   # REPLACE circle WITH FILENAME OF CHOICE







# DECLARATIONS THAT CANNOT BE CHANGED
multfactor = numpoints / radius
halfpoints = numpoints / 2
doublepoints = (numpoints) * 2


# ARRAY INITIALIZATIONS
xvalues = []; xvalues = resize(array(xvalues,'f'),(doublepoints+1)) 
yvalues = []; yvalues = resize(array(yvalues,'f'),(doublepoints+1))
lat = []; lat = resize(array(lat,'f'),(doublepoints+1)) 
lon = []; lon = resize(array(lon,'f'),(doublepoints+1)) 


# SET INITIAL VALUES FOR ARRAYS
yvalues[0] = -1.
xvalues[0] = 0.


# COMPUTE TOP HALF OF A CIRCLE NORMALIZED TO RADIUS 1
for p  in range(1,numpoints+1):
  yvalues[p] = -1. + (radius / numpoints) * ((multfactor * p) / halfpoints) 
  xvalues[p] = ((1. - yvalues[p]**2)**0.5)
    

# COMPUTE BOTTOM HALF OF CIRCLE NORMALIZED TO RADIUS 1
for p in range(numpoints+1, doublepoints+1):
  m = p - numpoints 
  yvalues[p] = (-1. + (radius / numpoints) * ((multfactor * m) / halfpoints)) * -1. 
  xvalues[p] = ((1. - yvalues[p]**2)**0.5) * -1.


############################################################################
# CONVERT NORMALIZED CIRCLE TO LAT LONS (BASED OFF OF LATCENT AND LONCENT)
# AND THEN WRITE OUT LAT LONS IN LINE FORMAT FOR USE IN GRLEVELX
############################################################################
f.write("; GRLevelX placefile creating range ring of %2ikm centered at %3f ,%3f\n\n" % (radius,latcent,loncent))
f.write("Title: RANGE RING\n")
f.write("Threshold: 999\n")
f.write("Color: 255 255 255\n\n")
f.write("Line: 1, 0,\"\"\n")
for p in range(0,doublepoints+1):
  lat[p] = latcent + (yvalues[p] * (radius / 111.325))
  lon[p] = loncent + (xvalues[p] * (radius / (numpy.cos(lat[p] * numpy.pi / 180) * 111.325)))
  f.write("  %3f, %3f\n" % (lat[p],lon[p]))
f.write("End:\n")
f.write(" ")
f.close