;+
; Project     : STEREO - PLASTIC
;
; Name        : plastic_plot_beacon
;
; Purpose     : create images of plastic beacon data
;
; Category    : Quicklook, Telemetry, Graphics, Time-series
;
; Explanation : Looks for PLASTIC beacon data for the last 
;               24 hours prior to date. Creates image files in .ps and
;               .png format for moments and heavy ion counts.
;
; Syntax      : IDL> plastic_plot_beacon
;
; Examples    : IDL> plastic_plot_beacon
;
; Inputs      : none
;
; Opt. Inputs : none
;
; Outputs     : none 
;
; Opt. Outputs: none
;
; Keywords    : date       Date to plot (together with previous day).  Defaults
;                          to today. ex: '2007-01-20'
;
;               cdf_filepath   Path to read files from.  Defaults to current
;                              directory.
;               image_filepath Path to write image files to. Defaults
;                              to current directory.
;
; Common      : none
;
; Restrictions: requires plastic_plot_stamp routine
;
; Side effects: creates .ps and .png files in current directory
;
; History     : version 1, 11/18/05, LB Ellis, written
;               version 2, 12/01/05, W. Thompson, added keywords date,filepath
;               version 3, 12/01/05, LB Ellis, changed Pressure to Temperature
;                    broke filepath into two keywords
;               version 4, 12/13/06, LB Ellis, new moments
;               version 6, 01/29/07, change in CDF version
;               version 7, 05/21/07, change in CDF version and update
;               velocity to _Spcrft instead of _Inst
;
; Contact     : Lorna Ellis (Lorna.Ellis@unh.edu)
;-


PRO plastic_plot_beacon, date = date, cdf_filepath = cdf_filepath, image_filepath = image_filepath

; find which files needed

if n_elements(date) eq 1 then utc = anytim2utc(date, /external) else $
  get_utc, utc, /external
today = anytim2cal(utc,form=8,/date)
today_filename_a = "STA_LB_PLASTIC_"+today+"_V08.cdf"
today_filename_b = "STB_LB_PLASTIC_"+today+"_V08.cdf"
today_st = strmid(today,4,2) + '/' + strmid(today,6,2) + '/' + $
  strmid(today,0,4)
utc.day = utc.day-1
check_ext_time, utc
yesterday = anytim2cal(utc,form=8,/date)
yesterday_filename_a = "STA_LB_PLASTIC_"+yesterday+"_V08.cdf"
yesterday_filename_b = "STB_LB_PLASTIC_"+yesterday+"_V08.cdf"
yesterday_st = strmid(yesterday,4,2) + '/' + strmid(yesterday,6,2) + '/' + $
  strmid(yesterday,0,4)
if n_elements(cdf_filepath) eq 1 then begin
    today_filename_a = concat_dir(cdf_filepath,today_filename_a)
    today_filename_b = concat_dir(cdf_filepath,today_filename_b)
    yesterday_filename_a = concat_dir(cdf_filepath,yesterday_filename_a)
    yesterday_filename_b = concat_dir(cdf_filepath,yesterday_filename_b)
endif

; set up variables
num_min = 2880 ; two days
density_out = replicate(!values.f_nan, num_min)
velocity_out = replicate(!values.f_nan, 4, num_min)
start_hour = utc.hour
IF start_hour EQ 24 THEN start_hour = 0

