;+
; NAME:
;    THM_ASI_CREATE_MOSAIC
;
; PURPOSE:
;    create mosaic with all THEMIS ASI
;
; CATEGORY:
;    None
;
; CALLING SEQUENCE:
;    THM_ASI_CREATE_MOSAIC,time
;
; INPUTS:
;    Time	like '2006-01-01/05:00:00'
;
; OPTIONAL INPUTS:
;    None
;
; KEYWORD PARAMETERS:
;    cal_files	calibration files if they do not need to be read
;    gif_out    create a gif-file
;    verbose    print some diagnostics
;    pgm_file   do not read CDF, but pgm-files
;    zbuffer    do in z-buffer, not on the screen
;    thumb      use thumbnails instead of full resolution images
;    exclude    string of station names that should not be plotted
;    show       string of station names that should only be plotted
;    top        top color to be used for polyfill
;    special    special treatment for not mapped thumbnail
;    scale               scale for map set
;    central_lon         geographic longitude of center of plot
;    central_lat         geographic latitude of center of plot
;    xsize               xsize of window
;    ysize               ysize of window
;    color_continent     shade of continent fill
;    color_background    shade of background
;    position=position   position of plot on window (normal coordinates)
;    noerase=noerase     do not erase current window (no effect if {x,y}size set
;    cursor     finish with cursor info
;    projection projection for map set, MAP_PROJ_INFO, PROJ_NAMES=names
;    window     set window number
;    rotation	rotate map
;    full_minute	full minute for pgm
;    minimum_elevation	minimum elevation to plot in degrees
;
; OUTPUTS:
;    None
;
; OPTIONAL OUTPUTS:
;    None
;
; COMMON BLOCKS:
;    None
;
; SIDE EFFECTS:
;    None
;
; RESTRICTIONS:
;    None
;
; EXAMPLE:
;    THM_ASI_CREATE_MOSAIC,'2006-01-01/05:00:00'
;    THM_ASI_CREATE_MOSAIC,'2007-03-01/04:00:00',/thumb,exclude='ekat'
; 
;
; MODIFICATION HISTORY:
;    Written by: Harald Frey, 02/06/2007
;                based on example from Donovan/Jackel
;                2007-03-15, hfrey, thumbnails, keyword exclude
;                2007-03-27, hfrey, special treatment for not mapped
;                thumbsnails
;                2007-12-21, jmm, added explicit set_plot,'z' for zbuffer
;                2008-07-21, jmm, added gif_dir, for output directory option
; NOTES:
;     
; VERSION:
;   $LastChangedBy: hfrey $
;   $LastChangedDate: 2008-10-24 13:43:00 -0700 (Fri, 24 Oct 2008) $
;   $LastChangedRevision: 3811 $
;   $URL: svn+ssh://thmsvn@ambrosia.ssl.berkeley.edu/repos/thmsoc/tags/tdas_5_1/idl/themis/ground/thm_asi_create_mosaic.pro $
;
;-

PRO THM_ASI_CREATE_MOSAIC,time,cal_files=cal_files,gif_out=gif_out,$
    verbose=verbose,pgm_file=pgm_file,thumb=thumb,$
    exclude=exclude,top=top,special=special,show=show,$
    scale=scale,$                      ; scale for map set
    central_lat=central_lat,$          ; geographic latitude of center of plot
    central_lon=central_lon,$          ; geographic longitude of center of plot
    color_continent=color_continent,$  ; shade of continent fill
    color_background=color_background,$; shade of background
    position=position,$                ; position of plot on window (normal coordinates)
    xsize=xsize,$                      ; xsize of window
    ysize=ysize,$                      ; ysize of window
    noerase=noerase,$                  ; do not erase current window (no effect if {x,y}size set
    zbuffer=zbuffer,$		       ; do it in z-buffer
    cursor=cursor,$		       ; finish with cursor info
    projection=projection,$	       ; map projection
    maxval=maxval,$		       ; brightness scaling of images
    minval=minval,$
    window=window,$      		; set window number
    rotation=rotation,$
    full_minute=full_minute,$
    minimum_elevation=minimum_elevation,$
    merge=merge,$			; merge full resolution and thumbnail images
    gif_dir=gif_dir                  ; An output directory for the gif output, default is the local working dir.

	; input check
if keyword_set(verbose) then verbose=systime(1)
if (strlen(time) ne 19) then begin
  print,'Wrong time input'
  print,'Correct value like 2006-01-01/05:00:00'
  return
  endif

	; check brightness scaling
if keyword_set(maxval) then begin
  if not keyword_set(show) then begin
     print,'Show has to be set with maxval'
     return
     endif
  if n_elements(show) ne n_elements(maxval) then begin
     print,'N_elements of show and maxval have to match'
     return
  endif
  endif  

	; strip time
year=strmid(time,0,4)
month=strmid(time,5,2)
day=strmid(time,8,2)
hour=strmid(time,11,2)
minute=strmid(time,14,2)
second=strmid(time,17,2)

	; setup
thm_init
timespan,year+'-'+month+'-'+day+'/'+hour,1,/hour
thm_asi_stations,site,location
if keyword_set(zbuffer) then map_scale=2.6e7 else map_scale=4.e7
if keyword_set(scale) then map_scale=scale
if not keyword_set(central_lon) then central_lon=255.
if not keyword_set(central_lat) then central_lat=63.
if not keyword_set(xsize) then xsize=700
if not keyword_set(ysize) then ysize=410
if not keyword_set(top) then top=254

	; characters
if keyword_set(zbuffer) then chars=1.15 else chars=1.5

	; some setup
if keyword_set(minimum_elevation) then minimum_elevation_to_plot=minimum_elevation else minimum_elevation_to_plot=8. ;degrees
if (keyword_set(thumb) and not keyword_set(special)) then n1=1024l else n1=256l*256l

     	; clean up before start
names=tnames('thg_asf*')
if (names[0] ne '') then store_data,delete=names
names=tnames('thg_ast*')
if (names[0] ne '') then store_data,delete=names

	; read available data
thm_mosaic_array,year,month,day,hour,minute,second,strlowcase(site),$
        image,corners,elevation,pixel_illuminated,n_sites,verbose=verbose,$
        cal_files=cal_files,pgm_file=pgm_file,thumb=thumb,special=special,$
        show=show,exclude=exclude,full_minute=full_minute,merge=merge

	; fill variables
bytval=fltarr(n_sites)+1.
bitval=fltarr(n_sites)
if keyword_set(maxval) then begin
  for i=0,n_elements(maxval)-1 do begin
    index=where(strlowcase(site) eq strlowcase(show[i]))
    bytval[index]=maxval[i]
    bitval[index]=minval[i]
    endfor
  for i=0,n_sites-1 do image[*,i]=bytscl(image[*,i],min=bitval[i],max=bytval[i])
  endif else begin
  if keyword_set(special) then for i=0,n_sites-1 do bytval[i]=(median(image[0:1023,i]) > 1) else begin
     if keyword_set(full_minute) then begin
        for i=0,n_sites-1 do bytval[i]=(median(image[*,i,*]) > 1) 
        for i=0,n_sites-1 do image[*,i,*]=((image[*,i,*]/bytval[i])*64.) < 254
        endif else begin ; prevent divide by zero
        for i=0,n_sites-1 do bytval[i]=(median(image[*,i]) > 1) ; prevent divide by zero
        mm_dumm=where(merge ne -1,count_mm_dumm)
        if (count_mm_dumm ne 0) then begin
           for dki=0,count_mm_dumm-1 do bytval[mm_dumm[dki]]=(median(image[0:1023,mm_dumm[dki]]) > 1)
           endif
        for i=0,n_sites-1 do image[*,i]=((image[*,i]/bytval[i])*64.) < 254
        endelse
     endelse
  endelse

	; no images found
if (max(bytval) eq 1.) then begin
  print,'No images for ',time
  return
  endif

	; exclude unwanted sites
if keyword_set(exclude) then begin
  for i=0,n_elements(exclude)-1 do begin
    not_site=where(strlowcase(site) eq strlowcase(exclude[i]),count)
    if (count eq 0) then print,'Not a valid site: ',exclude[i] else $
        corners[*,*,*,not_site]=!values.f_nan
    endfor
  endif

if keyword_set(full_minute) then begin
time_start=time

	; search for midnight file
	; does not change much over one minute
f=file_search('midnight.sav',count=count)
if (count eq 1) then begin
  midlons=fltarr(40)+!values.f_nan
  ut_hour=float(hour)+float(minute)/60.+float(second)/3600.
  restore,f
  for i=0,39 do begin
    lon=interpol(findgen(141)+start_longitude,reform(midnight[i,*]),ut_hour)
    midlons[i]=lon[0]
    endfor
  bad=where(midlons gt 360.,count)
  if (count gt 0) then midlons[bad]=!values.f_nan
  endif

	; run through images
for ikk=0,19 do begin

;zbuffer needs to be set before the loadct call in thm_map_set,
;otherwise this bombs the second time through because of reset to 'x'
;later in this program, jmm 21-dec-2007
if(keyword_set(zbuffer)) then set_plot, 'z'
; set up the map
thm_map_set,scale=map_scale,$
     central_lat=central_lat,$           ; geographic latitude of center of plot
     central_lon=central_lon,$           ; geographic longitude of center of plot
     color_continent=color_continent,$   ; shade of continent fill
     color_background=color_background,$ ; shade of background
     position=position,$                 ; position of plot on window (normal coordinates)
     xsize=xsize,$                       ; xsize of window
     ysize=ysize,$                       ; ysize of window
     noerase=noerase,$                   ; do not erase current window (no effect if {x,y}size set
     zbuffer=zbuffer,$
     projection=projection,$
     window=window,$
     rotation=rotation

	; normal filling
for pixel=0l,n1-1l do begin
  for i_site=0,n_sites-1 do begin
    if ((pixel_illuminated[pixel,i_site] eq 1) and $
           (elevation[pixel,i_site] gt minimum_elevation_to_plot)) then $
       polyfill,corners[pixel,[0,1,2,3,0],0,i_site],$
            corners[pixel,[0,1,2,3,0],1,i_site],color=image[pixel,i_site, ikk] < top
    if keyword_set(special) then if ((spec_illu[pixel*64 < (n1-1),i_site] eq 1) and $
           (elevation[pixel*64 < (n1-1),i_site] gt minimum_elevation_to_plot)) then begin
       for ijk=0,63 do polyfill,corners[pixel*64+ijk < (n1-1),[0,1,2,3,0],0,i_site],$
            corners[pixel*64+ijk < (n1-1),[0,1,2,3,0],1,i_site],color=image[pixel*64+ijk < (n1-1),i_site] < top
;       stop
       endif
    endfor
  endfor

	; finish map
return_lons=1
return_lats=1
thm_map_add,invariant_lats=findgen(4)*10.+50.,invariant_color=210,$
    invariant_linestyle=1,/invariant_lons,return_lons=return_lons,$
    return_lats=return_lats
time=time_string(time_double(time_start)+ikk*3)

xyouts,0.005,0.018,time,color=0,/normal,charsize=chars
xyouts,0.005,0.060,'THEMIS-GBO ASI',color=0,/normal,charsize=chars

	; plot midnight file
plots,smooth(midlons-360.,5),findgen(40)+40.,color=255,/data

	; gif output
if keyword_set(gif_out) then begin
   If(keyword_set(gif_dir)) Then gdir = gif_dir Else gdir = './';jmm, 21-jul-2008
   tvlct,r,g,b,/get
   tvlct,r,g,b,/get
   img=tvrd()
	; strip time
year=strmid(time,0,4)
month=strmid(time,5,2)
day=strmid(time,8,2)
hour=strmid(time,11,2)
minute=strmid(time,14,2)
second=strmid(time,17,2)
   out_name='MOSA.'+year+'.'+month+'.'+day+'.'+hour+'.'+minute+'.'+second+'.gif'
   write_gif,gdir+out_name,img,r,g,b
   print,'Output in ',out_name
   endif
   
if keyword_set(zbuffer) then begin
     zbuffer=tvrd()
     device,/close
     set_plot,'x'
   endif

endfor
endif else begin	; end of full_minute loop

;zbuffer needs to be set before the loadct call in thm_map_set,
;otherwise this bombs the second time through because of reset to 'x'
;later in this program, jmm 21-dec-2007
if(keyword_set(zbuffer)) then set_plot, 'z'
; set up the map
thm_map_set,scale=map_scale,$
     central_lat=central_lat,$           ; geographic latitude of center of plot
     central_lon=central_lon,$           ; geographic longitude of center of plot
     color_continent=color_continent,$   ; shade of continent fill
     color_background=color_background,$ ; shade of background
     position=position,$                 ; position of plot on window (normal coordinates)
     xsize=xsize,$                       ; xsize of window
     ysize=ysize,$                       ; ysize of window
     noerase=noerase,$                   ; do not erase current window (no effect if {x,y}size set
     zbuffer=zbuffer,$
     projection=projection,$
     window=window,$
     rotation=rotation
     
	; fill it
if keyword_set(special) then begin
  spec_illu=pixel_illuminated-2
  for i=0,n_elements(special)-1 do begin
    spec_stat=where(strlowcase(site) eq strlowcase(special[i]))
    spec_illu[*,spec_stat]=pixel_illuminated[*,spec_stat]
    endfor
;  for pixel=0l,n1-1l do begin
;    for i_site=0,n_sites-1 do begin
;      if ((spec_illu[pixel,i_site] eq 1) and $
;           (elevation[pixel,i_site] gt minimum_elevation_to_plot)) then $
;       polyfill,corners[pixel,[0,1,2,3,0],0,i_site],$
;            corners[pixel,[0,1,2,3,0],1,i_site],color=image[pixel,i_site] < top
;       endfor
;    endfor
;  for i=0,n_elements(special)-1 do begin
;    spec_stat=where(strlowcase(site) eq strlowcase(special[i]))
;    pixel_illuminated[*,spec_stat]=-1
;    endfor
  endif

	; normal filling
for pixel=0l,n1-1l do begin
  for i_site=0,n_sites-1 do begin
    if ((pixel_illuminated[pixel,i_site] eq 1) and $
           (elevation[pixel,i_site] gt minimum_elevation_to_plot)) then $
       polyfill,corners[pixel,[0,1,2,3,0],0,i_site],$
            corners[pixel,[0,1,2,3,0],1,i_site],color=image[pixel,i_site] < top
    if keyword_set(special) then if ((spec_illu[pixel*64 < (n1-1),i_site] eq 1) and $
           (elevation[pixel*64 < (n1-1),i_site] gt minimum_elevation_to_plot)) then begin
       for ijk=0,63 do polyfill,corners[pixel*64+ijk < (n1-1),[0,1,2,3,0],0,i_site],$
            corners[pixel*64+ijk < (n1-1),[0,1,2,3,0],1,i_site],color=image[pixel*64+ijk < (n1-1),i_site] < top
;       stop
       endif
    endfor
;if pixel eq 500 then stop
  if (keyword_set(special) and pixel ge 1024) then break
  endfor

	; finish map
return_lons=1
return_lats=1
thm_map_add,invariant_lats=findgen(4)*10.+50.,invariant_color=210,$
    invariant_linestyle=1,/invariant_lons,return_lons=return_lons,$
    return_lats=return_lats
xyouts,0.005,0.018,time,color=0,/normal,charsize=chars
xyouts,0.005,0.060,'THEMIS-GBO ASI',color=0,/normal,charsize=chars

if keyword_set(verbose) then print,'After map: ',systime(1)-verbose,$
   ' Seconds'

	; search for midnight file
f=file_search('midnight.sav',count=count)
if (count eq 1) then begin
  midlons=fltarr(40)+!values.f_nan
;  early=fltarr(40)+!values.f_nan
;  late=fltarr(40)+!values.f_nan
  ut_hour=float(hour)+float(minute)/60.+float(second)/3600.
  restore,f
  for i=0,39 do begin
    lon=interpol(findgen(141)+start_longitude,reform(midnight[i,*]),ut_hour)
    midlons[i]=lon[0]
;    lon=interpol(findgen(141)+start_longitude,reform(midnight[i,*]),ut_hour-2.)
;    early[i]=lon[0]
;    lon=interpol(findgen(141)+start_longitude,reform(midnight[i,*]),ut_hour+2.)
;    late[i]=lon[0]
    endfor
  bad=where(midlons gt 360.,count)
  if (count gt 0) then midlons[bad]=!values.f_nan
  plots,smooth(midlons-360.,5),findgen(40)+40.,color=255,/data
;  plots,smooth(early-360.,5),findgen(40)+40.,color=255,/data,line=2
;  plots,smooth(late-360.,5),findgen(40)+40.,color=255,/data,line=2
  endif

;thm_map_add,invariant_lats=[65,70],invariant_lons=0.

;stop
	; gif output
if keyword_set(gif_out) then begin
   If(keyword_set(gif_dir)) Then gdir = gif_dir Else gdir = './';jmm, 21-jul-2008
   tvlct,r,g,b,/get
   img=tvrd()
   out_name='MOSA.'+year+'.'+month+'.'+day+'.'+hour+'.'+minute+'.'+second+'.gif'
   write_gif,gdir+out_name,img,r,g,b
   print,'Output in ',out_name
   endif
   
if keyword_set(zbuffer) then begin
     zbuffer=tvrd()
     device,/close
     set_plot,'x'
   endif

if keyword_set(cursor) then begin
   print,'Point cursor on map!'
   cursor,x,y,/data
   print,'Location: ',x,y
   endif

endelse	; single time



if keyword_set(verbose) then print,'Calculation took ',systime(1)-verbose,$
   ' Seconds'

heap_gc
;stop
end