;********************************************
;  This plots a plan view (horizontal slice on a map)
;  Constant Latitude slice, constant longitude slice, and line segment slice
;  4 plot.
;**************************************************

;*****************************
;  Program to plot cband radar
;*****************************
pro cband_3panel,filenames
;filenames='/Volumes/mace-group5/field_programs/investigator/in2018_v01/cband/20180209/9776HUB-PPIVol-20180209-235943-0000.cfradial.pyart.81.401.nc'
;  Choose the horizonal surface height
surf=2.5  ;km
top_height=14 ;km

if filenames eq !NULL then begin
  ;  Choose radar site, date and time 
  ;program_name='nosmn'
  ;date_str='20200328' & time_str='*';'085004'

  ;program_name='normet'
  ;date_str='20200327' & time_str='*';'085004'

  program_name='in2018_v01'
  date_str='20180131' & time_str='*' 

  ;program_name='olympex'
  ;date_str='20151204' & time_str='163710' 

  ;  Get year month day strings
  syy=strmid(date_str,0,4)
  smm=strmid(date_str,4,2)
  sdd=strmid(date_str,6,2)

  shh=strmid(time_str,0,2)
  smi=strmid(time_str,2,2)
  sss=strmid(time_str,4,2)
  ;date_str=syy+smm+sdd+shh+smi+sss
  ;ymd_str=syy+smm+sdd

  ;  chpc or imac
  ;path_prefix='/Volumes/'  
  path_prefix='/uufs/chpc.utah.edu/common/home/'
  
  ;  directory
  if program_name eq 'nosmn' then begin  ;comble
    filedir=path_prefix+'mace-group5/field_programs/comble/cband_somna/'
    filestr='nosmn.ss.vol-dbz*noblock.'+syy+smm+sdd+'T'+time_str+'*pyart.41.201.nc'
    radar_name='cband'
  endif else if program_name eq 'in2018_v01' then begin
    ;filedir=path_prefix+'mace-group4/Capricorn2/CBand/'
    ;filedir=path_prefix+'mace-group4/Capricorn2/CBand/cfradial/'
    filedir=path_prefix+'mace-group5/field_programs/investigator/in2018_v01/cband/'+date_str+'/'
    filestr='9776HUB-PPIVol-'+syy+smm+sdd+'-'+time_str+'*pyart.81.401.nc'
    radar_name='cband'
  endif else if program_name eq 'in2019_v01' then begin
    filedir=path_prefix+'mace-group5/field_programs/investigator/in2019_v01_gpm_passes/out_vol/'
    filestr='9776HUB-PPIVol-'+syy+smm+sdd+'-'+time_str+'*pyart*.nc'
    radar_name='cband'
  endif else if program_name eq 'olympex' then begin
    filedir=path_prefix+'mace-group6/field_programs/olympex/sband/'+date_str+'/ppi/'
    filestr='olympex_NPOL1_'+date_str+'_'+time_str+'*pyart*.nc'
    radar_name='sband'
    ;  Get the flight path
    nav_date_str=date_str;'20151213'
    read_nav_olympex,path_prefix,nav_date_str,julian_day_plane,temp,tas,theta,pres,ppmvw,lat_plane,lon_plane,$
      alt_plane,wwind,uwind,vwind,cwc_kng,twc_nev,lwc_nev,iwc_nev,tot_s,cwc_s,rh,cwc_csi,tot_c,tot_p
    ;  Get the dc8 track - do this after reading in apr3
    ;read_apr3_olympex,path_prefix,date_str,lat_dc8,lon_dc8,alt_dc8,julian_day_dc8,$
    ;  kuz_dc8,kaz_dc8,alt3d_dc8
  endif else if program_name eq 'normet' then begin
    filedir=path_prefix+'mace-group5/field_programs/comble/cband/'+date_str+'/'
    filestr=date_str+time_str+'*pyart*.nc'
    radar_name='cband'
  endif
  
  ;  file name
  filenames=file_search(filedir+filestr,count=numfiles)
  print,'numfiles',numfiles
  
