c*************
c  Description:  This module contains primitive computational support routines.
c
c  Author:       JGS   4/2/90
c
c  MIB 12-Apil-2000                                            DR RRcl00788
c...The following use of XX and YY are a work around to  a complier problem
c for a new ES40 with VMS 7.1-2 and latest fortran complier EV7.2-180-4293Q
c per jgk: 16-Jul-1999
c apparently this compiler did not like (i+1:j-1) but (XX:YY) is ok MIB
c
c  22-Mar-2023 zwm Fixed lxor = '80808080'X to INT(Z'80808080')
c
c*************


	function M$WINDOW (buf,name,nout)
c
	include 'types.inc'
	include 'headers.inc'
	record /HEADER/ hcb
	character*(*) name
	real*4 buf(*)
	real*8 delta
c
c catch symmetric flag
	npts = nout
	if (nout.lt.0) npts = -nout - 1
c
c get window ID
	id = M$WINDOW_ID (name)
c
c standard window name
	if (id.gt.0) then
	  call M$WINDOW_BY_ID (buf,id,npts)
c
c non-standard window
	else
	  hcb.file_name = name
	  hcb.type_allow = '1000'
	  hcb.format_allow = 'SF'
	  hcb.switch = 'NOABORT'
	  call M$INOK (hcb)
	  if (.not.hcb.open) then
	    call M$ERROR('Couldnt open non-standard window file in M$WINDOW')
	  elseif (hcb.size .eq. npts) then
	    call M$GRAB (hcb,buf,1.d0,npts,-1)
	  else
	    call M$WARNING ('Interpolating window file in M$WINDOW')
	    nwind = hcb.size
	    delta = float(nwind)/float(npts)
	    call M$MALLOC (nwind*4+4,loc)
	    call M$GRAB (hcb,%val(loc),1.d0,nwind,-1)
	    call M$VMOV (%val(loc),1,%val(loc+nwind*4),1,1)
	    call M$INTERP (nwind+1,0.d0,1.d0,%val(loc),
     &			   npts,0.d0,delta,buf,1,3,ne)
	    call M$MFREE (nwind*4+4,loc)
	  endif
	endif
c
c add symmetric end point
	if (nout.lt.0) buf(npts+1) = buf(1)
c
	M$WINDOW = id
	return
	end



	function M$WINDOW_ID (wind)
	implicit none
	include 'constants.inc'
	include 'types.inc'
	character*4 wind, test
	integer*4   i, j
c
	i = 1
	do while (M$PARSE(standard_windows, i, test) .gt. 0)
	  if (test .eq. wind) then
	    M$WINDOW_ID = i
	    return
	  endif
	  i = i + 1
	enddo
	M$WINDOW_ID = 0
	return
	end



	subroutine M$WINDOW_BY_ID (buf,id,npts)
c
	include 'constants.inc'
	real*4 buf(*),data(4),c
c
c init variables for trig windows
	do n = 1,4
	  data(n) = 0.0
	enddo
c
c search for standard named window
	if (id.eq.2) then	! BART
	  nptsh = npts/2+1
	  do n = 1,nptsh
	    buf(n) = 2*float(n-1)/float(npts)
	  enddo
	  do n = nptsh+1, npts
	    buf(n) = 2*float(npts-n+1)/float(npts)
	  enddo
	  return
	elseif (id.eq.3) then	! HANN
	  data(1) = 0.5
	  data(2) = -.5
	elseif (id.eq.4) then	! HAMM
	  data(1) = 0.54
	  data(2) = -.46
	elseif (id.eq.5) then	! BLACK
	  data(1) = 0.42
	  data(2) = -.5
	  data(3) = 0.08
	elseif (id.eq.6) then	! BH61
	  data(1) = 0.44959
	  data(2) = -.49364
	  data(3) = 0.05677
	elseif (id.eq.7) then	! BH67
	  data(1) = 0.42323
	  data(2) = -.49755
	  data(3) = 0.07922
	elseif (id.eq.8) then	! BH74
	  data(1) = 0.40217
	  data(2) = -.49703
	  data(3) = 0.09897
	  data(4) = -.00183
	elseif (id.eq.9) then	! BH92
	  data(1) = 0.35875
	  data(2) = -.48829
	  data(3) = 0.14128
	  data(4) = -.01168
	else
	  if (id.ne.1) call M$WARNING
     &	  ('Unknown window ID in M$WINDOW_BY_ID - using rectangular')
	  do n = 1,npts
	    buf(n) = 1.0
	  enddo
	  return
	endif
c
c finish trigonometric windows
	if (M$GENRADIX(npts,m2,m3,m4,m5).eq.1) then	! use FFT
	  do n = 1,npts
	    buf(n) = 0.0
	  enddo
	  buf(1) = data(1)
	  do n = 2,4
	    buf(2*n-1) = data(n)/2.
	  enddo
	  call M$RFFT( buf, npts, -2)	! inverse packed FFT

	else						! use trig
	  ndo = npts
	  if (mod(npts,2).eq.0) ndo = npts/2+1	! even # points
	  do n = 1,ndo
	    c = twopi*(n-1)/float(npts)
	    buf(n) = data(4)*cos(3*c) + data(3)*cos(2*c) +
     &		 data(2)*cos(c) + data(1)
	  enddo
	  do n = ndo+1,npts
	    buf(n) = buf(npts-n+2)
	  enddo
	endif
c
	return
	end



c=======================================================================
c   SPLINE/LINEAR/HOLD INTERPOLATION ROUTINES                     JGS/891006
c
c      Interpolation from file in single precision
c      SUBROUTINE FINTERP (HCB, TVIN,
c                            NOUT, TSTRO, TDELO, TVOUT, NCOMP, NORDER, NERR)
c
c      Interpolation from file in double precision
c      SUBROUTINE FINTERPD (HCB, TVIN,
c                            NOUT, TSTRO, TDELO, TVOUT, NCOMP, NORDER, NERR)
c
c      Interpolation from single precision table
c      SUBROUTINE INTERP (NIN, TSTRI, TDELI, TVIN,
c                            NOUT, TSTRO, TDELO, TVOUT, NCOMP, NORDER, NERR)
c
c      Interpolation from double precision table
c      SUBROUTINE INTERPD (NIN, TSTRI, TDELI, TVIN,
c                             NOUT, TSTRO, TDELO, TVOUT, NCOMP, NORDER, NERR)
c
c  These subroutines perform multiple, hold, linear, or cubic-spline
c  interpolations on a vector in a tabulated table of vectors.   It
c  does not perform extrapolation.  If extrapolation is requested, the
c  output vector is equated to that at the nearest edge of the table.
c  Execution is not halted.
c
c  If interpolation immediately proximate to the beginning or end of
c  the table is requested, linear (vice cubic-spline interpolation) is
c  used and the error flag is set to note that fact.
c
c
c  INPUT PARAMETERS:
c  -----------------
c HCB    is the header control block of the file containing
c        vectors forming the table to be interpolated
c
c  NIN   is the number of vectors entered in the table to be
c        interpolated.  Since linear interpolation is used at the
c        extreme ends of the table, for maximum accuracy the table
c        would include at least two entries smaller than the smallest
c        abscissa for which interpolation is requested and two values
c        greater than the greatest abscissa for which interpolation is
c        requested.  If maximum accuracy is not achieved in all the
c        requested  interpolations, the error flag is set to indicate
c        that fact.
c
c TSTRI  is the start abscissa value corresponding to the table of
c        vectors to be interpolated.
c
c TDELI  is the abscissa delta value corresponding to the table of
c        vectors to be interpolated.
c
c  TVIN  is the table of tabulated vectors for which interpolation
c        will be done.  The table is arranged so that all the components
c        of each vector are stored contiguously in memory.
c
c  NOUT  is the number of abscissa values for which interpolation is
c        requested.
c
c TSTRO  is the start abscissa value corresponding to the table of
c        vectors to be output.
c
c TDELO  is the abscissa delta value corresponding to the table of
c        vectors to be output.
c
c NCOMP  is the length of the vector (i.e., the number of components).
c
c NORDER the highest order of interpolation to be used
c        (0-Hold, 1-Linear, or 3-Spline)
c
c
c  OUTPUT PARAMETERS:
c--------------------
c  TVOUT is the output table containing the NOUT interpolated vectors.
c        The table is arranged so that all the components of each
c        vector are stored contiguously in memory.  The i-th row of
c        the table contains the i-th components of the vectors;  the
c        j-th columns contains the j-th vector, corresponding to the
c        abscissa value TSTRO + TDELO*(j-1).  Thus the  dimensions of
c        this table are TVOUT(NCOMP, NOUT).
c
c  NERR  is the error flag.  If more than a single condition is
c        encountered during the execution of the subroutine, only
c        the last non-zero value is reported by the flag.
c        Its interpretation is as follows:
c
c     NERR = 0  => No error was encountered; maximum interpolation
c                  accuracy was achieved for all requested
c                  interpolations.
c
c     NERR = 1  => Linear interpolation (vice cubic-spline) was
c                  used for one or more of the requested
c                  interpolations.
c
c     NERR = 2  => Extrapolation was requested one or more times.
c                  That is, the value for which interpolation was
c                  requested fell outside the tabulated table.  In
c                  this case, the interpolated vector was set equal
c                  to that at the nearest edge of the table.
c
c=======================================================================
c  REVISIONS:
c  891006  -  Originated
c=======================================================================

      subroutine M$FINTERP (HCB,
     +             NOUT, TSTRO, TDELO, TVOUT, NCOMP, NORDER, NERR)
