program Pearson
    implicit none
    integer :: dim, i, j, k, l, m, Pminvalue, Subtnorm1value, Subtnorm2value
    character(80) ::  filename
    real(8) :: Smin, Bmin, Pmin, minCoeff, Coeff, DistM, time
    real(8), allocatable :: wavenumber(:), S(:), B(:), P(:), Sval(:), Bval(:), Pval(:), tempvector1(:),

tempvector2(:)
    real(8), allocatable :: Pnormvector(:), Samplenormvector(:)
    real(8), allocatable :: Svalmat(:,:), Bvalmat(:,:), Pvalmat(:,:), Snorm(:,:), Bnorm(:,:), Pnorm

(:,:), Subtract(:,:)
    real(8), allocatable :: Subtractnorm(:,:,:)

    write(*,*)'Specify the name of the file (format: data.txt):'
10  read(*,*)filename
    open(11,file=filename, status='old',err=20)
    goto 30
20  write(*,*) 'There is no file with this name in the current directory. Please check and try again'
    goto 10
30  continue
    dim=0
    do
	read(11,*,end=40)
	dim=dim+1
    end do
40  continue
    rewind 11
    allocate (wavenumber(dim), S(dim), B(dim), P(dim), Sval(dim), Bval(dim), Pval(dim), tempvector1

(dim), tempvector2(dim))
    allocate (Pnormvector(dim), Samplenormvector(dim))
    allocate (Svalmat(dim,1), Bvalmat(dim,1), Pvalmat(dim,1), Snorm(dim,dim), Bnorm(dim,dim), Pnorm

(dim,dim), Subtract(dim,dim))
    allocate (Subtractnorm(dim,dim,dim))
    do i=1,dim
	read(11,*,end=50) wavenumber(i), S(i), B(i), P(i)
    enddo
50  close(11)

    Smin=minval(S); Sval=S-Smin
    Bmin=minval(B); Bval=B-Bmin
    Pmin=minval(P); Pval=P-Pmin

    do i=1,dim
	do j=1,dim
		Snorm(i,j)=Sval(j)/Sval(i)
		Bnorm(i,j)=Bval(j)/Bval(i)
		Pnorm(i,j)=Pval(j)/Pval(i)
	enddo
    enddo

    Subtract=Snorm-Bnorm

    do i=1,dim
	do j=1,dim
		do k=1,dim
			Subtractnorm(i,j,k)=Subtract(j,k)/Subtract(j,i)
		enddo
	enddo
    enddo

    minCoeff=1.0; Pminvalue=0; Subtnorm1value=0; Subtnorm2value=0

    do i=1,dim
	do j=1,dim
		tempvector1(j)=Pnorm(i,j)
	enddo
	do k=1,dim
		do l=1,dim
			do m=1,dim
				tempvector2(m)=Subtractnorm(k,l,m)
			enddo

			Coeff=DistM(tempvector1,tempvector2, dim)
			if (Coeff.lt.minCoeff) then
				minCoeff=Coeff
				Pminvalue=i
				Subtnorm1value=k
				Subtnorm2value=l
			endif
		enddo
	enddo
    enddo

    print *, 'Minimal Correllation Coefficient: ', minCoeff
    print *, 'Pminvalue: ', Pminvalue
    print *, 'Subtnorm1value: ', Subtnorm1value
    print *, 'Subtnorm2value: ', Subtnorm2value
    call CPU_TIME(time)
    print *, 'Elapsed time: ', time, 'seconds'

    open(33,file='ResultsManh.auto')
    do i=1,dim
	Pnormvector(i)=Pnorm(Pminvalue,i)
        Samplenormvector(i)=Subtractnorm(Subtnorm1value,Subtnorm2value,i)
       	write(33,'(3F10.5)') wavenumber(i), Pnormvector(i),Samplenormvector(i)
    enddo
end

    real(8) function DistM(v, u, dim)
    integer dim
    real(8) v(dim), u(dim), Distance
    	Dif=0.0
	do i=1, dim
		Dif=Dif+((v(i)-u(i))*(v(i)-u(i)))
	enddo
    	Distance=Dif/dim
    	DistM=Distance
    end