;PRO plotaz2,cam,filt,exp
!order=0
username=''
SPAWN,'whoami',username
IF username(0) NE 'cnsr3' THEN BEGIN
	PRINT,'User is not cnsr3.'
	PRINT,'Please change users and try again'
	retall
ENDIF

PRINT,'enter the camera, filter, exposure time and day of month of exposure'
cam=0 & filt=0 & exp=0
read,cam,filt,exp,day
PRINT,'enter the LBS angle(-60,-40,-20,0,20,40,60)'
read,lbsang

;procedure to plot an image in terms of azimuth and elevation and to find 
;the azimuth and elevation of the centre of the LBS of a flat field image

;read in appropriate flat field images
path=STRING(cam,filt,FORM='("/jasper/cnsr3_data1/flat/c",I1,"/f",I1)')
cd,path
PRINT,'Enter the number of images to be averaged'
num=0
READ,num
print,'Enter the extension numbers of the flat field images you wish to use'
e=INTARR(num)
READ,e
itot=FLTARR(256,256)
;take average of all images
FOR i=0,num-1 DO BEGIN
	file=STRING(day,cam,filt,exp,e(i),FORM='("f",I2.2,2I1,I3.3,".",I3.3)')
	rdkimg,file,hb,img
	itot=itot+img
ENDFOR
kh=gethd(hb)
date=kh.misc.dt.(0)
img=itot/num
itot=FLTARR(256,256)

;read in dark frames (two)
path2=STRING(cam,FORM='("~cnsr3/dk/c",I1)')
cd,path2
;spawn,'ls d06*'

PRINT,'Enter the number of dark frames to be averaged'
num2=0
READ,num2

IF num2 GT 0 THEN BEGIN
    print,'Enter the extension numbers of the dark frame images you wish to use'
    e=INTARR(num2)
    READ,e
    tdk=INTARR(256,256)
;take mean of dark frames
    FOR i=0,num2-1 DO BEGIN
	file2=STRING(cam,day,cam,filt,exp,e(i),$
		FORM='("/jasper/cnsr3_data1/dk/c",I1,"/d",I2.2,2I1,I3.3,".",I3.3)')	
	rdkimg,file2,hbd,dk
	tdk=tdk+dk
        ENDFOR
    dk=tdk/num2
    tdk=''
    ENDIF ELSE rdmeandk,cam,exp,dk
;subtract mean dark frame from mean image
img=img-dk
dk=0
;create maps of the azimuth and elevation as done in plotorb.pro
l=256.
col=REBIN(LINDGEN(l),l,l)
row=TRANSPOSE(col)
meanf = STRING(cam,filt,FORMAT='("/wesson/user/steele/schee/dat/c",I1,"f",I1,".meanparms")')
;files '*.meanparms' contain factors to fit azimuth and zenith angle numbers
;to row and collumn numbers on an image
fl=findfile(meanf)
IF fl(0) EQ '' THEN BEGIN ;in case no .meanparms exists for this camera/filter.
	meanf=STRING(cam,FORM='("/wesson/user/steele/schee/dat/c",I1,"f*.meanparms")')
	fl=findfile(meanf)
	meanf=fl(0)
ENDIF
print,'meanf=',meanf
OPENR,u,meanf,/GET_LUN
c0=0. & r0=0. & r=0. & k=0. & camaz=0.
READF,u,c0,r0,r,k,camaz ;fitting factors
CLOSE,u & FREE_LUN,u

;c0=130.302
;r0=126.375
;r=199.895
;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
el=90.-!RADEG*(ASIN(rho/r)/k)
za=!RADEG*(ASIN(rho/r)/k)
az=!RADEG*ATAN(-y,x)
wn=WHERE(az LT 0.)
az(wn)=az(wn)+360.
wn=0
xang=!RADEG*ATAN(x,-y)
;save memory space by converting large arrays to small strings
rho='' & phi='' & yp='' & xp='' & x='' & y=''
dc='' & dr='' & col='' & row=''
loadct,13,/silent
;IF !D.NAME EQ 'X' THEN WINDOW,0

;create image only of pixels with LBS
;p is the part if the image with the LBS
p=WHERE(img GT MAX(img)/1.5)
;IF lbsang NE 0 THEN ip=WHERE(img0za GT MAX(img)/1.3)$
;	ELSE ip=p 

one2two,p,img,cp,rp
yang=ABS(lbsang*0.958863+0.338735)
IF lbsang LT 0 THEN zang=-123 ELSE zang=57



