C*********************************************************************
C  DEPAKE : depaking an NMR spectrum. 
C
C	This program is a modified version of the program written
C	by Edward Sternin as part of his M.Sc. thesis project in
C	1982 at the University of British Columbia.  The original
C	version was published in the Appendix to the thesis.  It 
C	was written for use on the Amdahl mainframe running MTS
C	(Michigan Terminal System) operating system.  Installation-
C	specific subroutines were restricted to the user I/O 
C	and the graphics support, and a number of versions
C	of the program exists (often unknown to the author) for
C	various computers.  
C
C	At UBC, the program has been used extensively for a number
C	of years and proved robust and reliable.  Later it  
C	became possible to move the program to the NMR lab computer,
C	DEC MicroVAX.  This version is a "quick and dirty" adaptation of the 
C	original program to the VAX/VMS (and FORTRAN77) environment. It
C	is not a "final" or a "released" version.  In the years since
C	the original was written, the level of sophistication of both the
C	author (he hopes) and of the hardware/support software available
C	to him increased dramatically, yet the limited time available
C	to polish up the program forced a minimalistic approach. 
C	One day...
C
C	The program does, however, reflects a certain philosophy that
C	the author feels strongly about: it is a highly interactive
C	package, guiding the user through the process, one step at a
C	time.  It would be very easy, for example to allow multiple
C	iterations to be performed, but this would result in the user
C	approaching the process of depaking in a "black box" fashion.
C	All my experience with depaking tells me it's the wrong approach.  
C
C	Depaking may appear to be "black magic", and it certainly can  
C	easily lead one astray.  To an attentive user, however, it does
C	give out signals when something is wrong.  Pay attention to these
C	signals!
C***********************************************************************
C
C	References:
C  - M.Bloom, J.H.Davis and A.L.MacKay, Chem.Phys.Lett. 80, p.198 (1981).
C  - E.Sternin, M.Sc. Thesis, Dept. of Physics, UBC (1982).
C  - E.Sternin, M.Bloom, and A.L.MacKay, J.Mag.Res. 55, p.274 (1983).
C
C	Author:		 Edward Sternin
C	Address:	 Department of Physics
C			 Brock University
C			 St.Catharines, Ontario, L2S 3A1 Canada
C	Phone:		 (905) 688-5550 ext. 3414
C	FAX		 (905) 682-9020
C	Electronic mail: edik@BrockU.ca
C
C	Copyright  (c)  Edward Sternin 1982, 1987, 1992, 1998.
C
C	This software is for non-commercial use only.  It is not to be
C	sold, as a whole or in part.  Unlimited duplication is permitted
C	provided this entire header is included. Use on any hardware is 
C	permitted.  
C		    * Modification is not discouraged; *
C	but, please send author a copy of your modifications if you make 
C	any. This way I could provide support, and there would be one
C	standard version to refer to.
C
C	Acknowledgements are appreciated.
C
C	History:
C	v.1	X.82	Thesis version published 
C	v.3	V.84	Last MTS revision, +0.5 bug fix in N0,
C			renormalized data output on unit 4,
C			pretty-plotting (DISSPLA v.8.2).
C	v.4	I.87	Conversion to VAX/VMS, Gplot (from
C			TRIUMF) graphics (slow!!), revised I/O,
C			some device support, allow smoothed output on
C			final plot and on file output. Tested under
C			VMS v.4.1M or later, FORTRAN v.3.5 or later.
C	v.4.1	V.91	small fixes, e.g. handling a header
C	Known bugs:
C  - Dumb TX terminals are badly supported, text may overwrite graphics
C  - Also there, the cursor to the top requires non-blank; a "." used.
C  - EOF during I/O is not trapped, a CTL-Z will crash the program.
C	v.5.0	XII.92  SGI version, complete reliance on GPLOT for graphics
C
C	Current version:	5.0
C	Date:			Dec-1992	
C
C**********************************************************************
C	A remnant from v.1 :
C========================================================================
C        Names of things used by the program
C 
C      NP= # of pts (channels) in the spectrum
C      SW= spectral width, kHz
C      N0= point, about which the first moment is zero
C      BR= Baseline Right
C      BL= Baseline Left
C      I = channel Index, counted from N0
C      P(I) / R(I) are intensities of Powder / Oriented spectra 
C      FAST= size of the averaging bin for the crude primary calculation
C      HEX=1  gives "edge on the left" (HEXagonal) configuration
C      HEX=-1 gives the bilayer configuration.
C
C========================================================================
       IMPLICIT REAL*4 (A-H,O-Z)
       COMPLEX*8 FT(8193)
       REAL*4 P(8193),R(8193),PF(8193),PFOUT(8193)
       REAL*4 Y(5000),X(5000),YINT1(2500),YINT2(2500)
       REAL*4 YOUT(5000),XOUT(5000)
       REAL*4 XSAVE(2500),YSAVE(2500),Y1(2500),Y2(2500)
       REAL*4 XSVOUT(2500),YSVOUT(2500)
       REAL*4 Z(8193),Z1(2),Z2(2),XAX(2),YAX(2)
       INTEGER BL,BR,BLSAVE,BRSAVE,FAST,HEX
       LOGICAL ASKIF,FANCY,LDUM,have_here,have_home
       CHARACTER*80 FNAME,TITLE

	real*4		r_ft(16386)
	integer*4	i_ft(16386)

	equivalence(i_ft,ft)
	equivalence(r_ft,ft)
