function kepler_mag,gmag,rmag,imag
; This routine returns an estimate of the Kepler magnitude of a star,
; given any two or all three of the SDSS g,r,i magnitudes.
; Invalid input magnitudes are indicated by any value larger than 40.
; The routine decides which formula to use based on which inputs are
; available.
; Formulas are accurate to about +/- 0.2 mag (relative to simulations)
; for stars hotter than 3500 K
; except if only gmag and rmag are valid, in which case errors for M stars
; may exceed 0.6 mag, in the sense that kepler_mag returns magnitudes that
; are too faint.
; The returned magnitudes have small discontinuous jumps depending
; on the input colors, so output should not be differentiated.

nst=n_elements(gmag)
gok=intarr(nst)
rok=intarr(nst)
iok=intarr(nst)
kp=fltarr(nst)

s=where(gmag gt 0. and gmag lt 40.,ns)
if(ns gt 0) then gok(s)=1
s=where(rmag gt 0. and rmag lt 40.,ns)
if(ns gt 0) then rok(s)=1
s=where(imag gt 0. and imag lt 40.,ns)
if(ns gt 0) then iok(s)=1

; check for insufficient data
sbad=where(gok+rok+iok eq 0,nsbad)
if(nsbad gt 0) then kp(sbad)=99.9
if(nsbad eq nst) then goto,fini

s1g=where(gok eq 1 and rok eq 0 and iok eq 0,ns1g)
if(ns1g gt 0) then kp(s1g)=gmag
s1r=where(gok eq 0 and rok eq 1 and iok eq 0,ns1r)
if(ns1r gt 0) then kp(s1r)=rmag
s1i=where(gok eq 0 and rok eq 0 and iok eq 1,ns1i)
if(ns1i gt 0) then kp(s1i)=imag

; all three inputs are okay
s3=where(gok eq 1 and rok eq 1 and iok eq 1,ns3)
if(ns3 gt 0) then begin
  gmr=gmag(s3)-rmag(s3)
  kp3=fltarr(ns3)
  gmag3=gmag(s3)
  rmag3=rmag(s3)
  imag3=imag(s3)
  sred=where(gmr gt 0.3,nsred) 
  if(nsred gt 0) then kp3(sred)=0.7*imag3(sred)+0.3*gmag3(sred) 
  sblue=where(gmr le 0.3,nsblue)
  if(nsblue gt 0) then kp3(sblue)=.75*rmag3(sblue)+.25*gmag3(sblue)
  kp(s3)=kp3
endif

; r and i are okay
sri=where(gok eq 0 and rok eq 1 and iok eq 1,nsri)
if(nsri gt 0) then begin
  rmi=rmag(sri)-imag(sri)
  kpri=fltarr(nsri)
  rri=rmag(sri)
  iri=imag(sri)
  sred=where(rmi gt .673,nsred) 
  sblue=where(rmi le .673,nsblue)
  if(nsred gt 0) then kpri(sred)=1.2*rri(sred)-0.2*iri(sred)
  if(nsblue gt 0) then kpri(sblue)=.65*rri(sblue)+.35*iri(sblue)
  kp(sri)=kpri
endif

; g and r are okay
sgr=where(gok eq 1 and rok eq 1 and iok eq 0,nsgr)
if(nsgr gt 0) then begin
  gmr=gmag(sgr)-rmag(sgr)
  kpgr=fltarr(nsgr)
  rgr=rmag(sgr)
  ggr=gmag(sgr)
  sred=where(gmr gt 0.8,nsred) 
  sblue=where(gmr le 0.8,nsblue)
  if(nsred gt 0) then kpgr(sred)=0.9*rgr(sred)+.1*ggr(sred)
  if(nsblue gt 0) then kpgr(sblue)=0.8*rgr(sblue)+0.2*ggr(sblue)
  kp(sgr)=kpgr
endif

; g and i are okay
sgi=where(gok eq 1 and rok eq 0 and iok eq 1,nsgi)
if(nsgi gt 0) then begin
  gmi=gmag(sgi)-imag(sgi)
  kpgi=fltarr(nsgi)
  igi=imag(sgi)
  ggi=gmag(sgi)
  sred=where(gmi gt (-.05),nsred) 
  sblue=where(gmi le (-.05),nsblue)
  if(nsred gt 0) then kpgi(sred)=0.7*igi(sred)+.3*ggi(sred)
  if(nsblue gt 0) then kpgi(sblue)=0.45*igi(sblue)+0.55*ggi(sblue)
  kp(sgi)=kpgi
endif

fini:
return,kp

end