c
      include 'headers.inc'
c
      record /HEADER/ HCB
      SAVE
      REAL*8      TVOUT(*), TSTRO, TDELO
      INTEGER*4   NCOMP,NORDER,NERR,NIN,NOUT
      character*2 format
      REAL*8      X_MIN,X_MAX,X_RAN,X_START,X_DELTA
c
c get necessary file specs and error check
      if (HCB.format.eq.'NH') call M$ERROR
     & ('M$FINTERP does not support non-homogeneous files')
      format = 'S'//HCB.format(2:2)
      NBPA  = M$BPA(format)
      NCOMP = HCB.bpe/NBPA
      HCB.buf_type = 'F'
c
c find range of output assuming monotonically increasing
      X_MIN = TSTRO
      X_MAX = TSTRO + (NOUT-1)*TDELO
      X_RAN = X_MAX-X_MIN
c
c set proper abscissa fields
      if (HCB.class .eq. 2) then
	X_START = HCB.ystart
	X_DELTA = HCB.ydelta
      else
	X_START = HCB.xstart
	X_DELTA = HCB.xdelta
      endif
c
c calculate read indices for grab and get data
      HCB.offset = dnint( (X_MIN-X_START)/X_DELTA - 2 )
      HCB.offset = max( 0.d0, min( HCB.size-1.d0, HCB.offset) )
      NGET = dnint( (X_MAX-X_START)/X_DELTA+3 ) - HCB.offset
      NGET = min( NGET, INT(HCB.size-HCB.offset) )
      X_START = X_START + HCB.offset*X_DELTA
      NTOT = NGET*NCOMP
      nbytes = max(4,NBPA) * NTOT
      call M$MALLOC ( nbytes, loc_TVIN )
      call M$GRAB (HCB, %val(loc_TVIN), HCB.offset+1.d0, NGET, -1)
c
c call interp routine
      call M$INTERP ( NGET, X_START, X_DELTA, %val(loc_TVIN),
     +             NOUT, TSTRO, TDELO, TVOUT, NCOMP, NORDER, NERR)
c
c free memory and report errors
      call M$MFREE ( nbytes, loc_TVIN )
c
      return
      end


      subroutine M$FINTERPD (HCB,
     +             NOUT, TSTRO, TDELO, TVOUT, NCOMP, NORDER, NERR)
c
      include 'headers.inc'
c
      record /HEADER/ HCB
      SAVE
      REAL*8      TVOUT(*), TSTRO, TDELO
      INTEGER*4   NCOMP,NORDER,NERR,NIN,NOUT
      character*2 format
      REAL*8      X_MIN,X_MAX,X_RAN,X_START,X_DELTA
c
c get necessary file specs and error check
      if (HCB.format.eq.'NH') call M$ERROR
     & ('M$FINTERPD does not support non-homogeneous files')
      format = 'S'//HCB.format(2:2)
      NBPA  = M$BPA(format)
      NCOMP = HCB.bpe/NBPA
      HCB.buf_type = 'D'
c
c find range of output assuming monotonically increasing
      X_MIN = TSTRO
      X_MAX = TSTRO + (NOUT-1)*TDELO
      X_RAN = X_MAX-X_MIN
c
c set proper abscissa fields
      if (HCB.class .eq. 2) then
	X_START = HCB.ystart
	X_DELTA = HCB.ydelta
      else
	X_START = HCB.xstart
	X_DELTA = HCB.xdelta
      endif
c
c calculate read indices for grab and get data
      HCB.offset = dnint( (X_MIN-X_START)/X_DELTA - 2 )
      HCB.offset = max( 0.d0, min( HCB.size-1.d0, HCB.offset) )
      NGET = dnint( (X_MAX-X_START)/X_DELTA+3 ) - HCB.offset
      NGET = min( NGET, INT(HCB.size-HCB.offset) )
      X_START = X_START + HCB.offset*X_DELTA
      NTOT = NGET*NCOMP
      nbytes = max(8,NBPA) * NTOT
      call M$MALLOC ( nbytes, loc_TVIN )
      call M$GRAB (HCB, %val(loc_TVIN), HCB.offset+1.d0, NGET, -1)
c
c call interp routine
      call M$INTERPD ( NGET, X_START, X_DELTA, %val(loc_TVIN),
     +             NOUT, TSTRO, TDELO, TVOUT, NCOMP, NORDER, NERR)
c
c free memory and report errors
      call M$MFREE ( nbytes, loc_TVIN )
c
      return
      end



      subroutine M$INTERP (NIN, TSTRI, TDELI, TVIN,
     +             NOUT, TSTRO, TDELO, TVOUT, NCOMP, NORDER, NERR)
c
      IMPLICIT REAL*4 (A-H,O-Z)
      SAVE
      PARAMETER (MAXL=100)
      REAL*8    TSTRI, TDELI, TSTRO, TDELO, TDELINV, TT
      INTEGER*4 K,KK,KLAST,N,NVEC,NCOMP,NORDER,NERR,NIN,NOUT
      DIMENSION TVIN(NCOMP,*), TVOUT(NCOMP,*)
      DIMENSION A0(MAXL),A1(MAXL),A2(MAXL),A3(MAXL)
      LOGICAL*4 HOLD,LINEAR
      KLAST = -1
c
      NERR = 0
      TDELINV = 1.D0/TDELI
c
c     Begin Looping over NOUT vectors
      do NVEC = 1,NOUT
c
c       Calculate appropriate interval
        TT = TSTRO + DBLE(NVEC-1)*TDELO
        K  = (TT-TSTRI)*TDELINV + 1.D0
        K  = max0( K, 0 )
        K  = min0( K, NIN )
	KK = max0( K, 1 )
        T  = ( TT - ((K-1)*TDELI+TSTRI) )*TDELINV
c
c       Calculate coefficients
	if (K.ne.KLAST) then
	  HOLD = NORDER.eq.0 .or. K.eq.0 .or. K.eq.NIN
	  LINEAR = NORDER.eq.1 .or. K.eq.1 .or. K.eq.NIN-1
	  do N = 1,NCOMP
	    aa0 = TVIN(N,KK)
	    A0(N) = aa0
	    if (.not. HOLD) THEN
	       yy1 = TVIN(N,K+1)
	       A1(N) = yy1 - aa0
	       if (.not. LINEAR) THEN
	    	  aa1 = 5.E-1*( yy1-TVIN(N,K-1) )
	    	  A1(N) = aa1
	    	  aa3 = 5.E-1*( TVIN(N,K+2)+3.E0*aa0 ) - yy1 - yy1 + aa1
	    	  A3(N) = aa3
	    	  A2(N) = yy1 - aa1 - aa3 - aa0
	       endif				!.not LINEAR
            endif				!.not HOLD
  	  enddo
	endif
