FUNCTION dustem_sed_extractor,maps, $
                              index_roi, $
                              filters, $

                              comments=comments, $
                              maps_order=maps_order, $
                              method=method, $
                              index_off=index_off, $
                              reference_filter=reference_filter, $
                              mask_I_undef=mask_I_undef, $
                              total_intensity_only=total_intensity_only, $

                              set_undefined=set_undefined, $

                              _extra=_extra

;+
; NAME:

;    dustem_sed_extractor

;

; PURPOSE:
;    Extract SEDs from a set of maps 

;

; CATEGORY:

;    DustEMWrap, Distributed, High-Level, Analysis Helper
;

; CALLING SEQUENCE:
;	sed=sed_extractor(maps,index,filters[,comments=][,maps_order=][,method=][,mask_I_undef=][,/total_intensity_only])

;

; INPUTS:

;    maps      : a set of Nmaps maps of same dimension [Nx,Ny,Nmaps,Ndata_sets]. Caution, maps are changed if set_undefined is set !

;    index_roi : spatial index in Nx,Ny where the SED is to be extracted from (ROI=region of interest)

;    filters   : DustEMWrap filters corresponding to each dataset [Ndata_sets]
;

; OPTIONAL INPUT PARAMETERS:
;    comments   = a vector of strings to be attached to the SED. Also an output.
;    maps_order = a vector of 9 integers giving the 0-based index for I,Q,U,II,QQ,UU,IQ,IU,QU in that order
;                 (default [0,1,2,3,4,5,6,7,8]). For total intensity only, use something like [0,-1,-1,1,-1,-1,-1,-1,-1].
;    method = extraction method (default = 'MEAN')/ Possible values : 'MEAN','MEAN_ONOFF','QU_CORREL'
;    index_off = same as index for the OFF region (needed only for method='MEAN_ONOFF')
;    reference_filter = reference filter to be correlated with each map (needed only for method='QU_CORREL')

;

; OUTPUTS:
;    SED  : a dustemwrap SED structure containing te extracted SED

;

; OPTIONAL OUTPUT PARAMETERS:
;    mask_I_undef = a map of saturated pixels found in Stokes I

;    comments   = a vector of strings to be written as comments into the SED.
;

; ACCEPTED KEY-WORDS:
;    total_intensity_only   = if set, data is assumed to contain only total intensity (no polarization)
;    comments               = comments to be added to the sed. sed_extractor adds comments to this

;

; COMMON BLOCKS:
;    None

;

; SIDE EFFECTS:
;    None

;

; RESTRICTIONS:
;    None

;

; PROCEDURE:
;    None

;

; EXAMPLES

;

; MODIFICATION HISTORY:

;    Written by JPB 2022
;    Evolution details on the DustEMWrap gitlab.
;    See http://dustemwrap.irap.omp.eu/ for FAQ and help.  

;-

IF keyword_set(help) THEN BEGIN
  doc_library,'dustem_sed_extractor'
  sed=0.
  goto,the_end
ENDIF

Nfilt=n_elements(filters)
sed=dustem_initialize_sed(Nfilt,comments=comments)
sed.filter=filters
sed.instru=dustem_filter2instru(sed.filter)
sed.wave=dustem_filter2wav(filters)
comments=[comments,'\Computed using sed_extractor on '+systime(0)]
order=sort(sed.wave)  ;wavelength order

count_roi=n_elements(index_roi)

use_method='MEAN'  ;This is the default extraction method
IF keyword_set(method) THEN BEGIN
	use_method=method
ENDIF
comments=[comments,'SED extracted using method '+use_method]
comments=[comments,'ROI has '+strtrim(count_roi,2)+' map pixels']
N1=(size(maps))[1]
N2=(size(maps))[2]
Nplanes=(size(maps))[3]

IF keyword_set(total_intensity_only) THEN BEGIN
	indI=0L & indII=1L
	IF keyword_set(maps_order) THEN BEGIN
		indI=maps_order[0]
		indII=maps_order[1]
	ENDIF
