Skip to content
Snippets Groups Projects
Commit 2d464b3b authored by mnb's avatar mnb Committed by afarrell
Browse files

Added RBSTFIT0 to l1_utils.f90 a variant of the robust line fitting routine

	Added RBSTFIT0 to l1_utils.f90 a variant of the robust line fitting routine
	RBSTFIT that assumes that line is defined by a single scale parameter. i.e.
	Y = SLOPE*X as opposed to Y = SHIFT + SLOPE*X
parent 2346f6e7
Branches
Tags 1.96
No related merge requests found
......@@ -117,6 +117,76 @@ END SUBROUTINE RSTFIT
! -----------------------------------------------------------------------
! -----------------------------------------------------------------------
SUBROUTINE RSTFIT0(X,Y,NDATA,SLOPE)
!+
! Purpose:
! Robust straight line fit using simple L1 algorithm method.
!
! Description:
! This routine fits the line Y = B*X to a set of points X,Y by the
! criterion of least absolute deviations.
!
! This version calls on the L1 subroutine as taken from the TOM archives
! (referred to as TOMS algorithm 478). This is a fast stable and simple
! routine which is satisfies our needs. See l1.F95 for further details.
!
! Notice the solution is not guaranteed to be unique and in such cases the
! answer is determined by the order of the data points given. However this
! case is sufficiently rare to not worry about.
!
! Arguments:
! X = REAL VECTOR(NDATA) (Given)
! X data points to be fitted
! Y = REAL VECTOR(NDATA) (Given)
! Y data points to be fitted
! NDATA = INTEGER (Given)
! Size of X and Y vectors
! SLOPE = REAL (Returned)
! Slope (or gradient)
!
!-
IMPLICIT NONE ! No implicit typing
! Arguments
INTEGER NDATA
REAL X(NDATA),Y(NDATA)
REAL SLOPE
! No of parameters to solve for
INTEGER,PARAMETER::N=1
! N2 = N+2
INTEGER,PARAMETER::N2=N+2
! Local variables :
INTEGER M,M2
REAL TOLER
REAL A(NDATA+2,N2) ! A has extra rows and cols to store return values
REAL XS(N) ! solution vector
REAL B(NDATA) ! Y data input vector
REAL E(NDATA) ! workspace vector
INTEGER S(NDATA)
! ???
M = NDATA
M2 = M+2
TOLER = 10**(-12*2/3) ! safe value below which is deemed zero
! fit function Yi=a*Xi+b*1.0, so data columns are Xi and 1.0 for i=1,ndata
! These are placed into A array whilst Y values are plced in B array
A(1:M,1) = X(1:M)
B(1:M) = Y(1:M)
! Call L1 routine, solution returned in XS vector
CALL L1(M,N,M2,N2,A,B,TOLER,XS,E,S)
SLOPE = XS(1)
END SUBROUTINE RSTFIT0
! -----------------------------------------------------------------------
! -----------------------------------------------------------------------
SUBROUTINE RSTFIT2(X,Y,NDATA,X2,Y2,NDATA2)
!+
! Purpose:
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment