pro lext, x, flux, fluxerr, z=z, extmodel=extmodel, fixav=fixav, fixrvinput=fixrvinput, $
    fixbeta=fixbeta, fixbump=fixbump, fixuv=fixuv, fixc1=fixc1, fixc2=fixc2, fixc3=fixc3, fixc4=fixc4, fixf0=fixf0, wv0=wv0, $
    fitredshift=fitredshift, galav=galav, rangebeta=rangebeta, title=title, labels=labels, nolegend=nolegend,$
    fitfilters=fitfilters, unextcurve=unextcurve, ext0fit=ext0fit, time=time, comment=comment, $
    morecomment=morecomment, xr=xr, yr=yr, centfilt=centfilt, shortcaption=shortcaption, plotfile=plotfile,$
    residuals=residuals, ext2model=ext2model, font=font, multi=multi, linestyle=linestyle, line2style=line2style, nofit=nofit, $
    hnudge=hnudge, hideigm=hideigm, wav=wav, freq=freq, bbplot=bbplot

; Wrapper for fitextinct, produces extinction plots fit to a set of data.
; Unfortunately you have to run it twice to get the fonts right.


if n_elements(multi) eq 2 then begin
  multi = [multi, 0, 0]
endif

if n_elements(nofit) eq 0 then begin
  nofit = 0
endif else begin
  nofit = 1
  extmodel = 'none'
endelse

if n_elements(ext2model) and (keyword_set(ext0fit) or keyword_set(unextcurve)) then shortcaption=1
if n_elements(font) eq 0 then font = 'times'
if n_elements(hnudge) eq 0 then hnudge = 0
if n_elements(unextcurve) eq 0 then unextcurve = 0

if n_elements(fixrvinput) gt 0 then begin
  if fixrvinput eq 0 then begin
     rvfixed = 0      ;fixrv = 0 FREES R_V
  endif else begin
     rvfixed = 1   
     fixrv = fixrvinput
  endelse
endif else begin
  if extmodel eq 'smc' then fixrv = 2.73
  if extmodel eq 'lmc' then fixrv = 3.2
  if extmodel eq 'mw'  then fixrv = 3.1
  if extmodel eq 'calz' then fixrv = 4.0
  rvfixed = 1
endelse
if keyword_set(fitredshift) then begin
  zfixed = 0         
endif else begin
  zfixed = 1
  fixz = z
endelse
if n_elements(fixbump) gt 0 then fixc3 = fixbump * 3.23 ; check this value
if n_elements(fixuv) gt 0 then fixc4 = fixuv * 0.41

betafixed = n_elements(fixbeta) gt 0
avfixed = n_elements(fixav) gt 0
f0fixed = n_elements(fixf0) gt 0
c1fixed = n_elements(fixc1) gt 0
c2fixed = n_elements(fixc2) gt 0
c3fixed = n_elements(fixc3) gt 0
c4fixed = n_elements(fixc4) gt 0
bumpfixed = c3fixed
uvfixed = c4fixed

if n_elements(ext0fit) eq 0 then ext0fit = 0
if n_elements(unextfit) eq 0 then unextfit = 0

storefluxerr = fluxerr
storex = x
if n_elements(fitfilters) gt 0 then storefitfilters = fitfilters

; filter/wavelength translation and Galactic extinction correction

if n_elements(wav) eq 0 and n_elements(freq) eq 0 then begin
   filters = x
   wveffobs = filtwv(filters, fwhm=wverrobs)
endif
if n_elements(wav) eq 1 then begin
   filters = strarr(n_elements(x))
   wveffobs = x
   wverrobs = x*0
endif
if n_elements(freq) eq 1 then begin
   filters = strarr(n_elements(x))
   wveffobs = 3e+18/x
   wverrobs = x*0
endif


fluxcor = flux
if n_elements(galav) gt 0 then begin  
  notxray = where(wveffobs gt 1000)
  extmag = extinctionfitz(wveffobs, 3.1, galav)
  bad = where(finite(extmag) eq 0, ct)
    if ct gt 0 then extmag[bad] = 0
  extfactor = 10.0^(extmag/2.5)
  fluxcor[notxray] = fluxcor[notxray] * extfactor
  fluxerr[notxray] = fluxerr[notxray] * extfactor
endif else begin
  print, 'Assuming Galactic A_V = 0, or already corrected.'
endelse
freqeff = 2.99792e18 / wveffobs
freqerr = 2.99721e18 * wverrobs / (wveffobs^2)

wverrobs = wverrobs / 2.0   ;2.35


; find the filters

if n_elements(fitfilters) eq 0 then fitfilters = indgen(n_elements(filters))
f = fitfilters
xx = 1 + intarr(n_elements(filters))
xx[f] = 0
nf = where(xx eq 1, nnf) ; count the number of usable filters in a weird way

if nofit eq 0 then begin
  fitextinct, wveffobs[f], fluxcor[f], fluxerr[f], exttype=extmodel, $
  par=sol, perr=parerr, pname=parname, fluxmodel=fluxmodel, $
  curvewav=curvewav, curveflux=curveflux, chisq=chisq, dof=dof, $ 
  rangebeta=rangebeta, fixf0=fixf0, wv0=wv0, fixbeta=fixbeta, fixav=fixav, fixrv=fixrv, fixz=fixz, $
  fixc1=fixc1, fixc2=fixc2, fixc3=fixc3, fixc4=fixc4


  storewveffobs = wveffobs[f]   ; I really should redefine to preserve these.
  storefluxmodel = fluxmodel


  betapar = (where(parname eq 'beta')) [0]
  beta = sol[betapar]
  betaerr = parerr[betapar]
  avpar = (where(parname eq 'a')) [0]
  av = sol[avpar]
  averr = parerr[avpar]
  rvpar = (where(parname eq 'r')) [0]
  rv = sol[rvpar]
  rverr = parerr[rvpar]
  bumppar = (where(parname eq 'c3')) [0]
  bump = sol[bumppar] 
  bumperr = parerr[bumppar]
  uvpar = (where(parname eq 'c4')) [0]
  uv = sol[uvpar] 
  uverr = parerr[uvpar]
  c1par = (where(parname eq 'c1')) [0]
  c1 = sol[c1par] 
  c1err = parerr[c1par]
  c2par = (where(parname eq 'c2')) [0]
  c2 = sol[c2par] 
  c2err = parerr[c2par]
  c3par = (where(parname eq 'c3')) [0]
  c3 = sol[c3par] 
  c3err = parerr[c3par]
  c4par = (where(parname eq 'c4')) [0]
  c4 = sol[c4par] 
  c4err = parerr[c4par]

  if extmodel eq 'fm1' then begin
      c1    =  2.030 - 3.007*c2
      c1err =  abs(2.030 - 3.007*(c2+c2err) - c1)
  endif
  if extmodel eq 'fm2' then begin
      c2    = -0.824 + 4.717/rv
      c2err = abs(-0.824 + 4.717/(rv+rverr) - c2)
      c1    =  2.030 - 3.007*c2
      c1err =  abs(2.030 - 3.007*(c2+c2err) - c1)
  endif
  if extmodel eq 'reichart' then begin
      c1 = -0.064 - 3.275*(c2 - 0.711)
      c1err = abs(-0.064 - 3.275*(c2+c2err - 0.711)  - c1)
      rv =  3.228 - 2.685*(c2 - 0.721) + 1.806*(c2 - 0.721)^2
      rverr = abs(3.228 - 2.685*(c2+c2err - 0.721) + 1.806*(c2+c2err - 0.721)^2 - rv)
  endif


  ; Intrinsic curve - best fit with extinction zeroed

  unextsol = sol
  unextsol[avpar] = 0.0     ;sol, but without extinction (A_V = 0)
  unextcurveflux = extpw(curvewav, unextsol, noigm=hideigm)
  unextmodelflux = extpw(wveffobs[f], unextsol)
  afilt = 2.5*alog10(unextmodelflux/fluxcor[f])

  openw, 2, 'extinfo.txt'
  printf, 2, extmodel
  printf, 2, 'beta=',+beta
  printf, 2, 'chi^2/dof = ',clip(chisq),'/',clip(dof)
  printf, 2, clip('filt',6), ' ', clip('wav_eff',8), ' ', clip('flux',11), ' ', clip('unc',10), ' ', clip('model',10), '', clip('intrins',10), ' ', clip('A_filt',10)
  for i = 0, n_elements(f)-1 do begin
    printf, 2, clip(filters[f[i]],6), ' ', fpr(wveffobs[f[i]],5.2), ' ', clip(fluxcor[f[i]],11), ' ', clip(fluxerr[f[i]],10), ' ', clip(fluxmodel[i],10), clip(unextmodelflux[i],10), ' ', clip(afilt[i],10)
  endfor
  close, 2

endif

; Secondary fit

if n_elements(ext2model) gt 0 then begin

  ; The behavior of the secondary fit is necessarily a bit weird because the lack of options to control it.

  if ext2model eq 'smc' then fixrv2 = 2.73
  if ext2model eq 'lmc' then fixrv2 = 3.2
  if ext2model eq 'mw'  then fixrv2 = 3.1
  if ext2model eq 'calz' then fixrv2 = 4.0

  fitextinct, wveffobs[f], fluxcor[f], fluxerr[f], exttype=ext2model, $
  par=sol2, perr=parerr2, fluxmodel=fluxmodel2, $
  curvewav=curvewav, curveflux=curveflux2, chisq=chisq2, dof=dof2, $ 
  rangebeta=rangebeta, fixf0=fixf0, fixbeta=fixbeta, fixz=fixz, fixrv=fixrv2
  ;fixf0=fixf0, fixbeta=fixbeta, fixav=fixav, fixrv2=fixrv2, fixz=fixz, fixc1=fixc1, fixc2=fixc2, fixc3=fixc3, $
  ;fixc4=fixc4
endif




if ext0fit ge 1 then begin
  ; No-extinction solution with no beta constraint
  fitextinct, wveffobs[f], fluxcor[f], fluxerr[f], exttype='none', fixz=fixz, $
  par=fitext0sol, perr=fitext0perr, pname=fitext0pname, fluxmodel=fitext0fluxmodel, $
  curvewav=curvewav, curveflux=fitext0curveflux, chisq=fitext0chisq, dof=fitext0dof
endif


xray = where(wveffobs lt 1000, ctxray)
if ctxray gt 0 then begin
   wveffobsx = wveffobs[xray]
   wverrobsx = wverrobs[xray]   ; preserve these
   fluxcorx = fluxcor[xray]
   fluxerrx = fluxerr[xray]
endif

notxray = where(wveffobs gt 1000)
if n_elements(notxray) ne n_elements(wveffobs) then begin
  filters = filters[notxray]
  fitfilters = fitfilters and notxray
  fluxcor = fluxcor[notxray]
  fluxerr = fluxerr[notxray]
  wveffobs = wveffobs[notxray]
  wverrobs = wverrobs[notxray]
endif



plotresiduals = keyword_set(residuals)
if n_elements(plotfile) eq 0 then plotfile = 'out.ps'

; open plot
if n_elements(multi) eq 0 then begin
  xs = 5.0
  ys = 3.8+1.5*plotresiduals
  initdev = !d.name
  set_plot, 'ps'
  if strpos(plotfile,'.eps') gt -1 then device, /encapsulated
  device, filename=plotfile
  !p.font = 0
  device, /symbol, font_index = 4
  if font eq 'helvetica' then begin
    device, /helvetica, font_index = 17
    device, /helvetica, /bold, font_index = 18
  endif else begin
    device, /times, font_index = 17
    device, /times, /bold, font_index = 18
  endelse
  device, /isolatin1, font_index = 19
  device, /inches, xsize=xs, ysize=ys
  device, /inches, yoffset=4.5+(5.5-ys)
  color = 1
  if (color) then begin
    device, /color
    loadfiltcolors
  endif
endif else begin
  !p.multi = [2,2,1]  ; never repage since !p.multi is crap
endelse

;minwv = (min(wveffobs[fitfilters])*0.85/1000)*1000
minwv = 912 * (1+z)
maxwv = (max(wveffobs[fitfilters])/0.85/5000)*5000
minflux = (min(fluxcor[fitfilters])*0.8)
maxflux = (max(fluxcor[fitfilters])/0.8)




xtv = [1000, 1500, 2000, 3000, 4000, 5000, 7500, 10000, 15000, 20000, 25000]
xtv = xtv[where(xtv ge minwv and xtv le maxwv)]

if n_elements(xr) eq 0 then xr = [max(wveffobs)*1.2, min(wveffobs)*0.8]


abmags = ujy2mag(fluxcor, 'AB')
abmagrange = ujy2mag([minflux,maxflux], 'AB')
abmagerr = abmags - ujy2mag(fluxcor+fluxerr,'AB')

ycent = mean(abmagrange)
if n_elements(centfilt) gt 0 then begin
  centindex = where(filters eq centfilt, ct) 
  if ct eq 1 then ycent = (ujy2mag(fluxcor[centindex],'AB')) [0]
endif

if n_elements(yr) eq 0 then yrange = abmagrange + [0.5, -0.5]
if n_elements(yr) eq 1 then yrange = ycent + [yr/2.,-yr/2.]
if n_elements(yr) eq 2 then yrange = yr
if yrange[0] lt yrange[1] then yrange = [yrange[1],yrange[0]]
minflux = mag2ujy(yrange[0],'AB')
maxflux = mag2ujy(yrange[1],'AB')

if maxflux/minflux lt 100 then ytv = [1, 2, 5] else ytv = [1]

ytv = [ytv*0.01, ytv*0.1, ytv, ytv*10, ytv*100, ytv*1000, ytv*1e4, ytv*1e5, ytv*1e6, ytv*1e7]
ytv = ytv[where(ytv ge minflux and ytv le maxflux)]

xtitle = '!4l!17!Deff!N (!19'+string(197B) + '!17)' ; !6!sA!r!u!9 %!6!n!4!17
x2title ='    !4l!17!Deff,rest!N (!19'+string(197B) + '!17)' ; !6!sA!r!u!9 %!6!n!4!17
ytitle = '!17magnitude (AB)'
y2title = '!17F!D!4n!17!N (!4m!17Jy)'

cols = getfiltcolors(filters)

if n_elements(multi) gt 0 then begin
  xtf =    (multi[3] eq multi[1]-1)       ? '(I)' : '(A1)'
  xtitle = (multi[3] eq multi[1]-1)       ? xtitle : ''
  ytitle = (multi[2] eq 0)                ? ytitle : ''
endif else begin
  xtf = '(I)'
endelse


xnorm_ch_size = 1.0*!d.x_ch_size / !d.x_vsize
ynorm_ch_size = 1.0*!d.y_ch_size / !d.y_vsize

left = 7. * xnorm_ch_size + hnudge
bot =  3.2 * ynorm_ch_size
right = 1 - 7. * xnorm_ch_size + hnudge
if n_elements(multi) gt 0 then begin
   if multi[0] gt 2 then right = 1 - 1.5 * xnorm_ch_size 
endif
top =   1 - 3. * ynorm_ch_size
if n_elements(multi) gt 0 then begin
  yinc =  (top-bot)/multi[1]
  xinc =  (right-left)/multi[0]
endif

if n_elements(multi) eq 0 then begin
  !p.position = [left, bot, right, top]
endif else begin
  !p.position = [left+xinc*(multi[2]), top-yinc*(multi[3]+1), left+xinc*(multi[2]+1), top-yinc*(multi[3])]
endelse

spacing = 1.4 * ynorm_ch_size * ((!p.position[3]-!p.position[1])/0.8)^0.2
xmargin = 3.  * xnorm_ch_size * (!p.position[2]-!p.position[0])/0.8
ymargin = 1.  * ynorm_ch_size * (!p.position[3]-!p.position[1])/0.8




