! @(#)spdblend.prg	17.1.1.1 (ES0-DMD) 01/25/02 17:56:03
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! PROCEDURE DBLEND
! Execute as @@ P2:DBLEND p1 p2 p3 p4 p5 p6 p7
! p1 : the name of a wavelength-calibrated, one-dimensional spectrum
! p2 : the name of the output file in which the fit is to be stored (def MDTMP1)
! p3 : name of output table file (def MDFIT)
! p4 : 1 to 3 letter fit option (def ABG)
!     available options:
!      first letter:
!       A: all fit parameters can be varied to get the best fit
!       W: center wavelengths held fixed at input values, heights and widths of
!          lines can vary
!       F: widths of lines held fixed at given values, centers and heights of
!          lines can vary
!       S: all fit parameters can vary, but all lines are required to have the
!          same width in the fit.
!      second letter:
!       B: non-weighted fit (ie. all points have equal weight)
!       P: Poisson-noise weighted fit (weight(i) = 1/datavalue(i))
!      third letter:
!       G: fit Gaussian curves: 
!            y(x) = max*exp(-0.5*((x-cent)/sigma)**2)
!       L: fit Lorentz curves: 
!            y(x) = max * (fwhm/2)**2 / ((x-cent)**2+(fwhm/2)**2)
!       E: fit power law curve (not available yet)
! p5 : continuum option (def P)
!     available options:
!       P: continuum is defined by two points
!       R: continuum is defined by two regions, bounded by four points; the
!          average pixel value in each region is calculated, and these two 
!          points are used to define the continuum
! p6 : input option (def C)
!     available options:
!       C: all initial parameters input interactively using cursor
!       T: all initial parameters read from an input table
!       B: available for A, W and S fit options.  For fit option = A, user
!          will be prompted for starting guesses at FWHM; centers are entered
!          interactively.  For fit option = W, input wavelengths are in table, 
!          remaining parameters are input interactively.  For fit option = S, 
!          input line full width half max is input in P7, remaining parameters 
!          are input interactively.
! p7 : input table name (def TTEMP, used for p6 = C); see also under B input
!       option
!
! This routine provides de-blending of spectral lines by means of a 
! simultaneous fit of up to six lines to a selected spectral region.
! Line profiles may be either Gaussian or Lorentzian in profile, or have
! a power law fall-off on either side of the maximum.
!
! Adapted from TVG program TVDBLE (command DBLEND) Jan. 1990 by Marylin Bell
! Added Lorentzian curve option Sep. 1990 (M.Bell)
!
! Modified to run on Nov 93 MIDAS at 4.94 
!------------------------------------------------------------------------------
CROSSREF INFILE FITIM FITPAR METHOD CONTIN INPUT INTAB
!
! set up inputs
!
DEFINE/PAR P1 ? IMAGE "Input spectrum :"
DEFINE/PAR P2 mdtmp1 IMAGE  !output fit file
DEFINE/PAR P3 mdfit TABLE
DEFINE/PAR P4 ABG CHARACTER
DEFINE/PAR P5 P CHARACTER
DEFINE/PAR P6 C CHARACTER
DEFINE/PAR P7 ttemp
!
DEFINE/LOCAL INFRAM/C/1/12 'P1'
DEFINE/LOCAL ZLIM/I/1/2 0,0            !limits for zoomed plot
DEFINE/LOCAL FITOPT/C/1/3 'P4'         !fit options chosen
DEFINE/LOCAL STATUS/I/1/1 0            !error status indicator
DEFINE/LOCAL N/I/1/1 0                 !counter
DEFINE/LOCAL CN/I/1/1 0                !counter for line changes
DEFINE/LOCAL FN/I/1/1 0                !counter for color changes
DEFINE/LOCAL LFILE/C/1/12              !individual line profile filename
DEFINE/LOCAL ISIT/I/1/1 0              !yes or no indicator for tests
DEFINE/LOCAL FWHM/R/1/6 0,0,0,0,0,0    !FWHM's to be entered at prompt
DEFINE/LOCAL NAX/I/1/1 0
DEFINE/LOCAL NP/I/2/1 0,0
DEFINE/LOCAL ST/R/2/1 0,0
!
! convert input 2D frame to 1D scratch frame
!
COPY/DK 'INFRAM' NAXIS NAX
IF NAX .NE. 1 THEN
  COPY/DK 'INFRAM' NPIX NP
  COPY/DK 'INFRAM' STEP ST
  IF NP(1) .EQ. 1 THEN
    EXTRACT/LINE dbtmp = 'INFRAM'[@1,@1:<,>] 'ST(2)'
  ELSE
    EXTRACT/LINE dbtmp = 'INFRAM'[<,>:@1,@1]
  ENDIF
  WRITE/KEY INFRAM dbtmp
ENDIF
!
! Note for input option = B, input table given will have to be transformed
!  to correct format for DBLEND; the transformed file will be names TTEMP.
!  Therefore, TTEMP is put in keyword IN_B and input table name given is
!  stored in BTAB.
!
IF P6(1:1) .EQ. "B" THEN
  WRITE/KEY IN_B "ttemp"
  DEFINE/LOCAL BTAB/C/1/12 'P7'
ELSE
  WRITE/KEY IN_B 'P7'
ENDIF
!
! write output filenames to be used to corresponding keywords
!
WRITE/KEY OUT_A 'P2'
WRITE/KEY OUT_B 'P3'
!
! initializations
!
SET/GRAP COLOUR=1 LTYPE=1 PMODE=1
CREATE/GRAPHIC ? 850,550,1,1
!
! delete old line files with default names, if they exist
!
DO N = 1 6 
LFILE = "mdtl{n}.bdf"
ISIT = M$EXIST(LFILE)
IF ISIT .NE. 0 DELETE/IMA 'LFILE(1:8)' NO
ENDDO
!
! for table input, skip plotting and interactive inputs
!
IF P6(1:1) .EQ. "T" GOTO DOIT
!
! clean up descriptor that may be left over by improper exit last time
!
ISIT = M$EXISTD(INFRAM,"COUT")
!write/out "isit = 'isit'"
IF ISIT .NE. 0 DELETE/DESCR 'INFRAM' COUT
!
! plot input spectrum to screen, get limits of area to zoom in on
!
FIND/MINMAX 'INFRAM'
CUTS/IMA 'INFRAM' 0,0
PLOT 'INFRAM'
!
! if line wavelengths are already in a table, overplot the given positions
!  for the user's benefit
!
IF FITOPT(1:1) .EQ. "W" THEN
  IF P6(1:1) .EQ. "B" OVERPLOT/IDENT 'P7'
ENDIF
WRITE/OUT
!
! get limits of area to zoom in on
!
WRITE/OUT " Use the cursor to enter beginning and end of region to zoom"
GET/GCURSOR COUT,DES N 2
COPY/DK 'INFRAM' COUT/R/3/1 ZLIM/I/1/1
COPY/DK 'INFRAM' COUT/R/10/1 ZLIM/I/2/1
DELETE/DESCR 'INFRAM' COUT
!
! plot chosen region
!
FIND/MINMAX 'INFRAM'[@'ZLIM(1)',@1:@'ZLIM(2)',@1]
PLOT 'INFRAM' @1 @'ZLIM(1)',@'ZLIM(2)'
IF FITOPT(1:1) .EQ. "W" THEN
  IF P6(1:1) .EQ. "B" OVERPLOT/IDENT 'P7'
ENDIF
WRITE/OUT
!
! do cursor parameter input
!
RUN STD_EXE:DBLGC
!
! do typed-in input
!
IF P4(1:1) .EQ. "A" THEN
  IF P6(1:1) .EQ. "B" INQUIRE/KEY FWHM " enter FWHM's of lines marked:"
ENDIF
!
DOIT:
!
! get continuum points into correct format for DBLEND
!
IF P5(1:1) .EQ. "R" RUN STD_EXE:AVECON
!
! get the rest of the table of parameters into correct format for DBLEND
!
IF P6(1:1) .EQ. "B" THEN
  RUN STD_EXE:MAKTAB
  WRITE/KEY IN_B "ttemp"
ENDIF
!
! error check
!
IF STATUS .NE. 0 GOTO END
!
! finally!
!
RUN STD_EXE:DBLEND
!
! (due to problems with CUNIT descriptor in MIDAS)
!
COPY/DD 'P1' CUNIT 'P2' CUNIT
!
! Plot fit on top of (already plotted) spectrum
!
FIND/MINMAX 'P2'
COPY/DD 'P2' LHCUTS/R/3/2 'P2' LHCUTS/R/1/2
IF P6(1:1) .EQ. "T" THEN
  PLOT 'P2'
  OVERPLOT 'P1'
ELSE
  OVERPLOT 'P2' 
ENDIF
!
! Overplot individual line profiles
!
!WRITE/OUT " Ignore any FILE NOT FOUND error messages you see next"
N = 1
LFILE = P2(1:3)//"l'N'.bdf"
ISIT = M$EXIST(LFILE)
IF ISIT .EQ. 0 GOTO NOTFOUND
!
FOUND:
CN = N+1
FN = CN
IF CN .GE. 5 FN = CN+1
SET/GRAP COLOUR='FN' LTYPE='CN'
FIND/MINMAX 'LFILE'
COPY/DD 'LFILE' LHCUTS/R/3/2 'LFILE' LHCUTS/R/1/2
OVER/ROW 'LFILE'
N = N+1
LFILE = P2(1:3)//"l'N'.bdf"
ISIT = M$EXIST(LFILE)
IF ISIT .NE. 0 GOTO FOUND
!
NOTFOUND:
!
! Write HISTORY descriptor to new frame
!
!WRITE/DESCR 'P2' HISTORY/C/-1/80 "DBLEND fit to frame 'P1' "
!
! Clean-up operation
!
FIND/MINMAX 'P1'
SET/GRAP
IF NAX .NE. 1 THEN
  DELETE/IMA dbtmp NO
ENDIF
!
WRITE/OUT
WRITE/OUT " fit image stored in file 'P2'"
WRITE/OUT " fit parameters stored in table 'P3'"
END:
WRITE/OUT "DBLEND finished"