c
c       Calculate interpolated values
	do N = 1,NCOMP
	  s = A0(N)
	  if (HOLD) goto 20
	  s = s + A1(N)*T
	  if (LINEAR) goto 20
	  s = s + T*T*( A2(N) + A3(N)*T )
  20	  TVOUT(N,NVEC) = s
	enddo
c
c	Update last and error parameters
	KLAST = K
	if (HOLD) NERR = 2
	if (LINEAR) NERR = max0( NERR, 1)
c
      enddo
c
      return
      end



      subroutine M$INTERPD (NIN, TSTRI, TDELI, TVIN,
     +             NOUT, TSTRO, TDELO, TVOUT, NCOMP, NORDER, NERR)
c
      IMPLICIT REAL*8 (A-H,O-Z)
      SAVE
      PARAMETER (MAXL=100)
      REAL*8    TSTRI, TDELI, TSTRO, TDELO, TDELINV, TT
      INTEGER*4 K,KK,KLAST,N,NVEC,NCOMP,NORDER,NERR,NIN,NOUT
      DIMENSION TVIN(NCOMP,*), TVOUT(NCOMP,*)
      DIMENSION A0(MAXL),A1(MAXL),A2(MAXL),A3(MAXL)
      LOGICAL*4 HOLD,LINEAR
      KLAST = -1
c
      NERR = 0
      TDELINV = 1.D0/TDELI
c
c     Begin Looping over NOUT vectors
      do NVEC = 1,NOUT
c
c       Calculate appropriate interval
        TT = TSTRO + DBLE(NVEC-1)*TDELO
        K  = (TT-TSTRI)*TDELINV + 1.D0
        K  = max0( K, 0 )
        K  = min0( K, NIN )
	KK = max0( K, 1 )
        T  = ( TT - ((K-1)*TDELI+TSTRI) )*TDELINV
c
c       Calculate coefficients
	if (K.ne.KLAST) then
	  HOLD = NORDER.eq.0 .or. K.eq.0 .or. K.eq.NIN
	  LINEAR = NORDER.eq.1 .or. K.eq.1 .or. K.eq.NIN-1
	  do N = 1,NCOMP
	    aa0 = TVIN(N,KK)
	    A0(N) = aa0
	    if (.not. HOLD) THEN
	    	yy1 = TVIN(N,K+1)
	    	A1(N) = yy1 - aa0
	    	if (.not. LINEAR) THEN
	    	    aa1 = 5.D-1*( yy1-TVIN(N,K-1) )
	    	    A1(N) = aa1
	    	    aa3 = 5.D-1*( TVIN(N,K+2)+3.D0*aa0 ) - yy1 - yy1 + aa1
	    	    A3(N) = aa3
	    	    A2(N) = yy1 - aa1 - aa3 - aa0
	     	endif			!.not. LINEAR
	    endif		!.not. HOLD
  	  enddo
	endif
c
c       Calculate interpolated values
	do N = 1,NCOMP
	  s = A0(N)
	  if (HOLD) goto 20
	  s = s + A1(N)*T
	  if (LINEAR) goto 20
	  s = s + T*T*( A2(N) + A3(N)*T )
  20	  TVOUT(N,NVEC) = s
	enddo
c
c	Update last and error parameters
	KLAST = K
	if (HOLD) NERR = 2
	if (LINEAR) NERR = max0( NERR, 1)
c
      enddo
c
      return
      end
c***********************************************************************
c  FUNCTION:
c    M$FFTSIZE_BEST()
c
c  SYNOPSIS:
c    integer*4 function M$FFTSIZE_BEST(nfft,type,search)
c
c  ARGUMENTS:
c    integer*4 nfft   --- desired FFT size
c    integer*4 type   --- FFT type (0=Complex, 1=Real)
c    integer*4 search --- direction to search for best FFT size
c                           Negative -> return size <= nfft
c                           Positive -> return size >= nfft
c                           Zero     -> return size closest to nfft
c
c  PURPOSE:
c    Returns nearest efficient FFT size (relative to the specified FFT
c    size) which can be computed by the MIDAS FFT library.
c
c  REVISIONS:
c    07/05/01    rdh     baseline
c    09/30/04    MGJ     DR 503888 - Corrected assignment of fftflags
c
        integer*4 function M$FFTSIZE_BEST(nfft,rcflag,searchdir)
        implicit none
        include 'libnmfft_f77.inc'

        integer*4 nfft          !Desired FFT size
        integer*4 rcflag        !1=RealFFT, 0=ComplexFFT
        integer*4 searchdir     !Negative -> return size <= nfft
                                !Positive -> return size >= nfft
                                !Zero     -> return size closest to nfft
        integer*4 nopt, fftflags

c       Create the appropriate flags
        fftflags = FFT_FLAGS_COMPLEX
        if (rcflag .eq. 1) fftflags = FFT_FLAGS_REAL

c       Initialize the optimal FFT size and call NM_FFTSIZE
        nopt = -1
        call NM_FFTSIZE_F77(nfft,fftflags,searchdir,nopt)
        M$FFTSIZE_BEST = nopt

        return
        end


c***********************************************************************
c  FUNCTION:
c    M$FFTSIZE()
c
c  SYNOPSIS:
c    integer*4 function M$FFTSIZE(nfft,type,search)
c
c  ARGUMENTS:
c    integer*4 nfft   --- desired FFT size
c    integer*4 type   --- FFT type (0=Complex, 1=Real)
c    integer*4 search --- direction to search for best FFT size
c                           Negative -> return size <= nfft
c                           Positive -> return size >= nfft
c                           Zero     -> return size closest to nfft
c
c  PURPOSE:
c    Returns nearest efficient FFT size (relative to the specified FFT
c    size) which can be computed by the MIDAS FFT library.
c
c  REVISIONS:
c    07/05/01    rdh     baseline
c    09/30/04    MGJ     DR 503888 - Corrected assignment of fftflags
c
        integer*4 function M$FFTSIZE(nfft,rcflag,searchdir)
        implicit none
c       include 'libnmfft_f77.inc'

	integer*4 nfft		!Desired FFT size
	integer*4 rcflag  	!1=RealFFT, 0=ComplexFFT
	integer*4 searchdir	!Negative -> return size <= nfft
				!Positive -> return size >= nfft
				!Zero     -> return size closest to nfft

	integer*4 nopt, fftflags, M$FFTSIZE_BEST

c       Create the appropriate flags
c       fftflags = FFT_FLAGS_COMPLEX
c       if (rcflag .eq. 1) fftflags = FFT_FLAGS_REAL

c       Initialize the optimal FFT size and call NM_FFTSIZE
        nopt = -1
c       call NM_FFTSIZE_F77(nfft,fftflags,searchdir,nopt)
        nopt = M$FFTSIZE_BEST(nfft,rcflag,searchdir)
        M$FFTSIZE = nopt

	return
	end


c***********************************************************************
c  FUNCTION:
c    M$FFTFLAGS()
c
c  SYNOPSIS:
c    integer*4 function M$FFTFLAGS(cflags)
c
c  ARGUMENTS:
c    character*(*) cflags --- string containing name list of desired flags
c
c  PURPOSE:
c    Returns a bitmapped integer which represents each of the specified
c    flag names. The bitmapped integer can then be used in the M$FFT()
c    subroutine call. This function is just a shorthand way to specify
c    the flags in human readable form.
c
c  REVISIONS:
c    07/05/01    rdh     baseline
c
        integer*4 function M$FFTFLAGS(cflags)
        implicit none
	include 'libnmfft_f77.inc'

	character*(*) cflags	!String of delimited input flags
	integer*4     iflags	!Bitmapped version of flags
	character*80  cwrk	!Working string
	integer*4     i,ilength
	integer*4     M$LENGTH, M$SEARCH

c	Create working copy of input string
	do i = 1,80
	  cwrk(i:i) = ' '
	enddo
	ilength = M$LENGTH(cflags)
	ilength = MIN(MAX(0,ilength),80)
	cwrk = cflags(1:ilength)
	call M$UPPERCASE(cwrk)