ENDIF ELSE BEGIN
	indI=0L & indQ=1L & indU=2L & indII=3L & indQQ=4L & indUU=5L &indIQ=6L & indIU=7L &indQU=8L
	IF keyword_set(maps_order) THEN BEGIN
		indI=maps_order[0]
		indQ=maps_order[1]
		indU=maps_order[2]
		indII=maps_order[3]
		indQQ=maps_order[4]
		indUU=maps_order[5]
		indIQ=maps_order[6]
		indIU=maps_order[7]
		indQU=maps_order[8]
	ENDIF
ENDELSE

;used_maps=maps

IF keyword_set(set_undefined) THEN BEGIN

	used_maps=maps

	;===== do a map of saturated stokesI pixels
	mask_I_undef=maps[*,*,0,0]
	mask_I_undef[*]=0

	FOR j=0L,Nfilt-1 DO BEGIN

		this_map=maps[*,*,indI,j]
		ind=where(this_map EQ la_undef(),count)
		IF count NE 0 AND count NE N1*N2 THEN BEGIN
			message,'Found '+strtrim(count,2)+' undefined values for map of '+filters[j],/continue
			;stop
			mask_I_undef[ind]=1
		ENDIF

	ENDFOR

	;set saturated stokesI pixels to la_undef in 
	indbad=where(mask_I_undef EQ 1,count)
	IF count NE 0 THEN BEGIN
		message,'Setting Stokes I saturated values to undef in the whole cube',/continue
		FOR j=0L,Nfilt-1 DO BEGIN
			FOR k=0L,Nplanes-1 DO BEGIN
				this_map=used_maps[*,*,k,j]
				this_map[indbad]=la_undef()
		  		used_maps[*,*,k,j]=this_map
		 	ENDFOR
		ENDFOR
	ENDIF

	maps=used_maps
	used_maps=0

ENDIF

