;+
; NAME:
;
;
;
; PURPOSE:
;
;
;
; CATEGORY:
;
;
;
; CALLING SEQUENCE:
;
;
;
; INPUTS:
;
;
;
; OPTIONAL INPUTS:
;
;
;
; KEYWORD PARAMETERS:
;
;
;
; OUTPUTS:
;
;
;
; OPTIONAL OUTPUTS:
;
;
;
; COMMON BLOCKS:
;
;
;
; SIDE EFFECTS:
;
;
;
; RESTRICTIONS:
;
;
;
; PROCEDURE:
;
;
;
; EXAMPLE:
;
;
;
; MODIFICATION HISTORY:
;
;-
FUNCTION get_thm_pxxm, start_time, $ ; this is done in the themis format
                       n_days, $     ; how many days to grab data for
                       probe, $      ; which probe
                       ARRAY = ARRAY, $ ; create an array or structures (should use this)
                       VERSION = version, $ ; override the current version number
                       CLEAN = clean, $     ; clean up all the tplot variables created by this
                       NO_ORBIT = no_orbit

on_error, 0
compile_opt strictarr

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; This is input data, can be made as a pass from a procedure
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
SetDefaultValue, Version, 2     ; this is the savefile version that will be used to restore
thm_save_dir = '/Users/balarsen/Documents/Research/Data/THEMIS/save'

IF n_params() NE 3 THEN $
   message, 'Must enter start_time, duration, and probe'

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; Get the data
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
timespan,start_time, n_days, /day
IF file_test(/read, THM_SAVE_DIR+'/'+thm_save_name(probe, version)) EQ 0 THEN BEGIN 
   thm_load_mom, probe=probe
   ;; need to figure out how to only load the s/c potential and not the whole momnets set as its huge
   
   thm_load_state, probe=probe,/get_support_data
   thm_cotrans, 'th'+probe+'_state_pos', out_coord='sm', out_suffix='_sm'
   
   
;; grab and separate the ROI data to find the eclipse times
;; from thm_roi_bar
   get_data, 'th'+probe+'_state_roi', tha_state_roi_time, tha_state_roi, data = d
   
   bit_mask = [1,1,0,0,0,0,0,1,1,1,1,1,1,1,1,0,0,0,0,1]
   
   in = d.y
   out = ulonarr(n_elements(in))
   
   j = 0
;packs all the requested bits to the left
   for i = 0,n_elements(bit_mask)-1 do begin
      out += ishft(in and ishft(bit_mask[i],i),-j)
      if ~bit_mask[i] then j++
   endfor
   t = total(bit_mask)
   if t le 16 then begin
      out = uint(out)
   ENDIF
   
   d = {x:d.x,y:out}
   size = expand_in_base(max(d.y), /get_size)
   bits = bytarr(size, n_elements(d.y)) 
   FOR i = 0UL, n_elements(d.y)-1 DO BEGIN
      bits[*, i] = expand_in_base(d.y[i], size = size)
   ENDFOR
;; drop all data that is in eclipse
;; Earth Shadow 1
;; Lunar Shadow 2
;; Atmospheric Absorption Zone 4
;; South Atlantic Anomaly 8
;; Northern Auroral Zone 16
;; Southern Auroral Zone 32
;; Perigee Passage 64
;; Inner & Outer Radiation Belts 128
;; Deep Plasma Sphere 256
;; Foreshock Solar Wind 512
;; Solar Wind Beam 1024
;; High Magnetic Field 2048
;; Average Plasma Sheet 4096
;; Bowshock Crossing 8192
;; Magnetopause Crossing 16384
;; Ground Based Observatories 32768
;; 2-Day Conjunctions 65536
;; 4-Day Conjunctions 131072
   
   ind_ecl = where(bits[0, *] EQ 1 OR bits[1, *] EQ 1, n_ind_ecl)
   IF n_ind_ecl NE 0 THEN BEGIN 
      regions = label_region(reform(bits[0,*]))
      h = HISTOGRAM(regions, REVERSE_INDICES=ri)
      ecl_times = dblarr(n_elements(h)-1, 2) 
      FOR i = 1UL, n_elements(h)-1 DO BEGIN 
         p = ri[ri[i]:ri[i+1]-1]
         min = min(p,min_i, max = max, subscript_max = max_i)
         ecl_times[i-1, *] = [tha_state_roi_time[min], tha_state_roi_time[max]] 
                                ; -1 is since I started the loop at 1
      ENDFOR 
   ENDIF 
   
   
