;
; photometry.pro
;   support functions for photometry
;

pro calcBandData,bandName,bandData
	bandName = ['B','V','R','I','J','H','Ks']
	bandData = [[0.4413,0.0959], $
	            [0.5512,0.0891], $
	            [0.6471,0.1436], $
	            [0.8060,0.1495], $
	            [1.2535,0.1531], $
	            [1.6343,0.2908], $
	            [2.1450,0.3173]]
end
; -- calculate the number of photons in a given filter passband, in photons/s/m^2 --
function photoCounts,band,m
	c = 3.0e8 ; m/s
	h = 6.63e-34 ; Js
	k = 1.38e-23 ; J/K
	microns = 1.e-6
	calcBandData,bandName,bandData
	i = (where(band eq bandName))[0]
	lambda = bandData[0,i]*microns
	dlambda = bandData[1,i]*microns
	dlambda0 = 0.09*microns
	nu = c/lambda
	F0 = 3.209e-9 ; W/m2 Vega in V band
	F_lambda_ratio = [0.07081, $ ; B
	                  0.03602, $ ; V
	                  0.02201, $ ; R
	                  0.01177, $ ; I
	                  0.00318, $ ; J
	                  0.00118, $ ; H
	                  0.000417] ; Ks
	F_lambda_ratio /= 0.0364 ; @ 5480 A
	F = F0*10^(-m/2.5)*F_lambda_ratio[i]*dlambda/dlambda0
	nPh = F/(h*nu) ; ph/s/m2
	return,nPh
end

; -- estimate the Plank star temperature given pass band data, degrees K --
function photoColor,band1,m1,band2,m2
	forward_function fluxRatio
	; - find the relative flux in each band -
	Flux_ratio = photoCounts(band1,m1)/photoCounts(band2,m2)
	dT = 1000
	nT = 15
	T = (findgen(nT)+1)*dT
	fr = fluxRatio(T,band1,band2) ; this generates a table of Plank flux ratios vs Temperature
	Tstar = interpol(T,fr,[Flux_ratio])
	return,Tstar
end

; -- generate Plank energy or photon flux ratio vs Temperature: flux(band1)/flux(band2) --
function fluxRatio,T,band1,band2,photons=photons
	c = 3.0e8 ; m/s
	h = 6.63e-34 ; Js
	k = 1.38e-23 ; J/K
	microns = 1.e-6
	calcBandData,bandName,bandData
	k1 = (where(band1 eq bandName))[0]
	k2 = (where(band2 eq bandName))[0]
	lambda1 = bandData[0,k1]*microns
	dlambda1 = bandData[1,k1]*microns
	lambda2 = bandData[0,k2]*microns
	dlambda2 = bandData[1,k2]*microns
	fluxRatio = (dlambda1/dlambda2)*(lambda2/lambda1)^5*(exp(h*c/(lambda2*k*T))-1)/(exp(h*c/(lambda1*k*T))-1)
	if (keyword_set(photons)) then fluxRatio *= (lambda1/lambda2)
	return,fluxRatio
end

; -- examples --
band = 'V'
mV = 0
print,photoCounts(band,mV)
band1 = 'B'
band2 = 'V'
n = 100
bv_max = 1.6
bv_min = 0.2
dbv = (bv_max-bv_min)/float(n-1)
m_band1_set = bv_min+findgen(n)*dbv
Tset = fltarr(n)
for k = 0,n-1 do begin
	m_band1 = m_band1_set[k]
	m_band2 = 0
	Tset[k] = photoColor(band1,m_band1,band2,m_band2)
endfor
plot,m_band1_set,Tset,xrange=[0,1.6],yrange=[3500,7000]
end