      SUBROUTINE CURFIT ( Y,IDIM,N,A,B,SIGMAA,SN,CHISQR,YFIT,NTERMS )
C   
C    
C     THIS ROUTINE COMPUTES A NON-LINEAR LEAST SQUARES FIT. 
C     TO A DOUBLE-GAUSSIAN ABSORBTION LINE.  
C     THE METHOD IS INTERMEDIATE BETWEEN A GRADIENT SEARCH AND   
C     A FIT TO THE CHI SQUARE FUNCTION. 
C    
C     FOR A DESCRIPTION OF THE ALGORYTHM, SEE
C    
C         DATA REDUCTION AND ERROR ANALYSIS FOR THE PHYSICAL SCIENCES.
C             PHILIP R. BEVINGTON, MCGRAW-HILL   1969  CHAPTER 11.    
C    
C     PARAMETERS.   
C     Y(I),I=-N,+N        ARRAY OF DATA POINTS.   
C     A(7) CONTINUUM INTENSITYF THE CENTER POSITION.   
C     A(1),A(4)           CENTRAL DEPTHS OF GAUSSAIN CENTERS.    
C     A(2),A(5)           CENTERS OF THE GUASSIANS.         
C     A(3),A(6)           WIDTHS OF THE GAUSSIANS.
C     SIGMAA              STANDARD DEVIATIONS OF THE PARAMETERS  
C     SN   STANDARD DEVIATION OF DATA        
C     CHISQR              CHI-SQUARE OF FIT. 
C     YFIT THE FITTED SPECTRUM
C     NTERMS              NUMBER OF TERMS IN THE FIT        
C            IF NTERMS = 4 THEN ONLY ONE GAUSSAIN IS FITTED. 
C    
C     OTHER REQUIREMENTS 
C     FUNCTN ( I, A ,NTERMS ) 
C             EVALUATES THE FITTING FUNCTION FOR THE I-TH TERM.  
C     FDERIV ( I,A,DERIV,NTERMS)        
C             EVALUATE DERIVATIVES OF THE FITTING FUNCTION. 
C     FCHISQ   ( Y,N0,N,NFREE,WEIGHT,YFIT )  
C         EVALUATES REDUCED CHI SQUARE FIT TO THE DATA.     
C     MATINV ( B,NTERMS )
C         INVERTS NTERMS BY NTERMS MATRIX B AND CALCULATES ITS   
C         ITS DETERMINENT DET.  THE MATRIX MUST BE SYMMETRIC.    
C    
C     THIS PARAMETER MUST BE MODIFIED IF ANOTHER FUNCTION IS TO BE FIT
C     NTERMS      NUMBER OF PARAMETERS IN FITTING FUNCTION  
C    
C     THE X COORDINATE IS TAKEN TO BE THE INDEX OF Y.       
C    
C     N0 IS THE CENTER OF THE FIT...TAKEN TO BE ZERO        
C    
C    
C     FLAMDA IS THE RATIO OF THE IMPORTANCE OF CURVATURE TO GRADIENT  
C             IN THE FITTING.
      REAL ARRAY(7,7)    
      DIMENSION Y(-IDIM:IDIM),YFIT(-IDIM:IDIM)    
      DIMENSION A(7),B(7),SIGMAA(7)     
      DIMENSION ALPHA(7,7),BETA(7),DERIV(7)  
      INTEGER         UOUT,URES         
      DATA UOUT,URES /  6 ,  9/         
      FLAMDA = .001 
      NFREE = 2*N + 1 - NTERMS
      IF ( NFREE .LE. 0 ) THEN
          CHISQR = 0.    
          WRITE(UOUT,*) ' NOT ENOUGH POINTS FOR A UNIQUE FIT. '  
          WRITE(URES,*) ' NOT ENOUGH POINTS FOR A UNIQUE FIT. '  
          RETURN    
      END IF        
C    
C     OUTPUT INITIAL PARAMETERS.        
      WRITE(URES,'(A6,7F10.3)' ) ' A =  ',A  
