PRO BPNECALC
;
; Numerically evaluates the Back Projection Normalization Error, as described in GH's 11-19-01 notebook for
; for all subcollimators, and generates a page of plots showing its dependence on radial offset.
;
; Initialize by defining vectors, etc.
pitch = 2.* 2.26 / 3600. * 3^ (0.5*INDGEN(9))						; degrees
fov = [1., 1., 1., 1., 1., 1., 4.56, 7.71, 2.81]				; degrees
ampl = 0.5
npt =5000
theta = FINDGEN(npt) * !PI / npt								; radians
nr = 500
phaseoff = !PI * 0
rtab = FINDGEN(nr) * 0.5 / nr									; degrees
z	= FLTARR(nr)
!P.MULTI = [0,3,3]
PRINT, '   SMECALC OUTPUT'
PRINT, 'SUBCOLL   MIN     MAX     +-     BIAS'
FOR nsc = 1,9 DO BEGIN
	FOR j = 0, nr-1 DO BEGIN
		r = rtab(j)													; degrees
		g = (1. - r * SIN(theta) / fov[nsc-1])^2
		phase = 2.* !PI * r * SIN(theta) / pitch(nsc-1)	+phaseoff
		p = 1. + ampl * COS(phase)
		c = g*p
		avp 	= MEAN(p)
		avg 	= MEAN(g)
		avp2 	= MEAN(P^2)
		varp = avp2-avp^2
		IF (varp EQ 0) THEN z[j] = 1. ELSE z[j]	= TOTAL(c*(p-avp)) / (avg*varp) / npt
	ENDFOR
	zmax = MAX(z[WHERE(rtab GT pitch(nsc-1))], MIN=zmin)
	zdiff = (zmax-zmin)/2.
	zbias = zmax - zdiff
	PRINT, nsc, zmin, zmax, zdiff, zbias, FORMAT='(I4, 3X, 4F8.3)'
	PLOT, rtab, z, YRANGE=[0.8, 1.2], XTITLE='SOURCE OFFSET(degrees)', YTITLE='RELATIVE NORMALIZATION',	$
			TITLE=STRING(nsc, FORMAT='("SUBCOLLIMATOR", I2)')
ENDFOR
END