;*****************************
;  Function to plot a circle
;*****************************
function circle, xcenter,ycenter,radius
  points = (2 * !PI / 99.0) * FINDGEN(100)
  x = xcenter + radius * COS(points )
  y = ycenter + radius * SIN(points )
  RETURN, TRANSPOSE([[x],[y]])
end

;*****************************
;  Program to plot cband radar
;*****************************
pro read_cband,julian_day_ship,dbz_grid,dbz_mask,x_out,y_out,$
  nbins,rstart,rscale,julian_day,elangle
  
;*** Color table ***
;  Top is the last color to scale 256 colors, 0-255
top_color=252
;  Colortable  0-252  253=white
mytable=colortable(39,ncolors=254)
;254=hot pink               ;gray=255
mytable=[mytable,transpose([238,18,137]),transpose([150,150,150])]
cbtable=mytable[0:top_color,*]  
  
;*******************  
;  Plots the first scan only, nelev=1
;  To plot other scans change nelev at line 69
;julian_day_ship=julday(1,25,2018,19,22,29)  ;test
;julian_day_ship=julday(1,28,2018,16,44,36)
;  directory
filedir='/uufs/chpc.utah.edu/common/home/mace-group4/Capricorn2/CBand/'
caldat,julian_day_ship,smm,sdd,syy,shh,smi,sss
syy=string(syy,format='(I4)')
smm=string(smm,format='(I02)')
sdd=string(sdd,format='(I02)')
shh=string(shh,format='(I02)')
smi=string(smi,format='(I02)')
sss=string(sss,format='(I02)')
date_str=syy+smm+sdd+shh+smi+sss
ymd_str=syy+smm+sdd

;  File search string
;filestr='9776HUB-PPIVol-'+syy+smm+sdd+'-'+shh+'*hdf'
filestr='9776HUB-PPIVol-'+syy+smm+sdd+'*hdf'
;filestr='9776HUB-PPIVol-*hdf'

;  file name
filenames=file_search(filedir+ymd_str+'/'+filestr,count=numfiles)
finfo=file_info(filenames)
result_size=where(finfo.size gt 2e7,count_size)
filenames=filenames[result_size]
numfiles=count_size

;  Get the times of the file start and find the closest to ship time
time_str=strmid(file_basename(filenames),24,6)
fhh=strmid(time_str,0,2)
fmi=strmid(time_str,2,2)
fss=strmid(time_str,4,2)
julian_day_files=julday(smm,sdd,syy,fhh,fmi,fss)
closest=min(abs(julian_day_files-julian_day_ship),idx)
;  Closest must be in 2 hours
if closest le 2D/24D then begin