endif else begin
 
  ;filenames='/uufs/chpc.utah.edu/common/home/mace-group6/field_programs/olympex/sband/20151130/olympex_NPOL1_20151130*.pyart.81.401.nc"

  filedir=file_dirname(filenames)+'/'

  ;  Get year month day strings
  if strmid(file_basename(filenames),0,7) eq 'olympex' then begin
    parts=strsplit(file_basename(filenames),'_',/extract)
  endif else begin
    parts=strsplit(file_basename(filenames),'-',/extract)
  endelse
  syy=strmid(parts[2],0,4)
  smm=strmid(parts[2],4,2)
  sdd=strmid(parts[2],6,2)

  shh=strmid(parts[3],0,2)
  smi=strmid(parts[3],2,2)
  sss=strmid(parts[3],4,2)
  
  date_str=syy+smm+sdd
  if syy eq '2015' or syy eq '2016' then radar_name='sband'
  if syy eq '2018' then radar_name='cband'

  numfiles=1
endelse

;  Loop through the cband files
for idx=0,numfiles-1 do begin

  filename=filenames[idx]
  print,filename
  
  if program_name eq 'olympex' then begin
    parts=strsplit(file_basename(filename),'_',/extract)
    shh=strmid(parts[3],0,2)
    smi=strmid(parts[3],2,2)
    sss=strmid(parts[3],4,2)
    time_str=shh+smi+sss
  endif else if program_name eq 'in2019_v01' then begin
    parts=strsplit(file_basename(filename),'-',/extract)
    shh=strmid(parts[3],0,2)
    smi=strmid(parts[3],2,2)
    sss=strmid(parts[3],4,2)
    time_str=shh+smi+sss   
  endif else if program_name eq 'nosmn' then begin
    parts=strsplit(file_basename(filename),'.',/extract)
    shh=strmid(parts[4],9,2)
    smi=strmid(parts[4],11,2)
    sss=strmid(parts[4],13,2)
    time_str=shh+smi+sss
  endif else if program_name eq 'normet' then begin
    parts=strsplit(file_basename(filename),'.',/extract)
    shh=strmid(parts[0],8,2)
    smi=strmid(parts[0],10,2)
    sss=strmid(parts[0],12,2)
    time_str=shh+smi+sss
  endif
  
  print,'start reading data'  
  if program_name eq 'nosmn' then begin
    read_nosmn_cband,filename,lat,lon,point_z,rho,kdp,zdr,dbz,$
      point_x,point_y,base_time,time_offset
  endif else begin
  fid=ncdf_open(filename)
  vid=ncdf_varid(fid,'base_time') & ncdf_varget,fid,vid,base_time
  vid=ncdf_varid(fid,'time_offset') & ncdf_varget,fid,vid,time_offset  
  ;vid=ncdf_varid(fid,'radar_time') & ncdf_varget,fid,vid,radar_time 
  ;ncdf_attget,fid,vid,'units',astr & astr=string(astr)
  vid=ncdf_varid(fid,'radar_latitude') & ncdf_varget,fid,vid,lat_radar
  vid=ncdf_varid(fid,'radar_longitude') & ncdf_varget,fid,vid,lon_radar 
  vid=ncdf_varid(fid,'point_latitude') & ncdf_varget,fid,vid,lat
  vid=ncdf_varid(fid,'point_longitude') & ncdf_varget,fid,vid,lon      
  vid=ncdf_varid(fid,'point_altitude') & ncdf_varget,fid,vid,alt   
  vid=ncdf_varid(fid,'point_x') & ncdf_varget,fid,vid,point_x
  vid=ncdf_varid(fid,'point_y') & ncdf_varget,fid,vid,point_y
  vid=ncdf_varid(fid,'point_z') & ncdf_varget,fid,vid,point_z
  ;vid=ncdf_varid(fid,'cross_correlation_ratio') & ncdf_varget,fid,vid,rho
  vid=ncdf_varid(fid,'specific_differential_phase')
  if vid eq -1 then vid=ncdf_varid(fid,'corrected_specific_differential_phase')
  if vid ne -1 then begin
    ncdf_varget,fid,vid,kdp
  endif else begin
    kdp=!null
  endelse
  vid=ncdf_varid(fid,'corrected_differential_reflectivity') 
  if vid ne -1 then begin
    ncdf_varget,fid,vid,zdr
  endif else begin
    zdr=!null
  endelse
  vid=ncdf_varid(fid,'corrected_reflectivity')  
  if vid eq -1 then vid=ncdf_varid(fid,'reflectivity')
  if vid eq -1 then vid=ncdf_varid(fid,'DBZH_gmm')
  ncdf_varget,fid,vid,dbz
  ncdf_close,fid
  print,'end reading data'
  endelse
  
  ;  convert point_x,y,z from m to km
  point_x=point_x/1000.0
  point_y=point_y/1000.0
  point_z=point_z/1000.0
    
  ;  convert time to julian_day
  ;time_offset=0L
  day1=julday(1,1,1970,0,0,0)
  sec_per_day=long(24L*60L*60L)
  julian_day=double(day1+((base_time+time_offset)/sec_per_day))    
  caldat,julian_day,mm,dd,yy,hh,mi,ss
  print,yy,mm,dd,hh,mi,ss
  
  ;  Calculate colongitude
  colon=lon
  result=where(lon lt 0,count)
  if count gt 0 then colon[result]=colon[result]+360.0

  ;  Zoom in on radar, dist_km amount (max distance is 123)
  dist_km=150.0  ;75
  x1=point_x[*,0,0]
  xr=where(x1 ge -1*dist_km and x1 le dist_km,xc)
  ;xr=where(x1 lt 0,xc)
  y1=point_y[0,*,0]
  yr=where(y1 ge -1*dist_km and y1 le dist_km,yc)    
  
  dbz=dbz[xr,*,*]
  lat=lat[xr,*,*]
  lon=lon[xr,*,*]
  dbz=dbz[*,yr,*]
  lat=lat[*,yr,*]
  lon=lon[*,yr,*]    
  
  if kdp ne !NULL then begin
    kdp=kdp[xr,*,*]
    kdp=kdp[*,yr,*]
  endif
  if zdr ne !NULL then begin
    zdr=zdr[xr,*,*]
    zdr=zdr[*,yr,*]
  endif
  
  ;  Height array
  ;height_cband=alt[0,0,*]/1000.0  ;msl
  height_cband=point_z[0,0,*]  ;agl
  
  ;  Index of surface
  r=where(height_cband ge surf,c)
  sidx=r[0]
  
  ;*************************************
  ;  Find flight track at this time
  ;*************************************
  if julian_day_plane ne !NULL then begin
    xtime=5D/(60D*24D)
    r=where(julian_day_plane ge julian_day-xtime and julian_day_plane le julian_day+xtime and $
          lat_plane ge min(lat) and lat_plane le max(lat) and $
          lon_plane ge min(lon) and lon_plane le max(lon) ,c)
    if c gt 1 then begin
      sub_lat_plane=lat_plane[r]
      sub_lon_plane=lon_plane[r]
      sub_alt_plane=alt_plane[r]/1000.0
      print,'found citation track'
    endif else begin
      print,'no citation track'
    endelse
  endif
  
  if syy eq '2015' or syy eq '2016' then begin
  ;if julian_day_dc8 ne !NULL then begin
    xtime=5D/(60D*24D)
    sjulian_day=julian_day-xtime
    ejulian_day=julian_day+xtime
    read_apr3_olympex,path_prefix,date_str,sjulian_day,ejulian_day,$
      lat_dc8,lon_dc8,alt_dc8,julian_day_dc8,kuz_dc8,kaz_dc8,alt3d_dc8
    r=where(julian_day_dc8 ge julian_day-xtime and julian_day_dc8 le julian_day+xtime and $
      lat_dc8 ge min(lat) and lat_dc8 le max(lat) and $
      lon_dc8 ge min(lon) and lon_dc8 le max(lon),numtimes_dc8)
    if numtimes_dc8 gt 1 then begin
      sub_lat_dc8=lat_dc8[r]
      sub_lon_dc8=lon_dc8[r]
      sub_alt_dc8=alt_dc8[r]/1000.0
      sub_julian_day_dc8=julian_day_dc8[r]
      sub_kuz_dc8=kuz_dc8[r,*]
      sub_kaz_dc8=kaz_dc8[r,*]
      sub_alt3d_dc8=alt3d_dc8[r,*]
      print,'found dc8 track'
      max_height=14000 ;m
      dh=30.0 ;m resolution
      numheight_grid=ceil(max_height/dh)+1
      height_grid=findgen(numheight_grid)*dh
      ;  Put the data on the regular grid
      print,'start grid'
      sub_kuz_grid=make_array(numtimes_dc8,numheight_grid,value=-9999,/float)
      sub_kaz_grid=make_array(numtimes_dc8,numheight_grid,value=-9999,/float)
      for i=0,numtimes_dc8-1 do begin
        for j=0,numheight_grid-1 do begin
          ku1=sub_kuz_dc8[i,*]
          ka1=sub_kaz_dc8[i,*]
          height1=sub_alt3d_dc8[i,*]
          r=where(height1 ge height_grid[j]-dh/2.0 and height1 lt height_grid[j]+dh/2.0,c)
          if c gt 0 then begin
            sub_kuz_grid[i,j]=ku1[r[0]]
            sub_kaz_grid[i,j]=ka1[r[0]]
          endif
        endfor
      endfor
      print,'end grid'
    endif else begin
      print,'no dc8 track'
    endelse
  endif

  ;  Filetag
  ;file_tag=strmid(file_basename(filename),16,19)   
  ;  Filetag
  if program_name eq 'olympex' then file_tag=strmid(file_basename(filename),0,35)
  if program_name eq 'in2018_v01' then file_tag=strmid(file_basename(filename),0,44)
  if program_name eq 'nosmn' then file_tag=strmid(file_basename(filename),0,75)
  if program_name eq 'normet' then file_tag=strmid(file_basename(filename),0,16)   
  
  ;  Set up plot size
  pxdim=1000 & pydim=1000
  ;  Position the plots
  xl=0.1 & xr=0.90
  yb=0.1 & yt=0.90
  sx=0.08
  sy=0.15
  if sub_lat_dc8 ne !NULL then begin
    numplots_x=2
    numplots_y=3
  endif else begin
    numplots_x=2
    numplots_y=2
  endelse
  position_plots,xl,xr,yb,yt,sx,sy,numplots_x,numplots_y,pos
  ;  Colorbar position
  cbpos=pos
  cbpos[*,0]=pos[*,2]+0.08  ;vertical
  cbpos[*,2]=cbpos[*,0]+0.01
  ;cbpos[*,1]=cbpos[*,3]+0.005 ;horizontal
  ;cbpos[*,3]=cbpos[*,1]+0.01
      
  ;  Top is the last color to scale 256 colors, 0-255
  top_color=252
  ;  Colortable  0-252  253=white
  ;mytable=colortable(39,ncolors=254)
  mytable=colortable(43,ncolors=253)
                   ;254=hot pink               ;gray=255
  mytable=[mytable,transpose([255,255,255]),transpose([238,18,137]),transpose([230,230,230])]
  mycbtable=mytable[0:top_color,*]
 
  ;  Top is the last color to scale 256 colors, 0-255
  top_color_dc8=252
  ;  Colortable  0-252  253=white
  ;mytable=colortable(39,ncolors=254)
  ;mytable=colortable(43,ncolors=253)
  mytable_dc8=colortable(33,ncolors=253)
  ;254=hot pink               ;gray=255
  mytable_dc8=[mytable_dc8,transpose([255,255,255]),transpose([238,18,137]),transpose([230,230,230])]
  mycbtable_dc8=mytable_dc8[0:top_color_dc8,*]

    
  ;top_color=254
  ;mytable=colortable(39,ncolors=255)
  ;mytable[0,*]=[250,250,250];  first color =gray
  ;mycbtable=mytable[1:-1,*]    

  ;  Strings
  dx=0.01 & dy=0.01
  fs1=12

  ;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])

  ;  radar min and max
  ;dmin_cband=-17 & dmax_cband=31
  ;dmin_cband=-5 & dmax_cband=25
  
  ;r=where(kdp ne -9999,c)
  ;dmin_cband=min(kdp[r])
  ;dmax_cband=max(kdp[r])
  
  ;r=where(reflect ne -9999,c)
  ;dmin_cband=min(reflect[r])
  ;dmax_cband=max(reflect[r])
  
  ;*** choose variable to plot***
  plot_var=dbz
  var_name='dbz'
  ;plot_var=zdr
  ;var_name='zdr'
  ;plot_var=kdp
  ;var_name='kdp'
  
  if var_name eq 'dbz' then begin
    dmin_cband=5
    dmax_cband=50
  endif else if var_name eq 'zdr' then begin
    r=where(plot_var ne -9999,c)
    dmin_cband=-2.0;min(plot_var[r])
    dmax_cband=2.0;max(plot_var[r])  
  endif else if var_name eq 'kdp' then begin
    r=where(plot_var ne -9999,c)
    dmin_cband=-0.5;min(plot_var[r])
    dmax_cband=1.0;max(plot_var[r])  
  endif
  
  ;***************************
  ;  Plot the map in position 1
  ;***************************
  ;  Choose a surface
  ;sidx=4  ;chosen at top
  data_var=reform(plot_var[*,*,sidx])
  lat_var=reform(lat[*,*,sidx])
  lon_var=reform(lon[*,*,sidx])
  llat=min(lat) & ulat=max(lat)
  llon=min(lon) & rlon=max(lon)
  data_image=bytscl(data_var,top=top_color,min=dmin_cband,max=dmax_cband)
  r=where(data_var eq -9999,c)
  if c gt 0 then data_image[r]=253
  r=where(data_var eq 0,c)
  if c gt 0 then data_image[r]=255
  p1=image(data_image,lon_var,lat_var,limit=[llat,llon,ulat,rlon],$
    /current,dimensions=[pxdim,pydim],grid_units='degrees',$
    ;  ;center_longitude=(rcolon_ovp-lcolon_ovp)/2.0,
    rgb_table=mytable,$
    label_format='MapGrid_Labels',/box_axes,label_position=0,$
    label_show=1,irregular=0,$
    map_projection='Mercator',position=pos[pnum,*],font_size=fs1)
  p2=mapcontinents(color='black',/hires)
  ;p0.mapgrid.box_axes=1
  ;p0.mapgrid.box_color='black'
  ;0.mapgrid.color='black'
  ;p0.mapgrid.label_color='black'
  p1.mapgrid.linestyle='dot'
  ;p0.mapgrid.label_position=0    
  if sub_lon_plane ne !NULL then p7=plot(sub_lon_plane,sub_lat_plane,color='blue',thick=2,/overplot)
  if sub_lon_dc8 ne !NULL then p7=plot(sub_lon_dc8,sub_lat_dc8,color='black',thick=2,/overplot)  
  t1=text(pos[pnum,0]-5*dx,pos[pnum,3]+2*dy,string(height_cband[0,0,sidx],format='(F4.2)')+'km AGL',$
    font_size=fs1)

