pro tn95,init,numim,itr,norm

; Extracting "clear" field unaffected by shifts and rotation,
; 2-D parabolic correction for intensity trends,
; normalization of images. 
; input: files stra[init:init+numim-1], output:  files stro........
; input parameters: init = first frame no.
;		    numim = number of frames
;		    itr  = option (1 .. correct for trends)
; output parameter: norm = array of means and standard deviations

; reading the input image
; !!!! TO BE CHANGED WITH RESPECT TO WHERE THE IMAGES ARE STORED !!!!

stra=string('/scratch/sob/lap95/ima/a_rs30Jun95.')
stro=string('/scratch/sob/lap95/ima/a_n30Jun95.')

; !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

;input parameters

x1= 0       ;boundaries of the "clear" field
y1= 0
x2= 1309
y2= 1031

;xx1=375	      ;boundaries of the field to be excluded fm normalization
;yy1=370	      ;  with resp. to the "clear" field
;xx2=719
;yy2=710

xx1=60		      ;boundaries of field to be used for normalization
yy1=50
xx2=560
yy2=330

ima1=intarr(1310,1032)  ;input image array
norm=fltarr(3,numim)

for i=0,numim-1 do begin

   dcn=strtrim(i+init,2)

   print,'Now reading image:'
   print,stra+dcn
   openr,1,stra+dcn
   readu,1,ima1
   close,1

; extraction and trend correction (optional)
   ima=float(ima1(x1:x2,y1:y2))

   if itr eq 1 then begin
      sima=sfit(congrid(ima,(x2-x1+1)/3,(y2-y1+1)/3,/cubic),2)
      sima=congrid(sima,x2-x1+1,y2-y1+1,/cubic)
      ima=ima/sima	;this is the "secondary flatfielding",
			; trends are not physical, but due to wrong
			; previous flatfielding
   endif

; statistics
; a case when the subfield is not excluded but used as the normalization
; area
   std=stdev(ima(xx1:xx2,yy1:yy2),mea)

  ; un con~azo con "excluded field" called exim
;	exim=double(ima(xx1:xx2,yy1:yy2))
;	nel=n_elements(ima)-n_elements(exim)  ; n_elements taken into account
;	mea=(total(double(ima))-total(exim))/nel   ; mean
;	std=total((double(ima)-mea)^2)-total((exim-mea)^2)
;       std=sqrt(std/(nel-1))
;	mea=float(mea)
;	std=float(std)

   print,'mean =',mea,'   norm.std =',std/mea
   norm(0,i)=[i+init,mea,std/mea]

; normalization

   ima=ima/mea*10000.
   ima=nint(ima)
 
   print,'Saving image:'
   print,stro+dcn
   openw,1,stro+dcn
   writeu,1,ima
   close,1

endfor

print,'normalization terminated ok'

end