caldat,julian_day_files[idx],xmm,xdd,xyy,xhh,xmi,xss
print,xyy,xmm,xdd,xhh,xmi,xss
numfiles=1
kk1=idx
;for kk=0,numfiles-1 do begin
for kk=kk1,kk1 do begin
  filename=filenames[kk]
  print,filename
  file_tag=strmid(file_basename(filename),16,19)
  
  ;  dump the hdf structure
  presult=h5_parse(filename)

  ;  Open an existing hdf5 file
  file_id=h5f_open(filename)

  ;  Open the top level (root)
  root_id=h5g_open(file_id,'/')

  ; Root / where (lon, lat, height)
  group_id=h5g_open(root_id,'where')
  attr_id=h5a_open_name(group_id,'lon')  ;154.38200
  lon=h5a_read(attr_id)
  attr_id=h5a_open_name(group_id,'lat')  ;-26.617001
  lat=h5a_read(attr_id)
  attr_id=h5a_open_name(group_id,'height')  ;22.00
  height=h5a_read(attr_id)
  h5a_close,attr_id
  h5g_close,group_id
  ;print, 'lon/lat/ht = ', lon, lat, height

  ; Root / what (object, version, date, time, source)
  group_id=h5g_open(root_id,'what')
  attr_id=h5a_open_name(group_id,'source')  ;WMO:00000,RAD:,PLC:CSIRO,NOD:9776HUB,ORG:0,CTY:000,CMT:
  csgn=h5a_read(attr_id)
  h5a_close,attr_id
  h5g_close,group_id
  ;print, 'station ID = ', csgn

  ; Root / how (wavelength, beamwidth)
  group_id=h5g_open(root_id,'how')
  attr_id=h5a_open_name(group_id,'wavelength')  ;5.3399911
  wavelength=h5a_read(attr_id)
  attr_id=h5a_open_name(group_id,'beamwH')  ;1.3
  beamwH=h5a_read(attr_id)
  h5a_close,attr_id
  h5g_close,group_id
  ;print, 'wavelength,beamwH = ', wavelength, beamwH

  ; Set up the arrays to hold the data
  ;nelev=17  ;how many scans to read, total number scans is 17
  ;nelev=12  ;new file total number of scans
  nelev=1  ; first level scan
  startdate=strarr(nelev) & starttime=strarr(nelev) & endtime=strarr(nelev)
  dataset_n=strarr(nelev) & elangle=dblarr(nelev) & nbins=intarr(nelev)
  rstart=dblarr(nelev) & rscale=dblarr(nelev) & nrays=intarr(nelev)
  a1gate=intarr(nelev) & ray_spacing=fltarr(nelev) & sjulian_day=dblarr(nelev)
  ejulian_day=dblarr(nelev) & delta_t=dblarr(nelev)

  ;  Initialize gate_start to -0.5
  gate_start=make_array(/double,nelev,value=-0.5)

  totdata=17  ;total number of data arrays in a dataset
  data_n=strarr(totdata)
  quantity=strarr(totdata)
  quality_n=strarr(totdata)
  gain=dblarr(totdata)
  offset=dblarr(totdata)
  nodata=dblarr(totdata)
  undetect=dblarr(totdata)

  ;  Loop through the datasets at the root level
  for i=0,nelev-1 do begin
  
    dnum=strcompress(i+1,/remove_all)
  
    ; Root / datasetX
    dataset_n[i]='dataset'+dnum
    ;print,'opening ',dataset_n[i]
    group_id=h5g_open(root_id,dataset_n[i])

    ; ... finding start/end date/time
    g_id=h5g_open(group_id,'what')
    attr_id=h5a_open_name(g_id,'startdate')  ;20170822
    startdate[i]=h5a_read(attr_id)
    ;print, 'startdate:',startdate[i]

    attr_id=h5a_open_name(g_id,'starttime')  ;205654
    starttime[i]=h5a_read(attr_id)
    ;print, 'starttime:',starttime[i]

    attr_id=h5a_open_name(g_id,'endtime')  ;205730
    endtime[i]=h5a_read(attr_id)
    ;print, 'endtime:',endtime[i]

    h5a_close,attr_id
    h5g_close,g_id
    h5g_close,group_id

    ; Constructing the time array 
    syear=strmid(startdate[i],0,4)
    smonth=strmid(startdate[i],4,2)
    sday=strmid(startdate[i],6,2)
    shour=strmid(starttime[i],0,2)
    smin=strmid(starttime[i],2,2)
    ssec=strmid(starttime[i],4,2)
    ehour=strmid(endtime[i],0,2)
    emin=strmid(endtime[i],2,2)
    esec=strmid(endtime[i],4,2)
    ;date_str=syear+smonth+sday+shour+smin+ssec
  
    sjulian_day[i]=julday(smonth,sday,syear,shour,smin,ssec)
    ejulian_day[i]=julday(smonth,sday,syear,ehour,smin,ssec)
  
    ;print,'date and start and end time:',startdate[i],' ',starttime[i],' ',endtime[i]

    ; Root / datasetX / where to read (elangle, nbins, rstart, rscale, nrays, a1gate)
    group_id=h5g_open(root_id,dataset_n[i])
    g_id=h5g_open(group_id,'where')
  
    attr_id=h5a_open_name(g_id,'elangle')  ;dataset1=0.5
    elangle[i]=h5a_read(attr_id)
    h5a_close,attr_id

    attr_id=h5a_open_name(g_id,'nbins')  ;dataset1=960
    nbins[i]=h5a_read(attr_id)
    h5a_close,attr_id

    attr_id=h5a_open_name(g_id,'rstart')  ;dataset1=0
    rstart[i]=h5a_read(attr_id)
    h5a_close,attr_id

    attr_id=h5a_open_name(g_id,'rscale')  ;dataset1=125
    rscale[i]=h5a_read(attr_id)
    h5a_close,attr_id

    attr_id=h5a_open_name(g_id,'nrays')  ;dataset1=360
    nrays[i]=h5a_read(attr_id)
    h5a_close,attr_id

    attr_id=h5a_open_name(g_id,'a1gate')  ;dataset1=256
    a1gate[i]=h5a_read(attr_id)
    h5a_close,attr_id

    delta_t[i]=(ejulian_day[i]-sjulian_day[i])/nrays[i]
    ;print,i,elangle[i],nbins[I],rstart[I],rscale[I],nrays[I],a1gate[I]
    ;print, 'dataset ', i, ' nbins,nrays = ', nbins[I],nrays[i]
    h5g_close,g_id
  
    ; Root / datasetX / where to read (elangle, nbins, rstart, rscale, nrays, a1gate)
    group_id=h5g_open(root_id,dataset_n[i])
    g_id=h5g_open(group_id,'how')

    attr_id=h5a_open_name(g_id,'startazA')  ;dataset1=0.5
    startaza=h5a_read(attr_id)
    h5a_close,attr_id
    xstr='startaza'+dnum+'=startaza'
    result=execute(xstr)
  
    attr_id=h5a_open_name(g_id,'stopazA')  ;dataset1=960
    stopaza=h5a_read(attr_id) 
    h5a_close,attr_id
    xstr='stopaza'+dnum+'=stopaza'
    result=execute(xstr)  
  
    h5g_close,g_id
 
    ;  loop through data arrays in dataset
    for j = 0,totdata-1 do begin  
    
      dnumj=strcompress(j+1,/remove_all)
      data_n[j]='data'+dnumj
      ;print,'OPENING data=',j,data_n[j],'dnumj',dnumj

      ; Root / datasetX / dataX to read the data
      data_id=h5g_open(group_id,data_n[j])

      ; Grab the quanity value of the datsetX.dataX array
      ; Root / datasetX / dataX / what (gain, offset, nodata, undetect, quantity)
      g_id=h5g_open(data_id,'what')
      attr_id=h5a_open_name(g_id,'gain')  ;0.5
      gain[j]=h5a_read(attr_id)
      attr_id=h5a_open_name(g_id,'offset')  ;-32
      offset[j]=h5a_read(attr_id)
      attr_id=h5a_open_name(g_id,'nodata')  ;255
      nodata[j]=h5a_read(attr_id)
      attr_id=h5a_open_name(g_id,'undetect')  ;0
      undetect[j]=h5a_read(attr_id)    
      attr_id=h5a_open_name(g_id,'quantity')
      quantity[j]=h5a_read(attr_id)
      ;print, 'quantity ',j,quantity[j]   
      h5a_close,attr_id
      h5g_close,g_id
         
      ;  Continue if this is a quantity we want to plot
      if quantity[j] eq 'DBZH' then begin 
        print, 'quantity ',i,j,' ',quantity[j]    
        ; read data
        dset_id=h5d_open(data_id,'data')
        dataX1=float(h5d_read(dset_id))
        h5d_close,dset_id ; data
        h5g_close,data_id ; dataX    
    
        ;  apply gain and offset
        result=where(dataX1 eq nodata[j],count)  ;no values in dataX1=missing
        if count gt 0 then dataX1[result] = -8888.
        result=where(dataX1 eq undetect[j],count)  ;0 values = undetect=no cloud
        if count gt 0 then dataX1[result] = -9999.
        result=where(dataX1 gt -8888,count)
        if count gt 0 then dataX1[result]=dataX1[result]*gain[j]+offset[j]
      
        ;  Separate the arrays by the quality string
        if strmid(quantity[j],0,5) eq 'VRADH' then begin
          xstr='dop'+dnum+'=dataX1'
          result=execute(xstr)
        endif else if strmid(quantity[j],0,4) eq 'DBZH' then begin
          xstr='ref'+dnum+'=dataX1'
          result=execute(xstr)
        endif else if strmid(quantity[j],0,5) eq 'CLASS' then begin
          xstr='radar_flag'+dnum+'=dataX1'
          result=execute(xstr)             
        endif else if strmid(quantity[j],0,5) eq 'WRADH' then begin
          xstr='sw'+dnum+'=dataX1'
          result=execute(xstr)
        endif ;else begin
        ;stop
      endif  ;end of this is a quantity I want
    endfor  ;end  j loop through data arrays in dataset
    h5g_close,group_id ;close dataX  
 
    ; ray spacing is 1 degree  
    ray_spacing[i] = 360. / nrays[i]
    ;print,'ray spacing:',ray_spacing[I]
  
  endfor  ;end of i loop  loop of data arrays in data set
  h5g_close,root_id
  ;h5f_close,file_id
  ;print,'finished reading in data'

  do_plot1='no'
  if do_plot1 eq 'yes' then begin
    ;  Dimensions of the arrays as pulled out of the hdf file
    s1_orig=size(ref1,/dimensions)
    s2_orig=size(ref2,/dimensions)

    ;  Plot and histogram dbz data straight from the file
    ;  Set up the plot
    pxdim=900 & pydim=900
    ;  Position the plots
    xl=0.10 & xr=0.88
    yb=0.06 & yt=0.95
    sx=0.09
    sy=0.07
    numplots_x=2
    numplots_y=3
    position_plots,xl,xr,yb,yt,sx,sy,numplots_x,numplots_y,pos
    ;Clear plot space
    pnum=0
    p0=plot([0,1],[0,1],position=pos[pnum,*],/buffer,dimensions=[pxdim,pydim],$
      axis_style=4,/nodata)
    d=p0.convertcoord([pos[pnum,0],pos[pnum,2]],[pos[pnum,1],pos[pnum,3]],/normal,/to_device)
    isx=(d[0,1,0]-d[0,0,0])
    isy=(d[1,1,0]-d[1,0,0])

    pnum=0
    data_var=ref1
    result=where(data_var ne -9999,count)
    start_bin=min(data_var[result])
    end_bin=max(data_var[result])
    dbin=(end_bin-start_bin)/100.
    hist_generic,data_var,start_bin,end_bin,dbin,bins1,data_freq1,data_counts1
    p1=plot(bins1,data_counts1,/current,color='black',/hist,position=pos[pnum,*],$
      ytitle='counts',xtitle='ref1')
    t1=text(pos[pnum,0],pos[pnum,3],'total cloudy='+string(total(data_counts1)))
    
    pnum=2
    data_image=bytscl(data_var,top=top_color,min=start_bin,max=end_bin)
    data_var=data_image
    result=where(data_var ne -9999,count)
    start_bin=3;min(data_var[result])  ;eliminate 0
    end_bin=max(data_var[result])
    dbin=(end_bin-start_bin)/100.
    hist_generic,data_var,start_bin,end_bin,dbin,bins1,data_freq1,data_counts1
    p1=plot(bins1,data_freq1,/current,color='black',/hist,position=pos[pnum,*],$
      ytitle='frequency',xtitle='bytscaled data')

    pnum=4
    result=where(ref1 eq -9999,count)
    data_image[result]=253
    result=where(ref1 ne -9999,count)
    data_image[result]=70
    p1=image(data_image,/current,image_dimensions=[isx,isy],position=pos[pnum,*],rgb_table=mytable)
    c1=contour(data_var,ytitle='Y (km)',xtitle='X (km)',$
      xstyle=1,ystyle=1,font_size=10,/nodata,/current,position=pos[pnum,*],$
      irregular=0,xtickdir=1,ytickdir=1,font_name='Helvetica',$
      axis_style=2)

    pnum=1
    data_var=ref2
    result=where(data_var ne -9999,count)
    start_bin=min(data_var[result])
    end_bin=max(data_var[result])
    dbin=(end_bin-start_bin)/100.
    hist_generic,data_var,start_bin,end_bin,dbin,bins1,data_freq1,data_counts1
    p1=plot(bins1,data_counts1,/current,color='black',/hist,position=pos[pnum,*],$
      ytitle='counts',xtitle='ref2')

    pnum=3
    data_image=bytscl(data_var,top=top_color,min=start_bin,max=end_bin)
    data_var=data_image
    result=where(data_var ne -9999,count)
    start_bin=3;min(data_var[result])
    end_bin=max(data_var[result])
    dbin=(end_bin-start_bin)/100.
    hist_generic,data_var,start_bin,end_bin,dbin,bins1,data_freq1,data_counts1
    p1=plot(bins1,data_freq1,/current,color='black',/hist,position=pos[pnum,*],$
      ytitle='frequency',xtitle='bytscaled data')

    pnum=5
    result=where(ref2 eq -9999,count)
    data_image[result]=253
    result=where(ref2 ne -9999,count)
    data_image[result]=70
    p1=image(data_image,/current,image_dimensions=[isx,isy],position=pos[pnum,*],rgb_table=mytable)
    c1=contour(data_var,ytitle='Y (km)',xtitle='X (km)',$
      xstyle=1,ystyle=1,font_size=10,/nodata,/current,position=pos[pnum,*],$
      irregular=0,xtickdir=1,ytickdir=1,font_name='Helvetica',$
      axis_style=2)
    p1.save,'dbz_read_'+file_tag+'.png',height=pydim
  endif
