module fit_interp_etc
contains

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!*******************************************************
!*   PROGRAM TO DEMONSTRATE ONE DIMENSIONAL OPERATION  *
!*     OF THE MULTI-NONLINEAR REGRESSION SUBROUTINE    *
!* --------------------------------------------------- *
!* Ref.: BASIC Scientific Subroutines Vol. II,         *
!*       By F.R. Ruckdeschel, Byte/McGRAW-HILL,1981    *
!*       [BIBLI 01].                                   *
!*                                                     *
!*                  F90 version by J-P Moreau, Paris   *
!* --------------------------------------------------- *
!* SAMPLE RUN:                                         *
!*                                                     *
!* MULTI-NON LINEAR REGRESSION                         *
!*                                                     *
!* Number of data points (mini 3, maxi 25)........: 11 *
!* Degree of the polynomial to be fitted (maxi 10): 4  *
!*                                                     *
!* Input the data in (x,y) pairs as prompted:          *
!*  1   X  Y = 1,1                                     *
!*  2   X  Y = 2,2                                     *
!*  3   X  Y = 3,3                                     *
!*  4   X  Y = 4,4                                     *
!*  5   X  Y = 5,5                                     *
!*  6   X  Y = 6,6                                     *
!*  7   X  Y = 7,7                                     *
!*  8   X  Y = 8,8                                     *
!*  9   X  Y = 9,9                                     *
!* 10   X  Y = 0,0                                     *
!* 11   X  Y = 1,1                                     *
!*                                                     *
!* The calculated coefficients are:                    *
!*  0    0.0000                                        *
!*  1    1.0000                                        *
!*  2    0.0000                                        *
!*  3    0.0000                                        *
!*  4    0.0000                                        *
!*                                                     *
!* Standard deviation:  0.000000                       *
!*                                                     *
!*******************************************************
subroutine LEASTSQR(x,y,m,n,D)
!x is x data, y is y data, m is number of data points, n is the order of the polynomial

implicit none
INTEGER, PARAMETER :: iknd=SELECTED_INT_KIND(8)  ! nobody does this
INTEGER, PARAMETER :: DP=SELECTED_REAL_KIND(14)  ! this is double precision
INTEGER, PARAMETER :: SP=SELECTED_REAL_KIND(6)  ! this is single precision   
    
integer :: mmax,nmax
parameter(MMAX=365,NMAX=2)

real(DP) ::  X(MMAX), Y(MMAX), X2(MMAX), Y2(MMAX) 
integer :: i,m,n,k,n2
real(DP) ::  dd
real(DP) :: D(3)
real(DP) ::  Z(MMAX,NMAX+1)          ! maxi m=25 data points
real(DP) ::  A(MMAX,MMAX)            ! maxi n=10 order of regression
real(DP) ::  B(MMAX,2*MMAX)
real(DP) ::  C(MMAX,MMAX),check(MMax)
real(DP),allocatable :: xuse(:),yuse(:)

!allocate(D(n))
allocate(xuse(m))
allocate(yuse(m))

xuse = x(1:m)
yuse = y(1:m)

  if(maxval(y) > 100)then
      k = 1
      do i = 1,365
        if(y(i) > 0)then
            y2(k) = y(i)
            x2(k) = x(i)
            k = k+1
        endif
      enddo
  else
    y2 = y
    x2 = x
    k = m+1
  endif

   ! Call coefficients generation subroutine
  call Gene_Coeffs(k-1,n,X2(1:k-1),Z)
  ! Call order n regression subroutine and print result
  n2 = n + 1
  call LstSqrN(m,n2,A,B,C,D,Y2(1:k-1),Z)

end subroutine leastsqr


! Matrix transpose B = Tr(A)
! A and B are n by m matrices
Subroutine Transpose(n,m,A,B)
  parameter(MMAX=365,NMAX=2)
  integer n,m
  real*8  A(MMAX,MMAX), B(MMAX,2*MMAX)
  integer i,j
  do i = 1, n
    do j = 1, m
      B(i,j) = A(j,i)
    end do
  end do
  return
end subroutine Transpose

! Matrix save (A in C) 
Subroutine A_IN_C(n1,n2,A,C) 
  parameter(MMAX=365,NMAX=2)
  integer n1,n2
  real*8  A(MMAX,MMAX), C(MMAX,MMAX)
  integer i1,i2
  if (n1 * n2.ne.0) then  
    do i1 = 1, n1
      do i2 = 1, n2
        C(i1,i2) = A(i1,i2)
      end do
    end do
  end if
  return
