C#######################################################################
C
C	PICLIB.FOR - Library for PICAX
C
C#######################################################################
C
	SUBROUTINE FILTER ( YIN, YOUT, NDATA, FLOW, FHIGH, A, COEF,
     #                  TEMP, NTERMS, IFLAG )
C
C**********************************************************************
C
C     DESCRIPTION:
C       THIS SUBROUTINE MAY BE USED AS A LOWPASS, HIGHPASS, BANDPASS
C       OR BANDSTOP NON-RECURSIVE FILTER FOR EVENLY SPACED DATA.
C       THE K-TH FILTERED DATA POINT IS CALCULATED FROM THE 
C       (K-NTERMS)-TH THROUGH (K+NTERMS)-TH DATA POINTS USING A
C       KAISER WINDOW.  NTERMS DATA POINTS AT EACH END OF THE DATA
C       REMAIN UNFILTERED DUE TO THE FILTER ALGORITHM.  FREQUENCIES
C       ARE EXPRESSED IN TERMS OF THE NYQUIST FREQUENCY, 1/2T, WHERE
C       T IS THE TIME BETWEEN DATA SAMPLES.
C
C     PARAMETERS:
C       YIN    - DATA ARRAY TO BE FILTERED
C       YOUT   - DATA ARRAY AFTER FILTERING (MAY BE THE SAME AS YIN)
C	NDATA  - NUMBER OF DATA POINTS IN YIN
C       FLOW   - LOWER FREQUENCY OF FILTER AS A FRACTION OF NYQUIST
C                FREQUENCY
C       FHIGH  - UPPER FREQUENCY OF FILTER AS A FRACTION OF NYQUIST
C                FREQUENCY
C       A      - SIZE OF GIBBS PHENOMENON WIGGLES IN -DB
C	         (50 IS A GOOD CHOICE)
C       COEF   - ARRAY OF CALCULATED FILTER COEFFICIENTS
C                (DIMENSION=NTERMS)
C       TEMP   - TEMPORARY ARRAY USED BY SUBROUTINE (DIMENSION=NTERMS)
C       NTERMS - NUMBER OF TERMS IN THE FILTER FORMULA
C       IFLAG  - FLAG INDICATING WHETHER FILTER COEFFICIENTS ARE TO
C                BE CALCULATED AGAIN.  IFLAG SHOULD BE SET EQUAL
C                TO ZERO ON THE FIRST CALL.
C
C     USAGE:
C       THE FOLLOWING CONDITIONS ARE NECESSARY FOR VARIOUS TYPES OF
C       FILTERS:
C         NO FILTERING - FLOW = 0,  FHIGH = 1
C         LOW  PASS    - FLOW = 0,  0 < FHIGH < 1
C         HIGH PASS    - 0 < FLOW < 1,  FHIGH = 1
C         BAND PASS    - 0 < FLOW < FHIGH < 1
C         BAND STOP    - 0 < FHIGH < FLOW < 1
C
C     PROGRAMMER:  ROBERT WALRAVEN     VERSION 2.1     2 APR 79
C
C**********************************************************************
C
	DIMENSION YIN(1), YOUT(1), COEF(1), TEMP(1)
	DATA PI/3.14159265/
C
C--------------------BRANCH IF COEFFICIENTS ALREADY CALCULATED---------
C
	IF (IFLAG .NE. 0) GO TO 30
C
C--------------------CALCULATE KAISER WEIGHTS--------------------------
C
	CALL KAISER (COEF, NTERMS, A)
C
C--------------------CALCULATE FILTER COEFFICIENTS---------------------
C
	DO 20 I=1,NTERMS
20	COEF(I) = COEF(I)*(SIN(PI*I*FHIGH)-SIN(PI*I*FLOW))/(PI*I)
C
C--------------------FILTER DATA---------------------------------------
C
30    IF (NDATA .LT. NTERMS*2+1) RETURN     !RETURN IF NDATA TOO SMALL
      STOP = 1.                             !IF BANDSTOP THEN STOP=1
      IF (FHIGH .GE. FLOW) STOP = 0.        !            ELSE STOP=0
      DO 40 I=1,NTERMS                      !SAVE INITIAL POINTS
40    TEMP(I) = YIN(I)
      DO 60 I=NTERMS+1,NDATA-NTERMS         !FILTER MIDDLE POINTS
      SUM = 0.0
      DO 50 J=1,NTERMS
50    SUM = SUM + COEF(J)*(YIN(I-J)+YIN(I+J))
60    YOUT(I-NTERMS) = SUM + (FHIGH-FLOW+STOP)*YIN(I)
      DO 70 I=NDATA-NTERMS+1,NDATA          !MOVE FINAL UNFILTERED PNTS
70    YOUT(I) = YIN(I)
      DO 80 I=NDATA-NTERMS,NTERMS+1,-1      !SHIFT MIDDLE POINTS
80    YOUT(I) = YOUT(I-NTERMS)
      DO 90 I=1,NTERMS                      !RESTORE INITIAL POINTS
90    YOUT(I) = TEMP(I)
      RETURN
      END
	SUBROUTINE KAISER (W, N, A)
C
C**********************************************************************
C
C     COMPUTES KAISER WEIGHTS W(N,K) FOR DIGITAL FILTERS.
C
C     PARAMETERS:
C       W - CALCULATED ARRAY OF KAISER WEIGHTS.
C       N - VALUE OF N IN W(N,K), I.E., NUMBER OF TERMS.
C       A - SIZE OF GIBBS PHENOMENON WIGGLES IN -DB.
C
C     PROGRAMMER:  ROBERT WALRAVEN   31 MARCH 79
C
C**********************************************************************
C
	DIMENSION W(1)
	IF (A .LE. 21.) ALPHA = 0.
	IF (A .GE. 50.) ALPHA = 0.1102*(A-8.7)
	IF (A .LT. 50.  .AND. A .GT. 21.)
     #     ALPHA = 0.5842*(A-21.)**0.4 + 0.07886*(A-21.)
	DENOM = BESI0(ALPHA)
	DO 10 K=1,N
	  ARG  = FLOAT(K)/FLOAT(N)
	  TEMP = BESI0(ALPHA*SQRT(1.-ARG*ARG))/DENOM
	  W(K) = TEMP
10	CONTINUE
	RETURN
	END
	FUNCTION BESI0 (X)
C
C----------------------------------------------------------------------
C
C     COMPUTES THE ZERO-TH ORDER MODIFIED BESSEL FUNCTION I(X)
C
C----------------------------------------------------------------------
C
	T = X/3.75
	IF (T.GT.1.) GO TO 10
	T = T*T
	BESI0=1.+T*(3.5156229+T*(3.089424+T*(1.2067492+T*(.2659732
     #      +T*(.0360768+T*.0045813)))))
	RETURN
10	T = 1./T
	BESI0=SQRT(1./X)*EXP(X)*(.39894228+T*(.01328592+T*.00225319
     #      +T*(-.00157565+T*(.00916281+T*(-.02057706
     #      +T*(.02635537+T*(-.01647633+T*.00392377)))))))
	RETURN
	END
      SUBROUTINE LOWPAS ( YIN, YOUT, NPTS, FCUTOF )
C
C----------------------------------------------------------------------
C
C     THIS SUBROUTINE REMOVES HIGH FREQUENCY COMPONENTS FROM EVENLY
C     SAMPLED DATA BY PASSING THE DATA THROUGH A 9-TERM LOW-PASS DIGITAL
C     FILTER WITH A 50 DB KAISER WINDOW.  THE CUTOFF FREQUENCY MAY BE
C     VARIED WITH THE VARIABLE 'FCUTOF'.   NINE DATA POINTS AT EACH
C     END OF THE DATA REMAIN UNFILTERED DUE TO THE FILTER ALGORITHM.
C
C     PARAMETERS:
C        YIN    - DATA ARRAY TO BE FILTERED
C        YOUT   - DATA ARRAY AFTER FILTERING (MAY BE THE SAME AS YIN)
C        NPTS   - NUMBER OF POINTS IN DATA ARRAY
C        FCUTOF - CUTOFF FREQUENCY OF FILTER AS A FRACTION OF
C                 THE NYQUIST FREQUENCY, 1/2T, WHERE T IS THE
C                 TIME BETWEEN DATA SAMPLES.
C
C     PROGRAMMER: ROBERT WALRAVEN, APPLIED SCIENCE, UCD
C                 VERSION 1.0:  8 MARCH 1979
C
C----------------------------------------------------------------------
C
      DIMENSION YIN(1), YOUT(1), W(9), Z(9), COEF(9)
      DATA PI /3.14159265/
      DATA W  /.97546667,.90469038,.79570758,.66053343,.51335543,
     #         .36854291,.23853482,.13233204,.05468065/
C
      IF (NPTS .LE. 18) RETURN            !RETURN IF TOO FEW PTS
      DO 10 I=1,9                         !MOVE INITIAL UNFILTERED PTS
10    Z(I) = YIN(I)
      DO 20 J=1,9                         !CALCULATE COEFFICIENTS
20    COEF(J) = W(J) * SIN(PI*J*FCUTOF) / (PI*J)
      DO 40 I=10,NPTS-9                   !LOOP ON DATA
      SUM=0.0
      DO 30 J=1,9                         !COMPUTE AND ADD 9 TERMS
30    SUM = SUM + COEF(J) * (YIN(I-J)+YIN(I+J))
40    YOUT(I-9) = SUM + FCUTOF * YIN(I)   !INSERT FIRST TERM
      DO 50 I=NPTS-8,NPTS                 !MOVE FINAL UNFILTERED PTS
50    YOUT(I) = YIN(I)
      DO 60 I=NPTS-9,10,-1                !SHIFT MIDDLE PTS BY 9
60    YOUT(I) = YOUT(I-9)
      DO 70 I=1,9                         !PUT IN FIRST 9 PTS
70    YOUT(I) = Z(I)
      RETURN
      END
	SUBROUTINE H OF Z NR (Y, NDATA, NTERMS, A, FLOW, FHIGH)
C
C**********************************************************************
C
C    COMPUTES TRANSFER FUNCTION OF NON-RECURSIVE FILTER.
C
C     FILTER DESIGN PARAMETERS:
C	Y      - ARRAY WHERE TRANSFER FUNCTION IS STORED
C	NDATA  - NUMBER OF ELEMENTS IN Y
C       NTERMS - NUMBER OF FILTER TERMS
C       A      - SIZE OF GIBBS PHENOMENON RIPPLES IN -DB
C       FLOW   - LOWER CUTOFF OF FILTER AS MULTIPLE OF NYQUIST FREQ.
C       FHIGH  - UPPER CUTOFF OF FILTER AS MULTIPLE OF NYQUIST FREQ.
C
C     PROGRAMMER:  ROBERT WALRAVEN     27 SEP 80
C
C**********************************************************************
C
	DIMENSION Y(1), COEF(50)
	DATA PI/3.14159265/
C
	CALL KAISER (COEF, NTERMS, A)
	DO 40 I=1,NTERMS
40	COEF(I) = COEF(I)*(SIN(PI*I*FHIGH)-SIN(PI*I*FLOW))/(PI*I)
	STOP = 0.
	IF (FHIGH .LT. FLOW) STOP = 1.
	CONST = PI / FLOAT(NDATA-1)
	DO 50 I=1,NDATA
	F = CONST * (I-1)
	Y(I) = FHIGH - FLOW + STOP
	DO 50 J=1,NTERMS
50	Y(I) = Y(I) + 2.*COEF(J)*COS(F*J)
	RETURN
	END
	SUBROUTINE COEFS (A, B, NORDER, FCUT1, FCUT2, ACHEB)
C
C	COMPUTES RECURSIVE DIGITAL FILTER COEFFICIENTS.
C
C	A,B ARE ARRAYS OF COEFFICIENTS TO BE RETURNED.  A AND B
C		SHOULD BE DIMENSIONED GREATER THAN OR EQUAL TO
C		NORDER+1 AND NORDER, RESPECTIVELY, IN THE MAIN PROGRAM.
C	NORDER IS THE ORDER OF THE FILTER.
C	FCUT1,FCUT2 ARE LOWER AND UPPER CUTOFFS OF FILTER, I.E., THE
C		FILTER PASSES FREQUENCIES BETWEEN FCUT1 AND FCUT2.
C		FOR DESIRED TYPE OF FILTER CHOOSE FCUT1 AND FCUT2
C		AS FOLLOWS:
C					FCUT1	FCUT2
C			LOWPASS		  0	0 TO 1
C			HIGHPASS	0 TO 1	  1
C			BANDPASS	0 TO 1	>FCUT1
C			BANDSTOP	0 TO 1	<FCUT1
C	IF ACHEB = 0, FILTER IS BUTTERWORTH
C	IF ACHEB > 0, ACHEB IS ATTENUATION FACTOR OF CHEBYCHEV FILTER
C	      IN DB WHEN ITYPE = 2
C
	DIMENSION A(1),B(1)
	COMMON /COEFS/NTYPE,CA,CK,PREAL,PIMAG,A0,B1,B0,C(5),D(5),NCOEF
C
	IF (NORDER .LT. 1) NORDER=1
	IF (FCUT1 .EQ. FCUT2) FCUT2 = FCUT1 + .2
	IF (FCUT1.LT.0. .OR. FCUT1.GT.1.) FCUT1=0.
	IF (FCUT2.LT.0. .OR. FCUT2.GT.1.) FCUT2=1.
	IF (ACHEB .LT. 0.) ACHEB = -ACHEB
	ITYPE = 1
	IF (ACHEB .GT. 0.) ITYPE = 2
C
	CALL CONST (NORDER,FCUT1,FCUT2,ITYPE,ACHEB)
	N = (NTYPE+1)/2
	N = N*NORDER
	DO 10 I=1,N
	A(I) = 0.
10	B(I) = 0.
	A(N+1) = 0.
	A(1) = 1
	B(1) = 1.
	N = (NORDER+1)/2
	DO 100 I=1,N
	II = I
	CALL POLE(II,NORDER,ITYPE,ACHEB)
	IF (2*I .LE. NORDER) CALL S TO Z
	CALL TRNSFM (II,NORDER)
	CALL NXT TRM (A,B,II,NORDER)
100	CONTINUE
	RETURN
	END
	SUBROUTINE CONST (NORDER,FCUT1,FCUT2,ITYPE,ACHEB)
	COMMON /COEFS/NTYPE,CA,CK,PREAL,PIMAG,A0,B1,B0,C(5),D(5),NCOEF
	DATA PI/3.14159265/
C
	GO TO (10,20),ITYPE
10	OMEGA = 3.**(0.5/FLOAT(NORDER))
	GO TO 100
20	X = SQRT(3./(10.**(0.1*ACHEB)-1.))
	X = ALOG(X+SQRT(X**2-1))/FLOAT(NORDER)
	OMEGA = 0.5*(EXP(X)+EXP(-X))
C
100	OMEGA = ATAN(OMEGA)
	D1 = PI*FCUT1/2.
	D2 = PI*FCUT2/2.
	IF (FCUT1 .GT. 0.) GO TO 200
	CA = SIN(OMEGA-D2)/SIN(OMEGA+D2)
	NTYPE = 1
	RETURN
C
200	IF (FCUT2 .LT. 1.) GO TO 300
	CA = COS(OMEGA-D1)/COS(OMEGA+D1)
	NTYPE = 2
	RETURN
C
300	CA = COS(D1+D2)/COS(D2-D1)
	IF (FCUT2 .LT. FCUT1) GO TO 400
	CK = SIN(OMEGA)*COS(D2-D1)/(COS(OMEGA)*SIN(D2-D1))
	NTYPE = 3
	RETURN
C
400	CK = SIN(OMEGA)*SIN(D1-D2)/(COS(OMEGA)*COS(D1-D2))
	NTYPE = 4
	RETURN
	END
	SUBROUTINE POLE (I, NORDER, ITYPE, ACHEB)
C
C	CALCULATES POLES OF RECURSIVE FILTER OF DESIRED TYPE
C
	COMMON /COEFS/NTYPE,CA,CK,PREAL,PIMAG,A0,B1,B0,C(5),D(5),NCOEF
	DATA PI,PCREAL,PCIMAG/3.14159265,0.,0./
C
	IF (ITYPE .EQ. 2) GO TO 10
	THETA = FLOAT(NORDER+2*I-1)*PI/FLOAT(2*NORDER)
	PCREAL = COS(THETA)
	PREAL = PCREAL
	PCIMAG = SIN(THETA)
	PIMAG = PCIMAG
	RETURN
C
10	IF (I .NE. 1) GO TO 20
	X = 1./SQRT(10.**(0.1*ACHEB)-1.)
	X = ALOG(X+SQRT(X**2+1.))/FLOAT(NORDER)
	PCREAL = -0.5*(EXP(X)-EXP(-X))
	PCIMAG =  0.5*(EXP(X)+EXP(-X))
C
20	THETA = FLOAT(2*I-1)*PI/FLOAT(2*NORDER)
	PREAL = PCREAL * SIN(THETA)
	PIMAG = PCIMAG * COS(THETA)
	RETURN
	END
	SUBROUTINE S TO Z
C
C	PERFORMS VARIABLE TRANSFORMATION FROM S TO Z
C
	COMMON/COEFS/NTYPE,CA,CK,PREAL,PIMAG,A0,B1,B0,C(5),D(5),NCOEF
C
	P2 = PREAL**2 + PIMAG**2
	DENOM = 1 + P2 -2.*PREAL
	QREAL = (1.-P2)/DENOM
	QIMAG = 2.*PIMAG/DENOM
	A0 = P2/DENOM
	B1 = -2.*QREAL
	B0 = QREAL**2 + QIMAG**2
	RETURN
	END
	SUBROUTINE TRNSFM (I,NORDER)
	COMMON /COEFS/NTYPE,CA,CK,PREAL,PIMAG,A0,B1,B0,C(5),D(5),NCOEF
	GO TO (10,20,30,40),NTYPE
10	C0=-CA
	C1 = 1.
	GO TO 25
20	C0 = CA
	C1 = -1.
25	C2 = 0.
	D1 = -CA
	D2 = 0.
	NCOEF = 3
	GO TO 100
30	C0 = (1.-CK)/(1.+CK)
	C1 = 2.*CA*CK/(1.+CK)
	C2 = -1.
	D1 = -C1
	D2 = -C0
	NCOEF = 5
	GO TO 100
40	C0 = (1.-CK)/(1.+CK)
	C1 = -2.*CA/(1.+CK)
	C2 = 1.
	D1 = C1
	D2 = C0
	NCOEF = 5
100	IF (2*I .GT. NORDER) GO TO 200
	C(5) = A0*(C2+D2)**2
	C(4) = 2.*A0*(C1+D1)*(C2+D2)
	C(3) = A0*((C1+D1)**2+2.*(1.+C0)*(C2+D2))
	C(2) = 2.*A0*(1.+C0)*(C1+D1)
	C(1) = A0*(1.+C0)**2
	D(5) = C2**2+B1*C2*D2+B0*D2**2
	D(4) = 2.*C1*C2+B1*(C1*D2+C2*D1)+2.*B0*D1*D2
	D(3) = 2.*C0*C2+C1**2+B1*(C2+C1*D1+C0*D2)+B0*(2.*D2+D1**2)
	D(2) = 2.*C0*C1+B1*(C1+C0*D1)+2.*B0*D1
	D(1) = C0**2+B1*C0+B0
	GO TO 300
200	NCOEF = (NCOEF+1)/2
	A0= PREAL/(PREAL-1.)
	B0 = (PREAL+1.)/(PREAL-1.)
	C(3) = A0*(C2+D2)
	C(2) = A0*(C1+D1)
	C(1) = A0*(C0+1.)
	D(3) = C2+B0*D2
	D(2) = C1+B0*D1
	D(1) = C0+B0
300	DO 400 J=1,NCOEF
	C(J) = C(J)/D(NCOEF)
400	D(J) = D(J)/D(NCOEF)
	RETURN
	END
	SUBROUTINE NXT TRM (A, B, I, NORDER)
	COMMON/COEFS/NTYPE,CA,CK,PREAL,PIMAG,A0,B1,B0,C(5),D(5),NCOEF
	DIMENSION A(1),B(1)
C
	N = NCOEF
	M = (NTYPE+1)/2
	IF (2*I .GT. NORDER) N=2*NCOEF-1
	J = (N-1)*I+1
	JOLD = J-N+1
	DO 20 K=J,1,-1
	AK = 0.
	BK = 0.
	DO 10 L=1,NCOEF
	L1 = K-L+1
	IF (L1 .GT. JOLD) GO TO 10
	IF (L1 .LT. 1) GO TO 10
	AK = AK + A(L1)*C(L)
	BK = BK + B(L1)*D(L)
10	CONTINUE
	A(K) = AK
	IF (K .LE. NORDER*M) B(K) = BK
20	CONTINUE
	RETURN
	END
	SUBROUTINE H OF Z R (A, B, NORDER, F, HABS, HPHASE)
C
C	COMPUTES TRANSFER FUNCTION FOR RECURSIVE DIGITAL FILTER.
C
C	A,B ARE ARRAYS OF FILTER COEFFICIENTS.
C	NORDER IS ORDER OF FILTER.
C	F IS FREQUENCY AT WHICH TRANSFER FUNCTION IS TO BE COMPUTED.
C		(F IS IN UNITS OF 1/2T, SO RANGE IS 0. TO 1.)
C	HABS,HPHASE ARE ABSOLUTE VALUE AND PHASE OF TRANSFER FUNCTION.
C
	DIMENSION A(1),B(1)
	DATA PI/3.14159265/
	ZR = COS(PI*F)
	ZI = SIN(PI*F)
	ZNR = 1.
	ZNI = 0.
	HNR = 0.
	HNI = 0.
	HDR = 0.
	HDI = 0.
	DO 10 I=1,NORDER
	HNR = HNR + A(I)*ZNR
	HNI = HNI + A(I)*ZNI
	HDR = HDR + B(I)*ZNR
	HDI = HDI + B(I)*ZNI
	ZTEMP = ZNR
	ZNR = ZNR*ZR - ZNI*ZI
	ZNI = ZTEMP*ZI + ZNI*ZR
10	CONTINUE
	HNR = HNR + A(NORDER+1)*ZNR
	HNI = HNI + A(NORDER+1)*ZNI
	HDR = HDR + ZNR
	HDI = HDI + ZNI
	T = HDR*HDR + HDI*HDI
	HR = (HNR*HDR + HNI*HDI)/T
	HI = (HNI*HDR - HNR*HDI)/T
	T = SQRT(HR*HR+HI*HI)
	HABS = T
	T = ATAN2(HI,HR) * 180./PI
	HPHASE = T
	RETURN
	END
C
C******************* F I T T I N G   R O U T I N E S ********************
C
C	SUBROUTINE CURFIT
C
C	PURPOSE
C	   MAKE A LEAST-SQUARES FIT TO A NON-LINEAR FUNCTION
C	     WITH A LINEARIZATION OF THE FITTING FUNCTION.
C
C	USAGE
C	   COMMON /FIT/ NPTS,NTERMS,A(10),SIGMAA(10),FLAMDA,CHISQR
C	   CALL CURFIT (S, Y, WEIGHT)
C
C	DESCRIPTION OF PARAMETERS
C	   NPTS   - NUMBER OF PAIRS OF DATA POINTS
C	   NTERMS - NUMBER OF PARAMETERS
C	   A      - ARRAY OF PARAMETERS
C	   SIGMAA - ARRAY OF STANDARD DEVIATIONS FOR PARAMETERS A
C	   FLAMDA - PROPORTION OF GRADIENT SEARCH INCLUDED
C	   CHISQR - REDUCED CHI SQUARE FOR FIT
C	   X      - ARRAY OF DATA POINTS FOR INDEPENDENT VARIABLE
C	   Y      - ARRAY OF DATA POINTS FOR DEPENDENT VARIABLE
C	   WEIGHT - ARRAY OF WEIGHTS FOR Y DATA POINTS
C
C	SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C	   FUNCTN (X, A, NTERMS)
C	     EVALUATES THE FITTING FUNCTION
C	   FCHISQ (X, Y, WEIGHT, PARAM)
C	     EVALUATES REDUCED CHI SQUARE FOR FIT TO DATA
C	   FDERIV (X, A, NTERMS)
C	     EVALUATES THE DERIVATIVES OF THE FITTING FUNCTION
C	     WITH RESPECT TO EACH PARAMETER
C	   MATINV
C	     INVERTS A SYMMETRIC TWO-DIMENSIONAL MATRIX OF DEGREE
C	     NTERMS AND CALCULATES ITS DETERMINANT
C
C	COMMENTS
C	   DIMENSION STATEMENT VALID FOR NTERMS UP TO 10
C	   SET FLAMDA = 0.001 AT BEGINNING OF SEARCH
C
C	VERSION 2.0     10-SEP-79
C
	SUBROUTINE CURFIT (X, Y, WEIGHT)
	DIMENSION X(1), Y(1), WEIGHT(1)
	DIMENSION ALPHA(10,10), BETA(10), B(10), C(10)
	COMMON /MATINV/ ARRAY(10,10), NORDER, DET
	COMMON /DERIV / DERIV(10),YMODEL
	COMMON /FIT   / NPTS,NTERMS,A(10),SIGMAA(10),FLAMDA,CHISQR
C
	NORDER = NTERMS
	IF (NORDER .GT. 10) RETURN
	CHISQR = 0.
	NFREE = NPTS - NTERMS
	IF (NPTS .LT. NTERMS) RETURN
	W1 = WEIGHT(1)
	DO 1 I=1,NPTS
1	WEIGHT(I) = WEIGHT(I)/W1
C
C	   EVALUATE ALPHA AND BETA MATRICES
C
10	DO 20 J=1, NTERMS
	BETA(J) = 0.
	DO 20 K=1, J
20	ALPHA(J,K) = 0.
	DO 40 I=1, NPTS
	CALL FDERIV (X(I), A, NTERMS)
	WI = WEIGHT(I)
	FACTOR = WI * (Y(I) - YMODEL)
	DO 30 J=1, NTERMS
	DJ = DERIV(J)
	BETA(J) = BETA(J) + FACTOR*DJ
	DO 30 K=1, J
30	ALPHA(J,K) = ALPHA(J,K) + WI*DJ*DERIV(K)
40	CONTINUE
	DO 50 J=2, NTERMS
	DO 50 K=1, J-1
50	ALPHA(K,J) = ALPHA(J,K)
C
C	   EVALUATE CHI SQUARE AT STARTING POINT
C
	CHISQ1 = FCHISQ (X, Y, WEIGHT, A)
C
C	   INVERT MODIFIED CURVATURE MATRIX TO FIND NEW PARAMETERS
C
70	DO 80 J=1, NTERMS
	TEMP = SQRT(ABS(ALPHA(J,J)))
80	C(J) = TEMP
	DO 90 J=1, NTERMS-1
	CJ = C(J)
	ARRAY(J,J) = 1.+FLAMDA
	DO 90 K=J+1, NTERMS
	ARRAY(J,K) = ALPHA(J,K) / (CJ*C(K))
90	ARRAY(K,J) = ARRAY(J,K)
	ARRAY(NTERMS,NTERMS) = 1. + FLAMDA
	CALL MATINV
	DO 100 J=1, NTERMS
	B(J) = A(J)
	CJ = C(J)
	DO 100 K=1, NTERMS
100	B(J) = B(J) + BETA(K)*ARRAY(J,K)/(CJ*C(K))
C
C	   IF CHI SQUARE INCREASED, INCREASE FLAMDA AND TRY AGAIN
C
	CHISQR = FCHISQ (X, Y, WEIGHT, B)
	IF (CHISQ1 .GE. CHISQR) GO TO 120
	FLAMDA = 10.*FLAMDA
	GO TO 70
C
C	   EVALUATE PARAMETERS AND UNCERTAINTIES
C
120	DO 130 J=1, NTERMS
	A(J) = B(J)
	TEMP = SQRT (ARRAY(J,J) / (W1*ALPHA(J,J)))
130	SIGMAA(J) = TEMP
	DO 140 I=1,NPTS
140	WEIGHT(I) = WEIGHT(I)*W1
	CHISQR = CHISQR * W1
	FLAMDA = FLAMDA/10.
	RETURN
	END
C	FUNCTION FCHISQ
C
C	PURPOSE
C	   EVALUATE REDUCED CHI SQUARE FOR FIT TO DATA
C	     FCHISQ = SUM (WEIGHT*(Y-YFIT)**2) / NFREE
C
C	USAGE
C	   COMMON /FIT/NPTS,NTERMS,A(10),SIGMAA(10),FLAMDA,CHISQR
C	   RESULT = FCHISQ (X, Y, WEIGHT, PARAM)
C
C	DESCRIPTION OF PARAMETERS
C	   NPTS   - NUMBER OF DATA POINTS
C	   NFREE  - NUMBER OF DEGREES OF FREEDOM
C	   X      - ARRAY OF DATA POINTS FOR INDEPENDENT VARIABLE
C	   Y      - ARRAY OF DATA POINTS FOR DEPENDENT VARIABLE
C	   WEIGHT - ARRAY OF WEIGHTS FOR DATA POINTS
C	   PARAM  - ARRAY OF FIT PARAMETERS
C
C	VERSION 2.0     10-SEP-79
C
	FUNCTION FCHISQ (X, Y, WEIGHT, PARAM)
	DIMENSION X(1), Y(1), WEIGHT(1), PARAM(1)
	COMMON /FIT/ NPTS,NTERMS,A(10),SIGMAA(10),FLAMDA,CHISQR
	CHISQ = 0.
	NFREE = NPTS - NTERMS
	IF (NFREE.GT.0) GO TO 10
	FCHISQ = 0.
	RETURN
C
C	   ACCUMULATE CHI SQUARE
C
10	DO 20 I=1, NPTS
20	CHISQ = CHISQ+WEIGHT(I)*(Y(I)-FUNCTN(X(I),PARAM,NTERMS))**2
C
C	   DIVIDE BY NUMBER OF DEGREES OF FREEDOM
C
	FCHISQ = CHISQ / FLOAT(NFREE)
	RETURN
	END
C	SUBROUTINE MATINV
C
C	PURPOSE
C	   INVERT A SYMMETRIC MATRIX AND CALCULATE ITS DETERMINANT
C
C	USAGE
C	   COMMON /MATINV/ ARRAY(10,10),NORDER,DET
C	   CALL MATINV
C
C	DESCRIPTION OF PARAMETERS
C	   ARRAY  - INPUT MATRIX WHICH IS REPLACED BY ITS INVERSE
C	   NORDER - DEGREE OF MATRIX (ORDER OF DETERMINANT)
C	   DET    - DETERMINANT OF INPUT MATRIX
C
C	COMMENTS
C	   DIMENSION STATEMENT VALID FOR NORDER UP TO 10
C
C	VERSION 2.0     10-SEP-79
C
	SUBROUTINE MATINV
	DIMENSION IK(10), JK(10)
	COMMON /MATINV/ ARRAY(10,10), NORDER, DET
C
	DET = 1.
	DO 130 K=1, NORDER
C
C	  FIND LARGEST ELEMENT ARRAY(I,J) IN REST OF MATRIX
C
	AMAX = 0.
10	DO 20 I=K, NORDER
	DO 20 J=K, NORDER
	IF (ABS(AMAX) .GT. ABS(ARRAY(I,J))) GO TO 20
	AMAX = ARRAY(I,J)
	IK(K) = I
	JK(K) = J
20	CONTINUE
C
C	   INTERCHANGE ROWS AND COLUMNS TO PUT AMAX IN ARRAY(K,K)
C
	IF (AMAX .NE. 0) GO TO 30
	DET = 0.
	RETURN
30	I = IK(K)
	IF (I - K) 10, 60, 40
40	DO 50 J=1, NORDER
	SAVE = ARRAY(K,J)
	ARRAY(K,J) = ARRAY(I,J)
50	ARRAY(I,J) = -SAVE
60	J = JK(K)
	IF (J - K) 10, 90, 70
70	DO 80 I=1, NORDER
	SAVE = ARRAY(I,K)
	ARRAY(I,K) = ARRAY(I,J)
80	ARRAY(I,J) = -SAVE
C
C	   ACCUMULATE ELEMENTS OF INVERSE MATRIX
C
90	DO 100 I=1, NORDER
	IF (I .NE. K) ARRAY(I,K) = -ARRAY(I,K) / AMAX
100	CONTINUE
	DO 110 I=1, NORDER
	DO 110 J=1, NORDER
	IF (I.NE.K .AND. J.NE.K)
     1		ARRAY(I,J) = ARRAY(I,J) + ARRAY(I,K)*ARRAY(K,J)
110	CONTINUE
	DO 120 J=1, NORDER
	IF (J .NE. K) ARRAY(K,J) = ARRAY(K,J) / AMAX
120	CONTINUE
	ARRAY(K,K) = 1./AMAX
130	DET = DET * AMAX
C
C	   RESTORE ORDERING OF MATRIX
C
	DO 170 L=1, NORDER
	K = NORDER - L + 1
	J = IK(K)
	IF (J .LE. K) GO TO 150
	DO 140 I=1, NORDER
	SAVE = ARRAY(I,K)
	ARRAY(I,K) = -ARRAY(I,J)
140	ARRAY(I,J) = SAVE
150	I = JK(K)
	IF (I .LE. K) GO TO 170
	DO 160 J=1, NORDER
	SAVE = ARRAY(K,J)
	ARRAY(K,J) = -ARRAY(I,J)
160	ARRAY(I,J) = SAVE
170	CONTINUE
	RETURN
	END
C
C********** F A S T   F O U R I E R   T R A N S F O R M ***************
C
C----------------------------------------------------------------------
C Fast Fourier Transform Package from
C	Programs for Digital Signal Processing
C	IEEE Press, 1979
C	Book number 0-87942-128-2
C
C Program Authors:  G.D. Bergland and M.T. Dolan
C                   Bell Laboratories, Murray Hill, New Jersey, 07974
C----------------------------------------------------------------------
C
C SUBROUTINE:  FAST
C REPLACES THE REAL VECTOR B(K), FOR K=1,2,...,N,
C WITH ITS FINITE DISCRETE FOURIER TRANSFORM
C----------------------------------------------------------------------
C
	SUBROUTINE FAST (B, N)
C
C THE DC TERM IS RETURNED IN LOCATION B(1) WITH B(2) SET TO 0.
C THEREAFTER THE JTH HARMONIC IS RETURNED AS A COMPLEX
C NUMBER STORED AS B(2*J+1) + I B(2*J+2).
C THE N/2 HARMONIC IS RETURNED IN B(N+1) WITH B(N+2) SET TO 0.
C HENCE, B MUST BE DIMENSIONED TO SIZE N+2.
C THE SUBROUTINE IS CALLED AS FAST(B,N) WHERE N=2**M AND
C B IS THE REAL ARRAY DESCRIBED ABOVE.
C
	DIMENSION B(2)
	COMMON /CONST/ PII,P7,P7TWO,C22,S22,PI2
C
	PII = 3.14159265
	PI8 = PII/8.
	P7 = 1./SQRT(2.)
	P7TWO = 2.*P7
	C22 = COS(PI8)
	S22 = SIN(PI8)
	PI2 = 2.*PII
	DO 10 I=1,15
	  M = I
	  NT = 2**I
	  IF (N.EQ.NT) GO TO 20
10	CONTINUE
	WRITE (5,9999)
9999	FORMAT(' N is not a power of two for FAST')
	STOP
20	N4POW = M/2
C
C DO A RADIX 2 ITERATION FIRST IF ONE IS REQUIRED.
C
	IF (M-N4POW*2) 40,40,30
30	NN = 2
	INT = N/NN
	CALL FR2TR (INT,B(1),B(INT+1))
	GO TO 50
40	NN = 1
C
C PERFORM RADIX 4 ITERATIONS.
C
50	IF (N4POW.EQ.0) GO TO 70
	DO 60 IT=1,N4POW
	NN = NN*4
	INT = N/NN
	CALL FR4TR(INT,NN,B(1),B(INT+1),B(2*INT+1),B(3*INT+1),
     1	B(1),B(INT+1),B(2*INT+1),B(3*INT+1))
60	CONTINUE
C
C PERFORM IN-PLACE REORDERING.
C
70	CALL FORD1 (M,B)
	CALL FORD2 (M,B)
	T = B(2)
	B(2) = 0.
	B(N+1) = T
	B(N+2) = 0.
	DO 80 IT=4,N,2
	  B(IT) = -B(IT)
80	CONTINUE
	RETURN
	END
C
C----------------------------------------------------------------------
C SUBROUTINE:  FSST
C FOURIER SYNTHESIS SUBROUTINE
C----------------------------------------------------------------------
C
	SUBROUTINE FSST (B,N)
C
C THIS SUBROUTINE SYNTHESIZES THE REAL VECTOR B(K), FOR
C K=1,2,...,N, FROM THE FOURIER COEFFICIENTS STORED IN THE
C B ARRAY OF SIZE N+2.  THE DC TERM IS IN B(1) WITH B(2) EQUAL
C TO 0.  THE JTH HARMONIC IS STORED AS B(2*J+1) + I B(2*J+2).
C THE N/2 HARMONIC IS IN B(N+1) WITH B(N+2) EQUAL TO 0.
C THE SUBROUTINE IS CALLED AS FSST(B,N) WHERE N=2**M AND
C B IS THE REAL ARRAY DISCUSSED ABOVE.
C
	DIMENSION B(2)
	COMMON /CONST/ PII,P7,P7TWO,C22,S22,PI2
C
	PII = 3.14159265
	PI8 = PII/8.
	P7 = 1./SQRT(2.)
	P7TWO = 2.*P7
	C22 = COS(PI8)
	S22 = SIN(PI8)
	PI2 = 2.*PII
	DO 10 I=1,15
	  M = I
	  NT = 2**I
	  IF (N.EQ.NT) GO TO 20
10	CONTINUE
	WRITE (5,9999)
9999	FORMAT(' N IS NOT A POWER OF TWO FOR FSST')
	STOP
20	B(2) = B(N+1)
	DO 30 I=4,N,2
	  B(I) = -B(I)
30	CONTINUE
C
C SCALE THE INPUT BY N
C
	DO 40 I=1,N
	  B(I) = B(I)/FLOAT(N)
40	CONTINUE
	N4POW = M/2
C
C SCRAMBLE THE INPUTS
C
	CALL FORD2(M,B)
	CALL FORD1(M,B)
C
	IF (N4POW.EQ.0) GO TO 60
	NN = 4*N
	DO 50 IT=1,N4POW
	  NN = NN/4
	  INT = N/NN
	  CALL FR4SYN(INT,NN,B(1),B(INT+1),B(2*INT+1),B(3*INT+1),
     1	B(1),B(INT+1),B(2*INT+1),B(3*INT+1))
50	CONTINUE
C
C DO A RADIX 2 ITERATION IF ONE IS REQUIRED
C
60	IF (M-N4POW*2) 80,80,70
70	INT = N/2
	CALL FR2TR(INT,B(1),B(INT+1))
80	RETURN
	END
C
C----------------------------------------------------------------------
C SUBROUTINE:  FR2TR
C RADIX 2 ITERATION SUBROUTINE
C----------------------------------------------------------------------
C
	SUBROUTINE FR2TR (INT, B0, B1)
	DIMENSION B0(2), B1(2)
	DO 10 K=1,INT
	  T = B0(K) + B1(K)
	  B1(K) = B0(K) - B1(K)
	  B0(K) = T
   10	CONTINUE
	RETURN
	END
C
C----------------------------------------------------------------------
C SUBROUTINE:  FR4TR
C RADIX 4 ITERATION SUBROUTINE
C----------------------------------------------------------------------
C
	SUBROUTINE FR4TR(INT,NN,B0,B1,B2,B3,B4,B5,B6,B7)
	DIMENSION L(15),B0(2),B1(2),B2(2),B3(2),B4(2),B5(2),B6(2),B7(2)
	COMMON /CONST/ PII,P7,P7TWO,C22,S22,PI2
	EQUIVALENCE (L15,L(1)),(L14,L(2)),(L13,L(3)),(L12,L(4)),
     *	(L11,L(5)),(L10,L(6)),(L9,L(7)),(L8,L(8)),(L7,L(9)),
     *	(L6,L(10)),(L5,L(11)),(L4,L(12)),(L3,L(13)),(L2,L(14)),
     *	(L1,L(15))
C
C JTHET IS A REVERSED BINARY COUNTER, JR STEPS TWO AT A TIME TO
C LOCATE THE REAL PARTS OF INTERMEDIATE RESULTS, AND JI LOCATES
C THE IMAGINARY PART CORRESPONDING TO JR.
C
	L(1) = NN/4
	DO 40 K=2,15
	  IF (L(K-1)-2) 10,20,30
10	L(K-1) = 2
20	L(K) = 2
	GO TO 40
30	L(K) = L(K-1)/2
40	CONTINUE
C
	PIOVN = PII/FLOAT(NN)
	JI = 3
	JL = 2
	JR = 2
C
	DO 120 J1=2,L1,2
	DO 120 J2=J1,L2,L1
	DO 120 J3=J2,L3,L2
	DO 120 J4=J3,L4,L3
	DO 120 J5=J4,L5,L4
	DO 120 J6=J5,L6,L5
	DO 120 J7=J6,L7,L6
	DO 120 J8=J7,L8,L7
	DO 120 J9=J8,L9,L8
	DO 120 J10=J9,L10,L9
	DO 120 J11=J10,L11,L10
	DO 120 J12=J11,L12,L11
	DO 120 J13=J12,L13,L12
	DO 120 J14=J13,L14,L13
	DO 120 JTHET=J14,L15,L14
	TH2 = JTHET - 2
	IF (TH2) 50,50,90
50	DO 60 K=1,INT
	  T0 = B0(K) + B2(K)
	  T1 = B1(K) + B3(K)
	  B2(K) = B0(K) - B2(K)
	  B3(K) = B1(K) - B3(K)
	  B0(K) = T0 + T1
	  B1(K) = T0 - T1
60	CONTINUE
C
	IF (NN-4) 120,120,70
70	K0 = INT*4+1
	KL = K0+INT-1
	DO 80 K=K0,KL
	  PR = P7*(B1(K)-B3(K))
	  PI = P7*(B1(K)+B3(K))
	  B3(K) = B2(K) + PI
	  B1(K) = PI - B2(K)
	  B2(K) = B0(K) - PR
	  B0(K) = B0(K) + PR
80	CONTINUE
	GO TO 120
C
90	ARG = TH2*PIOVN
	C1 = COS(ARG)
	S1 = SIN(ARG)
	C2 = C1**2 - S1**2
	S2 = C1*S1 + C1*S1
	C3 = C1*C2 - S1*S2
	S3 = C2*S1 + S2*C1
C
	INT4 = INT*4
	J0 = JR*INT4 + 1
	K0 = JI*INT4 + 1
	JLAST = J0 + INT - 1
	DO 100 J=J0,JLAST
	  K = K0 + J - J0
	  R1 = B1(J)*C1 - B5(K)*S1
	  R5 = B1(J)*S1 + B5(K)*C1
	  T2 = B2(J)*C2 - B6(K)*S2
	  T6 = B2(J)*S2 + B6(K)*C2
	  T3 = B3(J)*C3 - B7(K)*S3
	  T7 = B3(J)*S3 + B7(K)*C3
	  T0 = B0(J) + T2
	  T4 = B4(K) + T6
	  T2 = B0(J) - T2
	  T6 = B4(K) - T6
	  T1 = R1 + T3
	  T5 = R5 + T7
	  T3 = R1 - T3
	  T7 = R5 - T7
	  B0(J) = T0 + T1
	  B7(K) = T4 + T5
	  B6(K) = T0 - T1
	  B1(J) = T5 - T4
	  B2(J) = T2 - T7
	  B5(K) = T6 + T3
	  B4(K) = T2 + T7
	  B3(J) = T3 - T6
100	CONTINUE
C
	JR = JR + 2
	JI = JI - 2
	IF (JI-JL) 110,110,120
110	JI = 2*JR - 1
	JL = JR
120	CONTINUE
	RETURN
	END
C
C----------------------------------------------------------------------
C SUBROUTINE:  FR4SYN
C RADIX 4 SYNTHESIS
C----------------------------------------------------------------------
C
	SUBROUTINE FR4SYN (INT,NN,B0,B1,B2,B3,B4,B5,B6,B7)
	DIMENSION L(15),B0(2),B1(2),B2(2),B3(2),B4(2),B5(2),B6(2),B7(2)
	COMMON /CONST/ PII,P7,P7TWO,C22,S22,PI2
	EQUIVALENCE (L15,L(1)),(L14,L(2)),(L13,L(3)),(L12,L(4)),
     *	(L11,L(5)),(L10,L(6)),(L9,L(7)),(L8,L(8)),(L7,L(9)),
     *	(L6,L(10)),(L5,L(11)),(L4,L(12)),(L3,L(13)),(L2,L(14)),
     *	(L1,L(15))
C
	L(1) = NN/4
	DO 40 K=2,15
	  IF (L(K-1)-2) 10,20,30
10	  L(K-1) = 2
20	  L(K) = 2
	  GO TO 40
30	  L(K) = L(K-1)/2
40	CONTINUE
C
	PIOVN = PII/FLOAT(NN)
	JI = 3
	JL = 2
	JR = 2
C
	DO 120 J1=2,L1,2
	DO 120 J2=J1,L2,L1
	DO 120 J3=J2,L3,L2
	DO 120 J4=J3,L4,L3
	DO 120 J5=J4,L5,L4
	DO 120 J6=J5,L6,L5
	DO 120 J7=J6,L7,L6
	DO 120 J8=J7,L8,L7
	DO 120 J9=J8,L9,L8
	DO 120 J10=J9,L10,L9
	DO 120 J11=J10,L11,L10
	DO 120 J12=J11,L12,L11
	DO 120 J13=J12,L13,L12
	DO 120 J14=J13,L14,L13
	DO 120 JTHET=J14,L15,L14
	  TH2 = JTHET - 2
	  IF (TH2) 50,50,90
50	  DO 60 K=1,INT
	    T0 = B0(K) + B1(K)
	    T1 = B0(K) - B1(K)
	    T2 = B2(K)*2.0
	    T3 = B3(K)*2.0
	    B0(K) = T0 + T2
	    B2(K) = T0 - T2
	    B1(K) = T1 + T3
	    B3(K) = T1 - T3
60	CONTINUE
C
	IF (NN-4) 120,120,70
70	K0 = INT*4 + 1
	KL = K0 + INT - 1
	DO 80 K=K0,KL
	  T2 = B0(K) - B2(K)
	  T3 = B1(K) + B3(K)
	  B0(K) = (B0(K)+B2(K))*2.0
	  B2(K) = (B3(K)-B1(K))*2.0
	  B1(K) = (T2+T3)*P7TWO
	  B3(K) = (T3-T2)*P7TWO
80	CONTINUE
	GO TO 120
90	ARG = TH2*PIOVN
	C1 = COS(ARG)
	S1 = -SIN(ARG)
	C2 = C1**2 - S1**2
	S2 = C1*S1 + C1*S1
	C3 = C1*C2 - S1*S2
	S3 = C2*S1 + S2*C1
C
	INT4 = INT*4
	J0 = JR*INT4 + 1
	K0 = JI*INT4 + 1
	JLAST = J0 + JINT - 1
	DO 100 J=J0,JLAST
	  K = K0 + J - J0
	  T0 = B0(J) + B6(K)
	  T1 = B7(K) - B1(J)
	  T2 = B0(J) - B6(K)
	  T3 = B7(K) + B1(J)
	  T4 = B2(J) + B4(K)
	  T5 = B5(K) - B3(J)
	  T6 = B5(K) + B3(J)
	  T7 = B4(K) - B2(J)
	  B0(J) = T0 + T4
	  B4(K) = T1 + T5
	  B1(J) = (T2+T6)*C1 - (T3+T7)*S1
	  B5(K) = (T2+T6)*S1 + (T3+T7)*C1
	  B2(J) = (T0-T4)*C2 - (T1-T5)*S2
	  B6(K) = (T0-T4)*S2 + (T1-T5)*C2
	  B3(J) = (T2-T6)*C3 - (T3-T7)*S3
	  B7(K) = (T2-T6)*S3 + (T3-T7)*C3
100	CONTINUE
	JR = JR + 2
	JI = JI - 2
	IF (JI-JL) 110,110,120
110	JI = 2*JR - 1
	JL = JR
120	CONTINUE
	RETURN
	END
C
C----------------------------------------------------------------------
C SUBROUTINE:  FORD1
C IN-PLACE REORDERING SUBROUTINE
C----------------------------------------------------------------------
C
	SUBROUTINE FORD1 (M,B)
	DIMENSION B(2)
C
	K = 4
	KL = 2
	N = 2**M
	DO 40 J=4,N,2
	  IF (K-J) 20,20,10
10	  T = B(J)
	  B(J) = B(K)
	  B(K) = T
20	  K = K-2
	  IF (K-KL) 30,30,40
30	  K = 2*J
	  KL = J
40	CONTINUE
	RETURN
	END
C
C----------------------------------------------------------------------
C SUBROUTINE:  FORD2
C IN-PLACE REORDERING SUBROUTINE
C----------------------------------------------------------------------
	SUBROUTINE FORD2 (M,B)
	DIMENSION L(15), B(2)
	EQUIVALENCE (L15,L(1)),(L14,L(2)),(L13,L(3)),(L12,L(4)),
     *	(L11,L(5)),(L10,L(6)),(L9,L(7)),(L8,L(8)),(L7,L(9)),
     *	(L6,L(10)),(L5,L(11)),(L4,L(12)),(L3,L(13)),(L2,L(14)),
     *	(L1,L(15))
	N = 2**M
	L(1) = N
	DO 10 K=2,M
	  L(K) = L(K-1)/2
10	CONTINUE
	DO 20 K=M,14
	  L(K+1) = 2
20	CONTINUE
	IJ = 2
	DO 40 J1=2,L1,2
	DO 40 J2=J1,L2,L1
	DO 40 J3=J2,L3,L2
	DO 40 J4=J3,L4,L3
	DO 40 J5=J4,L5,L4
	DO 40 J6=J5,L6,L5
	DO 40 J7=J6,L7,L6
	DO 40 J8=J7,L8,L7
	DO 40 J9=J8,L9,L8
	DO 40 J10=J9,L10,L9
	DO 40 J11=J10,L11,L10
	DO 40 J12=J11,L12,L11
	DO 40 J13=J12,L13,L12
	DO 40 J14=J13,L14,L13
	DO 40 JI=J14,L15,L14
	  IF (IJ-JI) 30,40,40
30	  T = B(IJ-1)
	  B(IJ-1) = B(JI-1)
	  B(JI-1) = T
	  T = B(IJ)
	  B(IJ) = B(JI)
	  B(JI) = T
40	  IJ = IJ + 2
	RETURN
	END