;****************************
;  Done reading from HDF file, now do coordinate conversion
;****************************

  ;  Convert azimuth, elevation, range to x,y,z

  ;  Loop through the scans
  for i=0,nelev-1 do begin  ;scans
    dnum=strcompress(i+1,/remove_all)
    ;print,dnum,'dnum'
    n=0L  ;number of points in scan

    ;  Define arrays at this level
    numcount=long(nrays[i])*long(nbins[i])
    azimuth=fltarr(numcount)
    degree=fltarr(numcount)
    elevation=fltarr(numcount)
    range=fltarr(numcount)
    x=fltarr(numcount)
    y=fltarr(numcount)
    z=fltarr(numcount)
    dbz=fltarr(numcount)
    distance=fltarr(numcount)
    julian_day=dblarr(numcount)
 
    for j=0,nrays[i]-1 do begin  ;rays
      xstr='startaza=startaza'+dnum
      result=execute(xstr)
      xstr='stopaza=stopaza'+dnum 
      result=execute(xstr)
      if startaza[j] ge 350 and stopaza[j] le 2 then begin
        daza=(stopaza[j]+360.0-startaza[j])/nbins[i]  
        ;print,'in here'    
      endif else begin
        daza=(stopaza[j]-startaza[j])/nbins[i]
      endelse
      ;print, startaza[j],stopaza[j],nbins[i],daza

      for k=0,nbins[i]-1 do begin  ;bins
      
        ;  azimuth, 0degNorth
        ;azimuth[n] = 90.0-(float(j)*ray_spacing[i] + gate_start[i]) ;my old way
        ;azimuth[n] = (float(j)*ray_spacing[i] + gate_start[i]) ;degrees
        ;azimuth[n] = float(j)*ray_spacing[i]
        azimuth[n]=startaza[j]+float(k)*daza
     
        ;  Convert azimuth to math coordinate degree
        degree[n]=450.0-azimuth[n]
        if degree[n] gt 360.0 then degree[n]=degree[n]-360.0
 
        elevation[n]=elangle[i]   
        range[n]=rstart[i]+float(k)*rscale[i]  
        ;  height
        height=(float(k)*rscale[i])*sin(elangle[i]*(!pi/180.))+rstart[i]  ; msl
        z[n]=height      
        ;  distance      
        d_hor=(float(k)*rscale[i])*cos(elangle[i]*(!pi/180.))
        distance[n]=d_hor      
        x[n]=d_hor*sin(azimuth[n]*!pi/180.)
        y[n]=d_hor*cos(azimuth[n]*!pi/180.)    
        ;  dbz      
        xstr='dbz[n]=ref'+dnum+'[k,j]
        result=execute(xstr)
        ;  Time
        julian_day[n]=sjulian_day[i]+double(j)*delta_t[i]    
      
        n=n+1l
      endfor  ;end of loop through ranges k
    endfor  ;end of loop through azimuths j

    ;  plot these arrays
    do_plot2='no'
    if do_plot2 eq 'yes' then begin
      pxdim=700 & pydim=600
      ;p1idx=0 & p2idx=2000
      ;  Position the plots
      xl=0.10 & xr=0.95
      yb=0.05 & yt=0.90
      sx=0.14
      sy=0.03
      numplots_x=1
      numplots_y=3
      position_plots,xl,xr,yb,yt,sx,sy,numplots_x,numplots_y,pos
      ;azimuth[0:-1:50],degree[0:-1:50],dbz[0:-1:50]
      p0=plot(azimuth[0:500],/buffer,position=pos[0,*],dimensions=[pxdim,pydim],ytitle='azimuth')
      p1=plot(degree[0:500],/current,position=pos[1,*],ytitle='degree')
      ;p1=plot(elevation[0:-1:50],/current,position=pos[2,*],ytitle='elevation')
      p1=plot(dbz[0:500],/current,position=pos[2,*],ytitle='dbz',linestyle=6,symbol='o',sym_size=0.3,$
        yrange=[-50,50])
      ;p1=plot(distance[p1idx:p2idx],/current,position=pos[4,*],ytitle='distance')  
      ;p1=plot(x[p1idx:p2idx],/current,position=pos[5,*],ytitle='x')
      ;p1=plot(y[p1idx:p2idx],/current,position=pos[6,*],ytitle='y')
      ;p1=plot(z[p1idx:p2idx],/current,position=pos[7,*],ytitle='z')
      ddx=0.1
      ddy=0.02
      t1=text(xl+0.*ddx,yt+ddy,'scan '+dnum)
      t1=text(xl+1.*ddx,yt+ddy,'nrays '+strcompress(nrays[i],/remove_all)) 
      t1=text(xl+2.*ddx,yt+ddy,'nbins '+strcompress(nbins[i],/remove_all))    
      t1=text(xl+3.*ddx,yt+ddy,'ray spacing '+strcompress(ray_spacing[i],/remove_all)) 
      t1=text(xl+0.*ddx,yt+2.*ddy,'elevation '+strcompress(elangle[i],/remove_all))    
      t1=text(xl+2.*ddx,yt+2.*ddy,'rscale '+strcompress(rscale[i],/remove_all) )     
      t1=text(xl+4.*ddx,yt+2.*ddy,'rstart '+strcompress(rstart[i],/remove_all) )   
      p0.save,'coor_conv.'+dnum+'.aza.png',height=pydim
    endif  ;end of do plot

    ;  finished a scan.  Rename arrays
    xstr='azimuth'+dnum+'=azimuth'
    result=execute(xstr)
    xstr='degree'+dnum+'=degree'
    result=execute(xstr)
    xstr='elevation'+dnum+'=elevation'
    result=execute(xstr)
    xstr='range'+dnum+'=range'
    result=execute(xstr)
    xstr='x'+dnum+'=x'
    result=execute(xstr)
    xstr='y'+dnum+'=y'
    result=execute(xstr)
    xstr='z'+dnum+'=z'
    result=execute(xstr)
    xstr='dbz'+dnum+'=dbz'
    result=execute(xstr)
    xstr='distance'+dnum+'=distance'
    result=execute(xstr)
    xstr='julian_day'+dnum+'=julian_day' 
    result=execute(xstr)
  endfor  ;end of loop through scans i
  ;print,'finished coordinate conversion'

  do_plot3='no'
  if do_plot3 eq 'yes' then begin
  ;  Plot and histogram dbz data straight from the file
  ;  Set up the plot
  pxdim=900 & pydim=900
  ;  Position the plots
  xl=0.10 & xr=0.88
  yb=0.06 & yt=0.95
  sx=0.09
  sy=0.07
  numplots_x=2
  numplots_y=3
  position_plots,xl,xr,yb,yt,sx,sy,numplots_x,numplots_y,pos
  ;Clear plot space
  pnum=0
  p0=plot([0,1],[0,1],position=pos[pnum,*],/buffer,dimensions=[pxdim,pydim],$
    axis_style=4,/nodata)
  d=p0.convertcoord([pos[pnum,0],pos[pnum,2]],[pos[pnum,1],pos[pnum,3]],/normal,/to_device)
  isx=(d[0,1,0]-d[0,0,0])
  isy=(d[1,1,0]-d[1,0,0])

  dbz1r=reform(dbz1,s1_orig[0],s1_orig[1])  ;maybe this is nbins,nrays
  dbz2r=reform(dbz2,s2_orig[0],s2_orig[1])

  pnum=0
  data_var=dbz1r
  result=where(data_var ne -9999,count)
  start_bin=min(data_var[result])
  end_bin=max(data_var[result])
  dbin=(end_bin-start_bin)/100.
  hist_generic,data_var,start_bin,end_bin,dbin,bins1,data_freq1,data_counts1
  p1=plot(bins1,data_counts1,/buffer,color='black',/hist,position=pos[pnum,*],$
    ytitle='counts',xtitle='dbz1r1')
  t1=text(pos[pnum,0],pos[pnum,3],'total cloudy='+string(total(data_counts1)))

  pnum=2
  data_image=bytscl(data_var,top=top_color,min=start_bin,max=end_bin)
  data_var=data_image
  result=where(data_var ne -9999,count)
  start_bin=3;min(data_var[result])  ;eliminate 0
  end_bin=max(data_var[result])
  dbin=(end_bin-start_bin)/100.
  hist_generic,data_var,start_bin,end_bin,dbin,bins1,data_freq1,data_counts1
  p1=plot(bins1,data_freq1,/current,color='black',/hist,position=pos[pnum,*],$
    ytitle='frequency',xtitle='bytscaled data')

  pnum=4
  result=where(dbz1r eq -9999,count)
  data_image[result]=253
  result=where(dbz1r ne -9999,count)
  data_image[result]=70
  p1=image(data_image,/current,image_dimensions=[isx,isy],position=pos[pnum,*],rgb_table=mytable)
  c1=contour(data_var,ytitle='Y (km)',xtitle='X (km)',$
    xstyle=1,ystyle=1,font_size=10,/nodata,/current,position=pos[pnum,*],$
    irregular=0,xtickdir=1,ytickdir=1,font_name='Helvetica',$
    axis_style=2)

  pnum=1
  data_var=dbz2r
  result=where(data_var ne -9999,count)
  start_bin=min(data_var[result])
  end_bin=max(data_var[result])
  dbin=(end_bin-start_bin)/100.
  hist_generic,data_var,start_bin,end_bin,dbin,bins1,data_freq1,data_counts1
  p1=plot(bins1,data_counts1,/current,color='black',/hist,position=pos[pnum,*],$
    ytitle='counts',xtitle='dbz2r')

  pnum=3
  data_image=bytscl(data_var,top=top_color,min=start_bin,max=end_bin)
  data_var=data_image
  result=where(data_var ne -9999,count)
  start_bin=3;min(data_var[result])
  end_bin=max(data_var[result])
  dbin=(end_bin-start_bin)/100.
  hist_generic,data_var,start_bin,end_bin,dbin,bins1,data_freq1,data_counts1
  p1=plot(bins1,data_freq1,/current,color='black',/hist,position=pos[pnum,*],$
    ytitle='frequency',xtitle='bytscaled data')

  pnum=5
  result=where(dbz2r eq -9999,count)
  data_image[result]=253
  result=where(dbz2r ne -9999,count)
  data_image[result]=70
  p1=image(data_image,/current,image_dimensions=[isx,isy],position=pos[pnum,*],rgb_table=mytable)
  c1=contour(data_var,ytitle='Y (km)',xtitle='X (km)',$
    xstyle=1,ystyle=1,font_size=10,/nodata,/current,position=pos[pnum,*],$
    irregular=0,xtickdir=1,ytickdir=1,font_name='Helvetica',$
    axis_style=2)
  p1.save,'dbz_cconv_'+file_tag+'.png',height=pydim
  endif
  
  ;*******************
  ;  Now plot the scans
  ;*******************

  ;r=range
  ;theta=azm
  ;f=ref
  ;  Loop throught the elevations to plot
