module marq
 
   ! These are Koji Mukai's versions of the Marquardt routines.  
   ! For a good reference as to how it works, see Numerical Recipes.
   ! The conversion to Fortran 90 was originally a bit basic, but it 
   ! was improved by Tim Naylor and Lee Howells.

   implicit none

   ! Define large numbers.
   Real, Parameter, private :: n_Large=-1.0*huge(n_Large), p_Large=huge(P_Large)

   integer, private, parameter :: m_Par=32

   ! "The taming factor".
   Real, private :: tame_Factor = 1.0

   Logical, private :: cyclic_Data = .False.
   Real, private :: cycle_Data = 0.0

   Integer, dimension(m_Par), private :: list_Cyclic=0
   Real, dimension(m_Par), private :: tab_Period=0.0, tab_Bottom=0.0

   ! Initializes the free_list common block so that no parameter is fixed.
   Integer, dimension(m_Par), private :: Index = (/ 1, 2, 3, 4, 5, 6, 7, 8, &
   9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, &
   27, 28, 29, 30, 31, 32 /)
   Integer, dimension(m_Par), private :: Table_Free=0

   ! Initializes the link_status so that no parameter is linked.
   Integer, dimension(m_Par), private :: l_Index = (/ 1, 2, 3, 4, 5, 6, 7, 8, &
   9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, &
   27, 28, 29, 30, 31, 32 /)
   Real, dimension(m_Par), private :: l_Offset=0.0
   
   ! Set up the default boundary within which parameters are confined.
   Real, dimension(m_par), private :: lo_A = n_Large, hi_A = p_Large
     
   Logical Debug
   Common / bug / Debug
       
   contains

       Subroutine COV_SORT( Covar, n_Par, n_Term )
 
       ! Repacks the covariance matrix into convenient form

       Real, intent(inout), dimension(:,:) ::  Covar
       Integer, intent(in) :: n_Par, n_Term

       Real Work
       Integer I, J, K, L

       ! Zero all elemsnts below diagonal
       Do J = 1, n_Par - 1
         Do I = J + 1, n_Par
           Covar( I, J ) = 0.0
         End Do
       End Do
                             ! Repack off-diagnoal to correct location
       Do I = 1, n_Term - 1
         Do J = I + 1, n_Term
           K = Index( I )
           L = Index( J )
           If( L .gt. K ) Then
             Covar( L, K ) = Covar( I, J )
           Else
             Covar( K, L ) = Covar( I, J )
           End If
         End Do
       End Do
                             ! Swap diagnoals using the top row
       Work = Covar( 1, 1 )
       Do J = 1, n_Par
         Covar( 1, J ) = Covar( J, J )
         Covar( J, J ) = 0.0
       End Do
       K = Index( 1 )
       Covar( K, K ) = Work
       Do J = 2, n_Term
         K = Index( J )
         Covar( K, K ) = Covar( 1, J )
       End Do
                             ! Finally, fill in above diagonals
       Do J = 2, n_Par
         Do I = 1, J - 1
           Covar( I, J ) = Covar( J, I )
         End Do
       End Do
 
       End Subroutine COV_SORT



       Subroutine CURFIT( X, Y, W, Points, mode_W, A, dA, n_Par, &
                                  n_Term, Covar, d_Chi, b_Chi, &
                                 sig_Clip, m_Rep, Flag, FUNCTN, FDERIV )
 
       ! A new version of CURFIT; see write-up

       ! Number of data points
       Integer, intent(in) :: Points
       ! X-data, Y-data
       Real, dimension(Points), intent(in) :: X, Y
       ! Weight
       real, dimension(Points), intent(inout) ::  W
       ! Weight mode.  1 Poisson, 0 equal weight, -1 explicit
       Integer, intent(in) :: mode_W
       ! Number of parameters
       Integer, intent(in) :: n_Par
       ! The parameters
       Real, intent(inout) :: A( n_Par )
       ! Step for numerical differenciation
       Real, intent(in) :: dA( n_Par )
       ! Number of free parameters
       Integer, intent(in) :: n_Term
       ! The covariance matrix
       Real, intent(out) :: Covar(:,:)
       ! Fractional Chi change for convergence
       Real, intent(in) :: d_Chi
       ! The final (best) reduced Chi squared
       Real, intent(out) :: b_Chi
       ! Iterative clipping flag.  If >0, do iterative clipping
       Real, intent(in) :: sig_Clip
       ! Maximum number of trials
       Integer, intent(in) :: m_Rep
       ! Controls output; error handling.
       Integer, intent(inout) :: Flag
       !   0: ok; 
       !  -1: not converged; 
       !  -2: lambda too big;
       ! -51: Too many iterations in svd_comp
       ! -91: initial parameter out of range; 
       ! -92: no D.O.F.
       ! -93: chi-squared is zero/
       ! User-supplied model function and its derivatives.       
       real, external :: FUNCTN
       external FDERIV

       Integer :: n_Free, J, K, Reject, flag_mrq
       Real Lambda, Temp
       Real o_Chi, Chi, Thres, imp_Chi
       Logical Interactive

       Call SETBUG( )

       ! Set all the outputs to real numbers, so alternate returns don't
       ! cause NaN errrors.
       covar=1.0
       b_chi=0.0
       
       If( Flag .gt. 1000 ) Then
         Interactive = .False.
       Else
         Interactive = .True.
       End If
       Flag = 0

       Do K = 1, n_Par
         If( A( K ) .le. lo_A( K ) .or. A( K ) .ge. hi_A( K ) ) Then
           Flag = -91
           Return
         End If
       End Do

       Do K = 1, n_Par
         If( list_Cyclic( K ) .eq. 1 ) Then       ! actively cyclic
           If( A( K ) .lt. tab_Bottom( K ) ) Then
             Write( *, 90 )
 90             Format( ' WARNING:: Parameter #', I2, ' out of ''cyclic'' range' / &
                      '           It will be modified before fitting...' )
             A( K ) = tab_Bottom( K ) + tab_Period( K ) &
                           + Mod( A( K ) - tab_Bottom( K ), tab_Period( K ) )
           Else If( A( K ) .ge. tab_Bottom( K ) + tab_Period( K ) ) Then
             Write( *, 90 )
             A( K ) = tab_Bottom( K ) &
                          + Mod( A( K ) - tab_Bottom( K ), tab_Period( K ) )
           End If
         End If
       End Do

       Reject = 1
       If( mode_W .gt. 0 ) Then       ! Poisson weighting
         Do J = 1, Points
           If( Y( J ) .gt. 0 ) Then
             W( J ) = 1.0 / Y( J )
           Else                            ! Ignore negative points
             W( J ) = 0.0
           End If
         End Do
       Else If( mode_W .eq. 0 ) Then       ! Equal weighting
         Do J = 1, Points
           W( J ) = 1.0
         End Do
       ! Else                            ! Predetermined (W, not sigma)
       End If

       Do While( Reject .gt. 0 )       ! Iterative Clipping Loop
         imp_Chi = 1000.0
                                   ! First, degree of freedom
         n_Free = 0
         Do J = 1, Points
           If( abs(W(J)) > tiny(W(J)) ) Then
             n_Free = n_Free + 1
           End If
         End Do
         n_Free = n_Free - n_Term
         If( n_Free .le. 0 ) Then
           Write( *, 120 ) Points, n_Term, n_Free
120           Format( ' ERROR:: No degrees of freedom in CURFIT' / &
                    '         # of points =', I5, ' # of free parameters =', &
                    I2, ' DOF =', I3 )
           Flag = -92
           Return
         End if
         ! First call to MARQUARDT
         Lambda = -0.001
         ! Make sure there's no junk in covar.
         covar = 0.0
         Call MARQUARDT( X, Y, W, Points, Lambda, A, dA, Covar, &
                            o_Chi, n_Par, n_Term, FUNCTN, FDERIV, flag )
         if (flag /= 0) return
         o_Chi = o_Chi / n_Free
         If( Interactive ) Then
           Write( *, 520 ) 0, o_Chi, Lambda, ( A( J ), J = 1, n_Par )