C    
C     *****************WEIGHTS ARE TAKEN TO BE UNIFORM
C     FOR NON UNIFORM WEIGHTS INTRODUCE THE ARRAY 
C     WEIGHTS AND INITIALIZE IT HERE. SEE BEVINGTON         
C     FOR A DISCRIPTION. 
C     *********************************************         
C    
C     START ITERATIONS.  FINISHED WHEN THE CHANGE IS CHI SQUARE IS    
C     LESS THAN ONE PART IN 1.E6         
C    
      NITER = 0     
      CHICHI = 1.E30
   35 CONTINUE      
C    
C     EVALUATE ALPHA AND BETA MATRICES  
C    
      DO 34 J = 1,NTERMS 
          BETA(J) = 0.   
          DO 34 K = 1,NTERMS  
   34 ALPHA(J,K) = 0.    
      DO 50 I = -N,N
          CALL FDERIV ( I,A,DERIV,NTERMS )   
          FCT = FUNCTN(I,A,NTERMS)      
          DO 46 J = 1,NTERMS  
              BETA(J) = BETA(J) +   (Y(I)-FCT) * DERIV(J)   
              DO 46 K = 1,J   
   46     ALPHA(J,K) = ALPHA(J,K) +  DERIV(J)  * DERIV(K)   
   50 CONTINUE      
      DO 53 J = 1,NTERMS 
          DO 53 K = 1,J  
   53     ALPHA(K,J) = ALPHA(J,K)       
C    
C         EVALUATE CHI SQUARE AT STARTING POINT   
C    
      DO 62 I = -N,N
   62     YFIT(I) = FUNCTN(I,A,NTERMS)  
      CHISQ1 = FCHISQ (Y,IDIM,N,NFREE,SN,YFIT )   
C    
C     INVERT MODIFIED CURVATURE MATRIX TO FIND NEW PARAMETERS    
C    
C     MATRIX IS AN IMPROVEMENT OVER BEVINGTON'S MATINV.     
C    
   71 DO 74 J = 1,NTERMS 
          DO 73 K = 1,NTERMS  
   73         ARRAY(J,K) = ALPHA(J,K) / SQRT( ALPHA(J,J) * ALPHA(K,K) )         
   74 ARRAY(J,J) = ( 1. + FLAMDA )      
      CALL MATINV ( ARRAY,NTERMS )      
      DO 84 J = 1,NTERMS 
          B(J) = A(J)    
          DO 84 K = 1,NTERMS  
   84 B(J)=B(J)+.8* BETA(K)*ARRAY(J,K)/ SQRT( ALPHA(J,J)*ALPHA(K,K) ) 
C    
C     IF CHI SQUARE INCREASED, INCREASE FLAMDA AND TRY AGAIN
C    
      DO 92 I = -N,N
   92     YFIT(I) = FUNCTN ( I,B,NTERMS )    
      CHISQR = FCHISQ ( Y,IDIM,N,NFREE,SN,YFIT )  
      IF ( CHISQ1 - CHISQR ) 95,95,101  
   95 FLAMDA = 10. * FLAMDA   
          IF ( FLAMDA .GT. 1.E10) THEN  
              WRITE(UOUT,'(A46)')       
     &              ':I^* CANNOT FIND LINE...FITTING ROUTNE ABORTED'  
              WRITE(URES,'(A46)')       
     &              '     CANNOT FIND LINE...FITTING ROUTINE ABORTED' 
              RETURN
          END IF    
          GO TO 71  
C    
C     EVALUATE PARAMETERS AND UNCERTAINTIES. 
C    
  101 DO 103 J = 1,NTERMS
          A(J) = B(J)    
  103 SIGMAA(J) = SQRT( ABS(ARRAY(J,J)/ALPHA(J,J) ) )       
      FLAMDA = FLAMDA / 10.   
      NITER = NITER + 1  
          IF ( NITER .GT. 15 ) THEN     
              WRITE(UOUT,'(A33)') ':I^* LINE FIT DOES NOT CONVERGE..' 
              WRITE(URES,'(A33)') '     LINE FIT DOES NOT CONVERGE..' 
              RETURN
          END IF    
      CHANGE = ABS ( 1. - CHISQR/CHICHI )    
      WRITE(URES,'(A6,7F10.3)' ) ' A =  ',A  
      CHICHI = CHISQR    
