; Concepts illustrated:
;   - Specifying the central lat and lon
;   - Specifying the radial distances of points to be interpolated
;   - Use local function 'rgrid2geocircle' located in 'grid2geolocation'
;     to get the radial locations via 'nggcog'; use bilinear
;     interpolation (linint2_points) to get the values at various radii.
;   - Plotting the results     
;     Contour radial values         

begin
;-------------------------------------------------
load "./grid2geocircle.ncl"        ; radial interpolation library
;================================================================================================
;                               MAIN
;================================================================================================

;---Read in TRMM precipitation data

 f   = addfile("aila.nc", "r")
 prc = f->pcp(:,{0:30},{75:100})    ; prc[time | 32] x [latitude | 120] x [longitude | 100]
 
 ;printVarSummary(prc)
   
  dimprc  = dimsizes(prc)
  rankprc = dimsizes(dimprc) 
  ntim  = dimprc(0)
   
;---Specify center locations 

 filename = "aila_track.txt"
 ncols = numAsciiCol(filename)                   
 data_ = readAsciiTable(filename,ncols,"float",0) 
 
   clat=data_(:,0)  ; Latitude:  storm centers
   clon=data_(:,1)  ; Longitude: storm centers 
 
   clat!0= "center_lat"
   clon!0= "center_lon"
   nctr  = dimsizes(clat)

   ;rval = css2c(clat, clon)
   ;print(rval)
   ;printVarSummary(rval)
   ;x = rval(0,:)
   ;y = rval(1,:)
   ;print(x)
   ;print(y)
  
 
;---Specify radii: can be unequally spaced but must be monotonically increasing

   nrad             = 17
   radius           = fspan(0,300,nrad)  ; radii (kilometers)
   radius!0         = "radius"
   radius@units     = "kilometers"
   radius@long_name = "Radii"
   nrad             = dimsizes(radius) ; 17

;---Calculate the interpolated value and lat/lon locations of each circle at each radius
    
   N = 180
   prc_geocircle = rgrid2geocircle(prc, clat, clon, radius, N, False)
   ;printVarSummary(prc_geocircle)
   

;---Clarity: Explicitly extract the returned variables from the returned list variable

   prcrad = prc_geocircle[0]
   plat   = prc_geocircle[1]           ; radial [azimuthal] latitudes 
   plon   = prc_geocircle[2]           ;   "        :       longutudes
   ;printVarSummary(prcrad)            ; [time | 32] x [center | 22] x [radius | 17] x [circle | 180]
   ;printVarSummary(plat)              ; [center | 22] x [radius | 17] x [circle | 180]          
  
  
  plat_cen_avg = dim_avg_n_Wrap(plat,0)
  plon_cen_avg = dim_avg_n_Wrap(plon,0)
  

  time_avg = dim_avg_n_Wrap(prcrad,0)
  ;printVarSummary(time_avg)
  
  cen_avg  = dim_avg_n_Wrap(time_avg,0)
  ;printVarSummary(cen_avg)
   
  
  wks = gsn_open_wks("png","cartesian")           ; send graphics to PNG file

  cnres                   = True
  cnres@gsnMaximize       = True
  cnres@sfXArray          = plat_cen_avg
  cnres@sfYArray          = plon_cen_avg
  cnres@cnFillOn          = True
  cnres@cnFillPalette     = "rainbow"         ; set color map
  cnres@cnLinesOn         = False
  ;cnres@cnFillMode        = "RasterFill"      ; this mode is fastest
  ;cnres@trGridType        = "TriangularMesh"

  contour = gsn_csm_contour(wks,cen_avg,cnres)
end