;
;  This program makes the quicklooks for mpl data
;

pro plot_mpl,input_flag

;
;  Common to hold data
;
common mpl_data, base_time, time_offset, height,$
	backscat,mpl_base,km_flag
;
;  find date and site depending on input flag
;
input_flag=strcompress(input_flag,/remove_all)
if input_flag eq 'one' then begin
  site=' '
  read,site,prompt='Enter ARM site (nsa,sgp,twpc1,twpc2): '
  site=strcompress(site,/remove_all)
  curdate=' '
  read,curdate,prompt='Enter date (yyyymmdd): '
  curdate=strcompress(curdate,/remove_all)
endif else if input_flag eq 'eos' then begin
  fname=''
  openr,1,'/data/mace4/arm_data/eos/eos_netcdf_file.dat'
  readf,1,fname
  close,1
  fname=strcompress(fname,/remove_all)
  check=strmid(fname,0,29)
  if check eq '/data/mace2/arm_data/nsa/nsaC' then begin
    site='nsa'
    curdate=strmid(fname,69,15)
    openw,12,'/data/mace4/arm_data/eos/eos_page_date.dat'
    printf,12,strmid(fname,69,8)
    close,12
  endif else if check eq '/data/mace3/arm_data/sgp/sgpC' then begin
    site='sgp'
    curdate=strmid(fname,69,15)
    openw,12,'/data/mace4/arm_data/eos/eos_page_date.dat'
    printf,12,strmid(fname,69,8)
    close,12
  endif else if check eq '/data/mace/arm_data/twp/Manus' then begin
    site='twpc1'
    curdate=strmid(fname,77,15)
    openw,12,'/data/mace4/arm_data/eos/eos_page_date.dat'
    printf,12,strmid(fname,77,8)
    close,12
  endif else if check eq '/data/mace/arm_data/twp/Nauru' then begin
    site='twpc2'
    curdate=strmid(fname,77,15)
    openw,12,'/data/mace4/arm_data/eos/eos_page_date.dat'
    printf,12,strmid(fname,77,8)
    close,12
  endif
endif  
;
;  Format the date to process
;
print,'date  ',curdate,'site',site
curdate=strcompress(string(curdate),/remove_all)
date_string_fm_adtg, curdate, date_string
;
;  Get the data
;
get_mpl, curdate,site,datastream,input_flag

;*********************************
;  Sub set the data
;*********************************

;
;  Determine the height subsetting indicies
;
if site eq 'sgp' then eheight=15.0 else $
if site eq 'nsa' then eheight=10.0 else $
if site eq 'twpc1' or site eq 'twpc2' then eheight=20.0
;
;  Enter start height below
;
sheight=0.0
closeheight=height-sheight
closest=min(abs(closeheight),indsh)
sheight=height[indsh]
;
;  Enter end height below
;
;eheight=15.0
closeheight=height-eheight
closest=min(abs(closeheight),indeh)
eheight=height[indeh]
;
;  Subset the arrays for the height bounds
;
backscat=backscat[*,indsh:indeh]
height=height[indsh:indeh]

;***************************************
;  End of height subsetting
;***************************************

;
;  Calculate fractional hour and Julian day
;
ijday_ihrfrac_fm_numsec,base_time,time_offset,hrfrac,jday,yy,mm,dd,hh,mi,ss
;
;  Calculate IDL Julian day for plotting
;
julian_day=make_array(n_elements(jday),/double,value=-9999)
julian_day=julday(mm,dd,yy,hh,mi,ss)

;******************************************************
;  Calculate the background backscatter and subtract it
;  from the backscatter array
;******************************************************

;
;  Find the size of the arrays
;
numtimes=n_elements(time_offset)  ;number of times
numheights=n_elements(height)    ;number of heights
;
;  Make array to hold the background and new signal information
;
background=make_array(numtimes,2,/float,value=0)
signal=make_array(numtimes,numheights,/float,value=0)
;
;  Loop though the shots and calculate background
;
for t=0,numtimes-1 do begin
;
;  One shot
;
  shot=reform(backscat[t,*])
;
;  Calculate the average of the last 50 pts as background
;
  back_stat=moment(shot[numheights-50:numheights-2])
  background[t,0]=back_stat[0]  ;average
  background[t,1]=back_stat[1]  ;variance
;
;  Subtract the background from the shot
;
  signal[t,*]=backscat[t,*]-background[t,0]
endfor  ;end of loop through shots

;
;  Find where background is <= 0 and replace with 
;  the average backgroud value.
;
;if min(background) le 0.0 then begin
;  stat_back=moment(background)
;  result=where(background le 0.0,count)
;  if count gt 0 then background[result]=stat_back[0]
;endif
;
;  Now divide the shots by the shot background
;
;for t=0,numtimes-1 do begin
;  signal[t,*]=signal[t,*]/background[t,0]
;endfor

