      PROGRAM XPLOT
c	program to read information from a SF history file created 
c       by XSPEC for plotting using PGPLOT routines
c		version 1.0 fwj haberl   march 1988
c		version 1.1 fwj haberl   17-MAR-1989 
c		    pause command:    leave xplot with pause statement
c		    contour plotting: rdcont command to read steppar' package'
c				      plcont to make chi-square contours
c              version 2.0 fwj haberl   8-AUG-1990 
c                  uses PLT for plotting (not for contours)
c                  PLOT command returns to XPLOT prompt
c                  IPLOT command stays in PLT
c                  hardcopy command 17-AUG-1990 
c              version 2.1
c                  show command 1-NOV-1990
c              version 2.2
c                  on ALPHA/OSF
c
      INCLUDE '../inc/historysf.inc'
      INCLUDE '../inc/sitedef.inc'
      INCLUDE '../inc/writebuf.inc'
      INCLUDE 'xplot.inc'
      INTEGER*4 ierrsf, iunit, iparse, ios, lc, lseshd, lentmp
      INTEGER*4 len, lenact, ibeg, iend, ierr, i, iflag, idelim, nret
      INTEGER*4 maxhist, nhist
      PARAMETER (maxhist=60)
      CHARACTER*80 header, tmplte, history(maxhist)
      CHARACTER*128 string, tmpstr, mlabel
      CHARACTER*64 seshead, filenm
      CHARACTER*4 version
      INTEGER*4 dset_new, dset
      LOGICAL*4 qskip, qini, qmodel, qhist, qopen, qinter, qpause,
     &          qcontour
      INTEGER*4 savcom, nccom, nacom
      PARAMETER (savcom=20)
      CHARACTER*255 savstrg(savcom)
      REAL*4 chimin
      RECORD /his_steppar_sub1/stepsubx, stepsuby
c list of possible commands
      INTEGER*4 icom, maxcom
      PARAMETER (maxcom=42)
      CHARACTER*10 command(maxcom)
      CHARACTER*7 prompt
      CHARACTER*12 recfile
      DATA command/'chatter', 'reexecute', 'help', 'helpxspec',
     &     'history', 'clear', 'setplot', 'font', 'symbol', 'label',
     &     'plot', 'iplot', 'data', 'spectrum', 
     &     'kev_ang', 'rdmodel',
     &     'unfoldspec', 'divide', 'multiply', 'add', 
     &     'residuals','chiresid',
     &     'model', 'logarith', 'linear', 'rescale', 'rdcont', 
     &     'plcont', 'rdfit', 'distrib', 'fold', 'average', 
     &     'norm', 'wrqdp', 'exit', 'pause', 'viewport', 'newset', 
     &     'show', 'rdqdp','hardcopy','cpd'/
      DATA filenm/'none'/
      DATA prompt/'xplot> '/
      DATA version/'2.2'/
 
c initialize timer
      CALL gtinit

c write out the news file if one exists
      CALL trlog('XNEWS_DISABLE', 13, string, iparse)
      IF (iparse.EQ.0) THEN
