undef("create_grid_azimuthal")
function create_grid_azimuthal(xx:numeric, clat[*]:numeric, clon[*]:numeric, phi[*]:numeric, radius[*]:numeric, opt[1]:logical)
;
; A conceptual sketch ...; A pic:
; https://cs.stackexchange.com/questions/43744/fast-algorithm-for-interpolating-data-from-polar-coordinates-to-cartesian-coordi
; This is a 'polar' (not cylindrical)  view (z=constant;  p=isobaric_level)
;
; A possible issue is that the further away from the center the 'spatial' the locations are further apart.
;
;
local dimxx, rankxx, ntim, nlev, nlat, mlon, nctr, nrad, nphi, rlat, rlon, drad
;
; Nomenclature
; xx:          source rectininear grid: (time,lat,lon),  (time,lev,lat,lon)
; clat, clon:  center latitudes(s) / longitude(s); degrees_north , degrees_east
; phi:         angles   
; radius:      one or more radii (degrees)
; opt:         option. Currently, not used. Set to False
;
; Return
begin
  dimxx = dimsizes(xx)
  rankxx= dimsizes(dimxx)
  if (rankxx.lt.3 .or. rankxx.gt.4) then
      print("create_grid_azimuthal: Only ranks 3, 4 are allowed: rankxx="+rankxx)
      exit
  end if

  nctr = dimsizes(clat)
  nrad = dimsizes(radius)
  nphi = dimsizes(phi)

  dnam = getvardimnames(xx)
  ntim = dimxx(0)

  if (rankxx.eq.3) then
      nlev = 0
      nlat = dimxx(1)
      mlon = dimxx(2)
      lat  = xx&$dnam(1)$
      lon  = xx&$dnam(2)$
  else
      nlev = dimxx(1)
      nlat = dimxx(2)
      mlon = dimxx(3)
      plev = xx&$dnam(1)$
      lat  = xx&$dnam(2)$
      lon  = xx&$dnam(3)$
  end if
 
  rlat    = new( (/ntim,nctr,nrad,nphi/)  , "float")    ; lat at each radius, center and time
  rlon    = new( (/ntim,nctr,nrad,nphi/)  , "float")    ; lon

   if (rankxx.eq.3) then
       drad = new( (/ntim,nctr,     nrad,nphi/), "float")   ; xx=>3D
   else
       drad = new( (/ntim,nctr,nlev,nrad,nphi/), "float")   ; xx=>4D
   end if
   
;---Loop over time & radius
;---At each iteration use some builtin functions to interpolate. Could be slow!
 
   do nt=0,ntim-1
     do nc=0,nctr-1
       do nr=0,nrad-1
          nggcog(clat(nc), clon(nc), radius(nr), rlat(nt,nc,nr,:), rlon(nt,nc,nr,:))
          if (rankxx.eq.4) then
              do np=0,nlev-1
                 drad(nt,nc,np,nr,:) = linint2_points(lon, lat, xx(nt,np,:,:), False \
                                                     ,rlon(nt,nc,nr,:), rlat(nt,nc,nr,:), 0) ; 4D rectilinear
              end do                ; np
          else
              drad(nt,nc,nr,:) = linint2_points(lon, lat, xx(nt,:,:), False \
                                               ,rlon(nt,nc,nr,:), rlat(nt,nc,nr,:), 0) ; 3D rectilinear
          end if
       end do                  ; nr
     end do                    ; nc
   end do                      ; nt

  drad!0 = dnam(0)
  drad&$dnam(0)$ = xx&$dnam(0)$
  if (rankxx.eq.3) then
      drad!1 = "center"
      drad!2 = "radius"
      drad!3 = "phi"
      drad&radius = radius
      drad&phi    = phi
  else
      drad!1 = "center"
      drad!2 = dnam(1)
      drad!3 = "radius"
      drad!4 = "phi"
      drad&$dnam(1)$ = plev
      drad&radius    = radius
      drad&phi       = phi
  end if
  if (isatt(xx,"long_name")) then
      drad@long_name = xx@long_name
  else
      drad@long_name = "grid to azimuth"
  end if
  if (isatt(xx,"units")) then
      drad@units = xx@units
  end if
  drad@NCL_tag = "create_grid_azimuthal"

  return(drad)