for i=0,nelev-1 do begin

  dnum=strcompress(i+1,/remove_all)
  print,dnum,'dnum'

  ;  Rename arrays
  xstr='azimuth=azimuth'+dnum
  result=execute(xstr)
  xstr='degree=degree'+dnum
  result=execute(xstr)    
  xstr='elevation=elevation'+dnum
  result=execute(xstr)
  xstr='range=range'+dnum
  result=execute(xstr)
  xstr='x=x'+dnum
  result=execute(xstr)
  xstr='y=y'+dnum
  result=execute(xstr)
  xstr='z=z'+dnum
  result=execute(xstr)
  xstr='dbz=dbz'+dnum
  result=execute(xstr)
  xstr='distance=distance'+dnum
  result=execute(xstr)
  xstr='julian_day=julian_day'+dnum
  result=execute(xstr)
 
  ;  exclude -999 points for cloud area griding
  result=where(dbz ne -9999,count_good)
  azimuth_good=azimuth[result]  ;true azimuth
  degree_good=degree[result]    ;azimuth converted to degrees
  elevation_good=elevation[result]
  range_good=range[result]
  x_good=x[result]
  y_good=y[result]
  z_good=z[result]
  dbz_good=dbz[result]
  distance_good=distance[result]
  julian_day_good=julian_day[result]

  ;  Convert from r,theta to x,y
  ;  Conversion for mask - have to grid twice to find the non-scanned area
  grid_input,range,degree,dbz,x_new,y_new,dbz_new,/polar,/degrees
  ;  Conversion for cloud area
  grid_input,range_good,degree_good,dbz_good,x_good_new,y_good_new,dbz_good_new,/polar,/degrees  
    
  ;  decrease size of arrays so that it doesn't take too long to plot
  num=10
  x_new=x_new[0:-1:num]
  y_new=y_new[0:-1:num]
  dbz_new=dbz_new[0:-1:num]
  julian_day_new=julian_day[0:-1:num]
  num=2
  x_good_new=x_good_new[0:-1:num]
  y_good_new=y_good_new[0:-1:num]
  dbz_good_new=dbz_good_new[0:-1:num]

  ;  Define x,y grid to put the data on
  dist_max=120000.0
  x_out=-1.0*dist_max
  while max(x_out) lt dist_max do x_out=[x_out,max(x_out)+100.]
  y_out=-1.0*dist_max
  while max(y_out) lt dist_max do y_out=[y_out,max(y_out)+100.]

