;+
; NAME:
;	vu_synopticmap
; PURPOSE:
;	Plot a synoptic map extracted from a heliospheric tomography file
; CATEGORY:
;	Tomography
; CALLING SEQUENCE:
	PRO vu_synopticmap, hdr, ff, $
		ut0 		= ut0		, $
		ihdr		= ihdr		, $
		radius		= radius	, $
		type		= type		, $
		ylshow		= ylshow	, $
		degrees 	= degrees	, $
		earth		= earth 	, $
		forecast	= forecast	, $
		printer		= printer	, $
		timeaxis	= timeaxis	, $
		eastlimb	= eastlimb	, $
		westlimb	= westlimb	, $
		upto		= upto		, $
		xysize		= xysize	, $
		module		= module	, $
		diskcenter	= diskcenter, $
		fromcenter	= fromcenter, $
		ascii_file	= ascii_file, $
		mhd_string	= mhd_string, $
		mhd_unit	= mhd_unit	, $
		silent		= silent	, $
		_extra		= _extra
; INPUTS:
;	hdr				array[1]; type: structure
;						header information about the tomographic reconstruction
;						as returned from href=vu_select=. If 'hdr' does not exist,
;						vu_select is called.
;	ff				array[n,m,l,2]; type: float
;						3D volume data; if not present then the files
;						specified in the 'hdr' argument are read.
; OPTIONAL INPUT PARAMETERS:
;	ihdr=ihdr		scalar; type: integer: default: 0
;						If 'hdr' is an array (and 'ff' the corresponding data
;						arrays) then a synoptic map for data at index 'ihdr' is
;						displayed.
;
;	ut0=ut0			scalar or array[1]; type: time structure, float or integer
;						Time (converted to an Earth-based Carrington variable)
;						on which the map is centered. If not specified then the
;						time in 'hdr' is used.
;						If a sequence of files from a time-dependent tomography
;						run is used, and ihdr is not specified, then the hdr
;						closest to ut0 is used.
;
;	type=type		array[1]; type: structure
;						output structure from href=vu_type_insitu=
;
;	/timeaxis       if set, adds a time axis at the top of the map, specifying
;						the time at which the corresponding heliographic
;						longitude crossed the central meridian on the solar disk
;						(as seen from Earth).
;
;	radius=radius	scalar; type: float: default: 1 AU
;						heliocentric distance for synoptic map (AU)
;						The radius is rounded to the nearest grid distance in
;						the 3D heliospheric matrix.
;
;	/earth			sets keyword 'radius' to the heliocentric distance of Earth
;						at the time corresponding to the center of the synoptic
;						map plotted on the screen (see PROCEDURE). Also plots
;						the Earth symbol at the heliographic location of Earth
;						at that time. If /timeaxis is set, also forces plotting
;						of a time series at Earth along the bottom of the map.
;
;	/forecast		if /timeaxis is set, plots vertical dashed line in center of
;						map; and the word 'FORECAST' inside left part of map.
;
;	ylshow=ylshow	array[2]; type: float; default: [-60,60] degrees
;						range of latitudes shown in the synoptic map in units
;						determined by keyword /degrees
;	/degrees		if set, then ylshow is assumed to be in degrees (default: radians)
;	/printer		if set, print to print device (or EPS) instead of to z-buffer
;	ascii_file=ascii_file
;					if set to valid file name the map array will be written to this
;					file (at the resolution of the tomography file).
;					One record for each latitude is written (from bottom to
;					top of the map)
;
;	Keywords used only when a time-ordered sequence is processed:
;
;	/eastlimb		same as fromcenter = -90 (degrees)
;	/westlimb		same as fromcenter = +90 (degrees)
;	/diskcenter 	same as fromcenter =   0 (degrees)
;	fromcenter=fromcenter
;					scalar; type: float; default: none
;						offset in heliographic longitude -180 < fromcenter < 180 degrees)
;
;					If none of the above four keywords is set than a synoptic map
;					at a single time is plotted. The time is specified through keywords
;					ut0 or ihdr.
;
;					If a time-ordered sequence is processed and one of these keywords is
;					set then a synoptic map is build from narrow longitudinal strips
;					with each strip extracted at a different time.
;
;					E.g. if /diskcenter is SET then strips perpendicular to the solar
;					equator at the sub-Earth point (at diskcenter) are extracted from
;					each of sequence members.
;					If /eastlimb or /westlimb are set the strips are extracted above
;					the east or west limb, respectively.