if plotresiduals then begin
  !p.multi = [0,1,2] ; safe generally
  !p.position = [0.12, 0.369, 0.89, 0.93]
  respos = [0.12, 0.086, 0.89, 0.369]
endif




  ; set up plot frame

  if plotresiduals then xtf1 = '(A1)' else xtf1 = xtf
  if plotresiduals then xtitle1 = '' else xtitle1 = xtitle

  ytitle1 = ytitle
  if n_elements(yrange) eq 2 then begin
     if max(yrange)-min(yrange) gt 3.5 then ytf = '(I)'
  endif
  if (minflux gt 1 and maxflux lt 100000.) then rytf = '(I)'

  if n_elements(multi) gt 0 then begin
     hugewhitespace = '                                                       '

     ystyle = 1
     ; left axis
     if multi[2] gt 0 then ytitle1=''                          ;no title unless on left side
     if multi[1] eq 3 and multi[3] ne 1 then ytitle1=''        ;if 3 vertical plots, title only on middle plot
     if multi[0] eq 2 and multi[2] eq 1 then ytf='(A1)'        ;if 2 horizontal plots, no tick labels on right plots 

     ;bottom axis
     if multi[0] eq 2 and multi[2] gt 0 then xtitle1=''                              ;if 2 horizontal plots, only one title
     if multi[0] eq 2 and multi[2] eq 0 then xtitle1=hugewhitespace+'      '+xtitle1 ;space out that title to put it in the center

     ;top axis
     xtf     = (multi[3] eq 0) ? '(I)' : '(A1)'                                 ;tick labels only on top plots
     x2title = (multi[3] eq 0) ? x2title : ''                                   ;title only on top plots
     if multi[0] eq 2 and multi[2] gt 0 then x2title=''                         ;if 2 horizontal plots, only one title
     if multi[0] eq 2 and multi[2] eq 0 then x2title=hugewhitespace+x2title     ;space out that title to put it in the center

     ;right axis
     rytf = '(A1)'                                            ;blank by default
     y2title = (multi[2] eq multi[0]-1) ? y2title : ''        ;no title unless on right side
     if multi[0] eq 2 and multi[2] eq 1 then rytf='(I)'       ;if right of 2 horizontal plots, label right side
     if multi[0] eq 2 and multi[2] eq 1 then rytitle=ytitle   ;if right of 2 horizontal plots, title right side
  endif

  megaplot, xrange=xr, yrange=yrange, xstyle=1+8, ystyle=8+1, /xlog, $
            xtitle=xtitle1, ytitle=ytitle1, xtickformat=xtf1, xminor=2, ytickformat=ytf

  ;plot measurements
  megaplot, wveffobs[f], abmags[f], yerr=abmagerr[f], xerr=wverrobs[f], /noerase
  if nnf gt 0 then $
   megaplot, wveffobs[nf], abmags[nf], yerr=abmagerr[nf], xerr=wverrobs[nf], errcolor=15, /noerase

  ; plot top x-axis
  if max(xr/(1+z)) lt 8000. then begin
     xtickv = [1500.,2000.,3000.,4000.,6000.]
     xticks = 3
     xticklabelv = [1000,2000,4000,6000]
     xticklabels = 3
  endif
  if max(xr/(1+z)) ge 8000. and max(xr/(1+z)) le 10000. then begin
    xticklabelv = [1000.,2000.,4000.,6000.,8000.]
    xticklabels = 4
    xtickv = [1,2,4,6,8,10]*1000.
  endif
  if max(xr/(1+z)) gt 10000. then begin
    xticklabelv = [1000.,2000.,4000.,6000.,10000.,20000.]
    xticklabels = 5
    xtickv = [1,2,4,6,8,10,20]*1000.
;    xticks = 7
  endif
  axis, xaxis=1, xrange=xr/(1.+z), xstyle=1, xtitle=x2title, /xlog, $
        xticks=xticklabels, xtickv=xticklabelv, $
        xtickformat=xtf, xminor=1
  axis, xaxis=1, xrange=xr/(1.+z), xstyle=1, /xlog, $
        xticks=n_elements(xtickv), xtickv=xtickv, $
        xtickformat='(A1)', xminor=2, /noerase
    
  ; plot right y-axis
  if n_elements(multi) eq 0 then begin
    ;flux axis

    ytickv=ytv
  
    axis, yaxis=1, yrange=[minflux,maxflux], ystyle=1, yminor=9, ytickformat='(A1)', /ylog
    
    axis, yaxis=1, yrange=[minflux,maxflux], ystyle=1, ytitle='!17'+y2title,$
        /ylog, yticks=n_elements(ytickv), ytickv=ytickv, ytickformat=rytf,yminor=0
    ;axis, yaxis=1, yrange=[minflux,maxflux], ystyle=1, $
    ;    /ylog, yticks=n_elements(ytickv)+1, ytickv=ytickv, ytickformat='(I)',yminor=10
     
    ;if n_elements(!p.multi[1] > 1 or !p.multi[2] > 1) then device, font_size=9
  

    ;plot, wveffobs, fluxcor, xrange=[maxwv,minwv], xstyle=1, yrange=[minflux,maxflux], ystyle=1, /xlog, /ylog, xtitle='!  17'+xtitle, ytitle='!17'+ytitle, /nodata,$
    ;  title = '!17'+title+'!17', $
    ;  xticks = n_elements(xtv)-1, $
    ;  xtickv = xtv, $
    ;  xtickformat = '(I)', $
    ;  yticks = n_elements(ytv)-1, $
    ;  ytickv = ytv
  endif else begin
    if multi[0] gt 2 then $
      axis, yaxis=1, ytickformat='(A1)', yticks=1, yminor=1 ;null axis if more than 3 horizontal plots
      if multi[0] eq 2 then axis, yaxis=1, yrange=yrange, ytickformat=rytf, yminor=10, ytitle=rytitle, ystyle=1
  endelse