end subroutine A_IN_C

! Matrix save (B in A) 
Subroutine B_IN_A(n1,n2,A,B) 
  parameter(MMAX=365,NMAX=2)
  integer n1,n2
  real*8  A(MMAX,MMAX), B(MMAX,2*MMAX)
  integer i1,i2
  if (n1 * n2.ne.0) then  
    do i1 = 1, n1
      do i2 = 1, n2
        A(i1,i2) = B(i1,i2)
      end do
    end do
  end if
  return
end subroutine B_IN_A

! Matrix save (C in B)
Subroutine C_IN_B(n1,n2,B,C) 
  parameter(MMAX=365,NMAX=2)
  integer n1,n2
  real*8  C(MMAX,MMAX), B(MMAX,2*MMAX)
  integer i1,i2
  if (n1 * n2.ne.0) then  
    do i1 = 1, n1
      do i2 = 1, n2
        B(i1,i2) = C(i1,i2)
      end do
    end do
  end if
  return
end Subroutine C_IN_B

! Matrix save (C in A)
Subroutine C_IN_A(n1,n2,A,C) 
  parameter(MMAX=365,NMAX=2)
  integer n1,n2
  real*8  A(MMAX,MMAX), C(MMAX,MMAX)
  integer i1,i2
  if (n1 * n2.ne.0) then  
    do i1 = 1, n1
      do i2 = 1, n2
        A(i1,i2) = C(i1,i2)
      end do
    end do
  end if
  return
end Subroutine C_IN_A

! Matrix multiplication C = A * B
! A is m1 by n1, B is m2 by n2 ( m2 = n1 )
! Result C is m1 by n2.
Subroutine Matmult(m1,n1,n2,A,B,C) 
  parameter(MMAX=365,NMAX=2)
  integer m1,n1,n2
  real*8  A(MMAX,MMAX), B(MMAX,2*MMAX), C(MMAX,MMAX)  
  integer i,j,k
  do i = 1, m1
    do j = 1, n2	 
      C(i,j) = 0
      do k = 1, n1
        C(i,j) = C(i,j) + A(i,k) * B(k,j)
      end do
    end do
  end do
  return
end Subroutine Matmult

! Matrix inversion: B = Inv(A) by Gauss-Jordan method
! A and B are n by n matrices
Subroutine Matinv(n,A,B)  
  ! Labels: 10, 20, 30
  parameter(MMAX=365,NMAX=2)
  integer n
  real*8  A(MMAX,MMAX), B(MMAX,2*MMAX)  
  integer i,j,k
  real*8 bb

  do i = 1, n
    do j = 1, n
      B(i,j + n) = 0.d0
      B(i,j) = A(i,j)
    end do
    B(i,i + n) = 1.d0
  end do  

  do k = 1, n
    if (k.eq.n) goto 10
    m = k
    do i = k+1, n
      if (abs(B(i,k)) > abs(B(m,k)))  m = i
    end do
    if (m == k) goto 10
    do j = k, 2*n
      bb = B(k,j)
      B(k,j) = B(m,j)
      B(m,j) = bb
    end do
10  do j = k+1, 2*n 
      B(k,j) = B(k,j) / B(k,k)
    end do
    if (k.eq.1) goto 20
    do i = 1, k-1
      do j = k+1, 2*n
        B(i,j) = B(i,j) - B(i,k) * B(k,j)
      end do
    end do
    if (k.eq.n) goto 30
20  do i = k+1, n
      do j = k+1, 2*n
        B(i,j) = B(i,j) - B(i,k) * B(k,j)
      end do
    end do
  end do    ! k loop

30  do i = 1, n
      do j = 1, n   
        B(i,j) = B(i,j + n)
      end do
    end do
  return
end Subroutine Matinv ! Matinv()

! Coefficient matrix generation
Subroutine Gene_Coeffs(m,n,X,Z)
  parameter(MMAX=365,NMAX=2)
  integer m, n
  real*8  X(MMAX), Z(MMAX,NMAX+1)   
  integer i,j
  real*8 b
  do i = 1, m
    b = 1.d0
    do j = 1, n+1
      Z(i,j) = b
      b = b * X(i)
    end do
  end do
  return
end Subroutine Gene_Coeffs