CASE use_method OF
	'MEAN':BEGIN

		   FOR j=0L,Nfilt-1 DO BEGIN
                   this_map=maps[*,*,indI,j]

                   sed[j].stokesI=la_mean(this_map[index_roi])

                   this_map=maps[*,*,indII,j]

                   sed[j].sigmaII=la_div(la_mean(this_map[index_roi]),count_roi)
                   IF not keyword_set(total_intensity_only) THEN BEGIN

                      this_map=maps[*,*,indQ,j]

                      sed[j].stokesQ=la_mean(this_map[index_roi])

                      this_map=maps[*,*,indU,j]

                      sed[j].stokesU=la_mean(this_map[index_roi])

                      this_map=maps[*,*,indQQ,j]

                      sed[j].sigmaQQ=la_div(la_mean(this_map[index_roi]),count_roi)

                      this_map=maps[*,*,indUU,j]

                      sed[j].sigmaUU=la_div(la_mean(this_map[index_roi]),count_roi)
                      

                      this_map=maps[*,*,indIQ,j]

                      sed[j].sigmaIQ=la_div(la_mean(this_map[index_roi]),count_roi)

                      this_map=maps[*,*,indIU,j]

                      sed[j].sigmaIU=la_div(la_mean(this_map[index_roi]),count_roi)

                      this_map=maps[*,*,indQU,j]

                      sed[j].sigmaQU=la_div(la_mean(this_map[index_roi]),count_roi)
                   ENDIF

            ENDFOR
         END

        'MEAN_ONOFF': BEGIN     ;This is mean ON - mean MEAN_ONOFF
           IF not keyword_set(index_off) THEN BEGIN
              message,'index_off keyword must be set with method '+use_method,/continue
              stop
           ENDIF
           count_off=n_elements(index_off)
           comments=[comments,'OFF has '+strtrim(count_off,2)+' map pixels']
           FOR j=0L,Nfilt-1 DO BEGIN

              this_map=maps[*,*,indI,j]

              sed[j].stokesI=la_sub(la_mean(this_map[index_roi]),la_mean(this_map[index_off]))

              this_map=maps[*,*,indII,j]

              sig_on=la_div(la_mean(this_map[index_roi]),count_roi)
              sig_off=la_div(la_mean(this_map[index_off]),count_off)
              sed[j].sigmaII=la_add(sig_on,sig_off)
              IF not keyword_set(total_intensity_only) THEN BEGIN

                 this_map=maps[*,*,indQ,j]

                 sed[j].stokesQ=la_sub(la_mean(this_map[index_roi]),la_mean(this_map[index_off]))

                 this_map=maps[*,*,indU,j]

                 sed[j].stokesU=la_sub(la_mean(this_map[index_roi]),la_mean(this_map[index_off]))
                                ;variances

                 this_map=maps[*,*,indQQ,j]

                 sig_on=la_div(la_mean(this_map[index_roi]),count_roi)
                 sig_off=la_div(la_mean(this_map[index_off]),count_off)
                 sed[j].sigmaQQ=la_add(sig_on,sig_off)

                 this_map=maps[*,*,indUU,j]

                 sig_on=la_div(la_mean(this_map[index_roi]),count_roi)
                 sig_off=la_div(la_mean(this_map[index_off]),count_off)			
                 sed[j].sigmaUU=la_add(sig_on,sig_off)
                                ;co-variances

                 this_map=maps[*,*,indIQ,j]

                 sig_on=la_div(la_mean(this_map[index_roi]),count_roi)
                 sig_off=la_div(la_mean(this_map[index_off]),count_off)			
                 sed[j].sigmaIQ=la_add(sig_on,sig_off)

                 this_map=maps[*,*,indIU,j]

                 sig_on=la_div(la_mean(this_map[index_roi]),count_roi)
                 sig_off=la_div(la_mean(this_map[index_off]),count_off)			
                 sed[j].sigmaIU=la_add(sig_on,sig_off)

                 this_map=maps[*,*,indQU,j]

                 sig_on=la_div(la_mean(this_map[index_roi]),count_roi)
                 sig_off=la_div(la_mean(this_map[index_off]),count_off)			
                 sed[j].sigmaQU=la_add(sig_on,sig_off)
              ENDIF
           ENDFOR	
                                ;stop

	END
	'QU_CORREL': BEGIN    ;This is to compute SED of polarization correlated to a reference
        IF keyword_set(total_intensity_only) THEN BEGIN
        	message,'mode QU_CORREL can only be used with polarization data'
        ENDIF
	    IF not keyword_set(reference_filter) THEN BEGIN
	    	message,'reference_filter keyword must be set with method '+use_method,/continue
	    	stop
	    ENDIF
	    indf=where(filters EQ reference_filter,countf)
	    IF countf EQ 0 THEN BEGIN
	    	message,'reference filter '+reference_filter+' not found in provided filter list',/continue
	    	stop
	    ENDIF

        mapI_ref=maps[*,*,indI,indf]
        mapQ_ref=maps[*,*,indQ,indf]
        mapU_ref=maps[*,*,indU,indf]
        mapII_ref=maps[*,*,indII,indf]
        mapQQ_ref=maps[*,*,indQQ,indf]
        mapUU_ref=maps[*,*,indUU,indf]

		Iref=la_mean(mapI_ref[index_roi])
		Qref=la_mean(mapQ_ref[index_roi])
		Uref=la_mean(mapU_ref[index_roi])
		FOR j=0L,Nfilt-1 DO BEGIN
			;==== do the correlation in Stokes I
			x=mapI_ref[index_roi]

			y=(maps[*,*,indI,j])[index_roi]

			sx=la_power(mapII_ref[index_roi],0.5)

			sy=la_power((maps[*,*,indI,j])[index_roi],0.5)

			ind_good=where(x NE la_undef() AND y NE la_undef() AND sx NE la_undef() and sy NE la_undef(),count_good)
			comments=[comments,'I corelation has '+strtrim(count_good,2)+' map pixels']
			IF count_good NE 0 THEN BEGIN
	            x=x[ind_good]
	            y=y[ind_good]
	            sx=sx[ind_good]
	            sy=sy[ind_good]
	            start=[1.,1.]
		    	try=linear_mpfit(x,y,sx,sy,start,status=status $
	                      ,perror=perror,bestnorm=bestnorm,_extra=_extra)
		    	;stop
	    	ENDIF ELSE BEGIN
	    		status=0
	    		try=[la_undef(),la_undef()]
	    	ENDELSE
	    	IF status NE 0 AND status NE -16 THEN BEGIN
				sed[j].stokesI=try[0]*Iref
				sed[j].sigmaII=la_power(la_mul(perror[0],Iref),2.)
			ENDIF
			;Do the correlation between Q maps
			x=[mapQ_ref[index_roi],mapU_ref[index_roi]]

			y=[(maps[*,*,indQ,j])[index_roi],(maps[*,*,indU,j])[index_roi]]

			sx=la_power([mapQQ_ref[index_roi],mapUU_ref[index_roi]],0.5)

			sy=la_power([(maps[*,*,indQ,j])[index_roi],(maps[*,*,indU,j])[index_roi]],0.5)

			ind_good=where(x NE la_undef() AND y NE la_undef() AND sx NE la_undef() and sy NE la_undef(),count_good)
			comments=[comments,'QU corelation has '+strtrim(count_good,2)+' map pixels']
			IF count_good NE 0 THEN BEGIN
	            x=x[ind_good]
	            y=y[ind_good]
	            sx=sx[ind_good]
	            sy=sy[ind_good]
	            start=[1.,1.]
		    	try=linear_mpfit(x,y,sx,sy,start,status=status $
                      ,perror=perror,bestnorm=bestnorm,_extra=_extra)
		    	;stop
	    	ENDIF ELSE BEGIN
	    		status=0
	    		try=[la_undef(),la_undef()]
	    	ENDELSE
	    	IF status NE 0 AND status NE -16 THEN BEGIN
				sed[j].stokesQ=la_mul(try[0],Qref)
				sed[j].stokesU=la_mul(try[0],Uref)
				sed[j].sigmaQQ=la_power(la_mul(perror[0],Qref),2.)
				sed[j].sigmaUU=la_power(la_mul(perror[0],Uref),2.)
	            ;for now assume no co-variances. But we could aclculate this, right ?
				sed[j].sigmaIQ=0.
				sed[j].sigmaIU=0.
				sed[j].sigmaQU=0.
				;IF filters[j] EQ reference_filter THEN stop
			ENDIF
		ENDFOR
		comments=[comments,'Reference filter for QU correlation '+reference_filter]     
	END
    ELSE: BEGIN
    	message,'Method '+use_method+' not implemented',/continue
    	sed=0L
    	goto,the_end
    END
ENDCASE

IF not keyword_set(total_intensity_only) THEN BEGIN
	polar_iqu2ippsi,sed.stokesI,sed.stokesQ,sed.stokesU,p_values,psi_values,/lac,largep=largep
	sed.smallp=p_values
	sed.psi=psi_values
	sed.largep=largep
	polar_variance_iqu2ippsi,sed.stokesI,sed.stokesQ,sed.stokesU,sed.sigmaII,sed.sigmaQQ,sed.sigmaUU,sed.sigmaIQ,sed.sigmaIU,sed.sigmaQU,variances_p,variances_psi,variances_largep,/lac
	sed.sigma_smallp=variances_p
	sed.sigma_largep=variances_largep
	sed.sigma_psi=variances_psi
	ssed=sed[order]
	sed=ssed
ENDIF

;stop

the_end:
RETURN,sed

END