;+
; 
;             NAME : stereo_pointing
; 
;          PURPOSE : to use a point selected on a stereo image to estimate when a region
;                    will be observable with spacecraft located at Earth.
;
; CALLING SEQUENCE : stereo_pointing - click on "Get Stereo Browser" to download a recent
;                    browser image or "Read My Stereo File" to read a file that already
;                    exists on your disk.
; 
;           INPUTS : none
; 
;  OPTIONAL INPUTS : none
; 
;         KEYWORDS : none
; 
;          OUTPUTS : Pointings! Positions are displayed on the widget and written to a
;                    file (stereo_pointing.txt) in the currect directory. You can also
;                    make a PostScript file.
;
; OPTIONAL OUTPUTS : none
;
;    COMMON BLOCKS : none
; 
;          EXAMPLE : stereo_pointing
;
;    SPECIAL CALLS : NOTE: USES CURL TO DOWNLOAD THE FILE. Also uses get_stereo_browse,
;                    stereo_pointing_display, and several VERY useful WCS conversion
;                    routines.
;
;          WRITTEN : HPW. 10-MAR-2010. Appears to work!
;                    HPW. 24-AUG-2010. Removed call to eit_colors. Added ability to read
;                                      AIA data.
;                    HPW. 23-AUG-2010. Fixed bug with OBSRVTRY tag (from Jian Sun)
;                    HPW. 21-FEB-2012. lengend -> ssw_legend
;                    PRY. 03-MAR-2012. legend -> ssw_legend  :-)
;                    HPW. 05-JUN-2013. Fixed longitude bug so that B images project forward. 
; 
;-

;; --- svn check 10-FEB-2012

;; #######################################################################################
;;
;; Read the fits file and store the data and the header
;;

pro wt_read_file,WidgetData

  ;; --- read the file

  file = WidgetData.file

  head = headfits(file)
  INSTRUME = strtrim(fxpar(head,'INSTRUME'),2)
  OBSRVTRY = strtrim(fxpar(head,'OBSRVTRY'),2)
  if strpos(INSTRUME,'AIA') ne -1 then INSTRUME = 'AIA'

  case INSTRUME of
    'SECCHI': begin
      box_message,'READING STEREO FILE'
      Sdata = sccreadfits(file,Shead,/silent)
    end
    'AIA' : begin
      box_message,'READING AIA FILE'
      mreadfits,file,Shead,Sdata
      Shead = add_tag(Shead,'AIA','OBSRVTRY')
    end
    else : begin
      box_message,'READING UNKNOWN FILE'
      mreadfits,file,Shead,Sdata
      if not(tag_exist(Shead,'OBSRVTRY')) then Shead = add_tag(Shead,INSTRUME,'OBSRVTRY')
    end
  endcase

  print,'INSTRUME = '+Shead.INSTRUME
  print,'OBSRVTRY = '+Shead.OBSRVTRY

  wcs = fitshead2wcs(Shead)

  ;; --- some browse data has negative intensties

  image = sigrange(Sdata,range=range,f=0.99)
  if range[1] lt 0 then begin
    message,'reversing intensities',/info
    Sdata = abs(Sdata)
  endif

  WidgetData.Sdata = ptr_new(Sdata)
  WidgetData.Shead = ptr_new(Shead)
  WidgetData.wcs   = ptr_new(wcs)

  ;; --- there is no point

  WidgetData.point_defined = 0

  ;; --- erase some of the displays

  wset,WidgetData.windowIDs[1]
  erase
  wset,WidgetData.windowIDs[2]
  erase

  wset,WidgetData.windowIDs[0]
  
return
end

;; #######################################################################################
;;
;; Convert from device coordinates to data coordinats for the Sun display. Use saved
;; coordinates because other plots may have been made since the Sun was displayed.
;;

function wt_convert_coord,WidgetData,x0,y0,window=window

  Xs = !x
  Ys = !y
  Zs = !z
  Ps = !P

  if not(keyword_set(window)) then window = 0  
  draw = WidgetData.draw_coord[window]

  !x = *draw.X
  !y = *draw.Y
  !z = *draw.Z
  !p = *draw.P
  pt = convert_coord([x0,y0],/device,/to_data)

  !x = Xs
  !y = Ys
  !z = Zs
  !p = Ps