520           Format( ' Iteration # ', I2, ' Red. Chi-squ = ', 1PE10.3, &
                          ' Lambda = ', E8.1 / ( 5( ' ', E10.3 ) ) )
         End If

         ! Iteration
         Do K = 1, m_Rep
           Call MARQUARDT( X, Y, W, Points, Lambda, A, dA, &
                                 Covar, Chi, n_Par, n_Term, &
                                                      FUNCTN, FDERIV, flag )
           if (flag /= 0) return
           Chi = Chi / n_Free
           If( Interactive ) Then
             Write( *, 520 ) K, Chi, Lambda, ( A( J ), J = 1, n_Par )
           End If
           if (chi <= tiny(chi)) then
             print*, ' WARNING:: Chi-squared is zero.'
             flag = -93
             goto 300
           end if
           ! Small bug (which has a huge effect) spotted by Timn and Tariq.
           ! This line used to read .gt., which meant that if the fit had
           ! converged so well that o_chi was equal to chi you never tested
           ! against the convergence condition.
           ! If( o_Chi .ge. Chi ) Then
           ! But even this caused the g95 compiler problems, so I had to
           ! use this slightly laxer condition.
           if (o_chi - chi >= -1.0*epsilon(chi)) then
             imp_Chi = ( o_Chi - Chi ) / Chi
             If( imp_Chi .lt. d_Chi ) Goto 300
             o_Chi = Chi
           End If
           If( Lambda .ge. 1.0E+30 ) Then
             print*, o_chi, chi, o_chi-chi
             Write( *, 530 ) K, imp_Chi
530             Format( ' WARNING:: After ', I3, 'th iteration, lambda', &
                      ' is becoming too big' / &
                      '           The last fractional improvement in reduced', &
                      ' chi-squared is ', 1PE8.2 )
             If( imp_Chi .lt. 1.0E-03 ) Then
              Write( *, 531 )
531              Format( '           I would consider it as having converged!' )
             Else If( imp_Chi .lt. 0.1 ) Then
              Write( *, 532 )
532              Format( '           I might consider it as having converged' )
             Else
              Write( *, 533 )
533              Format( ' ERROR:: Curfit should not behave like this!' )
             End If
             Flag = -2
             Goto 300
           End If
         End Do
         ! .not. converged
         Write( *, 540 ) m_Rep
540         Format( ' WARNING:: Not converged after ', I2, ' iterations!' )
         Flag = -1

300      Continue
         if (debug) print*, '@ Converged, with flag ', flag
         b_Chi = Chi
         ! Check for clipping
         Reject = 0
         If( sig_Clip .gt. 0.0 ) Then
           Thres = Chi * sig_Clip * sig_Clip
           Do J = 1, Points
             Temp = Y( J ) - FUNCTN( X( J ), A, n_Par )
             Temp = Temp * Temp * W( J )
             If( Temp .gt. Thres ) Then
              Reject = Reject + 1
              W( J ) = 0.0
             End If
           End Do
         End If
         If( Reject .gt. 0 ) Then
           Write( *, 550 ) Reject
550           Format( ' * ', I3, ' points are rejected as outliers.' / &
                          '   Restarting fit again...' )
         End If
       End Do
       ! Final call...calculate covariances
       if (debug) print*, '@ Calculating co-variances.'
       Lambda = 0.0
       Call MARQUARDT( X, Y, W, Points, Lambda, A, dA, Covar,&
                             Chi, n_Par, n_Term, FUNCTN, FDERIV, flag_mrq )
       if (debug) print*, '@ done'
       if (flag == 0) flag = flag_mrq
       Do K = 1, n_Par
         If( list_Cyclic( K ) .ne. 0 ) Then       ! Cyclic
           If( A( K ) .lt. tab_Bottom( K ) ) Then
             A( K ) = tab_Bottom( K ) + tab_Period( K ) &
                           + Mod( A( K ) - tab_Bottom( K ), tab_Period( K ) )
           Else
             A( K ) = tab_Bottom( K ) &
                          + Mod( A( K ) - tab_Bottom( K ), tab_Period( K ) )
           End If
         End If
       End Do

       End Subroutine CURFIT



       Integer Function CYC_DAT( Cycle )

       Real Cycle

       If( Cycle .le. 0.0 ) Then
         Write( *, 100 )
100         Format( ' ERROR:: unacceptable periodicity' )
         cyclic_Data = .False.
         cycle_Data = 0.0
         CYC_DAT = -1
       Else
         cyclic_Data = .True.
         cycle_Data = Cycle
         CYC_DAT = 1
       End If

       End Function CYC_DAT



       Integer Function CYC_PAR( N, Lo, Hi )

       ! Select the relevant parameter
       Integer N
       ! Specify the range for cyclic behaviour
       Real Lo, Hi       

       Integer K

       If( Hi .le. Lo ) Then
         Write( *, 100 )
100         Format( ' ERROR:: Range of cyclic behaviour is not positve' )
         CYC_PAR = -1
       Else If( N .le. 0 .or. N .gt. m_Par ) Then
         Write( *, 110 )
110         Format( ' ERROR:: A non-existent parameter has bee specified' )
         CYC_PAR = -1
       Else If( list_Cyclic( N ) .ne. 0 ) Then
         Write( *, 120 ) N
120         Format( ' ERROR:: Parameter ', I2, ' is already declared cyclic' )
         CYC_PAR = -1
       Else If( l_Index( N ) .ne. N ) Then
         Write( *, 130 ) N
130         Format( ' ERROR:: Parameter ', I2, ' is linked' )
         CYC_PAR = -1
       Else If( lo_A( N ) .ne. n_Large .or. hi_A( N ) .ne. p_Large ) Then
         If( hi_A( N ) - lo_A( N ) .ge. Hi - Lo ) Then
           Write( *, 140 )
140           Format( ' ERROR:: Conflict between LIMIT_PAR and CYC_PAR' )
           CYC_PAR = -1
         Else
           Do K = 1, m_Par
             If( l_Index( K ) .eq. -N ) Then
              Write( *, 150 )
150               Format( ' ERROR:: Illegal attempt ', &
                                         'to use CYC_PAR on a LINKed parameter' )
              CYC_PAR = -1
              Return
             Else If( l_Index( K ) .eq. N .and. K .ne. N ) Then
              list_Cyclic( K ) = -1
              tab_Bottom( K ) = Lo
              tab_Period( K ) = Hi - Lo
             End If
           End Do
           Write( *, 160 )
160           Format( ' WARNING:: This parameter has already been LIMITed' )
           list_Cyclic( N ) = -1
           tab_Bottom( N ) = Lo
           tab_Period( N ) = Hi - Lo
           CYC_PAR = 0
         End If
       Else
         Do K = 1, m_Par
           If( l_Index( K ) .eq. -N ) Then
             Write( *, 160 )
             CYC_PAR = -1
             Return
           Else If( l_Index( K ) .eq. N .and. K .ne. N ) Then
             list_Cyclic( K ) = -1
             tab_Bottom( K ) = Lo
             tab_Period( K ) = Hi - Lo
             lo_A( K ) = Lo - tab_Period( K ) * 0.25 + l_Offset( K )
             hi_A( K ) = Lo + tab_Period( K ) * 1.25 + l_Offset( K )
           End If
         End Do
         list_Cyclic( N ) = 1
         tab_Bottom( N ) = Lo
         tab_Period( N ) = Hi - Lo
         lo_A( N ) = Lo - tab_Period( N ) * 0.25
         hi_A( N ) = Lo + tab_Period( N ) * 1.25
         CYC_PAR = 1
       End If

       End Function CYC_PAR



       Subroutine D_SCALE( Factor )

       Real Factor

       If( Factor .le. 0.0 ) Then
         tame_Factor = 1.0
       Else
         tame_Factor = Factor
       End If

       End Subroutine D_SCALE


       Integer Function FIX_PAR( N, n_Par, n_Term )

       ! Used to fix parameter #N throughout a fit
       !                          28 July 1987, KM

       ! Slave parameter
       Integer N
       ! Total number of parameters
       Integer n_Par
       ! Number of free parameters
       Integer n_Term              

       Integer I, K

       ! First, make sure it's not a linked parameter
       If( l_Index( N ) .ne. N ) Then
         Write( *, 100 ) N
100         Format( ' ERROR:: Parameter #', I2, ' is linked!' / &
                          '           This call to FIX_PAR is ignored' )
         FIX_PAR = -1              ! Error status
         Return
       End If
       Do K = 1, n_Par
         If( Abs( l_Index( K ) ) .eq. N .and. K .ne. N ) Then
           Write( *, 100 ) N
           FIX_PAR = -1       ! Error status
           Return
         End If
       End Do
!              Is it already fixed?
       If( Table_Free( N ) .eq. 1 ) Then
         Write( *, 200 ) N