;**************************
;  Plot the longitude constant slice
;**************************
pnum=1
lat1=reform(lat[0,*,0])
lon1=reform(lon[*,0,0])
;  Print along constant longitude line
;xlon=138
xlon=((rlon-llon)/3.0)+llon
closest=min(abs(xlon-lon1),idx_lon)
data_var=reform(plot_var[idx_lon,*,*])
;  lat[0,*,0]=all lats   lat[*,0,0]=constant lat
lat_var=reform(lat[idx_lon,*,*])
;  lon[0,*,0]=constant lon   lon[*,0,0]=all lon
lon_var=reform(lon[idx_lon,*,*])

;p0=image(data_var,$
;  /current,dimensions=[pxdim,pydim],$
;  ;  ;center_longitude=(rcolon_ovp-lcolon_ovp)/2.0,
;  rgb_table=mytable,min_value=dmin_cband,max_value=dmax_cband,$
;  irregular=0,$
;  position=pos[pnum,*],font_size=fs1)

;  Choose height range to plot
hidx=where(height_cband lt top_height)
data_image=bytscl(data_var,top=top_color,min=dmin_cband,max=dmax_cband)
r=where(data_var eq -9999,c)
if c gt 0 then data_image[r]=253
r=where(data_var eq 0,c)
if c gt 0 then data_image[r]=255
p4=image(data_image[*,hidx],/current,image_dimensions=[isx,isy],$
  ;min_value=dmin_cband,max_value=dmax_cband,$
  position=pos[pnum,*],rgb_table=mytable)
  
