PRO flatnorm,az,za,img,file,ccol,crow

;Procedure to produce an image which will normalize an LBS image to a uniform 
;brightness based on the centre (2 degrees of zenith angle) of the LBS

;option to create Post Script file
pso='' & psfile='' 
print,'Do you wish to open a PS file? (y/[n])'
READ,pso
IF STRUPCASE(pso) EQ 'Y' THEN BEGIN
	PRINT,'Enter the name of the file'
	READ,psfile
	PSOPEN,psfile
ENDIF

;take the point found by fitting a circle 
;over the original image to find the centre row and collumn.
	c0=ccol
	r0=crow
l=256.
;create azimuth and zenith angle images where the centre is the specified 
;collumn and row
col=REBIN(LINDGEN(l),l,l)
row=TRANSPOSE(col)
r=199.895	;+((ABS(r0-126.375))<(ABS(c0-130.302)))
k=0.638308
ca=0.
dc=col-REPLICATE(c0,l,l)
dr=row-REPLICATE(r0,l,l)
rho=SQRT(dc^2+dr^2)
phi=ATAN(dc,-dr)
xp=COS(phi)
yp=SIN(phi)
cc=COS(ca)
sc=SIN(ca)
x=xp*cc+yp*sc
y=-xp*sc+yp*cc
zacent=!RADEG*(ASIN(rho/r)/k)
azcent=!RADEG*ATAN(-y,x)

;take 2 degrees of zenith angle from centre as brightest pixels
small=(zacent LT 2)*img
;mean of brightest pixels
smean=mean(small(where(small)))
;create normalizing map of LBS
p=where(img GT max(img)/1.5)
inorm=smean/img(p)

;create normalising image
normimg=FLTARR(256,256)
normimg(p)=inorm
;plot contours of normalizing image
p0=[!X.MARGIN(0)*!D.X_CH_SIZE,!Y.MARGIN(0)*!D.Y_CH_SIZE]
dpy=!D.Y_SIZE-!D.Y_CH_SIZE*TOTAL(!Y.MARGIN)
pp=[p0,p0+[dpy,dpy]]
CONTOUR,normimg,LEVELS=((FINDGEN(8)+20)/20.),/FOLLOW $
       ,TITLE=psfile,C_LABELS=BINDGEN(8),POSITION=pp,/DEVICE $
       ,XRANGE=[84,171],YRANGE=[84,171],XSTYLE=1,YSTYLE=1 $
       ,XTITLE='Column',YTITLE='Row'

openw,1,'azzatest.dat'
printf,1,n_elements(p)
printf,1,az(p)
printf,1,za(p)
printf,1,col(p)
printf,1,row(p)
close,1

;save data
PRINT,'Do you wish to save the normalizing map? ([Y]/N)'
ans=''
READ,ans
IF STRLOWCASE(ans) NE 'n' THEN BEGIN
    sp=STRPOS(file,'.')
    name=STRMID(file,sp-8,12)
    flatfile=STRING(name,FORM='("~cnsr3/flat/fnorm.",A,".dat")')
    OPENW,unit,flatfile,/GET_LUN
    WRITEU,unit,N_ELEMENTS(p)
    WRITEU,unit,az(p),za(p),inorm
    close,unit & FREE_LUN,unit
    ENDIF
IF STRUPCASE(pso) EQ 'Y' THEN psclose
PRINT,smean,FORM='("mean over 2 degrees is: ",F12.5)'

RETURN
END