;	Additional keywords can be passed to vu_get_page if /printer is not set
; OUTPUTS:
;	(plot to screen)
; OPTIONAL OUTPUT PARAMETERS:
;	hdr				array[1]; type: structure
;	ff				array[n,m,l,2]; type: float
;						Returns headers and arrays read by vu_select and vu_read
; INCLUDE:
	@compile_opt.pro					; On error, return to caller
	@vu_fnc_include.pro 				;
; CALLS:
;	InitVar, IsTime, Carrington, big_eph 
;	ToRadians, PlotSynopticMap, TimeGet, TimeUnit, GetColors, vu_fnc
;	vu_TimeSeries, set_page, vu_get_page, vu_getdata, gridgen, destroyvar
;	vu_type_insitu, vu_mean, vu_gettime, vu_atlocation, vu_get, IsType, SuperArray
; PROCEDURE:
; >	A t3d file contains a 3D matrix defined on a regular grid in heliographic
;	coordinates. The array processed here covers one Carrington rotation [xc,xc+1].
;	Typically, the array refers to time xc+0.5.
; >	A 2D matrix is extracted at the distance set by keyword 'radius', or implied by
;	keyword /earth.
; >	If argument 'ut' is not specified then the synoptic map plotted covers [xc, xc+1].
; >	If argument 'ut' is specified the plot covers the Carrington variable at time 'ut
;	+/- 0.5 rotation.
; MODIFICATION HISTORY:
;	AUG-1999, Paul Hick (UCSD/CASS)
;	DEC-2000, Paul Hick (UCSD/CASS)
;		Added keyword degree
;		Added keyword ylshow to set latitude range shown in the synoptic map.
;	JAN-2001, Paul Hick (UCSD/CASS)
;		Added keyword /earth.
;		Added capability to plot synoptic map at any radius inside the range
;		covered by the t3d array by linear interpolation on the 3D matrix,
;	MAR-2001, Paul Hick (UCSD/CASS)
;		Rearrange argument list. Arguments 'hdr' and 'type' are now used to
;		set up the display.
;	SEP-2001, Paul Hick (UCSD/CASS)
;		Fixed bug in call to vu_type_insitu
;	JUL-2002, Paul Hick (UCSD/CASS)
;		Added interpolation between neigbouring times of ut0, instead of
;		truncating to nearest grid time.
;	AUG-2002, Paul Hick (UCSD/CASS)
;		Activated option to generate synoptic maps simulating the way some
;		types of synoptic data are made: assembling the map from strips of
;		data (limb scans, slit along central meridian) taken at different times.
;	NOV-2006, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Added ascii_file keyword.
;-

InitVar, silent 	, 0

InitVar, timeaxis	, /key
InitVar, earth		, /key
InitVar, forecast	, /key
InitVar, printer	, /key

InitVar, eastlimb	, /key					; Only one of these three should be set
InitVar, westlimb	, /key
InitVar, diskcenter	, /key

rpm = ToRadians(degrees=degrees)

; 'fromcenter' can only be used when a time-ordered sequence is processed.

IF eastlimb+diskcenter+westlimb NE 0 THEN fromcenter = (-1*eastlimb+0*diskcenter+1*westlimb)*!pi/2/rpm

ihdr_set = IsType(ihdr,/defined)

utt = vu_getdata(hdr, ff, ut0=ut0, ihdr=ihdr, count=count, silent=silent+1, _extra=_extra)
IF count EQ 0 THEN RETURN

utc = Carrington(utt, /get_variable)
utt = Carrington(utt, /get_time    )

; We assume that all 3D matrices have the same basic parameters.
; The first header is used to set up a few constants

hdri = hdr[0]

xcdisplay = vu_get(hdri, /array_range)
xcdisplay = utc+(xcdisplay-mean(xcdisplay))	; Carrington var at central meridian on time utt +/- 0.5 rotation

nTim  = n_elements(hdr)

nLat  = vu_get(hdri, /nlat)
nRad  = vu_get(hdri, /nrad)
R3D	  = vu_get(hdri, /start_distance)
dR3D  = vu_get(hdri, /distance_step)
ylmap = vu_get(hdri, /latitude_range, degrees=degrees)

CASE earth OF
0: InitVar, radius, R3D
1: BEGIN

	; If synoptic map at earth is required, calculate location of Earth in
	; heliographic coordinates. This is used to set 'radius'.
	;'rr_earth' may still be needed later on to plot a time
	; axis at the top of the synoptic map.

	rr_earth = big_eph(utt,	$
		body=jpl_body(/earth,/string),	$
		/to_heliographic	, $
		/to_sphere			, $
		/precess			, $
		degrees=degrees		, $
		/silent)
	rr_earth[0] = utc

	radius = rr_earth[2]
END
ENDCASE


sequence = IsType(fromcenter,/defined) and vu_is_sequence(hdr,silent=silent)