200         Format( ' WARNING:: Parameter #', I2, ' is already fixed!' / &
                          '           This call to FIX_PAR is ignored' )
         FIX_PAR = 0              ! Warning status
         Return
       End If
       Table_Free( N ) = 1
       I = 1
       Do While( I .le. n_Term .and. Index( I ) .ne. N )
         I = I + 1
       End Do
       n_Term = n_Term - 1       ! One less free parameter
       Do K = I, n_Term
         Index( K ) = Index( K + 1 )
       End Do

       FIX_PAR = 1              ! OK status

       End Function FIX_PAR



       Integer Function FREE_PAR( N, n_Par, n_Term )

       ! Parameter N (>0) becomes free again!
       ! Or if N=0, all parameters will become free

       ! Fixed ---> Free
       Integer N
       ! Number of parameters
       Integer n_Par
       ! Number of free parameters
       Integer n_Term       

       Integer K

       If( N .ge. 1 .and. N .le. n_Par ) Then       ! Single mode
         If( l_Index( N ) .ne. N ) Then
           Write( *, 100 ) N
100           Format( ' ERROR:: Parameter #', I2, ' is linked!' /  &
                           '           Use SEPARATE_PAR to separate' )
           FREE_PAR = -1
           Return
         Else If( Table_Free( N ) .eq. 0 ) Then
           Write( *, 200 ) N
200           Format( ' WARNING:: Parameter #', I2, ' is free!' /  &
                   '           This call to FREE_PAR is ignored' )
           FREE_PAR = 0
           Return
         End If
         n_Term = n_Term + 1       ! One more free parameter
         Table_Free( N ) = 0
         Index( n_Term ) = N       ! N must now become free
         FREE_PAR = 1
       Else If( N .eq. 0 ) Then              ! All-clear mode
         Do K = 1, n_Par
           If( Table_Free( K ) .eq. 1 ) Then
             n_Term = n_Term + 1
             Table_Free( K ) = 0
             Index( n_Term ) = K
           End If
         End Do
         FREE_PAR = 1
       Else                                   ! Uninown mode
         Write( *, 300 )
300         Format( ' ERROR:: invalid argument!',  &
                           ' This call to FREE_PAR is ignored' )
         FREE_PAR = -1
       End If

       End Function FREE_PAR



       Subroutine GAUSS_JORDAN( Matrix, Order, Vector, Flag )

!       Solves the linear equation Matrix.Answer=Vector
!       and replaces Vector by Answer, and Matrix by its inverse.
!
!                                          7/July/1987, KM
!                     Modified a bit              August 1 1989, KM


       Integer, intent(in) :: Order                     ! Used size
       Real, dimension(:,:), intent(inout) :: matrix    ! Matrix of interest
       Real, dimension(:), intent(inout) :: vector        ! Vector of interest
       Integer, intent(out) :: Flag                     ! -1 if singular

       Integer Pivot( Order )       ! Bookkeeping
       Integer r_Index( Order )       ! variables for
       Integer c_Index( Order )       ! pivoting

       Real Big, inv_Pivot, Work
       Integer Row, Col, I, J, K, L, M

       Flag = 0
       
       Pivot=0

       ! Main loop over the columns to be reduced
       Do I = 1, Order
         Big = 0.0
         ! Outer loop to search for a pivot element
         Do J = 1, Order
           If( Pivot( J ) .ne. 1 ) Then
             Do K = 1, Order
              If( Pivot( K ) .eq. 0 ) Then
                If( Abs( Matrix( J, K ) ) .ge. Big ) Then
                  Big = Abs( Matrix( J, K ) )
                  Row = J
                  Col = K
                End If
              Else If( Pivot( K ) .gt. 1 ) Then
                Flag = -1
                Return
              End If
             End Do
           End If
         End Do
         Pivot( Col ) = Pivot( Col ) + 1

!         We now have the pivot element, so exchange rows to put the
!         pivot element on the diagonal.  Changes are labeled using
!         *_Index arrays (not actually interchanged).

         If( Row .ne. Col ) Then
           Do L = 1, Order
             Work = Matrix( Row, L )
             Matrix( Row, L ) = Matrix( Col, L )
             Matrix( Col, L ) = Work
           End Do
           Work = Vector( Row )
           Vector( Row ) = Vector( Col )
           Vector( Col ) = Work
         End If
         r_Index( I ) = Row       ! Now ready to divide the pivot row
         c_Index( I ) = Col       ! by the pivot element
         If( Matrix( Col, Col ) <= tiny(matrix(col,col))) Then       
           ! New - protect.  
           ! Used to have a number in there.
           Flag = -1
           Return
         End If
         inv_Pivot = 1.0 / Matrix( Col, Col )
         Matrix( Col, Col ) = 1.0
         Do L = 1, Order
           Matrix( Col, L ) = Matrix( Col, L ) * inv_Pivot
         End Do
         Vector( Col ) = Vector( Col ) * inv_Pivot
                             ! Next, reduce the rows
         Do M = 1, Order
           If( M .ne. Col ) Then       ! ...except for the pivot one
             Work = Matrix( M, Col )
             Matrix( M, Col ) = 0.0
             Do L = 1, Order
              Matrix( M, L ) &
                           = Matrix( M, L ) - Matrix( Col, L ) * Work
             End Do
             Vector( M ) = Vector( M ) - Vector( Col ) * Work
           End If
         End Do
       End Do              ! Endof the main loop over columns
                      ! Now unscramble the solution
       Do L = Order, 1, -1
         If( r_Index( L ) .ne. c_Index( L ) ) Then
           Do K = 1, Order
             Work = Matrix( K, r_Index( L ) )
             Matrix( K, r_Index( L ) ) = Matrix( K, c_Index( L ) )
             Matrix( K, c_Index( L ) ) = Work
           End Do
         End If
       End Do

       End Subroutine GAUSS_JORDAN




       Integer Function LIMIT_PAR( N, Lo, Hi )

       ! Set up the limit within which parameter A( N ) must change
       ! It does NOT check if the current value of A( N ) is within the bound.
       ! The default values are +/-1.0E+36.

       Integer N       ! Parameter of interest
       Real Lo              ! Lower boundary
       Real Hi              ! Upper boundary

       Integer K

       If( Lo .ge. Hi ) Then
         Write( *, 100 )
100         Format( ' ERROR:: Upper limit is smaller than lower limit' / &
                           '           This call to LIMIT_PAR is ignored' )
         LIMIT_PAR = -1
       Else If( N .ge. 1 .and. N .le. m_Par ) Then
         If( list_Cyclic( N ) .ne. 0 ) Then       ! Cyclic
           If( l_Index( N ) .ne. N ) Then
             Write( *, 105 )
105             Format( ' ERROR:: Sorry, you cannot do that' )
             LIMIT_PAR = -1
             Return
           Else If( Lo .ge. tab_Bottom( N ) &
                           .and. Hi .lt. tab_Bottom( N ) + tab_Period( N ) ) Then
             Write( *, 110 )
110             Format( ' WARNING:: LIMITing a cyclic parameter' )
             lo_A( N ) = Lo
             hi_A( N ) = Hi
             list_Cyclic( N ) = -1
             LIMIT_PAR = 0
           Else
             If( Hi - Lo .lt. tab_Period( N ) ) Then
              Write( *, 110 )
              Lo_A( N ) = Lo
              hi_A( N ) = Hi
              list_Cyclic( N ) = -1
               LIMIT_PAR = 0
             Else
              Write( *, 120 )
120              Format( ' ERROR:: Illegal attempt to LIMIT a cyclic parameter' )
              LIMIT_PAR = -1
             End If
           End If
           Do K = 1, N
             If( l_Index( K ) .eq. N .and. K .ne. N ) Then
              lo_A( K ) = lo_A( N ) + l_Offset( K )
              hi_A( K ) = hi_A( N ) + l_Offset( K )
             End If
           End Do
         Else
           Do K = 1, N
             If( l_Index( K ) .eq. -N ) Then
              lo_A( K ) = lo_A( N ) * l_Offset( K )
              hi_A( K ) = hi_A( N ) * l_Offset( K )
             Else If( l_Index( K ) .eq. N .and. K .ne. N ) Then
              lo_A( K ) = lo_A( N ) + l_Offset( K )
              hi_A( K ) = hi_A( N ) + l_Offset( K )
             End If
           End Do
           lo_A( N ) = Lo
           hi_A( N ) = Hi
           LIMIT_PAR = 1
         End If
       Else
         Write( *, 200 ) N