return,pt[0:1]
end

;; #######################################################################################
;;
;; plot a series of Stonyhurst grids showing the Sun at different times during the disk
;; passage. Dates and pointings are displayed
;;

pro wt_plot_grid,WidgetData,ps=ps,pos=pos

  if not(keyword_set(ps)) then wset,WidgetData.windowIDs[2]
  sp_setup_dots,size=1
  linecolors
  !p.multi=-1
  !p.position=-1

  lon0     = WidgetData.lon0       ;; these variabiles correspond to the STEREO pointing
  lat0     = WidgetData.lat0
  tim0     = WidgetData.tim0
  lon_list = *WidgetData.lon_list  ;; these variables correspond to the disk passage
  tim_list = *WidgetData.tim_list
  nlist    = n_elements(lon_list) 

  if keyword_set(ps) then begin
    noerase = 1 
    color   = 256
    cs      = 1.0
  endif else begin
    noerase  = 0
    !p.multi = [0,nlist,1,0,0]
    color    = 0
    cs       = 1.25
  endelse

  for i=0,nlist-1 do begin

    if keyword_set(ps) then begin
      x1 = i*0.2
      x2 = (i+1)*0.2
      y1 = pos[1]
      y2 = pos[3]
      !p.position = [x1,y1,x2,y2]
    endif

    plot,[-1200,1200],[-1200,1200],/nodata,noerase=noerase,xstyle=4,ystyle=4
    plot_helio,tim_list[i],/limb,grid=20,/over

    a = 60*hel2arcmin(lat0,lon_list[i],date=tim_list[i])
    plots,a[0],a[1],psym=8,color=2

    title = anytim(tim_list[i],/vms,/trunc)
    if keyword_set(ps) then y0 = -1350 else y0 = -1550
    xyouts,0,y0,title,align=0.5,charsize=cs

    title = trim(a[0],'(f10.1)')+','+trim(a[1],'(f10.1)')
    if keyword_set(ps) then y0 = 1100 else y0 = 1300
    xyouts,0,y0,title,align=0.5,charsize=cs

  endfor

  !p.multi    = -1
  !p.region   = -1
  !p.position = -1

return
end

;; #######################################################################################
;; 
;; --- display a Stonyhurst plot of the Sun with the disk passage track outlined on it
;;

