;+
; NAME:
;     ARLOCPLOT
; PURPOSE:
;     Plot active region and phase center locations on the solar disk
;     for a solar scan, and allow readout of coordinates using the mouse.
; CATEGORY:
;     OVRO APC DATA ANALYSIS
; CALLING SEQUENCE:
;     arlocplot,filename,hrec
; INPUTS:
;     filename   the name of the file containing the data.
;     hrec       the record number of the scan header for a solar scan.
; OPTIONAL (KEYWORD) INPUT PARAMETERS:
; ROUTINES CALLED:
;     openarc, newscan, freescan, ovsa_get_index, findseg, getdata, decode,
;     cvdoy, cnv_xy2h, cnv_h2xy, stonyhurst
; OUTPUTS:
; COMMENTS:
;     Crude exit indicator: Click on left side within 10 pixels of edge
;     of window to exit.
; SIDE EFFECTS:
; RESTRICTIONS:
; MODIFICATION HISTORY:
;     Written 31-May-2002 by Dale Gary
;     08-Jun-2002  DG
;       Cleaned up coordinate display, and made it exit gracefully.
;     11-Jan-2004  DG
;       Added plot of position angle of fringes.  Currently only does 12
;       baseline, but easily expanded to do any baseline--I just need a
;       way for the user to specify the desired baseline.
;-
pro arlocplot,filename,hrec

   color1 = 255
   color2 = 200

   lun = openarc(filename,a,nrec)
   if (lun eq 0) then begin
      print,'ARLOCPLOT: Error opening file ',filename
      return
   endif

   ; Get index segment
   idx = ovsa_get_index(a,nrec,/create)

   ; Read header info, etc.
   rec = newscan(a,nrec,hrec,header,obseq,cfg,traj,refcal,gparm,geometry,pRecent,pOld,cyc_str) - 1

   ; Calculate uv for a set of times during the day (will be used to plot
   ; the position angles of the fringes for the different baselines)
   uv = dblarr(header.nant,header.nant,9)  ; do for 9 times (one per hour)
   for i = 0, 8 do begin   ; Loop over hour

      time = hms2sec(160000L+i*10000L)*1000.D       ; Ranges from 16 UT to 24 UT
      tnow = time + 64000.D  ; Add approx ET - UT1

      ; Calculate the RA, Dec, and HA at the current time, in absence of
      ; horizontal parallax.  Multiply powers of tnow by fit coefficients and
      ; use TOTAL() function to obtain sum.

      tarray = [1D,tnow,tnow^2]
      ranow = total(geometry.a_ra*tarray)

      decnow = total(geometry.a_dec*tarray)

      ; Convert date and time to single MJD double-precision number, for use in
      ; OVSA_LST() routine.
      mjd = cvdoy(geometry.year,geometry.doy[0],/JUL) + (time+geometry.dut)/8.64D7
      hnow = ovsa_lst(mjd) - ranow

      ; Time delay [msec] relative to current time
      taunow = bdots(hnow,decnow,geometry,uvnow)
      uv[*,*,i] = uvnow
   endfor

   ; Get MJD from TLS of header
   mon = strmid(header.tls.date,5,2)
   day = strmid(header.tls.date,8,2)
   yr = header.tls.year
   hr = strmid(header.tls.timstr,1,2)
   mn = strmid(header.tls.timstr,4,2)
   sc = strmid(header.tls.timstr,7,2)
   mjd = julday(mon,day,yr,hr,mn,sc) - 2400000D

   ; Read the EPHEMERIS segment from the file
   rec = findseg(a,hrec,nrec,!SEGM.EPHEM,idx=idx)
   if (rec le 0) then begin
      print,'ARLOCPLOT: Error reading EPHEM segment from file ',filename
      free_lun,lun
      return
   endif else begin
      data = getdata(rec-1,a)
      ephem = decode(data,!SEGM.EPHEM)
   endelse

   tls = tl_decode(data)

   ; Extract the REGION structure from the EPHEM structure
   regn = ephem.regn
   nregn = regn.nregn

   ; Determine difference in time between EPHEM time and current scan time
   res = cvdoy(ephem.regn.epoch.yr,ephem.regn.epoch.doy)
   hour = fix(ephem.regn.epoch.sec/3600.)
   regn_mjd = julday(res[1],res[2],res[0],24,0,0) - 2400000D
   now_mjd = double(mjd)

   ; plot the grid lines and active regions
   ; Determine the P-angle, B0-angle, and the solar radius in arcsec for
   ; date of EPHEM.REC
 ;  rotangle = 0  ; Equatorial rotation rate (degrees)
   rotnow = 13.83*(regn_mjd-now_mjd)
   pandb,regn_mjd,p,b,r
   loadct,3
   out = [256,256,200]
   nx = 512
   ny = 512
   window,/free,xsiz=nx,ysiz=ny
   wid = !d.window

   ; Plot a circle representing the disk
   th = findgen(361)*!dtor
   plots,/dev,out[2]*cos(th)+out[0],out[2]*sin(th)+out[1],color=200

   ; Overplot the solar grid lines
   psty = !p.linestyle
   pcol = !p.color
   !p.color = color1
   !p.linestyle = 1
   stonyhurst,b,p,out[2],center=[out[0],out[1]],deg=15
   !p.linestyle = psty
   !p.color = pcol

   for i = 0, nregn-1 do begin
      arstr = string(regn.info(i).noa,format='("AR",I4)')
      lat = regn.info(i).lat
      lng = regn.info(i).lng

      ; Get new longitude, correcting rotation rate for dependence
      ; on latitude (from Newton-Nunn formula) [NB: -rotangle because
      ; of definition of longitude as +E]
