PRO POStats,R,NumAcq,NumData,outfile
;
;  This program was written by Laura Ellen Dafoe on 24 November 1993.
;  It calculates the six statistical algorithms sent by Paris Observatory
;  to us 30 September 1993.  These are the same algorithms used in their
;  GSE.  It is coded here for comparison with data in the full DISR.
;
;	U	Float of the data
;	M	Array of the mean for each acquisition
;	P	Array of the mean over the acquistions for each pixel
;	PS	Array of the std dev over the acquisitions for each pixel
;
;	Tested at PM Delivery against PO data
;  Modified 29 June 1994 by LED and Lyn Doose to incorporate flight software
;  normalization FOR 6-12-6 ONLY!

U=double(R)
M=DBLARR(NumAcq)
P=DBLARR(NumData)
PS=DBLARR(NumData)
printf,outfile,'Standard Deviation of All Data Values ',stdev(U)
for k=0,NumAcq-1 do  M(k)=total(U(*,k))/double(NumData)
;
; PO Number One
;
printf,outfile,'Mean over the number of acqs of the means for a given acq, PO#1'
POOne=double(total(M))/double(NumAcq)
printf,outfile,POOne*3.
;
; PO Number Two
;
if NumData GT 1 then begin
sum=0.D
for k=0,NumAcq-1 do sum = sum + stdev(U(*,k))
POTwo = sum/double(NumAcq)
printf,outfile,'Mean over the number of acqs of the stddevs for a given acq, PO#2'
printf,outfile,POTwo*3.
endif
;
; PO Number Three
;
if NumAcq eq 1 then STOP
sum=0.D
for i=0,NumData-1 do sum = sum + stdev(U(i,*))
POThree = sum/double(NumData)
print,'Mean IR Noise Figure of Merit ',POThree*3.
printf,outfile,'Mean over the num of pixels of the stddev over the num of acqs, PO#3'
printf,outfile,POThree*3.
;
; PO Number Four
;
if NumData GT 1 then begin
sum=0.D
for i=0,NumData-1 do sum = sum + stdev(U(i,*))^2 - POThree^2
POFour = sqrt(sum/double(NumData-1))
printf,outfile,'Standard deviation of the standard deviations over the acqs, PO#4'
printf,outfile,POFour*3.
endif
printf,outfile,' '
printf,outfile,' '
;
; PO Number Six, the average of the pixels over the number of acquisitions 
;
for i=0,NumData-1 do P(i) = total(U(i,*))*3./NumAcq
;
; PO Number Five, the standard deviation of the pixels over the number of acq's
;
for i=0,NumData-1 do PS(i) = stdev(U(i,*))*3.
;
printf,outfile,'Statistics of the Pixels over the Acquisitions'
printf,outfile,'Pix Num       Mean          Standard Deviation'
for i=0,NumData-1 do begin
printf,outfile,format='(i6,4x,f10.2,10x,f10.4)',i,P(i),PS(i)
endfor

;WINDOW,1,$
;TITLE='Average of the Pixels over the Acquisitions'
;plot,P,XTITLE='Pixel Number',YTITLE='Avg',CHARSIZE=2,$
;yrange=[min(P),max(P)],ystyle=1
;WINDOW,2,$
;TITLE='Standard Deviation of the Pixels over the Acquisitions'
;plot,PS,XTITLE='Pixel Number',YTITLE='Std Dev',CHARSIZE=2,$
;yrange=[min(PS),max(PS)],ystyle=1
;
RETURN
END