pro wt_plot_sun,WidgetData,plot_time=plot_time,ps=ps,pos=pos

  if not(keyword_set(ps)) then begin
    wset,WidgetData.windowIDs[1]
    erase
  endif
  if keyword_set(ps) then size = 1.0 else size = 1.25
  sp_setup_dots,size=1.5
  !p.multi=-1
  linecolors

  lon0     = WidgetData.lon0         ;; these variabiles correspond to the STEREO pointing
  lat0     = WidgetData.lat0
  tim0     = WidgetData.tim0 
  lon_list = *WidgetData.lon_list    ;; these variables correspond to the disk passage
  tim_list = *WidgetData.tim_list
  nlist    = n_elements(lon_list)

  ;; --- time interval for disk passage
  time_start = anytim(tim_list[0])
  time_end   = anytim(tim_list[nlist-1])
  dt_all     = time_end-time_start

  ;; --- time of interest based on slider or field
  if keyword_set(plot_time) then begin
    this_time = anytim(plot_time)
    if this_time gt time_start and this_time le time_end then begin
      f = (this_time - time_start)/dt_all
      WIDGET_CONTROL,WidgetData.sliderIDs[0],set_value=f
    endif else begin
      message,'INPUT TIME OUTSIDE OF DISK PASSAGE',/info
      return
    endelse
  endif else begin
    WIDGET_CONTROL,WidgetData.sliderIDs[0],get_value=f
    this_time = time_start + f*dt_all
    WIDGET_CONTROL,WidgetData.fieldIDs[0],set_value=anytim(this_time,/vms,/trunc)
  endelse
  
  ;; --- plot a track across the Sun
  if keyword_set(pos) then !p.position = pos
  if keyword_set(ps) then color = 255 else color = 0
  plot,[-1200,1200],[-1200,1200],/nodata,/noerase,xstyle=4,ystyle=4,$
      xtickname=replicate(' ',20),ytickname=replicate(' ',20)
  plot_helio,anytim(this_time,/ccsds),/over,/limb

  ngrid = 20
  x = fltarr(ngrid)
  y = fltarr(ngrid)
  for i=0,ngrid-1 do begin

    ;; --- fraction of the time interval
    frac = i/float(ngrid-1)
    dt   = frac*dt_all/86400.
    t    = anytim(time_start + dt*86400.,/ccsds)

    ;; --- time since STEREO point
    dt_s = (anytim(t) - anytim(tim0))/86400.
    lon  = lon0 + diff_rot(dt_s,lat0,synodic=1,sidereal=0,rigid=0)

    ;; --- convert to arcsec
    a    = 60*hel2arcmin(lat0,lon,date=t)
    x[i] = a[0]
    y[i] = a[1]

  endfor
  if keyword_set(ps) then color = 8 else color = 5
  if keyword_set(ps) then thick = 6 else thick = 2
  plots,x,y,color=color,thick=thick

  ;; --- plot the point of interest

  ;; --- fraction of the time interval
  dt   = f*dt_all/86400.
  t    = anytim(time_start + dt*86400.,/ccsds)

  ;; --- time since the STEREO point
  dt_s = (anytim(t) - anytim(tim0))/86400.
  lon  = lon0 + diff_rot(dt_s,lat0,synodic=1,sidereal=0,rigid=0)

  a    = 60*hel2arcmin(lat0,lon,date=t)
  r    = sqrt(a[0]^2+a[1]^2)
  plots,a[0],a[1],psym=8,/data,color=2

  if keyword_set(ps) then cs = 1.0 else cs = 1.25

  title = anytim(t,/vms,/trunc)
  if keyword_set(ps) then y0 = -1250 else y0 = -1200
  xyouts,0,y0,title,align=0.5,charsize=cs
  
  title = trim(a[0],'(f10.1)')+','+trim(a[1],'(f10.1)')
  if keyword_set(ps) then y0 = 1100 else y0 = 1200
  xyouts,0,y0,title,align=0.5,charsize=cs

  !p.position = -1

  wait,0.1

return
end

;; #######################################################################################
;;
;; --- compute the disk passage of the point selected in the STEREO image
;;

pro wt_compute_hg,WidgetData

  wcs   = *WidgetData.wcs
  xpix  = WidgetData.point[0]
  ypix  = WidgetData.point[1]
  Shead = *WidgetData.Shead

  ;; --- peform the conversion from wcs to heliographic coordinates
  coord = wcs_get_coord(wcs,[xpix,ypix])
  wcs_convert_from_coord,wcs,coord,'HG',lon0,lat0

  ;; --- map 180 degrees to -180 degrees
  WIDGET_CONTROL,WidgetData.bgroupIDs[0],get_value=v
  this_rotation = v[0]
  if this_rotation eq 0 then begin
    if lon0 gt 0 then lon0 = lon0 - 360
  endif

  ;; --- compute the rotation of the selected position for the next 28 days
  time_start    = anytim(Shead.date_obs)
  time_end      = time_start + 28*86400.  
  time_grid_abs = anytim(timegrid(time_start,time_end,hours=6))
  time_grid_rel = time_grid_abs - time_grid_abs[0]
  ntime         = n_elements(time_grid_abs) 

  lon = fltarr(ntime)
  for i=0,ntime-1 do begin
    dt_s   = time_grid_rel[i]/86400.
    lon[i] = lon0 + diff_rot(dt_s,lat0,synodic=1,sidereal=0,rigid=0)
  endfor

  ;;plot,time_grid_rel,lon

  ;; --- find times of interest for disk observations
  lon_list = [-90,-45,0,45,90]
  nlist    = n_elements(lon_list) 
  tim_list = strarr(nlist)
  openw,unit,'stereo_pointing.txt',/get
  printf,unit,'# '
  printf,unit,'# POINT IN STEREO IMAGE '
  printf,unit,'# TIME          = '+Shead.date_obs
  printf,unit,'# X (PIXELS)    = '+trim(xpix)
  printf,unit,'# Y (PIXELS)    = '+trim(ypix)
  printf,unit,'# lon (degrees) = '+trim(lon0,'(f10.1)')
  printf,unit,'# lat (degrees) = '+trim(lat0,'(f10.1)')
  printf,unit,'#','LON','DATE','X"','Y"','R','Rsun',format='(a1,a4,a30,4a10)'
  print,'#','LON','DATE','X"','Y"','R','Rsun',format='(a1,a4,a30,4a10)'
  for i=0,nlist-1 do begin
    this_time   = interpol(time_grid_rel,lon,lon_list[i])
    tim_list[i] = anytim(time_grid_abs[0]+this_time,/ccsds)
    a  = 60.0*hel2arcmin(lat0,lon_list[i],date=tim_list[i])
    r  = sqrt(a[0]^2+a[1]^2)
    Rsun = (get_sun(tim_list[i]))[1]
    print,lon_list[i],anytim(tim_list[i],/vms,/trunc),a[0],a[1],r,Rsun,$
        format='(i5,a30,4f10.1)'
    printf,unit,lon_list[i],anytim(tim_list[i],/vms,/trunc),a[0],a[1],r,Rsun,$
        format='(i5,a30,4f10.1)'
  endfor
  free_lun,unit

  WidgetData.lon0     = lon0
  WidgetData.lat0     = lat0
  WidgetData.tim0     = anytim(Shead.date_obs,/ccsds)
  WidgetData.lon_list = ptr_new(lon_list)
  WidgetData.tim_list = ptr_new(tim_list)

  WIDGET_CONTROL,WidgetData.fieldIDs[0],set_value = anytim(tim_list[0],/vms,/trunc)