C    
C     IF LINE CENTER IS NO LONGER IN REGION QUIT. 
      IF ( MAX( ABS(A(2)),ABS(A(5)) ) .GT.   FLOAT(N)  )  THEN   
          DO 104 J = 1,7 
              A(J) = 0.  
  104     CONTINUE  
          CHISQR=0. 
          WRITE(UOUT,'(A37)') ':I^* LINE HAS MIGRATED OUT OF REGION ' 
          WRITE(URES,'(A37)') '     LINE HAS MIGRATED OUT OF REGION ' 
          RETURN    
      END IF        
C    
C     CHECK TO SEE IF PARAMETERS ARE GIVING TROUBLE.        
C    
      IF (( NTERMS .EQ. 7 ).OR.(NTERMS.EQ.6) ) THEN         
          IF ( ( A(3)*A(6) ) .LT. 0.)  THEN  
              NTERMS = NTERMS-3         
              WRITE(UOUT,'(A35)')':I^* A SINGLE GAUSSIAN WILL BE USED'
              WRITE(URES,'(A35)')'     A SINGLE GAUSSIAN WILL BE USED'
              IF ( A(3) .LT. 0. ) THEN  
                  A(1) = A(4) 
                  A(2) = A(5) 
                  A(3) = A(6) 
              END IF
              A(4) = A(7)
          END IF    
      END IF        
      IF ( A(3) .LT. 0. ) THEN
          WRITE(UOUT,'(A33)') ' WIDTH OF GAUSSIAN IS NEGATIVE...'
          WRITE(URES,'(A33)') ' WIDTH OF GAUSSIAN IS NEGATIVE...'
          WRITE(UOUT,'(A17)') ':I^* FIT ABORTED ' 
          WRITE(URES,'(A17)') '     FIT ABORTED ' 
          RETURN    
      END IF        
      IF ( CHANGE .GT. 1.E-4 ) GO TO 35 
      WRITE(UOUT,*) NITER ,' ITERATIONS REQUIRED '
      WRITE(URES,*) NITER ,' ITERATIONS REQUIRED '
      RETURN        
      END 
      FUNCTION FCHISQ ( Y,IDIM,N,NFREE,SN,YFIT )  
C    
C     EVALUATES REDUCED CHI SQUARE FIT TO DATA    
C    
C         FCHISQ = SUM ( (Y-YFIT)/SIGMA ) **2 / NFREE       
C    
C     PARAMETERS    
C     YARRAY OF DATA POINTS   
C     IDIM            DIMENSION FOR Y.  
C     NNUMBER OF POINTS ON EACH SIDE OF CENTER    
C     NFREE           NUMBER OF DEGREES OF FREEDOM
C     SN              THE WEIGHT OF EACH DATA POINT.        
C     YFIT            ARRAY OF CALCULATED VALUES OF Y       
C    
      DIMENSION  Y(-IDIM:IDIM),YFIT(-IDIM:IDIM)   
      FCHISQ = 0.   
      IF ( NFREE .GT. 0 ) THEN
C    
C         ACCUMULATE CHI SQUARE         
C         ALWAYS USE NO WEIGHTING       
          DO 30 I = -N,N 
   30         FCHISQ = FCHISQ +  (Y(I)-YFIT(I))*( Y(I)-YFIT(I) ) 
          FCHISQ = FCHISQ / ( FLOAT(NFREE) * SN*SN )        
      END IF        
   40 RETURN        
      END 
      FUNCTION FUNCTN ( I,A,NTERMS )    