c1=contour(data_var[*,hidx],lat1,height_cband[hidx],$
  ytitle='height (km)',xtitle='latitude',irregular=0,xtickdir=1,ytickdir=1,$
  xstyle=1,ystyle=1,/nodata,/current,position=pos[pnum,*],font_size=fs1,$
  axis_style=2,$
  yrange=[height_cband[hidx[0]],height_cband[hidx[-1]]]);,xrange=[lat_var[0,0],lat_var[-1,0]])
p6=plot(lat1,make_array(n_elements(lat1),value=height_cband[sidx]),$
        /overplot,color='gray',thick=3,linestyle=2)     
if sub_lat_plane ne !NULL then p7=plot(sub_lat_plane,sub_alt_plane,/overplot,color='blue',thick=2)        
if sub_lat_dc8 ne !NULL then p7=plot(sub_lat_dc8,sub_alt_dc8,/overplot,color='black',thick=2)
             
p5=plot(lon_var[*,sidx],lat1,overplot=p1,color='gray',thick=3,linestyle=2)
t1=text(pos[pnum,0]+dx,pos[pnum,3]+2*dy,'Constant Longitude '+string(xlon,format='(F6.1)'),font_size=fs1)

;**************************
;  Plot the latitude constant slice
;**************************
pnum=2
lat1=reform(lat[0,*,0])
lon1=reform(lon[*,0,0])
;  Print along constant longitude line
;xlat=-65.0
xlat=((ulat-llat)/3.0)+llat
closest=min(abs(xlat-lat1),idx_lat)

