; Given type (string, e.g, "O"), subtype (float, e.g. 3)
; and class (integer), return H-K from
; AQ 4th ed
;

function aq_hk, type, subtype, class
  functionname = 'aq_hk'
  CASE class OF
    5: begin ; Main sequence
	    CASE type OF
		'O': begin 
			x = [9,9.5] & y = [-0.04,-0.04]
		     end
		'B': begin 
			x = [0,1,2,3,4,5,6,7,8,9] 
			y = [-0.04,-0.03,-0.03,-0.02,-0.02,-0.01,-0.01,-0.01,0.00,0.00]
		     end
		'A': begin 
			x = [0,2,5,7] & y = [0.00,0.01,0.02,0.03]
		     end
		'F': begin 
			x = [0,2,5,7] & y = [0.03,0.04,0.04,0.05]
		     end
		'G': begin 
			x = [0,2,4,6] & y = [0.05,0.05,0.06,0.06]
		     end
		'K': begin 
			x = [0,2,4,5,7]
			y = [0.08,0.09,0.11,0.11,0.15]
		     end
		'M': begin 
			x = [0,1,2,3,4,5,6] 
			y = [0.17,0.18,0.20,0.23,0.27,0.33,0.38]
		     end
	  ENDCASE
       end
    3: begin ; giants. Teff for O-F are from xx
	    CASE type OF
		'G': begin 
			x = [0,4,6,8] & y = [0.07,0.08,0.09,0.09]
		     end
		'K': begin 
			x = [0,1,2,3,4,5] & y = [0.10,0.10,0.12,0.14,0.15,0.17]
		     end
		'M': begin 
			x = [0,1,2,3,4,5,6,7]
			y = [0.19,0.21,0.22,0.24,0.25,0.29,0.30,0.31]
		     end
		 ELSE: begin ; out of bounds! 
		    print, functionname+': spectral type out of bounds'
		    return, -100
		    end
	  ENDCASE
       end
    1: begin ; super giants. 
	    CASE type OF
		'B': begin 
 			   x = [5,6,7,8,9] 
			   y = [0.00,-0.02,-0.02,-0.02,-0.02]
		     end
		'A': begin 
			x = [0,1,2,5] & y = [-0.02,-0.01,-0.01,0.02]
		     end
		'F': begin 
			x = [0,2,5,8] & y = [0.04,0.05,0.06,0.07]
		     end
		'G': begin 
			x = [0,3,8] & y = [0.08,0.09,0.11]
		     end
		'K': begin 
			x = [0,1,2,3,5] & y = [0.12,0.13,0.13,0.13,0.14]
		     end
		'M': begin 
			x = [0,1,2,3,4] & y = [0.18,0.20,0.22,0.24,0.26]
		     end
		  ELSE: begin ; out of bounds! 
		    print, functionname+': spectral type out of bounds'
		    return, -100
		    end
	  ENDCASE
       end
  ELSE: begin ; out of bounds! 
    print, functionname+': spectral type out of bounds'
    return, -100
    end
  ENDCASE

;  res = poly_fit(x, y, 1) 
;  color = poly(subtype,res)
  color = interpol(y,x,subtype)
  return, color
end