;
; Copyright (c) 1995, Research Systems, Inc.  All rights reserved.
;       Unauthorized reproduction prohibited.
;+
; NAME:
;       POLAR_CONTOUR
;
; PURPOSE:
;       Produce a contour plot from data in polar coordinates.
;       Data may be in a regular or scattered grid.
; CATEGORY:
;       Graphics.
;
; CALLING SEQUENCE:
;       POLAR_CONTOUR, Z, Theta, R
; INPUTS:
;       Z = data values.  If regulary gridded, Z must have dimensions of
;               (NTheta, NR).
;       Theta = values of theta in radians.  For the regular grid, Theta
;               must have the same number of elements as the first dimension
;               of Z.  For the scattered grid, Theta must have the same number
;               of elements as Z.
;       R = values of radius.  For the regular grid, R
;               must have the same number of elements as the second dimension
;               of Z.  For the scattered grid, R must have the same number
;               of elements as Z.
; KEYWORD PARAMETERS:
;       Any of the keywords accepted by CONTOUR may be specified.
;       In addition to:
;       SHOW_TRIANGULATION = color.  If set, the triangulation connecting
;               the data points is overplotted in the designated color.
;   nrings=number of rings for the polar style overlay
;   nspokes= number of spokes for polar overlay
;   ring_thick= thickness of rings
;       spoke_thick= thickness of spokes
;   labels= flag for placing angle and radius labels on overlay
;   rotate= rotation angle for the data to be contoured
;       charsize= size of characters on polar overlay
;   charthick= thickness of characters on polar overlay
;   ps = switch to turn postscript output on and off
; OUTPUTS:
;       A contour plot is produced.
; COMMON BLOCKS:
;       None.
; SIDE EFFECTS:
;       Plot is produced on the current graphics output device.
; RESTRICTIONS:
;       None.
; PROCEDURE:
;       The cartesian coordinates of each point are calculated.
;       TRIANGULATE is called to grid the data into triangles.
;       CONTOUR is called to produce the contours from the triangles.
; EXAMPLE:
;               Example with regular grid:
;       nr = 12         ;# of radii
;       nt = 18         ;# of thetas
;       r = findgen(nr) / (nr-1)   ;Form r and Theta vectors
;       theta = 2 * !pi * findgen(nt) / (nt-1)
;       z = cos(theta*3) # (r-.5)^2   ;Fake function value
;       tek_color
;               Create filled contours:
;       Polar_Contour, z, theta, r, /fill, c_color=[2,3,4,5]
;
;               Example with random (scattered) grid:
;       n = 200
;       r = randomu(seed, n)            ;N random r's and Theta's
;       theta = 2*!pi * randomu(seed, n)
;       x = cos(theta) * r              ;Make a function to plot
;       y = sin(theta) * r
;       z = x^2 - 2*y
;       Polar_Contour, z, theta, r, Nlevels=10, xrange=[0,1], yrange=[0,1]
;
; MODIFICATION HISTORY:
;       Written by:     Your name here, Date.
;       January, 1995.  DMS, RSI.
;       August, 1997 RLK
;-
pro plotSpokeTicks,finalRing,angle,color=spokeColor, ticklen=ticklen, spokeTicks=spokeTicks

;Utility routine to plot the ticks on the spokes

if finalRing eq 0 then return ;with no rings, cannot figure out ticks
if not keyword_set(ticklen) then ticklen=.01

intervalBetweenTicks = finalRing/float(spokeTicks)
for i=1,spokeTicks-1 do begin
  xcenter = (i * intervalBetweenTicks)*cos(angle)
  ycenter = (i * intervalBetweenTicks)*sin(angle)
  xTickBottom = ticklen * cos(angle-(!pi/2.))
  yTickBottom = ticklen * sin(angle-(!pi/2.))
  xTickTop = ticklen * cos(angle+(!pi/2.))
  yTickTop = ticklen * sin(angle+(!pi/2.))
  coord = convert_coord(xcenter,ycenter,/data,/to_normal)
  xTick = [coord[0]+xTickTop,coord[0]+xTickBottom]
  yTick = [coord[1]+yTickTop,coord[1]+yTickBottom]
  plots,xTick,yTick,/normal
endfor
return
end