! Least squares multidimensional fitting subroutine
Subroutine LstSqrN(m,n,A,B,C,D,Y,Z)
  parameter(MMAX=365,NMAX=2)
  integer m,n,m2
  real*8 A(MMAX,MMAX),B(MMAX,2*MMAX),C(MMAX,MMAX)
  real*8 D(NMAX+1), Y(MMAX), Z(MMAX,NMAX+1) 
  integer i,j,m4,n4
  m4 = m
  n4 = n
  do i = 1, m
    do j = 1, n 
      A(i,j) = Z(i,j)
    end do
  end do
  call Transpose(n,m,A,B)                ! B=Transpose(A) 
  call A_IN_C(m,n,A,C)                   ! move A to C
  call B_IN_A(n,m,A,B)                   ! move B to A
  call C_IN_B(m,n,B,C)                   ! move C to B
  call Matmult(n,m,n,A,B,C)              ! multiply A and B, result in C
  call C_IN_A(n,n,A,C)                   ! move C to A
  call Matinv(n,A,B)                     ! B=Inverse(A) 
  m2 = m4                                 ! restore m
  call B_IN_A(n,n,A,B)                   ! move B to A
  do i = 1, m2
    do j = 1, n
      B(j,i) = Z(i,j)
    end do
  end do
  call Matmult(n,n,m2,A,B,C)              ! multiply A and B, result in C 
  call C_IN_A(n,m2,A,C)                   ! move C to A
  do i = 1, m2
    B(i,1) = Y(i)
  end do
  call Matmult(n,m2,1,A,B,C)              ! multiply A and B

  ! Product C is n by 1 - Regression coefficients are in C(i,1)
  ! and put for convenience into vector D(i).
  do i = 1, n
    D(i) = C(i,1)
  end do
  return
end Subroutine LstSqrN ! LstSqrN()

! Standard deviation calculation for a polynomial fit
Subroutine Standard_deviation(m,n,dd,D,X,Y) 
  parameter(MMAX=365,NMAX=2)
  integer m,n
  real*8  dd,D(NMAX+1),X(MMAX),Y(MMAX) 
  integer i, j
  real*8 b, yy
  dd = 0.d0
  do i = 1, m
    yy = 0.d0 
    b = 1.d0
    do j = 1, n+1
      yy = yy + D(j) * b
      b = b * X(i)
    end do
    dd = dd + (yy - Y(i)) * (yy - Y(i))
  end do
  if (m - n - 1 > 0) then
    dd = dd / (m - n - 1)
  else
    dd = 0.d0
  end if
  dd = dsqrt(dd)
  return
end Subroutine Standard_deviation

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!linear curve fit
subroutine fit(Tmean,Twbmean,m,c)
!linear curve fit, m is the number of data points, a and b are the coefficients with
!y = ax + b
use global, only : DP
!adapted from Numerical Mathematics and Computing by Cheney and Kincaid, pg 497-498

implicit none
real(DP) :: Tmean(365),Twbmean(365)
real(DP) :: a,b,c(2)
real(DP) :: p,q,r,s,d
integer :: m

p = sum(Tmean)
q = sum(Twbmean)
r = sum(Tmean*Twbmean)
s = sum(Tmean*Tmean)
d = (m+1)*s - p*p

a = 1/d*( (m+1)*r - p*q )
b = 1/d*(s*q - p*r)

c(2) = a
c(1) = b
end subroutine fit

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!*****************************************************
!* Sorts an array ARR of length N in ascending order *
!* by straight insertion.                            *
!* ------------------------------------------------- *
!* INPUTS:                                           *
!*	    N	  size of table ARR                  *
!*          ARR	  table to be sorted                 *
!* OUTPUT:                                           *
!*	    ARR   table sorted in ascending order    *
!*                                                   *
!* NOTE: Straight insertion is a Nē routine and      *
!*       should only be used for relatively small    *
!*       arrays (N<100).                             *
!*****************************************************         
SUBROUTINE PIKSRT(N,ARR,out)
use global, only: DP
implicit none
real(DP) :: ARR(N),a,out(N,2)
integer :: i,j,n,rvect(n),b
    do i = 1,n
        rvect(i) = i
    enddo
  do j=2,N
    a = ARR(j)
    b = rvect(j)
    do i=j-1,1,-1
      if (ARR(i)<=a) goto 10
      ARR(i+1)=ARR(i)
      rvect(i+1) = rvect(i)
    end do
	i=0
10  ARR(i+1) = a
    rvect(i+1) = b
  end do
  out(:,1) = arr
  out(:,2) = rvect
  return
END subroutine PIKSRT

end module fit_interp_etc