if nofit ne 0 then begin
  goto, printlabels ; gratuitous hack
endif

;axis, xaxis = 1, xrange=[25000, 5000]/(1.0+z), xtickv=xtv/(1.0+z), /xlog

!p.psym = 0

if n_elements(linestyle) gt 0 then extfitlinestyle = linestyle else extfitlinestyle = 0
if keyword_set(overplot) then extfitlinestyle = 1
unextlinestyle = 3
bestunextlinestyle = 2
if n_elements(line2style) gt 0 then model2linestyle = line2style else model2linestyle = 1

;plot fitted curve

valid = where(curvewav/(1.+z) gt 1000.)

oplot, curvewav[valid], ujy2mag(curveflux[valid],'AB'), psym=0, linestyle=extfitlinestyle
if unextcurve and av gt 0.01 then $
  oplot, curvewav, ujy2mag(unextcurveflux,'AB'), psym=0, linestyle=unextlinestyle
if ext0fit then begin
  notdamped = where(curvewav/(1+z) gt 1217)
  oplot, curvewav[notdamped], ujy2mag(fitext0curveflux[notdamped],'AB'), psym=0, linestyle=bestunextlinestyle
endif
if n_elements(ext2model) gt 0 then begin
  oplot, curvewav, ujy2mag(curveflux2,'AB'), psym=0, linestyle=model2linestyle
endif




fixedstr = '' ; ' (fixed)'

pm = ' !19'+string(177B) + '!17 '

betastr = fpr(beta+0.0001, 2.2) + $
                      (betafixed ? fixedstr : pm +fpr(betaerr,1.2))
avstr   = fpr(av,1.2) + (avfixed ? fixedstr : pm +fpr(averr,1.2))
if extmodel eq 'none' then avstr = '0 (fixed)'
rvstr   = clip(rv,4) + (rvfixed ? fixedstr : pm +fpr(rverr,1.2))
bumpstr = clip(3.23*bump,4)+(bumpfixed ? fixedstr : pm +fpr(3.23*bumperr,1.2))
uvstr   = clip(0.41*uv,4) + (uvfixed ? fixedstr : pm +fpr(0.41*uverr,1.2))


if n_elements(ext2model) gt 0 then begin
  beta2    = sol2[where(parname eq 'beta')]
  beta2err = parerr2[where(parname eq 'beta')]
  av2      = sol2[where(parname eq 'a')]
  av2err   = parerr2[where(parname eq 'a')]
  av2str   = fpr(av2,1.2) + (avfixed ? fixedstr : pm +fpr(av2err,1.2))


  beta2str = fpr(beta2+0.0001, 2.2) + $
                        (betafixed ? fixedstr : pm +fpr(beta2err,1.2))
  av2str   = fpr(av2,1.2) + (avfixed ? fixedstr : pm +fpr(av2err,1.2))
endif


  ;averrstr = '!S+ ' + fpr(averr,1.3) + ' !R_'
  ; rvstr = fpr(rv,1.3)+ '!S+ ' + fpr(rverr,1.3) + ' !R_'

chisqstr = clip(chisq,4) + ' / ' + clip(dof)
if chisq lt 0.01 then chisqstr = '0' + ' / ' + clip(dof)