C     EVALUATES A DOUBLE GAUSSIAN ABSORBTON ON A FLAT BACKGROUND.
C           I     THE X VALUE WHERE THE GAUSSIANS WILL BE EVALUATED   
C         THE PARAMETERS OF THE FUNCTION ARE GIVEN BY       
C                     A(J)    J = 1,7 
      DIMENSION A(7)
      A0 = A(1)     
      X0 = A(2)     
      D0 = A(3)     
      A1 = A(4)     
      X1 = A(5)     
      D1 = A(6)     
      IF ( ( NTERMS .EQ. 3 ) .OR. ( NTERMS .EQ. 4 ) ) THEN  
          Y0 = A(4) 
      ELSE
          Y0 = A(7) 
      END IF        
      FUNCTN = Y0   
      IF ( D0 .GT. .02 ) THEN 
          Z  = ( FLOAT(I) - X0 ) / D0   
          Z2 = 0.5 * Z * Z    
          IF ( Z2.LT.50. ) FUNCTN = FUNCTN - Y0 * A0 * EXP( -Z2 )
      END IF        
      IF ( ( NTERMS .EQ. 4 ) .OR. ( NTERMS .EQ. 3 ) ) RETURN
      IF ( D1 .GT. .02 ) THEN 
          Z  = ( FLOAT(I) - X1 ) / D1   
          Z2 = 0.5 * Z * Z    
          IF ( Z2.LT.50. ) FUNCTN = FUNCTN - Y0 * A1 * EXP ( -Z2 )    
      END IF        
      RETURN        
      END 
      SUBROUTINE FDERIV ( I,A,DERIV,NTERMS)  
C    
C     ANALYTICALLY EVALUATES THE DERIVATIVES OF A 
C     A DOUBLE GAUSSIAN ABSORBTION LINE.
C     LINE ON A FLAT BACKGROUND AT THE POINT X(I) 
C    
C     I       THE X VALUE FOR THE FUNCTION.  
C     A       ARRAY OF PARAMETERS  FOR THE FUNCTION.        
C     DERIV   DERIVATIVES OF THE FUNCTION    
C    
      DIMENSION A(7),DERIV(7) 
      A0 = A(1)     
      X0 = A(2)     
      D0 = A(3)     
      A1 = A(4)     
      X1 = A(5)     
      D1 = A(6)     
      IF ( ( NTERMS .EQ. 3 ) .OR. ( NTERMS .EQ. 4 ) ) THEN  
          Y0 = A(4) 
      ELSE
          Y0 = A(7) 
      END IF        
      DO 10 J = 1,7 
          DERIV(J) = 0.  
   10 CONTINUE      
      DERIV(4) =1.  
      IF ( D0 .GT. .02 ) THEN 
          Z  = (FLOAT(I) - X0) / D0     
          Z2 = 0.5 * Z*Z 
          IF ( Z2 .LT. 50. ) THEN       
              DERIV(1) = - Y0 * EXP ( -Z2 )  
              DERIV(4) = DERIV(4) - A0 * DERIV(1) / Y0      
              DERIV(2) = A0  * Z * DERIV(1) / D0  
              DERIV(3) =       Z * DERIV(2)  
          END IF    
      END IF        
      IF ( ( NTERMS .EQ. 3 ) .OR. ( NTERMS .EQ. 4 ) ) RETURN
C    
C     DERIV(7) = DERIV(4)
      DERIV(4) = 0. 
      IF ( D1 .GT. .02 ) THEN 
          Z  = (FLOAT(I) - X1 ) / D1    
          Z2 = 0.5 * Z*Z 
          IF ( Z2. LT. 50. ) THEN       
              DERIV(4) = - Y0 * EXP ( -Z2 )  
              DERIV(7) = DERIV(7) - A1 * DERIV(4) / Y0      
              DERIV(5) = A1 * Z * DERIV(4) / D1   
              DERIV(6) =      Z * DERIV(5)   
          END IF    
      END IF        
      RETURN        
      END 
      SUBROUTINE MATINV ( A,M )         