end  
;================================================================================================
;                               MAIN
;================================================================================================

;---Bogus input
;---Add minimal meta data to facilitate information tracking

   clat  = (/ 25.0, 40.0 /)    ; center lat 
   clon  = (/ 45.0, 65.0 /)    ;        lon
   clat!0= "center_lat"
   clon!0= "center_lon"
   nctr  = dimsizes(clat)
   printVarSummary(clat)       ; [center_lat | 2]
   print("-----")
   printVarSummary(clon)       ; [center_lon | 2]
   print("-----")

   radius= (/ 10,20,25,35,50/)*1.0   ; radii (degrees)
   radius!0 = "radius"
   radius@units = "degrees"
   nrad  = dimsizes(radius)
   printVarSummary(radius)     ; [radius | 5]
   print("-----")

   ntim  = 6
   time  = fspan(0,100,ntim)   ; [time | 6]
   time!0= "time"
   time@units = "fantasy"

   nlev    = 11                  ; pressure levels
   plev  = fspan(1000,100,nlev)
   plev!0= "plev"
   plev@long_name = "pressure"
   plev@units = "hPa"

   phi   = ispan(0,315,45)*1.0 ; nphi = [0,45,90,135,180,225,270,315] 
   nphi  = dimsizes(phi)       ; # of points at each radii

   Nlat  = 9                   ; rectilinear grid sizes
   Nlon  = 13                  ;                     
   lat   = fspan(20, 60, Nlat) ; rectilinear grid latitudes
   lon   = fspan(30,150, Nlon) ;                  longitudes
   lat!0 = "lat"
   lon!0 = "lon"

   xx    = random_uniform(-8,100, (/ntim,nlev,Nlat,Nlon/))  ; data: rectilinear grid
   xx!0  = "time"
   xx!1  = "plev"
   xx!2  = "lat"
   xx!3  = "lon"
   xx&time = time
   xx&plev = plev
   xx&lat  = lat 
   xx&lon  = lon 
   xx@long_name = "Bogus Data"
   xx@units     = "Whatever"
   printVarSummary(xx)              ; [time | 6] x [plev | 11] x [lat | 9] x [lon | 13]
   printMinMax(xx,0)
   print("----")

   drad = create_grid_azimuthal(xx,clat,clon,phi,radius,False)

   printVarSummary(drad)
   printMinMax(drad,0)
   print("----")


;---Average all 'phi' (azimuth) at each time, each center, pressure levels and radius

   rankxx = dimsizes(dimsizes(xx))
   dradAvg = dim_avg_n_Wrap(drad, rankxx)
   printVarSummary(dradAvg)   
   printMinMax(dradAvg,0)

   wks  = gsn_open_wks("png","azimuthal_avg")   ; send graphics to PNG file

  res                     = True          ; plot mods desired
  res@cnFillOn            = True          ; turn on color
  res@gsnMaximize         = True          ; maximize plot size
  res@trYReverse          = True
  res@tiYAxisString      = "Pressure (hPa)"                   
  res@tiXAxisString      = "Radius"  
  res@lbOrientation      = "vertical"

  nt = 0
  nc = 0
  do nt=0,0              ; ntim-1
    do nc=0,0            ; nctr-1
       res@tiMainString = "create_grid_azimuthal: nt="+time(nt)+": ["+sprintf("%5.1f",clat(nc)) \
                                                                 +","+sprintf("%5.1f",clon(nc))+"]"
       plot = gsn_csm_contour(wks,dradAvg(nt,nc,:,:),res)       ; contour the variable
    end do
  end do