! @(#)irstarflux.prg	17.1.1.1 (ES0-DMD) 01/25/02 17:53:01
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!.COPYRIGHT   (C) 1992 European Southern Observatory
!.IDENT       .prg
!.AUTHOR      E. Oliva,  Firenze, Arcetri
!.KEYWORDS    Spectroscopy, IRSPEC
!.PURPOSE     Execute command  STANDARD/IRSPEC
!
!                used to create a table file with :wl and :flux
!                from an ASCII file containing the star flux at
!                a few wavelengths. Interpolate using either a BB
!                or a polynomial in log-log.
!                Needs irsbbfit.for
!
!.VERSION     1.0    Creation     02.09.1992   E. Oliva
!
!-------------------------------------------------------
!
! Tested O.K. 23-sep-92 Tino
!
crossref in out method degree step limits units plot
def/param p1 ? cha "Name of ASCII file containing star fluxes ? : "
def/param p2 ? tab "Name of output table ? : "
def/param p3 ? cha "Method for interpolation (bb, poly, spline) ? : "
def/param p4 2 num "Polynomial degree "  1,10
def/param p5 .01 num "Wavelength step in out_table " 
def/param p6 0,1e7 num "Optional limits start,end-wl " 
def/param p7 m cha "Units (m=micron a=Angstrom)"
def/param p8 1 num "Plot option (0,1)"
!
def/loc incmd/c/1/132 " " 
copy/key MID$LINE/C/1/132 incmd/C/1/132
!def/loc incmd/c/1/132 "{MID$LINE}"
!
def/loc met/c/1/80
def/loc deg/i/1/1 0
def/loc nwl/i/1/1 0 ? +lower
def/loc step/r/1/1 0.
def/loc tlim/r/1/2 0,0
def/loc i/i/1/1 0
def/loc j/i/1/1 0
def/loc table/c/1/60
def/loc temp/r/1/1 0.
def/loc flux0/r/1/1 0.
def/loc limits/r/1/2 0,1e7
def/loc units/i/1/1 0

!VERIFY/IRSPEC 'p1'
@ creifnot 1

met = m$lower(p3)
if met(1:1) .ne. "s" then
  if met(1:1) .ne. "b" then
    if met(1:1) .ne. "p" then
      write/out Invalid method. Only bb, poly and spline are available.
      return/exit
    endif
  endif
endif
deg = 'p4'
step = 'p5'
write/keyw limits/r/1/2 'p6'
p7 = m$lower(p7)
if p7(1:1) .eq. "a" units = 1
crea/table &a 2 10 'p1' 
name/column &a #1 :wl
name/column &a #2 :flux
sort/table &a :wl
comp/table &a :log_wl = log10(:wl)
comp/table &a :log_flux = log10(:flux)
copy/dk &a.tbl tblcontr/i/4/1 nwl/i/1/1
tlim(1) = '&a,:wl,@1'
tlim(2) = '&a,:wl,@{nwl}'
if limits(1) .lt. tlim(1) limits(1) = tlim(1)
if limits(2) .gt. tlim(2) limits(2) = tlim(2)
i = (limits(2)-limits(1))/step+1
if i .gt. 3000 then
  write/out "The output table is going to be very large ({i} rows)"
  write/out "Make sure that you properly chose the limits and step"
  write/out "(expecially in you are working in Angstrom)"
  inqu/keyw j "Hit Return to continue, enter 1 to stop here : "
  if j .ge. 1 then
    return/end
  endif
endif
crea/table 'p2'  2 'i' null
comp/table 'p2' :wl = 'limits(1)'+(seq-1)*'step'
comp/table 'p2' :log_wl = log10(:wl)
if met(1:1) .eq. "s" then
  if nwl .le. 3 then
    write/out Need at least 4 input points for spline interpolation
    write/out .....Sorry.....
    return/end
  endif
  copy/table &a &b
  interp/tt &a :log_wl,:log_fit &b :log_wl,:log_flux 1 'deg'
  comp/table &a :fit = 10.**(:log_fit)
  comp/table &a :err = (:fit/:flux-1.)*100.
  interp/tt 'p2' :log_wl,:log_flux &b :log_wl,:log_flux 1 'deg'
  goto done
endif
if met(1:1) .eq. "p" then
  if nwl .le. deg then
    write/out Cannot fit a 'deg'-polynon to 'nwl' data points
    return/end
  endif
  regr/poly &a :log_flux :log_wl 'deg'
  save/regr &a poly
  comp/regr &a :log_fit = poly
  comp/table &a :fit = 10.**(:log_fit)
  comp/table &a :err = (:fit/:flux-1.)*100.
  comp/table 'p2' :log_flux = 'outputd(1)'+'outputd(2)'*:log_wl
  deg = deg+1
  do i = 3 'deg' 
    comp/table 'p2' :log_flux = :log_flux+'outputd({i})'*:log_wl**({i}-1)
  enddo
  goto done
endif
if met(1:1) .eq. "b" then
  if nwl .le. 2 then
    write/out Cannot fit a black-body to only 'nwl' data points
    return/end
  endif
  table = "&a"
  run STD_EXE:irsbbfit.exe
  temp = outputr(1)
  flux0 = outputr(2)
  write/out Fitted temperature='temp'
  comp/table &a :log_fit = log10('flux0')
  if units .eq. 0 then
    comp/table &a :log_fit = log10('flux0')-5*:log_wl
    comp/table &a :log_fit = :log_fit-log10((exp(14388/'temp'/:wl)-1))
  else
    comp/table &a :log_fit = log10('flux0')-5*:log_wl+20.
    comp/table &a :log_fit = :log_fit-log10((exp(1.4388e8/'temp'/:wl)-1))
  endif 
  comp/table &a :fit = 10.**(:log_fit)
  comp/table &a :err = (:fit/:flux-1.)*100.
  if units .eq. 0 then 
    comp/table 'p2' :log_flux = log10('flux0')-5*:log_wl
    comp/table 'p2' :log_flux = :log_flux-log10((exp(14388/'temp'/:wl)-1))
  else
    comp/table 'p2' :log_flux = log10('flux0')-5*:log_wl+20.
    comp/table 'p2' :log_flux = :log_flux-log10((exp(1.4388e8/'temp'/:wl)-1))
  endif 
  goto done
endif
done:
comp/table 'p2' :flux = 10.**(:log_flux)
copy/tt &a :wl 'p2' :wl_tab
copy/tt &a :flux 'p2' :flux_tab
copy/tt &a :fit 'p2' :flux_int
copy/tt &a :err 'p2' :err_percent
name/column 'p2' :err_percent f11.1
sel/table 'p2' seq.le.'nwl'
read/table 'p2' :wl_tab :flux_tab :flux_int :err_percent
sel/table 'p2' all
if 'p8' .gt. 0  then
  set/grap xaxis=auto yaxis=auto stype=2
  plot/table &a :log_wl :log_flux
  set/grap stype=0
  overplot/table 'p2' :log_wl :log_flux 
endif
copy/kd incmd/c/1/132 'p2'.tbl history/c/-1/132
return