C    
C     REPLACES THE M BY M MATRIX A BY ITS INVERSE.
C    
      REAL A(7,7),AA(7,7)
      REAL MAX,PIVOT
      INTEGER C(7),R(7),RM,CM 
      DO 20 J = 1,M 
          C(J) = J  
          R(J) = J  
          DO 10 I = 1,M  
              AA(I,J) = A(I,J)
              A(I,J) = 0.    
   10     CONTINUE  
          A(J,J) = 1.    
   20 CONTINUE      
C    
      DO 100 I = 1,M
          MAX = 0.  
          DO 50 J = I,M  
C    
C             FIND THE MAXIMIAL PIVOT   
              DO 40 K = I,M   
                 IF ( ABS( AA(R(J),C(K)) ) .GE. MAX ) THEN 
                     MAX = AA( R(J),C(K) )  
                     RM = J  
                     CM = K  
                 END IF 
   40         CONTINUE   
   50     CONTINUE  
C    
          IC    = C(I)   
          C(I)  = C(CM)  
          C(CM) = IC
          IR    = R(I)   
          R(I)  = R(RM)  
          R(RM) = IR
          RM    = R(I)   
          CM    = C(I)   
          PIVOT = AA(RM,CM)   
C         DIVIDE THE PIVOT ROW BY THE PIVOT  
          DO 60 J = 1,M  
              AA( RM,C(J) ) = AA( RM,C(J) ) / PIVOT         
              A ( RM,C(J) ) = A ( RM,C(J) ) / PIVOT         
   60     CONTINUE  
          DO 90 J = 1,M  
              IF ( J .NE. RM ) THEN     
                  COEF = AA( J,CM )     
                  DO 70 K = 1,I-1       
                      A ( J,C(K) ) = A ( J,C(K) ) - A ( RM,C(K) ) * COEF
   70             CONTINUE              
                  DO 80 K = I,M         
                      AA( J,C(K) ) = AA( J,C(K) ) - AA( RM,C(K) ) * COEF
                      A ( J,C(K) ) = A ( J,C(K) ) - A ( RM,C(K) ) * COEF
   80             CONTINUE
              END IF
   90     CONTINUE  
  100 CONTINUE
C
C     ERROR CHECKING SHOULD GO HERE.    
      RETURN        
      END 
      FUNCTION FUNCTN ( I,A,NTERMS )    
C     EVALUATES A DOUBLE GAUSSIAN ABSORBTON ON A FLAT BACKGROUND.
C           I     THE X VALUE WHERE THE GAUSSIANS WILL BE EVALUATED   
C         THE PARAMETERS OF THE FUNCTION ARE GIVEN BY       
C                           A(J)    J = 1,7 
      DIMENSION A(7)
      A0 = A(1)     
      X0 = A(2)     
      D0 = A(3)     
      A1 = A(4)     
      X1 = A(5)     
      D1 = A(6)     
      IF ( ( NTERMS .EQ. 3 ) .OR. ( NTERMS .EQ. 6 ) ) THEN  
          Y0 = 1.   
      ELSE
          Y0 = A(NTERMS) 
      END IF        
      FUNCTN = Y0   
      IF ( D0 .GT. .02 ) THEN 
          Z  = ( FLOAT(I) - X0 ) / D0   
          Z2 = 0.5 * Z * Z    
          IF ( Z2.LT.50. ) FUNCTN = FUNCTN - Y0 * A0 * EXP( -Z2 )
      END IF        
      IF ( ( NTERMS .EQ. 4 ) .OR. ( NTERMS .EQ. 3 ) ) RETURN
      IF ( D1 .GT. .02 ) THEN 
          Z  = ( FLOAT(I) - X1 ) / D1   
          Z2 = 0.5 * Z * Z    
          IF ( Z2.LT.50. ) FUNCTN = FUNCTN - Y0 * A1 * EXP ( -Z2 )
      END IF        
      RETURN        
      END 