if not(keyword_set(nolegend)) then begin
if keyword_set(shortcaption) then begin

   ; "Short" caption
 
   lx = !p.position[0] + xmargin
   ly = !p.position[1] + ymargin + spacing * 3
   if ext0fit then ly = ly + spacing
   if n_elements(ext2model) gt 0 then ly = ly + spacing
   lsp = spacing
   
   astr = 'A!DV!N = '
   a2str = astr
   extstring = '' ; general - FM2, etc.
   ext2string = ''
   if extmodel eq 'none' then extstring = ''
   if extmodel eq 'smc' then extstring = '(SMC)'
   if extmodel eq 'lmc' then extstring = '(LMC)'
   if extmodel eq 'mw' then extstring =  '(MW)'
   if extmodel eq 'calz' then extstring = '(Calz)'
   if extmodel eq 'maio' then extstring = '(SN)'
   if extmodel eq 'maio' then astr = 'A!D3000!N = '
   if n_elements(fixrv) gt 0 then begin
     if extmodel eq 'smc' and fixrv ne 2.73 then extstring = '(SMC, R!DV!N='+fpr(fixrv,1.1)+')'
     if extmodel eq 'lmc' and fixrv ne 3.2 then extstring = '(LMC, R!DV!N='+fpr(fixrv,1.1)+')'
     if extmodel eq 'mw' and fixrv ne 3.1 then  extstring =  '(MW, R!DV!N='+fpr(fixrv,1.1)+')'
   if extmodel eq 'calz' and fixrv ne 4.0 then extstring = '(Calz, R!DV!N='+fpr(fixrv,1.1)+')'
   endif
   if n_elements(ext2model) gt 0  then begin
    if ext2model eq 'none' then extstring = ''
    if ext2model eq 'smc' then ext2string = '(SMC)'
    if ext2model eq 'lmc' then ext2string = '(LMC)'
    if ext2model eq 'mw' then  ext2string =  '(MW)'
    if ext2model eq 'maio' then ext2string = '(SN)'
    if ext2model eq 'maio' then a2str = 'A!D3000!N = '
    if ext2model eq 'calz' then ext2string = '(Calz)'
   endif 
   

   if beta lt 0 then betastr = clip(fpr(beta,2.1)) else betastr = clip(fpr(beta,1.2)) 
   if n_elements(ext2model) gt 0  then begin
     if beta2 lt 0 then beta2str = clip(fpr(beta2,2.1)) else beta2str = clip(fpr(beta2,1.2)) 
   endif

   plots, [lx,lx+0.05], [ly,ly]+0.01,       /norm, linestyle=extfitlinestyle
   ;xyouts, lx+0.33, ly,       /norm, extstring
   xyouts, lx+0.19, ly,       /norm, astr+fpr(av,1.2) + ' ' + extstring
   xyouts, lx+0.06, ly,       /norm, '!4b!17 = '+betastr


   if n_elements(ext2model) gt 0 then begin
     plots, [lx,lx+0.05], [ly,ly]+0.01+lsp*1, /norm, linestyle=model2linestyle
     ;xyouts, lx+0.33, ly+lsp*1, /norm, ext2string
     xyouts, lx+0.19, ly+lsp*1, /norm, a2str+fpr(av2,1.2) + ' ' + ext2string
     xyouts, lx+0.06, ly+lsp*1, /norm, '!4b!17 = '+beta2str 
   endif

   if ext0fit then begin
     plots, [lx,lx+0.05], [ly,ly]+0.01-lsp*1, /norm, linestyle=bestunextlinestyle
     xyouts, lx+0.19, ly-lsp*1, /norm, 'A!DV!N = 0'
     xyouts, lx+0.06, ly-lsp*1, /norm, '!4b!17 = '+fpr(fitext0sol[betapar]+0.0001,1.2)
   endif

   if unextcurve and av gt 0.01 then begin 
     lx = !p.position[0] + (0.58-0.13)
     ly = !p.position[3] - 0.06 * (3.8/ys)
     plots, [lx,lx+0.05], [ly,ly]+0.01-lsp*0, /norm, linestyle=unextlinestyle
     xyouts, lx+0.19, ly-lsp*0, /norm, 'A!DV!N = 0'
     xyouts, lx+0.06, ly-lsp*0, /norm, '!4b!17 = '+clip(fpr(beta,1.2))
   endif


endif else begin

   ; Standard caption

   lux = !p.position[0] + xmargin

   lx = !p.position[2] - xmargin - 15 * xnorm_ch_size
   ly = !p.position[3] - ymargin - ynorm_ch_size * 1.25
   lsp = spacing
   nl = 0

   if extmodel eq 'none' and keyword_set(ext2model) then begin
     hold = lux
     lux = lx     ; swap locations
     lx = hold
   end

   xyouts, lx, ly-lsp*nl, /norm, '!4b!17 = ' + betastr  & nl = nl + 1
   if extmodel ne 'maio' then begin
       xyouts, lx, ly-lsp*nl, /norm, 'A!DV!N = ' + avstr    & nl = nl + 1
   endif else begin
       xyouts, lx, ly-lsp*nl, /norm, 'A!D3000!N = ' + avstr    & nl = nl + 1
   endelse
   if (((avfixed eq 0) or av gt 0.0) and extmodel ne 'maio' and extmodel ne 'none') then begin
     xyouts,lx,ly-lsp*nl, /norm, 'R!DV!N = ' + rvstr    & nl = nl + 1
   endif
   if extmodel eq 'fm' then begin
     xyouts, lx, ly-lsp*nl, /norm, 'C!D3!N = ' + bumpstr  & nl = nl + 1
     xyouts, lx, ly-lsp*nl, /norm, 'C!D4!N = ' + uvstr    & nl = nl + 1
   endif else begin
    if n_elements(multi) eq 0 and extmodel ne 'none' then begin
     if extmodel eq 'smc' then xyouts, lx, ly-lsp*nl, /norm, 'SMC extinction'
     if extmodel eq 'lmc' then xyouts, lx, ly-lsp*nl, /norm, 'LMC extinction'
     if extmodel eq 'calz' then xyouts, lx, ly-lsp*nl, /norm, 'Calzetti extinction'
     if extmodel eq 'maio' then xyouts, lx, ly-lsp*nl, /norm, 'Maiolino extinction'
     nl = nl + 1
    endif
   endelse
   xyouts, lx, ly-lsp*nl, /norm, '!4c!17!U2!N/dof = ' + chisqstr


   if keyword_set(ext0fit) then begin
      if beta gt 2 then lux = lux - 0.11
      ubetastr = fpr(fitext0sol[betapar]+0.0001,2.2) + pm +fpr(fitext0perr[betapar],1.2)
      uavstr = '0.0 (fixed)'
      uchisqstr = clip(fitext0chisq,4) + ' / ' + clip(fitext0dof)
      xyouts, lux, ly-lsp*0, /norm, '!4b!17 = ' + ubetastr
      xyouts, lux, ly-lsp*1, /norm, 'A!DV!N = ' + uavstr
      xyouts, lux, ly-lsp*2, /norm, '!4c!17!U2!N/dof = ' + uchisqstr
      plots, [lux+0.03,lux+0.18], [ly+spacing,ly+spacing], /norm, linestyle=bestunextlinestyle

      plots, [lx+0.03,lx+0.18], [ly+spacing,ly+spacing], /norm, linestyle=extfitlinestyle
   endif
   if n_elements(ext2model) then begin
      chisq2str = clip(chisq2,4) + ' / ' + clip(dof2)
      xyouts, lux, ly-lsp*0, /norm, '!4b!17 = ' + beta2str
      xyouts, lux, ly-lsp*1, /norm, 'A!DV!N = ' + av2str
      if ext2model eq 'smc' then  xyouts, lux, ly-lsp*2, /norm, 'SMC extinction'
      if ext2model eq 'lmc' then  xyouts, lux, ly-lsp*2, /norm, 'LMC extinction'
      if ext2model eq 'calz' then xyouts, lux, ly-lsp*2, /norm, 'Calzetti extinction'
      if ext2model eq 'maio' then xyouts, lux, ly-lsp*2, /norm, 'Maiolino extinction'
      xyouts, lux, ly-lsp*3, /norm, '!4c!17!U2!N/dof = ' + chisq2str
      plots, [lux+0.03,lux+0.18], [ly+spacing,ly+spacing], /norm, linestyle=model2linestyle

      plots, [lx+0.03,lx+0.18], [ly+spacing,ly+spacing], /norm, linestyle=extfitlinestyle
   endif


