;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

PRO aia_psf_correction_jr, data_raw, aia_psf_fft_channel, data_cor, error=error

; correcting the RAW_DATA AIA image for the already prepared PSF influence
; depending on the selected channel and data size

; INPUT: data_raw  - RAW_DATA AIA image
;                    RESTRICTION: raw AIA image data to be deconvolved for the PSF has to
;                                 be larger in both dimensions than 801 points
;        aia_psf_fft_channel - FFT of PSF (complex)
;                    
; OUTPUT: data_cor - AIA DATA corrected for the PSF (real)
;         error    - 0=creation OK, 1=too small data size

; PSFs of SDO AIA derived by the procedure 'aia_psf_fft_creation_jr' based on article and data of
; title  : POINT-SPREAD FUNCTIONS FOR THE EXTREME-ULTRAVIOLET CHANNELS OF SDO/AIA TELESCOPES'
; autors : B. Poduval, C. E. DeForest, J. T. Schmelz, and S. Pathak
; journal: The Astrophysical Journal, 765:144 (15pp), 2013 March 10

; OK
error=0

; too small data size
data_size_tmp=size(data_raw)
if min([data_size_tmp(1), data_size_tmp(2)]) lt 801 then begin
     print
     print,' INFO : too small data size, it has to be >= 801 !'
     error=1 & data_cor=0
     return
endif

; parameters
data_raw_size=size(data_raw)
data_raw_size_max=max([data_raw_size(1),data_raw_size(2)])
data_raw_size_max_half=data_raw_size_max/2
; data array for the result
data_cor=fltarr(data_raw_size_max,data_raw_size_max)

; cases of the non-squared data array
if data_raw_size(1) ne data_raw_size(2) then begin 
   data_raw_tmp=fltarr(data_raw_size_max,data_raw_size_max)
   x1=0 & x2=data_raw_size_max-1 & y1=0 & y2=data_raw_size_max-1
   if data_raw_size(1) lt data_raw_size(2) then begin ; a tall rectangle
      x1=fix((data_raw_size_max-data_raw_size(1))/2.)
      x2=x1+data_raw_size(1)-1
   endif
   if data_raw_size(1) gt data_raw_size(2) then begin ; a long rectangle
      y1=fix((data_raw_size_max-data_raw_size(2))/2.)
      y2=y1+data_raw_size(2)-1
   endif   
   data_raw_tmp(x1:x2,y1:y2)=data_raw
   ; carrousel - temporal data keeping action
   data_raw_stored=data_raw
   data_raw=data_raw_tmp
endif

; FFT of the raw_ data
fft_data=fft(data_raw,-1)
; division of the FFT DATA by FFT PSF = deconcolution of the DATA and the PSF
div=fft_data/aia_psf_fft_channel
; IFFT of the deconcolved data 
res=fft(div,+1)
; real part of the deconvolved data
res_real=real_part(res)

; replacement of the resulting data segmets to a structure of the
; original AIA image
; 1/ BL - bottom left corner
data_cor(0:data_raw_size_max_half-1,0:data_raw_size_max_half-1)=$
   res_real(data_raw_size_max_half:*,data_raw_size_max_half:*)
; 2/ BR - bottom right corner
data_cor(data_raw_size_max_half:*,0:data_raw_size_max_half-1)=$
   res_real(0:data_raw_size_max_half-1,data_raw_size_max_half:*)
; 3/ UL - upper left corner
data_cor(0:data_raw_size_max_half-1,data_raw_size_max_half:*)=$
   res_real(data_raw_size_max_half:*,0:data_raw_size_max_half-1)
; 4/ UR - upper right corner
data_cor(data_raw_size_max_half:*,data_raw_size_max_half:*)=$
   res_real(0:data_raw_size_max_half-1,0:data_raw_size_max_half-1)

; normalization of the cedoncolved data
data_cor=data_cor/float(data_raw_size_max)/float(data_raw_size_max)

; temporal data keeping inverse action
if data_raw_size(1) ne data_raw_size(2) then begin
   ; carrousel
   data_raw=data_raw_stored
   ; trucation (case of an rectangular data)
   data_cor=data_cor(x1:x2,y1:y2)
endif

return
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;