; open and read files 
; yesterday's data
loadct, 39
FOR loop = 0, 1 DO BEGIN 
    IF loop EQ 0 THEN spacecraft = 'a' ELSE spacecraft = 'b'
    ; yesterday's data
    IF loop EQ 0 THEN temp_name = yesterday_filename_a ELSE temp_name = yesterday_filename_b
    IF loop EQ 1 THEN BEGIN ; reset variables
        density_out[*] = !values.f_nan
        velocity_out[*, *] = !values.f_nan
    ENDIF 

    IF file_test(temp_name) EQ 1 THEN BEGIN ; file exists
                                ; read cdf
        cdf_id = cdf_open(temp_name)
        cdf_control, cdf_id, get_var_info = info, variable = 'Epoch1'
        cdf_varget, cdf_id, 'Epoch1', epoch_cdf, rec_count = info.maxrec
        cdf_varget, cdf_id, 'Density', density_cdf, rec_count = info.maxrec
        cdf_varget, cdf_id, 'Velocity_HGRTN', velocity_cdf, rec_count = info.maxrec
        cdf_close, cdf_id
                                ; find which data we need
        cdf_i = 0
        out_i = 0
        cdf_epoch, epoch_cdf[cdf_i], epoch_year, epoch_month, epoch_day, epoch_hour, epoch_min, epoch_sec, /breakdown_epoch
                                ; get to start_hour
        WHILE epoch_hour LT start_hour AND cdf_i LT info.maxrec DO BEGIN 
            cdf_i = cdf_i + 1
            IF cdf_i LT info.maxrec THEN $
              cdf_epoch, epoch_cdf[cdf_i], epoch_year, epoch_month, epoch_day, epoch_hour, epoch_min, epoch_sec, /breakdown_epoch
        ENDWHILE 
                                ; get data from first day
        WHILE cdf_i LT info.maxrec DO BEGIN 
            density_out[out_i] = density_cdf[cdf_i]
            velocity_out[0:3, out_i] = velocity_cdf[0:3, cdf_i]
            old_time = epoch_hour + (epoch_min * 60) + (epoch_sec * 3600)
            cdf_i = cdf_i + 1
            out_i = out_i + 1
            IF cdf_i LT info.maxrec THEN BEGIN 
                cdf_epoch, epoch_cdf[cdf_i], epoch_year, epoch_month, epoch_day, epoch_hour, epoch_min, epoch_sec, /breakdown_epoch
                new_time = epoch_hour + (epoch_min * 60) + (epoch_sec * 3600)
                WHILE new_time GT old_time + 75 AND cdf_i LT info.maxrec DO BEGIN ; check that within 1 min and 15 seconds
                    old_time = new_time
                    cdf_i = cdf_i + 1
                    out_i = out_i + 1
                    IF cdf_i LT info.maxrec THEN $
                      cdf_epoch, epoch_cdf[cdf_i], epoch_year, epoch_month, epoch_day, epoch_hour, epoch_min, epoch_sec, /breakdown_epoch
                    new_time = epoch_hour + (epoch_min * 60) + (epoch_sec * 3600)
                ENDWHILE 
            ENDIF 
        ENDWHILE 
    ENDIF 

    hour = 0
    out_i = (24-start_hour)*60
    
; today's data
    IF loop EQ 0 THEN temp_name = today_filename_a ELSE temp_name = today_filename_b
    IF file_test(temp_name) EQ 1 THEN BEGIN ; file exists
                                ; read cdf
        cdf_i = 0
        cdf_id = cdf_open(temp_name)
        cdf_control, cdf_id, get_var_info = info, variable = 'Epoch1'
        cdf_varget, cdf_id, 'Epoch1', epoch_cdf, rec_count = info.maxrec
        cdf_varget, cdf_id, 'Density', density_cdf, rec_count = info.maxrec
        cdf_attget, cdf_id, 'UNITS', 'Density', density_units
        cdf_varget, cdf_id, 'Velocity_HGRTN', velocity_cdf, rec_count = info.maxrec
        cdf_attget, cdf_id, 'UNITS', 'Velocity_HGRTN', velocity_units
        cdf_close, cdf_id
                                ; find which data we need
        cdf_epoch, epoch_cdf[cdf_i], epoch_year, epoch_month, epoch_day, epoch_hour, epoch_min, epoch_sec, /breakdown_epoch
        WHILE cdf_i LT info.maxrec AND out_i LT num_min DO BEGIN 
            density_out[out_i] = density_cdf[cdf_i]
            velocity_out[0:3, out_i] = velocity_cdf[0:3, cdf_i]
            old_time = epoch_hour + (epoch_min * 60) + (epoch_sec * 3600)
            cdf_i = cdf_i + 1
            out_i = out_i + 1
            IF cdf_i LT info.maxrec THEN BEGIN 
                cdf_epoch, epoch_cdf[cdf_i], epoch_year, epoch_month, epoch_day, epoch_hour, epoch_min, epoch_sec, /breakdown_epoch
                new_time = epoch_hour + (epoch_min * 60) + (epoch_sec * 3600)
                WHILE new_time GT old_time + 75 AND cdf_i LT info.maxrec DO BEGIN ; check that within 1 min and 15 seconds
                    old_time = new_time
                    cdf_i = cdf_i + 1
                    out_i = out_i + 1
                    IF cdf_i LT info.maxrec THEN $
                      cdf_epoch, epoch_cdf[cdf_i], epoch_year, epoch_month, epoch_day, epoch_hour, epoch_min, epoch_sec, /breakdown_epoch
                    new_time = epoch_hour + (epoch_min * 60) + (epoch_sec * 3600)
                ENDWHILE 
            ENDIF 
        ENDWHILE 
    ENDIF