C
	COMMON SW,N0,BL,BR,NP,NITER
	COMMON/PLTPRM/ XMINSV,XMAXSV,Y1,Y2
C=======================================================================
C	Using GPLOT package for plotting: default to the device types
C	defined in the environment variables (setenv in UNIX)
C		TRIUMF_TERMINAL_TYPE
C		TRIUMF_PLOTTER_TYPE
C	and the default bitmap device = HP LaserJet at 150dpi.
C	If there is a file 'plotdata.ini' in the current directory, 
C	or in the user's home directory, the settings defined there
C	will override the default GPLOT parameter values.

C	call find_unit(lun)		!find an unused i/o unit

  	write(6,'(A,$)') 
     1 	   ' Looking for .gplot.dpk ...'

        inquire(FILE = '.gplot.dpk', EXIST = have_here)
        call getenv('HOME',fname)       !find user's home dir
        n = index(fname,' ')
        fname(n:) = '/.gplot.dpk'
        inquire(FILE = fname(1:n+10), EXIST = have_home)
        if (have_here) then
	  write(6,'(A)') ' found in current directory!'
          call GPLOT_SETUP('.gplot.dpk')
        else if (have_home) then
	  write(6,'(A)') ' found in home directory!'
          call GPLOT_SETUP(fname(1:n)//'/.gplot.dpk')
        else
	  write(6,'(/,A)') 
     1     ' *** Initialization file not found, using GPLOT defaults'
          call GPLOT_SETUP(' ')
        end if


C	fname='gplot_dpk.ini'		!in the current directory
C	open(lun,FILE=fname(1:13),STATUS='OLD',ERR=190)
C	close(lun)
C	call gplot_setup(fname(1:13))
C	goto 192

C  190	call getenv('HOME',fname)	!find user's home dir
C	len=index(fname,' ')
C	fname(len:)='/gplot_dpk.ini'	!look there
C	write(6,'(1X,A,$)') fname(1:len+13)
C	open(lun,FILE=fname(1:len+13),STATUS='OLD',ERR=191)
C	close(lun)
C	write(6,'(A)') '  found!'
C	call gplot_setup(fname(1:len+13))
C	goto 192

C  191		call gplot_setup(' ')

  192   call setnam('NSYINC',2.0)	!a personal preference, default=5.0
	call setnam('NSXINC',2.0)
	call dlinescale(3,11./640.)	!adjust dashes for [Postscript et al.]
	call dlinescale(7,11./640.)	!page size of 11x8.5 in
	call dlinescale(9,11./640.)
	call transparent_mode(1)

	write(6,'(/,A,/,A,/)')
     1	   '               Welcome  to  interactive  De-Pake-ing !',
     2	   '                (c) Ed Sternin, 1982,1987,1992 v.5.0'
C
C========================================================================
C
C	Part 0:	Read in the data. We expect here D.Hare's FTNMR format,
C
  10  WRITE(6,'(A,$)') '  depk-inpt> Data file name : '
      READ(5,FMT='(A)') FNAME
      call find_unit(lun)		!find an unused i/o unit
      OPEN(unit=lun, file=fname, status='OLD', form='UNFORMATTED',
     1	iostat=ios, err=9001)
      sw_in = 0.
  11  read(lun) Npts,(r_FT(I),I=1,abs(2*Npts))
      if(Npts.lt.0) then
C...interpret the header, if any
	    if (i_FT(3).eq.0) then
		write(6,'(a)') '  depk-ERR> need FT data !'
		close(lun)
		goto 10
		end if
	    sw_in = r_FT(17)*0.001	!in kHz
	    NP = i_FT(1)
	    if (i_FT(2).eq.0) then
		jstep = 1		!real data
	    else
		jstep = 2		!complex data
		end if
	    goto 11			!go read another block from the file
      else
	    NP = Npts			!assume complex
	    jstep = 2
	    end if
      close (lun)

C...done: r_FT contains our data, select the real part only
      DO 103 I=1,NP
	P(I)=r_FT( (I-1)*jstep + 1 )
        PF(I)=P(I)
        X(I)=FLOAT(I)
  103 CONTINUE
  151 CONTINUE
      NPSAVE=NP
C====================================================================
C	Part 1:	Parameters.
C	Here we interactively determine a number of settings:
C	 - baseline, i.e. which is the data to be depaked
C	 - binning, i.e. if we want a quick pre-view by reducing the
C	   number of points in the spectrum
C	 - IMPORTANT: determine the precise center of weight of the 
C	   spectrum, a.k.a. zero of the first moment
C
C 	Default the baseline to 0-25% and 75-100% of the full spectrum
C
       BLSAVE=NPSAVE/4
       BRSAVE=NPSAVE*3/4
       	if (sw_in.eq.0.0) then 		!if don't know any better,
	   SWSAVE=NPSAVE		!assume 1 kHz/data point
	   CALL SETNAM('NLXINC',8.0)	!the # of data pts is a multiple of 2
	else
	  SWSAVE = sw_in		!if a known SW, use the normal
          CALL SETNAM('NLXINC',10.0)	!10.-division scale on X axis
	end if
	XMIN=0.0
	XMAX=NPSAVE			!plot against channel # initially

C--------------------------------------------------------------------
C...Do a quick plot here to show the user what he's got 

	IF(.NOT.ASKIF('  depk-inqr> Want to see the data ? @'))
     2	  GOTO 107
       CALL CLEAR_PLOT
        CALL SETLAB('XLABEL','Channel Numbers')
	CALL SETNAM('XMIN',XMIN)
	CALL SETNAM('XMAX',XMAX)
	CALL SETNAM('XAUTO',0.0)
C...to do: change to 10., 2. if DW is read in from the data file
	CALL NARGSI(4)
        CALL GPLOT(X,PF,NPSAVE,1)
       CALL TRANSPARENT_MODE(0)
  107  WRITE(6,'(A,I8,/,A)') 
     1	'  depk-info> No of data points read: ',NPSAVE,
     2	'  depk-WARN> Only the real part will be used.'
C---------------------------------------------------------------------

  101  BL=BLSAVE
       BR=BRSAVE
       NP=NPSAVE
       SW=SWSAVE
       WRITE(6,160)BL,BR,NP,SW
  160  FORMAT(/'  depk-info>                       x'/
     1         '                                  x x'/
     1         '                                 x  x'/
     1         '                             x x    x'/
     1         '                   IxxxxxxIxx        xxxIxxxxxxI'/
     1         '                   0      BL            BR     NP'//
     1 '             BL = ',I6,'   BR = ',I6,'   NP = ',I6,/
     1 '             SW = spectral width = ',G10.4,' kHz'/)
       IF(ASKIF('  depk-inqr> Change these values ? @')) THEN
	  WRITE(6,'(A,$)') '  depk-inpt> Enter  BL,BR,NP,SW : '
          READ(5,*) BLSAVE,BRSAVE,NPSAVE,SWSAVE
	  end if
C---------------------------------------------------------------------
C...To sketch or not to sketch...   Bin averaging.

  105  BL=BLSAVE
       BR=BRSAVE
       NP=NPSAVE
       SW=SWSAVE
       DO 154 I=1,NP
	PF(I)=0
  154  CONTINUE
       FAST=1
       IF(.NOT.ASKIF('  depk-inqr> Want a crude calculation ? @'))
     1	GO TO 104
       WRITE(6,164)
  164  FORMAT(
     1  '  depk-info>',/
     2  '                     *** note ***',/
     3  '     For a quick pre-view we can average over a few ',/
     4  '     data points.The resulting averages are treated as',/
     5  '     single points.This speeds up the calculations but',/
     6  '     should not be relied upon. Set the averaging',/
     7	'     parameter so that the total number of points is a',/
     8  '     multiple of it.',/)
	write(6,'(A,$)') 
     1	   '  depk-inpt> Number of points to be averaged into one bin : '

        READ(5,*)FAST

  104   CONTINUE
       L=0
       PERPT=SW/FLOAT(NP)
       PMAX=0
       IFIN=NP-FAST+1
       DO 111 I=1,IFIN,FAST
        L=L+1
        JFIN=I+FAST-1
        DO 131 J=I,JFIN
        PF(L)=PF(L)+P(J)
  131   CONTINUE
        PF(L)=PF(L)/FAST
        IF(PMAX.LT.PF(L))PMAX=PF(L) 
        X(L)=0.5*FLOAT(2*I+FAST)*PERPT 
  111  CONTINUE
       NP=NP/FAST
       BL=BL/FAST
       BR=BR/FAST
       PERPT=SW/FLOAT(NP)
C---------------------------------------------------------------------
C...Baseline correction

       B1=0
       DO 112 I=1,BL
         B1=B1+PF(I)
  112  CONTINUE
       DO 113 I=BR,NP
       B1=B1+PF(I)
  113  CONTINUE
       B1=B1/(NP-BR+BL+1)
       DO 114 I=1,NP
         PF(I)=PF(I)-B1
  114  CONTINUE

C...Determine zero of the first moment
       AR1=0
       AR2=0
       DO 115 I=BL,BR
         AR1=AR1+PF(I)
         AR2=AR2+I*PF(I)
  115  CONTINUE
       N0=AINT(AR2/AR1+.5)

C...Primary plotting
       Z2(1)=FLOAT(N0)*PERPT
       Z2(2)=Z2(1)
       Z1(1)=-0.1*PMAX
       Z1(2)=1.05*PMAX
       NBIG=MAX0(N0-BL,BR-N0)
       XMIN=SW*FLOAT(N0-NBIG)/FLOAT(NP)
       XMAX=SW*FLOAT(N0+NBIG)/FLOAT(NP)
       XMINSV=XMIN
       XMAXSV=XMAX
       IF(ASKIF('  depk-inqr> Plot what is to be depaked ? @')) THEN
       	   CALL CLEAR_PLOT
           CALL SETLAB('XLABEL','Frequency (kHz)')
           CALL SETLAB('YLABEL','Intensity (arb. units)')
	   CALL SETNAM('XAUTO',0.0)
           CALL SETNAM('XMIN',XMIN)
           CALL SETNAM('XMAX',XMAX)
           CALL SETNAM('YMIN',-0.1*PMAX)
           CALL SETNAM('YMAX',1.05*PMAX)
	   CALL NARGSI(4)
           CALL GPLOT(X,PF,L,1)
           CALL SETNAM('LINTYP',3.)
	   CALL NARGSI(4)
	   CALL NARGSI(4)
           CALL GPLOT(Z2,Z1,2,2)
           CALL TRANSPARENT_MODE(0)
	   LDUM=ASKIF('  depk-paus> Hit <CR> to continue : @')
	   end if

C...Calculation of the moments
       WM0=0.0
       WM1=PF(N0)
       WM2=0.0
       WM3=0.0
       WM4=0.0
       KM=BR-N0
       IF(KM.LT.(N0-BL))KM=N0-BL
       DO 141 K=1,KM
         WM1MIN=PF(N0-K)-PF(N0+K)
         WM1PLU=PF(N0-K)+PF(N0+K)
C...! do not use integer arithmetics, it saves time, but may overflow
         WM1=WM1+FLOAT(K)*WM1MIN
         WM2=WM2+FLOAT(K)**2*WM1PLU
         WM3=WM3+FLOAT(K)**3*WM1MIN
         WM4=WM4+FLOAT(K)**4*WM1PLU
         WM0=WM0+WM1PLU
         Y(K)=WM1/WM0
         X(K)=FLOAT(K)
  141  CONTINUE
C
C What we have done here is to calculate the zero of the first
C moment from the given baseline limits, and then re-calculated
C the first moment about that point extending the area of integration
C from +/- 0 to +/- edges of the baseline.  While the limits of
C integration are within the spectrum, the result not =0 
C (WARNING: unless the spectrum has been artificially symmetrized by
C zeroing the out of phase channel), but as the limits extend into
C the baseline, the result should tail off at zero, if the zero
C of the first moment has been determined correctly.  If there is no
C tail-off, then either the spectrum extends beyond the given
C baseline limits, i.e. we will miss some systematic intensity if
C we depake now, or there is a tilt in the baseline.  If, on the 
C other hand, the tail is too long, we are depaking noise and wasting
C CPU time.
C--------------------------------------------------------------------
       CALL CLEAR_PLOT
       CALL SETNAM('XAUTO',2.)
       CALL SETNAM('YAUTO',2.)
       CALL SETLAB('XLABEL','DISTANCE from N<_>0')
       CALL SETLAB('YLABEL','FIRST MOMENT minus N<_>0')
       CALL SETNAM('LINTYP',1.)
       CALL NARGSI(4)
       CALL GPLOT(X,Y,KM,1)
       CALL TRANSPARENT_MODE(0)  

C...Normalize the moments
       TWOPI=4.0*ASIN(1.0)
       WM2=WM2/WM0*TWOPI**2
       WM3=WM3/WM0*TWOPI**3
       WM4=WM4/WM0*TWOPI**4
       WRITE(6,168)WM2,WM3,WM4
  168  FORMAT('  depk-info> Moments:  M2 = ',G16.8,/
     1        '                       M3 = ',G16.8,/
     2        '                       M4 = ',G16.8,/)
       WRITE(6,166)BL,N0,BR,NP,SW
  166  FORMAT('  depk-info> BL = ',I4,' N0 = ',I4,' BR = ',I4,
     1 ' NP = ',I4,' SW = ',G10.4,/)

C...This is to be used by an experienced depaker only !!
       WRITE(6,'(A)')
     2	'  depk-WARN> To be used by an experienced depaker only !'
       IF(ASKIF('  depk-inqr> Want to modify N0 ? @')) THEN
          WRITE(6,'(A,$)') '  depk-inpt> Enter new N0 : '
          READ(5,*)N0
	  end if
  169  IF(.NOT.ASKIF('  depk-inqr> OK to depake now ? @')) GO TO 101 
C=======================================================================
C 
C      PART 2 : De-Pake-ing
C
C     First time around we know that the contribution
C     from the other side is zero ( since we haven't
C     been there yet ).
C
 200  NMAX=2*AINT(MAX0(BR-N0,N0-BL)/2.)
C         =even #

 210  HEX=1
      IF(ASKIF('  depk-inqr> Expecting an edge on the right ? @'))HEX=-1
      IPLOT=0
      YMIN=1.E10
      YMAX=-1.E10
      DO 201 I=1,NP
 201  R(I)=0.
C
      DO 213 I=1,NMAX
      N=NMAX-I+1
      NN=2*N
      NSGND=HEX*N
      SUM=0.
      IF(I.EQ.1)GO TO 218
      N21=NN+2
      N22=2*NMAX
      JFIN=N22-N21+1
      IF(HEX.EQ.1)GO TO 217
      N21=-2*NMAX
 217  DO 212 J=1,JFIN,2
      M=N21-1+J
 212  SUM=SUM+R(N0+M)*RAT1(M,-NSGND)
 218  NOR=N0+2*NSGND
      ORIENT=(1.732050808*PF(N0-NSGND)-SUM)/RAT2(-NSGND)
      R(NOR)=ORIENT
C
      IPLOT=IPLOT+1
      Y(IPLOT)=ORIENT
      X(IPLOT)=PERPT*NOR
C
 213  CONTINUE
C
      NITER=0
      NS=0
      NPSAV=0
 220  CALL DEPLOT(X,Y,IPLOT,NS)
      NITER=NITER+1
 221  WRITE(6,'(A,I3,A,$)') '  depk-inqr>',NITER,' iterations done,'
      IF(ASKIF('  enough ? @')) GO TO 300
      NPSAV=IPLOT
      DO 224 I=1,IPLOT
       XSAVE(I)=X(I)
       YSAVE(I)=Y(I)
 224  CONTINUE

C...Iterative de-Pake-ing
      IPLOT=0
      IHEX=HEX*(-1)**NITER
      DO 231 I=1,NMAX
      N=NMAX+1-I
      NSGND=IHEX*N
      SUM1=0.
      SUM2=0.
      IF(I.EQ.1)GO TO 238
      N31=2*AINT(N/2.)+2
      N32=-2*NMAX
      N33=-2*N-2
      IF(IHEX.EQ.-1)GO TO 237
      N31=2*N+2
      N33=-N-1
C              note that the last to contribute is (-N-2); +(-n)
 237  J1FIN=2*NMAX-N31+1
      J2FIN=N33-N32+1
      DO 232 J=1,J1FIN,2
      M=N31-1+J
 232  SUM1=SUM1+R(N0+M)*RAT1(M,-NSGND)
      DO 233 J=1,J2FIN,2
      M=N32-1+J
 233  SUM2=SUM2+R(N0+M)*RAT1(M,-NSGND)
 238  NOR=N0+2*NSGND
      ORIENT=1.732050808*PF(N0-NSGND)-SUM1-SUM2-R(N0-NSGND)*RAT3(-NSGND)
      ORIENT=ORIENT/RAT2(NSGND)
      R(NOR)=ORIENT
C
      IPLOT=IPLOT+1
      Y(IPLOT)=ORIENT
      X(IPLOT)=PERPT*NOR
 231  CONTINUE
C 
      GO TO 220
C
C  End of de-Pake-ing
C  At this stage, X(IPLOT),Y(IPLOT) and XSAVE(NPSAVE),YSAVE(NPSAVE)
C  contain the results of the (Niter)'th and (Niter-1)'th 
C  iterations, respectively 
C
C===================================================================
C
C
C      PART 3 : The time of judgement 
C
  300 FANCY=.FALSE.
      CALL TRANSPARENT_MODE(1)
      WRITE(6,309)
  309 FORMAT('  depk-info> The following options are available :'/
     1       '             (1)  start from scratch'/
     1       '             (2)  change the averaging bin size'/
     1       '             (3)  start dePakeing from the other side'/
     1       '             (4)  continue dePakeing (more iterations)'/
     1       '             (5)  fancy plot of the results'/
     1       '             (6)  save output in a file'/
     1       '             (7)  quit')
      WRITE(6,'(A,$)') '  depk-inpt> Enter your choice : '
      READ(5,*)NCHOIC
      GOTO (101,105,210,221,310,320,600),NCHOIC
      WRITE(6,308)
  308 FORMAT('  depk-ERR-> Not understood')
      GOTO 300

C...Set this flag if the plot is wanted
  310 FANCY=.TRUE.

C...Calculate the integral starting from both ends towards
C   the mid-point ( = N0 = center-of-weight ).
C
  320 CONTINUE
      ORMAX=Y(1)
      ORMIN=ORMAX
      YINT1(1)=ORMAX
      YIMAX=ORMAX
      DO 321 I=2,IPLOT
      YDUM=Y(I)
      IF(YDUM.GT.ORMAX)ORMAX=YDUM
      IF(YDUM.LT.ORMIN)ORMIN=YDUM
      YINT1(I)=YINT1(I-1)+YDUM
      IF(YINT1(I).GT.YIMAX)YIMAX=YINT1(I)
  321 CONTINUE
C
      YINT2(1)=YSAVE(1)
      DO 322 I=2,NPSAV
      YDUM=YSAVE(I)
      IF(YDUM.GT.ORMAX)ORMAX=YDUM
      IF(YDUM.LT.ORMIN)ORMIN=YDUM
      YINT2(I)=YINT2(I-1)+YDUM
      IF(YINT2(I).GT.YIMAX)YIMAX=YINT2(I)
  322 CONTINUE

C...Normalize all curves to unit max intensity 
      SW2=0.5*SW
      IF(ASKIF(
     1	'  depk-inqr> Use filtered data on the final plot ? @')) THEN
	  DO 341 I=1,IPLOT
	    XOUT(I)=X(I)-SW2
	    YOUT(I)=Y1(I)/ORMAX
  341       YINT1(I)=YINT1(I)/YIMAX
	  DO 342 I=1,NPSAV
	    XSVOUT(I)=XSAVE(I)-SW2
	    YSVOUT(I)=Y2(I)/ORMAX
  342 	    YINT2(I)=YINT2(I)/YIMAX
        ELSE
	  DO 344 I=1,IPLOT
	    XOUT(I)=X(I)-SW2
	    YOUT(I)=Y(I)/ORMAX
  344       YINT1(I)=YINT1(I)/YIMAX
	  DO 345 I=1,NPSAV
	    XSVOUT(I)=XSAVE(I)-SW2
	    YSVOUT(I)=YSAVE(I)/ORMAX
  345 	    YINT2(I)=YINT2(I)/YIMAX
       ENDIF	  
       ORMIN=ORMIN/ORMAX
       IF(-.20.LT.ORMIN)ORMIN=-.20
       DO 343 I=1,NP
         Z(I)=FLOAT(I)*PERPT-SW2
  343    PFOUT(I)=PF(I)/PMAX
C
      IF(.NOT.FANCY) GOTO 500
C==================================================================
C	Do some fancy plotting -- access GPLOT in an interactive 
C	mode, allow to change parameters and produce hardcopy, etc.
C
      CALL CLEAR_PLOT
      CALL SETNAM('XAUTO',0.)
      CALL SETNAM('YAUTO',0.)
      CALL SETNAM('XMIN',XMINSV-SW2)
      CALL SETNAM('XMAX',XMAXSV-SW2)
      CALL SETNAM('YMIN',ORMIN)
      CALL SETNAM('YMAX',1.0)
      CALL SETLAB('XLABEL','Frequency (kHz)')
      CALL PLOT_COLOR(7,1)		!start off in white
      CALL SETNAM('CURSOR',2.)		!interactive text, abut to LC

 465  call clear_plot
      CALL SETNAM('LINTYP',1.)
      CALL NARGSI(4)
      CALL GPLOT(XOUT,YOUT,IPLOT,4)     !draw axes also
C Make an artificial x,y cross in the middle
      XAX(1)=0.0
      XAX(2)=0.0
      YAX(1)=ORMIN
      YAX(2)=1.0
      CALL SETNAM('LINTYP',3.)		!dashed line
      CALL NARGSI(4)
      CALL GPLOT(XAX,YAX,2,2)
      XAX(1)=-SW2
      XAX(2)=SW2
      YAX(1)=0.0
      YAX(2)=0.0
      CALL NARGSI(4)
      CALL GPLOT(XAX,YAX,2,2)
      CALL PLOT_COLOR(1,2)   		!Make depaked results red
      CALL SETNAM('LINTYP',1.)
      CALL NARGSI(4)
      CALL GPLOT(XOUT,YOUT,IPLOT,2)
      CALL NARGSI(4)
      CALL GPLOT(XSVOUT,YSVOUT,NPSAV,2)
      CALL PLOT_COLOR(6,1)  		!Make powder blue (cyan)
      CALL SETNAM('LINTYP',9.)		! and dotted
      CALL NARGSI(4)
      CALL GPLOT(Z,PFOUT,NP,2)
      CALL PLOT_COLOR(4,1) 		!Make integrals  green
      CALL SETNAM('LINTYP',7.)		! and short-dashed
      CALL NARGSI(4)
      CALL GPLOT(XOUT,YINT1,IPLOT,2)
      CALL NARGSI(4)
      CALL GPLOT(XSVOUT,YINT2,NPSAV,2)
      CALL PLOT_COLOR(7,1)		!restore white
      call transparent_mode(0)
      call gplot_control(' depk-gplot> ',&465)

  409 call gploti			!next time start with defaults

C     IF(ASKIF('  depk-inqr> Do you want a hardcopy ? @')) THEN
C	CALL NARGSI(1)
C	CALL GRAPHICS_HARDCOPY(0)
C     END IF
C
      IF(.NOT.ASKIF('  depk-inqr> List the results to a file ? @'))
     2	GOTO 600

C...Output Xori, Yori, Yintegral, and Yori-smoothed
C   Do not bother sorting the data in x, plotting program can do it.
  500 WRITE(6,FMT='(A,$)')'  depk-inpt> File name : '
      READ(5,FMT='(A)') FNAME
      call find_unit(lun)		!find an unused i/o unit      
      OPEN(UNIT=lun,FILE=FNAME,STATUS='NEW',IOSTAT=IOS,
     1	CARRIAGECONTROL='LIST',ERR=9002)
C...always write out the unfiltered data
      WRITE(lun,508)(XOUT(I),Y(I),YINT1(I),I=1,IPLOT)
      WRITE(lun,508)(XSVOUT(I),YSAVE(I),YINT2(I),I=1,NPSAV)
      CLOSE(lun)
  508 FORMAT(3G16.8)

  600 IF(.NOT.ASKIF('  depk-inqr> Last chance: ready to quit ? @'))
     1	GO TO 300

       CALL TRANSPARENT_MODE(1)
       WRITE(6,'(A)') '  depk-stop> Well... See ya.'
       STOP

C...Error handling
 9001  WRITE(6,9901) IOS
 9901  FORMAT('  depk-ERR-> Input data file error, I/O status = ',I8)
       GO TO 10
C
 9002  WRITE(6,9902) IOS
 9902  FORMAT('  depk-ERR-> Listing file error, I/O status = ',I8)
       GO TO 500

       END
C==================================================================
C
C        Functions, subroutines etc.
C
C==================================================================
C
      FUNCTION RAT1(M,N)
      MP=M+1
      MM=M-1
      RAT=SQRT(MP*(MP+2.*N))+MP+N
      RAT=RAT/(SQRT(MM*(MM+2.*N))+MM+N)
      RAT1=ALOG(ABS(RAT))
      RETURN
      END
C
      FUNCTION RAT2(N)
      RAT=SQRT(2.*IABS(N)+1.)+N+1.
      IF(N)421,422,422
 422  RAT2=ALOG(ABS(RAT/N))
      RETURN
 421  RAT2=ALOG(ABS(N/(RAT-2)))
      RETURN
      END
C
      FUNCTION RAT3(N)
      NABS=IABS(N)
      RAT=SQRT((NABS+1.)*(3.*NABS+1.))+2.*N+1.
      IF(N)431,432,432
 432  RAT3=ALOG(ABS(RAT/3.732050808/NABS))
      RETURN
 431  RAT3=ALOG(ABS(-.267949192*NABS/(RAT-2.)))
      RETURN
      END
C
C====================================================================
C
      SUBROUTINE DEPLOT(X,Y,IPLOT,NS)
      COMMON SW,N0,BL,BR,NP,NITER
      COMMON/PLTPRM/ XMINSV,XMAXSV,Y1,Y2
      LOGICAL ASKIF
      INTEGER ATT,BINO,BL,BR
      CHARACTER*2 NUMB
      REAL*4 YMINSV/1.0/,YMAXSV/0.0/
      REAL*4 XMINSV,XMAXSV
      DIMENSION X(2500),Y(2500),Y1(2500),Y2(2500)
      DIMENSION X2(2),X3(2),Z2(2),Z3(2)
      DIMENSION ATT(36),NATT(6)
      DIMENSION XPSAVE(2500),YPSAVE(2500)
C
      DATA ATT/1,0,0,0,0,0,2,1,0,0,0,0,6,4,1,0,0,0,20,15,6,1,0,0,
     1 70,56,28,8,1,0,252,210,120,45,10,1/
      DATA NATT/1,4,16,64,256,1024/
C
      YMIN=YMINSV
      YMAX=YMAXSV
      XMIN=XMINSV
      XMAX=XMAXSV
      DO 1011 I=1,IPLOT
      ORIENT=Y(I)
      IF(YMIN.GT.ORIENT) YMIN=ORIENT
      IF(YMAX.LT.ORIENT) YMAX=ORIENT
1011  Y1(I)=Y(I)
      DO 1012 I=1,NPSAV
      ORIENT=YPSAVE(I)
      IF(YMIN.GT.ORIENT) YMIN=ORIENT
      IF(YMAX.LT.ORIENT) YMAX=ORIENT
      XDUM=XPSAVE(I)
      IF(XMINSV.GT.XDUM) XMINSV=XDUM
      IF(XMAXSV.LT.XDUM) XMAXSV=XDUM
 1012 Y2(I)=YPSAVE(I)
      IF(YMIN.GT.-0.2*YMAX)YMIN=-0.2*YMAX
      X3(1)=XMIN
      X3(2)=XMAX
      Z3(1)=0.0
      Z3(2)=0.0
      X2(1)=FLOAT(N0)*SW/FLOAT(NP)
      X2(2)=X2(1)
      Z2(1)=YMIN
      Z2(2)=1.05*YMAX
 1002 CONTINUE
C
C	Plot out the current and the previous iteration
C
      CALL CLEAR_PLOT
      CALL SETNAM('XAUTO',0.)
      CALL SETNAM('YAUTO',0.)
      CALL SETNAM('PMODE',0.)
      CALL SETNAM('XMIN',XMIN)
      CALL SETNAM('XMAX',XMAX)
      CALL SETNAM('YMIN',YMIN)
      CALL SETNAM('YMAX',1.05*YMAX)
C
      CALL SETNAM('LINTYP',1.)  	!Solid line
      CALL PLOT_COLOR(1,2)		!Depaked data in red
      CALL NARGSI(4)
      CALL GPLOT(X,Y1,IPLOT,2)  	!Just data line
      CALL PLOT_COLOR(5,1)  !add iteration count in yellow, black
      YDUM=85.
      IF(X(2).LT.X(1)) then
           XDUM=55.
      	   CALL SETNAM('CURSOR',-1.)
      else
           XDUM=45.
      	   CALL SETNAM('CURSOR',-3.)
	   end if
      CALL SETNAM('%XLOC',XDUM)
      CALL SETNAM('%YLOC',YDUM)
      WRITE(NUMB,'(I2)')NITER+1		!encode # iterations in a string
      CALL SETLAB('TEXT','Iteration No '//NUMB)
      CALL PLOT_COLOR(1,2)		!prev iteration in red
      IF(NITER.GE.1) THEN
	   CALL NARGSI(4)
	   CALL GPLOT(XPSAVE,Y2,NPSAV,2)
	   END IF
      CALL SETNAM('LINTYP',3.)  	!Dotted line
      CALL PLOT_COLOR(7,1)		!Axes in white,black
      CALL NARGSI(4)
      CALL GPLOT(X2,Z2,2,2)
      CALL NARGSI(4)
      CALL GPLOT(X3,Z3,2,2)
      CALL SETLAB('XLABEL','Frequency (kHz)')
      CALL SETLAB('YLABEL','Intensity (arb. units)')
      CALL NARGSI(4)
      CALL GPLOT(X,Y1,IPLOT,4)  	!Just axes
      CALL TRANSPARENT_MODE(0)
C
      IF(.NOT.ASKIF('  depk-inqr> Want to filter the output ? @'))
     1	GO TO 1006
      WRITE(6,1009)
 1009 FORMAT('  depk-inpt> Enter correlation +/- 1/2 length',/
     2       '             for binomial filtering, (max 5) : ',$)
C
      READ(5,*)BINO
      JFIN=BINO*2+1
      IBEG=BINO+1
      IFIN=IPLOT-BINO
       DO 1004 I=IBEG,IFIN
       Y1(I)=0
       DO 1041 J=1,JFIN
       N=I-BINO+J-1
       IDUM=6*BINO+IABS(BINO-J+1)+1
       Y1(I)=Y1(I)+Y(N)*ATT(IDUM)
 1041  CONTINUE
 1004  Y1(I)=Y1(I)/NATT(BINO+1)
C
      IF(NITER.LT.1)GO TO 1002
      IFIN=NPSAV-BINO
       DO 1005 I=IBEG,IFIN
       Y2(I)=0
       DO 1051 J=1,JFIN
       N=I-BINO+J-1
       IDUM=6*BINO+IABS(BINO-J+1)+1
       Y2(I)=Y2(I)+YPSAVE(N)*ATT(IDUM)
 1051  CONTINUE
 1005  Y2(I)=Y2(I)/NATT(BINO+1)
      GO TO 1002
C
 1006 DO 1061 I=1,IPLOT
      XPSAVE(I)=X(I)
      YPSAVE(I)=Y(I)
 1061 CONTINUE
      NPSAV=IPLOT
      YMINSV=YMIN
      YMAXSV=YMAX
      RETURN
      END
C=======================================================================
C
C...ASKIF -- a short Logical Function to get .TRUE. or .FALSE. result as 
C            a Y or N answer to the Character QUESTION, terminated by '@'    
C
	logical function askif(string) 
	character*80	string
        character*1	char/' '/,eol/'@'/
C--------------------------------------------------------------------<     
	nch = index(string,eol)
	if(nch.gt.0) write(6,'(A,$)') string(1:nch-1)

 1000   read(5,'(A1)') char
        
	if (char .eq. 'y' .or. char .eq. 'Y') then
	     askif = .true.  
	     return 
	else if (char.eq.' ' .or. char.eq.'n' .or. char.eq.'N') then
	     askif = .false. 
 	     return  
	else
	     write (6,'(A,$)') 'Please answer Y or N: '	
	end if
     
	go to 1000	
	end