FUNCTION Go_On
Print,'Press any key to go on'
k=Get_Kbrd(1)
Return,(StrUpCase(k) NE 'Q')
END

;+
; NAME:
;	CalComp
; PURPOSE:
;	Comparison of calibration images from different dates
; CATEGORY:
;	Data validation
; CALLING SEQUENCE:
;	CalComp, File1, File2 [, Dark1 = dfile1, Dark2 = dfile2 ]
; INPUTS:
;	File1:  The name (only; no path required) of the (presumed valid)
;               calibration image file.
;       File2:  The name (only; as above) of another calibration image
;               file acquired with the same camera, filter, and exposure
;               time.
; OPTIONAL INPUTS:
;	None.
; KEYWORD PARAMETERS:
;	Dark1:  If supplied, gives the name of the appropriate dark
;               frame to subtract from the image of File1.
;       Dark2:  As above, for File2.
; OUTPUTS:
;	None.
; OPTIONAL OUTPUTS:
;	None.
; COMMON BLOCKS:
;	None.
; SIDE EFFECTS:
;	The specified files are opened and read.  Various plots are
;       created on the screen.
; RESTRICTIONS:
;	The specified files must exist, and contain valid data.
; PROCEDURE:
;	Fairly straightforward, except the statistical calculations at
;       the end, for which references are given.
; EXAMPLE:
;	
; SEE ALSO:
;	
; MODIFICATION HISTORY:
; 	Written by:	DPSteele, ISAS, June 4, 1996.
;-
PRO CalComp,file1,file2,dark1=dk1,dark2=dk2

@isitdos

IF !D.Window NE -1 THEN WDelete,!D.Window

path=iroot+dd+imagetype(2)+dd
cam=StrMid(file1,3,1)
filt=StrMid(file1,4,1)
time=StrMid(file1,5,3)
mismatch=(cam NE StrMid(file2,3,1)) OR $
         (filt NE StrMid(file2,4,1)) OR $
         (time NE StrMid(file2,5,3))
IF mismatch THEN Message,'Image file mismatch'
path=path+cc+cam+dd+ff+filt+dd
Help,path
IF NOT Go_On() THEN Return

IF (N_Elements(dk1) GT 0) OR (N_Elements(dk2) GT 0) THEN BEGIN
    dpath=iroot+dd+imagetype(3)+dd+cc+cam+dd
    Help,dpath
    IF NOT Go_On() THEN Return
    ENDIF
    
RdKImg,path+file1,hb1,im1
RdKImg,path+file2,hb2,im2
h1=GetHd(hb1)
h2=GetHd(hb2)
Print,'Image 1:'
Help,/St,h1.misc.dt
Help,/St,h1.misc.tm
Print,'Image 2:'
Help,/St,h2.misc.dt
Help,/St,h2.misc.tm
IF NOT Go_On() THEN Return
sim1=StdSize(im1)
sim2=StdSize(im2)

IF N_Elements(dpath) GT 0 THEN BEGIN
    IF N_Elements(dk1) EQ 0 $
        THEN RdMeanDk,Fix(cam),Fix(time),dim1,/NoInteractive $
        ELSE RdKImg,dpath+dk1,dhb1,dim1
    dim1=StdSize(dim1)
    sim1=(sim1-dim1)>0
    IF N_Elements(dk2) EQ 0 $
        THEN RdMeanDk,Fix(cam),Fix(time),dim2,/NoInteractive $
        ELSE RdKImg,dpath+dk2,dhb2,dim2
    dim2=StdSize(dim2)
    sim2=(sim2-dim2)>0
    ENDIF
    
pc1=Pcentile(sim1,[0.01,0.99])
pc2=Pcentile(sim2,[0.01,0.99])
Print,'Image 1 in standard size: 1st and 99th percentiles'
Print,pc1
Print,'Image 2 in standard size: 1st and 99th percentiles'
Print,pc2
Window,0,XSize=512,YSize=512
TVScl,sim1>pc1(0)<pc1(1),0
TVScl,sim2>pc2(0)<pc2(1),1
IF NOT Go_On() THEN Return

MeanCent,sim1,cc1,cr1,tot1,Radius=rad1,/NoTV
MeanCent,sim2,cc2,cr2,tot2,Radius=rad2,/NoTV
csh=Round(cc1-cc2)
rsh=Round(cr1-cr2)
Print,csh,rsh $
     ,Format='("Shift image 2 ",I0," columns and ",I0,' $
            +'" rows to register PLBS images?")'
IF NOT Go_On() THEN Return
sim2=Shift(sim2,csh,rsh)

resp=''
Read,'Rotate image 2 to match image 1 better? ([Y]/N) ',resp
IF NOT (StrUpCase(resp) EQ 'N') THEN BEGIN
    Print,'Rotating in 30 deg steps ...'
    ang=0
    REPEAT BEGIN
        ang=ang+30
        TVScl,Rot(sim2,ang,1,cc1,cr1,/Pivot)>pc2(0)<pc2(1),2
        Read,'Is this satisfactory? (Y/[N]) ',resp
        ENDREP UNTIL (StrUpCase(resp) EQ 'Y')
    sim2=Rot(sim2,ang,1,cc1,cr1,/Pivot)
    TVScl,sim2>pc2(0)<pc2(1),1
    ENDIF
IF NOT Go_On() THEN Return