data_var=reform(plot_var[*,idx_lat,*])
;  lat[0,*,0]=all lats   lat[*,0,0]=constant lat
lat_var=reform(lat[*,idx_lat,*])
;  lon[0,*,0]=constant lon   lon[*,0,0]=all lon
lon_var=reform(lon[*,idx_lat,*])

;p0=image(data_var,$
;  /current,dimensions=[pxdim,pydim],$
;  ;  ;center_longitude=(rcolon_ovp-lcolon_ovp)/2.0,
;  rgb_table=mytable,min_value=dmin_cband,max_value=dmax_cband,$
;  irregular=0,$
;  position=pos[pnum,*],font_size=fs1)

;  Choose height range to plot
hidx=where(height_cband lt top_height)
data_image=bytscl(data_var,top=top_color,min=dmin_cband,max=dmax_cband)
r=where(data_var eq -9999,c)
if c gt 0 then data_image[r]=253
r=where(data_var eq 0,c)
if c gt 0 then data_image[r]=255
p7=image(data_image[*,hidx],/current,image_dimensions=[isx,isy],$
  ;min_value=dmin_cband,max_value=dmax_cband,$
  position=pos[pnum,*],rgb_table=mytable)

c1=contour(data_var[*,hidx],lon1,height_cband[hidx],$
  ytitle='height (km)',xtitle='longitude',irregular=0,xtickdir=1,ytickdir=1,$
  xstyle=1,ystyle=1,/nodata,/current,position=pos[pnum,*],font_size=fs1,$
  axis_style=2,$
  yrange=[height_cband[hidx[0]],height_cband[hidx[-1]]]);,xrange=[lat_var[0,0],lat_var[-1,0]])