return
end

;; #######################################################################################
;;
;; --- Display the selected EUVI Image

pro wt_display_euvi,WidgetData,ps=ps,pos=pos

  if not(keyword_set(ps)) then wset,WidgetData.windowIDs[0]
  if not(keyword_set(ps)) then pos = sp_plot_pos(1,1,/nolabels)

  Sdata = *WidgetData.Sdata
  Sdata = float(Sdata)
  Shead = *WidgetData.Shead
  wave  = Shead.WAVELNTH

  ;; --- scale the image
  image = sigrange(Sdata,range=range,f=0.999)
  Sdata = Sdata - range[0]
  image = sigrange(Sdata,range=range,f=0.999)
  range = range[1]*[1.0E-3,1.0]
  image = eis_scale_image(Sdata,'logarithmic',range)
  IMGnx = Shead.naxis1
  IMGny = Shead.naxis2

  ;; --- display the image
  ;;eit_colors,wave
  loadct,0,/silent
  stereo_pointing_display,[0,IMGnx-1],[0,IMGny-1],image,/nox,/noy,pos=pos,/black

  ;; --- the limb as seen from Earth
  linecolors
  if keyword_set(ps) then size = 0.25 else size = 0.25
  sp_setup_dots,size=size
  if keyword_set(ps) then symsize = 10 else symsize = 1
  color = 255
  nlat = 45 
  latL = -90 + 180*findgen(nlat)/(nlat-1)
  wcs  = *WidgetData.wcs
  for k=1,2 do begin
    if k eq 1 then lonL = 90 else lonL = -90
    for i=0,nlat-1 do begin
      wcs_convert_to_coord,wcs,coord,'hg',lonL,latL[i]
      pix = wcs_get_pixel(wcs,coord)
      plots,pix[0],pix[1],psym=8,color=color,/data
    endfor
  endfor

  ;; --- the selected point
  if WidgetData.point_defined eq 1 then begin
    if keyword_set(ps) then size = 1.0 else size = 1.25
    sp_setup_dots,size=size
    plots,WidgetData.point[0],WidgetData.point[1],psym=8,color=2
  endif

  title = str_replace(Shead.OBSRVTRY,'_',' ')
  title = title+' '+trim(wave)+' '+anytim(Shead.date_obs,/vms,/trunc)
  if keyword_set(ps) then cs = 1.0 else cs = 1.5
  ssw_legend,title,box=0,/bottom,charsize=xs,textcolor=255

  if keyword_set(ps) then color = 255 else color = 0
  ssw_legend,'Limb as Seen From Earth',box=0,/top,charsize=xs,textcolor=255,linestyle=1,color=255


  if not(keyword_set(ps)) then begin
    WidgetData.draw_coord[0].X = ptr_new(!x)
    WidgetData.draw_coord[0].Y = ptr_new(!y)
    WidgetData.draw_coord[0].Z = ptr_new(!z)
    WidgetData.draw_coord[0].P = ptr_new(!p)
  endif

