;  Open trajectory file
;  Loop through each point
;  Compare the point to cloudsat orbit track
;  Print out the closest overpass in time and space 

pro trajectory_satellite_ovp_cs

;*************************************
;  Input
;*************************************
;  Directories
;path_prefix='/Volumes/'   ;imac
path_prefix='/uufs/chpc.utah.edu/common/home/'  ;chpc

;  Trajectory file directory
tdir='u0079358/public_html/realm/mod06/lo_traj/'
traj_files=file_search(path_prefix+tdir+'traj_mod06*txt',count=ncount)

;  Loop through all the trajectory files in the directory
for n=0,ncount-1 do begin
  ;  Trajectory file name
  tfile=traj_files[n]
  tname=file_basename(tfile)
;tname='traj_mod06_20170112_07_-64.02_100.09.txt'
;tname='traj_mod06_20170113_07_-65.44_119.37.txt'
;tname='traj_mod06_20170109_05_-65.08_135.15.txt'
;tname='traj_mod06_20170109_05_-64.50_138.98.txt'
;tname='traj_mod06_20170121_09_-65.54_82.86.txt'
;tname='traj_mod06_20170117_08_-64.59_91.33.txt'
;tname='traj_mod06_20170109_05_-64.48_137.00.txt'

;  Absolute path to trajectory file
;tfile=path_prefix+tdir+tname

;*************************************
;  Read the trajectory file
;*************************************
read_hysplit_traj,tfile,julian_day_traj,age_traj,$
  lat_traj,lon_traj,ht_traj,pres_traj,lat_start,lon_start

;*************************************
;  Plotting region
;*************************************
;llat=10 & ulat=40 & llon=-179.0 & rlon=-140.0 ;kilauea
llat=-76 & ulat=-45 & llon=40 & rlon=152 ;antarctic region
  
;  Colortable
ct=['red','orange','green','blue','saddle brown','dodger blue','spring green',$
    'gray','crimson','olive','purple','magenta','hot pink','navy','purple','teal',$
    'dark goldenrod','tan','aqua','forest green']

;  Map the trajectory file
pxdim=900  ;  pixel size of the image in the x direction
pydim=1000   ;  pixel size of the image in the y direction
;  These values divide up the plot space using the subroutine
;  position plots.  The location of the plots will be defined by
;  pos
xl=0.10 & xr=0.95   
yb=0.05 & yt=0.45
sx=0.05
sy=0.07
numplots_x=1
numplots_y=1
position_plots,xl,xr,yb,yt,sx,sy,numplots_x,numplots_y,pos
dx=0.01 & dy=0.01
map_pos=pos[0,*]
m0=map('mercator',$
  limit=[llat,llon,ulat,rlon],/buffer,$
  box_axes=1,label_position=0,linestyle='dotted',$
  position=map_pos,font_size=5,clip=0)
mc=mapcontinents(/countries,thick=1.0)
mc=mapcontinents(/USA,thick=1.0)
;  Plot the track projection
p1=plot(lon_traj,lat_traj,/overplot)
;  Plot the track initial trajectory
p0=plot(lon_start,lat_start,/overplot,color='black',symbol='diamond',thick=1.5,sym_size=1)
parts=strsplit(tname,'txt',/extract)
tfile_str='t'+parts[0]
;m0.save,tfile_str+'.png',height=pydim

;*************************************************
;  Loop through all the points in the trajectory
;*************************************************
k=-1
prev_geop_file=''
;  Loop through each point in the trajectory
for t=0,n_elements(lat_traj)-1 do begin
  
;  Pick out the point we are looking at  
julian_day_t=julian_day_traj[t]
lat_t=lat_traj[t]
lon_t=lon_traj[t]

;  Small region around the point
ulat_st=lat_t+5.0
llat_st=lat_t-5.0
rlon_st=lon_t+5.0
llon_st=lon_t-5.0
print,llat_st,ulat_st,llon_st,rlon_st
  
;  Get the date strings
caldat,julian_day_t,tmonth,tday,tyear,thour,tmin,tsec

;  Calculate day of year 0-365
julian_date,tyear,tmonth,tday,tdoy

;  strings
tyear_str=string(tyear,format='(I4)')
tmonth_str=string(tmonth,format='(I02)')
tday_str=string(tday,format='(I02)')
thour_str=string(thour,format='(I02)')
tmin_str=string(tmin,format='(I02)')
tsec_str=string(tsec,format='(I02)')
tdoy_str=string(tdoy,format='(I03)')
;  yyyymmdd
tdate_str=tyear_str+tmonth_str+tday_str
;  hhmmss
ttime_str=thour_str+tmin_str+tsec_str

;******************************************************************************************
;  List of 3 consecutive days
;******************************************************************************************
;  Day before
julian_day_tb=julian_day_t-1D
caldat,julian_day_tb,tbmonth,tbday,tbyear,tbhour,tbmin,tbsec
julian_date,tbyear,tbmonth,tbday,tbdoy
tbyear_str=string(tbyear,format='(I4)')
tbdoy_str=string(tbdoy,format='(I03)')
;  Day after
julian_day_ta=julian_day_t+1D
caldat,julian_day_ta,tamonth,taday,tayear,tahour,tamin,tasec
julian_date,tayear,tamonth,taday,tadoy
tayear_str=string(tayear,format='(I4)')
tadoy_str=string(tadoy,format='(I03)')
;  Day after
julian_day_ta2=julian_day_t+2D
caldat,julian_day_ta2,ta2month,ta2day,ta2year,ta2hour,ta2min,ta2sec
julian_date,ta2year,ta2month,ta2day,ta2doy
ta2year_str=string(ta2year,format='(I4)')
ta2doy_str=string(ta2doy,format='(I03)')
;  Array of the 3 consecutive doys and years
year_str=[tbyear_str,tyear_str,tayear_str,ta2year_str]
doy_str=[tbdoy_str,tdoy_str,tadoy_str,ta2doy_str]