p8=plot(lon1,make_array(n_elements(lon1),value=height_cband[sidx]),$
  /overplot,color='gray',thick=3,linestyle=3)
if sub_lon_plane ne !NULL then p7=plot(sub_lon_plane,sub_alt_plane,/overplot,color='blue',thick=2)  
if sub_lon_dc8 ne !NULL then p7=plot(sub_lon_dc8,sub_alt_dc8,/overplot,color='black',thick=2)

p9=plot(lon1,lat_var[*,sidx],overplot=p1,color='gray',thick=3,linestyle=3)
t1=text(pos[pnum,0]+dx,pos[pnum,3]+2*dy,'Constant Latitude '+string(xlat,format='(F6.1)'),font_size=fs1)

;**************************
;  Plot a line
;**************************
if 1 eq 1 then begin
pnum=3
;lon_a=136.50 & lat_a=-65.50
;lon_b=139.50 & lat_b=-64.50
lon_a=min(lon) & lat_a=min(lat)
lon_b=max(lon) & lat_b=max(lat)

x_a=lon_a
x_b=lon_b
y_a=lat_a
y_b=lat_b

slope_m=(y_a-y_b)/(x_a-x_b)
intercept=y_a-(slope_m*x_a)
;y=(slope_m*x)+intercept
y=(slope_m*lon1)+intercept
r=where(lon1 ge lon_a and lon1 le lon_b and $
        y ge lat_a and y le lat_b,c)