200         Format( ' ERROR:: No such parameter #', I2 / &
                           '           This call to LIMIT_PAR is ignored' )
         LIMIT_PAR = -1
       End If

       End Function LIMIT_PAR




       Integer Function NOLIMIT_PAR( N )

       Integer N       ! Parameter of interest

       Integer K

       If( N .ge. 1 .and. N .le. m_Par ) Then
         If( list_Cyclic( N ) .eq. 1 ) Then
           Write( *, 110 ) N
110           Format( ' ERROR:: No limit has been set for cyclic parameter', I3 )
           NOLIMIT_PAR = -1
         Else If( list_Cyclic( N ) .eq. 0 ) Then
           lo_A( N ) = n_Large
           hi_A( N ) = p_Large
           Do K = 1, m_Par
             If( l_Index( K ) .eq. N .and. K .ne. N ) Then
              lo_A( K ) = lo_A( N ) + l_Offset( K )
              hi_A( K ) = hi_A( N ) + l_Offset( K )
             Else If( l_Index( K ) .eq. -N ) Then
              lo_A( K ) = lo_A( N ) * l_Offset( K )
              hi_A( K ) = hi_A( N ) * l_Offset( K )
             End If
           End Do
           NOLIMIT_PAR = 1
         Else              ! Limited cyclic parameter
           If( l_Index( N ) .eq. N ) Then
             lo_A( N ) = tab_Bottom( N ) - tab_Period( N ) * 0.25
             hi_A( N ) = tab_Bottom( N ) + tab_Period( N ) * 1.25
             Do K = 1, m_Par
              If( l_Index( K ) .eq. N .and. K .ne. N ) Then
                lo_A( K ) = lo_A( N ) + l_Offset( K )
                hi_A( K ) = hi_A( N ) + l_Offset( K )
              End If
             End Do
           Else
             Write( *, 120 )
120             Format( ' ERROR:: Illegal attempt on a cyclic slave parameter' )
             NOLIMIT_PAR = -1
             Return
           End If
           list_Cyclic( N ) = 1
           NOLIMIT_PAR = 1
         End If
       Else If( N .eq. 0 ) Then
         NOLIMIT_PAR = 0
         Do K = 1, m_Par
           If( list_Cyclic( K ) .eq. 1 ) Then
             Write( *, 110 ) K
             NOLIMIT_PAR = -1
           Else If( list_Cyclic( K ) .eq. 0 ) Then
             If( l_Index( K ) .eq. K ) Then
              lo_A( K ) = n_Large
              hi_A( K ) = p_Large
             Else If( l_Index( K ) .gt. 0 ) Then
              lo_A( K ) = lo_A( l_Index( K ) ) + l_Offset( K )
              hi_A( K ) = hi_A( l_Index( K ) ) + l_Offset( K )
             Else
              lo_A( K ) = lo_A( -l_Index( K ) ) * l_Offset( K )
              hi_A( K ) = hi_A( -l_Index( K ) ) * l_Offset( K )
             End If
           Else
             If( l_Index( K ) .eq. K ) Then
              lo_A( K ) = tab_Bottom( K ) - tab_Period( K ) * 0.25
              hi_A( K ) = tab_Bottom( K ) + tab_Period( K ) * 1.25
             Else
              lo_A( K ) = tab_Bottom( l_Index( K ) ) &
                           - tab_Period( l_Index( K ) ) * 0.25 + l_Offset( K )
              hi_A( K ) = tab_Bottom( l_Index( K ) ) &
                           + tab_Period( l_Index( K ) ) * 1.25 + l_Offset( K )
             End If
             list_Cyclic( K ) = 1
           End If
         End Do
       Else
         Write( *, 200 ) N