c	Now, form bitmapped flag
	iflags = FFT_FLAGS_NULL
	if (M$SEARCH(cwrk,'REAL') .gt. 0)
     -	  iflags = IOR(iflags,FFT_FLAGS_REAL)
	if (M$SEARCH(cwrk,'COMP') .gt. 0)
     -	  iflags = IOR(iflags,FFT_FLAGS_COMPLEX)
	if (M$SEARCH(cwrk,'FOR') .gt. 0)
     -	  iflags = IOR(iflags,FFT_FLAGS_FORWARD)
	if (M$SEARCH(cwrk,'INV') .gt. 0)
     -	  iflags = IOR(iflags,FFT_FLAGS_INVERSE)
	if (M$SEARCH(cwrk,'NONYQ') .gt. 0)
     -	  iflags = IOR(iflags,FFT_FLAGS_NONYQUIST)
	if (M$SEARCH(cwrk,'SCALE') .gt. 0)
     -	  iflags = IOR(iflags,FFT_FLAGS_SCALE)
	if (M$SEARCH(cwrk,'UNPACK') .gt. 0) then
     	  iflags = IOR(iflags,FFT_FLAGS_UNPACKED)
	else
	  if (M$SEARCH(cwrk,'PACK') .gt. 0)
     -	    iflags = IOR(iflags,FFT_FLAGS_PACKED)
	endif

c	These flags should NOT be used except, perhaps, by XMLIB
	if (M$SEARCH(cwrk,'INITPLANS') .gt. 0)
     -	  iflags = IOR(iflags,FFT_FLAGS_INITPLANS)
	if (M$SEARCH(cwrk,'FREEPLANS') .gt. 0)
     -	  iflags = IOR(iflags,FFT_FLAGS_FREEPLANS)
	if (M$SEARCH(cwrk,'VALIDPLAN') .gt. 0)
     -	  iflags = IOR(iflags,FFT_FLAGS_VALIDPLAN)

	M$FFTFLAGS = iflags

        return
        end

#if _NMFFT
c***********************************************************************
c  SUBROUTINE:
c    M$RFFT_NM()
c
c  SYNOPSIS:
c    M$RFFT_NM(afft,nfft,mflag)
c
c  ARGUMENTS:
c    real*4    afft(*) --- buffer containing data to be transformed
c    integer*4 nfft    --- desired real FFT size
c    integer*4 mflag   --- direction flag (0,1=Forward unpacked)
c                                         ( -1=Inverse unpacked)
c                                         (  2=Forward packed)
c                                         ( -2=Inverse packed)
c                                         (  3=Forward without Nyquist bin)
c
c  PURPOSE:
c    subroutine call to NM_FFT_F77 "C" wrapper to NextMidas FFT routine
c    Drop-in replacement for M$RFFT
c
c  REVISIONS:
c    07/05/01    rdh     baseline
c
        subroutine M$RFFT(afft,nfft,mflag)
        implicit none
	include 'libnmfft_f77.inc'

	real*4    afft(*)	!Vector of data to transform
	integer*4 nfft		!Desired FFT size
	integer*4 mflag		!Original MIDAS calling flag
	integer*4 iflags	!Bitmapped FFT flags for NM_FFT
	integer*4 irtnval	!Return status value (0=Success, -1=Failure)

c	Set the flags
	iflags = FFT_FLAGS_NULL
	iflags = IOR(iflags,FFT_FLAGS_REAL)
	if (mflag .ge. 0) then
	  iflags = IOR(iflags,FFT_FLAGS_FORWARD)
	  if (mflag .eq. 1) then
	    iflags = IOR(iflags,FFT_FLAGS_UNPACKED)
	  elseif (mflag .eq. 2) then
	    iflags = IOR(iflags,FFT_FLAGS_PACKED)
	  elseif (mflag .eq. 3) then
	    iflags = IOR(iflags,FFT_FLAGS_NONYQUIST)
	  else
	    iflags = IOR(iflags,FFT_FLAGS_UNPACKED)
	  endif
	else
	  iflags = IOR(iflags,FFT_FLAGS_INVERSE)
	  if (mflag .eq. -1) then
	    iflags = IOR(iflags,FFT_FLAGS_UNPACKED)
	  elseif (mflag .eq. -2) then
	    iflags = IOR(iflags,FFT_FLAGS_PACKED)
	  elseif (mflag .eq. -3) then
	    iflags = IOR(iflags,FFT_FLAGS_NONYQUIST)
	  else
	    iflags = IOR(iflags,FFT_FLAGS_UNPACKED)
	  endif
	endif

	irtnval = 0
	call NM_FFT_F77(afft,nfft,iflags,irtnval)
	if (irtnval .ne. 0) then
	  if (irtnval .eq. -2) then
	    call M$ERROR('Insufficient memory for FFT size')
	  else
	    call M$ERROR('Invalid FFT size/parameters')
	  endif
	endif

	return
	end
#else
	subroutine M$RFFT (C,N,F)
c
	real*4 C(*)
	integer*4 N,F,FF,M2,M3,M4,M5,loc_tab,loc_scr
c
c build the twiddle table and calc RADIX factors (if not already done)
	call M$FFTPREP (1,N,M2,M3,M4,M5,loc_tab,loc_scr)
c
c pack data for inverse?
	if (F.eq.-1) c(2) = c(n+1)
c
c call the FPS RFFT routine - assumes packed data
	FF = min(0,F)	! forward or inverse
	call M$RFFTM (C,%val(loc_scr),%val(loc_tab),0,M2,M3,M4,M5,N,FF,-1)
c
c unpack data for forward?
	if (F.eq.0 .or. F.eq.1) then
	  c(n+1) = c(2)	! set nyquist real
	  c(n+2) = 0.0	! zero nyquist imaginary
	  c(2) = 0.0	! zero DC imaginary
	elseif (F.eq.3) then
	  c(2) = 0.0	! just lose the nyquist
	endif
c
	return
	end
#endif

#if _NMFFT
c***********************************************************************
c  SUBROUTINE:
c    M$CFFT_NM()
c
c  SYNOPSIS:
c    M$CFFT_NM(afft,nfft,mflag)
c
c  ARGUMENTS:
c    real*4    afft(*) --- buffer containing data to be transformed
c    integer*4 nfft    --- desired complex FFT size
c    integer*4 mflag   --- direction flag (0,1=Forward, -1=Inverse)
c
c  PURPOSE:
c    subroutine call to NM_FFT_F77 "C" wrapper to NextMidas FFT routine
c    Drop-in replacement for M$CFFT
c
c  REVISIONS:
c    07/05/01    rdh     baseline
c
        subroutine M$CFFT(afft,nfft,mflag)
        implicit none
	include 'libnmfft_f77.inc'

	real*4    afft(*)	!Vector of data to transform
	integer*4 nfft		!Desired FFT size
	integer*4 mflag		!Original MIDAS calling flag
	integer*4 iflags	!Bitmapped FFT flags for NM_FFT
	integer*4 irtnval	!Return status value (0=Success, -1=Failure)

c	Set the flags
	iflags = FFT_FLAGS_NULL
	iflags = IOR(iflags,FFT_FLAGS_COMPLEX)
	if (mflag .ge. 0) then
	  iflags = IOR(iflags,FFT_FLAGS_FORWARD)
	else
	  iflags = IOR(iflags,FFT_FLAGS_INVERSE)
	endif

	irtnval = 0
	call NM_FFT_F77(afft,nfft,iflags,irtnval)
	if (irtnval .ne. 0) then
	  if (irtnval .eq. -2) then
	    call M$ERROR('Insufficient memory for FFT size')
	  else
	    call M$ERROR('Invalid FFT size/parameters')
	  endif
	endif

	return
	end
#else
	subroutine M$CFFT (C,N,F)
c
	real*4 C(*)
	integer*4 N,F,FF,M2,M3,M4,M5,loc_tab,loc_scr
c
c build the twiddle table and calc RADIX factors (if not already done)
	call M$FFTPREP (2,N,M2,M3,M4,M5,loc_tab,loc_scr)
c
c call the FPS CFFT routine
	FF = min(0,F)
	call M$CFFTM (C,%val(loc_scr),%val(loc_tab),0,M2,M3,M4,M5,N,FF,-1)
c
	return
	end
#endif