;; now onto the pxxm data
   get_data, 'th'+probe+'_pxxm_pot', tha_pxxm_pot_time, tha_pxxm_pot, data = d
   IF n_ind_ecl NE 0 THEN BEGIN 
      FOR i = 0UL, n_elements(ecl_times[*, 0])-1 DO BEGIN  
         ind = where(tha_pxxm_pot_time GE ecl_times[i, 0]-30 AND $
                     tha_pxxm_pot_time LE ecl_times[i, 1]+30, n_ind)
         ;; the -30 and +30 are 30 secodes on both sides of eclipse
         IF n_ind NE 0 THEN $
            tha_pxxm_pot[ind] = !values.f_nan   
      ENDFOR
   ENDIF 
   
;; store the corrected potential data
   store_data, 'th'+probe+'_pxxm_pot_trim', data ={x:tha_pxxm_pot_time, y:tha_pxxm_pot}
   
;; get the sm coordinates
   get_data, 'th'+probe+'_state_pos_sm', a, b, v, data=d
   
;; change units to Re and store it back
   d.y /= earthrad('km')
   store_data,  'th'+probe+'_state_pos_sm', data = d
   
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; This is an orbit diagnostic plot, might want this
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   plot, d.y[*, 0], d.y[*, 1], /iso,  $
         title = 'th'+probe+' ' +start_time + ' -> ' + trim(n_days) + ' days',  charsize = 1.2,  $
         xtitle = 'SM-x', ytitle = 'SM-y'
   plot_earth, color = 0   
   
   
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;  calculate L - centered dipole
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; ok, calculate R 
   R = sqrt(d.y[*, 0]*d.y[*, 0] + d.y[*, 1]*d.y[*, 1] + d.y[*, 2]*d.y[*, 2])
   maglat = atan(d.y[*, 2], sqrt(d.y[*, 0]*d.y[*, 0] + d.y[*, 1]*d.y[*, 1]))
   L = R/(cos(maglat)*cos(maglat))
;; store the l_shell and maglat
   store_data,  'th'+probe+'_maglat', data = {x:d.x, y:maglat}
   store_data,  'th'+probe+'_l_shell', data = {x:d.x, y:l}
;tplot, ['th'+probe+'_pxxm_pot', 'th'+probe+'_l_shell']         
   
;;;;;;;;;;;;;;;;;;;;;;;;
;;  calculate MLT
;;;;;;;;;;;;;;;;;;;;;;;;
;; start in GSE coords and then angle in x/y plane
   
   thm_cotrans, 'th'+probe+'_state_pos', out_coord='gse', out_suffix='_gse'
   get_data, 'th'+probe+'_state_pos_gse',  data = lmlt_data
   angle = atan(lmlt_data.y[*, 1], lmlt_data.y[*, 0])
   mlt = rad2mlt(angle)         ; assumes that 0 degrees is 1200 mlt
   
   store_data, 'th'+probe+'_l_mlt', data = {x:lmlt_data.x, y:[[l], [mlt]]}
   
;;;;;;;;;;;;;;;;;;;;;;;;;;
;; interpolate L and MLT to the same times as the pxxm data
;;;;;;;;;;;;;;;;;;;;;;;;;;
   