CASE sequence OF

0: BEGIN

	say, 'instant '+TimeGet(utt,/ymd,upto=TimeUnit(/minute))+	$
		 ' (CR'+strcompress(utc)+') at'+strcompress(radius)+' AU', threshold=0, silent=silent

	; Determine the radius indices low and high that bracket 'radius'.
	; Also determine the weight 'rad' needed for linear interpolation between the two levels..

	rad  = (radius-R3D)/dR3D
	rad  = (rad > 0) < (nRad-1)

	low  = floor(rad)
	high = (low+1) < (nRad-1)
	rad  -= low


	; Make a synoptic map using a single density/velocity matrix

	CASE ihdr_set OF
	0: hdri = vu_gettime(hdr, ff, ut=utt, ff_ut=ffi); Interpolate at time utt
	1: hdri = vu_gettime(hdr, ff, ut=vu_get(hdr[ihdr], /uttime), ff_ut=ffi, /fixgrid); ihdr set explicitly
	ENDCASE

	; Extract requested data (density, velocity, etc.). Interpolate on radius

	mapi  = (1-rad)*ffi[*,*,low,*]+rad*ffi[*,*,high,*]

	xcmap = vu_get(hdri, /array_range)				; Only need this to center the rotation on utt

END

1: BEGIN

	say, 'scan '+TimeGet(utt,/ymd,upto=TimeUnit(/minute))+	$
		 ' (CR'+strcompress(utc)+') at'+strcompress(radius)+' AU', threshold=0, silent=silent

	; Make a synoptic map using the whole time-ordered sequence.
	; First get heliographic coordinates at all times.

	tt = SuperArray(Carrington(vu_get(hdr, /uttime)), nlat, /trail)

	map = fltarr(3,ntim,nlat)

	map[0,*,*] = tt-fromcenter*rpm/(2*!pi)				; Carrington var at extraction point
	map[1,*,*] = SuperArray(gridgen(nlat,range=ylmap), ntim, /lead) ; Latitude
	map[2,*,*] = Radius									; Heliocentric distance

	mapi = vu_atlocation(vu_get(hdr,/array), vu_get(hdr,/start_distance),	$
		vu_get(hdr,/distance_step), ff, sequence=vu_get(hdr,/uttime),		$
		location=map, times=tt, degrees=degrees, /periodic, /is_carrington)

	; For now we assume that the times hdr.time are equally spaced

	xcmap = [min(tt),max(tt)]-fromcenter*rpm/(2*!pi)

END

ENDCASE

CASE n_elements(ylshow) OF
0: ylshow = !pi/3/rpm*[-1,1]
1: ylshow = abs(ylshow)*[-1,1]
ELSE:
ENDCASE

;==============================
; Pick up 'type' based on input

type = vu_type_insitu('synoptic', xysize=xysize, hdr=hdr, type=type, _extra=_extra)
destroyvar, module

