;+
; NAME:
;	nrh_map.pro
;
; PURPOSE:
;	General purpose routine to retrieve maps from the NRH.
;
; CALLING SEQUENCE:
;	
;	maps=nrh_map(time_intv, ...)
;
; INPUT:
;	time_intv : a 2-element array of times in any of the ANYTIM formats
;	
; KEYWORDS (INPUT):
;
;	dir	: directory where files are located
;	MEUDON	: if set, assumes we are working in Meudon, on the server where the data are. 'dir' is not used.
;	archive_path	: assumes files are located in a $archive_path/yyyy/mm/dd/ directory structure
;
;	QUICKLOOK: if set, will look for '2q' files instead of default '2i' files. Combination with /MEUDON keyword not implemented yet.
;
;	STOKESV	: if set, will retrieve Stokes V maps instead of default Stokes I maps. Looks like data before 1999/04/20 does not have Stokes V...
;	SFU	: if set, returns map pixel values in SFU instead of brightness temperature [Kelvin] (default).
;
;	size	: map pixel length. Allowed values are 128 (default), 256, and 512.
;
;	CLEAN	: if set, maps are CLEANed. Possible clean modes:
;			0: no clean
;			1: usual clean (with residuals added)
;			2: clean without residuals
;			3: clean components outside of clean box are returned
;			4: ...
;		Later, routine should be modified to propagate the whole set of CLEAN parameters...
;
;	xrange, yrange: sub maps' x- and y- intervals to be returned.
;
;	ALIAS_MIN_DISTANCE_ESTIMATE: if set, will add a tag in each map structure containing the estimated minimal distance (in arcsececonds) between a source and its alias.
;					This estimation is based on the 100m (the longest) baseline of the E-W array, and should be exact for sources and alias aligned with Earth's E-W axis.
;					2010/06/11: this feature was disabled, as it was not working properly. Now replaced by the much better ALIAS_VECTOR keyword.
;
;	ALIAS_VECTOR: if set, yields the [x,y] displacement (in arcsecs) on a solar map between aliases due to the EW and NS network of antennae.
;
;	STARTENDTIMES : if set, output will be start/end times of daily data file instead of maps.
;	CATCH_ERRORS=CATCH_ERRORS : if set, catches problems and returns -1.
;
;	/PSF: is set, will add PSF information to maps (HWHM of major axis & minor axis, tilt angle from x-axis in radians)
;		(HWHM = 1.177*sigma). The routine ellipse_maille.pro --which might be outdated-- is used.
;
; OUTPUT:
;	If keyword /STARTENDTIMES is set, will return a 2-element array of doubles: these are the ANYTIM times of the beginning and end of data in daily file.
;	Otherwise, returns an array [n,m] of map structures, where n is the number of time intervals, and m the number of frequencies. 
;
;
; EXAMPLES:
;
;	maps=nrh_map('2003/03/27 14:56:56', size=512 ,/MEUDON)
;	plot_map, maps[0,2],/LIMB
;
;	nrh_data_intv=nrh_map2('2007/01/23',/STARTENDTIMES ,/MEUDON)
;	ptim,  nrh_data_intv
;		
;	maps=nrh_map('2002/04/16 '+['13:11:30','13:11:41'],/MEUDON,/ALIAS,CLEAN=3)
;	plot_map, maps[0,2],/LIMB
;
;	maps=nrh_map('2002/02/20 '+['10:50','10:50:10'],/MEUDON,/CLEAN)
;	plot_map, maps[0,4],/LIMB
;
;	maps=nrh_map('2002/02/20 09:00', archive_path='/home/shilaire/data/nrh/',/PSF)
;	i=0 & j=1
;	plot_map, maps[i,j],/LIMB
;	tvellipse, maps[i,j].psf[0], maps[i,j].psf[1], 0.,0., maps[i,j].psf[2]/!DPI*180, 255, /DATA,/MAJOR,/MINOR
;
;
; RESTRICTIONS:
;	-It is assumed that time_intv does not span more than one day!!
;
; HISTORY:
;	Pascal Saint-Hilaire (shilaire@ssl.berkeley.edu), 2008/09/22, adapted from previous software.
;	PSH, 2010/04/19: added psf stuff (needs newer version of ellipse_maille.pro).
;	PSH, 2010/06/11: deactivated keyword ALIAS_MIN_DISTANCE_ESTIMATE, because it was wrong! Replaced by /ALIAS_VECTOR (which necessitates PSH's nrh_alias_vector.pro)
;
;-
FUNCTION nrh_map, time_intv, STOKESV=STOKESV, size=size, QUICKLOOK=QUICKLOOK, xrange=xrange, yrange=yrange, $
	STARTENDTIMES=STARTENDTIMES, CATCH_ERRORS=CATCH_ERRORS, PSF=PSF, $
	dir=dir, MEUDON=MEUDON, archive_path=archive_path, $
	SFU=SFU, ALIAS_VECTOR=ALIAS_VECTOR, $
	CLEAN=CLEAN


	;CATCH+++
	IF KEYWORD_SET(CATCH_ERRORS) THEN BEGIN
		err_status=1
		CATCH,err_status
		IF err_status NE 0 THEN BEGIN
			CLOSE,/ALL
			RETURN,-1
		ENDIF
	ENDIF
	;CATCH---

	time_intv2=anytim(time_intv)
	IF N_ELEMENTS(time_intv2) EQ 1 THEN time_intv2+=[-10,10]
	IF KEYWORD_SET(QUICKLOOK) THEN filetype='q' ELSE filetype='i'
	IF KEYWORD_SET(MEUDON) THEN BEGIN
		IF time_intv2[0] LT anytim('1999/04/20') THEN basedir='/juke/10secbis/' ELSE basedir='/juke/10sec/'
		pattdir=basedir+'INT'+STRMID(time2file(anytim(time_intv[0]),/year2digit,/date_only),0,4)+'?'
		files=FINDFILE(pattdir+'/2'+filetype+time2file(time_intv[0],/date_only, /year2digit)+'.??')
	ENDIF ELSE BEGIN
		IF KEYWORD_SET(archive_path) THEN dir=archive_path+psh_time2dir(time_intv[0]) ELSE default, dir, './'
		files=FINDFILE(dir+'/2'+filetype+time2file(time_intv[0],/date_only, /year2digit)+'.??')
	ENDELSE

	IF files[0] EQ '' THEN BEGIN
		MESSAGE,/INFO,'NRH file for '+anytim(time_intv[0],/ECS)+' not found. Not present locally?...'
		RETURN,-1
	ENDIF

	FOR j=0, N_ELEMENTS(files)-1 DO BEGIN
		;PSH, 2008/09/19: better define HBEG and HEND at the start of every j loop, as read_nrh.pro modifies HBEG and HEND !!!
		HBEG=anytim(time_intv[0],/time_only,/ECS)
		IF N_ELEMENTS(time_intv) GT 1 THEN HEND=anytim(time_intv[1],/time_only,/ECS) ELSE HEND=HBEG

		break_file, files[j], DISK_LOG, DIR, FILNAM, EXT, FVERSION, NODE
		file=FILNAM+EXT
		IF KEYWORD_SET(STARTENDTIMES) THEN BEGIN
			read_nrh, file, index, dir=dir
			curstart=anytim(index.DATE_OBS)
			read_nrh, file, index, dir=dir, HBEG='23:00:00'
			curend=anytim(index.DATE_OBS)
			IF exist(nrh_data_intv) THEN BEGIN
				curstart= nrh_data_intv[0] < curstart
				curend=   nrh_data_intv[1] > curend
			ENDIF			
			nrh_data_intv=[curstart,curend]
		ENDIF ELSE BEGIN
			nfreq=1 & i=0
			WHILE i LT nfreq DO BEGIN
				read_nrh, file, index, data, dir=dir, HBEG=HBEG, HEND=HEND, freq=i, STOKES=STOKESV, size=size, nfrq=nfrq, SFU=SFU, CLEAN=CLEAN
				nfreq=nfrq
				index2map,index,data,newmaps
				t=anytim(newmaps.time)
				ss=WHERE((t GE  anytim(time_intv2[0])) AND (t LE anytim(time_intv2[1])))	;could have done this earlier, and spare the computing power required in index2map for useless maps...
														;also, potential pitfall, as higher freqs can have times ~20 ms later than lower freqs....
				IF ss[0] NE -1 THEN BEGIN
					newmaps=newmaps[ss]
					
					IF KEYWORD_SET(xrange) OR KEYWORD_SET(yrange) THEN BEGIN
						sub_map, newmaps, smaps, xrange=xrange, yrange=yrange
						newmaps=smaps
					ENDIF
					
					newmaps=add_tag(newmaps,index[0].FREQ,'freq')
					newmaps.ID='NRH '+strn(index[0].FREQ, format='(f10.1)')+' MHz'
					
					IF KEYWORD_SET(ALIAS_VECTOR) THEN BEGIN
						newmaps=add_tag(newmaps,[0.,0.],'alias_vector_ew')
						newmaps=add_tag(newmaps,[0.,0.],'alias_vector_ns')
						FOR ii=0, N_ELEMENTS(ss)-1 DO BEGIN   
							nrh_alias_vector, anytim(newmaps[ii].TIME), newmaps[ii].FREQ, v1,v2
							newmaps[ii].alias_vector_ew=v1	;;in arcsecs
							newmaps[ii].alias_vector_ns=v2	;;in arcsecs
						ENDFOR;ii
					ENDIF

					IF KEYWORD_SET(PSF) THEN BEGIN
						IF NOT tag_exist('psf',newmaps,/Q) THEN newmaps=add_tag(newmaps, [0.,0,0],'psf',/NO_PAIR)
						FOR ii=0, N_ELEMENTS(ss)-1 DO newmaps[ii].psf=ellipse_maille(struct2fitshead(index[ss[ii]]),1d3*anytim(index[ss[ii]].DATE_OBS,/TIME_ONLY))
						newmaps.psf[0:1]*=psh_solrad(index[0].DATE_OBS)	;; convert from solar radius to arcseconds
						newmaps.psf[2]*= !DPI/180.				;; convert from degrees to radians.
					ENDIF
					
					IF exist(maps) THEN maps=[[maps],[newmaps]] ELSE maps=newmaps
				ENDIF
				i+=1
			ENDWHILE;i
		ENDELSE
	ENDFOR;j
	IF KEYWORD_SET(STARTENDTIMES) THEN RETURN, nrh_data_intv
	IF exist(maps) THEN RETURN, maps ELSE RETURN,-1
END