PRINT,'enter the camera and filter numbers'
READ,cam,filt

img= FLTARR(256,256,31)
i=FLTARR(256,256)

file1=STRING(cam,filt,$
	FORM='("/jasper/cnsr3_data1/flat/fn??",2I1,"*")')
fl1=FINDFILE(file1,count=nfl1)

; Make sure these files are up-to-date, i.e., created in October 93.
ok=BYTARR(nfl1)
FOR j=0,nfl1-1 DO BEGIN
    SPAWN,'ls -l '+fl1(j),result
    ok(j)=(STRPOS(result,'Oct') GT 0)
    ENDFOR
fl1=fl1(WHERE(ok))
nfl1=N_ELEMENTS(fl1)

FOR k=0,nfl1-1 DO BEGIN
	rdnflat,fl1(k),im
	img(*,*,k)=im
ENDFOR

FOR j=nfl1-1,0,-1 DO BEGIN
w=WHERE(img(*,*,j))
i(w)=img(w+65536*j)
ENDFOR

WINDOW,2,XSIZE=256,YSIZE=256
loadct,13
TVSCL,i

;dc0=profile(i,xx,yy)
;IF xx(0) GT xx(50) THEN xx=ROTATE(xx,2)
;IF yy(0) GT yy(50) THEN BEGIN 
;	yy=ROTATE(yy,2)
;	xx=ROTATE(xx,2)
;ENDIF
;two2one,xx,yy,i0,w

;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
OPENR,u,meanf,/GET_LUN
k=0.
READF,u,c0,r0,r,k,camaz ;fitting factors
CLOSE,u & FREE_LUN,u
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=90.-el
az=!RADEG*ATAN(-y,x)
az(WHERE(az LT 0))=az(WHERE(az LT 0))+360

eln=el*((az GT 90) AND (az LE 270))
elp=el*((az LE 90) OR (az GT 270))
elp(WHERE(elp))=180-elp(WHERE(elp))
el(WHERE(elp))=elp(WHERE(elp))


mk,prof
prof(127:129,*)=1
w=WHERE(prof)

PRINT,'Would you like the y-axis to go to zero?'
zer=''
READ,zer

PRINT,'Would you like to create a PS file ? (y/[n])
ans=''
READ,ans
IF STRLOWCASE(ans) EQ 'y' THEN BEGIN
    PRINT,'Name the PS file'
    psname=''
    READ,psname
    psfile="/jasper/cnsr3_data1/flat/"+STRTRIM(psname,2)
    psopen,psfile,/land
    !P.MULTI=[0,4,3,0,0]
    ;device,yoffset=1.905
    ENDIF ELSE BEGIN
    WINDOW,0,XSIZE=975,YSIZE=850
    psname=STRING(cam,filt,FORM='("c",I1,"/f",I1)')
    WSET,0
    IF !p.multi(1) NE 2 THEN !P.MULTI=[0,2,2,0,0]
    ENDELSE

;making profiles of the fov along lines of constant azimuth, every 15 degrees.
maxdn=FLTARR(180)
mindn=maxdn+50000
numtot=maxdn & totaldn=maxdn
elarr=FLTARR(1800,12) & dnarr=elarr