;tplot, ['th'+probe+'_pxxm_pot_trim', 'th'+probe+'_l_mlt'] 
   get_data,  'th'+probe+'_pxxm_pot_trim',  data = d1
   get_data, 'th'+probe+'_l_mlt',  data = d2
   
   l_int=interpol(d2.y[*,0],d2.x,d1.x)
   mlt_int=interpol(d2.y[*,1],d2.x,d1.x)
   
;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;  calculate and plot up 2d histograms
;;;;;;;;;;;;;;;;;;;;;;;;;;;
   
;;  Change to easier variables
   
   pxxm = d1.y
   time = d1.x
   lshell = l_int
   mlt = mlt_int

;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; calculate orbits and break them at apogee
;;;;;;;;;;;;;;;;;;;;;;;;;;;
   IF NOT keyword_set(no_orbit) THEN BEGIN 
      ind = where(lshell LT 8)
      ltmp = lshell
      ltmp[ind] = 0
      
      obrks = dblarr(10*n_days)
      reg=label_region(ltmp)
      FOR i = 1UL, max(reg) DO BEGIN
         ind = where(reg EQ i, n_ind)
         IF n_ind NE 0 THEN BEGIN 
            max = max(ltmp[ind], max_ind)
            tm = time[ind[max_ind]]
            plots, tm, max, psym = 1
            print, tm, max,  max_ind
            obrks[i-1] = tm
         ENDIF 
      ENDFOR 
      
      ind = where(obrks NE 0., n_ind)
      IF n_ind NE 0 THEN obrks = obrks[ind]
      
;; break the data into orbits at apogee
      o_num = uintarr(n_elements(time))
      ind = where(time LT obrks[0], n_ind)
      IF n_ind NE 0 THEN o_num[ind] = 1
      FOR i = 0UL, n_elements(obrks)-2 DO BEGIN  
         ind = where(time GE obrks[i] AND time LT obrks[i+1], n_ind)
         IF n_ind NE 0 THEN o_num[ind] = max(o_num)+1
      ENDFOR 
      ind = where(time GE obrks[i], n_ind)
      IF n_ind NE 0 THEN o_num[ind] = max(o_num)+1
      
;; break the orbits into inbound and outbound
      inbound = bytarr(n_elements(time)) 
      FOR i = 1UL, max(o_num) DO BEGIN 
         ind = where(o_num EQ i, n_ind)
         IF n_ind GT 1 THEN BEGIN 
            tmp = lshell[ind] - lshell[ind[1:*]]
            ind2 = where(tmp LT 0, n_ind2)
            IF n_ind2 NE 0 THEN BEGIN 
               inbound[ind[ind2]] = 1
            ENDIF 
         ENDIF 
      ENDFOR 
   ENDIF 
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; save the vars
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   
;; Make a savefile of the info we want
   save, pxxm, time, lshell, mlt,  $
         filename = THM_SAVE_DIR+'/'+thm_save_name(probe, version)
   
   
;; if there is a savefile just load it
ENDIF ELSE restore, THM_SAVE_DIR+'/'+thm_save_name(probe, version)

IF keyword_set(array) THEN BEGIN
   dat =  {pxxm:0., $
           time:0.d, $
           lshell:0.d, $
           mlt:0.d,  $
           START_TIME:'', $
           N_DAYS:0UL, $
           PROBE:''}

   dat = replicate(dat, n_elements(time))
   dat.pxxm = pxxm
   dat.time = time
   dat.lshell = lshell
   dat.mlt = mlt
   dat.START_TIME = start_time
   dat.PROBE = probe
ENDIF ELSE $ 
   dat = {pxxm:pxxm, $
          time:time, $
          lshell:lshell, $
          mlt:mlt,  $
          START_TIME:start_time, $
          N_DAYS:n_days, $
          PROBE:probe}