;  radial_basis_function does not work
  print,'start gridding'
  qhull,x_new,y_new,triangle,/delaunay
  dbz_mask=griddata(x_new,y_new,dbz_new,xout=x_out,yout=y_out,/grid,triangles=triangle,$
  ;/kriging,min_points=8,sectors=8,empty_sectors=4,missing=-8888)
  ;/inverse_distance,min_points=10,sectors=8,empty_sectors=3,missing=-8888) 
  /kriging,min_points=10,sectors=8,empty_sectors=5,missing=-8888)  
  qhull,x_good_new,y_good_new,triangle,/delaunay
  ;dbz_grid=griddata(x_good_new,y_good_new,dbz_good_new,xout=x_out,yout=y_out,/grid,triangles=triangle,$
  ;/kriging,min_points=10,sectors=8,empty_sectors=4,missing=-9999)
  dbz_grid=griddata(x_good_new,y_good_new,dbz_good_new,xout=x_out,yout=y_out,/grid,triangles=triangle,$
  ;/nearest_neighbor,missing=-9999)
  /kriging,min_points=10,sectors=8,empty_sectors=5,missing=-9999)  
  ;/kriging,min_points=8,sectors=8,empty_sectors=4,missing=-9999)
  ;/inverse_distance,min_points=10,sectors=8,empty_sectors=3,missing=-9999) 
    
  print,'finish griddata'  