FOR n=3,168,15 DO BEGIN
    eltot=0 & dntot=0
    rprofile=rot(prof,n,missing=0)
    w=WHERE(rprofile)
    ;one2two,w,profile,xx,yy
    ;yf=1
    ;fit=svdfit(xx,yy,2,yfit=yf)
    ;xx=INDGEN(256)
    ;yy=xx*fit(1)+fit(0)
    ;wy=WHERE((yy GE 0) AND (yy LE 255))
    ;xx=xx(wy)
    ;yy=yy(wy)
    ;yy=ROUND(yy)
    ;mk,profile
    ;profile(xx,yy)=1
    FOR m=0,nfl1-1 DO BEGIN
        profile2=img(*,*,m)*rprofile
        wp=WHERE(profile2)
        IF wp(0) NE -1 THEN BEGIN
            eltot=[eltot,el(wp)]
            dntot=[dntot,img(wp+m*65536)]
            smaldrop=WHERE(ABS(dntot-SHIFT(dntot,1))/$
                          (ABS(eltot-SHIFT(eltot,1))>1) lt 500)
            et=eltot(smaldrop)
            dt=dntot(smaldrop)
            smaldrop=WHERE(ABS(dt-SHIFT(dt,1))/$
                          (ABS(eltot-SHIFT(eltot,1))>1) lt 500)
            eltot=et(smaldrop)
            dntot=dt(smaldrop)

        ENDIF
    ENDFOR
    dntot=dntot(1:*) & eltot=eltot(1:*)

    ;ptitle=STRING(cam,filt,n,$
       ;FORM='("c",i1,"/f",i1,"  ",i3.3,"degrees azimuth")')
    ;plot,eltot,dntot,psym=3,XTITLE='elevation',YTITLE='DN',TITLE=ptitle

    esort=sort(eltot)
    elsort=eltot(esort)
    dnsort=dntot(esort)
    smaldrop=WHERE(ABS(dnsort-SHIFT(dnsort,1))/$
                  (ABS(elsort-SHIFT(elsort,1))>1) lt 1000)
    eltot=elsort(smaldrop) & dntot=dnsort(smaldrop)
    ;plotting a minimum and maximum envelope onto each graph
    FOR j=FIX(MIN(eltot)),ROUND(MAX(eltot)) DO BEGIN
        wel=WHERE(FIX(eltot) EQ j)
        IF wel(0) NE -1 THEN BEGIN
	    jj=j-FIX(MIN(eltot))
            IF MIN(dntot(wel)) LT mindn(jj) THEN mindn(jj)=MIN(dntot(wel))
            IF MAX(dntot(wel)) GT maxdn(jj) THEN maxdn(jj)=MAX(dntot(wel))

            totaldn(jj)=totaldn(jj)+TOTAL(dntot(wel))

            numtot(jj)=numtot(jj) + N_ELEMENTS(wel)
            ENDIF
        ENDFOR
    elarr(0:N_ELEMENTS(eltot)-1,n/15)=eltot
    dnarr(0:N_ELEMENTS(dntot)-1,n/15)=dntot
    ENDFOR
avdn=totaldn(WHERE(numtot))/numtot(WHERE(numtot))

FOR n=0,11 DO BEGIN
    elevation=elarr(*,n)
    elevation=elevation(WHERE(elevation))
    dn=dnarr(*,n)
    dn=dn(WHERE(dn))
    ptitle=STRING(cam,filt,n*15+3 $
                 ,FORMAT='("c",i1,"/f",i1,"  ",i3.3," degrees azimuth")')
    IF STRLOWCASE(zer) EQ 'n' THEN $
        PLOT,elevation,dn,PSYM=3,XTITLE='elevation',YTITLE='DN',TITLE=ptitle $
    	    ,yrange=[MIN(dnarr(WHERE(dnarr)),MAX=mxdn),mxdn]

    IF STRLOWCASE(zer) EQ 'y' THEN $
        PLOT,elevation,dn,PSYM=3,XTITLE='elevation',YTITLE='DN',TITLE=ptitle
    xx=INDGEN(ROUND(MAX(eltot))-FIX(MIN(eltot)))+FIX(MIN(eltot))
    wm=WHERE(mindn NE 50000) & wmx=WHERE(maxdn)
    mindn=mindn(wm) & xx1=xx(wm) & maxdn=maxdn(wmx) & xx2=xx(wmx)
    OPLOT,xx1,mindn
    OPLOT,xx2,maxdn
    OPLOT,xx1,avdn
    IF !D.(0) EQ 'X' THEN BEGIN
	PRINT,'Press any key to continue'
	g=GET_KBRD(1)
	ENDIF
    ENDFOR
IF STRLOWCASE(ans) EQ 'y' THEN psclose
!P.MULTI=0

; try to generate a 'mean' image from the 31 overlapping pieces
mimg=FLTARR(256,256) & nimg=mimg
FOR j=0,nfl1-1 DO BEGIN
    mimg=mimg+img(*,*,j)
    nimg=nimg+(img(*,*,j) GT 0)
    ENDFOR
mimg=mimg/(nimg>1.)
WSET,2
TVSCL,mimg

; save 'mean' image along with corresponding azimuth and elevation maps
el=90.-za
outname=STRING(cam,filt $
	      ,FORMAT='("~cnsr3/flat/ffimg.",2I1,".raw")')
PRINT,'Saving mean image with azimuth and elevation maps in '+outname
OPENW,out,outname,/GET_LUN
sm=SIZE(mimg)
WRITEU,out,sm(1),sm(2)
WRITEU,out,mimg
WRITEU,out,az
WRITEU,out,el
CLOSE,out
FREE_LUN,out

END