return
end

;; #######################################################################################

pro stereo_pointing_event,event

  ;; --- get information on the event
  WIDGET_CONTROL,event.id,GET_UVALUE=event_type

  ;; --- get WidgetData
  BaseWidget = event.handler
  WIDGET_CONTROL,BaseWidget,GET_UVALUE=unseen
  WidgetData = *unseen
  
  case event_type of 
    
    "QUIT": begin
      WIDGET_CONTROL,/DESTROY,event.top
    end

    "DATA1": begin
      file = get_stereo_browse(group_leader=WidgetData.base)
      if file_exist(file) then begin
        WidgetData.file = file
        wt_read_file,WidgetData
        wt_display_euvi,WidgetData
      endif
    end

    "DATA2": begin
      file = dialog_pickfile(/read,filter=['*.fits','*.fts','*.fits.gz'])
      if file_exist(file) then begin
        WidgetData.file = file
        wt_read_file,WidgetData
        wt_display_euvi,WidgetData
      endif
    end

    "DRAW1": begin
      if event.release eq 1 then begin
        point = wt_convert_coord(WidgetData,event.x,event.y,window=0)
        WidgetData.point = point
        WidgetData.point_defined = 1
        wt_display_euvi,WidgetData
        wt_compute_hg,WidgetData
        wt_plot_grid,WidgetData
        WIDGET_CONTROL,WidgetData.sliderIDs[0],set_value=0
        wt_plot_sun,WidgetData
      endif
    end

    "PLOT": begin
      if WidgetData.point_defined then begin
        WIDGET_CONTROL,WidgetData.fieldIDs[0],get_value=t
        t = anytim(t,/ccsds)
        wt_plot_sun,WidgetData,plot_time=t
      endif
    end

    "REPLOT": begin
      if WidgetData.point_defined eq 1 then begin
        wt_display_euvi,WidgetData
        wt_compute_hg,WidgetData
        wt_plot_grid,WidgetData
        WIDGET_CONTROL,WidgetData.sliderIDs[0],set_value=0
        wt_plot_sun,WidgetData
      endif
    end

    "SLIDER": begin
      if WidgetData.point_defined then begin
        wt_plot_sun,WidgetData
      endif
    end

    "PS": begin
      if WidgetData.point_defined then begin

        break_file,WidgetData.file,disk,dir,name,ext
        TOTALx = 1024
        TOTALy = 512+205
        PSw    = 10.0
        scale  = PSw/TOTALx
        PSh    = scale*TOTALy
        sp_setup_ps,/ps,w=PSw,h=PSh,/color,/land,font=10,file=name

        pos = [0,0,1024,205]/float([TOTALx,TOTALy,TOTALx,TOTALy])
        wt_plot_grid,WidgetData,/ps,pos=pos

        pos = [512,200,1024,717]/float([TOTALx,TOTALy,TOTALx,TOTALy])
        wt_plot_sun,WidgetData,/ps,pos=pos

        pos = [0,200,512,717]/float([TOTALx,TOTALy,TOTALx,TOTALy])
        wt_display_euvi,WidgetData,/ps,pos=pos

        sp_setup_ps,/ps,/close
        !p.region = -1
        !p.position = -1
        !p.multi = -1
      endif
    end

    ELSE: begin
      (*unseen) = WidgetData
      return
    end
    
  endcase
  
  (*unseen) = WidgetData
  
return
end

;; #######################################################################################