;  pass out these arrays for plotting
;dbz_grid
;x_out/1000.0,y_out/1000.0
;nbins

do_plot4='no'
if do_plot4 eq 'yes' then begin  
  ;*** Color table ***
  ;  Top is the last color to scale 256 colors, 0-255
  top_color=252
  ;  Colortable  0-252  253=white
  mytable=colortable(39,ncolors=254)
  ;254=hot pink               ;gray=255
  mytable=[mytable,transpose([238,18,137]),transpose([150,150,150])]
  cbtable=mytable[0:top_color,*]

  ;  Quick plot of dbz_grid
  pxdim=700 & pydim=600
  ;  Position the plots
  xl=0.15 & xr=0.80
  yb=0.15 & yt=0.90
  sx=0.14
  sy=0.03
  numplots_x=1
  numplots_y=1
  position_plots,xl,xr,yb,yt,sx,sy,numplots_x,numplots_y,pos
  ;  Colorbar position
  cbpos=pos
  cbpos[*,0]=pos[*,2]+0.14
  cbpos[*,2]=cbpos[*,0]+0.02
  ;  set plot area
  pnum=0
  p0=plot([0,1],[0,1],/buffer,/nodata,axis_style=4,dimensions=[pxdim,pydim],position=pos[pnum,*])
  d=p0.convertcoord([pos[pnum,0],pos[pnum,2]],[pos[pnum,1],pos[pnum,3]],/normal,/to_device)
  isx=(d[0,1,0]-d[0,0,0])
  isy=(d[1,1,0]-d[1,0,0])
  px0=d[0,0,0]
  py0=d[1,0,0]
  ;  Plot reflectivity
  pnum=0
  
  ;data_var=dbz_mask
  ;dmin=min(data_var)
  ;dmax=max(data_var)
  ;data_image=bytscl(data_var,min=dmin,max=dmax)
  
  data_var=dbz_grid
  ;result=where(data_var ne -9999,count)
  dmin=-20.0;min(data_var[result])
  dmax=50.0;max(data_var[result])
  data_image=bytscl(data_var,top=top_color,min=dmin,max=dmax)
  ;result=where(data_var eq -999.,count)
  result=where(data_var eq -9999,count)
  if count gt 0 then data_image[result]=253
  ;result=where(dbz_mask eq -8888,count)
  ;if count gt 0 then data_image[result]=255  
  p1=image(data_image,/current,image_dimensions=[isx,isy],position=pos[pnum,*],rgb_table=mytable)
  c1=contour(data_var,x_out/1000.0,y_out/1000.0,ytitle='Y (km)',xtitle='X (km)',$
    xstyle=1,ystyle=1,font_size=14,/nodata,/current,position=pos[pnum,*],$
    irregular=0,xtickdir=1,ytickdir=1,font_name='Helvetica',$
    axis_style=2)
  ;  plot range rings
  for r=0,nbins[i]-1,150 do begin
    ;print,'range ring',r
    ;  radius in km
    range_ring=(rstart[i]+float(r)*rscale[i])/1000.0
    result=circle(0,0,range_ring)
    p2=plot(result,/current,/overplot)
    ;stop
  endfor
  c2=colorbar(orientation=1,rgb_table=cbtable,range=[dmin,dmax],title='Z (dBZ)',$
    position=reform(cbpos[pnum,*]),font_size=14,font_name='Helvetica',tickdir=1,border_on=1)
  caldat,julian_day1[0],mm,dd,yy,hh,mi,ss
  date_time=string(yy,format='(I04)')+string(mm,format='(I02)')+string(dd,format='(I02)')+$
    ' '+string(hh,format='(I02)')+':'+string(mi,format='(I02)')+':'+string(ss,format='(I02)')
  scan_str=string(elangle[i],format='(F05.2)')
  t1=text(xl,0.95,'el='+scan_str,font_size=14)
  t1=text(0.5,0.95,date_time,font_size=14)
  imagename='dbz_grid_'+file_tag+'.'+dnum+'.png'
  p0.save,imagename,height=pydim
  ;p0.save,'dbz_grid.el'+scan_str+'.'+startdate[i]+'.'+starttime[i]+'.mp08.png',height=pydim
  print,imagename
endif  ; do plot
  
endfor  ;  end of loop through scans to plot
;stop
endfor  ;end of loop through files
endif  ;scan within 2 hours
end