y=y[r]  ;lat
x=lon1[r]

data_var=make_array(n_elements(x),n_elements(height_cband))
for i=0,n_elements(y)-1 do begin
  closest=min(abs(y[i]-lat1),idx_lat)
  closest=min(abs(x[i]-lon1),idx_lon) 
  data_var[i,*]=plot_var[idx_lon,idx_lat,*]
endfor

;  Choose height range to plot
hidx=where(height_cband lt top_height)
data_image=bytscl(data_var,top=top_color,min=dmin_cband,max=dmax_cband)
r=where(data_var eq -9999,c)
if c gt 0 then data_image[r]=253
r=where(data_var eq 0,c)
if c gt 0 then data_image[r]=255
p11=image(data_image[*,hidx],/current,image_dimensions=[isx,isy],$
  ;min_value=dmin_cband,max_value=dmax_cband,$
  position=pos[pnum,*],rgb_table=mytable)

c1=contour(data_var[*,hidx],x,height_cband[hidx],$
  ytitle='height (km)',xtitle='Longitude',irregular=0,xtickdir=1,ytickdir=1,$
  xstyle=1,ystyle=1,/nodata,/current,position=pos[pnum,*],font_size=fs1,$
  axis_style=2,$
  yrange=[height_cband[hidx[0]],height_cband[hidx[-1]]]);,xrange=[lat_var[0,0],lat_var[-1,0]])
xaxis=axis('X',location=-2,title='Latitude',target=c1,$
  coord_transform=[intercept,slope_m],tickfont_size=fs1)

p8=plot(x,make_array(n_elements(x),value=height_cband[sidx]),$
  /overplot,color='gray',thick=3,linestyle=4)

if sub_lon_plane ne !NULL then p7=plot(sub_lon_plane,sub_alt_plane,/overplot,color='blue',thick=2) 
if sub_lon_dc8 ne !NULL then p7=plot(sub_lon_dc8,sub_alt_dc8,/overplot,color='black',thick=2)



p10=plot(x,y,overplot=p1,color='gray',thick=3,linestyle=4)
;t1=text(pos[pnum,0]+dx,pos[pnum,3]+2*dy,'Constant Latitude '+string(xlat,format='(F6.1)'),font_size=fs1)
t1=text(pos[pnum,0]+dx,pos[pnum,3]+2*dy,'line segment',font_size=fs1)