;
;  Replace negatives and zeros with really small number
;
result=where(signal le 0.0,count)
if count gt 0 then signal[result]=1.e-5
;
;  Perform an r squared correction
;
for t=0,numtimes-1 do begin
;  signal[t,*]=(alog10(signal[t,*]))+(2.*alog10(height))
  signal[t,*]=(10.*alog10(signal[t,*]))+(20.*alog10(height))
endfor
;
;  Data max and min
;
print,min(signal)
if min(signal) le -20.0 then begin
  result=where(signal lt -20.0,count)
  if count gt 0 then signal[result]=-20.0
endif
print,max(signal)
if max(signal) gt 20.0 then begin
  result=where(signal gt 20.0,count)
  if count gt 0 then signal[result]=20.0
endif
;
;  Convert mpl_bases to km
;
mpl_base=mpl_base*km_flag


;**************************
;  Set up the plot
;**************************

!p.background=!d.n_colors-1
loadct,5
gamma_ct,1.65
!p.charsize=1.5
!p.font=-1
window, 0, xsize=800, ysize=400, TITLE='MPL',/pixmap
!p.multi = [0,1,1]
!x.range=[julian_day[0],julian_day[numtimes-1]]
!x.style=1
!y.style=1
if input_flag eq 'eos' then $
  dummy=label_date(date_format=['%N/%D/%Z!C%H:%I:%S']) $
else $
  dummy=label_date(date_format=['%H:%I'])

;xl=0.075 & yl=0.30 & xr=0.85 & yr=0.75
position1=[0.1,0.15,0.85,0.85]
;
;  Plot 1  -  Backscatter
;

;
;  Find the data max and min
;
dmax=max(signal)
dmin=min(signal)
dmax=20.0
dmin=-20.0
;
;  Set up the plot coordinates
;
contour,signal,julian_day,height,/nodata,xstyle=4,ystyle=4,$
  position=position1
px=!x.window*!d.x_vsize
py=!y.window*!d.y_vsize
sx=px[1]-px[0]+1
sy=py[1]-py[0]+1
;
;  Scale and create the image
;
image=bytscl(signal)
tv,congrid(image,sx,sy),px[0],py[0]
;
;  Set up the axes
;
contour,signal,julian_day,height,/noerase,/nodata,$
  xtickformat='label_date',xtickunits='time',$
  color=0,xticklen=-0.02,yticklen=-0.02,position=position1,$
  yrange=[height[0],height[n_elements(height)-1]],$
  ytitle='Height (km)',xtitle='Time (UTC)'
;
;  Add a color bar
; 
colorbar,image,dmax,dmin,'Normalized Signal',px,py,sx,sy
;
;  Add title to plot page
;
if site eq 'twpc1' then capsite='TWP C1'
if site eq 'twpc2' then capsite='TWP C2'
if site eq 'nsa' then capsite='NSA'
if site eq 'sgp' then capsite='SGP'
a_string=capsite+' MicroPulse Lidar Observations, '+date_string
xyouts, 0.5, 0.93, alignment=0.5,a_string, /normal, color=0, charsize=2.2
xyouts, 0.5, 0.88, alignment=0.5,datastream, /normal, color=0, charsize=1.9
;
;  Create the gif file
;

if input_flag eq 'eos' then begin

  gif_file='/home/meteor/y/mace/homepages/research/eos/images/'+$
    site+'/'+datastream+'.'+curdate
  write_gif,gif_file+'.gif',tvrd()
  unix_command='rm '+gif_file+'.png'
  spawn,unix_command

  openw,13,'/data/mace4/arm_data/eos/eos_image_name.dat'
  printf,13,datastream+'.'+curdate+'.gif'
  close,13

endif else $
  write_gif,'/data/mace4/arm_data/temp_image/'+$
  datastream+'.'+curdate+'.gif',tvrd()

end

;****************************************************
;  colorbar.pro creates a colorbar for the DBZ images
;****************************************************

pro colorbar, image,dmax,dmin,cb_title,px,py,sx,sy

min_image=0;min(image)
max_image=255;max(image)

cb_arr=fltarr(2,max_image-min_image+1)
for j=0,max_image-min_image do $
  cb_arr(*,j)=dmin+((dmax-dmin)/(max_image-min_image))*j

cb_img=bindgen(2,max_image-min_image+1)

for j=0,max_image-min_image do cb_img(*,j)=min_image+j

tv, congrid(cb_img,(sx/35),sy),px[1]+sx/7,py[0]
axis,px[1]+sx/7,py[0],yaxis=0,color=0,$
  yrange=[cb_arr[0,0],cb_arr[0,max_image-min_image]],$
  ystyle=1,/device,charsize=1.2,ytitle=cb_title,ticklen=-0.02

return 
end