200         Format( ' ERROR:: No such parameter #', I3 / &
                          '         This call to NOLIMIT_PAR is ignored' )
         NOLIMIT_PAR = -1
       End If

       End Function NOLIMIT_PAR



       Integer Function LINK_PAR( Master, Slave, Cnst, Mode, A, n_Par, n_Term )

       ! Used to `link' parameter #N to parameter #Master such that
       ! A( Slave ) = A( Master ) + Cnst (Mode = 0)
       ! or A( Slave ) = Cnst * A( Master ) (Mode = 1)
       ! at any time during the fit
       ! e.g. fiiting a resolved doublet whose separation is known.
       ! Both Master and Slave must be on the list of free parameters (Index)
       ! prior to the call.  Designed to be called as the final part
       ! of parameter initialization.
       !                                   12 August 1987, KM

       Integer Master              ! Master parameter
       Integer Slave              ! Slave parameter
       Real Cnst              ! Difference or Ratio
       Integer Mode              ! Selects difference or ratio
       Integer n_Par              ! Total number of parameters
       Real A( n_Par )              ! Parameters
       Integer n_Term              ! Number of free parameters

       Integer J, K

       ! First, look for Master and Slave in Index
       LINK_PAR = 1
       If( Master .le. 0 .or. Master .gt. n_Par ) Then
         Write( *, 110 ) Master
110         Format( ' ERROR::', I3, ' is out of range; ', &
                                  'this call to LINK_PAR is ignored' )
         LINK_PAR = -1
         Return
       Else If( Slave .le. 0 .or. Slave .gt. n_Par ) Then
         Write( *, 110 ) Slave
         LINK_PAR = -1
         Return
       Else If( Table_Free( Master ) .eq. 1 ) Then
         Write( *, 120 ) Master
120         Format( ' ERROR:: Parameter #', I3, ' is fixed' / &
                           '         This call to LINK_PAR is ignored' )
         LINK_PAR = -1
         Return
       Else If( Table_Free( Slave ) .eq. 1 ) Then
         Write( *, 120 ) Slave
         LINK_PAR = -1
         Return
       Else If( Table_Free( Master ) .eq. -1 ) Then
         If( Mode .eq. 0 ) Then
           Write( *, 130 ) Master, l_Index( Master )
130           Format( ' WARNING:: Parameter #', I3, ' is linked to', &
                           ' parameter #', I3 )
           Master = l_Index( Master )
           Cnst = Cnst + l_Offset( Master )
           LINK_PAR = 0
         Else
           Write( *, 140 )
140           Format( ' ERROR:: Link mode mismatch' / &
                           '         This call to LINK_PAR is ignored' )
           LINK_PAR = -1
           Return
         End If
       Else If( Table_Free( Master ) .eq. -2 ) Then
         If( Mode .eq. 1 ) Then
           Write( *, 130 ) Master, -l_Index( Master )
           Master = -l_Index( Master )
           Cnst = Cnst * l_Offset( Master )
           LINK_PAR = 0
         Else
           Write( *, 140 )
           LINK_PAR = -1
           Return
         End If
       Else If( Table_Free( Slave ) .lt. 0 ) Then
         Write( *, 150 ) Slave
150         Format( ' ERROR:: Parameter #', I3, ' is already linked' / &
                           '         This call to LINK_PAR is ignored' )
         LINK_PAR = -1
         Return
       Else If( list_Cyclic( Master ) .ne. 0 ) Then
         If( list_Cyclic( Slave ) .eq. 0 ) Then
           Write( *, 160 )
160           Format( ' ERROR:: Illegal attempt to link ', &
                                  'a cyclic parameter with a non-cyclic one' )
           LINK_PAR = -1
           Return
         Else If( list_Cyclic( Slave ) .eq. 1 ) Then
           If( Mode .eq. 1 ) Then
             Write( *, 170 )
170             Format( ' ERROR:: Illegal attempt ', &
                                  'to fix the ratio of cyclic parameters' )
             LINK_PAR = -1
             Return
           End If
           lo_A( Slave ) = lo_A( Master ) + Cnst
           hi_A( Slave ) = hi_A( MAster ) + Cnst
           list_Cyclic( Slave ) = -1
         Else
           Write( *, 180 )
180           Format( ' ERROR:: Sorry, you cannot do that' )
           LINK_PAR = -1
           Return
         End If
       Else If( list_Cyclic( Slave ) .ne. 0 ) Then
         Write( *, 160 )
         LINK_PAR = -1
         Return
       End If

       If( Mode .eq. 0 ) Then
         l_Index( Slave ) = Master
         l_Offset( Slave ) = Cnst
         A( Slave ) = A( Master ) + Cnst
         Table_Free( Slave ) = -1
       Else If( Mode .eq. 1 ) Then
         l_Index( Slave ) = -Master
         l_Offset( Slave ) = Cnst
         A( Slave ) = A( Master ) * Cnst
         Table_Free( Slave ) = -2
       End If
       n_Term = n_Term - 1       ! One less free parameter
       J = 1
       Do While( J .le. n_Term .and. Index( J ) .ne. Slave )
         J = J + 1
       End Do
       Do K = J, n_Term
         Index( K ) = Index( K + 1 )       ! Adjust Index
       End Do

       End Function LINK_PAR




       Integer Function SEPARATE_PAR( Slave, n_Par, n_Term )

       ! Parameter Slave (>0) becomes independent again!
       ! Or if Slave=0, all parameters will become independent

       Integer Slave       ! Slave ---> Free Man
       Integer n_Par       ! Number of parameters
       Integer n_Term       ! Number of free parameters

       Integer K

       SEPARATE_PAR = 1
       If( Slave .ge. 1 .and. Slave .le. n_Par ) Then       ! Single mode
         If( l_Index( Slave ) .eq. Slave ) Then
           Write( *, 100 ) Slave
100           Format( ' ERROR:: Parameter #', I2, ' is not linked!'/ &
                                  '         Call to SEPARATE_PAR ignored' )
           SEPARATE_PAR = -1
           Return
         End If
         l_Index( Slave ) = Slave       ! Now un-linked
         n_Term = n_Term + 1              ! One more free parameter
         Index( n_Term ) = Slave       ! Slave must now become free
         Table_Free( Slave ) = 0
         If( list_Cyclic( Slave ) .eq. -1 ) Then
           lo_A( Slave ) = tab_Bottom( Slave ) - tab_Period( Slave ) * 0.25
           hi_A( Slave ) = tab_Bottom( Slave ) + tab_Period( Slave ) * 1.25
           list_Cyclic( Slave ) = 1
         End If
       Else If( Slave .eq. 0 ) Then              ! All-clear mode
         Do K = 1, n_Par
           If( l_Index( K ) .ne. K ) Then       ! Found a linked par
             l_Index( K ) = K
             n_Term = n_Term + 1
             Index( n_Term ) = K
             Table_Free( K ) = 0
             If( list_Cyclic( K ) .eq. -1 ) Then
              lo_A( K ) = tab_Bottom( K ) - tab_Period( K ) * 0.25
              hi_A( K ) = tab_Bottom( K ) + tab_Period( K ) * 1.25
              list_Cyclic( K ) = 1
             End If
           End If
         End Do
       Else                                   ! Unknown mode
         Write( *, 200 )
200         Format( ' ERROR:: invalid argument!', &
                                  ' This call to SEPARATE_PAR ignored' )
         SEPARATE_PAR = -1
       End If

       End Function SEPARATE_PAR



       Subroutine MARQUARDT( X, Y, W, Points, Lambda, A, dA, &
         Covar, Chi, n_Par, n_Term, FUNCTN, FDERIV, flag_marq )


       ! Number of data points
       Integer, intent(in) :: Points
       ! Data and weights
       Real, dimension(Points), intent(in) :: X, Y, W
       ! Fit control variable
       Real, intent(inout) :: Lambda
       ! # of params of model
       Integer, intent(in) :: n_Par
       ! Parameters and increments
       Real, dimension(n_Par), intent(inout) :: A
       Real, dimension(n_Par), intent(in) :: dA
       ! Covariance matrix
       Real, dimension(:,:), intent(out) :: Covar
       ! Chi squared
       Real, intent(out) :: Chi
       ! # of free paras
       Integer, intent(in) :: n_Term
       ! Model Function and its derivative
       Real, external :: FUNCTN              
       External FDERIV
       ! A flag, normally zero.  Inherited from svd_comp.
       integer, intent(out) :: flag_marq                    

       Real, save :: Alpha( m_Par, m_Par )       ! Curvature matrix
       Real, save :: try_A( m_Par )              ! Trial paras
       Real, save :: Beta( m_Par )
       Real, save :: incr_A( m_Par )
       Real, save :: svd_U( m_Par, m_Par )
       Real, save :: svd_W( m_Par )
       Real, save :: svd_V( m_Par, m_Par )
       Integer, save :: I, J, K, Flag
       Real, save :: o_Chi

       flag_marq=0

       If( Lambda .lt. 0.0 ) Then
         ! First call: Initialize and return
         Lambda = 0.001       ! Initialize Lambda
         Call MRQ_COF( X, Y, W, Points, A, dA, n_Par, n_Term, &
                           Alpha, Beta, Chi, FUNCTN, FDERIV )
         o_Chi = Chi
         Do K = 1, n_Par
           try_A( K ) = A( K )
         End Do
         Return
       End If
       ! Change diagonal elements
       Do I = 1, n_Term
         Do J = 1, n_Term
           Covar( I, J ) = Alpha( I, J )
         End Do
         Covar( I, I ) = Alpha( I, I ) * ( 1.0 + Lambda )
         incr_A( I ) = Beta( I )
       End Do
       Call GAUSS_JORDAN( Covar, n_Term, incr_A, Flag )

       If( Lambda .eq. 0.0 ) Then       ! Final call...error calc only
         If( Flag .ne. 0 ) Then
           Do I = 1, n_Term              ! Recalculate COVAR
             Do J = 1, n_Term
              svd_U( I, J ) = Alpha( I, J )
             End Do
             svd_U( I, I ) = Alpha( I, I ) * ( 1.0 + Lambda )
             incr_A( I ) = Beta( I )
           End Do
           If( Debug ) &
           Print *, '@ Singular Matrix trying Singular Value Decomposition.'
           ! You may like to try svdcmp instead.
           Call SVD_COMP( svd_U, n_Term, n_Term, svd_W, &
           svd_V, flag_marq )
           if (flag_marq /= 0) return
           Call SVD_REC( svd_U, svd_W, svd_V, n_Term, n_Term )
           Do I = 1, n_Term
             Do J = 1, n_Term
              Covar( I, J ) = svd_U( I, J )
             End Do
           End Do
         End If
         Call COV_SORT( Covar, n_Par, n_Term )
         Return
       End If

       If( Flag .ne. 0 ) Then              ! Singular Matrix
         If( Debug ) Then
           Print *, '@ Singular Matrix...'
           Do I = 1, n_Term
             Print *, ( Alpha( I, J ), J = 1, n_Term )
           End Do
           Print *, ' ...trying Singular Value Decomposition'
         End If
         Do I = 1, n_Term               ! Recalculate COVAR
           Do J = 1, n_Term
             svd_U( I, J ) = Alpha( I, J )
           End Do
           svd_U( I, I ) = Alpha( I, I ) * ( 1.0 + Lambda )
           incr_A( I ) = Beta( I )
         End Do
         ! You may like to try svdcmp instead.
         Call SVD_COMP( svd_U, n_Term, n_Term, svd_W, &
         svd_V, flag_marq )
         if (flag_marq /= 0) return
         Call SV_BCSB ( svd_U, svd_W, svd_V, n_Term, n_Term, incr_A )
       End If

       Do I = 1, n_Term
         incr_A( I ) = incr_A( I ) * tame_Factor
       End Do

       Call UPDATE_PARAM( try_A, incr_A, n_Par, n_Term )

       Call MRQ_COF( X, Y, W, Points, try_A, dA, n_Par, n_Term, &
                          Covar, incr_A, Chi, FUNCTN, FDERIV )

       ! This line used to read...
       ! If( Chi .lt. o_Chi ) Then       
       ! But is clearly better done like this.
       if (o_chi - chi > -1.0*epsilon(chi)) then
         ! Success
         Lambda = Lambda * 0.1
         o_Chi = Chi
         Do I = 1, n_Term
           Do J = 1, n_Term
             Alpha( I, J ) = Covar( I, J )
           End Do
           Beta( I ) = incr_A( I )
         End Do
         Do J = 1, n_Par
           A( J ) = try_A( J )
         End Do
       Else                            ! Failure
         Lambda = Lambda * 100.0
         Chi = o_Chi
         Do I = 1, n_Term
           J = Index( I )
           try_A( J ) = A( J )
         End Do
       End If

       End Subroutine MARQUARDT



       Subroutine MRQ_COF( X, Y, W, Points, A, dA, n_Par, n_Term, &
                                 Alpha, Beta, Chi, FUNCTN, FDERIV )

       ! Calculates Alpha, Beta and Chi.
       !                 8 July 1987, KM

       ! Parameters
       Integer, intent(in) :: Points
       Real, dimension(Points), intent(in) :: X, Y, W
       Integer, intent(in) :: n_Par
       Real, dimension(n_par), intent(in) :: A, dA
       Integer, intent(in) :: n_Term
       Real, dimension(:,:), intent(out) :: Alpha 
       Real, dimension(n_Par), intent(out) :: Beta
       Real, intent(out) :: Chi
       Real, external :: FUNCTN
       External FDERIV

       ! Locals
       Real, allocatable, dimension(:) :: dF_dA
       Real :: dY, Work
       Integer :: I, J, K
       
       ! Initialize (symmetric) Alpha, Beta
       Do I = 1, n_Term
         Do J = 1, I
           Alpha( I, J ) = 0.0
         End Do
         Beta( I ) = 0.0
       End Do
       Chi = 0.0

       ! Create array of differentials.
       allocate (df_da(n_par))
       
       ! Summation loop over all data
       Do K = 1, Points
         If( W( K ) .gt. 0.0 ) Then
           dY = Y( K ) - FUNCTN( X( K ), A, n_Par )
           If( cyclic_Data ) Then
             dY = dY - Nint( dY / cycle_Data ) * cycle_Data
           End If
           Call FDERIV( X( K ), A, dA, n_Par, dF_dA )
           Do I = 1, n_Par
             If( Table_Free( I ) .eq. -1 ) Then
              J = l_Index( I )
              dF_dA( J ) = dF_dA( J ) + dF_dA( I )
             Else If( Table_Free( I ) .eq. -2 ) Then
              J = -l_Index( I )
              dF_dA( J ) = dF_dA( J ) + dF_dA( I ) * l_Offset( I )
             End If
           End Do
           Do I = 1, n_Term
             Work = dF_dA( Index( I ) ) * W( K )
             Do J = 1, I
              Alpha( I, J ) = Alpha( I, J ) + Work * dF_dA( Index( J ) )
             End Do
             Beta( I ) = Beta( I ) + dY * Work
           End Do
           Chi = Chi + dY * dY * W( K )
         End If
       End Do
       ! Finished with df_da
       deallocate(df_da)
       
       ! Fill in the symmetric side
       Do J = 2, n_Term
         Do I = 1, J - 1
           Alpha( I, J ) = Alpha( J, I )
         End Do
       End Do

       End Subroutine MRQ_COF



       Subroutine RESET( )

       Integer K

       Do K = 1, m_Par
         Index( K ) = K
         Table_Free( K ) = 0
         l_Index( K ) = K
         l_Offset( K ) = 0.0
         lo_A( K ) = n_Large
         hi_A( K ) = p_Large
         list_Cyclic( K ) = 0
         tab_Period( K ) = 0.0
         tab_Bottom( K ) = 0.0
       End Do
       cyclic_Data = .False.
       cycle_Data = 0.0
       tame_Factor = 1.0

     End subroutine reset



     Subroutine SVD_COMP( Matrix, size_X, size_Y, W, V, flag_svd )

       ! Singular Value Decomposition.

       ! Given an arbitrary matrix A, computes its singular value decomposition
       ! A = U W Vt; the matrix U replaces A on output.  The diagonal matrix
       ! W is output as a vector, and the matrix V (not its transpose Vt) is
       ! calculated.  size_X must be greater or equal to dim_Y; if it is smaller,
       ! then A should be filled up to square with zero rows.

       ! The Numerical Recipes used to look to see if certain numbers were
       ! close to zero, by adding them to a large number, and looking for
       ! equality with that large number.  Fortran 95 allows you to use the
       ! tiny function.

       ! Matrix to be decomposed / U
       real, dimension(:,:), intent(inout) :: Matrix
       ! Actual size of the matrix
       integer, intent(in) :: size_X, size_Y              
       ! Diagonal matrix
       real, dimension(:), intent(out) :: W
       ! V as in A = U.W.Vt
       real, dimension(:,:), intent(out) :: V
       ! A flag, set to -51 if 30 iterations are tried.
       integer, intent(out) :: flag_svd

       Real RV1( m_Par )
       Real Scale
       Real C, G, H, S, X, Y, Z
       double precision :: F
       Integer I, J, K, L, N, Loop
       real :: anorm

       flag_svd=0

       If( size_X .lt. size_Y ) Stop 'ERROR:: unacceptable size'
       If( size_Y .gt. m_Par ) Stop 'ERROR:: matrix too big'

       ! Householder reduction to bidiagonal form

       G = 0.0
       Scale = 0.0
       anorm=0.0
       Do I = 1, size_Y
         L = I + 1
         RV1( I ) = Scale * G
         G = 0.0
         S = 0.0
         Scale = 0.0
         Do K = I, size_X
           Scale = Scale + Abs( Matrix( K, I ) )
         End Do
         If( abs(Scale) > tiny(Scale) ) Then
           Do K = I, size_X
             Matrix( K, I ) = Matrix( K, I ) / Scale
             S = S + Matrix( K, I ) * Matrix( K, I )
           End Do
           F = Matrix( I, I )
           G = -Sign( Sqrt( S ), real(F) )
           H = F * G - S
           Matrix( I, I ) = F - G
           If( I .ne. size_Y ) Then
             Do J = L, size_Y
              S = 0.0
              Do K = I, size_X
                S = S + Matrix( K, I ) * Matrix( K, J )
              End Do
              F = S / H
              Do K = I, size_X
                Matrix( K, J ) = Matrix( K, J ) + F * Matrix( K, I )
              End Do
             End Do
           End If
           Do K = I, size_X
             Matrix( K, I ) = Scale * Matrix( K, I )
           End Do
         End If
         W( I ) = Scale * G
         G = 0.0
         S = 0.0
         Scale = 0.0
         If( I .ne. size_Y ) Then
           Do K = L, size_Y
             Scale = Scale + Abs( Matrix( I, K ) )
          End Do
           If( abs(Scale) > tiny(Scale) ) Then
             Do K = L, size_Y
               Matrix( I, K ) = Matrix( I, K ) / Scale
               S = S + Matrix( I, K ) * Matrix( I, K )
             End Do
             F = Matrix( I, L )
             G = -Sign( Sqrt( S ), real(F) )
             H = F * G - S
             Matrix( I, L ) = F - G
             Do K = L, size_Y
              RV1( K ) = Matrix( I, K ) / H
             End Do
             If( I .ne. size_X ) Then
              Do J = L, size_X
                S = 0.0
                Do K = L, size_Y
                  S = S + Matrix( J, K ) * Matrix( I, K )
                End Do
                Do K = L, size_Y
                  Matrix( J, K ) = Matrix( J, K ) + S * RV1( K )
                End Do
              End Do
             End If
             Do K = L, size_Y
               Matrix( I, K ) = Scale * Matrix( I, K )
             End Do
           End If
         End If
         anorm=max(anorm,(abs(w(i))+abs(rv1(i))))
       End Do

       ! Accumulation of right-hand transformations
       Do I = size_Y, 1, -1
         If( I .lt. size_Y ) Then
           if (abs(G) > tiny(G))then
             Do J = L, size_Y       
               ! Double division to avoid possible underflow
               V( J, I ) = ( Matrix( I, J ) / Matrix( I, L ) ) / G
             End Do
             Do J = L, size_Y
               S = 0.0
               Do K = L, size_Y
                 S = S + Matrix( I, K ) * V( K, J )
               End Do
               Do K = L, size_Y
                 V( K, J ) = V( K, J ) + S * V( K, I )
               End Do
             End Do
           End If
           Do J = L, size_Y
             V( I, J ) = 0.0
             V( J, I ) = 0.0
           End Do
         End If
         V( I, I ) = 1.0
         G = RV1( I )
         L = I
       End Do

       ! Accumulation of left-hand transformations

       Do I = size_Y, 1, -1
         L = I + 1
         G = W( I )
         If( I .lt. size_Y ) Then
           Do J = L, size_Y
             Matrix( I, J ) = 0.0
           End Do
         End If
         if (abs(G) > tiny(G))then
           G = 1.0 / G
           If( I .ne. size_Y ) Then
              Do J = L, size_Y
                 S = 0.0
              Do K = L, size_X
                S = S + Matrix( K, I ) * Matrix( K, J )
              End Do
              F = ( S / Matrix( I, I ) ) * G
              Do K = I, size_X
                Matrix( K, J ) = Matrix( K, J ) + F * Matrix( K, I )
              End Do
             End Do
           End If
           Do J = I, size_X
             Matrix( J, I ) = Matrix( J, I ) * G
           End Do
         Else
           Do J = I, size_X
             Matrix( J, I ) = 0.0
           End Do
         End If
         Matrix( I, I ) = Matrix( I, I ) + 1.0
       End Do

       ! Diagonalization of the bidiagonal form
       Do K = size_Y, 1, -1
         Do Loop = 1, 30
           Do L = K, 1, -1
             N = L - 1
             if (abs(rv1(l)) < epsilon(anorm))  goto 200
             if (abs(w(n))   < epsilon(anorm))  goto 100
           End Do
100          Continue
           C = 0.0
           S = 1.0
           Do I = L, K
             F = S * RV1( I )
             if (abs(f) > epsilon(anorm)) then
               G = W( I )
               H = Sqrt( F * F + G * G )
               W( I ) = H
               H = 1.0 / H
               C = G * H
               S = -F * H
               Do J = 1, size_X
                 Y = Matrix( J, N )
                 Z = Matrix( J, I )
                 Matrix( J, N ) = Y * C + Z * S
                 Matrix( J, I ) = -Y * S + Z * C
               End Do
             End If
           End Do
200          Continue
           Z = W( K )
           If( L .eq. K ) Then              
             ! Convergence
             If( Z .lt. 0.0 ) Then       
               ! Singular value is made nonnegative
               W( K ) = -Z
               Do J = 1, size_Y
                 V( J, K ) = -V( J, K )
               End Do
             End If
             Goto 300
           End If
           If ( Loop == 30 ) then
             if (debug) print*, '@ ERROR:: No convergence in 30 iterations'
             !flag_svd = -51
             !return
           end if
           ! Shift from bottom 2-by-2 minor
           X = W( L )               
           N = K - 1
           Y = W( N )
           G = RV1( N )
           H = RV1( K )
           F = ( ( Y - Z ) * ( Y + Z ) + ( G - H ) * ( G + H ) ) &
                                                      / ( 2.0 * H * Y )
           if (abs(F) < Sqrt(huge(G)) ) then 
             G = Sqrt( F * F + 1.0 )
           else 
             G=F
           End if 
           F = ( ( X - Z ) * ( X + Z ) &
                          + H * ( Y / ( F + Sign( G, real(F) ) ) - H ) ) / X
           ! Next QR transformation
           C = 1.0
           S = 1.0
           Do J = L, N
             I = J + 1
             G = RV1( I )
             Y = W( I )
             H = S * G
             G = C * G
             Z = Sqrt( F * F + H * H )
             RV1( J ) = Z
             C = F / Z
             S = H / Z
             F = X * C + G * S
             G = -X * S + G * C
             H = Y * S
             Y = Y * C
             Do N = 1, size_Y
              X = V( N, J )
              Z = V( N, I )
              V( N, J ) = X * C + Z * S
              V( N, I ) = -X * S + Z * C
             End Do
             Z = Sqrt( F * F + H * H )
             W( J ) = Z
             ! Rotation can be arbitrary if Z=0
             if (abs(Z) > tiny(Z)) then
               Z = 1.0 / Z
               C = F * Z
               S = H * Z
             End If
             F = C * G + S * Y
             X = -S * G + C * Y
             Do N = 1, size_X
              Y = Matrix( N, J )
              Z = Matrix( N, I )
              Matrix( N, J ) = Y * C + Z * S
              Matrix( N, I ) = -Y * S + Z * C
             End Do
           End Do
           RV1( L ) = 0.0
           RV1( K ) = F
           W( K ) = X
         End Do
300      Continue
       End Do

       End Subroutine SVD_COMP



       Subroutine SV_BCSB( U, W, V, size_X, size_Y, B )

       ! U, W, V as returned by SVD_COMP, B is A.X on input, X on output.
       real, dimension(:,:) :: U, V
       real, dimension(:) :: W, B
       Integer size_X, size_Y       ! Actual sizes used

       Real Sum, Work( m_Par ), max_W, min_W
       Integer I, J, K
!                            ZERO the small W's
       max_W = W( 1 )
       Do J = 2, size_Y
         max_W = Max( max_W, W( J ) )
       End Do
       min_W = max_W * 1.0E-06
       Do J = 1, size_Y
         If( W( J ) .lt. min_W ) W( J ) = 0.0
       End Do

       Do J = 1, size_Y
         Sum = 0.0
         If( W( J ) .ne. 0.0 ) Then
           Do I = 1, size_X
             Sum = Sum + U( I, J ) * B( I )
           End Do
           Sum = Sum /  W( J )
         End If
         Work( J ) = Sum
       End Do

       Do J = 1, size_Y
         Sum = 0.0
         Do K = 1, size_Y
           Sum = Sum + V( J, K ) * Work( K )
         End Do
         B( J ) = Sum
       End Do

       End Subroutine SV_BCSB



      Subroutine SVD_REC( U, W, V, size_X, size_Y )

       ! U, W, V as returned by SVD_COMP
       ! U is replaced by the product Ut.W-1.Vt
       real, dimension(:,:) :: U, V       
       real, dimension(:) :: W
       Integer size_X, size_Y       ! Actual sizes used
       Real Work( m_Par, m_Par ), max_W, min_W
       Integer I, J, K
       
       max_W = maxval(W(1:size_Y))
       ! Handle the small W's
       min_W = max(max_W * 1.0E-06,tiny(max_W))
       Do J = 1, size_Y
         W( J ) = 1.0 / Max( W( J ), min_W )
       End Do

       Do I = 1, size_X
         Do J = 1, size_Y
           Work( I, J ) = 0.0
           Do K = 1, size_Y
             Work( I, J ) = Work( I, J ) + V( J, K ) * W( K ) * U( I, K )
           End Do
         End Do
       End Do

       Do I = 1, size_X
         Do J = 1, size_Y
           U( I, J ) = Work( I, J )
         End Do
       End Do

!       now answer(1...dimy)=sum(u(i,j)*question(i)),i=1...dimx

       End Subroutine SVD_REC



       Subroutine UPDATE_PARAM( try_A, incr_A, n_Par, n_Term )

       ! Update trial parameters within FIT2D. Checks for boundaries, and
       ! scales incr_A down if necessary to keep try_A within a sensible range.
       !                                   13 July 1987, KM
       !                            revised 23 Dec  1987

       ! ARGUMENTS
       Integer n_Par              ! Number of parameters
       Integer n_Term              ! Number of free parameters
       Real try_A( n_Par )       ! The current parameter values
       Real incr_A( n_Par )       ! To be changed by...

       Real Ratio, m_Ratio, Scale
       Real Labour
       Integer J, K, L
       Integer r_Index( m_Par )

       ! Is any parameter going off-boundary?
       m_Ratio = 0.0
       Do K = 1, n_Term
         ! Un-linked parameters first
         J = Index( K )
         r_Index( J ) = K
         If( incr_A( K ) .gt. 0.0 ) Then
           If( hi_A( J ) .le. try_A( J ) ) Then
             If( Debug ) Then
              Write( *, 99 ) J
 99              Format( ' @ Parameter #', I2, ' was about to go off-limit' / &
                           '    NOT UPDATING PARAMETERS DURING THIS ITERATION' )
             End If
             Return
           End If
           Ratio = incr_A( K ) / ( hi_A( J ) - try_A( J ) )
         Else If( incr_A( K ) .eq. 0.0 ) Then
           Ratio = 0.0
         Else
           If( lo_A( J ) .ge. try_A( J ) ) Then
             If( Debug ) Then
              Write( *, 99 ) J
             End If
             Return
           End If
           Ratio = incr_A( K ) / ( lo_A( J ) - try_A( J ) )
         End If
         If( Ratio .ge. 1.0 ) Then
           m_Ratio = Max( m_Ratio, Ratio )
           If( Debug ) Then
             Write( *, 100 ) J
100             Format( ' @ Parameter #', I2, ' was about to go off-limit' )
           End If
         End If
       End Do
       Do K = 1, n_Par                            ! Linked parameters next
         If( Table_Free( K ) .lt. 0 ) Then       ! Linked
           If( Table_Free( K ) .eq. -1 ) Then       ! Mode = 0
             Labour = incr_A( r_Index( l_Index( K ) ) )
           Else                            ! Mode = 1
             Labour = incr_A( r_Index( -l_Index( K ) ) ) * l_Offset( K )
           End If
           If( Labour .gt. 0.0 ) Then
             Ratio = Labour / ( hi_A( K ) - try_A( K ) )
           Else If( Labour .eq. 0.0 ) Then
             Ratio = 0.0
           Else
             Ratio = Labour / ( lo_A( K ) - try_A( K ) )
           End If
           If( Ratio .ge. 1.0 ) Then
             If( Debug ) Then
              Write( *, 100 ) K
             End If
             m_Ratio = Max( m_Ratio, Ratio )
           End If
         End If
       End Do
                             ! Rescale incr_A if necessary
       If( m_Ratio .ge. 1.0 ) Then
         Scale = 0.95 / m_Ratio
         If( Debug ) Then
           Write( *, 110 ) Scale
110           Format( ' @ Changes in trial parameters are scaled down by ', &
                                                              1PE10.3 )
         End If
       Else
         Scale = 1.0
       End If
                             ! Now calculate the new trial parameters
       Do K = 1, n_Term
         J = Index( K )
         try_A( J ) = try_A( J ) + incr_A( K ) * Scale
       End Do
                      ! Linked changes
       Do K = 1, n_Par
         If( Table_Free( K ) .eq. -1 ) Then
           L = l_Index( K )
           try_A( K ) = try_A( L ) + l_Offset( K )
         Else  If( Table_Free( K ) .eq. -2 ) Then
           L = -l_Index( K )
           try_A( K ) = try_A( L ) * l_Offset( K )
         End If
       End Do
                      ! Actively cyclic
       Do K = 1, n_Par
         If( list_Cyclic( K ) .eq. 1 ) Then
           If( try_A( K ) .lt. tab_Bottom( K ) ) Then
             try_A( K ) = tab_Bottom( K ) + tab_Period( K ) &
                           + Mod( try_A( K ) - tab_Bottom( K ), tab_Period( K ) )
           Else
             try_A( K ) = tab_Bottom( K ) &
                           + Mod( try_A( K ) - tab_Bottom( K ), tab_Period( K ) )
           End If
         End If
       End Do

       End Subroutine UPDATE_PARAM


       ! Here is a lightly modified version of the Numerical Recipes 
       ! svdcmp, so you can call it directly instead of the ARK one.
       ! A futher improvement would be to change the statements with
       ! anorm in them to use the epsilon function. 

      SUBROUTINE svdcmp(a, m, n, w, v, flag_svd)

      INTEGER m,n
      real, dimension(:,:) :: a, v
      real, dimension(:) :: w
      ! A flag, set to -51 if 30 iterations are tried.
      integer, intent(out) :: flag_svd

      INTEGER i,its,j,jj,k,l,nm
      REAL anorm,c,f,g,h,s,scale,x,y,z,rv1(m_par)

      flag_svd=0

      g=0.0
      scale=0.0
      anorm=0.0
      do 25 i=1,n
        l=i+1
        rv1(i)=scale*g
        g=0.0
        s=0.0
        scale=0.0
        if(i.le.m)then
          do 11 k=i,m
            scale=scale+abs(a(k,i))
11        continue
          if(abs(scale) > tiny(scale))then
            do 12 k=i,m
              a(k,i)=a(k,i)/scale
              s=s+a(k,i)*a(k,i)
12          continue
            f=a(i,i)
            g=-sign(sqrt(s),f)
            h=f*g-s
            a(i,i)=f-g
            do 15 j=l,n
              s=0.0
              do 13 k=i,m
                s=s+a(k,i)*a(k,j)
13            continue
              f=s/h
              do 14 k=i,m
                a(k,j)=a(k,j)+f*a(k,i)
14            continue
15          continue
            do 16 k=i,m
              a(k,i)=scale*a(k,i)
16          continue
          endif
        endif
        w(i)=scale *g
        g=0.0
        s=0.0
        scale=0.0
        if((i.le.m).and.(i.ne.n))then
          do 17 k=l,n
            scale=scale+abs(a(i,k))
17        continue
          if(abs(scale) > tiny(scale))then
            do 18 k=l,n
              a(i,k)=a(i,k)/scale
              s=s+a(i,k)*a(i,k)
18          continue
            f=a(i,l)
            g=-sign(sqrt(s),f)
            h=f*g-s
            a(i,l)=f-g
            do 19 k=l,n
              rv1(k)=a(i,k)/h
19          continue
            do 23 j=l,m
              s=0.0
              do 21 k=l,n
                s=s+a(j,k)*a(i,k)
21            continue
              do 22 k=l,n
                a(j,k)=a(j,k)+s*rv1(k)
22            continue
23          continue
            do 24 k=l,n
              a(i,k)=scale*a(i,k)
24          continue
          endif
        endif
        anorm=max(anorm,(abs(w(i))+abs(rv1(i))))
25    continue
      do 32 i=n,1,-1
        if(i.lt.n)then
          if(abs(g) > tiny(g))then
            do 26 j=l,n
              v(j,i)=(a(i,j)/a(i,l))/g
26          continue
            do 29 j=l,n
              s=0.0
              do 27 k=l,n
                s=s+a(i,k)*v(k,j)
27            continue
              do 28 k=l,n
                v(k,j)=v(k,j)+s*v(k,i)
28            continue
29          continue
          endif
          do 31 j=l,n
            v(i,j)=0.0
            v(j,i)=0.0
31        continue
        endif
        v(i,i)=1.0
        g=rv1(i)
        l=i
32    continue
      do 39 i=min(m,n),1,-1
        l=i+1
        g=w(i)
        do 33 j=l,n
          a(i,j)=0.0
33      continue
        if(abs(g) > tiny(g))then
          g=1.0/g
          do 36 j=l,n
            s=0.0
            do 34 k=l,m
              s=s+a(k,i)*a(k,j)
34          continue
            f=(s/a(i,i))*g
            do 35 k=i,m
              a(k,j)=a(k,j)+f*a(k,i)
35          continue
36        continue
          do 37 j=i,m
            a(j,i)=a(j,i)*g
37        continue
        else
          do 38 j= i,m
            a(j,i)=0.0
38        continue
        endif
        a(i,i)=a(i,i)+1.0
39    continue
      do 49 k=n,1,-1
        do 48 its=1,30
          do 41 l=k,1,-1
            nm=l-1
            if((abs(rv1(l))+anorm).eq.anorm)  goto 2
            if((abs(w(nm))+anorm).eq.anorm)  goto 1
41        continue
1         c=0.0
          s=1.0
          do 43 i=l,k
            f=s*rv1(i)
            rv1(i)=c*rv1(i)
            if((abs(f)+anorm).eq.anorm) goto 2
            g=w(i)
            h=pythag(f,g)
            w(i)=h
            h=1.0/h
            c= (g*h)
            s=-(f*h)
            do 42 j=1,m
              y=a(j,nm)
              z=a(j,i)
              a(j,nm)=(y*c)+(z*s)
              a(j,i)=-(y*s)+(z*c)
42          continue
43        continue
2         z=w(k)
          if(l.eq.k)then
            if(z.lt.0.0)then
              w(k)=-z
              do 44 j=1,n
                v(j,k)=-v(j,k)
44            continue
            endif
            goto 3
          endif
          if (its == 30) then
            if (debug) print*, '@ ERROR:: No convergence in 30 iterations'
            flag_svd = -51
            return
          end if
          x=w(l)
          nm=k-1
          y=w(nm)
          g=rv1(nm)
          h=rv1(k)
          f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y)
          g=pythag(f,1.0)
          f=((x-z)*(x+z)+h*((y/(f+sign(g,f)))-h))/x
          c=1.0
          s=1.0
          do 47 j=l,nm
            i=j+1
            g=rv1(i)
            y=w(i)
            h=s*g
            g=c*g
            z=pythag(f,h)
            rv1(j)=z
            c=f/z
            s=h/z
            f= (x*c)+(g*s)
            g=-(x*s)+(g*c)
            h=y*s
            y=y*c
            do 45 jj=1,n
              x=v(jj,j)
              z=v(jj,i)
              v(jj,j)= (x*c)+(z*s)
              v(jj,i)=-(x*s)+(z*c)
45          continue
            z=pythag(f,h)
            w(j)=z
            if(abs(z) > tiny(z))then
              z=1.0/z
              c=f*z
              s=h*z
            endif
            f= (c*g)+(s*y)
            x=-(s*g)+(c*y)
            do 46 jj=1,m
              y=a(jj,j)
              z=a(jj,i)
              a(jj,j)= (y*c)+(z*s)
              a(jj,i)=-(y*s)+(z*c)
46          continue
47        continue
          rv1(l)=0.0
          rv1(k)=f
          w(k)=x
48      continue
3       continue
49    continue
      return
      END subroutine svdcmp

      real FUNCTION pythag(a,b)

      REAL a,b
      REAL absa,absb
      absa=abs(a)
      absb=abs(b)
      if(absa.gt.absb)then
        pythag=absa*sqrt(1.+(absb/absa)**2)
      else
        if(absb.eq.0.)then
          pythag=0.
        else
          pythag=absb*sqrt(1.+(absa/absb)**2)
        endif
      endif
      return
      END function pythag

  end module marq