;  loop through the days and get a list of files and their start times
geop_jdays=!NULL & geop_files=!NULL
for i=0,n_elements(doy_str)-1 do begin
  cdir=path_prefix+'mace-group6/cloudsat/2B-GEOPROF.P1_R05/'+year_str[i]+'/'+doy_str[i]+'/'
  cstr='*.hdf'
  ;  Get a list of geoprof files
  cfiles=file_search(cdir+cstr,count=num_c)
  ;  This makes an array of files and an array of file starting times in julian days
  for j=0,num_c-1 do begin
    ;  Extract the file date and time from filename
    parts = strsplit(file_basename(cfiles[j]),'.',/extract)
    year=fix(strmid(parts[0],0,4))
    doy=fix(strmid(parts[0],4,3))
    hh=fix(strmid(parts[0],7,2))
    mm=fix(strmid(parts[0],9,2))
    ss=fix(strmid(parts[0],11,2))
    day_of_year_to_date,doy,year,month,day
    jday=julday(month,day,year,hh,mm,ss)
    ;print,file_basename(cfiles[j])
    ;print,year,month,day,hh,mm,ss
     ;  Accumulate these start times and filenames in an array
    if geop_jdays eq !NULL then begin
      geop_jdays=jday
      geop_files=cfiles[j]
    endif else begin
      geop_jdays = [geop_jdays,jday]
      geop_files = [geop_files,cfiles[j]]
    endelse
  endfor   ;end of loop through files
endfor  ;end of loop through consecutive days

;  sort the files by closest in time
abs_jday=abs(julian_day_t-geop_jdays)
sidx=sort(abs_jday)
sort_geop_files=geop_files[sidx]
sort_geop_times=geop_jdays[sidx]

;  Starting with the closest in time file, look for files with 
;  satellite tracks in the region of interest
m=0
count_ovp=0
;  Loop until we find a track in our region
while count_ovp eq 0 and m lt n_elements(sort_geop_files) do begin

  close_geop_file=sort_geop_files[m]
  close_geop_time=sort_geop_times[m]
  
  ;  Open geoprof and get lat and lon
  lat_cs=!NULL & lon_cs=!NULL
  get_cloudsat,'',close_geop_file,'Latitude',lat_cs
  get_cloudsat,'',close_geop_file,'Longitude',lon_cs

  ;  Caculate colongitude because the pacific crosses the 180 line
  colon_cs=lon_cs
  result=where(lon_cs lt 0,count)
  if count gt 0 then colon_cs[result]=colon_cs[result]+360.0

  ;  Open geoprof to get time 
  get_cloudsat,'',close_geop_file,'TAI_start',base_time
  get_cloudsat,'',close_geop_file,'Profile_time',time_offset

  ;  Julian day is calculated
  sec_per_day=long(24L*60L*60L)
  day1=julday(1,1,1993,0,0,0)
  julian_day_cs=double(day1+((base_time+time_offset)/sec_per_day))

  ;  Look for cloudsat track points in the region
  bidx=where(lat_cs ge llat and lat_cs le ulat and $  ;big map- need this for plotting
    lon_cs ge llon and lon_cs le rlon,bcount_ovp)
  cidx=where(lat_cs ge llat_st and lat_cs le ulat_st and $  ;small region around traj point
    lon_cs ge llon_st and lon_cs le rlon_st,count_ovp)
  
  print,m,' ',count_ovp,close_geop_file  
  
  ;  If there are no points, go get the next closest file
  if count_ovp eq 0 then m=m+1
endwhile

;  If there are points in the region, plot it
if close_geop_file ne prev_geop_file and count_ovp gt 0 then begin
  prev_geop_file=close_geop_file
  print,'plot this one'

  ;index for my color table
  k=k+1
  
  parts=strsplit(file_basename(close_geop_file),'_',/extract)
  orbit=parts[1]
  sym_str=tyear_str+tmonth_str+tday_str+' '+thour_str+':'+tmin_str
  p1=symbol(lon_t,lat_t,'star',/data,sym_color=ct[k],$;label_string=sym_str,$
    label_font_size=4,label_color=ct[k])
  t1=text(0.1,0.97-k*2.0*dy,'trajectory point '+sym_str+' '+string(lat_t,format='(F6.2)')+'lat'+$
    ' '+string(lon_t,format='(F7.2)')+' lon',color=ct[k],font_size=6)
  ovp_time=julian_day_cs[cidx[0]]
  caldat,ovp_time,omonth,oday,oyear,ohour,omin,osec
  oyear_str=string(oyear,format='(I4)')
  omonth_str=string(omonth,format='(I02)')
  oday_str=string(oday,format='(I02)')
  ohour_str=string(ohour,format='(I02)')
  omin_str=string(omin,format='(I02)')
  osec_str=string(osec,format='(I02)')
  ovp_str=oyear_str+omonth_str+oday_str+' '+ohour_str+':'+omin_str+' '+string(orbit)
  p1=plot(lon_cs[bidx],lat_cs[bidx],/overplot,color=ct[k])   
  ;t1=text(min(lon_cs[cidx])+1,min(lat_cs[cidx])+1,ovp_str,color=ct[k],font_size=5,/data,orientation=90)
  t1=text(0.5,0.97-k*2.0*dy,'sat ovp '+ovp_str,color=ct[k],font_size=6)
  m0.save,tfile_str+'cloudsat.png',height=pydim
  print,tfile_str+'cloudsat.png'
endif
endfor  ;end of loop through trajectory points    
  
endfor  ;end of loop through trajectory files    
                         

end