#if _NMFFT
c***********************************************************************
c  SUBROUTINE:
c    M$FFT()
c
c  SYNOPSIS:
c    M$FFT(afft,nfft,flags)
c
c  ARGUMENTS:
c    real*4    afft(*) --- buffer containing data to be transformed
c    integer*4 nfft    --- desired FFT size
c    integer*4 flags   --- bitmapped integer of FFT flags
c
c  PURPOSE:
c    subroutine call to NM_FFT_F77 "C" wrapper to NextMidas FFT routine
c
c  REVISIONS:
c    07/05/01    rdh     baseline
c
        subroutine M$FFT(afft,nfft,flags)
        implicit none

	real*4    afft(*)	!Vector of data to transform
	integer*4 nfft		!Desired FFT size
	integer*4 flags		!Bitmapped FFT flags (real/cmplx/etc.)
	integer*4 irtnval	!Return status value (0=Success, -1=Failure)

	irtnval = 0
	call NM_FFT_F77(afft,nfft,flags,irtnval)
	if (irtnval .ne. 0) then
	  if (irtnval .eq. -2) then
	    call M$ERROR('Insufficient memory for FFT size')
	  else
	    call M$ERROR('Invalid FFT size/parameters')
	  endif
	endif

	return
	end
#else

	subroutine M$FFTPREP (MODE,N,M2,M3,M4,M5,ltab,lscr)
c
	parameter (L_setups=32)
	integer*4 setup(L_setups,8)
	save loc_scr,ns,nscr,setup
	data ns,nscr/0,0/
        logical*4       pwr2
        integer*4       ier
c
c release dynamic memory associated with scratch and twiddle table buffers
	if (MODE.eq.0) then
	  call M$MFREE (NSCR,loc_scr)
	  do i = 1,NS
	    call M$MFREE (setup(i,8), setup(i,7) )
	  enddo
	  NS = 0
	  NSCR = 0
	  return
c
c build the twiddle table and calc RADIX factors
	else
	  do i = 1,NS
	    if (setup(i,1).eq.N .and. setup(i,6).eq.MODE) goto 10
	  enddo
	  if (NS.ge.L_setups) call M$WARNING
     &	  ('Too many different FFT sizes to store twiddle tables for')
	  i = min(NS+1,L_setups)

	  m = M$GENRADIX (N,setup(i,2),setup(i,3),setup(i,4),setup(i,5))
	  if (m.eq.-1) call M$ERROR ('Could not factor FFT radix by 2,3,4,5')
	  if (m.eq.-2) then
	  write(*,*) '** Mixed-radix FFTs must have at least one radix-2/2 **'
	  write(*,*) '** pass to maintain the MIDAS interpretation of the  **'
	  write(*,*) '** nyquist bin for backwards compatibility  DR#219   **'
	  call M$ERROR ('Must have minimum of 1 radix-2/4 pass in M$GENRADIX')
	  endif
          pwr2 = setup(i,3) + setup(i,5) .eq. 0
	  setup(i,1) = N
	  setup(i,6) = MODE
#if _CNVX
          if (pwr2) then
             NSCRN = 0
          else
             NSCRN = N*4*MODE
          endif
#else
	  NSCRN = N*4*MODE
#endif
	  if (NSCRN.gt.NSCR) then
	    call M$MFREE ( NSCR, loc_scr )
	    NSCR = NSCRN
	    call M$MALLOC ( NSCR, loc_scr )
	  endif
#if _CNVX
          if (.not. pwr2) then
             setup(i,8) = N*8
          else
            if (MODE.eq.1) setup(i,8) = (3*N/2)*4
            if (MODE.eq.2) setup(i,8) = (5*N/2)*4
          endif
#else
	  setup(i,8) = N*8
#endif
	  call M$MALLOC ( setup(i,8), setup(i,7) )
#if _CNVX
          if (.not.pwr2) then
             call M$GENTAB ( %val(setup(i,7)), N, 2-MODE)
          else
             if (MODE.eq.1) call SRC1FT (0, N, %val(setup(i,7)), -3, ier)
             if (MODE.eq.2) call C1DFFT (0, N, %val(setup(i,7)), -3, ier)
             if (ier.ne.0) then
                type *,'M$FFTPREP - SRC1FT/C1DFFT error = ', ier
                call M$ERROR ('M$FFTPREP - FFT initialization failed')
             endif
          endif
#else
	  call M$GENTAB ( %val(setup(i,7)), N, 2-MODE)
#endif
	  NS = i
	endif
c
c recall known FFT size/mode setup
  10	M2 = SETUP(i,2)
	M3 = SETUP(i,3)
	M4 = SETUP(i,4)
	M5 = SETUP(i,5)
	ltab = SETUP(i,7)
	lscr = loc_scr
c
	return
	end

#endif


	function M$GENRADIX (N,M2,M3,M4,M5)
c
c computes Radix-2,3,4,5 factorization of FFT size - JGS
	M2 = 0
	M3 = 0
	M4 = 0
	M5 = 0
	NN = N
	do while (NN.gt.1 .and. MOD(NN,4).eq.0)
	  M4 = M4+1
	  NN = NN/4
	enddo
	do while (NN.gt.1 .and. MOD(NN,2).eq.0)
	  M2 = M2+1
	  NN = NN/2
	enddo
	do while (NN.gt.1 .and. MOD(NN,3).eq.0)
	  M3 = M3+1
	  NN = NN/3
	enddo
	do while (NN.gt.1 .and. MOD(NN,5).eq.0)
	  M5 = M5+1
	  NN = NN/5
	enddo
	if (NN.ne.1) then
	  M$GENRADIX = -1
	elseif (M2+M4.eq.0) then
	  M$GENRADIX = -2
	else
	  M$GENRADIX = 1
	endif
	return
	end



	function M$POWER2GE (iarg)
c
c find next power of 2 greater than or equal to arg
	if (iarg .le. 1) then
	  M$POWER2GE = 1
	else if (iarg.lt.2**10) then
	  M$POWER2GE = 2
	  do while (M$POWER2GE.lt.iarg)
	    M$POWER2GE = M$POWER2GE + M$POWER2GE
	  enddo
	else if (iarg.le.2**30) then
	  M$POWER2GE = 2**10
	  do while (M$POWER2GE.lt.iarg)
	    M$POWER2GE = M$POWER2GE + M$POWER2GE
	  enddo
	else
	  call M$ERROR ('Argument out of range in M$POWER2GE')
	endif
c
	return
	end



	function M$POWER2LE (iarg)
c
c find last power of 2 less than or equal to arg
	if (iarg .le. 1) then
	  M$POWER2LE = 1
	else if (iarg.lt.2**10) then
	  M$POWER2LE = 2
	  do while (M$POWER2LE.le.iarg)
	    M$POWER2LE = M$POWER2LE + M$POWER2LE
	  enddo
	  M$POWER2LE = M$POWER2LE/2
	else if (iarg.le.2**30) then
	  M$POWER2LE = 2**10
	  do while (M$POWER2LE.le.iarg)
	    M$POWER2LE = M$POWER2LE + M$POWER2LE
	  enddo
	  M$POWER2LE = M$POWER2LE/2
	else
	  call M$ERROR ('Argument out of range in M$POWER2LE')
	endif
c
	return
	end



	subroutine M$REFORMAT (A,formin,B,typeout,natom)
c
	implicit none
	include 'types.inc'
	include 'parameters.inc'
	include 'bounds.inc'
	character formin*(*),typeout
	integer*4 L_bbuf
	parameter (L_bbuf=16)
	include 'buffers.inc'
	record /BUFFER/ A,B
	integer*4 i,j,k,n,natom,nscalar,nbytes
	integer*4 it1,it2,itypein,itypeout,ibyte
	byte p(4)
	equivalence (p,j)
#if _IEEE
	parameter (ibyte=4)
#else
	parameter (ibyte=1)
#endif
c
c convert format input to type and SPA
	if (LEN(formin).eq.1) then
	  it1 = ICHAR(formin)
	  nscalar = natom
	else
	  it1 = ICHAR(formin(2:2))
	  nscalar = natom * M$SPA(formin)
	endif
	it2 = ICHAR(typeout)
	goto 1
c
c just format type entry point
	entry M$RETYPE (A,itypein,B,itypeout,natom)
c
	it1 = itypein
	it2 = itypeout
	nscalar = natom
