pro persefona,num,sym,right

; Eclipse observation: determination of the PSF from the
;                      profile of the Moon's limb.
; Philosophy: the derivative of the Moon's limb profile
;             is equal to the "compressed" PSF.
;
; Creation of array of derivatives in selected regions of
; eclipse image. The input image is first rotated, to have
; the moon limb vertical. Then a derivative image is computed
; (function "derivatives" from Oslo library), and
; the regions of interest are interactively selected. An array
; of selected regions is constructed, and each row is shifted
; with sub-pixel precision (Fourier interpolation) to have a
; sample just at the maximum. Then each maximum is shifted to
; the position of the maximum of the first row, so that the
; Moon's limb irregularities are compensated. Finally, this
; "rectified" derivative image is averaged in rows to obtain
; the mean "PSF". Result: the symmetrized "PSF", normalized to
; its area (an array of 1024 elements).
;
; Parameters: 1. serial number of the observed image (input)
;                image form: fltarr(1315,1035)
;             2. symmetrized "PSF" (output)
;             3. right wing of the "PSF", containing the
;                granular noise, for comparison (output).
; Note: The Moon's limb must cover at least 527 pixels
;       in the rotated image!
;       Some routines from ASP and Oslo libs are used.
;

if n_params() ne 3 then begin
   print,'usage: persefona, imag.number, output PSF name, output right part'
   stop
endif

; !!!! TO BE CHANGED WITH RESPECT TO WHERE THE IMAGES ARE STORED !!!!
stra=string('/scratch/sob/')
strb=string('imc.')
; !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

ysel=intarr(2,20)  ; array containing low. and upp. limits of selected areas

window,0,xsize=657,ysize=517
           
; reading the input image

   ima=fltarr(1315,1035)
   dcn=strtrim(num,2)
   if strlen(dcn) eq 3 then str='0'+dcn
   if strlen(dcn) eq 2 then str='00'+dcn
   if strlen(dcn) eq 1 then str='000'+dcn
   print,'Now reading image:'
   print,strb+str
   openr,1,stra+strb+str
   readu,1,ima
   close,1

; rotating the image

   read,'rotation angle? ',alpha
   ima=rot_int(ima,alpha)
   im2=rebin(ima(0:1313,0:1033),657,517)
   tvscl,im2

; selection  of area to compute derivatives

again:
   wset,0
   wdelete,1
   xloadct
   print,'set vertical limits of area to determine the PSF'
cur1:
   cursor,x1,y1,/device,/down
   print,y1
   print,'accept this selection? (y/n)'
   if get_kbrd(1) eq 'n' then goto,cur1
cur2:
   cursor,x2,y2,/device,/down
   print,y2
   print,'accept this selection? (y/n)'
   if get_kbrd(1) eq 'n' then goto,cur2

   y1=y1*2
   y2=y2*2
   x=x1+x2              ; (x1+x2)/2 *2

; check if we have the Moon long enough
   if x lt 527 then begin
	print,' SORRY, but the distance covered by Moon is too short.'
	stop
   endif
   print,' OK'

   if y2 lt y1 then begin   ; changing order if y2<y1
	y=y2
	y2=y1
	y1=y
   endif

; "cleaning" the edges of the image

   im=ima(*,y1:y2)
   for i=0,y2-y1 do begin
	me=total(im(5:14,i))/10.	; left edge	
	im(0:4,i)=me
	me=total(im(1300:1309,i))/10.	; right edge
	im(1310:1314,i)=me
   endfor

; computing derivatives

   dero=derivatives(im)
   window,1,xsize=1315,ysize=y2-y1+1,xpos=10,ypos=500
   tvscl,dero

; selection of areas with PSF unaffected by granulation
   
   n=0     ;counter of selected areas