c       writing news file not disabled
         INQUIRE (FILE=docdir//'xplot.new', IOSTAT=ios)
         CALL getlun(iunit)
         IF (ios.EQ.0) THEN
            CALL openwr(iunit, docdir//'xplot.new', 'old', ' ', ' ',
     &                   0, 1, ios)
            IF (ios.NE.0) THEN
               WRITE (*, *) ' *WARNING*: Unable to open News file'
            ELSE
               DO WHILE (ios.EQ.0)
                  READ (iunit, '(a)', IOSTAT=ios) string
                  IF (ios.EQ.0) 
     &               WRITE (*, '(1x,a)') string(:lenact(string))
               ENDDO
               CLOSE (iunit)
            ENDIF
         ENDIF
         CALL frelun(iunit)
      ENDIF
c entry in user record file
      recfile = 'xplot'
      CALL usrinf(recfile)
c create the startup banner
      string = 'XPLOT '//version
      lc = lenact(string) + 1
      string(lc:lc) = ' '
      lc = lc + 1
      CALL gettim(string(lc:))
      lc = lc + 8
      string(lc:lc) = ' '
      lc = lc + 1
      CALL getdat(string(lc:))
      lc = lc + 8
      seshead = string(:lc)
      lseshd = lc
      WRITE (*, '(/15x,a/)') string(:lc)
c default start values:
      indent = '      > '
      qmodel = .TRUE.
      qhist = .TRUE.
      qcontour = .FALSE.
      iucom = mxucmd
      CALL clear(dset, dset_new)

c check for the logical/environment variable PGPLOT_TYPE
c and set the plotting device. If no variable then set to none.

      CALL trlog('PGPLOT_TYPE', 11, tmpstr, lentmp)
      IF (lentmp.EQ.0) tmpstr = '/nu'
      string = 'device '//tmpstr(:lenact(tmpstr))
      iparse = 0
      CALL stplot(string, iparse, qini, dset, dset_new)
      if(tmpstr.ne.'/nu' .AND. qini)then
         tmpstr= ' Current plot device '//tmpstr(:lentmp)
      else
         tmpstr= ' Plot device not set, use "setplot dev" to set it'
      endif
      call xwrite(tmpstr,5)
      call xwrite(' ',5)

c initialize the indirect file suffix
      CALL xsetin('XPC')

c the initial command line is taken from the invocation line.  
      CALL xinird(' ', string, iparse)
c there was something on that command line:
      IF (iparse.EQ.0) GOTO 210

c general command input
 200  CONTINUE
 
c still more commands to reexecute
      IF (nacom.GT.1) THEN
         GOTO 215
      ENDIF
 
      CALL xnxtrd(prompt, string, iparse, ierr, *950)
 
      nacom = 0
c come from for reexecuting of command
 201  CONTINUE
 
      IF (ierr.NE.0) THEN
         WRITE (*, '(a)')
     &                ' *ERROR*:XPLOT: An error during the command read'
         iparse = len(string)
         GOTO 200
      ENDIF
c mis-spelled command name
 210  CONTINUE
      CALL xgtarg(string, iparse, ibeg, iend, qskip, iflag, idelim,
     &            *200, *200, *950)
      IF (.NOT.qskip) THEN
         CALL xmatch(string(ibeg:iend), command, maxcom, icom)
         IF (icom.EQ.0) THEN
            CALL xunids(string(ibeg:iend), command, maxcom, icom,
     &                  'XPLOT top level commands')
            CALL xinfix('Insert correct command:', string, iparse, ierr,
     &                  *200, *200)
            IF (ierr.LT.0) THEN
               GOTO 200
            ENDIF
            GOTO 210
         ENDIF
      ENDIF
 
      IF (icom.NE.1) CALL xcmsav(string, iparse, savstrg, savcom, nccom)
 
      GOTO (212, 215, 225, 226, 250, 255, 300, 325, 326, 327, 350, 351,
     &      400, 405, 410, 415, 417, 420, 421, 425, 430, 431, 
     &      450, 550, 560,
     &      600, 630, 631, 650, 654, 655, 670, 750, 940, 950, 951, 800,
     &      850, 880, 941, 900, 910), icom
      WRITE (*, '(a)') ' *ERROR*:XPLOT: command index out of range'
      GOTO 200

c chatter: set the chattyness level used for output verbosity
 212  CONTINUE
      CALL setcht(string, iparse)
      GOTO 200

c reexecute: execute previous command(s)
 215  CONTINUE
      CALL xcmget(string, iparse, savstrg, savcom, nccom, nacom, ierr)
      IF (nacom.GT.0) THEN
         GOTO 201
      ELSE
         GOTO 200
      ENDIF

c help: get help on the available commands
 225  CONTINUE
      string = string(iparse+1:)
      CALL f_ihf(docdir//'xplot.hlp', string)
      GOTO 200

c helpxspec: get help on the available commands in XSPEC
 226  CONTINUE
      string = string(iparse+1:)
      CALL f_ihf(docdir//'xspec.hlp', string)
      GOTO 200

c history: read in history filename
 250  CONTINUE
      CALL xgtstr(string, iparse, 1, 'file', 1, filenm, nret, iflag, -1)
      IF (iflag.LT.0) THEN
         GOTO 200
      ENDIF
      IF (filenm(1:4).EQ.'none') THEN
         IF (qopen) THEN
            CLOSE (iunit)
            CALL frelun(iunit)
            qopen = .FALSE.
         ELSE
            WRITE (*, '(2a)') indent, ' No history file open yet'
         ENDIF
         GOTO 200
      ENDIF
      CALL xtend(filenm, 'xhs')
      IF (qopen) THEN
         CLOSE (iunit)
         CALL frelun(iunit)
         qopen = .FALSE.
      ENDIF
      ierrsf = 0
      CALL getlun(iunit)
      CALL opnrsf(filenm, iunit, 'XSPEC histry', header, tmplte, nhist,
     &            ierrsf)
      IF (ierrsf.NE.0) THEN
         CALL frelun(iunit)
         GOTO 200
      ELSE
         qopen = .TRUE.
      ENDIF
      WRITE (*, '(3a)') indent, 'History file header:',
     &                  header(27:lenact(header))
c Print out SF History Records onto terminal.
      IF (nhist.NE.0) THEN
         WRITE (*, '(2a)') indent, 'File history:'
         CALL rstrsf(iunit, history, maxhist, nhist, ierrsf)
         DO i = 1, nhist
            WRITE (*, '(2a)') indent, history(i)
         ENDDO
      ENDIF
      GOTO 200

c clear: reset data and plot environment to start values
 255  CONTINUE
      CALL clear(dset, dset_new)
      qcontour = .FALSE.
      GOTO 200

c setplot: command with subcommands to set up plot environment
 300  CONTINUE
      CALL stplot(string, iparse, qini, dset, dset_new)
      GOTO 200

c font: change character font
 325  CONTINUE
      CALL stscf(string, iparse)
      GOTO 200

c symbol: change symbol for data points
 326  CONTINUE
      CALL stsym(string, iparse, dset, dset_new, *200)
      GOTO 200

c label: to write text along axes and on top of plot
 327  CONTINUE
      CALL stlab(string, iparse, mlabel, dset, dset_new)
      GOTO 200

c plot: plot current data
 350  CONTINUE
      qinter = .FALSE.
      GOTO 352

c iplot: plot current data, stay in PLT
 351  CONTINUE
      qinter = .TRUE.
 352  CONTINUE
      IF (qcontour) THEN
         GOTO 631
      ENDIF
      CALL plotter(string, iparse, dset_new, qmodel, qhist, qini, dset,
     &             qinter)
      GOTO 200

c data: read next 'data file' package
 400  CONTINUE
      CALL rddata(string, iparse, iunit)
      GOTO 200

c spectrum: read next 'spectrum' package
 405  CONTINUE
      CALL rdspec(string, iparse, iunit, dset_new, dset, mlabel)
      GOTO 200

c kev_ang: convert keV into angstrom
 410  CONTINUE
      CALL kevang(dset_new, dset)
      GOTO 200

c rdmodel: read next 'photon' package (photon model spectrum)
 415  CONTINUE
      CALL rdmod(iunit, dset_new, dset)
      GOTO 200

c unfoldspec: unfold PHA spectrum using photon model spectrum
 417  CONTINUE
      CALL ufspec(string, iparse, dset, dset_new)
      GOTO 200

c divide: divide two data sets
 420  CONTINUE
      CALL divdat(string, iparse, dset_new, dset)
      GOTO 200

c multiply: multiply two data sets
 421  CONTINUE
      CALL muldat(string, iparse, dset_new, dset)
      GOTO 200

c add: add two data sets
 425  CONTINUE
      CALL addat(string, iparse, dset_new, dset)
      GOTO 200

c residuals: subtract model
 430  CONTINUE
      CALL resid(string, iparse, dset_new, dset)
      GOTO 200

c chi2 residuals: subtract model and divide by error
 431  CONTINUE
      CALL chires(string, iparse, dset_new, dset)
      GOTO 200

c model: plot with or without model
 450  CONTINUE
      CALL pmodel(string, iparse, qmodel, qhist)
      GOTO 200

c logarith: logarithmic axes
 550  CONTINUE
      CALL stlog(string, iparse)
      GOTO 200

c linear: linear axes
 560  CONTINUE
      CALL stlin(string, iparse)
      GOTO 200

c rescale: change scaling on x and y
 600  CONTINUE
      CALL stlim(string, iparse)
      GOTO 200

c rdcont: read steppar package
 630  CONTINUE
      CALL contrd(iunit, chimin, stepsubx, stepsuby, dset, dset_new)
      IF (dset_new.GT.0) qcontour = .TRUE.
      GOTO 200

c plcont: plot chi square contours
 631  CONTINUE
      CALL contpl(string, iparse, chimin, stepsubx, stepsuby, dset,
     &            dset_new, qini)
      GOTO 200

c rdfit: read fit results
 650  CONTINUE
      CALL rdfit(iunit, string, iparse, dset, dset_new)
      IF (dset_new.GT.0) THEN
         CLOSE (iunit)
         CALL opnrsf(filenm, iunit, 'XSPEC histry', header, tmplte,
     &               nhist, ierrsf)
         WRITE (*, '(3a)') indent, 'History file header:',
     &                     header(27:lenact(header))
      ENDIF
      qmodel = .FALSE.
      GOTO 200

c distrib: plot data distribution
 654  CONTINUE
      CALL distri(string, iparse, dset, dset_new)
      GOTO 200

c fold: fold data (time)
 655  CONTINUE
      CALL tfold(string, iparse, dset, dset_new, qmodel)
      GOTO 200

c average: calculate average values
 670  CONTINUE
      CALL averag(string, iparse, dset)
      GOTO 200

c norm: normalize data_set
 750  CONTINUE
      CALL norm(string, iparse, dset, dset_new)
      GOTO 200

c viewport: reset size of viewports
 800  CONTINUE
      CALL newfrm(string, iparse)
      GOTO 200

c newset: rearrange data sets
 850  CONTINUE
      CALL newset(string, iparse, dset, dset_new)
      qmodel = .FALSE.
      GOTO 200

c show: show information
 880  CONTINUE
      CALL xplinf(string, iparse, dset, dset_new, filenm)
      GOTO 200

c hardcopy: make a hardcopy plot
 900  CONTINUE
      CALL ldbuf('@hardcopy.xcm', 1, ierr)
      IF (ierr.NE.0) CALL xwrite('Failed to make hardcopy', 10)
      GOTO 200

c cpd: change the plot device
 910  CONTINUE
      if(iparse.eq.len(string))then
c       no further argument on command
         tmpstr = ' '
         iparse = 0
      else
         tmpstr = string(iparse+1:)
      endif
      string(iparse+1:) = 'device '//tmpstr
      CALL stplot(string, iparse, qini, dset, dset_new)
      GOTO 200

c wrqdp: create qdp file and write data into it
 940  CONTINUE
      CALL mkqdp(string, iparse, dset)
      GOTO 200

c rdqdp: read data from qdp file
 941  CONTINUE
      CALL xrdqdp(string, iparse, dset, dset_new)
      GOTO 200

c exit: exit
 950  CONTINUE
      qpause = .FALSE.
      GOTO 999

c pause: exit with pause
 951  CONTINUE
      qpause = .TRUE.
      GOTO 999

c the end
 999  CONTINUE
      IF (qopen) CLOSE (iunit)
      CALL frelun(iunit)
      IF (qpause) THEN
         PAUSE 'XPLOT: exit with pause'
         GOTO 200
      ENDIF
      WRITE (*, *) ' XPLOT: exit'
      END