FOR itype=0,n_elements(type)-1 DO BEGIN

	id = (where( type[itype].data ))[0]

	; Extract fields from 'type' structure

	ctable	= type[itype].ctable

	Fmin	= type[itype].Fmin
	Fmax	= type[itype].Fmax
	display	= type[itype].display
	label	= type[itype].label

	format	= type[itype].format
	charsize= type[itype].charsize
	xysize  = type[itype].xysize

	BreakVal= type[itype].breakval


	i = where(finite(BreakVal), nbreak)
	IF nbreak GT 0 THEN BEGIN
		BreakVal = type[itype].breakval[i]
		IF format EQ '' THEN BreakVal = round(BreakVal)
	ENDIF

	; Done with 'type' structure
	;===========================

	map = vu_fnc(id, mapi, dimension=4-sequence, /t3d_array)

	; Optional: write map to file at the resolution of the tomography files

	IF IsType(ascii_file,/string) THEN BEGIN
		say, 'write '+ascii_file, threshold=-1, silent=silent
		openw, /get_lun, iu, ascii_file, width=strlen(string(map[0]))*(size(map,/dim))[0],append=itype NE 0
		printf, iu, map
		free_lun, iu
	ENDIF

	set_page, $
		zbuffer		= 1-printer	, $
		printer		= printer	, $
		xysize		= xysize	, $
		fullsize	= printer	, $
		_extra		= _extra	, $
		old_device	= old_device, $
		ctable		= ctable	, $
		/silent

	; Plot the synoptic map.
	; The center of the window is xysize/2. The map is shifted to the right to move it
	; away from the vertical color bar. The size of the window is controlled by plotsize.

	PlotSynopticMap, xcdisplay, map, 				$
		xcmap		= xcmap,						$
		ylmap		= ylmap,						$
		ylplot		= ylplot,						$
		ylshow		= ylshow, degrees=degrees,		$
		breakval	= BreakVal,						$
		format		= format,						$
		plotcenter	= xysize/2+[0.04*xysize[0],0],	$	; Center plot right of window center
		plotsize	= 0.38*xysize[0],				$
		charsize	= charsize,						$
		timeaxis	= timeaxis,						$
		_extra		= _extra

	; Title above map

	mhd_s = IsType(mhd_string,/defined) ? mhd_string : vu_fnc(id,/symbol,/norm)
	mhd_u = IsType(mhd_unit  ,/defined) ? mhd_unit   : vu_fnc(id,/unit)
	xyouts, 0.05, 0.95, /norm, mhd_s+' ('+mhd_u+')'+' at '+string(format='(F4.2)',radius)+' AU', charsize=charsize, _extra=_extra

	IF timeaxis AND forecast THEN BEGIN
		n = 20
		cc = gridgen(n, range=!y.crange)
		xc = mean(XCdisplay)*[1,1]				; Vertical dashed line in center
		FOR i=0,n-1,2 DO oplot, xc, cc[i:i+1];, _extra=_extra
												; 'FORECAST' label inside map
		xyouts, total(([0,1]+0.75*[1,-1])*!x.crange), total(([0,1]+.1*[1,-1])*!y.crange), 'FORECAST', align=0.5, charsize=charsize;, _extra=_extra
	ENDIF

	; UT time above map

	xyouts, 0.99, 0.95, /norm, align=1, TimeGet(utt, /ymd, upto=upto, /round)+' UT', charsize=charsize, _extra=_extra

	IF earth THEN BEGIN

		IF timeaxis THEN BEGIN					; Plot earth time series along bottom

			; Get the time series at Earth.

			destroyvar, T1						; Safety belt: make sure T1 does not exist
			F1 = vu_TimeSeries(xcmap, R3D, dR3D, vu_fnc(id,ffi,/t3d_array), range=xcdisplay, T=T1, radius=radius)

			IF NOT finite(Fmin) THEN Fmin = min(F1,/nan)
			IF NOT finite(Fmax) THEN Fmax = max(F1,/nan)
			say, 'minimum='+strcompress(Fmin)+'  maximum='+strcompress(Fmax), threshold=-1, silent=silent

			; Plot the timeseries along the bottom of the map

			T1 = Carrington(T1)

			Yb = ylplot[0]
			Yt = 0.25*(ylplot[1]-ylplot[0])

			oplot, T1, Yb+(F1-Fmin)/(Fmax-Fmin)*Yt;, _extra=_extra

			dxc = 0.015*(xcdisplay[0]-xcdisplay[1])
			FOR i=0,3 DO plots, xcdisplay[0]+[0,dxc], Yb+Yt*(i/3.)+[0,0]

			xyouts, xcdisplay[0]+1.5*dxc, Yb   , strcompress(round(Fmin), /rem), charsize=charsize;, _extra=_extra
			xyouts, xcdisplay[0]+1.5*dxc, Yb+Yt, strcompress(round(Fmax), /rem), charsize=charsize;, _extra=_extra

			; Get map value at Earth location by interpolation on the time series

			F1 = interpol(F1,T1,rr_earth[0])

			CASE finite(F1) OF
			0: F1 = !p.color
			1: BEGIN
				Yb += (F1-Fmin)/(Fmax-Fmin)*Yt

				plots , xcdisplay[0]-[0,dxc], Yb+[0,0];, _extra=_extra
				xyouts, xcdisplay[0]-1.5*dxc, Yb, strcompress(round(F1), /rem), align=1.0, charsize=charsize;, _extra=_extra

				; Use the value at Earth to find a contrasting color with the map.

				F1 = GetColors(F1, BreakVal, /open, /badbackground, _extra=_extra)
				F1 = (F1 LT !d.n_colors/2)*(!d.n_colors-1)
			END
			ENDCASE

		ENDIF

		; Plot Earth symbol at location Earth

		cc = gridgen(33, range=[0,2*!pi])
		usersym, [cos(cc),-1,0,0,0], [sin(cc),0,0,1,-1], _extra=_extra	; Circle with cross inside
		oplot, [rr_earth[0]], [rr_earth[1]], color=F1, psym=8, thick=0.5, symsize=1.3*charsize;, _extra=_extra

	ENDIF

	boost, module, display+label
	vu_get_page, hdri, title='', module[itype], ut=utt, device=old_device, upto=upto, silent=silent, new=itype ne 0, _extra=_extra

ENDFOR

RETURN  &  END