pro stereo_pointing

  BASE = WIDGET_BASE(/ROW,TITLE='STEREO Pointing',/FRAME)

  BASE1 = WIDGET_BASE(BASE,/COLUMN,/FRAME)

    BUTTON  = WIDGET_BUTTON(BASE1,VALUE='   Quit   ',UVALUE='QUIT')
    BUTTON  = WIDGET_BUTTON(BASE1,VALUE='  Get Stereo Browse  ',UVALUE='DATA1')
    BUTTON  = WIDGET_BUTTON(BASE1,VALUE=' Read My Fits File ',UVALUE='DATA2')
    SLIDER1 = CW_FSLIDER(BASE1,MIN=0,MAX=1.0,UVALUE='SLIDER',/SUPPRESS_VALUE,/DRAG,$
                         title='Rotate Region',xsize=160)
    FIELD1  = CW_FIELD(BASE1,title='Observation Time',/COLUMN,xsize=24)
    BUTTON  = WIDGET_BUTTON(BASE1,VALUE=' Plot Observation Time ',UVALUE='PLOT')

    BUTTON  = WIDGET_BUTTON(BASE1,VALUE=' Print to PS file ',UVALUE='PS')

    LABEL   = WIDGET_LABEL(BASE1,value='')
    BGROUP1 = CW_BGROUP(BASE1,'Compute for this rotation',/nonexclusive,UVALUE='NULL')
    BUTTON  = WIDGET_BUTTON(BASE1,VALUE=' Re-Plot ',UVALUE='REPLOT')

  BASE2 = WIDGET_BASE(BASE,/COLUMN)

  BASE2a = WIDGET_BASE(BASE2,/ROW,xpad=0,ypad=0)
  
    DRAW1 = WIDGET_DRAW(BASE2a,XSIZE=512,YSIZE=512,/MOTION,/BUTTON,$
                        UVALUE='DRAW1')
    DRAW2 = WIDGET_DRAW(BASE2a,XSIZE=512,YSIZE=512,UVALUE='DRAW2')

    DRAW3 = WIDGET_DRAW(BASE2,XSIZE=1027,YSIZE=200)

  WIDGET_CONTROL,BASE,/REALIZE

  ;; -------------------------------------------------------------------------------------
  ;; --- WidgetData ----------------------------------------------------------------------
  ;; -------------------------------------------------------------------------------------

  ;; --- get window IDs
  WIDGET_CONTROL,DRAW1,GET_VALUE=window1
  WIDGET_CONTROL,DRAW2,GET_VALUE=window2
  WIDGET_CONTROL,DRAW3,GET_VALUE=window3
  
  windowIDs = [window1,window2,window3]
  sliderIDs = [SLIDER1]
  fieldIDs  = [FIELD1]
  bgroupIDs = [BGROUP1]

  draw_coord = {X:  ptr_new(),$
                Y:  ptr_new(),$
                Z:  ptr_new(),$
                P:  ptr_new()}
  draw_coord = replicate(draw_coord,n_elements(windowIDs))
  
  ;; --- structure to hold program info 
  WidgetData = {base:           BASE, $
                windowIDs:      windowIDs,$
                sliderIDs:      sliderIDs,$
                fieldIDs:       fieldIDs,$
                bgroupIDs:      bgroupIDs,$
                file:           '',$
                Sdata:          ptr_new(),$  ;; the STEREO data
                Shead:          ptr_new(),$  ;; the STEREO index
                wcs:            ptr_new(),$  ;; the wcs structure for the data
                draw_coord:     draw_coord,$ ;; saved coordinates for device conversion
                point_defined:  0,$          
                point:          [0.0,0.0],$  ;; stereo point in stereo pixel coordinates
                lon0:           0.0,$        ;; stereo point in helographic coordinates
                lat0:           0.0,$        ;; stereo point in helographic coordinates
                tim0:           '',$         ;; stereo point in helographic coordinates
                lon_list:       ptr_new(),$  ;; longitudes of interest
                tim_list:       ptr_new()}   ;; times for these longitudes

  ;; --- store WidgetData as pointer
  unseen = ptr_new(WidgetData)
  WIDGET_CONTROL,BASE,SET_UVALUE=unseen

  ;; -------------------------------------------------------------------------------------
  ;; --- Register widget with XMANAGER ---------------------------------------------------
  ;; -------------------------------------------------------------------------------------

  XMANAGER,'stereo_pointing',BASE
  
  WidgetData = *unseen

return
end