;+
; NAME:
;	vu_losmap
; PURPOSE:
;	Plots lines of sight stored in v_los.txt or g_los.txt output
;	from tomography programs.
; CATEGORY:
;	sat/idl/util/vupack
; CALLING SEQUENCE:
	PRO vu_losmap, ff,	$
		carrington_nr	= carrington_nr	, $
		xrange			= xrange		, $
		timenr			= timenr		, $
		losnrs			= losnrs		, $
		all				= all			, $
		_extra			= _extra
; INPUTS:
;	ff				scalar; type: string
;						fully-qualified name of file with los info
;						(typically v_los.txt or g_los.txt)
; OPTIONAL INPUT PARAMETERS:
;	/carrington_nr	if set, plot Carrington rotation number on x-axis
;	xrange=xrange	array[2]; type: float
;						controls range of x-axis (this is passed
;						directly to the IDL plot function).
;						Note that to get Carrington variable increasing
;						left to right xrange[0] must be larger than xrange[1]
;	timenr=timenr	array; type: integer
;						list of times to be plotted
;						(selected from 3rd column of input file)
;	losnrs=losnrs	array; type: integer
;						list of lines of sight to be plotted
;						(selected from 1st column in input file)
; OUTPUTS:
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; CALLS:
;	InitVar
; PROCEDURE:
;	The format of the los text file is:
;    los nr    seg nr   time nr         cr var            lat
;         2        40         4   2055.9404297      9.0538921
;
;	If /carrington_nr NOT set then the heliographic longitude
;	is used for the x-axis. All Carrington variables will be
;	mapped back into the range [0,360], i.e. a single rotation
;	is plotted with longitude running from [0,360].
;
;	If /carrington_nr is SET then the Carrington variable (4th
;	column in input file) is used as x-axis. Since the range
;	of Carrington variables will generally not be an integer
;	number of Carrington rotations, the keyword xrange is
;	usually needed to get the x-axis to look good.
; MODIFICATION HISTORY:
;	DEC-2008, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;-

InitVar, ff, filepath(root=getenv('TUB'),'v_los.txt')

IF IsType(all,/undefine) THEN	$
	IF NOT flt_read(ff,all) THEN RETURN

InitVar, carrington_nr, /key

times = round(reform(all[2,*]))
CASE IsType(timenr,/defined) OF
0: BEGIN
	i = lindgen(n_elements(times))	; Plot all times
	message, /info, 'plotting all times'
END
1: BEGIN
	CASE n_elements(timenr) EQ 1 OF
	0: i = where_common(times,timenr)
	1: i = where(times EQ timenr[0])
	ENDCASE
	message, /info, 'plotting time(s) '+strjoin(strcompress(timenr,/rem),',')
END
ENDCASE

IF i[0] EQ -1 THEN RETURN

los = round(reform(all[0,i]))		; LOS identifier
seg = round(reform(all[1,i]))		; LOS segment numbers
lng = reform(all[3,i])			; Carrington number
lat = reform(all[4,i])			; Heliographic latitude (deg)

CASE carrington_nr OF
0: BEGIN				; Heliographic longitude on x-axis
	lng = Carrington(lng,/get_longitude,/degrees)
	InitVar, xrange, [0,360]

	plot, lng, lat, /nodata, xrange=xrange, yrange=[-90,90], xstyle=1, ystyle=1, $
		xtitle	= 'Longitude (deg)' , $
		ytitle	= 'Latitude (deg)'	, $
		charsize= 2					, $
		xticks	= 6
END
1: BEGIN				; Carrington var on x-axis
 	InitVar, xrange, [ceil(max(lng)),floor(min(lng))]

	xticks = round(abs(xrange[0]-xrange[1])/0.5)
 	plot, lng, lat, /nodata, xrange=xrange, yrange=[-90,90], xstyle=1, ystyle=1, $
		xtitle	= 'Carrington Number' , $
		ytitle	= 'Latitude (deg)'	, $
		charsize= 2					, $
		xticks	= xticks,	$
		yticks  = 6	, $
		xmargin = [10,5]
END
ENDCASE

nsegments = 80

i = uniq(los,sort(los))					; List of unique los numbers (array index into los)
message, /info, 'nr of lines of sight:'+strcompress(n_elements(i))

nseg_total = 0

FOR u=0,n_elements(i)-1 DO BEGIN

	losu = los[i[u]]

	IF IsType(losnrs,/defined) THEN BEGIN
		IF (where(losu EQ losnrs))[0] EQ -1 THEN	$
			continue
		message, /info, 'plotting los '+strcompress(losu)
	ENDIF

	useg = where(los EQ losu, nseg)		; Pick up all los segments for losu
	nseg_total += nseg

	IF nseg NE 0 THEN BEGIN				; nseg LOS segments found

		ulng = lng[useg]
		ulat = lat[useg]
		useg = seg[useg]

		IF max(useg) GT nsegments THEN message, 'too many los segments'
		useg -= 1

		qlng = replicate(!values.f_nan, nsegments)
		qlat = replicate(!values.f_nan, nsegments)
		qseg = replicate(0, nsegments)

		qlng[useg] = ulng
		qlat[useg] = ulat
		qseg[useg] = useg

		;FOR s=0,n_elements(qlng)-1 DO print, qlng[s],qlat[s],qseg[s]

		qlng = qlng[where(finite(qlng))]
		qlat = qlat[where(finite(qlat))]

		nlng = n_elements(qlng)

		CASE carrington_nr OF

		0: BEGIN

			WHILE nlng GT 1 DO BEGIN

				j = where(abs((qlng-shift(qlng,-1))[0:nlng-2]) GT 180)
				IF j[0] EQ -1 THEN BEGIN
					PlotCurve, qlng, qlat, /oplot, /silent, _extra=_extra
					break
				ENDIF

				PlotCurve, qlng[0:j[0]], qlat[0:j[0]], /oplot, /silent, _extra=_extra

				qlng = qlng[j[0]+1:*]
				qlat = qlat[j[0]+1:*]

				nlng = n_elements(qlng)

			ENDWHILE

		END

		1: PlotCurve, qlng, qlat, /oplot, /silent, _extra=_extra

		ENDCASE

	ENDIF

ENDFOR

message, /info, 'nr of segments:'+strcompress(nseg_total)

RETURN  &  END