newbox:
   print,'set the center position for zoom'
   cursor,x0,y0,/device,/down
   print,x0,y0

   dx=20       ; setting parameters of zoom window
   dy=30
   fact=10.

   uly=y0+dy   ; the following block assures that the zoom is...
   lly=y0-dy   ; ...always inside the array 
   if uly gt y2-y1 then begin
	y0=(y2-y1)-dy
        print,'top limit: new center position for zoom ',y0
   endif
   if lly lt 0     then begin
	y0=dy
	print,'bottom limit: new center position for zoom ',y0
   endif

   box=dero(x0-dx:x0+dx,y0-dy:y0+dy)   ; zoom
   fact=10.
   box=rebin(box,fact*(2*dx+1),fact*(2*dy+1),/sample)
   window,2,xsize=fact*(2*dx+1),ysize=fact*(2*dy+1)
   tvscl,box

;      interactive inspection/selection
   ss=-1
   ses=0
curs:
   print,'profiles - left; selection - middle; quit - right (mouse)
   cursor,x,y,/device,/down
   case !err of
	1:begin
		window,3,xsize=400,ysize=400,xpos=700
		plot,box(*,y),yrange=[-5,30],/yst
		wset,2
		goto,curs
	  end
	2:begin
		ss=ss*(-1)
		case ss of
			1:begin
				ys1=y0-dy+fix(y/fact) ; coord. in dero
				print,'1st position: ',ys1
			  end
		       -1:begin
				ys2=y0-dy+fix(y/fact) ; coord in dero
				print,'2nd position: ',ys2
			  end
		endcase
		ses=ses+1
		goto,curs
	  end
	4:begin
		if ses eq 0 then begin
			print,'you did not select anything!'
			goto,curs
		endif
		goto,continue
	  end
   endcase

continue:
   if ys2 lt ys1 then begin
	y=ys2
	ys2=ys1
	ys1=y
   endif
   print,'definitive limits of selected range in derivative image:'
   print,ys1,ys2

   ysel(0,n)=ys1
   ysel(1,n)=ys2
   n=n+1          ; counter of selected areas (max 20 areas)
   if n gt 19 then goto,finish
   print,'new box? (y/n)'
   if get_kbrd(1) eq 'y' then begin
        wset,1
	wdelete,2
	goto,newbox
   endif

; end of selection
  
finish:

   arr=fltarr(1315,y2-y1+1)   ; array of selected regions
   nsh=0   		      ; counter of rows in arr

; filling the array

   for i=0,n-1 do begin				;loop in areas
	for j=ysel(0,i),ysel(1,i) do begin	;loop in rows
		arr(*,nsh)=dero(*,j)
		nsh=nsh+1	
	endfor
   endfor

   nsh=nsh-1		      ; last row in arr
   arr=arr(*,0:nsh)

   wset,1
   wdelete,2
   wdelete,3
   erase
   tvscl,arr

   print,'Are you satisfied with your selections?'
   print,'If not, press "a" to repeat it. If yes, press any other key.'
   pau=get_kbrd(1)
   if pau eq 'a' then goto,again
   print,'OK'

; ---------------------------------------

out=arr

; sub-pixel Fourier shift to have a sample in the maximum position

for i=0,nsh do begin	
	mx=maxpos(arr(*,i),imax)
	del=imax-mx
	out(0:1023,i)=fftshift(arr(0:1023,i),del,25,998)
endfor

ma=max(out(0:1023,0),ref)          ; reference "ref"
print,'reference position of maximum:',ref

; shift in pixels to the reference position

for i=1,nsh do begin
	ma=max(out(0:1023,i),pos)
	sh=ref-pos
	out(*,i)=shift(out(*,i),sh)
endfor

; normalization to unit area (only left part considered)

for i=0,nsh do begin
	norm=out(ref,i)+2*total(out(ref-200:ref-1,i))
	out(*,i)=out(*,i)/norm
endfor
tvscl,out
print,'This is the shifted derivative image. Press any key to continue'
pau=get_kbrd(1)

sum=total(out,2)/(nsh+1)
right=sum(ref:1023)
sym=[reverse(sum(ref-511:ref)),sum(ref-512:ref-1)]

plot,sum(ref:ref+20),linestyle=5
oplot,sym(0:20)
print,'This is the result of the symmetrization'

;final normalization of sym

sym=sym/total(sym)

end