!Create X summation matrix for Dth polynomial
DO i = 1,d+1
  DO j = 1,d+1
     !Go through length of point matrix for summation
	sum1 = 0.0
         DO k = 1,n
	  IF (j+i-2 == 0) THEN
	      sum1 = sum1 + 1
	   ELSEIF (p(k,1) == 0.0) THEN
	       sum1 = sum1 + 0
	    ELSE
	       sum1 = sum1 + p(k,1)**(j+i-2)
	    ENDIF
        ENDDO
        A(i,j) = sum1
    ENDDO
ENDDO
	
!Create Y summation matrix of Fn and Xn variables
	DO i = 1,d+1
		sum1 = 0.0
		DO k = 1,n
			IF (i-1 == 0) THEN
			    sum1 = sum1 + p(k,2)
			ELSEIF (p(k,1) == 0.0) THEN
			    sum1 = sum1 + 0
			ELSE
			     sum1 = sum1 + p(k,2)*(p(k,1)**(i-1))
			ENDIF
		ENDDO
		B(i) = sum1
	ENDDO
	
	!Inverse A
	CALL inverse(d+1,A,AI)
	
	!Multiply A^-1 by B to get the coefficients
	coef = MATMUL(AI, B)
	
	!Now check standard deviation versus standard error
	!If SE < SD, then it’s a good fit
	
	!Standard Deviation
	favg = SUM(p(1:n,2))/real(n) !F average / F bar
	sum1 = 0.0
	DO i = 1,n
		sum1 = sum1 + (p(i,2) - favg)**2
	ENDDO
	sd = SQRT(sum1/(real(n)-1.))
	
	!Standard Error
	sum1 = 0.0
	sum2 = 0.0
	DO i = 1,n
		sum2 = 0.0
		DO j = 1,d+1
			IF (j-1 == 0) THEN
			  sum2 = sum2 + coef(j)
			ELSEIF (p(i,1) == 0) THEN
			   sum2 = sum2
			ELSE
			   sum2 = sum2 + coef(j)*(p(i,1)**(j-1))
			ENDIF
		ENDDO
		sum1 = sum1 + (p(i,2) - sum2)**2
	ENDDO
	se = SQRT(sum1/(n-2))
	
	!Other computations
	!Norm Error
	sum1 = 0.0
	DO i = 1,n
		sum2 = 0.0
		DO j = 1,d+1
			IF (j-1 == 0) THEN
			   sum2 = sum2 + coef(j)
			ELSEIF (p(i,1) == 0) THEN
			    sum2 = sum2
			ELSE
			   sum2 = sum2 + coef(j)*(p(i,1)**(j-1))
			ENDIF
		ENDDO
		sum1 = sum1 + (p(i,2) - sum2)**2
	ENDDO
	ne = SQRT(sum1)