endelse
endif


if not(keyword_set(nolegend)) then begin
  l2x = !p.position[0] + xmargin*0.85
  l2sp = spacing*0.85
  l2y = !p.position[1] + ymargin + l2sp*0.25;*0.75
  ypos = l2y
  if n_elements(multi) eq 0 then begin
    xyouts, l2x, ypos, /norm, 'Host redshift z = ' + clipzero(z), size=0.85
    ypos = ypos + l2sp
  endif
  if n_elements(galav) gt 0 and n_elements(multi) eq 0 then begin
    xyouts, l2x, ypos, /norm, 'Corrected for Galactic A!DV!N =' + string(galav,format='(F5.2)'), size=0.85
    ypos = ypos + l2sp
  endif
  if n_elements(time) gt 0 then begin
    if time gt 0.5 then timestr = clipzero(time) + ' day' $
                   else timestr = clip(round(time*24.*3600.))+' s'
    xyouts, l2x, ypos, /norm, 't = '+timestr, size=0.85
    ypos = ypos + l2sp
  endif
  if n_elements(comment) gt 0 then begin
    xyouts, l2x, ypos, /norm, comment
    ypos = ypos + l2sp
  endif
  if n_elements(morecomment) gt 0 then begin
    xyouts, l2x, ypos, /norm, morecomment
    ypos = ypos + l2sp
  endif

endif

printlabels:
labels = 1
if keyword_set(labels) then begin
  for i = 0, n_elements(filters)-1 do begin
    col = cols[i]
    xpos = wveffobs[i]/1.03
    ypos = abmags[i]-0.65*spacing*abs(yrange[1]-yrange[0])
 ;   if xpos/xr[1] lt 2 and filters[i] ne 'clear' then begin
 ;     xpos = xpos * 1.1
 ;   ;  ypos = ypos - 0.04
 ;   endif
 ;   nearcurve = (where(    abs(wveffobs[i] - curvewav) eq $
 ;                     min(abs(wveffobs[i] - curvewav)))) [0]
 ;     ;print, '-', filters[i],  abmags[i] - ujy2mag(curveflux[nearcurve],'AB')
 ;   if abmags[i] - ujy2mag(curveflux[nearcurve],'AB') ge 0.15 then begin
 ;     xpos = xpos * 1.10
 ;   endif  
  ;  if xpos gt xr[1]*1.15 and xpos lt xr[0] then $
      xyouts, xpos, ypos, filters[i], color=col
  endfor
endif

if nofit then begin
  psclose 
  set_plot, initdev
  return
endif


; RESIDUALS PLOT