;      lng = lng - rotangle*(1-0.2058*(sin(lat*!dtor)^2))
      ; Advance to 2000 UT today
      lngnow = lng + rotnow*(1-0.2058*(sin(lat*!dtor)^2))

      ; Repeat for current region location
      cnv_h2xy,p,b,r,lat,-lngnow,x,y
      ; Convert offsets in arcsec to device coords
      xdev1 = (x/r)*out[2]+out[0]
      ydev1 = (y/r)*out[2]+out[1]
      plots,/dev,xdev1,ydev1,psym=6
      xyouts,/dev,xdev1,ydev1+10,arstr,align=0.5
;      plots,/dev,[xdev,xdev1],[ydev,ydev1],linesty=2

   endfor

   ; Overplot the phase center
   hao = (*traj.pha)[0]/1000.
   deco = (*traj.pdec)[0]/1000.
   plots,hao*(3600/r)*out[2] + out[1],deco*(3600/r)*out[2] + out[1],/dev,psym=2,symsize=2

   ; Overplot HAO and DECO lines
   psty = !p.linestyle
   pcol = !p.color
   !p.color = color2
   !p.linestyle = 1
   for i = -3,3 do begin
      plots,[i,i]*(360./r)*out[2] + out[0],[0,ny],/dev
      plots,[0,nx],[i,i]*(360./r)*out(2) + out[1],/dev
   endfor

   ; Determine the position angles of baselines, and plot them in multiple
   ; pixmap windows
   kwin = 4
   nant = header.nant
   ; Loop over all baselines
   for iant = 0, nant-2 do begin
      for jant = iant+1, nant-1 do begin
         kwin = kwin+1
         window,kwin,/pixmap,xsiz=nx,ysiz=ny
         wset,kwin
         device,copy=[0,0,nx,ny,0,0,wid]
         pa = reform(atan(uv[iant,jant,*],uv[jant,iant,*]))  ; PA of fringes
         x = out[2]*[1.025,0.975]#cos(pa)
         y = out[2]*[1.025,0.975]#sin(pa)
         !p.linestyle=0
         txt = ['16','17','18','19','20','21','22','23','24']
         for i = 0, 8 do begin
            plots,x[*,i]+out[0],y[*,i]+out[1],/dev
            xyouts,x[0,i]+out[0],y[0,i]+out[1],/dev,txt[i],align=0,color=255
         endfor
         for i = 0, 8 do begin
            plots,-x[*,i]+out[0],-y[*,i]+out[1],/dev
            xyouts,-x[0,i]+out[0],-y[0,i]+out[1],/dev,txt[i],align=1,color=255
         endfor
      endfor
   endfor
   wset,wid

   !p.linestyle = psty
   !p.color = pcol
   free_lun,lun
   ; Free memory associated with structures
   freescan,obseq,traj,refcal,gparm,geometry,pRecent,pOld,cyc_str

   pchar = !p.charsize
   !p.charsize=1.5
   ; Make exit "button"
   polyfill,/dev,[0,50,50,0],[0,0,30,30]
   xyouts,/dev,25,5,'!6Exit',align=0.5,color=0,charthick=2
   for i = 0, kwin-5 do begin
      polyfill,/dev,[0,19,19,0]+i*20+60,[1,1,20,20]
      xyouts,/dev,10+i*20+60,5,align=0.5,color=0,charthick=2,string(i,format='(I2)')
   endfor

   xd = 100
   while (xd gt 10) do begin
      cursor,xd,yd,/dev,/down
      if (!d.window ne wid) then begin
         !p.charsize = pchar
         for k = 5, kwin do wdelete,k  ; Delete all the pixmap windows
         return  ; Window was closed by user
      endif
      if (xd le 50 and yd le 30) then begin
         wdelete,wid
         for k = 5, kwin do wdelete,k  ; Delete all the pixmap windows
         !p.charsize = pchar
         return
      endif

      if (xd gt 60 and yd le 20) then begin
         device,copy=[0,0,nx,ny,0,0,5+(xd-60)/20]
         polyfill,/dev,[0,50,50,0],[0,0,30,30]
         xyouts,/dev,25,5,'!6Exit',align=0.5,color=0,charthick=2
         for i = 0, kwin-5 do begin
            polyfill,/dev,[0,19,19,0]+i*20+60,[1,1,20,20]
            xyouts,/dev,10+i*20+60,5,align=0.5,color=0,charthick=2,string(i,format='(I2)')
         endfor
      endif else begin

         ha = (xd-out[0])*(r/3600)/out[2]
         dec = (yd-out[1])*(r/3600)/out[2]

         x = (xd-out[0])*r/out[2]
         y = (yd-out[1])*r/out[2]
         cnv_xy2h,p,b,r,lat,lng,x,y
         if (lat gt 0) then slat = 'N ' else slat = 'S '
         if (lng gt 0) then slng = 'W' else slng = 'E'

         ; Write out results to window
         polyfill,/dev,[0,250,250,0],[ny-60,ny-60,ny-1,ny-1],color=color2
         xyouts,/dev,5,ny-15,'!6'+string(ha,dec,format='("Geo: ",2(f7.3,"d "))'),color=0
         xyouts,/dev,5,ny-35,'Helio: '+string(abs(lat),format='(f5.2)')+slat+' '+string(abs(lng),format='(f5.2)')+slng+'!3',color=0
         xyouts,/dev,5,ny-55,'!7D!6r: '+string((ha-hao)*3600,(dec-deco)*3600, format='(2(f8.1,''"''))'),color=0
      endelse
   endwhile

return
end