Print,'Generating weighting arrays for images ...'
spc1=(1.-Float(tot1)/N_Elements(sim1))  ; fraction of image pixels 
spc2=(1.-Float(tot2)/N_Elements(sim2))  ; not showing PLBS screen
wt1=(Float(sim1)/Pcentile(sim1,spc1))<1
wt2=(Float(sim2)/Pcentile(sim2,spc2))<1
TVScl,wt1,2
TVScl,wt2,3
k=Where((wt1 GT 0.5) AND (wt2 GT 0.5),nk)
Print,nk,Format='(I0," key points common to both images")'
IF NOT Go_On() THEN Return

Print,'Evaluating statistics ...'

; Note: Formulae used below are from Bevington & Robinson, Data
; Reduction and Error Analysis for the Physical Sciences, 2nd ed.,
; McGraw-Hill, 1992, pp. 50, 60.

mean1=Total(wt1(k)*sim1(k))/Total(wt1(k))
mean2=Total(wt2(k)*sim2(k))/Total(wt2(k))
var1=((Total(wt1(k)*sim1(k)*sim1(k))/Total(wt1(k)))-mean1*mean1) $
    *(nk/(nk-1))
var2=((Total(wt2(k)*sim2(k)*sim2(k))/Total(wt2(k)))-mean2*mean2) $
    *(nk/(nk-1))
cov=(Total(Sqrt(wt1(k))*(sim1(k)-mean1)*Sqrt(wt2(k))*(sim2(k)-mean2)) $
    /Total(Sqrt(wt1(k)*wt2(k))))*(nk/(nk-1))
Help,mean1,var1,mean2,var2,cov

ss=Screen_Size()
goldrat=(1+Sqrt(5))/2
Window,4,XSize=ss(1)/goldrat,YSize=ss(1)
!P.Multi=[0,0,2,0,0]
rng1=[Min(sim1(k),Max=mx),mx]
rng2=[Min(sim2(k),Max=mx),mx]

bin1=NiceNumber(0.01*(rng1(1)-rng1(0)))
x1=MakeX(rng1(0)/bin1,rng1(1)/bin1+1,1)*bin1+bin1/2
nx1=N_Elements(x1)
h1=Histogram(sim1(k),Min=rng1(0),Max=rng1(1)+bin1,Bin=bin1)
nh1=N_Elements(h1)
CASE 1 OF
nx1 LT nh1: BEGIN
        x1=[x1,(FIndGen(nh1-nx1)+1)*bin1+x1(nx1-1)]
        nx1=N_Elements(x1)
        END
nx1 GT nh1: BEGIN
        x1=x1(0:nh1-1)
        nx1=N_Elements(x1)
        END
ELSE:
    ENDCASE
; x1=x1(0:(nx1<nh1)-1)
Plot,x1,h1,PSym=10
g1=Exp(-0.5*(x1-mean1)*(x1-mean1)/var1)
f1=LADFit(g1,h1)
OPlot,x1,f1(0)+f1(1)*g1
; Print,'Best-fit Gaussian for sim1:'
; Print,f1,Format='("DC: ",F5.1,"  Peak Height: ",F6.1)'

bin2=NiceNumber(0.01*(rng2(1)-rng2(0)))
x2=MakeX(rng2(0)/bin2,rng2(1)/bin2+1,1)*bin2+bin2/2
nx2=N_Elements(x2)
h2=Histogram(sim2(k),Min=rng2(0),Max=rng2(1)+bin2,Bin=bin2)
nh2=N_Elements(h2)
CASE 1 OF
nx2 LT nh2: BEGIN
        x2=[x2,(FIndGen(nh2-nx2)+1)*bin2+x2(nx2-1)]
        nx2=N_Elements(x2)
        END
nx2 GT nh2: BEGIN
        x2=x2(0:nh2-1)
        nx2=N_Elements(x2)
        END
ELSE:
    ENDCASE
; x2=x2(0:(nx2<nh2)-1)
Plot,x2,h2,PSym=10
g2=Exp(-0.5*(x2-mean2)*(x2-mean2)/var2)

; save_multi=!P.Multi
; !P.Multi=0
; Window,/Free,XSize=512,YSize=384
; Plot,g2,h2,PSym=3

f2=LADFit(g2,h2)

; OPlot,g2,f2(0)+f2(1)*g2
; Hlp,x2
; Print,f2
; IF NOT Go_On() THEN Stop
; WDelete,!D.Window
; !P.Multi=save_multi

OPlot,x2,f2(0)+f2(1)*g2
; Print,'Best-fit Gaussian for sim2:'
; Print,f2,Format='("DC: ",F5.1,"  Peak Height: ",F6.1)'
Wait,5
!P.Multi=0
WDelete,!D.Window
IF NOT Go_On() THEN Return

Print,'Generating ratio of key points of sim2/sim1 ...'

; Note: formulae used below are from Mood, Graybill, & Boes,
; Introduction to the Theory of Statistics, 3rd ed., McGraw-Hill, 1974,
; p. 181 (Theorem 4).

meanrat=(mean2/mean1)-(cov/(mean1*mean1))+(mean2/mean1^3)*var1
tmp=((var2/(mean2^2))+(var1/(mean1^2))+(2*cov)/(mean1*mean2))
varrat=(mean2/mean1)^2*tmp
sdrat=Sqrt(varrat)
Print,meanrat,sdrat $
     ,Format='("Mean ratio is ",F5.3/"standard deviation is ",F5.3)'
TVSI,(sim2/Float(sim1))>(meanrat-sdrat)<(meanrat+sdrat) $
    ,80,80,175,175,CTable=0
; WDelete,0
Return
END