c
c check for degenerate cases or memory move
  1	if (nscalar.le.0) return
	if (it1.eq.it2) then
	  if (%loc(A) .eq. %loc(B)) return
	  nbytes = nscalar * M$NPT(CHAR(it1))
	  if (nbytes.lt.0) nbytes = (-nbytes+7)/8
	  call M$MOVE (A,B,nbytes)
	  return
	endif
c
c branch on float output
	if (it2.eq.MP_FLOAT) then

	  if (it1.eq.MP_BYTE) then
	    do n = nscalar,1,-1
	      B.fbuf(n) = A.bbuf(n)
	    enddo
	  else if (it1.eq.MP_OFFSET) then
	    call M$REFORMAT_OFF (A,B,nscalar)
	    do n = nscalar,1,-1
	      B.fbuf(n) = B.bbuf(n)
	    enddo
	  else if (it1.eq.MP_INT) then
	    do n = nscalar,1,-1
	      B.fbuf(n) = A.ibuf(n)
	    enddo
	  else if (it1.eq.MP_LONG) then
	    do n = 1,nscalar
	      B.fbuf(n) = A.lbuf(n)
	    enddo
	  else if (it1.eq.MP_DOUBLE) then
	    do n = 1,nscalar
	      B.fbuf(n) = A.dbuf(n)
	    enddo
	  else if (it1.eq.MP_PACKED) then
	    call M$BIT2BYTE (A,B,nscalar)
	    do n = nscalar,1,-1
	      B.fbuf(n) = B.bbuf(n)
	    enddo
	  else
	    call M$ERROR('Unsupported input type in M$REFORMAT: '
     +			 //CHAR(it1))
	  endif
c
c branch on double output
	else if (it2.eq.MP_DOUBLE) then

	  if (it1.eq.MP_BYTE) then
	    do n = nscalar,1,-1
	      B.dbuf(n) = A.bbuf(n)
	    enddo
	  else if (it1.eq.MP_OFFSET) then
	    call M$REFORMAT_OFF (A,B,nscalar)
	    do n = nscalar,1,-1
	      B.dbuf(n) = B.bbuf(n)
	    enddo
	  else if (it1.eq.MP_INT) then
	    do n = nscalar,1,-1
	      B.dbuf(n) = A.ibuf(n)
	    enddo
	  else if (it1.eq.MP_LONG) then
	    do n = nscalar,1,-1
	      B.dbuf(n) = A.lbuf(n)
	    enddo
	  else if (it1.eq.MP_FLOAT) then
	    do n = nscalar,1,-1
	      B.dbuf(n) = A.fbuf(n)
	    enddo
	  else if (it1.eq.MP_PACKED) then
	    call M$BIT2BYTE (A,B,nscalar)
	    do n = nscalar,1,-1
	      B.dbuf(n) = B.bbuf(n)
	    enddo
	  else
	    call M$ERROR('Unsupported input type in M$REFORMAT: '
     +			 //CHAR(it1))
	  endif
c
c branch on long integer output
	else if (it2.eq.MP_LONG) then

	  if (it1.eq.MP_BYTE) then
	    do n = nscalar,1,-1
	      B.lbuf(n) = A.bbuf(n)
	    enddo
	  else if (it1.eq.MP_OFFSET) then
	    call M$REFORMAT_OFF (A,B,nscalar)
	    do n = nscalar,1,-1
	      B.lbuf(n) = B.bbuf(n)
	    enddo
	  else if (it1.eq.MP_INT) then
	    do n = nscalar,1,-1
	      B.lbuf(n) = A.ibuf(n)
	    enddo
	  else if (it1.eq.MP_FLOAT) then
	    do n = 1,nscalar
#if _SUN|!_SATURATE
	      B.lbuf(n) = nint(A.fbuf(n))
#else
	      B.lbuf(n) = nint(min(P_FMAXL,max(P_FMINL,A.fbuf(n))))
#endif
	    enddo
	  else if (it1.eq.MP_DOUBLE) then
	    do n = 1,nscalar
#if _SUN|!_SATURATE
	      B.lbuf(n) = nint(A.dbuf(n))
#else
	      B.lbuf(n) = nint(min(P_DMAXL,max(P_DMINL,A.dbuf(n))))
#endif
	    enddo
	  else if (it1.eq.MP_PACKED) then
	    call M$BIT2BYTE (A,B,nscalar)
	    do n = nscalar,1,-1
	      B.lbuf(n) = B.bbuf(n)
	    enddo
	  else
	    call M$ERROR('Unsupported input type in M$REFORMAT: '
     +			 //CHAR(it1))
	  endif
c
c branch on short integer output
	else if (it2.eq.MP_INT) then

	  if (it1.eq.MP_BYTE) then
	    do n = nscalar,1,-1
	      B.ibuf(n) = A.bbuf(n)
	    enddo
	  else if (it1.eq.MP_OFFSET) then
	    call M$REFORMAT_OFF (A,B,nscalar)
	    do n = nscalar,1,-1
	      B.ibuf(n) = B.bbuf(n)
	    enddo
	  else if (it1.eq.MP_LONG) then
	    do n = 1,nscalar
	      B.ibuf(n) = min(P_LMAXI,max(P_LMINI,A.lbuf(n)))
	    enddo
	  else if (it1.eq.MP_FLOAT) then
	    do n = 1,nscalar
	      B.ibuf(n) = nint(min(P_FMAXI,max(P_FMINI,A.fbuf(n))))
	    enddo
	  else if (it1.eq.MP_DOUBLE) then
	    do n = 1,nscalar
	      B.ibuf(n) = nint(min(P_DMAXI,max(P_DMINI,A.dbuf(n))))
	    enddo
	  else if (it1.eq.MP_PACKED) then
	    call M$BIT2BYTE (A,B,nscalar)
	    do n = nscalar,1,-1
	      B.ibuf(n) = B.bbuf(n)
	    enddo
	  else
	    call M$ERROR('Unsupported input type in M$REFORMAT: '
     +			 //CHAR(it1))
	  endif
c
c branch on byte or offset byte output
	else if (it2.eq.MP_BYTE .or. it2.eq.MP_OFFSET) then

	  if (it1.eq.MP_BYTE .or. it1.eq.MP_OFFSET) then
	    call M$REFORMAT_OFF (A,B,nscalar)
	    return
	  else if (it1.eq.MP_INT) then
	    do n = 1,nscalar
	      B.bbuf(n) = min(P_IMAXB,max(P_IMINB,A.ibuf(n)))
	    enddo
	  else if (it1.eq.MP_LONG) then
	    do n = 1,nscalar
	      B.bbuf(n) = min(P_LMAXB,max(P_LMINB,A.lbuf(n)))
	    enddo
	  else if (it1.eq.MP_FLOAT) then
	    do n = 1,nscalar
	      B.bbuf(n) = nint(min(P_FMAXB,max(P_FMINB,A.fbuf(n))))
	    enddo
	  else if (it1.eq.MP_DOUBLE) then
	    do n = 1,nscalar
	      B.bbuf(n) = nint(min(P_DMAXB,max(P_DMINB,A.dbuf(n))))
	    enddo
	  else if (it1.eq.MP_PACKED) then
	    call M$BIT2BYTE (A,B,nscalar)
	  else
	    call M$ERROR('Unsupported input type in M$REFORMAT: '
     +			 //CHAR(it1))
	  endif
	  if (it2.eq.MP_OFFSET)
     +	    call M$REFORMAT_OFF (B,B,nscalar)
c
c branch on packed bit output
	else if (it2.eq.MP_PACKED) then

	  nscalar = (nscalar+7)/8
	  if (it1.eq.MP_BYTE) then
	    do n = 1,nscalar
	      k = n*8
	      j = 0
	      do i = 7,0,-1
		if (A.bbuf(k-i).gt.0) j = JIBSET(j,i)
	      enddo
	      B.bbuf(n) = p(ibyte)
	    enddo
	  else if (it1.eq.MP_OFFSET) then
	    do n = 1,nscalar
	      k = n*8
	      j = 0
	      do i = 7,0,-1
		if (A.bbuf(k-i).le.0) j = JIBSET(j,i)
	      enddo
	      B.bbuf(n) = p(ibyte)
	    enddo
	  else if (it1.eq.MP_INT) then
	    do n = 1,nscalar
	      k = n*8
	      j = 0
	      do i = 7,0,-1
		if (A.ibuf(k-i).gt.0) j = JIBSET(j,i)
	      enddo
	      B.bbuf(n) = p(ibyte)
	    enddo
	  else if (it1.eq.MP_LONG) then
	    do n = 1,nscalar
	      k = n*8
	      j = 0
	      do i = 7,0,-1
		if (A.lbuf(k-i).gt.0) j = JIBSET(j,i)
	      enddo
	      B.bbuf(n) = p(ibyte)
	    enddo
	  else if (it1.eq.MP_FLOAT) then
	    do n = 1,nscalar
	      k = n*8
	      j = 0
	      do i = 7,0,-1
		if (A.fbuf(k-i).gt.0.) j = JIBSET(j,i)
	      enddo
	      B.bbuf(n) = p(ibyte)
	    enddo
	  else if (it1.eq.MP_DOUBLE) then
	    do n = 1,nscalar
	      k = n*8
	      j = 0
	      do i = 7,0,-1
		if (A.dbuf(k-i).gt.0.d0) j = JIBSET(j,i)
	      enddo
	      B.bbuf(n) = p(ibyte)
	    enddo
	  else
	    call M$ERROR('Unsupported input type in M$REFORMAT: '
     +			 //CHAR(it1))
	  endif

	else
	  call M$ERROR('Unsupported output type in M$REFORMAT: '
     +			 //CHAR(it2))
	endif
c
	return
	end



	subroutine M$REFORMAT_OFF (bufi,bufo,N)
	integer*4 bufi(*),bufo(*),N,buft,lxor
c
c convert offset byte to byte
c (compliment the sign bit of byte data) (in groups of 4 for speed)
	lxor = INT(Z'80808080')
	n4 = N/4
	do i = 1,n4
	  bufo(i) = IEOR(bufi(i),lxor)
	enddo
	i = N - n4*4
	if (i.gt.0) then
	  call M$MOVE (bufi(n4+1),buft,i)
	  buft = IEOR(buft,lxor)
	  call M$MOVE (buft,bufo(n4+1),i)
	endif
	return
	end


	subroutine M$BIT2BYTE (buf1,buf2,nbits)
	include 'parameters.inc'
	byte buf1(*),buf2(*)
	nbytes = nbits/8
	nleft = nbits - nbytes*8
	if (nleft.gt.0) then
	  j = buf1(nbytes+1)
	  k = (nbytes+1)*8
	  do i = 7,8-nleft,-1
	    buf2(k-i) = JIBITS (j,i,1)
	  enddo
	endif
	do n = nbytes,1,-1
	  j = buf1(n)
	  k = n*8
	  do i = 7,0,-1
	    buf2(k-i) = JIBITS (j,i,1)
	  enddo
	enddo
	return
	end



	subroutine M$BITSHIFT (buf,nbits,nbs,insert)
	integer*4 j,ibyte,nbits,nbs
	byte buf(*),p(4),insert
	equivalence (j,p)
#if _IEEE
	parameter (ibyte=4)
#else
	parameter (ibyte=1)
#endif
	if (nbs.gt.0) then
	  nbq = (nbits+nbs+7)/8
	  no = 8 - nbs
	  j1 = buf(1)
	  do n = 1,nbq
	    j = j1
	    j1 = buf(n+1)
	    j = JIOR( JISHFT(j,nbs), JIBITS(j1,no,nbs))
	    buf(n) = p(ibyte)
	  enddo
	elseif (nbs.lt.0) then
	  nbq = (nbits-nbs+7)/8
	  no = 8 + nbs
	  nbsi = -nbs
	  j1 = insert
	  j1 = JISHFT (j1,-no)
	  do n = 1,nbq
	    j = j1
	    j1 = buf(n)
	    j = JIOR( JISHFT(j,no), JIBITS(j1,nbsi,no))
	    buf(n) = p(ibyte)
	  enddo
	endif
	return
	end



	function M$RECAST (bin,tin,lin,bout,tout,lout)
c
	character tin,tout,type
	integer*4 lin,lout,ld,lt,ndo,npt,la
	real*8    dval
	include 'types.inc'
	include 'values.inc'
	record /VALUE/ bin,bout
        INTEGER*4 XX,YY                                    !MIB 12-Apil-2000

#if _PGF77
#else
	automatic type,ld,lt,ndo,npt,dval
#endif
c
	M$RECAST = -1
C	TYPE *, '...recasting '''//tin//''' to '''//tout//'''...'
c
	if (tout.eq.' ' .OR. tout.eq.'_') then
	  ! Output type is unknown
	  if (tin.le.'A') then	! (tin.lt.'A')==(tin.eq.' ')?
	    ! Input type is also unknown
	    if (bin.a(1:1).eq.'(' .and. bin.a(lin:lin).eq.')' .and.
     &		INDEX(bin.a(1:lin),',').gt.0 .and.
     &		INDEX(bin.a(1:lin),' ').eq.0 ) then
	      ! It's an atomic results expression
	      tout = 'D'
	      goto 7
	    endif
	    ld = M$GET_RESULT(bin.a(1:lin), bout, lout, 0, type, la)
	    if (ld .lt. 0) then
	      ! It's not a results parameter - just some value
	      if (M$A2D(bin.a(1:lin),dval).gt.0) then
		bout.d = dval
		tout = 'D'
		ld = 8
	      else
		tout = 'A'
		ld = min(lin,lout)
		call M$MOVE (bin,bout,ld)
	      endif
	    elseif (tout.eq.'_' .AND. type.eq.'A') then
	      ! It's an ASCII results parameter
	      if (M$A2D(bout.a(1:ld), dval) .gt. 0) then
		bout.d = dval
		tout = 'D'
		ld = 8
	      endif
	    else
	      if (type.eq.'T') type = 'L'
	      tout = type
	    endif
	  else
	    ! Input type is known - use that for the output type
	    tout = tin
	    ld = min(lin,lout)
	    call M$MOVE (bin,bout,ld)
	  endif

	elseif (tout.eq.tin) then
	  ! Input and output types match - no cast is necessary
	  ! (maybe the memory move isn't even necessary)
	  if (%loc(bin) .ne. %loc(bout)) then
	    ld = min(lin,lout)
	    call M$MOVE (bin,bout,ld)
	  else
	    ld = lin
	  endif

	elseif (tout.eq.'A') then
	  ! Convert to ASCII
	  if (tin.le.'A') then
	    ! Input type is unknown
	    ld = M$GET_RESULT( bin.a(1:lin), bout, lout, 0, type, la)
	    if (type.ne.'A') then
	      ld = lin
	      call M$MOVE (bin,bout,ld)
	    endif
	  else
	    ! Numeric input
	    npt = M$NPT(tin)
	    if (npt.ne.lin) then		! it's an array
	      i = 1
	      n = 0
	      ld = lin/npt
	      nn = min(ld,lout/2+1)
	      call M$MALLOC(nn*npt,lt)
	      call M$MOVE(bin,%val(lt),nn*npt)
	      bout.a(1:1) = '('
	      do while (i.lt.lout-8 .and. n.lt.nn)
		if (tin.eq.'D') then
		  j = M$D2A (%val(lt+n*8), bout.a(i+1:min(i+20,lout-2)) )
		elseif (tin.eq.'F') then
		  j = M$F2A (%val(lt+n*4), bout.a(i+1:min(i+20,lout-2)) )
		elseif (tin.eq.'L') then
		  j = M$L2A (%val(lt+n*4), bout.a(i+1:min(i+20,lout-2)) )
		else
		  call M$REFORMAT (%val(lt+n*npt),tin,dval,'L',1)
		  j = M$L2A (dval, bout.a(i+1:min(i+20,lout-2)) )
		endif
		n = n+1
		i = i+j+1
		bout.a(i:i) = ','
	      enddo
	      call M$FREE(lt)
	      if (n.eq.ld) then
		bout.a(i:i) = ')'
		ld = i
	      else
		ld = i+1
		bout.a(ld:ld) = '+'
	      endif
	    elseif (tin.eq.'D') then
	      ld = M$D2A (bin, bout.a(1:min(20,lout)))
	    elseif (tin.eq.'F') then
	      ld = M$F2A (bin, bout.a(1:min(20,lout)))
	    else
	      call M$REFORMAT (bin, tin, dval, 'L', 1)
	      ld = M$L2A (dval, bout.a(1:min(20,lout)))
	    endif
	  endif

	elseif (tin.le.'A') then
	  ! Output type is numeric, and input type is either
	  ! ASCII or unknown...
	  if (M$TYPE(bin.a(1:1)) .eq. 2) then	! Alphabetic
	    ld = M$GET_RESULT( bin.a(1:lin), bout, lout, 0, type, la)
	    if (ld.gt.0 .and. type.gt.'A') then
	      if (type.eq.'T') type = 'L'
	      if (tout.ne.type) then
		n = ld / M$NPT(type)
	        call M$REFORMAT (bout,type, bout,tout, n)
	        ld = n * M$NPT(tout)
	      endif
	    else
	      ! Either not a results parameter, or an ASCII results parameter
	      if (M$A2D(bin.a(1:lin),dval).lt.0) then
		call M$WARNING
     &		  ('Problem converting ascii to numeric: '//bin.a(1:lin))
		return
	      endif
	      call M$REFORMAT (dval,'D',bout,tout,1)
	      ld = M$NPT(tout)
	    endif
	  elseif (INDEX(bin.a(1:lin),',').gt.0) then
	    goto 7

	  elseif (M$A2D(bin.a(1:lin),dval).lt.0) then
	    call M$WARNING ('Problem converting ascii to numeric: '//
     &				bin.a(1:lin))
	    return
	  else
	    call M$REFORMAT (dval,'D',bout,tout,1)
	    ld = M$NPT(tout)
	  endif

	else
	  ! Input and ouput types are numeric and different
	  n = 1
	  type = tin
	  if (tin .eq. 'T') type = 'L'
	  npt = M$NPT(type)
	  if (npt.ne.lin) n = lin/npt
	  call M$REFORMAT (bin,type,bout,tout,n)
	  ld = n*M$NPT(tout)

	endif
	goto 9000

 7	j = INDEX(bin.a(1:lin),',')
	npt = M$NPT(tout)
	i = 0
	k = lin
	if (bin.a(1:1).eq.'(' .and. bin.a(k:k).eq.')') then
	  i = 1
	  k = lin-1
	endif
	n = 0
	nn = min(lin/2+1,lout/npt)
	call M$MALLOC(nn*npt,lt)
	do while (j.gt.i .and. n.lt.nn)
C...The following use of XX and YY are a work around to  a complier problem
C for a new ES40 with VMS 7.1-2 and latest fortran complier EV7.2-180-4293Q
C per jgk: 16-Jul-1999
          XX=i+1                             !MIB 12-Apil-2000
          YY=j-1                             !MIB 12-Apil-2000
          if (M$A2D(bin.a(XX:YY),dval).lt.0) then         !MIB 12-Apil-2000
		call M$WARNING ('Problem converting ascii to numeric: '
     &				//bin.a(XX:YY))     !MIB 12-Apil-200
		return
	  endif
	  call M$REFORMAT(dval,'D',%val(lt+n*npt),tout,1)
	  n = n+1
	  i = j
	  if (i.lt.k) then
	    j = INDEX(bin.a(i+1:k),',') + i
	    if (j.eq.i) j=k+1
	  endif
	enddo
	ld = n*npt
	call M$MOVE (%val(lt),bout,ld)


 9000	continue
	M$RECAST = ld
	return
	end



	function M$DEQT(d1, d2, tol)

	implicit none

	include 'types.inc'
c Arguments
	real*8		d1		! (IN)	first number to compare
	real*8		d2		! (IN)	second number to compare
	real*4		tol		! (IN)  tolerance within which to
					!	compare them
c
c Code
c
	if (d1 .eq. d2) then
	  M$DEQT = .TRUE.
	else if ((ABS(d1 - d2)) / MAX(ABS(d1), ABS(d2)) .le. tol) then
	  M$DEQT = .TRUE.
	else
	  M$DEQT = .FALSE.
	endif

	return

	end


	function M$DEQ(d1, d2)

	implicit none

	include 'midas.inc'
	include 'types.inc'
c Arguments
	real*8		d1, d2
c
c Code
c
	M$DEQ = M$DEQT(d1, d2, Mc.tolerance)
	return
	end



	function M$CALC (string, buf, nbuf, type, mode)

	implicit none
	include 'types.inc'
	include 'midas.inc'
	character string*(*), type
	integer*4 nbuf, mode
	real*8 buf(*)

	include 'aop.inc'
	record /AOPSTRUCT/ op

	integer*4 L_stk,ns,nss,i,ip,lp,ls,istat,rstat,npt,savetype
	parameter (L_stk=16)
	character modes*(L_stk), tmp*24, rtype
	real*8 dval

	structure /STKSTRUCT/
	  real*8    val(3)
	endstructure
	record /STKSTRUCT/ stk(L_stk)

	ns = 0
	ip = 0
	ls = LEN(string)
	npt = M$NPT(type)

	do while (ip.lt.ls)
	  lp = M$LENGTHD (string(ip+1:))

	  istat = M$AOP_FIND (op, string(ip+1:ip+lp), modes(1:ns), type)

	  if (istat.eq.P_AOP_VALID) then		! operator
	    nss =  ns - op.inputs
	    ! RANDOM needs a special exception to remain backwards compatible.
	    ! We need to call M$VRAND using Mc.seed, and subtract 0.5; the
	    ! M$AOP_RAND routine does neither. We're also relying on there
	    ! only being one output value.  Note that this means different
	    ! behavior for FCALC and CALC.
	    if (string(ip+1:ip+lp) .eq. 'RANDOM') then
	      call MD$VTYPEQ(savetype)
	      call M$VTYPE('D')
	      call M$VRAND(Mc.seed, dval, 1)
	      dval = dval - 0.5d0
	      call M$REFORMAT(dval, 'D', stk(nss+1).val, type, 1)
	      call MD$VTYPE(savetype)
	    else
	      istat = M$AOPF (%val(op.func.p), 1,
     &		stk(nss+1).val,1, stk(nss+2).val,1,
     &		stk(nss+3).val,1, stk(nss+4).val,1 )
	      if (istat.lt.0) call M$AOP_ERROR (istat,1)
	    endif
	    do i = 1, op.outputs
	      stk(nss+i) = stk(ns+i)
	    enddo
	    ns = ns - op.inputs + op.outputs
	    modes(ns:ns) = op.mode
	  elseif (istat.eq.P_AOP_BADTYPE) then
	    call M$ERROR('Data type not supported for op: '//
     &				string(ip+1:ip+lp))
	  else
	    ns = ns + 1
	    rtype = ' '
	    rstat = M$RECAST (%val(%loc(string)+ip),' ',lp,
     &				stk(ns).val,rtype,24)
	    if (rtype.eq.type .and. rstat.eq.npt) then
	      modes(ns:ns) = 'S'
	    elseif (rtype.gt.'A') then
	      i = rstat / M$NPT(rtype)
	      modes(ns:ns) = 'S'
	      if (i.eq.2) modes(ns:ns) = 'C'
	      if (i.eq.3) modes(ns:ns) = 'V'
	      if (rtype.ne.type) call M$REFORMAT
     &		(stk(ns).val,rtype, stk(ns).val,type, i)
	    else
	      tmp = string(ip+1:ip+lp)
	      if (istat.eq.P_AOP_BADMODE) call M$WARNING
     &   ('Data modes not supported by operator: '//tmp//' or,')
	      if (istat.eq.P_AOP_BADSTK) call M$WARNING
     &   ('Insufficient stack to support operator: '//tmp//' or,')
	      call M$ERROR ('Couldnt find result: '//tmp)
	    endif
	  endif
	  ip = ip + lp + 1
	enddo
c
c check for valid termination
	if (ns.ne.1) call M$ERROR ('Extra parameters left on stack')
	i = M$NPT(type) * M$NPT(modes(1:1))
	call M$MOVE (stk(1).val, buf, min(i,nbuf))
	M$CALC = i
	return
	end