;create sub-image from flat field image of only LBS
si=img(MIN(cp,MAX=cm):cm,MIN(rp,MAX=rm):rm)
cp=0 & rp=0
;conversion of original image into azimuth and elevation coords.
;TRIANGULATE,az(p),el(p),triangles
;azmin=min(az(p))-10 & azmax=max(az(p))+10
;elmin=min(el(p))-10 & elmax=max(el(p))+10
;part of image is lost if the range exceeds 360 Degrees - set azimuth range to
;360 degrees if it presently exceeds 360 degrees
;IF (azmax-azmin) GT 360. THEN BEGIN
	;azmin=0D 
	;azmax=360D
;ENDIF

;triangulated image
;ti=TRIGRID(az(p),el(p),img(p),triangles,[.25,.25],[azmin,elmin,azmax,elmax])
;MAP_SET,(elmin+elmax)/2.,(azmin+azmax)/2.,0.,/ORTHOGRAPHIC,$
	;LIMIT=[elmin-05,azmin-05,elmax+05,azmax+05]

;mapped image in terms of azimuth and zenith angle, plotted on a latitude and
;longitude grid
;mi=MAP_IMAGE(ti,sx,sy,xsiz,ysiz,LATMIN=elmin-05,LATMAX=elmax+05,$
	;LONMIN=azmin-05,LONMAX=azmax+05,MISSING=255B)
;bmi=BYTSCL(mi,MIN=0,TOP=254)
mi=''


;make background white for a Post Script file
IF !D.NAME EQ 'PS' THEN BEGIN
    DEVICE,BITS_PER_PIXEL=8
    TVLCT,255B,255B,255B,!D.N_COLORS-1
    bmi(WHERE(bmi LE 9B))=255B
ENDIF
;MAP_GRID,LATDEL=5,LONDEL=5,/LABEL
;centre azimuth and elevation coords.
;cel=(elmax+elmin)/2. & caz=(azmax+azmin)/2.
;r=20.
;find centre row and collumn of LBS on original image
WSET,0
IF lbsang EQ 0 THEN meancent,img,x0,y0
;WINDOW,0
;find centre azimuth and elevation from azimuth/elevation image
;cirplot,cel,caz,r
;arcs,100,0,360,128,128,/dev
;PRINT,'CIRCLE IN GOOD POSITION? ([Y]/N)'
;res=GET_KBRD(1)
;WHILE(STRUPCASE(res) EQ 'N') DO BEGIN
;	PRINT,'Enter new center azimuth, elevation, and radius'
;	PRINT,'CURRENT: (DEGREES)'
;	PRINT,'azimuth = ',caz
;	PRINT,'elevation = ',cel
;	PRINT,'radius = ',r
;	READ,caz,cel,r
;	TVSCL,bmi,sx,sy,XSIZE=xsiz,YSIZE=ysiz
;	MAP_GRID,LATDEL=5,LONDEL=5,/LABEL
;	cirplot,cel,caz,r
;	PRINT,'CIRCLE IN GOOD POSITION? ([Y]/N)'
;	res=GET_KBRD(1)
;ENDWHILE
;MAP_SET,cel,caz,0.,/STEREOGRAPHIC,LIMIT=[cel-r-5,caz-r-5,cel+r-5,caz+r-5]
;mi2=MAP_IMAGE(ti,sx,sy,xsiz,ysiz,LATMIN=cel-r-5,LATMAX=cel+r-5,$
;	LONMIN=caz-r-5,LONMAX=caz+r-5,MISSING=255B)
;bmi2=BYTSCL(mi2,MIN=0,TOP=254)
;TVSCL,bmi2,sx,sy,XSIZE=xsiz,YSIZE=ysiz
;MAP_GRID,LATDEL=5,LONDEL=5,/LABEL
;cirplot,cel,caz,r

; display LBS image inside border
fimg=img
mp=MAX(img)
fimg(0,*)=mp
fimg(255,*)=mp
fimg(*,0)=mp
fimg(*,255)=mp
TVSCL,fimg,0
TVSCL,fimg,0
XYOUTS,128,490,"Dark-corrected image",/DEVICE,ALIGNMENT=0.5
fimg=''

;create image to normalize the LBS image only at 0 zenith angle
IF lbsang EQ 0 THEN BEGIN
    ; create a normalizing map
    flatnorm,az,za,img,file,x0,y0
    yang=0.
    zang=0.
    ENDIF ELSE $	; find out where the off-zenith LBS is
    rotflat,img,x0,y0,az,za,xang,zang,yang
normrot,img,p,az,el,zang,yang,cam,filt,lbsang,date

END