;  Colorbar
cb0=colorbar(rgb_table=mycbtable,range=[dmin_cband,dmax_cband],$
  orientation=0,position=[0.55,0.97,0.9,0.98],font_size=fs1,$
  tickdir=1,border_on=1,title='dBZ')
endif

;t1=text(0.1,0.97,'CBAND '+date_str+' '+shh+':'+smi+':'+sss+' KDP',font_size=fs1)
t1=text(0.1,0.97,radar_name+' '+date_str+' '+shh+':'+smi+':'+sss+' '+var_name,font_size=fs1)
if sub_lon_plane ne !NULL then t1=text(0.1,0.95,'+/- 5min citation track',color='blue',font_size=fs1)
if sub_lon_dc8 ne !NULL then t1=text(0.1,0.95,'+/- 5min dc8 track',color='black',font_size=fs1)

;**********Now plot curtains****
if sub_lat_dc8 ne !NULL then begin
  pnum=4
  dmin_ku=-40.0 & dmax_ku=60.0
  data_var=sub_kuz_grid
  height_var=height_grid
  time_var=sub_julian_day_dc8
  data_image=bytscl(data_var,top=top_color_dc8,min=dmin_ku,max=dmax_ku)
  r=where(data_var le -99.88,c)
  if c gt 0 then data_image[r]=253
  ;r=where(data_var eq 0,c)
  ;if c gt 0 then data_image[r]=255
  p4=image(data_image,/current,image_dimensions=[isx,isy],$
    ;min_value=dmin_cband,max_value=dmax_cband,$
    position=pos[pnum,*],rgb_table=mytable_dc8)
  c1=contour(data_image,time_var,height_var/1000.0,$
    ytitle='height (km)',irregular=0,xtickdir=1,ytickdir=1,$
    xstyle=1,ystyle=1,/nodata,/current,position=pos[pnum,*],font_size=fs1,$
    xtickunits='time',xtickformat='label_date',$
    axis_style=2);,$
    ;yrange=[height_var[0]/1000.0,height_var[-1]/1000.0]);,xrange=[lat_var[0,0],lat_var[-1,0]])
  ;  Colorbar
  cb0=colorbar(rgb_table=mycbtable_dc8,range=[dmin_ku,dmax_ku],$
    orientation=1,position=cbpos[pnum,*],font_size=fs1,$
    tickdir=1,border_on=1,title='dBZ')

  pnum=5
  dmin_ka=-40.0 & dmax_ka=60.0
  data_var=sub_kaz_grid
  height_var=height_grid
  time_var=sub_julian_day_dc8
  data_image=bytscl(data_var,top=top_color_dc8,min=dmin_ka,max=dmax_ka)
  r=where(data_var le -99.88,c)
  if c gt 0 then data_image[r]=253
  ;r=where(data_var eq 0,c)
  ;if c gt 0 then data_image[r]=255
  p4=image(data_image,/current,image_dimensions=[isx,isy],$
    ;min_value=dmin_cband,max_value=dmax_cband,$
    position=pos[pnum,*],rgb_table=mytable_dc8)
  c1=contour(data_image,time_var,height_var/1000.0,$
    ytitle='height (km)',irregular=0,xtickdir=1,ytickdir=1,$
    xstyle=1,ystyle=1,/nodata,/current,position=pos[pnum,*],font_size=fs1,$
    xtickunits='time',xtickformat='label_date',$
    axis_style=2);,$
  ;yrange=[height_var[0]/1000.0,height_var[-1]/1000.0]);,xrange=[lat_var[0,0],lat_var[-1,0]])
  ;  Colorbar
  cb0=colorbar(rgb_table=mycbtable_dc8,range=[dmin_ka,dmax_ka],$
    orientation=1,position=cbpos[pnum,*],font_size=fs1,$
    tickdir=1,border_on=1,title='dBZ')


  
endif




imagename=file_tag+'.'+string(surf,format='(I1)')+'km.'+var_name+'.png'
print,imagename
p0.save,filedir+imagename,height=pydim

endfor

end