CHARACTER*40 INFILE1,INFILE2,FILENAME
	CHARACTER*70 ASTRING
	LOGICAL EFLG
	real*8 mat(10000,10),Y(10000),YTM(10),P(10)
	real*8 MTM(4,4)
	integer*4 I,j,N,M
	CHARACTER*4   A
c 	M is the number of observations, number of rows in mat()
c 	N is the number of variables, number of columns in mat()
c	form i,j th element of mtm by dotting ith into j'th column of mat.
c	form ytm by dotting y unto each column of mat
c	multiply ytm by inverse of mtm, gives best coefficients.

c. Put application here to load mat() and Y(), or 
c   use this as a callable subroutine

C read atoms coordinates in PDB format and find coeficients of best-plane
C	passing through them. See notes on windows machine
C grep HEM an2ref.pdb | awk '$5~/C/ && $6~/501/' > ! blo.pdb
C awk '$3!~/^CB|^CG|^O/' blo.pdb >! temp.pdb

C final touch- subtract residual from each point to get projection on plane
C Write that out as new pdbfile. (not done yet)

6	WRITE(6,*) 'ENTER NAME OF pdb file:'
	READ(5,'(A)') INFILE1
c1200	FORMAT(A)
	IF (INFILE1(:5).EQ.'     ') STOP
	  INQUIRE(FILE=INFILE1, NAME=FILENAME,EXIST=EFLG)
	  IF (.NOT.EFLG) GOTO 6
        OPEN (UNIT=2,FILE=INFILE1,STATUS='OLD')
	
c 	read coord 
	I=1        
50	read(2,'(A)',end=70) ASTRING
	if ((ASTRING(1:6).ne.'ATOM  ').and.(ASTRING(1:6).ne.'HETATM')) 
     . goto 50
     	if (I.gt.10000) stop 'dimensioned for 10000 atoms max!'
	read (ASTRING,51) A, A,J,(mat(i,j),j=1,3)
51 	format (12x,A4,A6,I4,4x,3d8.3,2f6.2)
	Y(I)=1.0
60	I=I+1
	goto 50
70	close (unit=2)
	M=I-1
	N=3
	Write(6,*) 'finding best plane through',M,' atoms'



c	form i,j th element of mtm by dotting ith into j'th column of mat.
	do 100 i=1,n
	do 100 j=1,i
	x=0
	do 90 k=1,m
90	x=x+mat(k,i)*mat(k,j)
	mtm(i,j)=x
100	mtm(j,i)=x
c	form ytm by dotting y into each column of mat
	do 200 i=1,n
	x=0
	do 190 k=1,m
190	x=x+mat(k,i)*y(k)
200	ytm(i)=x
	
c	multiply ytm by inverse of mtm, gives best coefficients.
	call matinv(mtm,n,x)
c	do 15 i=1,3
c15	write(6,'(3F16.10)') (mtm(i,j),j=1,3)	
	do 300 i=1,n
	x=0
	do 290 j=1,n
290	x=x + mtm(i,j)*ytm(j)
300	p(i)=x


350	Write(6,*) 'best-fit coefficients:'
	Write(6,'(3F16.10)') (p(i),i=1,n)
	write(6,*) 'equation of plane:', P(1),'*X + ', P(2),' * Y + ', P(3), ' * Z = 1.0'
	x=0.0
	do 360 i=1,n
360	x=x+p(i)*p(i)
	x=sqrt(x)
	dist=1/x
	Write(6,*) 'distance of plane from origin:',dist
	do 370 i=1,n
370	p(i)=p(i)/x
	Write(6,*) 'normalized vector orthogonal to bf plane:'
	Write(6,'(3F16.10)') (p(i),i=1,n)
	Write(6,*) 'coordinates of closest point in plane to origin'
	Write(6,'(3F16.10)') (p(i)/x,i=1,n)
	pi=4*atan(1.)
	write(6,*) 'angle of orth vector to Z axis:', acos(p(3))*180/pi,'deg' 
	 
c  Calculate residuals
	sum=0.0
	do 400 k=1,m
	x=0
	do 390 i=1,3
390	x=x+mat(k,i)*p(i)
	x=x-dist
CCC	write(6,'(I6,7f10.2)') k,(mat(k,i),i=1,3), x, sum
400	sum=sum+x*x
	write(6,*) m,' Atoms have RMS deviation from plane:',sqrt(sum/m)


	
	end



	SUBROUTINE MATINV (ARRAY, NORDER, DET)
	DOUBLE PRECISION ARRAY, AMAX, SAVE
2	DIMENSION ARRAY(4,4), IK (10), JK (10)
10	DET = 1
11 	DO 100  K=1, NORDER
C
C	FIND LARGEST ELEMENT ARRAY(I,J) IN REST OF MATRIX
C
	AMAX=0
21 	DO 30 I=K, NORDER
     	DO 30 J=K, NORDER
23 	IF (DABS (AMAX) - DABS (ARRAY(I,J) ) ) 24, 24, 30
24 	AMAX=ARRAY(I,J)
     	IK(K) =I
     	JK(K) =J
30 	CONTINUE
C
C INTERCHANGE ROWS AND COLUMNS TO PUT AMAX IN ARRAY(K,K)
C
31 	IF (AMAX) 41, 32, 41
32 	DET=0.
     	GO TO 140
41 	I=IK (K)
     	IF (I-K) 21, 51, 43
43 	DO 50 J=1, NORDER
     	SAVE= ARRAY(K,J)
      	ARRAY(K,J)= ARRAY(I,J)
50 	ARRAY(I,J)= -SAVE
51 	J=JK (K)
     	IF (J-K) 21, 61, 53
53 	DO 60 I=1, NORDER
     	SAVE= ARRAY(I,K)
     	ARRAY(I,K)= ARRAY(I,J)
60 	ARRAY(I,J)= -SAVE
C
C	ACCUMULATE ELEMENTS OF INVERSE MATRIX
C
61 	DO 70 I=1, NORDER
     	IF(I-K) 63, 70, 63
63 	ARRAY(I,K) = -ARRAY(I,K)  / AMAX
70 	CONTINUE
71 	DO 80 I=1, NORDER
     	DO 80 J=1, NORDER
     	IF (I-K) 74, 80, 74
74 	IF (J-K) 75, 80, 75
75 	ARRAY(I,J) = ARRAY(I,J) + ARRAY(I,K) *ARRAY(K,J)
80 	CONTINUE
81 	DO 90 J=1, NORDER
     	IF (J-K) 83, 90, 83
83 	ARRAY(K,J) = ARRAY(K,J) / AMAX
90 	CONTINUE
     	ARRAY(K,K) = 1. / AMAX
100 	DET=DET *AMAX
c	write(6,*)'DET=',det
C
C 	RESTORE ORDERING OF MATRIX
C 
101 	DO 130 L=1, NORDER
       	K= NORDER  - L +1
       	J= IK (K)
       	IF (J-K) 111, 111, 105
105 	DO 110  I=1, NORDER
       	SAVE= ARRAY(I,K)
        ARRAY(I,K)= -ARRAY(I,J)
110 	ARRAY(I,J)= SAVE
111 	I= JK(K)
       	IF (I-K) 130, 130, 113
113 	DO 120  J=1, NORDER
       	SAVE= ARRAY(K,J)
       	ARRAY(K,J) = -ARRAY(I,J)
120 	ARRAY(I,J)= SAVE
130 	CONTINUE
140  	RETURN
	END