; get ready for plots
    dummy = label_date(date_format = "%H:%I")    
    mytimes = timegen(num_min, units = 'Minutes', start = julday(utc.month, utc.day, utc.year, start_hour))
    mydate = "Hours: "+yesterday_st+" - "+today_st
; plot moments
; first page (Density, Velocity)
    mydevice = !D.NAME
    SET_PLOT, 'PS'
    out_name = 'moment_' + spacecraft
    IF n_elements(image_filepath) EQ 1 THEN out_name = concat_dir(image_filepath, out_name)
    DEVICE, FILENAME = out_name+'.ps', /COLOR, PORTRAIT = 1, /INCHES, xsize = 8.5, ysize = 11, XOFFSET = 0, YOFFSET = 0
    !p.noerase = 1
    header = 'PLASTIC Beacon Moments: Density & Velocity, Spacecraft '+strupcase(spacecraft)
    xyouts, .15, .96, header, /normal, size = 1.5
    dummy = label_date(date_format = " ")
    ; density
    temp = where(finite(density_out) EQ 1, count)
    ytitle = 'Density!c!c'+density_units
    IF count GT 0 THEN $
      ;plot, mytimes, density_out, xtickunits = "time", xtickformat = 'label_date', ytitle = ytitle, /ylog, yrange = [.1, 100], $
       ;     color = 0, background = 255, position = [.15, .53, .9, .94] $
      plot, mytimes, density_out, xtickunits = "time", xtickformat = 'label_date', ytitle = ytitle, /ylog, $
            color = 0, background = 255, position = [.15, .53, .9, .94] $
    ELSE $
      plot, mytimes, /nodata, xtickunits = "time", xtickformat = 'label_date', ytitle = ytitle, /ylog, yrange = [.1, 100], $
            color = 0, background = 255, position = [.15, .53, .9, .94] 

    ; total velocity
    dummy = label_date(date_format = "%H:%I")
    temp = where(finite(velocity_out[3, *]) EQ 1, count)
    ytitle = 'Velocity!c!c'+velocity_units
    IF count GT 0 THEN $
      ;plot, mytimes, velocity_out[3, *], xtickunits = "time", xtickformat = 'label_date', yrange = [0, 2000], $
       ;     xtitle = mydate, ytitle = ytitle, color = 0, background = 255, position = [.15, .1, .9, .51] $
      plot, mytimes, velocity_out[3, *], xtickunits = "time", xtickformat = 'label_date', $
            xtitle = mydate, ytitle = ytitle, color = 0, background = 255, position = [.15, .1, .9, .51] $
    ELSE $
      plot, mytimes, /nodata, xtickunits = "time", xtickformat = 'label_date', yrange = [0, 2000], $
            xtitle = mydate, ytitle = ytitle, color = 0, background = 255, position = [.15, .1, .9, .51] 
    plastic_plot_stamp
    device, /close
    set_plot, mydevice
    !p.noerase = 0
    temp_name = out_name
    command_string = 'convert -density 288 -geometry 25% -sharpen 3 ' + out_name + '.ps ' + temp_name + '.png'
    spawn, command_string

ENDFOR 

END