pro cphot,fname,filtlist,teff,logg,logz,colors,names,ext,rv,erng,rvrng,ab=ab,$
atmos=atmos,castelli=castelli
; This routine reads the BaSeL flux distribuition file fname, and normalizes
; the fluxes so that the total from 500 to 550 nm is unity.
; If argument erng is present, the routine also constructs reddened models,
; with total extinction (magnitudes) given by the values in array erng.
; If erng is absent, extinction = 0.
; If argument rvrng is present, for each extinction value, fluxes are
; also computed for the range of R_v values given in array rvrng.
; If rvrng is absent, R_v defaults to 3.1 (typical for diffuse ISM).
; If erng is absent, rvrng must be absent also.
; It reads the ascii file filtlist, which contains the lower and upper
; wavelengths of a list of nfilt flat-topped transmission profiles.
; It integrates the normalized fluxes over each of these bandpasses,
; and converts the resulting fluxes to instrumental magnitudes (ie, -2.5*log).
; It also reads from file filtlist a list of color indices (filter i - filter j)
; It creates a table color(nindex,nmodel) of color indices for each model
; and returns these, along with string names for the colors.
; It also returns the model parameters teff(nmodel), logg(nmodel), logz(nmodel).
; Format for input file filtlist:
; If keyword ab is set, the colors returned are from AB magnitudes, not
; from those normalized to Vega.
; If keyword atmos is a nonzero value and keyword ab is set, then the
; AB magnitudes are computed at airmass=atmos, rather than airmass=0.
; nfilt = number of filters (i4)
; If keyword castelli is set, input flux data are read from a file containing
; castelli (not basel) models.
; for each filter:
; filtname = string, eg Z or g'
; lam_bot, lam_top = min, max wavelengths in bandpass (nm) (2 floating values)
; ftype, fwid = type of filter profile, width
; if ftype=0, boxcar filter spanning lambot:lamtop
; if ftype=1, gaussian centered at (lambot+lamtop)/2, with
; fwhm = fwid nm. Sampled only in range lambot:lamtop
; if ftype=2, read filter profile from file ftname (which is read from
; next line), and interpolate this onto range lambot:lamtop
; In this case, fwid is ignored.
; ftnam = name of ascii file containing filter profile. Read and ignored
; unless ftype=2
; Then for the color indices,
; nindex = number of color indices to compute (i4)
; for each index:
; i, j where color = mag_i - mag_j, 0 <= i,j <= nindex-1
;
; for example, this defines UBV filters, colors U-B and B-V:
; 3
; U
; 350. 400.
; 0 50.
; dummy string for a boxcar filter
; B
; 400. 500.
; 1 40.
; dummy line for a gaussian filter
; V
; 500. 620.
; 2 0.
; /name/of/V/filter/ascii/file
; 2
; 0 1
; 1 2
; constants
lammin=500.
lammax=550.
extinct_file='~/mosaic/extinct.asc'
; count the reddening cases to be considered
if(n_params() lt 10) then erng=[0.]
if(n_params() lt 11) then rvrng=[3.1]
nex=n_elements(erng)
nrv=n_elements(rvrng)
; read the flux file, normalize fluxes
if(keyword_set(castelli)) then begin
rd_castelli,fname,lambda,nmod,teff,logg,logz,fluxi
endif else begin
rd_basel,fname,lambda,nmod,mnum,teff,logg,logz,fluxi
endelse
nmodel=n_elements(teff)
nlambda=n_elements(lambda)
s=where(lambda ge lammin and lambda le lammax,ns)
for i=0,nmodel-1 do begin
fluxi(*,i)=fluxi(*,i)/total(fluxi(s,i))
endfor
; make the reddened fluxes
flux=fltarr(nlambda,nmodel,nex,nrv)
for i=0,nmodel-1 do begin
for j=0,nex-1 do begin
for k=0,nrv-1 do begin
flux(*,i,j,k)=redden(lambda,fluxi(*,i),erng(j),rvrng(k))
endfor
endfor
endfor
nmodtot=nmodel*nex*nrv
flux=reform(flux,nlambda,nmodtot)
; if keywords ab and atmos are set, further filter the fluxes with
; terrestrial atmospheric extinction.
if(keyword_set(ab) and keyword_set(atmos)) then begin
openr,iun,extinct_file,/get_lun
ss=''
readf,iun,ss
readf,iun,ss
lame=fltarr(35)
exte=fltarr(35)
for i=0,34 do begin
readf,iun,v1,v2
lame(i)=v1
exte(i)=v2
endfor
close,iun
free_lun,iun
exte=exte*atmos
extei=interpol(exte,lame,lambda)
extei=(extei < 30.) > 0.
atmtran=10.^(-0.4*extei)
atmtran=rebin(atmtran,nlambda,nmodel,nex,nrv)
flux=flux*atmtran
endif
; make all the other parameter lists concordant
teff=rebin(teff,nmodel,nex*nrv)
teff=reform(teff,nmodtot)
logg=rebin(logg,nmodel,nex*nrv)
logg=reform(logg,nmodtot)
logz=rebin(logz,nmodel,nex*nrv)
logz=reform(logz,nmodtot)
ext=reform(erng,1,nex,1)
ext=rebin(ext,nmodel,nex,nrv)
ext=reform(ext,nmodtot)
rv=reform(rvrng,1,1,nrv)
rv=rebin(rv,nmodel,nex,nrv)
rv=reform(rv,nmodtot)
; read the filtlist
openr,iun,filtlist,/get_lun
readf,iun,nfilt
filnames=strarr(nfilt)
lbot=fltarr(nfilt)
ltop=fltarr(nfilt)
ftype=intarr(nfilt)
fwid=fltarr(nfilt)
ftnam=strarr(nfilt)
ss=''
for i=0,nfilt-1 do begin
readf,iun,ss
filnames(i)=strtrim(ss,2)
readf,iun,wb,wt
lbot(i)=wb
ltop(i)=wt
readf,iun,it,wid
ftype(i)=it
fwid(i)=wid
readf,iun,ss
ftnam(i)=strtrim(ss,2)
endfor
readf,iun,nindex
names=strarr(nindex)
ija=lonarr(2,nindex)
for i=0,nindex-1 do begin
readf,iun,ix,jx
ija(0,i)=ix
ija(1,i)=jx
names(i)=filnames(ix)+'-'+filnames(jx)
endfor
; make zero-point offsets from width values, for ftype=2.
zpoff=fwid
sz=where(ftype ne 2,nsz)
if(nsz gt 0) then zpoff(sz)=0.
close,iun
free_lun,iun
; make magnitudes in the various filters
; photon flux = energy flux * nu/(c*h)
; dnu = -c*lambda^-2*dlambda
; so fe*dnu/nu = fp*dlambda
mags=fltarr(nfilt,nmodtot)
for i=0,nfilt-1 do begin
sl=where(lambda ge lbot(i) and lambda le ltop(i),nsl)
if(nsl gt 0) then begin
ccd_filtprofile,lambda(sl),ftype(i),fwid(i),ftnam(i),prof
dlam=deriv(lambda)
dlam=dlam(sl)
lams=lambda(sl)
for j=0,nmodtot-1 do begin
if(keyword_set(ab)) then begin
num=total(flux(sl,j)*prof*dlam/lams) > 1.e-8
den=total(prof*dlam/lams)
tf=num/den > 1.e-8
endif else begin
tf=total(flux(sl,j)*prof*dlam) > 1.e-6
endelse
mags(i,j)=-2.5*alog10(tf) + zpoff(i)
endfor
endif
endfor
; make the colors
colors=fltarr(nindex,nmodtot)
for i=0,nindex-1 do begin
for j=0,nmodtot-1 do begin
colors(i,j)=mags(ija(0,i),j)-mags(ija(1,i),j)
endfor
endfor
;plotfigs,lambda,flux,teff,logg,logz,ext,colors
;stop
end