pro phot_models,inlist,redprm,filtlist,savo,teff,logg,logz,ext,rv,names,$
    colors,ab=ab,atmos=atmos,castelli=castelli
; This routine computes a set of photometric indices (colors) for all of
; the BaSeL models in the files listed in the ascii file inlist.  Moreover,
; for each BaSeL model it computes a grid of reddened model spectra,
; with 2-parameter reddening values givin in input ascii file redprm.
; The full list of BaSeL models runs to about 8000 entries, so the total
; table of colors will have as many as ntot=8000*nred*nrv entries.
; Output consists of ntot-long vectors teff, logg, logz, red, rv, which give 
; the model parameters, string array names(nc), 
; which contains the names of the nc filter differences (colors) for which 
; values are computed, and array colors(nc,ntot), which contains the color
; information.
; Colors are computed as differences between photon fluxes integrated over
; the specified filter bands, converted to magnitudes, with zero points
; determined by specifying that for the spectrum lcb97_vega.ori, all colors
; are zero.
; Before returning, all output data are written as an IDL save file to savo.
; Reddening is defined by ext = A(V) = total extinction in the V band, and
; rv = A(V)/E(B-V), as in Cardelli et al. ApJ 345, 245 (1989).
; If keyword ab is set, then colors are computed from AB magnitudes,
; (ie, energy fluxes per unit frequency, not photon fluxes per unit lambda)
; and are not normalized to Vega.
; If keywords ab and atmos are set, then fluxes are computed as seen with
; airmass=atmos from Mt. Palomar, a la Hayes & Latham.
; If keyword castelli is set, the routines assume castelli models.

; constants
veganame='spectra/Original/lcb97_vega.ori'

; construct colors for vega, get number of colors
cphot,veganame,filtlist,vt,tg,vz,vega_c,names,vext,vrv,ab=ab,atmos=atmos
  
ncols=n_elements(names)

; read input file list, count them, make an array of their names
openr,iun,inlist,/get_lun
nfiles=0
ss=''
while(not eof(iun)) do begin
  readf,iun,ss
  ss=strtrim(ss,2)
  if(nfiles eq 0) then fnames=[ss] else fnames=[fnames,ss]
  nfiles=nfiles+1
endwhile
close,iun
free_lun,iun

; read redprm file, count number of extinction values, make oversized
; colors, parms arrays.  Also make array to count total models as they come in
openr,iun,redprm,/get_lun
readf,iun,nex
if(nex gt 0) then erng=fltarr(nex) else erng=[0.]
if(nex gt 0) then readf,iun,erng
readf,iun,nrv
if(nrv gt 0) then rvrng=fltarr(nrv) else rvrng=[3.1]
if(nrv gt 0) then readf,iun,rvrng
nred=nex*nrv > 1
colo=fltarr(ncols,500*nred,nfiles)        ; assumes <= 500 models per file
ngood=lonarr(nfiles)
teffo=fltarr(500*nred,nfiles)
loggo=fltarr(500*nred,nfiles)
logzo=fltarr(500*nred,nfiles)
exto=fltarr(500*nred,nfiles)
rvo=fltarr(500*nred,nfiles)

; loop over models, accumulating colors, etc.
for i=0,nfiles-1 do begin
  fname=fnames(i)
  cphot,fname,filtlist,teff,logg,logz,colors,names,ext,rv,erng,rvrng,ab=ab,$
    atmos=atmos,/castelli
  nmodtot=n_elements(teff)
  nm=nmodtot-1
  ngood(i)=nmodtot
  vega_t=rebin(vega_c,ncols,nmodtot,1)
  if(not keyword_set(ab)) then begin
    colo(*,0:nm,i)=colors-vega_t           ; normalization to Vega done here
  endif else begin
    colo(*,0:nm,i)=colors
  endelse
  teffo(0:nm,i)=teff
  loggo(0:nm,i)=logg
  logzo(0:nm,i)=logz
  exto(0:nm,i)=ext
  rvo(0:nm,i)=rv
endfor

; make output colors, etc. arrays, straining out all of the empty cells
nn=total(ngood)
colors=fltarr(ncols,nn)
teff=fltarr(nn)
logg=fltarr(nn)
logz=fltarr(nn)
ext=fltarr(nn)
rv=fltarr(nn)
ncur=0L
for i=0,nfiles-1 do begin
  ntop=ncur+ngood(i)-1L
  ng=ngood(i)-1L
  colors(*,ncur:ntop)=colo(*,0:ng,i)
  teff(ncur:ntop)=teffo(0:ng,i)
  logg(ncur:ntop)=loggo(0:ng,i)
  logz(ncur:ntop)=logzo(0:ng,i)
  ext(ncur:ntop)=exto(0:ng,i)
  rv(ncur:ntop)=rvo(0:ng,i)
  ncur=ncur+ng+1
endfor

; save the results
save,teff,logg,logz,ext,rv,names,colors,file=savo

end