;+
;NAME:
; 	rapp_get_spectrogram.pro
;PROJECT:
; 	ETHZ Radio Astronomy
;CATEGORY:
; 	Utility
;PURPOSE:
; This routine, given a time interval (anytim format) and a frequency interval (optional),
; returns the associated Phoenix data array, or a string if a problem occured (such as
; no data available). The intensity data is returned as SFU values, unless
; the /POL keyword is used, in which case the values are in % polarization 
; (negative=left).
;
;
;CALLING SEQUENCE:
; rapp_get_spectrogram,time_interval,freq_interval,ELIM=ELIM,BACK=BACK,POL=POL, $
;			xaxis=xaxis,yaxis=yaxis,avgback=avgback
;
;INPUT:
; time_interval : a 2-element array,in a format accepted by ANYTIM
;
;OPTIONAL INPUT:
; freq_interval : a 2-element int array (MHz). If not specified, the whole 
;				frequency range is returned.
;
;KEYWORD INPUT:
;	/NODECOMPRESS: if set, will NOT linearize the data. Normally, no point in using that keyword simultaneously with /LOG...
;			Also, when this keyword is used, the data type ("bunit") is not checked.
;	/LOG  : if set, will output 45*alog( (data > (-9.)) +10.) . Normally, no point in using that keyword simultaneously with /NODECOMPRESS...
;	/ELIM: if set, will automatically remove bad channels
;	/POL : if set, will look for polarization data instead of intensity/L1/L2	
;	/DESPIKE : if set, will use tracedespike.pro (from the TRACE part of the SSW) on the spectrogram. CAUTION: radio 'spikes' would surely be removed...
;	/FAST: if set, will not check each file for time interval (and assumes they are all 15 min. long -- works 99% of the time...)
;	BACK: if set, will automatically remove background 
;			if set to 1 (/BACK): will take the least 3% of the data as background
;			if set to a single number (ex: 0.05): is the least fraction of the data to be considered as background (/BACK or BACK=0.03 does the same)
;			if set to a 2-element array, assumes those are ANYTIM time intervals to be taken as background
;	/ANYTIME: if set, xaxis will return times in anytim format, instead of seconds since midnight...
;	REMOVED!!!: /CALLISTO: if set, will look for data in /ftp/pub/rag/callisto/observations/
;	valid_file_identifiers:
;	archivedir: ex: /ftp/pub/rag/callisto/observations/
;
;OUTPUT:
; a Phoenix-2 spectrogram image, in LINEAR SFU scale (for intensities) or % circular polarization (negative = left)
;
;KEYWORD OUTPUT:
;	xaxis : returns the xaxis, as per RagFitsRead.pro
;	yaxis : returns the yaxis, as per RagFitsRead.pro
;	avgback: array of averaged background values for each frequencies. Useful only if the BACK keyword has also been used!
;
;CALLS
;	has_overlap.pro
;	int2str.pro
;	eventually (keyword /DESPIKE): phoenix_despike.pro
;
;RESTRICTIONS:
;	For now, only works inside ETH (because of file locations)
;
;EXAMPLES:
;	time_intv='2001/09/01 '+['11:50:00','11:55:00']
;	data=rapp_get_spectrogram(time_intv,xaxis=xaxis,yaxis=yaxis,/ELIM,/BACK)
;	data=rapp_get_spectrogram(time_intv,xaxis=xaxis,yaxis=yaxis,/ELIM,BACK=0.05)
;	data=rapp_get_spectrogram(time_intv,xaxis=xaxis,yaxis=yaxis,/ELIM,BACK='2001/09/01 '+['11:45:00','11:46:00'])
;
;	time_intv='2001/09/20 '+['11:35:00','11:36:00']
;	data=rapp_get_spectrogram(time_intv,xaxis=xaxis,yaxis=yaxis,/POL)
;
;
;HISTORY:
;	Pascal Saint-Hilaire, 2001/11/16 created
;		shilaire@astro.phys.ethz.ch
;	PSH 2002/01/21: corrected a bug where a break error occured if no data files were found.
;			Now, simply returns -1 if a problem occured
;
;	PSH 2002/01/28: returns a STRING if something went wrong, instead of -1. This string 
;					usually contains an error message.
;
;	PSH 2002/01/28: major changes, renamed from rapp_get_array to rapp_get_spectrogram
;		and included most of the called subroutines (some of them heavily inspired by 
;		previous work from A. Csillaghy and/or P. Messmer)
;	
;	PSH 2002/05/03: Phoenix file duration is no longer always assumed to be 15 min...
;	PSH 2002/05/04: Added /DESPIKE keyword
;			Added /LOG keyword
;	PSH 2002/05/29: modified background subtraction (accepts a number, or an anytim time array)
;	PSH 2002/07/23: added an environmental variable: RAG_ARCHIVE
;	PSH 2002/07/24: replaced .gz by * for valid_file_identifiers (as some of the files aren't gzipped, for some reason...)
;	PSH 2002/08/05: added the /FAST keyword
;	PSH 2002/09/05: added the avgback keyword output
;	PSH 2002/10/09: corrected a bug where a freq_intv bigger than physical limits would crash the routine.
;	PSH 2002/11/18: corrected a bug where a low limit for the freq interval did not function properly.
;	PSH 2003/10/06: added keyword /ANYTIME
;	PSH 2004/04/26: added keyword /CALLISTO
;	PSH 2004/05/13: -further mods, so that the use of the keyword /CALLISTO also works for background-subtraction.
;			-added KEYWORD /NODECOMPRESS
;	PSH 2004/09/17: added keyword 'valid_file_identifiers' and 'archivedir', for greater flexibility.
;			Removed keyword /CALLISTO.
;			Uses phoenix_despike instead of tracedespike
;	PSH 2006/08/10: added recognition of the "quasi % " format.
;	PSH 2007/05/15: added recognition of the "% polarization (negative = LHCP " format
;-

;==============================================================================================
FUNCTION rapp_return_file_interval,ragfile,NOCHECK=NOCHECK

	IF KEYWORD_SET(NOCHECK) THEN BEGIN
		PRINT,'....not checking file intervals: should be a lot faster!'
		tmp=rsi_strsplit(ragfile,'/',/EXTRACT)
		timeintv=tmp[N_ELEMENTS(tmp)-4]+'/'+tmp[N_ELEMENTS(tmp)-3]+'/'+tmp[N_ELEMENTS(tmp)-2]+' '
		filname=tmp[N_ELEMENTS(tmp)-1]	
		timeintv=anytim(timeintv+STRMID(filname,8,2)+':'+STRMID(filname,10,2)+':'+STRMID(filname,12,2))+[0.,900.]
	ENDIF ELSE BEGIN	
		ragfitsread,ragfile,dateobs=dateobs,timeobs=timeobs,dateend=dateend,timeend=timeend,/SILENT
		timeintv=[dateobs+' '+timeobs,dateend+' '+timeend]
	ENDELSE

RETURN,anytim(timeintv,/ECS)
END

;===================================================================================================================================================================================================================
FUNCTION rapp_get_file_list,time_intv, valid_file_identifiers=valid_file_identifiers,NOCHECK=NOCHECK, archivedir=archivedir
; this function should work for any rag fits after 1978 (excl.)
if not exist(valid_file_identifiers) then valid_file_identifiers='[0-9]'+['*l1.fit*','*l2.fit*','*i.fit*','*p.fit*']
	; OK, I did not include the 'd' and 's' ones... as they were used only the first two days or so...

time_interval=anytim(time_intv)

;some constants
dir_time_interval=3600.*24.	;	the lowest possible directory levels are 1 day long
tmp=GETENV('RAG_ARCHIVE')
IF tmp NE '' THEN rag_root=get_logenv('RAG_ARCHIVE') ELSE rag_root='/ftp/pub/rag/phoenix-2/observations/'	;/global/hercules/data3/rag/observations/'
IF KEYWORD_SET(archivedir) THEN rag_root=archivedir	;'/ftp/pub/rag/callisto/observations/'

;initialize output
ragfile=-1

;get the list of absolute directory paths to be checked (chronological)
dirlist=-1
currentdatetime=time_interval[0]
enddatetime=time_interval(1)
while (currentdatetime le enddatetime) do begin
	bla=anytim(currentdatetime,/EX)
	bla=rag_root+int2str(bla[6],4)+'/'+int2str(bla[5],2)+'/'+int2str(bla[4],2)+'/'
	if datatype(dirlist) eq 'INT' then dirlist=bla else dirlist=[dirlist,bla]
	currentdatetime=currentdatetime+dir_time_interval
endwhile

if datatype(dirlist) eq 'INT' then return,-1

;now, for EACH directory in dirlist, get filelist
filelist=-1
for i=0,n_elements(dirlist)-1 do begin
	bla=findfile(dirlist[i]+'*')

	; now, for each filename in bla, check if ends with correct possible identifiers,
	; and append to filelist

	for j=0,n_elements(valid_file_identifiers)-1 do begin
		temp=WHERE(STRMATCH(bla,valid_file_identifiers[j],/FOLD_CASE) eq 1)
		if temp[0] ne -1 then begin
			temp=bla(temp)
			if datatype(filelist) eq 'INT' then filelist=temp else filelist=[filelist,temp]
		endif
	endfor
endfor

; now, check which ones are in the correct time interval, and return them...
if datatype(filelist) eq 'INT' then RETURN,-1

for i=0,n_elements(filelist)-1 do begin
	if has_overlap(time_interval,anytim(rapp_return_file_interval(filelist[i],NOCHECK=NOCHECK))) then begin
		if datatype(ragfile) eq 'INT' then ragfile=filelist[i] ELSE ragfile=[ragfile,filelist[i]]
	endif
endfor

return,ragfile
END
;===================================================================================================================================================================================================================
PRO rapp_elimwrong, image, xaxis, yaxis 
 print, "Eliminating wrong channels"
 image_old = image
 xaxis_old = xaxis
 yaxis_old = yaxis
 ElimWrongChannels, image, xaxis, yaxis
 if n_elements(image(0,*)) le 1 then begin 
    print, "Restoring original image"
    image = image_old
    xaxis = xaxis_old
    yaxis = yaxis_old
 endif
END
;===================================================================================================================================================================================================================
PRO rapp_merge, im1, im2, xAxis1, xAxis2, yAxis1, yAxis2, $
	im, xAxis, yAxis, ADJUST = adjust

  nx1 = N_Elements( xAxis1 ) & nx2 = N_Elements( xAxis2 )
  ny1 = N_Elements( yAxis1 ) & ny2 = N_Elements( yAxis2 )

  IF nx1 GT nx2 THEN BEGIN
    xAxis = xAxis1
    nxBoth = nx1
    IF Keyword_Set( ADJUST ) THEN $
      im2 = ConGrid1D( im2, nxBoth, ny2 ) 
  ENDIF ELSE  BEGIN
    xAxis = xAxis2
    nxBoth = nx2
    IF Keyword_Set( ADJUST)  THEN $
      im1 = ConGrid1D( im1, nxBoth, ny1 )
  ENDELSE

  yAxis = [ yAxis1, yAxis2 ]
  nyBoth = ny1 + ny2

  im = FltArr( nxBoth, nyBoth )
  im(0, 0) = im1
  im(0, ny1) = im2

END ; Merge images  
;===================================================================================================================================================================================================================
FUNCTION rapp_get_spectrogram,time_intv, freq_intv, ELIM=ELIM, BACK=BACK, POL=POL, DESPIKE=DESPIKE, FAST=FAST, $	;INPUT KEYWORDS
						ANYTIME=ANYTIME, NODECOMPRESS=NODECOMPRESS,LOG=LOG, valid_file_identifiers=valid_file_identifiers, archivedir=archivedir, $
						xaxis=xaxis, yaxis=yaxis,avgback=backavg				;OUTPUT KEYWORDS
	
	IF NOT KEYWORD_SET(valid_file_identifiers) THEN BEGIN
		IF KEYWORD_SET(POL) THEN valid_file_identifiers='*p.fit*' ELSE valid_file_identifiers=['*l1.fit*','*l2.fit*','*i.fit*']
	ENDIF
	ragfiles=rapp_get_file_list(time_intv,valid_file_identifiers=valid_file_identifiers,NOCHECK=FAST, archivedir=archivedir) 

	IF datatype(ragfiles) EQ 'INT' THEN RETURN,'No fits data'
	nb_files=N_ELEMENTS(ragfiles)
	FOR j=0,nb_files-1 DO BEGIN
		ragfitsread,ragfiles[j],img,xax,yax,/SILENT,bunit=bunit
			IF NOT KEYWORD_SET(NODECOMPRESS) THEN BEGIN					
				CASE STRLOWCASE(bunit) OF
					'solar flux units (sfu)': PRINT,bunit
					'45*ln10(solar flux units (sfu) + 10)': BEGIN
						img=10^(img/45.) - 10
						PRINT,'Decompressing '+bunit
					END
					'45*alog(solar flux units (sfu)+10)': BEGIN
						img=exp((img/45.)) - 10
						PRINT,'Decompressing '+bunit
					END						
					'45*log(sfu+10)': BEGIN
						img=10^((img/45.)) - 10
						PRINT,'Decompressing '+bunit
					END						
					'% polarization (negative = left)': PRINT,bunit
					'% polarization (negative = LHCP ': PRINT,bunit
					'quasi % ': PRINT,bunit
					ELSE: BEGIN
						PRINT,'......................The encoding format: "'+bunit+'" is NOT recognized!!!'
						RETURN,'Encoding not recognized'
					END
				ENDCASE
			ENDIF										
				
		IF j GE 1 THEN BEGIN
			rapp_merge, Transpose(image), Transpose(img), $
				yaxis, yax, xaxis, xax, $
				newimage, newyaxis, newxaxis,/ADJUST 
			image = Transpose(newimage)
			xaxis=newxaxis
			yaxis=newyaxis
		ENDIF ELSE BEGIN
			image=img
			xaxis=xax
			yaxis=yax
		ENDELSE
	ENDFOR
	
	;at this point, image is always an image in LINEAR SFU SCALE (or % polarization, for polarization files) !!! [Unless keyword /NODECOMPRESS was used]

	; the following takes only the wanted time interval:
	ss_start=WHERE(xaxis GE (anytim(time_intv[0],/time_only)))	
	IF ss_start[0] EQ -1 THEN ss_start=0  ELSE ss_start=ss_start[0]
	ss_end=WHERE(xaxis LT (anytim(time_intv[1],/time_only)))
	IF ss_end[0] EQ -1 THEN ss_end=N_ELEMENTS(xaxis)-1 ELSE ss_end=ss_end[N_ELEMENTS(ss_end)-1]
	image=image[ss_start:ss_end,*]
	xaxis=xaxis[ss_start:ss_end]		
	IF KEYWORD_SET(ANYTIME) THEN xaxis=anytim(time_intv[0],/date_only)+anytim(xaxis)
	
;DATA PROCESSING:
;==================================================================================================================================================================
	; if a backgnd time interval was entered, should be looked into BEFORE removing bad channels, in order to be sure one has the same freq. channels...
	IF N_ELEMENTS(BACK) EQ 2 THEN BEGIN
		backdata=rapp_get_spectrogram(BACK, freq_intv, POL=POL, valid_file_identifiers=valid_file_identifiers, archivedir=archivedir)
		IF datatype(backdata) NE 'STR' THEN BEGIN
			backavg = avg(backdata,0)
			FOR i=0,N_ELEMENTS(image[0,*])-1 DO image[*,i]=image[*,i]-backavg[i]
		ENDIF ELSE BEGIN
			PRINT,'..........No background data found!!!! Using constbacksub instead...!!!'
			image = constbacksub_v21(image, 0, n_elements(image[*,0])-1,'X', automatic=0.03, background=backavg )
		ENDELSE
	ENDIF

	IF KEYWORD_SET(ELIM) THEN rapp_elimwrong, image, xaxis, yaxis

	IF N_ELEMENTS(BACK) EQ 1 THEN BEGIN
		IF BACK EQ 1 THEN BACK=0.03
		image = constbacksub_v21(image, 0, n_elements(image[*,0])-1,'X', automatic=BACK, background=backavg )
	ENDIF
		
	IF KEYWORD_SET(DESPIKE) THEN image=phoenix_despike(image)
;==================================================================================================================================================================

	IF EXIST(freq_intv) THEN BEGIN	
		ss1=WHERE(yaxis GE freq_intv[0])
		IF ss1[0] EQ -1 THEN RETURN,'Specified frequency interval outside what is available !'
		ss2=WHERE(yaxis LE freq_intv[1])
		IF ss2[0] EQ -1 THEN RETURN,'Specified frequency interval outside what is available !'
		
		tmp=LINDGEN(N_ELEMENTS(yaxis))
		tmp[ss1]=1
		tmp[ss2]=tmp[ss2]+1	
		ss=WHERE(tmp EQ 2)
		IF ss[0] EQ -1 THEN RETURN,'No data in specified frequency range.'

		image=image[*,ss]
		yaxis=yaxis[ss]
		IF EXIST(backavg) THEN backavg=backavg[ss]
	ENDIF
	
	IF KEYWORD_SET(LOG) THEN image=45*ALOG10( (image > (-9.)) +10.)
	
	RETURN,image
END