;{{:|{{:|{{:|{{:|{{:|{{:|{{:|{{:|{{:|{{:|{{:|{{:|{{:|{{:|

PRO Polar_Contour, z, theta, r, SHOW_TRIANGULATION=show, _Extra = e,$
        nrings=nrings, nspokes=nspokes, ring_thick=ring_thick, $
        spoke_thick=spoke_thick,labels=labels,rotate=rotate,$
        charsize=charsize,charthick=charthick,ps=ps,ringColor=ringColor,$
        spokeColor = spokeColor, degrees=degrees, spokeStyle=spokeStyle, $
        spokeTicks=spokeTicks, ticklen=ticklen, labelColor=labelColor, $
        r_max=r_max, ringStyle=ringStyle, contourOnTop=contourOnTop

on_error, 2             ;Return to caller...
OS = !d.name
if keyword_set(ps) then begin
  file = dialog_pickfile(/write)
  set_plot,'ps'
  device,/landscape,filename=file
  device,/color
endif

s = size(z)
nz = n_elements(z)
r_max = 0 ;use as a kind of flag
twopi = !pi*2.

origTheta = theta

if not keyword_set(ring_thick) then ring_thick=1
if not keyword_set(spoke_thick) then spoke_thick=1
if not keyword_set(rotate) then rotate = 0.0
if n_elements(spokeColor) eq 0 then spokeColor = !d.n_colors-1
if n_elements(ringColor) eq 0 then ringColor = !d.n_colors-1
if n_elements(labelColor) eq 0 then labelColor = !d.n_colors-1
if keyword_set(degrees) then theta = theta *!dtor

theta = theta + (rotate*!dtor)

if (s(0) eq 2) and (n_elements(theta) eq s(1)) and $   ;Regular grid case?
   (n_elements(r) eq s(2)) then begin
  tt = theta # replicate(1,n_elements(r))
  rr = replicate(1,n_elements(theta)) # r
  x = rr * cos(tt)
  y = rr * sin(tt)
  irregular = 0
endif else begin                ;Irregular grid
  if n_elements(r) ne nz then message, $
        'R array has wrong number of elements'
  if n_elements(theta) ne nz then message, $
        'Theta array has wrong number of elements'
  x = r * cos(theta)
  y = r * sin(theta)
  irregular = 1
endelse

;old way FOR REALLY LARGE NUMBER OF POINTS THIS MAY STILL BE BEST.
;triangulate, x, y, tr
;newz = trigrid(x,y,z,tr,xgrid=xgrid,ygrid=ygrid)
;contour, newz, xgrid, ygrid, _EXTRA=e
contour, z, x, y, irregular=irregular, _extra=e, nodata=keyword_set(contourOnTop)

if n_elements(show) eq 1 then begin     ;Show the triangulation?
    triangulate, x, y, tr
    newz = trigrid(x,y,z,tr,xgrid=xgrid,ygrid=ygrid)
    oplot, x, y, /psym, color=show
    for i=0, n_elements(tr)/3-1 do begin
      t = [tr(*,i), tr(0,i)]
      plots, x(t), y(t), color=show
      endfor
    endif

if keyword_set(nrings) then begin
  ring = (findgen(360)/359.)*twopi
  if not keyword_set(r_max) then r_max = max(r)
  r_space = fix(r_max)/nrings
  if r_space le 0 then r_space=1
  r_max = fix(r_max/r_space) * r_space
  for i=0,r_max,r_space do begin
    oplot,i*cos(ring),i*sin(ring),thick=ring_thick,color=ringColor
    if ( keyword_set(labels) and (i gt 0.)) then begin
      xyouts,0,-i+(r_space/10.),strcompress(string(i,format='(f4.1)')),$
               charthick=charthick,charsize=charsize, color=labelColor
    endif
    finalRing = i
  endfor
endif

if keyword_set(nspokes) then begin
  theta_step = twopi/nspokes
  if not keyword_set(r_max) then r_max = max(r)
  for i=0.0,nspokes-1 do begin
    index = i*theta_step
    oplot,[0,r_max*cos(index)],[0,r_max*sin(index)],thick=spoke_thick,$
          color=spokeColor, linestyle=spokstyle
    if keyword_set(spokeTicks) then plotSpokeTicks,finalRing,index,color=spokeColor, $
                                    ticklen=ticklen, spokeTicks=spokeTicks
    if keyword_set(labels) then begin
      textval = ((index/!dtor)-rotate) mod 360.
      textval = strcompress(string(textval,format='(i5.1)'))
      just = 0.0
      if ( (index ge !pi/2.) and (index le !pi*1.5)) then begin
        angle = index + !pi
        just = 1.0
        textval = textval + ' '
      endif else angle = index
      xyouts,r_max*cos(index),r_max*sin(index),$
             orientation=angle/!dtor,textval,align=just,$
             charthick=charthick,charsize=charsize
    endif
  endfor
endif

if keyword_set(contourOnTop) then contour, z, x, y, /irregular, _extra=e,/over

if keyword_set(ps) then begin
  device,/close
  set_plot,OS
endif

theta = origTheta

return
end