pro IPSedit, reference=REFERENCE

@compile_opt.pro		; On error, return to caller

;+
; NAME:
;	IPSedit
; PURPOSE:
;	Display Cambridge IPS data (one day) in rectangular format.
;	Edit data using mouse
; CALLING SEQUENCE:
;	IPSedit, reference=REFERENCE
; OPTIONAL INPUTS: (only in X-windows; otherwise REFERENCE=0 is used)
;	REFERENCE	= 0 (default), displays only one dayly G-map
;			= 1 displays 1 extra (i.e. total of two) G-maps
;			    (user is prompted for both MJD and type of the G-map)
;			= 2 displays 2 extra (i.e. total of three) G-maps.
;			    (user is prompted only for MJD of the center day
;			    and data type of the first map; the other two will
;			    be for MJD-1 and MJD+1 and will be the same type.
;			= 3 displays 3 extra (i.e. total of four) G-maps.
;			    (user is prompted only for MJD and data type of the
;			    1st map; the other three will be for MJD+1,MJD+2,
;			    MJD+3 and will be the same type)
; CALLS:
;	call_external, ColorBox, Echo, GetColors, IPS_Gadgets,
;	IPS_OnGadgets, Pcursor, Read_GE, Read_GG, Julian, Sphere_Circle
;	twin, TimeSet, TimeGet, TimeUnit
; RESTRICTIONS:
;	Unedited data files are read using the Fortran subroutine Read_GG.
;	This module is called by call_external. The logical EXEIDL should be
;	assigned to the executable file for Read_GG. (The logical is best
;	defined in your login file; the executable can be created by the
;	command procedure X$COM:IDLEXE.COM).
; COMMON BLOCKS:
	common IPSedit_save, ilong, olong, mjd_save
; PROCEDURE:
;   >	User prompts:
;	- data directory (is used to set logical $DAT)
;	- modified Julian day
;	- edited or unedited data files
;   >   Edited data can be saved in new file
;   >	Unedited data files are identified by the name 'Gxxxxx.' (where xxxxx
;	stands for the 5-digit modified Julian day (these are the original
;	Cambridge files in hex format)
;   >	Edited data files are identified by the name 'Exxxxx.DAT', (xxxxx is
;	again the modified Julian day)
;   >	If more than one G-map is displayed, only the map in the window
;	title 'IPS edit' can be modified.
; MODIFICATION HISTORY:
;	SEP-1992, Tom Davidson, Paul Hick (UCSD)
;	JAN-1994, Paul Hick (UCSD); added the REFERENCE option to display
;		multiple files in separate windows (X-windows only).
;-

InitVar, ilong, 30.
InitVar, olong, 90.
InitVar, mjd_save, 48794l

InitVar, REFERENCE, 0
REFERENCE = (REFERENCE > 0) < 3

; The first call to twin will create the graphics window (TWIN window) in
; X-windows. For some reason this affects the setting of !d.n_colors (changes
; it from 256 to 246).					(PLH, 12/27/93)

if (!d.NAME eq 'X' or !d.NAME eq 'WIN') and REFERENCE ne 0 then begin
    if REFERENCE ge 2 then	$
        twin, /quarter, /hide, /bl, title='IPS reference'
    if REFERENCE ge 3 then	$
        twin, /quarter, /hide, /br, title='IPS reference'
        twin, /quarter, /hide, /tr, title='IPS reference'
	twin, /quarter, /hide, /tl, title='IPS edit'      , load = 5
endif else begin
    if REFERENCE then message, /info, 'Keyword REFERENCE only permitted in X-windows
    REFERENCE = 0
    twin, /hide, title='IPS edit'	, load = 5, aspectratio=1.2
endelse

Nt = 72				; # times
Xt = findgen(Nt+1)/3.		; edges of bins in hour angle
widXt = Xt[1]-Xt[0]		; bin width in hour angle

; Declination list (degrees)
Dec = [-7.8, 1.9, 9.8, 16.9, 23.4, 29.5,35.4,41.,46.6,52.1,57.7,63.2,68.9, 74.8]
TopDec = 82.			; Phony declination to mark top of plot
Nd  = 14			; # declinations
Yd  = [Dec,TopDec]			; edges of bins in declination
widYd = Yd[1:Nd]-Dec		; bin width in declination

mjd = mjd_save+indgen(REFERENCE+1)
xg  = reform( fltarr(Nt*Nd*(REFERENCE+1)) , Nt,Nd, REFERENCE+1 )
xgr = fltarr(Nt,Nd)		; Used to read the IPS files

Format = '(F8.1)'

ipsEdir = filepath(root=getenv('SSW_SMEI_DAT'),subdir='cambridge','')
echo, 'Data directory ?', ipsEdir
ipsEdir = replicate(ipsEdir,REFERENCE+1)

dum = ipsEdir[0]
if REFERENCE eq 1 then begin
	echo, 'Reference directory ?', dum
	ipsEdir[1] = dum
endif

echo, 'Inner Elongation', ilong, 0., 90.
echo, 'Outer Elongation', olong, 0., 90.

Qed = 1B+bytarr(REFERENCE+1)  &  iok = 0

GIF = 0B	; &  echo,'GIF Output?',GIF,/log

PROMPTS = ['editing','reference']

START:

if REFERENCE eq 1 then irmax = 1 else irmax = 0

for ir=0,irmax do begin
	dum = mjd[0]  &  echo, 'MJD for '+PROMPTS[ir]+' (-1 to Quit)', dum
	if dum lt 0 then begin
		twin, /delete, /all		; Delete all windows
		return
	endif
	mjd[ir] = dum
	dum = Qed[0]  &  echo, 'Edited data ?', dum, /logical  &  Qed[ir] = dum
endfor

if REFERENCE ge 2 then begin
	Qed[*] = Qed[0]
	if REFERENCE eq 2 then begin
		mjd[1] = mjd[0]-1  &  mjd[2] = mjd[0]+1
	endif else	$
		mjd = mjd[0]+indgen(REFERENCE+1)
endif

for ir=0,REFERENCE do begin

	case Qed[ir] of
	0: begin			; Read unedited file
		status = call_external('EXEIDL','read_gg',mjd[0],iok,xgr)  &  print
		file = ipsEdir[ir]+'G'+strcompress(mjd[ir],/remove_all)+'.'
	end
	1: begin			; Read edited data file
		read_ge, mjd[ir], iok, xgr, dir=ipsEdir[ir]
		file = filepath(root=ipsEdir[ir],'e'+strcompress(mjd[ir],/remove_all)+'.dat')
	end
	endcase

	if not iok then begin
		message, 'Data file '+file+' not read!', /info
		if ir eq 0 then goto, START		; 1st file not read: cancel
		file = ''
		xg[*,*,ir] = -1.
	endif else begin
		if ir eq 0 then mjd_save = mjd[0]	; Store in common block IPSedit_save
		file = (file_search(file))[0]		; Full file name
		xg[*,*,ir] = xgr
	endelse

	; Convert modified julian day to actual JD
	; Calculate right ascension, declination, and sidereal time of sun

	T = TimeSet(mjd1=mjd[ir]+0.5)

	gst   = TimeGST(T)
	tmp   = NewcombSun(T, /long, /lat, /degrees)
	slong = tmp[0]
	tmp   = CvSky(T, fromecliptic=tmp, /toequatorial, /degrees)
	srasn = tmp[0]
	sdec  = tmp[1]

	ha = gst-srasn			; ha is hour angle of Sun in degrees
	ha = ha+(ha lt 0)*360.

	case ir eq 0 of
	0: twin, /show, tr=ir eq 1, bl=ir eq 2, br=ir eq 3
	1: if REFERENCE ne 0 then twin, /show, /tl else twin, 1
	endcase

	erase

	POS = [.15,.1,.93,.95]		; set up the proper scaling
	plot, Xt, Yd,	/nodata, /noerase, xrange=[Xt(Nt),Xt(0)],	$
		ystyle=5,xstyle=5,position=POS

	if GIF then !X.S[1] = !X.S[1]*.852944	; Calibration for GIF output

	BreakVal = 0.+.25*16/!d.n_colors*findgen(!d.n_colors-1)
	if iok then	ColorBox , Xt[0:Nt-1] , Yd[0:Nd-1] , widXt, widYd ,	$
		GetColors ( xg[*,*,ir] , Breakval , /legend , format=Format, load=5 ) $
	else xyouts, 18., 70.,charsize=2.,	$
		'NO DATA FOR MJD '+strcompress(mjd[ir],/remove_all)

	plot, Xt, Yd,	$	; Draw the axes surrounding the color plot
		title = file,			$
		xtitle= 'Hour',			$
		ytitle= 'Declination (DEG)',	$
		/nodata, /noerase,		$
		xrange= !x.crange,		$
		xticks= 6,			$
		xstyle= 1,  ystyle=1,		$
		yticks= 14, ytickv=dec,		$
		position= POS

; Draw max and min elongation circles

	SPHERE_CIRCLE,latitude=sdec,elongation=ilong,dlon,dlat,angmin=-180,angmax=180
	oplot, (dlon+ha)/15., dlat

	SPHERE_CIRCLE,latitude=sdec,elongation=olong,dlon,dlat,angmin=-179,angmax=-1
	oplot, (dlon+ha)/15., dlat

	SPHERE_CIRCLE,latitude=sdec,elongation=olong,dlon,dlat,angmin=1,angmax=179
	oplot, (dlon+ha)/15., dlat

endfor

if REFERENCE ne 0 then twin, /show, /tl		; Show edit window

IPS_gadgets					; draw gadgets on right

stop = 0  &  mode = 0
while (not stop) do begin			; MAIN EDITING LOOP

	pcursor, xp, yp, /data, /kludge		; if button pushed, where?
	IPS_ongadget, xp, yp, result		; is crosshair on a gadget?

	case result of
	0 : begin					; not on gadget
		i = where ( xp-Xt ge 0 , cnt )
		if cnt gt 0 and cnt lt Nt+1 then i = i[cnt-1] else i = -1
		j = where ( yp-Yd ge 0 , cnt )
		if cnt gt 0 and cnt lt Nd+1 then j = j[cnt-1] else j = -1
						; erase only if mode ne 0
		if mode ne 0 and i ne -1 and j ne -1 then begin
			if mode eq 1 then j = indgen(Nd)	; delete column
			xg[i,j,0] = -1.
			ColorBox, Xt[i], Yd[j], widXt, widYd[j], 0, /black
		endif
	end
	3: mode = 1					; delete column (time) mode
	4: mode = 2					; delete one box mode
	5: begin					; delete entire day
		xg[*,*,0] = -1.
		ColorBox , Xt[0:Nt-1] , Yd[0:Nd-1] , widXt, widYd , 0, /black
    end
    else: stop = 1				; quit or save (result 1,2)
	endcase

	if not stop then begin
		DrawAxis = 0
		if result eq 0 then begin
			if mode eq 1 then DrawAxis = 1	else	$
			if mode eq 2 then if i eq Nt-1 or j eq Nd-1 then DrawAxis = 1
		endif else if result eq 5 then DrawAxis = 1

		if DrawAxis then	$		; Redraw axes
			plot, Xt, Yd,	$
				/nodata, /noerase,		$
				xrange= !x.crange,		$
				xticks= 6,			$
				xstyle= 1,  ystyle=1,		$
				yticks= 14, ytickv=dec,		$
				position= POS
	endif

endwhile

twin, /hide, /all			; Hide all windows

if result eq 2 then begin	; save data
	file = filepath(root=ipsEdir[0],'e'+strcompress(mjd[0],/remove_all))
	message, 'Saving data in file '+file, /info
	openw, /get_lun, IU, file, /fortran
	printf, IU, mjd[0]
	savefmt = '('+strcompress(Nd,/remove_all)+'f8.3)'
	for i=0,Nt-1 do printf, IU, format=savefmt,xg[i,*,0]
	free_lun, IU
endif

mjd = mjd+1  &  goto, START

end