if plotresiduals then begin
   !p.position = respos


   rightytitle = 'flux ratio'
   residual = ujy2mag(fluxcor,'AB') - ujy2mag(extpw(wveffobs,sol),'AB')

   resyrange = [0.5,-0.5]
   resyrangeright = mag2ujy(resyrange, 'AB') / mag2ujy(0.0, 'AB')
   resytickv = [-0.5,0,0.5]
   resyminor = 5

   resrightyminor = 5
   resrightytickv = [0.5,1.0,1.5,2.0]

  errcolor = intarr(n_elements(residual))
  if nnf gt 0 then errcolor[nf] = 15

  megaplot, wveffobs, residual, yerr=abmagerr, xerr=wverrobs, psym=-1, errcolor=errcolor, $
            xrange=xr, xstyle=1, ystyle=1, /xlog, $
            xtitle=xtitle, xtickformat=xtf, xminor=2, $ 
            ytitle='residual (mag)',yrange=resyrange, ytickv=resytickv, yminor=resyminor, $
            rightynormpos=resyrange,rightyrange=resyrangeright,rightystyle=1,/rightylog, $
            rightytitle=rightytitle, ryminor=resrightyminor, rytickv=resrightytickv, xticklen=0.04
  ;axis, xaxis=1, xrange=xr, xstyle=1, xtitle='', /xlog, $
  ;          xtickformat='(A1)', xminor=2


  oplot, xr, [0,0]

  if ext0fit then $
      oplot, curvewav[notdamped], (ujy2mag(fitext0curveflux,'AB')-ujy2mag(curveflux,'AB')) [notdamped], linestyle=bestunextlinestyle
  if n_elements(ext2model) gt 0 then begin
      valid = where(curveflux2 gt 0)
      oplot, curvewav[valid], ujy2mag(curveflux2[valid],'AB')-ujy2mag(curveflux[valid],'AB'), linestyle=model2linestyle
  endif

  ;diffcurve = ujy2mag(curveflux2,'AB')-ujy2mag(curveflux,'AB')
  ;oplot, curvewav, diffcurve, linestyle=1


  ;replot axes again
  plot, [0], [0],  xminor=1, xticks=1, xtickformat='(A1)', yminor=1, yticks=1, ytickformat='(A1)', /noerase



endif

if n_elements(multi) eq 0 then begin   ; close the plot if not in multi mode.
  psclose 
  set_plot, initdev
endif else begin
  multi[3] = multi[3] + 1
  if multi[3] eq multi[1] then begin
    multi[3] = 0
    multi[2] = multi[2] + 1
  endif
endelse


; Broadband plot for X-rays...
if ctxray gt 0 then begin
  initdev = !d.name
  set_plot, 'ps'
  !p.multi = 0
  !p.position = 0
  device, filename='broadbandplot.ps'
  device, /inches, xsize=xs, ysize=4
  device, /inches, yoffset=4.5+(5.5-ys)

    bbcurvewav = [curvewav, 1000, 100, 10]
    bbunextcurveflux = extpw(bbcurvewav, unextsol, /noigm)

    megaplot, [wveffobs, wveffobsx], [fluxcor, fluxcorx], yerr=[fluxerr, fluxerrx], xerr=[wverrobs, filtwv('keV')/2+fltarr(n_elements(fluxerrx))], xrange=[30000,5],   xstyle=1, /xlog, /ylog, xtitle=xtitle, ytitle=y2title
    oplot, bbcurvewav, bbunextcurveflux
    oplot, curvewav, unextcurveflux, linestyle=1
    oplot, curvewav, curveflux, linestyle=2
  psclose
  set_plot, initdev
endif 



;  INFORMATION PRINTED TO CONSOLE

;print, '             F0 = '+clip(f0,5) + ' +/- ' + clip(ef0,5)
print, '            beta = '+clip(beta,5) + ' +/-' + fpr(betaerr,2.3)
if av ne 0.0 then begin
print, '             A_V ='+fpr(av,2.3)+ ' +/-' + fpr(averr,2.3)
endif else begin
  if n_elements(fixav) eq 0 then $
    print, '             A_V <'+fpr(averr,2.3) $
  else $
    print, '             A_V ='+fpr(av,2.3)+ ' (fixed)'
endelse
if rvfixed eq 0 or strmid(extmodel,0,2) eq 'fm' or extmodel eq 'reichart' or extmodel eq 'reichartg' or extmodel eq 'ccm' or extmodel eq 'lmc2' then $
print, '             R_V ='+fpr(rv,2.3)+ ' +/-' + fpr(rverr,2.3)
;if av gt 1e-5 and extmodel ne 'maio' then $
;print, '          E(B-V) ='+fpr(av/rv,2.3)+ ' +/-' + fpr(averr/rv,2.3)
if strmid(extmodel,0,2) eq 'fm' or extmodel eq 'reichart' or extmodel eq 'reichartg' then begin
  print, '             C1 ='+fpr(c1,2.3)+ ' +/-' + fpr(c1err,2.3)
  print, '             C2 ='+fpr(c2,2.3)+ ' +/-' + fpr(c2err,2.3)
  print, '             C3 ='+fpr(c3,2.3)+ ' +/-' + fpr(c3err,2.3)
  print, '             C4 ='+fpr(c4,2.3)+ ' +/-' + fpr(c4err,2.3)
 if extmodel eq 'fm2g' or extmodel eq 'fm1g' or extmodel eq 'reichartg' then $
  print, '          gamma =',fpr(gamma,2.3)+ ' +/-' + fpr(gammaerr,2.3)
endif
if strmid(extmodel,0,2) eq 'li' then begin
  print, '             C1 ='+fpr(c1,3.2)+ ' +/-' + fpr(c1err,3.2)
  print, '             C2 ='+fpr(c2,3.2)+ ' +/-' + fpr(c2err,3.2)
  print, '             C3 ='+fpr(c3,3.2)+ ' +/-' + fpr(c3err,3.2)
  print, '             C4 ='+fpr(c4,3.2)+ ' +/-' + fpr(c4err,3.2)
endif
print, '           chi^2/dof = '+chisqstr

if n_elements(multi) eq 0 then begin
  print, 'Plotted to '+plotfile
endif

if n_elements(storefitfilters) gt 0 then fitfilters = storefitfilters
x = storex
fluxerr = storefluxerr


end