IF keyword_set(clean) THEN BEGIN
store_data, 'th'+probe+'_peim_density', /delete      
store_data, 'th'+probe+'_peim_flux', /delete 
store_data, 'th'+probe+'_peim_mftens', /delete
store_data, 'th'+probe+'_peim_eflux', /delete
store_data, 'th'+probe+'_peim_velocity', /delete
store_data, 'th'+probe+'_peim_ptens', /delete
store_data, 'th'+probe+'_peim_ptot', /delete
store_data, 'th'+probe+'_peem_density', /delete
store_data, 'th'+probe+'_peem_flux', /delete
store_data, 'th'+probe+'_peem_mftens', /delete
store_data, 'th'+probe+'_peem_eflux', /delete
store_data, 'th'+probe+'_peem_velocity', /delete
store_data, 'th'+probe+'_peem_ptens', /delete
store_data, 'th'+probe+'_peem_ptot', /delete
store_data, 'th'+probe+'_psim_density', /delete
store_data, 'th'+probe+'_psim_flux', /delete
store_data, 'th'+probe+'_psim_mftens', /delete
store_data, 'th'+probe+'_psim_eflux', /delete
store_data, 'th'+probe+'_psim_velocity', /delete
store_data, 'th'+probe+'_psim_ptens', /delete
store_data, 'th'+probe+'_psim_ptot', /delete
store_data, 'th'+probe+'_psem_density', /delete
store_data, 'th'+probe+'_psem_flux', /delete
store_data, 'th'+probe+'_psem_mftens', /delete
store_data, 'th'+probe+'_psem_eflux', /delete
store_data, 'th'+probe+'_psem_velocity', /delete
store_data, 'th'+probe+'_psem_ptens', /delete
store_data, 'th'+probe+'_psem_ptot', /delete
store_data, 'th'+probe+'_ptim_density', /delete
store_data, 'th'+probe+'_ptim_flux', /delete
store_data, 'th'+probe+'_ptem_mftens', /delete
store_data, 'th'+probe+'_ptim_eflux', /delete
store_data, 'th'+probe+'_ptim_velocity', /delete
store_data, 'th'+probe+'_ptim_ptens', /delete
store_data, 'th'+probe+'_ptim_ptot', /delete
store_data, 'th'+probe+'_ptem_density', /delete
store_data, 'th'+probe+'_ptem_flux', /delete
store_data, 'th'+probe+'_ptim_mftens', /delete
store_data, 'th'+probe+'_ptem_eflux', /delete
store_data, 'th'+probe+'_ptem_velocity', /delete
store_data, 'th'+probe+'_ptem_ptens', /delete
store_data, 'th'+probe+'_ptem_ptot', /delete
store_data, 'th'+probe+'_pxxm_pot', /delete
store_data, 'th'+probe+'_pxxm_qf', /delete
store_data, 'th'+probe+'_pxxm_shft', /delete
store_data, 'th'+probe+'_state_pos', /delete
store_data, 'th'+probe+'_state_vel', /delete
store_data, 'th'+probe+'_state_man', /delete
store_data, 'th'+probe+'_state_roi', /delete
store_data, 'th'+probe+'_state_spinras', /delete
store_data, 'th'+probe+'_state_spindec', /delete
store_data, 'th'+probe+'_state_spinalpha', /delete
store_data, 'th'+probe+'_state_spinbeta', /delete
store_data, 'th'+probe+'_state_spinper', /delete
store_data, 'th'+probe+'_state_spinphase', /delete
store_data, 'th'+probe+'_state_spin_spinper', /delete
store_data, 'th'+probe+'_state_pos_sm', /delete
store_data, 'th'+probe+'_pxxm_pot_trim', /delete
store_data, 'th'+probe+'_maglat', /delete
store_data, 'th'+probe+'_l_shell', /delete
store_data, 'th'+probe+'_state_pos_gse', /delete
store_data, 'th'+probe+'_l_mlt', /delete
ENDIF

return, dat

END