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

;array of 31 normalized flat field images to be read in
img= fltarr(256,256,31)
i=fltarr(256,256)

;files to read in. Extention number notpresent because all files with this base
;filename are needed. If there are not 31 files, the tiled image will be 
;incomplete or an error will occur in the program
file1=STRING(cam,filt,$
	FORM='("/jasper/cnsr3_data1/flat/fn??",2I1,"*")')
fl1=findfile(file1,count=nfl1)

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

IF (!D.WINDOW EQ 0) AND (!D.X_SIZE NE 975) THEN WINDOW,0,XSIZE=975,YSIZE=850
IF !D.WINDOW EQ -1 THEN WINDOW,0,XSIZE=975,YSIZE=850

WINDOW,2,XSIZE=300,YSIZE=300
loadct,13
tvscl,i

;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)
az=!RADEG*ATAN(-y,x)
az(where(az LT 0))=az(where(az LT 0))+360

;define elevations on one side of elevation map from 90-180 so that results 
;can be plotted from 0-180 degrees elevation
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)) ;el now defined 0-180


mk,prof
prof(127:129,0:128)=1 ;256X256 image with a radial line to identify the profile
w=where(prof)

psname=STRING(cam,filt,FORM='("c",I1,"/f",I1)') ;future title of graphs
WSET,0
IF !p.multi(1) NE 2 THEN !P.MULTI=[0,2,2,0,0] ;plot graphs 4 on a page

maxdn=FLTARR(180)
mindn=maxdn+50000 ;make minimum start excessively large so that later values 
		  ;to replace these will be smaller since a minimum is being
		  ;found

numtot=maxdn & totaldn=maxdn
elarr=fltarr(1800,72) & dnarr=elarr
avdn=FLTARR(180,72)
elev=avdn

;making profiles of the fov along lines of constant azimuth, every 5 degrees.
FOR n=3,358,5 DO BEGIN
	eltot=0 & dntot=0
	rprofile=ROT(prof,n,missing=0)
	w=where(rprofile)
	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)]

	      ;take away points where the slope is very steep. These points are
	      ;at the edges of the normalized LBS where things don't quite
	      ;line up
	      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
	;first element of dntot & eltot was 0 (when the were defined)
	;take away first element
	dntot=dntot(1:*) & eltot=eltot(1:*)

	;sort elevations into ascending order, then take away the sharp slopes
	;between points again
	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

	;loop from closest integer below minimum elevation seen so far, to 
	;integer closest to max elevation
	FOR j=FIX(min(eltot)),ROUND(max(eltot)) DO BEGIN 
	    wel=where(FIX(eltot) EQ j)
	    IF wel(0) NE -1 THEN BEGIN
		;find minimum and maximum dn at elevation 'j'
		IF (min(dntot(wel)) LT mindn(j-FIX(min(eltot)))) THEN $
	    	   mindn(j-FIX(min(eltot)))=min(dntot(wel))
		IF (max(dntot(wel)) GT maxdn(j-FIX(min(eltot)))) THEN $
	    	   maxdn(j-FIX(min(eltot)))=max(dntot(wel))

		;sum of all dn at this elevation - used to find mean later
		totaldn(j)=totaldn(j)+$
			TOTAL(dntot(wel))

		;total number of points at this elevation - used to find mean
		numtot(j)=numtot(j)+$
			N_ELEMENTS(wel)
	    ENDIF
	ENDFOR

	;mean DN at each elevation
	avdn(WHERE(numtot)+(n/5)*180)=$
		totaldn(where(numtot))/numtot(where(numtot))

	;place the array of dn and elevations at this profile line azimuth
	;into a 3-dimensional array, with the 3rd dimension representing each
	;azimuth angle of the profile
	elarr(0:N_ELEMENTS(eltot)-1,n/5)=eltot
	dnarr(0:N_ELEMENTS(dntot)-1,n/5)=dntot

	;re-initialize variables so they can be re-summed for the next azimuth
	totaldn=FLTARR(180)
	numtot=totaldn
ENDFOR

avtot=FLTARR(180) & sdev=avtot
angl=intarr(180)
FOR n=0,179 DO BEGIN
	wav=WHERE(avdn(n,*))
	sz=SIZE(wav)

	;create an average dn at every elevation where any points were
	;present on the profile
	IF wav(0) NE -1 THEN BEGIN
	    avtot(n)=TOTAL(avdn(n,wav))/$
		    N_elements(avdn(n,wav))
	    IF sz(1) GT 1 THEN sdev(n)=STDEV(avdn(n,wav))
	ENDIF
ENDFOR

;angl= values from 3 to 358 in 5 degree increments
angl='az = '+STRTRIM(INDGEN(72)*5+3,2)


PRINT,'Do you want a PS plot of the total mean vs. elevation (one graph only)?'
ans=''
READ,ans
;plot of average dn and standard dev. vs. elevation
IF strlowcase(ans) EQ 'y' THEN BEGIN
	filename=STRING(cam,filt,$
		FORM='("/wesson/user/steele/schee/dat/meanstdev.",2I1)')
	psname=STRTRIM(filename,2)+'.ps'
	psopen,psname
	PRINT,'File is called: ',psname	
	tname=STRING(cam,filt,FORM='("c",I1,"/f",I1)')
	!P.MULTI=0
	plot,avtot,title=tname,xtitle='elevation',$
		ytitle='DN, standard dev. * 10'
	oplot,sdev*10,line=1
	legend,['mean','stdev*10'],linestyle=[0,1] ;used to identify lines on plot
	psclose
	dname=STRTRIM(filename,2)+'.dat'
ENDIF

;print average dn and standard deviation values into a file
OPENW,unit,dname,/GET_LUN
printf,unit,180
printf,unit,avtot
printf,unit,sdev
close,unit & FREE_LUN,unit
For b=0,20 do bell


WDELETE,2

END