CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C									C
C	PROGRAM   DEFINE_Structure					C
C									C
C	A fortran program to determine the secondary and first level	C
C	supersecondary structure of a protein from alpha carbon atom	C
C	coordinates.							C
C									C
C	Copyright (c) 1988, 1994 Yale University			C
C									C
C	Authors: F.M. Richards and C.E.Kundrot				C
C									C
C	Reference for Version 1: F.M.Richards and C.E.Kundrot Proteins	C
C	     3: 71,84 (1988)  Identification of Structural Motifs from	C
C	Protein Coordinate Data: Secondary Structure and First-Level	C
C	Supersecondary Structure.					C
C									C
C	The program DEFINE_Structure was developed in the Center for	C
C	Structural Biology at Yale University. Its use is subject to 	C
C	the following conditions.					C
C									C
C     	a) Use of the program is restricted to the laboratory or	C
C	organization to which it has been supplied. The program may not	C
C	be sold nor may copies be distributed to third parties without	C
C	the express permission of Yale University.			C
C									C
C	b) This program is provided on an "as is" basis. Yale in no	C
C	way warrants this program or its results.			C
C									C
C     	c) Yale is under no obligation to provide any services by	C
C	way of maintenance, update, or otherwise. However, if a cost	C
C	recovery fee has been collected, future updates, if any, will   C
C	be supplied free of charge.					C
C									C
C     	d) Reports or publications based on this program should		C
C	contain an acknowledgement in a form commonly used in academic	C
C	research.							C
C									C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       Modifications
C       Version 2 April 1994 Dennis Xu, Art Perlo
C       (1) Convert graphics from UGS to Post-Script
C       (2) Common code for VAX/VMS and SGI/IRIX
C       (3) Fixed minor bugs
C
C	Further modifications by P.J. Fleming, April 1995
C	to fix Post_Script output to reflect different line widths
C	and other minor bugs
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C====================    PROGRAM DEFINE_STRUCTURE   ==================
C=======================    SHORT TITLE =  MASK     ==================
C=======================  (updated on 8 Mar 88)     ==================
C=====================================================================
C	This is Version 1.1 of this program. It differs from the 
C	first version described in Richards and Kundrot, J. Mol. Biol.
C	(1988) only by the addition of the secondary unit FIVETRN 
C       which roughly fills the gap between SHARPTRN and LOOPS.
C=====================================================================

C******		Program to produce secondary and supersecondary
C******		structure lists and diagonal maps of character matrix
C******		showing the structural motifs. The tables and maps are
C******		based on difference distances derived from CA distances
C******		in the coordinate files and a set of masks for distances
C******		is various standard structures.

C*********************************************************************

C******		INPUT :  XXXX.PDB -  a Brookhaven Data Bank File
C******			The PDB file must have residual numbered consequently
C******		        CHECK.FIL -  an auxiliary file named which 
C******                              contains the allowed errors
C******				     against which the distance
C******				     differences are evaluated.

C*********************  SAMPLE CHECK.FIL  ********

C**** 1.0,		CHKHA  FOR AHELIX-alpha-max rms error
C**** 1.0,		CHKH3  FOR AHELIX-310  -max rms error
C**** 1.0,		CHKSTR FOR SHARPTRN    -max rms error
C**** 1.0,		CHKFIV FOR FIVETRN     -max rms error
C**** 1.0,		CHKTRN FOR BSTRAND     -max rms error
C****10.0,              CHKLP1 FOR LOOPS       -max end-end dist
C**** 0.667,            CHKLP2 FOR    LOOPS    -end-end/max dist
C**** 1.5,		CHKSHT FOR BSHEET      -max rms error
C**** 6.0,		CHKSDT FOR    BSHEET   -max dist close pairs
C****10.0,		CHKHSD FOR SUPERSS     -max DMIN-helix/strand
C****10.0,		CHKSHD FOR SUPERSS     -max DMIN-strand/helix
C****15.0,		CHKHHD FOR SUPERSS     -max DMIN-helix-helix
C****		SEPMAX FOR SUPERSS     -max allowed atom sepn for inclusion
C****		SEPMIN FOR SUPERSS     -min     "     "    "   "       "
C****		AXDMAX FOR SUPERSS     -max allowed axis sepn for inclusion
C****		AXDMIN FOR SUPERSS     -min     "     "    "   "       "
C****		PT1MAX FOR SUPERSS     -range of positions 	
C****		PT1MIN FOR SUPERSS     - for Point 1 for inclusion
C****		PT2MAX FOR SUPERSS     -range of positions
C****		PT2MIN FOR SUPERSS     - for Point 2 for inclusion

C*********************************************************************

C		          SAMPLE COMMAND FILE
C			  -------------------
C	$ !
C	$ ! Command File DEFINE_S.COM
C	$ ! 
C	$ ! @DEFINE_S XXXX      where XXXX is the PDB file
C	$ ! 			idenfitcation symbol
C	$ ! For batch submission
C	$ ! SUBMIT/QUEUE=SYS$BATCH/NOTIFY/PARAMETER="XXXX" DEFINE_S 
C	$ !
C	$ SET DEFAULT [RICHARDS.MASK]	
C	$ ! 
C	$ ! Set up file names for input and output
C	$ !     
C	$ PDBFIL :='P1'
C	$ INPDB =PDBFIL+".PDB"
C	$ OUTDAT=PDBFIL+".DAT"
C	$ OUTHDR=PDBFIL+".HDR"
C	$ OUTSSS=PDBFIL+".SSS"
C	$ !
C	$ ! open and write file for data input to program MASK
C	$ !
C	$ OPEN/WRITE FTMP [RICHARDS.MASK]FILE.DAT
C	$ WRITE FTMP INPDB		! PDB file for primary input
C	$ WRITE FTMP OUTDAT		! Data output file for printing
C	$ WRITE FTMP OUTHDR		! PDB header notes for printing
C	$ WRITE FTMP OUTSSS		! Data output file for storage
C	$ WRITE FTMP "N"	! Print distance matrix ?
C	$ WRITE FTMP "N"	! Print character matrix ?
C	$ WRITE FTMP "Y"	! Print SS structure ?
C	$ WRITE FTMP "Y"	! Print helix-helix boxes ?
C	$ WRITE FTMP "Y"	! Print hlx-str & str-hlx boxes ?
C	$ WRITE FTMP "Y"	! Print strand-strand boxes ?
C	$ WRITE FTMP "IMA3320"	! Printer indentification for UGSYS package
C	$ CLOSE FTMP
C	$ DEFINE/USER_MODE SYS$INPUT FILE.DAT
C	$ RUN MASK
C	$ DEFINE/USER_MODE SYS$INPUT SYS$COMMAND
C	$ IMPRINT/LANDSCAPE/TWO_UP/DELETE 'OUTDAT'
C	$ IMPRINT/LANDSCAPE/TWO_UP/DELETE 'OUTHDR'
C	$ IMPRINT/DELETE SYS$SCRATCH:UGIM.IMP

C*********************************************************************

C******		Functions or subroutines called: 

C******		Library	  Library    Library   Library   Library
C******		SSUTIL	  SECSTR     AXISCHK   SUPERSS   MXMAP
C******		-------   -------    -------   -------   -------
C******		CPUTIM	  CISPEP     AXISCHK   BSHEET    MXMAP
C******		DATE(SYS) AHELIX     HAXIS     LISTBSH	 UGXERR
C******		TIME(SYS) SHARPTRN   SAXIS     SUPERSS
C******		ASKFIL	  FIVETRN    VDOTV     RELATE
C******		LOADCRD	  BSTRAND    VXV       MAPBOX
C******		PRDIST    LOOPS      POINT     CODE
C******		PRDIS2	  LIST2ND	       DET3R

C*********************************************************************

	IMPLICIT   INTEGER (I-N)
	PARAMETER  (MAXLEN=360 ) 
	PARAMETER  (LINS  =1600)
	PARAMETER  (ICP   =25  )
C*********************************************************************

	COMMON /FILES/FNAME

C********************************************************************

	DIMENSION XA(MAXLEN),YA(MAXLEN),ZA(MAXLEN),LI(MAXLEN)
	DIMENSION DIST(MAXLEN,MAXLEN),DMIN(MAXLEN),NMIN(MAXLEN)
	DIMENSION ALPHA(50), BETA(20), TURN(3,5),FIVE(4,6),BB(13)
	DIMENSION IFIRST(100), ILAST(100), CURV(100)
	DIMENSION FE(175),C3B(175),C3C(175)
	DIMENSION HA3(50,50),A310(3,50)
	DIMENSION HLXWT(30),STRWT(30)
	DIMENSION HLXCOS(ICP,3),HLXPT(ICP,3),STRCOS(ICP,3),STRPT(ICP,3)
	DIMENSION IHSTRT(ICP),IHEND(ICP),ISSTRT(ICP),ISEND(ICP)
	DIMENSION HRAD(ICP),HDIST(ICP),SRAD(ICP),SDIST(ICP)
	DIMENSION HTHA(ICP),HTHMAX(ICP),STHA(ICP),STHMAX(ICP)

	REAL*4    SIGMA(30)

	CHARACTER*1  LINE(MAXLEN),ASTER
	CHARACTER*6  KEYWRD
	CHARACTER*30 FNAME,GNAME,DEVICE
	CHARACTER*9  PRDATE
	CHARACTER*8  PRTIME
	CHARACTER*7  COLR(3,LINS)
	CHARACTER*4  ATM(MAXLEN)

	INTEGER*4  IX(3,LINS),IY(3,LINS),BIT(3,LINS),NCNT(3)

	BYTE        INT(12),DIST2(MAXLEN,MAXLEN),PT1,PT2
	BYTE	    HELIX,H310,TRN,TRN5,STRAND,SHORT,SHEET,LOOP
	BYTE	    TT(3), SYMBOL(100)
	BYTE	    ABC(26),SSCALC,LPDIS2
	BYTE	    HHPRNT,HSPRNT,SSPRNT

	INTEGER     SEQNUM, OFFSET, PAIR(100,6)

C******************************************************************

	DATA INT/'0','1','2','3','4','5','6','7','8','9','-',' '/

	DATA ABC/'A','B','C','D','E','F','G','H','I','J','K','L','M',
	1	 'N','O','P','Q','R','S','T','U','V','W','X','Y','Z'/

C******************************************************************

	DATA ALPHA/
	1   0.00, 3.75, 5.36, 5.02, 6.11, 8.53, 9.75,10.43,12.18,14.09,
	2  15.15,16.32,18.19,19.77,20.85,22.33,24.13,25.49,26.70,28.36,
	3  30.01,31.27,32.65,34.35,35.84,37.11,38.64,40.30,41.67,43.06,
	4  44.64,46.20,47.52,48.97,50.61,52.07,53.41,54.95,56.55,57.93,
	5  59.33,60.93,62.45,63.81,65.29,66.89,68.33,69.72,71.27,72.82/

	DATA BETA/      0.00, 3.75, 6.47, 9.89,12.94,
	1	       16.28,19.40,22.72,25.87,29.17,
	2	       32.34,35.62,38.81,42.09,45.28,
	3	       48.55,51.74,55.01,58.21,61.48/

	DATA TURN /     0.00, 0.00, 0.00,
	1		3.70, 0.00, 0.00, 
	2		5.60, 3.70, 0.00, 
	3		5.00, 5.10, 3.70, 
	4		5.40, 6.90, 6.10/

	DATA FIVE/	0.00, 0.00, 0.00, 0.00,
	1		3.70, 0.00, 0.00, 0.00,
	2		6.50, 3.70, 0.00, 0.00,
	3		8.40, 6.30, 3.70, 0.00,
	4		6.60, 5.80, 5.50, 3.70,
	5		8.70, 9.10, 8.70, 6.30/
 
	DATA BB/        4.95, 6.20, 8.45,11.45,14.55,
	1	       17.65,20.90,24.15,27.45,30.65,
	2	       34.00,37.25,40.55/

	DATA A310/
     1   -0.042,-1.388, 1.834, 1.472,-1.616,-1.673, 3.013, 1.872,-1.334,
     2    4.464, 1.004, 2.097, 5.939,-2.228, 0.679, 7.427,-0.325,-2.287,
     3    8.917, 2.295, 0.069,10.402,-0.459, 2.268,11.884,-2.174,-0.812,
     4   13.360, 1.146,-2.011,14.844, 1.792, 1.456,16.334,-1.730, 1.525,
     5   17.815,-1.221,-1.965,19.289, 2.143,-0.899,20.777, 0.548, 2.244,
     6   22.266,-2.289, 0.155,23.743, 0.238,-2.309,25.222, 2.260, 0.581,
     7   26.712,-0.916, 2.104,28.197,-1.898,-1.281,29.673, 1.605,-1.704,
     8   31.158, 1.456, 1.818,32.647,-1.992, 1.096,34.126,-0.717,-2.193,
     9   35.607, 2.325,-0.409,37.096, 0.079, 2.304,38.572,-2.247,-0.349,
     1   40.053, 0.737,-2.226,41.558, 2.124, 1.009,43.051,-1.299, 1.839,
     2   44.986,-0.755,-1.696,46.921, 1.843, 0.220,48.857,-1.132, 1.471,
     3   50.792,-0.684,-1.726,52.728, 1.833, 0.295,54.663,-1.192, 1.424,
     4   56.598,-0.613,-1.752,58.534, 1.819, 0.370,60.469,-1.249, 1.373,
     5   62.404,-0.541,-1.776,64.340, 1.803, 0.444,66.275,-1.304, 1.321,
     6   68.210,-0.467,-1.797,70.146, 1.783, 0.518,72.081,-1.357, 1.266,
     7   74.017,-0.393,-1.814,75.952, 1.760, 0.591,77.887,-1.408, 1.210,
     8   79.823,-0.318,-1.829,81.758, 1.734, 0.662/

C******  SYMBOL TABLE  ************************************************

	ASTER  = '*'
	HELIX  = '@'
	H310   = '3'
	TRN    = 'T'
	TRN5   = '5'
	STRAND = 'B'
	LOOP   = '<'
	SHORT  = 'S'
	SHEET  = '#'

C**********************************************************************

C* Do initialization to all the arrays involved
	
	DO I =1, MAXLEN
		XA(I)=0
		YA(I)=0
		ZA(I)=0
		LI(I)=0
		DMIN(I)=0
		NMIN(I)=0
		DO J=1, MAXLEN
			DIST(I,J)=0
		ENDDO
	ENDDO

	DO I=1,100
		IFIRST(I)=0
		ILAST(I)=0
		CURV(I)=0
	ENDDO

	DO I=1,175
		FE(I)=0
		C3B(I)=0
		C3C(I)=0
	ENDDO
	
	DO I=1,ICP
		IHSTRT(I)=0
		IHEND(I)=0
		ISSTRT(I)=0
		ISEND(I)=0
		HRAD(I)=0
		HDIST(I)=0
		SRAD(I)=0
		SDIST(I)=0
		HTHA(I)=0
		HTHMAX(I)=0
		STHA(I)=0
		STHMAX(I)=0
		DO J=1,3
			HLXCOS(I,J)=0
			HLXPT(I,J)=0
			STRCOS(I,J)=0
			STRPT(I,J)=0
		ENDDO
	ENDDO

	DO I=1,30
		HLXWT(I)=0
		STRWT(I)=0
		SIGMA(I)=0
	ENDDO

	DO I=1,50
		DO J=1,50
			HA3(I,J)=0
		ENDDO
	ENDDO
 
	DO I=1,3
		NCNT(I)=0
		DO J=1,LINS
			IX(I,J)=0
			IY(I,J)=0
			BIT(I,J)=0
		ENDDO
	ENDDO

	DO I=1,100
		DO J=1,6
			PAIR(I,J)=0
		ENDDO
	ENDDO

	SEQNUM=0
	OFFSET=0
	
C****** Assign input,output file names, and open files.
C****** Unit #1 = input file ( BNL Protein Data Bank format)
C****** Unit #2 = output file for data to be printed
C****** Unit #3 = output file for PDB header information
C****** Unit #4 = output file for structure data (for use in other pgms)
C****** Unit #7 = CHECK.FIL = list of distances or standard deviations
C******  	  used in accepting or rejecting mask fits to the 
C******		  distance matrix.

	TZERO = CPUTIM(0.)

	TYPE 5001

	CALL DATE(PRDATE)
	CALL TIME(PRTIME)

	CALL ASKFIL(1,'  WHAT IS THE NAME OF THE INPUT FILE?')
		GNAME=FNAME
	CALL ASKFIL(2,'  WHAT IS THE NAME OF THE OUTPUT FILE?')

C******  DISTANCE MATCH LIMITS  ***************************************

	OPEN(UNIT=7,NAME='check.fil',  FORM='FORMATTED',STATUS='OLD')

4900	FORMAT(F8.3)
	READ(7,4900) CHKHA		!for AHELIX-alpha max rms error
C							for last atom
	READ(7,4900) CHKH3		!for AHELIX-310   max rms error
C							for last atom
	READ(7,4900) CHKTRN		!for SHARPTRN max sigma for mask
	READ(7,4900) CHKFIV		!for FIVETRN  max sigma for mask
	READ(7,4900) CHKSTR		!for BSTRAND  max error each atom
	READ(7,4900) CHKLP1		!for LOOPS    max loop neck
	READ(7,4900) CHKLP2		!for   "      max loop dimension
	READ(7,4900) CHKSHT		!for BSHEET   
	READ(7,4900) CHKSDT		!for   "      max separation
	READ(7,4900) CHKHSD		!for   "
	READ(7,4900) CHKSHD		!for   "
	READ(7,4900) CHKHHD		!for HLX_HLX
	READ(7,4900) SEPMAX
	READ(7,4900) SEPMIN
	READ(7,4900) AXDMAX
	READ(7,4900) AXDMIN
	READ(7,4900) PT1MAX
	READ(7,4900) PT1MIN
	READ(7,4900) PT2MAX
	READ(7,4900) PT2MIN

C******  XXXXXD = minimum distances from atoms in one secondary unit
C******  to the other.(i.e. not the distances to the helix axes.)

	CLOSE(UNIT=7)

C*********************************************************************

C******  Load array HA3 (distance mask) from CA data in A310 
C******  CA atoms 1-30 are alpha carbon atoms, atoms 31-50 are in the
C******  310 conformation.

	DO I=1,50
		DO J=I,50
			HA3(I,J)=SQRT((A310(1,J)-A310(1,I))**2+
	1				      (A310(2,J)-A310(2,I))**2+
	2				      (A310(3,J)-A310(3,I))**2)
		ENDDO
	ENDDO

C**********************************************************************

	WRITE(2,5000)
	WRITE(2,4999) PRDATE,PRTIME
4999	FORMAT(' RUN   :',5X,A9,10X,A8)
	WRITE(2,5000)
5000 	FORMAT(' --------------------------------------------------')

	WRITE(2,5001)
5001	FORMAT(/,' PROGRAM:     DEFINE_Structure',/,
	1	 ' AUTHORS:	F.M.Richards & C.E.Kundrot',/,
	2	 ' PURPOSE:	Specification of secondary and first',/,
	3	 '		level supersecondary structure from',/,
	4	 '		alpha carbon coordinate list.',/,
	5	 ' COPYWRITE (C) Yale University,  1988,  1994',/)

	WRITE(2,5000)

	WRITE(2,5025) CHKHA,CHKH3,CHKTRN,CHKFIV,CHKSTR,CHKLP1,CHKLP2,
	1		CHKSHT,CHKSDT,CHKHSD,CHKHSD,CHKHHD
5025	FORMAT(/,' PARAM : **CHECK DISTANCES**   (FROM CHECK.FIL)',//,
	1	 ' PARAM :  HA =',F7.3,' rms limit for fit to ahelix mask',/,
  	2	 ' PARAM :  H3 =',F7.3,' rms limit for fit to H310 mask',/, 
	3	 ' PARAM : TRN =',F7.3,' rms limit for fit to TRN mask',/,
	3	 ' PARAM : FIV =',F7.3,' rms limit for fit to FIV mask',/,
	4	 ' PARAM : STR =',F7.3,' rms limit for fit to STRAND mask',/,
	5	 ' PARAM : LP1 =',F7.3,' max end to end distance for loop',/,
	6	 ' PARAM : LP2 =',F7.3,' fraction of max loop dimension',/,
	7	 ' PARAM ; SHT =',F7.3,'                               	',/,
	8	 ' PARAM : SDT =',F7.3,'			       	',/,
	9	 ' PARAM : HSD =',F7.3,'    				',/,
	1	 ' PARAM : SHD =',F7.3,'				',/,
	2	 ' PARAM : HHD =',F7.3,'				',/)
	WRITE(2,5026) SEPMAX,SEPMIN,AXDMAX,AXDMIN,
	2		PT1MAX,PT1MIN,PT2MAX,PT2MIN
5026	FORMAT(  ' PARAM : SEPMAX =',F7.3,'				',/,
	1	 ' PARAM : SEPMIN =',F7.3,'				',/,
	2	 ' PARAM : AXDMAX =',F7.3,'				',/,
	3	 ' PARAM : AXDMIN =',F7.3,'				',/,
	4	 ' PARAM : PT1MAX =',F7.3,' 				',/,
	5	 ' PARAM : PT1MIN =',F7.3,'				',/,
	6	 ' PARAM : PT2MAX =',F7.3,'				',/,
	7	 ' PARAM : PT2MIN =',F7.3,/)
	WRITE(2,5000)
	WRITE(2,5005) GNAME
	WRITE(2,5010) FNAME
5005	FORMAT(/,' INFIL :           INPUT FILE NAME IS :  ',A30,/)
5010	FORMAT(  ' OUTFIL:    DATA  OUTPUT FILE NAME IS :  ',A30,/)
	GNAME=FNAME 
	CALL ASKFIL(3,'  WHAT IS THE NAME OF THE HEADER FILE?')
	CALL ASKFIL(4,'  NAME OF THE FINAL STRUCTURE FILE?')
	WRITE(2,5011) FNAME
5011	FORMAT(  ' OUTFIL: FINAL STRUCTURE FILE NAME IS :  ',A30,/)
	WRITE(2,5000)
	WRITE(4,5012) FNAME
5012	FORMAT('REMARK  ',A60)
	WRITE(4,5013)
5013	FORMAT('REMARK  FORMATS FOR THE FOLLOWING LISTS')
	WRITE(4,5014)
5014	FORMAT('REMARK  REMARK = (A8,A60) from MASK')
	WRITE(4,5021)
5021	FORMAT('REMARK  NCAATS = (A6,I5) from MASK')
	WRITE(4,5015)
5015	FORMAT('REMARK  CISPEP = RES1,RES2  (A6,2I5)  from LIST2ND'
	1      ' in SECSTR')
	WRITE(4,5022)
5022	FORMAT('REMARK  ELEMNT = SYMBOL,HP(X,1),APT(X,1),IS,IEE',/,
	1      '                 (A6,1X,A1,6F10.3,2I5) from AXISCHK')
	WRITE(4,5016)
5016	FORMAT('REMARK  SECND  = IFIRST,ILAST,LENGTH,SYMBOL,CURV',/,
	1      '                 (A6,3I5,A1,F5.0) from AXISCHK')
	WRITE(4,5017)
5017	FORMAT('REMARK  DIMENH = I,IHSTRT,IHEND,ALEN,HRAD,HDIST'
	1      ',HTHA,HTHMAX',/,
	2      'REMARK  DIMENS = I,ISSTRT,ISEND,SLEN,SRAD,SDIST'
	3      ',STHA,STHMAX'/,
	4      '         both  = (A6,3I5,F6.1,2F6.3,2F6.1) from'
	5      ' AXISCHK')
	WRITE(4,5018)
5018	FORMAT('REMARK  PAIRS  = ATOM1, ADJACENT ATOMS  (A6,6I5)',/,
	1      '                 from LISTBSH in SUPERSS')
	WRITE(4,5019)
5019	FORMAT('REMARK  HLXHLX = IS,IE,JS,JE,DMIN,IM,JM,ANGL,D6,D1',
	1      ',D2',/,
	2      '                 (A6,4I5,F6.1,2I5,4F6.0)  from',
	3      ' SUPERSS',/,
	4      'REMARK  STRHLX    same format as HLXHLX',/,
	5      'REMARK  HLXSTR    same format as HLXHLX',/,
	6      'REMARK  STRSTR	  same format as HLXHLX')

C******************************************************************

	LPDIST='N'
	TYPE *,' PRINT DISTANCE MATRIX ? Y/N'
	ACCEPT 5100, LPDIST
5100	FORMAT(A1)

	LPDIS2='N'
	TYPE *,' PRINT CHARACTER MATRIX ? Y/N'
	ACCEPT 5100, LPDIS2

	SSCALC='Y'
	TYPE *, ' CALCULATE AND PRINT SUPERSECONDARY STRUCTURE ? Y/N '
	ACCEPT 5100, SSCALC

	HHPRNT='Y'
	TYPE *,' PRINT HELIX-HELIX BOXES ? Y/N '
	ACCEPT 5100, HHPRNT

	HSPRNT='Y'
	TYPE *,' PRINT HELIX-STRAND AND STRAND-HELIX BOXES ? Y/N '
	ACCEPT 5100, HSPRNT

	SSPRNT='Y'
	TYPE *,' PRINT STRAND-STRAND BOXES ? Y/N '
	ACCEPT 5100, SSPRNT

C******************************************************************

C****** Load atom coordinates from input file
C****** Select alpha carbon atoms

100	CALL LOADCRD (LSTATM,MAXLEN,XA,YA,ZA,JMAX,IPEND,ATM,
	1			PRDATE,PRTIME)

		IF(IPEND.EQ.0) IPEND=JMAX

	WRITE(4,5110) JMAX
5110	FORMAT('NCAATS',I5)
	
C********************************************************************

C******	Load upper half of distance matrix DIST.

	DO I=1,JMAX
		DO J=I,JMAX
			TEMP=(XA(I)-XA(J))**2+(YA(I)-YA(J))**2+
	1			(ZA(I)-ZA(J))**2
			DIST(I,J)=SQRT(TEMP)
		ENDDO
	ENDDO

C****** Initialize array DIST2 with blanks, except for '-' on diagonal

	DO I=1,MAXLEN
		DO J=1,MAXLEN
			DIST2(I,J)=' '
			IF(J.EQ.I) DIST2(I,J)=INT(11)
		ENDDO
	ENDDO

	TSTART = CPUTIM(TZERO)
	TLOAD=TSTART

C*********************************************************************

C****** Look for cis peptide bonds and check that all trans bonds are 
C****** within reasonable distance limits.

150	CALL CISPEP (DIST,DIST2,MAXLEN,JMAX)
	
C*********************************************************************
C****** Check helix mask (ALPHA) against each line. Enter symbol,
C****** HELIX, in each consecutive position where the difference 
C****** between DIST and the mask ALPHA is less than CHKHLX.

200	CALL AHELIX (DIST,DIST2,MAXLEN,JMAX,ALPHA,HA3,HELIX,H310,
	1		CHKHA,CHKH3,HLXWT)

C*********************************************************************

C****** Check in non-helical segments for tight beta turns using the 
C****** mask array TURN. The rms match must be less than or equal to
C****** CHKTRN.  

300	CALL SHARPTRN( DIST,DIST2,MAXLEN,JMAX,TURN,TRN,HELIX,H310,CHKTRN)

C********************************************************************* 

C******  Check in non-helical and non-tight turn segments for turns
C******  involving 5 residues using the mask array FIVE. The rms match
C******  must be less than or equal to CHKFIV.

350	CALL FIVETRN( DIST,DIST2,MAXLEN,JMAX,FIVE,TRN5,TRN,HELIX,
	1		H310,CHKFIV)

C*********************************************************************

C******  Check for extended chains with the mask BETA.
C******  This mask has distances that are equivalent to those found in
C******  parallel sheets, i.e. somewhat shorter than a fully extended
C******  chain. The symbol,STRAND, is entered in each consecutive 
C******  position where the difference between DIST and the mask BETA
C******  is less than CHKSTR.

400	CALL BSTRAND (DIST,DIST2,MAXLEN,JMAX,BETA,HELIX,H310,TRN,TRN5,
	1		STRAND,CHKSTR,STRWT)

C*******************************************************************

C****** Write to output file the beginning and ending alpha carbon
C****** numbers and length for each type of secondary structure:
C****** helices, turns, and strands.

500	CALL LIST2ND (DIST,DIST2,MAXLEN,IPEND,JMAX,IFIRST,ILAST,
	1		SYMBOL,HELIX,H310,TRN,TRN5,STRAND)

c	TEND=CPUTIM(TZERO)
c	T2A=TEND - TSTART
c	TSTART=TEND

C*******************************************************************

C******  Check for curvature in helix and strand segments. Link 
C******  helices or strands that seem to be reasonably uniformly bent.

	CALL AXISCHK(DIST,DIST2,MAXLEN,JMAX,IPEND,
	1	IFIRST,ILAST,SYMBOL,CURV,ABC,INT,
	2	CHKHHD,CHKHSD,CHKSHD,XA,YA,ZA,
	3	HELIX,H310,STRAND,TRN,TRN5,
	4	HLXCOS,HLXPT,STRCOS,STRPT,IHSTRT,IHEND,ISSTRT,
	5	ISEND,HRAD,HDIST,SRAD,SDIST,ICP,MAXH,MAXS,
	6	HTHA,HTHMAX,STHA,STHMAX)

c	TEND=CPUTIM(TZERO)
c	T3A=TEND - TSTART
c	TSTART=TEND

C********************************************************************

C******  Check for beta strand pairs either parallel or antiparallel

	IF(SSCALC.EQ.'N') GOTO 700
        IF(SSCALC.EQ.'n') GOTO 700


600	CALL BSHEET (DIST,DIST2,MAXLEN,JMAX,STRAND,SHORT,SHEET,BB,
	1		HELIX,H310,TRN,TRN5,CHKSHT,CHKSDT)

C********************************************************************

C******  Check for omega LOOPS, lengths of chain between 6 and 16 
C******  residues inclusive as described by Leszczynski and Rose,
C******  Science  234: 849-855 (1986). Output to DIST2 does not 
C******  overwrite any other symbols previously entered.

700	CALL LOOPS (DIST,DIST2,MAXLEN,JMAX,HELIX,H310,TRN,TRN5,
	1		STRAND,LOOP,CHKLP1,CHKLP2)

C********************************************************************

C******  Write out paired residues in the beta sheet

	IF(SSCALC.EQ.'N') GOTO 1100

800	CALL LISTBSH (DIST,DIST2,MAXLEN,IPEND,JMAX,
	1	IFIRST,ILAST,SYMBOL,PAIR,HELIX,H310,TRN,TRN5,
	2	STRAND,LOOP,SHORT,SHEET)

c	TEND=CPUTIM(TZERO)
c	T4A=TEND - TSTART
c	TSTART=TEND

C********************************************************************

C******	  Determine tertiary structure relations. Angles and distances 
C******   between the axes of the various segments. 
C******   Helix-helix,  helix-strand, strand-helix.

	DO I=1,3
		NCNT(I)=0
	ENDDO
900	CALL SUPERSS(DIST,DIST2,MAXLEN,JMAX,IPEND,
	1	IFIRST,ILAST,SYMBOL,ABC,INT,
	2	CHKHHD,CHKHSD,CHKSHD,XA,YA,ZA,
	3	LINS,IX,IY,BIT,COLR,NCNT,HELIX,H310,STRAND,
	4	SEPMAX,SEPMIN,AXDMAX,AXDMIN,
	5	PT1MAX,PT1MIN,PT2MAX,PT2MIN,
	6	HHPRNT,HSPRNT,SSPRNT,
	7	HLXCOS,HLXPT,STRCOS,STRPT,IHSTRT,IHEND,ISSTRT,
	8	ISEND,HRAD,HDIST,SRAD,SDIST,ICP,MAXH,MAXS)

c	TEND=CPUTIM(TZERO)
c	T5A=TEND - TSTART
c	TSTART=TEND

C*************************************************************

1000	IF(LPDIST.EQ.'N') GOTO 1100
	IF(LPDIST.EQ.'n') GOTO 1100
	CALL PRDIST(JMAX,ASTER,DIST,LI,MAXLEN)

1100	IF(LPDIS2.EQ.'N') GOTO 1150
	IF(LPDIS2.EQ.'n') GOTO 1150
	CALL PRDIS2(JMAX,ASTER,DIST2,MAXLEN)

1150	IF(NCNT(1).GT.LINS.OR.NCNT(2).GT.LINS.OR.NCNT(3).GT.LINS) THEN
		MAXCNT=MAX(NCNT(1),NCNT(2),NCNT(3))
		TYPE *,' ERROR FOR MXMAP: NCNT.GE.LINS',MAXCNT
		GOTO 9000
	ENDIF

1200	CONTINUE

C******  Enter the name for the PS file
5130	FORMAT(A30)

	TYPE *, 'WHAT IS THE NAME OF PS FILE (xxxxx.ps) ?'
	ACCEPT 5130, DEVICE

	OPEN(UNIT=8,NAME=DEVICE, carriagecontrol='list',
	1	FORM='FORMATTED',STATUS='NEW')

	write (8,1105)
	write (8,1106)
	write (8,1107)
	write (8,1110)
1105	FORMAT('%!PS-Adobe-3.0')
1106	FORMAT('25 25 translate')
1107	FORMAT('0.2 0.2 scale')
1110    FORMAT('.25 setlinewidth')
	
	TYPE *, SEPMAX,SEPMIN,AXDMAX,AXDMIN
	TYPE *, PT1MAX,PT1MIN,PT2MAX,PT2MIN

	NQ=1
	CALL MXMAP(DIST2,MAXLEN,JMAX,NQ,GNAME,HELIX,H310,
	1	STRAND,TRN,TRN5,LOOP,SHEET,SHORT,
	2	LINS,IX,IY,BIT,COLR,NCNT,DEVICE,
	3	SEPMAX,SEPMIN,AXDMAX,AXDMIN,PT1MAX,PT1MIN,
	4	PT2MAX,PT2MIN)

	IF(JMAX.GT.260) THEN
		DO NQ=2,3
			CALL MXMAP(DIST2,MAXLEN,JMAX,NQ,GNAME,HELIX,H310,
	1		STRAND,TRN,TRN5,LOOP,SHEET,SHORT,
	2		LINS,IX,IY,BIT,COLR,NCNT,DEVICE,
	3		SEPMAX,SEPMIN,AXDMAX,AXDMIN,PT1MAX,PT1MIN,
	4		PT2MAX,PT2MIN)
		ENDDO
	ENDIF
	
	CLOSE (unit = 8)
C*	CALL UGCLOS( 'NOCLEAR' )
9000	continue
c	TEND=CPUTIM(TZERO)
c	TDISP=TEND - TSTART

C*******************************************************************

c	WRITE(2,5139)
c5139	FORMAT(' RUN TIMES IN SECONDS',/)

c	WRITE(2,5140) TLOAD,T2A,T3A,T4A,T5A,TDISP

c5140	FORMAT(/,'    TLOAD =',F10.2,'  THRU LOADCRD',/,
c	1	 '    T2A   =',F10.2,'  THRU LIST2ND',/,
c	2	 '    T3A   =',F10.2,'  THRU AXISCHK',/,
c	3	 '    T4A   =',F10.2,'  THRU LISTBSH',/,
c	4	 '    T5A   =',F10.2,'  THRU SUPERSS',/,
c	5	 '    TDISP =',F10.2,'  THRU ALL DISPLAYS',/)

	STOP

	END

C********************************************************************
C=============       SUBROUTINE LIBRARY        =========================
C=================        SSUTIL           =============================
C================= (updated on 8 Mar 88)   =============================
C=======================================================================
C		CPUTIM			PRDIST
C		(DATE from system)	PRDIS2
C		(TIME from system)	LOADCRD
C		ASKFIL
C=======================================================================

C=======================  CPUTIM  ======================================

	FUNCTION CPUTIM (SECS)
C
C Rerurns number of seconds of CPU time since SECS
C
c	INCLUDE		'($JPIDEF)'
c	INCLUDE		'($SSDEF)'

c	INTEGER*4	SYS$GETJPI

c	INTEGER*2	ITMLST	(8)
c	INTEGER*4	JTMLST	(4)
c	EQUIVALENCE	(ITMLST,JTMLST)
c	INTEGER*4	NCPU

c	ITMLST(1) = 4			!length in bytes
c	ITMLST(2) = JPI$_CPUTIM		!item code
c	JTMLST(2) = %LOC (NCPU)
c	JTMLST(3) = 0			!no return list
c	JTMLST(4) = 0			!end of list

C $GETJPI Returns accumulated CPU time in 10-ms tics (1 centisecond)
C
c	ICODE = SYS$GETJPI (, , ,ITMLST,,,)

c	IF (ICODE .NE. SS$_NORMAL) CALL EXIT (ICODE)
c	CPUTIM = FLOAT(NCPU)/100. - SECS
	cputim = 0.0
	RETURN 
	END

C=======================   ASKFIL   =================================

	SUBROUTINE ASKFIL(LUNIT,CHAR)
C
	COMMON/FILES/ FNAME
C
	CHARACTER*(*)  CHAR
	CHARACTER*30   FNAME
C
C		THE SIZE OF THE MESSAGE BUFFER IS DEFINED AT RUN TIME
C		SINCE IT IS USED AS A PASSED LENGTH STRING. FNAME IS 
C		A LOCALLY DEFINED CONSTANT CURRENTLY SET AT 30.
C
	TYPE 10, CHAR
10	FORMAT(5X,A,'---->'$)
c	TYPE*,'--------------------------------'
C
	ACCEPT 100, FNAME
100	FORMAT(A30)
 	TYPE*,' '
C
	OPEN(UNIT=LUNIT,NAME=FNAME, STATUS='UNKNOWN')
C
	RETURN
C
	END

C========================   PRDIST   ===============================

C******  Call in program anywhere that a hard copy of actual distance
C******  matrix DIST is needed. Printed numbers are integers 
C******  representing truncated DIST(I<J)*10.0

	SUBROUTINE PRDIST(JMAX,ASTER,DISTT,LI,MAXLEN)

	CHARACTER*1 ASTER

	DIMENSION DISTT(MAXLEN,MAXLEN),LI(MAXLEN)

	INTEGER = OFFSET

	OFFSET=0

2000	IS=1
	JS=IS+OFFSET
	JE=JS+24
	IF(JE.GT.JMAX) JE=JMAX
	IE = JE

	WRITE(2,5299)
5299	FORMAT('1')

	WRITE(2,5300) (I,I=JS+4,JE,5)
5300	FORMAT(11X,10I20)

	WRITE(2,5310) (ASTER,I=1,5)
5310	FORMAT(11X,5(19X,A1),/)

	DO 2100 I=1,IE

		DO 2050 J=JS,JE
2050			LI(J)=DISTT(I,J)*10.

		WRITE(2,5320) I, (LI(L),L=JS,JE),I
5320		FORMAT(1X,I5,5X,25I4,I10)
2100	CONTINUE

	OFFSET = OFFSET + 25
	IF(OFFSET.GE.JMAX) RETURN

	GOTO 2000

	END

C======================   PRDIS2   ====================================

C******  Prints the matrix DIST2, the character representation derived
C******  from DIST through the mask checks.

	SUBROUTINE PRDIS2(JMAX,ASTER,DIST2,MAXLEN)

	CHARACTER*1 ASTER

	BYTE    DIST2(MAXLEN,MAXLEN)
	BYTE	LINE(500)
	
	JINCR=100	! Change depending on printer page width
	JLS=0		! See also FORMAT  5320
	DO IL=0,10
		IOFSET=IL*50
		IF(IOFSET.GT.JMAX) RETURN
		IS=IOFSET+1
		IE=IS+49
		DO JL=JLS,3
			JOFSET=JL*JINCR 
			JS=JOFSET + 1
			JE=JS+JINCR-1
			IF(JOFSET.GT.JMAX) GOTO 1000
			WRITE(2,5299)
5299			FORMAT('1')  !Print page indicator

			WRITE(2,5300) (I,I=JS+9,JE,10)
5300			FORMAT(6X,12I10)  !Print residue number 

			WRITE(2,5310) (ASTER,I=1,JINCR/10)
5310			FORMAT(6X,12(9X,A1),/)  !Print asterisks 

			DO  I=IS,IE
				DO L=JS,JE
					LINE(L)=' '
				ENDDO
				DO L=JS+9,JE,10
					LINE(L)='!'
				ENDDO
				IF(JMOD(I,10).EQ.0) THEN !- on diagonal
					DO L=JS,JE
						LINE(L)='-'
					ENDDO
				ENDIF
C******  JJE = last blank position starting from current page, right margin
				JJE=JE
				DO J=JE,JS,-1
					IF(DIST2(I,J).EQ.' ') THEN
						JJE=J
						IF(JJE.EQ.JS) JJE=JS-1
					ELSE
						GOTO 950
					ENDIF
				ENDDO
950				CONTINUE
C******  IIS = diagonal position on line or left margin if away
C******        from diagonal
				IIS=I
				IF(IE.LT.JS) IIS=JS
				IF(JJE.LE.JS) GOTO 960

				DO J=IIS,JJE
					LINE(J)=DIST2(I,J)
				ENDDO
960		IF (IIS.GT.1)	LINE(IIS-1) = ' '
		IF (IIS.GT.2)	LINE(IIS-2) = ' '
				LINE(JJE+1) = ' '
				LINE(JJE+2) = ' '	

				IF(JMAX.LT.JE.AND.JMAX.GT.JS)
	1					LINE(JMAX+1)='*'
				IF(I.EQ.JMAX+1) THEN
					DO L=JS,JE
						LINE(L)='*'
					ENDDO
				ENDIF


				WRITE(2,5320) I,(LINE(L),L=JS,JE)
5320				FORMAT(1X,I4,1X,<JINCR>A1)
			ENDDO
		ENDDO

1000		IF(IE.GT.JS) THEN
			JLS=JLS+1
		ENDIF

	ENDDO

	RETURN
		
	END

C========================   LOADCRD   ==================================

C****** Subroutine to load coordinates from input tape and set
C******	up alpha carbon arrays.

	SUBROUTINE LOADCRD (LSTATM,MAXLEN,XA,YA,ZA,JMAX,IPEND,ATM,
	1			PRDATE,PRTIME)
	
	DIMENSION XA(MAXLEN),YA(MAXLEN),ZA(MAXLEN)

	INTEGER	     SEQNUM,RESCNT

	CHARACTER*9  PRDATE
	CHARACTER*8  PRTIME
	CHARACTER*6  KEYWRD
	CHARACTER*4  ATM(MAXLEN),AATM

	BYTE         BLINE(74),ASTER

	DATA ASTER /'*'/

	WRITE(3,4900) (ASTER,I=1,60)
4900	FORMAT(1X,60A1)
	WRITE(3,4910) (ASTER,I=1,10),(ASTER,I=1,10),
	1	      (ASTER,I=1,10),PRDATE,PRTIME,(ASTER,I=1,10)
4910	FORMAT(1X,10A1,' Protein Data Bank Header File ',10A1,/,
	1      1X,10A1,' Printed on: ',A9,4X,A8,10A1)     
	WRITE(3,4920) (ASTER,I=1,60)
4920	FORMAT(1X,60A1,//)

100	READ(1,5001,END=5200) KEYWRD, BLINE
5001	FORMAT(A6,74A1)
5002	FORMAT(1X,A6,74A1)
	IF(KEYWRD.NE.'ATOM  ') THEN
		WRITE(3,5002) KEYWRD,BLINE
		GOTO 100
	ELSE
		GO TO 150
	ENDIF

5200	TYPE*,' ERROR IN FILE FORMAT. NO "ATOM" FOUND'
	STOP

150	LSTATM=0
	LAST=0
	NEXT=0
	J=0
	BACKSPACE 1
	DO 200 I=1,10000
		READ(1,5000,END=300) KEYWRD,AATM,RES3,SEQNUM,X,Y,Z
5000		FORMAT(A6,7X,2A4,1X,I4,4X,3F8.3)
		IF(KEYWRD.NE.'ATOM  ') GOTO 300
		LSTATM=LSTATM+1
		IF(I.EQ.1) THEN		! get correct residue number 
			RESCNT=SEQNUM-1	! if file does not start with
			J=RESCNT	! SEQNUM=1
			IF(J.LT.0) J=0	! correction if SEQNUM=0
		ENDIF
C****** Select alpha carbon atoms

		IF(AATM.EQ.'CA  ') THEN
			J=J+1
			XA(J)=X
			YA(J)=Y
			ZA(J)=Z
			ATM(J)=AATM
		ENDIF
		IF(AATM.EQ.'FE  ') IPEND=J
		IF(RES3.EQ.'HEM ') THEN
			J=J+1
			XA(J)=X
			YA(J)=Y
			ZA(J)=Z
			ATM(J)=AATM
		ENDIF
200	CONTINUE

300	JMAX = J
	JCA=J-RESCNT
C******	all arrays have now been loaded from input file

	WRITE(2,5104) LSTATM
5104	FORMAT(/,' ATOMS : TOTAL NUMBER OF PROTEIN ATOMS',/,
	1	 ' ATOMS :              IN INPUT LIST IS:',I5)

	WRITE(2,5290) JCA
5290	FORMAT(  ' ATOMS :    NUMBER OF C ALPHA ATOMS IS:',I5,/)
	WRITE(2,5003)
5003	FORMAT(' --------------------------------------------------')

	RETURN

	END

C======================================================================

C=============       SUBROUTINE LIBRARY        =========================
C=================        SECSTR           =============================
C================= (Updated on 29 Dec 87)  =============================
C=======================================================================
C		CISPEP			BSTRAND
C		AHELIX			LOOPS
C		SHARPTRN		LIST2ND
C		FIVETRN
C=======================================================================

C=========================   CISPEP   ==================================

	SUBROUTINE CISPEP (DIST,DIST2,MAXLEN,JMAX)

	DIMENSION    DIST(MAXLEN,MAXLEN)

	BYTE         DIST2(MAXLEN,MAXLEN),HELIX,H310,CHK1,CHK2

	DO I=1,JMAX-1
		IF(DIST(I,I+1).GE.2.7.AND.
	1	   DIST(I,I+1).LE.3.1) DIST2(I,I)='C'
	ENDDO

	DO I=1,JMAX-1
		IF(DIST2(I,I).NE.'C') THEN
			IF(DIST(I,I+1).LT.3.6.OR.DIST(I,I+1).GT.4.0) 
	1			THEN
				WRITE(2,5000) I,I+1,DIST(I,I+1)
5000				FORMAT(' BAD DISTANCE CA ATOMS',
	1				I4,' -',I4,' =',F8.2)
			ENDIF
		ELSEIF(DIST2(I,I).EQ.'C') THEN
			WRITE(4,5100)I,I+1
5100			FORMAT('CISPEP',2I5)
		ENDIF
	ENDDO

	RETURN

	END

C===========================     AHELIX     ============================

C******  Match of consectuve lines, in the form of a triangle, to the
C******  helix mask ALPHA is checked for agreement within the error limit
C******  specified by CHKHLX. The total number of misses, NO, is 
C******  accumulated except in the first line. When NO is greater than or
C******  equal to the number of lines tested, the helix is terminated.

	SUBROUTINE AHELIX (DIST,DIST2,MAXLEN,JMAX,ALPHA,HA3,
	1			HELIX,H310,CHKHA,CHKH3,HLXWT)

	DIMENSION    DIST(MAXLEN,MAXLEN),ALPHA(50)
	DIMENSION    HA3(50,50),HLXWT(30)

	BYTE         DIST2(MAXLEN,MAXLEN),HELIX,H310,CHK1,CHK2

	CHKSQA=CHKHA**2
	CHKSQ3=CHKH3**2
	I=1

1000	IF(I.GE.JMAX) RETURN
1100	CONTINUE
	JE=I+49
		IF(JE.GT.JMAX) JE=JMAX
	LH=0
	DO J=I+1,JE
		LH=LH+1
		DDK=0
		NDDK=0
		DO K= I,J-1
			JN=J-K+1
C******  Omit diagonal and first off-diagonal positions from sum
			IF(JN.GT.2.AND.K.LT.J-1) THEN
				DK=(DIST(K,J)-ALPHA(JN))**2
C******  Weight matrix elements for curvature according to length of helix
				DK=DK*((1/(1+HLXWT(J-K)))**2)	
				DDK=DDK+DK
				NDDK=NDDK+1		
			ENDIF
			DIST2(K,J)=HELIX
		ENDDO
		IF(NDDK.EQ.0) GOTO 1150
C******  Check for rms match in each successive J line.
		IF(DDK/NDDK.GT.CHKSQA) THEN
			DO K=I,J-1		!Erase vertical line J
				DIST2(K,J)=' '
			ENDDO
			LH=LH-1
			GOTO 1200
		ENDIF
1150	ENDDO

C******  Alpha helix terminated. Check to see if there is a 310 helix
C******  at the C terminus.

1200	JJS=J
	JJE=JJS+19
C******  Maximum length of alpha helix = JJS-I
	DO J=JJS,JJE
		DDK=0
		NDDK=0
		DO K=I,J-2
			DDK=DDK+(DIST(K,J)-HA3(31+K-JJS,31+J-JJS))**2
			NDDK=NDDK+1
			DIST2(K,J)=H310
		ENDDO
		IF(DDK/NDDK.GT.CHKSQ3) THEN
			DO K=I,J-2
				DIST2(K,J)=' '
			ENDDO
			GOTO 1300	! Terminate 310 helix
		ENDIF
	ENDDO
1300	JJE=J-1
	IF(JJE.LT.JJS)  GOTO 1500

C******  Fill in bottom of triangle with H310 if any 310 helix was 
C******  added to the end of the alpha helix.
	KS=I+LH
	DO J=JJS,JJE
		LH=LH+1
		DO K=KS,J-1
			DIST2(K,J)=H310
		ENDDO
	ENDDO
1500	CONTINUE

C******  Save all helices longer than 4 residues

	IF(LH.GT.4) THEN
		I=I+LH
		GOTO 1000

C******  Erase helical bits that are too short

	ELSE
		DO J=I+1,I+LH
			DO K=I,J-1
				DIST2(K,J)=' '
			ENDDO
		ENDDO
		I=I+1
		GOTO 1000
	ENDIF

	END

C========================     SHARPTRN     =============================

C****** Check for turns using the mask array TURN. 

	SUBROUTINE SHARPTRN (DIST,DIST2,MAXLEN,JMAX,TURN,TRN,
	1			HELIX,H310,CHKTRN)

	DIMENSION DIST(MAXLEN,MAXLEN),TURN(3,5)

	BYTE      DIST2(MAXLEN,MAXLEN),TRN, HELIX, H310

	ICHK=1
1000	IF(ICHK.GT.JMAX-4) RETURN

C******  Determine fit of turn mask to lines ICHK,+1,+2
	
	SCORE = 0.0
	DO  L=ICHK,ICHK+2
		DO M=L+1,ICHK+4
			DX2=(DIST(L,M) - TURN(L-ICHK+1,M-ICHK+1))**2
			SCORE = SCORE + DX2
		ENDDO
	ENDDO
	SCORE=SQRT(SCORE/9.)

	IF(SCORE.GT.CHKTRN) THEN
		ICHK=ICHK+1
		GOTO 1000
	ENDIF
C******  If match is better than CHKTRN, see if turn should be put in.
C******  Check whether or not current line has a HELIX symbol
	IF(DIST2(ICHK,ICHK+1).EQ.HELIX.OR.
	1  DIST2(ICHK,ICHK+1).EQ.H310) THEN
		N=1
C****** Check ahead for length of helix beyond ICHK
		DO I=ICHK+1,JMAX
			IF(DIST2(I,I+1).EQ.HELIX.OR.
	1		   DIST2(I,I+1).EQ.H310) THEN
				N=N+1
			ELSE
				GOTO 1100
			ENDIF
		ENDDO
C******  Start next check beyond helix
1100		ICHK=ICHK+N
		N=0
		GOTO 1000

	ELSEIF(DIST2(ICHK,ICHK+1).NE.HELIX.AND.
	1      DIST2(ICHK,ICHK+1).NE.H310) THEN
C******  Check ahead to see if space is clear to insert turn or
C******  whether a new potential helix is starting.
		N=0
		DO I=ICHK+1,ICHK+3
			IF(DIST2(I,I+1).EQ.HELIX.OR.
	1		   DIST2(I,I+1).EQ.H310) N=N+1
		ENDDO
C******  If N=0 then OK to insert turn
		IF(N.LE.0) THEN
			GOTO 1300	! Insert turn
		ELSE		! i.e. N>0
C******  If N>0 then next helix is starting just ahead.
			N=0
			ICHK=ICHK+1
			GOTO 1000
		ENDIF
	ENDIF

1300	CONTINUE
	DO L=ICHK,ICHK+2
		DO M=L+1,ICHK+3
			DIST2(L,M) = TRN
		ENDDO
	ENDDO

	ICHK=ICHK+3
	N=0
	GOTO 1000

	END


C========================     FIVETRN     =============================

C****** Check for turns using the mask array FIVE. 

	SUBROUTINE FIVETRN (DIST,DIST2,MAXLEN,JMAX,FIVE,TRN5,TRN,
	1			HELIX,H310,CHKFIV)

	DIMENSION DIST(MAXLEN,MAXLEN),FIVE(4,6)

	BYTE      DIST2(MAXLEN,MAXLEN),TRN, TRN5, HELIX, H310

	ICHK=1
1000	IF(ICHK.GT.JMAX-5) RETURN

C******  Determine fit of turn mask to lines ICHK,+1,+2,+3
	
	SCORE = 0.0
	DO  L=ICHK,ICHK+3
		DO M=L+1,ICHK+5
			DX2=(DIST(L,M) - FIVE(L-ICHK+1,M-ICHK+1))**2
			SCORE = SCORE + DX2
		ENDDO
	ENDDO
	SCORE=SQRT(SCORE/14.)

	IF(SCORE.GT.CHKFIV) THEN
		ICHK=ICHK+1
		GOTO 1000
	ENDIF
C******  If match is better than CHKFIV, see if turn should be put in.
C******  Check whether or not current line has a HELIX or TRN symbol
	IF(DIST2(ICHK,ICHK+1).EQ.HELIX.OR.
	1  DIST2(ICHK,ICHK+1).EQ.H310 .OR. 
	2  DIST2(ICHK,ICHK+1).EQ.TRN) THEN
		N=1
C****** Check ahead for length of HELIX or TRN beyond ICHK
		DO I=ICHK+1,JMAX
			IF(DIST2(I,I+1).EQ.HELIX.OR.
	1		   DIST2(I,I+1).EQ.H310 .OR. 
	2		   DIST2(I,I+1).EQ.TRN) THEN
				N=N+1
			ELSE
				GOTO 1100
			ENDIF
		ENDDO
C******  Start next check beyond HELIX or TRN
1100		ICHK=ICHK+N
		N=0
		GOTO 1000

	ELSE
C******  Check ahead to see if space is clear to insert turn or
C******  whether a new potential helix is starting.
		N=0
		DO I=ICHK+1,ICHK+4
			IF(DIST2(I,I+1).EQ.HELIX.OR.
	1		   DIST2(I,I+1).EQ.H310 .OR.
	2		   DIST2(I,I+1).EQ.TRN) N=N+1
		ENDDO
C******  If N=0 then OK to insert turn
		IF(N.LE.0) THEN
			GOTO 1300	! Insert turn
		ELSE		! i.e. N>0
C******  If N>0 then next helix is starting just ahead.
			N=0
			ICHK=ICHK+1
			GOTO 1000
		ENDIF
	ENDIF

1300	CONTINUE		! Insert "five" turn symbol
	DO L=ICHK,ICHK+3
		DO M=L+1,ICHK+4
			DIST2(L,M) = TRN5
		ENDDO
	ENDDO

	ICHK=ICHK+4
	N=0
	GOTO 1000

	END

C=============================     BSTRAND     ========================

C******  Match of consectuve lines, in the form of a triangle, to the
C******  strand mask BETA is checked. Total number of misses, NO, is 
C******  accumulated. When NO is greater than or equal to one half of 
C******  the number of lines tested, the strand is terminated. If the 
C******  strand runs into a helix or a turn before such a termination 
C******  then the helix or turn is given preference and the strand is 
C******  terminated before the first residue of the helix or turn.

	SUBROUTINE BSTRAND (DIST,DIST2,MAXLEN,JMAX,BETA,
	1		HELIX,H310,TRN,TRN5,STRAND,CHKSTR,STRWT)

	DIMENSION DIST(MAXLEN,MAXLEN),BETA(20),STRWT(30)

	BYTE      DIST2(MAXLEN,MAXLEN)
	BYTE      HELIX,H310,TRN,TRN5,STRAND

	CHKSQ=CHKSTR**2
	I=1
1000	IF(I.GE.JMAX) RETURN

C****** Do not check positions already identified as HELIX, TRN, or TRN5.

		IF(DIST2(I,I+1).EQ.HELIX.OR.
	1	   DIST2(I,I+1).EQ.H310 .OR.
	2	   DIST2(I,I+1).EQ.TRN  .OR.
	3	   DIST2(I,I+1).EQ.TRN5) THEN
			I=I+1
			GOTO 1000
		ENDIF

1100	NO=0
	JE=I+49
		IF(JE.GT.JMAX) JE=JMAX
	LH=0
	DO J=I+1,JE
		LH=LH+1
		DDK=0
		NDDK=0
		DO K= I,J-1
C******  If the strand runs into a helix or turn, these get preference
			IF(DIST2(K,J).EQ.HELIX.OR.
	1		   DIST2(K,J).EQ.H310 .OR.
	2		   DIST2(K,J).EQ.TRN  .OR.
	3		   DIST2(K,J).EQ.TRN5) GOTO 1200
			JN=J-K+1
C******  Omit diagonal and first off-diagonal positions from sum
			IF(JN.GT.2.AND.K.LT.J-1) THEN
				DK=(DIST(K,J)-BETA(JN))**2
				DK=DK*((1/(1+STRWT(J-K)))**2)
				DDK=DDK+DK
				NDDK=NDDK+1
			ENDIF
			DIST2(K,J)=STRAND
		ENDDO
		IF(NDDK.EQ.0) GOTO 1150
C******  Check for rms match in each successive J line.
		IF(DDK/NDDK.GT.CHKSQ) THEN
			DO K=I,J-1		!Erase vertical line J
				DIST2(K,J)=' '
			ENDDO
			LH=LH-1
			GOTO 1500
		ENDIF
1150	ENDDO
	GOTO 1500
1200	DO M=I,K-1
		DIST2(M,J)=' '
	ENDDO
	LH=LH-1
1500	IF(LH.GT.2) THEN
		I=I+LH
		GOTO 1000
	ELSE
		IF(J.GE.JMAX) THEN
			I=I+1
			GOTO 1000
		ENDIF
		DO J=I+1,I+LH
			DO K=I,J-1
				DIST2(K,J)=' '
			ENDDO
		ENDDO
		I=I+1
		GOTO 1000
	ENDIF

	END

C==========================   LOOPS   ===============================

C******  This subroutine identifies the omega loops as defined by 
C******  Leszczynski and Rose, Science 234:849-855 (1986) with regard to
C******  end to end distance constraints and maximum loop dimension. No
C******  account is taken of residue packing in this identification.

	SUBROUTINE LOOPS(DIST,DIST2,MAXLEN,JMAX,HELIX,H310,TRN,TRN5,
	1		STRAND,LOOP,CHKLP1,CHKLP2)

	DIMENSION DIST(MAXLEN,MAXLEN)

	BYTE  DIST2(MAXLEN,MAXLEN)
	BYTE  HELIX,H310,TRN,TRN5,STRAND,LOOP

	LOGICAL  LP(16)

C-----------------------------------------------------------------------

	WRITE(2,4910)
4910	FORMAT(/, ' SECND : POSSIBLE OMEGA  LOOPS',/,
	1	  ' SECND :   FIRST    LAST',/,
	2	  ' SECND : RESIDUE RESIDUE LENGTH',/)

	I=1
1000	IF(I.GE.JMAX) THEN
		WRITE(2,5001)
5001	FORMAT(' ---------------------------------------------------')
		RETURN
	ENDIF

C******  Start a loop check region if the current line I is not in a 
C******  helix or in a beta strand of total length equal to or greater
C******  than 5 residues.

	IF(DIST2(I,I+1).EQ.HELIX.OR.
	1  DIST2(I,I+1).EQ.H310) THEN
		I=I+1
		GOTO 1000
	ENDIF
	IF (DIST2(I,I+1).EQ.STRAND) THEN
		N=0
		DO J=I+1,JMAX
			IF(DIST2(I,J).EQ.STRAND) N=N+1
			IF(DIST2(I,J).NE.STRAND) GOTO 1100
		ENDDO
1100		IF(N.GE.4) THEN
			I=I+N
			GOTO 1000
		ENDIF
	ENDIF

C****** Start at this value of I to check ahead for the lengths of 
C****** possible loops. Stop at a helix or at a strand section that is 
C****** 5 residues in total length or longer.

	DO K=I+1,JMAX
		IF(DIST2(K,K+1).EQ.HELIX.OR.
	1	   DIST2(K,K+1).EQ.H310) GOTO 1300  ! stop search at K
		IF(DIST2(K,K+1).EQ.STRAND) THEN
			N=0
			DO J=K+1,JMAX
				IF(DIST2(K,J).EQ.STRAND) N=N+1
				IF(DIST2(K,J).NE.STRAND) GOTO 1200
			ENDDO
1200			IF(N.GE.4) GOTO 1300	     ! stop search at K
		ENDIF
	ENDDO
	KMAX=K-1
	GOTO 1400
1300	KMAX=K

C****** Start search for loops in the range I thru KMAX

1400	IF((KMAX-I).LT.5) THEN	! Region not long enough for any loop
		I=KMAX+1
		GOTO 1000
	ENDIF
	DO II=I,KMAX-5
		DO L=1,16
			LP(L)=.FALSE.
		ENDDO
		DMAX=0
		JE=II+15
		IF(JE.GT.KMAX) JE=KMAX
		DO J=II+1,JE
			DO M=II,J-1
				IF(DMAX.LT.DIST(M,J)) DMAX=DIST(M,J)
			ENDDO
			IF(J.GE.(II+5)) THEN
				IF(DIST(II,J).LT.CHKLP1.AND.DIST(II,J)
	1			.LT.(DMAX*CHKLP2)) LP(J-II+1)=.TRUE.
			ENDIF
		ENDDO

C****** Array LP now identifies all sizes of LOOPS permitted in this 
C****** segment of maximum length 16 that also match the requirements
C****** of the loop definition.
C****** Load DIST2 with all allowed loops each as a line of symbols, LOOP,
C****** of appropriate length on line I and then down to just above the 
C****** above the diagonal, to give an open triangle if nothing else 
C****** is in the way. Do not overwrite any other symbols.

		LCNT=0
		DO L=6,16
			IF(LP(L).EQ..TRUE.) THEN
				DO J=II+1,II+L-1
					IF(DIST2(II,J).EQ.' ') THEN
						DIST2(II,J)=LOOP
					ENDIF
				ENDDO
				JL=II+L-1
				DO M=II,J-1
					IF(DIST2(M,JL).EQ.' ') THEN
						DIST2(M,JL)=LOOP
					ENDIF
				ENDDO
				WRITE(2,5500) II,JL,L
5500				FORMAT(' SECND : ',I4,4X,I4,6X,I2)
				LCNT=LCNT+1
				WRITE(4,5520) II,JL,L
5520				FORMAT('LOOPS ',3I5)
			ENDIF
		ENDDO
		IF(LCNT.GT.0) THEN
			WRITE(2,5510)
5510			FORMAT(/)
		ENDIF
	ENDDO

	I=KMAX
	GOTO 1000

	END

C=========================     LIST2ND     =============================

C******  At this point the array DIST2 has been filled with characters 
C******  which correspond to the original DIST array modified by the 
C******  identification of alpha helix, sharp turns, five residue turns 
C******  and beta strands.
C******  These are contiguous stretches in the primary structure,
C******  and most obviously fit the definition of secondary structure.
C******  (In this program the beta sheet is considered part of the next
C******  level of structure.)

C******  Write to output file the beginning and ending alpha carbon
C******  numbers for each type of secondary structure. 
C******  A residue can be part of two consecutive secondary structure 
C******  units in sequence. In this case it is listed both as the last
C******  residue in the leading unit of secondary structure and as the 
C******  first residue in the following unit. ( The summed length of 
C******  such units, without correction, would then be longer than the
C******  actual peptide chain segment.

	SUBROUTINE LIST2ND (DIST,DIST2,MAXLEN,IPEND,JMAX,
	1		IFIRST,ILAST,SYMBOL,HELIX,H310,TRN,TRN5,STRAND)

	DIMENSION DIST(MAXLEN,MAXLEN)
	DIMENSION IFIRST(100),ILAST(100),LENGTH(100)
	
	INTEGER   HCNT,SCNT,TCNT,DCNT,UCNT,IZZ
	
	BYTE        DIST2(MAXLEN,MAXLEN),SYMBOL(100),ALAST,TEST
	BYTE	    HELIX,H310,TRN,TRN5,STRAND

C********************

	HCNT=0
	SCNT=0
	TCNT=0
	DCNT=0
	UCNT=0
	IZZ=0

	NCIS=0
	
	DO I=1,JMAX-1			! Check for cis peptide bonds
		IF(DIST2(I,I).EQ.'C') THEN
			WRITE(2,5030) I,I+1
5030			FORMAT(/,' PRIMRY: CIS PEPTIDE BOND BETWEEN ATOMS',
	1			I5,' -',I5)
			NCIS=NCIS+1
		ENDIF
	ENDDO
	IF(NCIS.EQ.0) THEN
		WRITE(2,5031)
5031		FORMAT(/,' PRIMRY: NO CIS PEPTIDE BONDS')
		WRITE(4,5032) IZZ,IZZ
5032		FORMAT('CISPEP',2I5)
	ENDIF
	WRITE(2,5001)
5001	FORMAT(/,' ---------------------------------------------------')
C********************
	WRITE(2,4999)
4999	FORMAT(1H1)
	WRITE(2,5001)

	DO I=1,100
		SYMBOL(I)=' '
	ENDDO

	N=1
	I=1
	J1CNT=0
	J2CNT=0

1000	IF(I.EQ.1.AND.J1CNT.EQ.0.AND.J2CNT.EQ.0) GOTO 1100

C******  Terminate segment at IEND

	ILAST(N)=IEND
	SYMBOL(N)=TEST
	IF(SYMBOL(N).EQ.H310) IFIRST(N)=IFIRST(N)+1
	LENGTH(N)=ILAST(N)-IFIRST(N)+1
	N=N+1			! increment counter for symbol list
	I=IEND			! reset start point
	IF(I.GT.JMAX-1) GOTO 2000	! check for completion

C******  Start search for next segment

1100	IF(DIST2(I,I+1).EQ.' '.AND.I.NE.1) I=I+1
	IFIRST(N)=I
	TEST=DIST2(I,I+1)
	J1CNT=0
	J2CNT=0
	IF(TEST.EQ.' ') GOTO 1400
	DO J=I+1,JMAX		! count # of TEST symbols on line I
		IF(DIST2(I,J).EQ.TEST) J1CNT=J1CNT+1
		IF(DIST2(I,J).EQ.' ')  GOTO 1200
	ENDDO

1200	IEND=I+J1CNT
	IEND1=IEND-1

1250	DO J=IEND,JMAX		! count # of TEST symbols on line IEND1
		IF(DIST2(IEND1,J).EQ.TEST) J2CNT=J2CNT+1
		IF(DIST2(IEND1,J).EQ.' ') GOTO 1300
	ENDDO

1300	IF(J2CNT.LE.1) GOTO 1000	! segment terminated OK at IEND

	IF(J2CNT.GT.1) THEN		! reset tentative segment end at
		IEND=IEND1+J2CNT	! IEND and try again
		IEND1=IEND-1
		J2CNT=0
		GOTO 1250
	ENDIF

C******  Special case of TEST = ' '.

1400	DO II=I,JMAX
		IF(DIST2(II,II+1).EQ.' ') THEN
			J1CNT=J1CNT+1
		ELSE		! Terminate blank segment
			ILAST(N)=IFIRST(N) + J1CNT - 1
			SYMBOL(N)='?'
			LENGTH(N)=J1CNT
			N=N+1
			I=I+J1CNT	! Reset start to next line
			IF(I.GT.JMAX-1) GOTO 2000  ! Search over terminate
			GOTO 1100
		ENDIF
	ENDDO

	ILAST(N)=JMAX
	LENGTH(N)=ILAST(N)-IFIRST(N)+1
	SYMBOL(N)=TEST
	IF(SYMBOL(N).EQ.' ') SYMBOL(N)='?'

2000	WRITE(2,5040)
5040	FORMAT(/,' SECND :  FIRST PASS - STRAIGHT SEGMENTS ONLY',/,
	1	 ' SECND :   SEQUENCE  SEQUENCE SEGMENT SECONDARY',/,
	1        ' SECND :      START       END  LENGTH STRUCTURE',/)
	DO I=1,100
		IF(IFIRST(I).EQ.0) GOTO 3000
		IF(SYMBOL(I).EQ.HELIX.OR.
	1	   SYMBOL(I).EQ.H310) THEN
			WRITE(2,5050) IFIRST(I),ILAST(I),LENGTH(I),SYMBOL(I)
5050			FORMAT(' SECND : ',2I10,I8,2X,A1)
		ELSEIF(SYMBOL(I).EQ.STRAND) THEN
			WRITE(2,5060) IFIRST(I),ILAST(I),LENGTH(I),SYMBOL(I)
5060			FORMAT(' SECND : ',2I10,I8,4X,A1)
		ELSEIF(SYMBOL(I).EQ.TRN.OR.
	1	       SYMBOL(I).EQ.TRN5) THEN
			WRITE(2,5070) IFIRST(I),ILAST(I),LENGTH(I),SYMBOL(I)
5070			FORMAT(' SECND : ',2I10,I8,6X,A1)
		ELSE
			WRITE(2,5080) IFIRST(I),ILAST(I),LENGTH(I),SYMBOL(I)
5080			FORMAT(' SECND : ',2I10,I8,8X,A1)
		ENDIF

	ENDDO

3000	DO I=1,100
		IF(IFIRST(I).EQ.0) GOTO 3100

		IF(SYMBOL(I).EQ.HELIX.OR.
	1	   SYMBOL(I).EQ.H310) THEN
			HCNT=HCNT+ILAST(I)-IFIRST(I)+1
		ELSEIF(SYMBOL(I).EQ.STRAND) THEN
			SCNT=SCNT+ILAST(I)-IFIRST(I)+1
		ELSEIF(SYMBOL(I).EQ.TRN.OR.
	1	       SYMBOL(I).EQ.TRN5) THEN
			TCNT=TCNT+ILAST(I)-IFIRST(I)+1
		ELSEIF(SYMBOL(I).EQ.'?') THEN
			UCNT=UCNT+ILAST(I)-IFIRST(I)+1
		ENDIF
		
		IF(I.GT.1) THEN
			IF(ILAST(I-1).EQ.IFIRST(I)) DCNT=DCNT-1
		ENDIF
	ENDDO

3100	WRITE(2,5100) HCNT,SCNT,TCNT,UCNT,DCNT
5100	FORMAT(//,' SECND : TOTAL NUMBER OF RESIDUES IN:',/,
	1	  ' SECND :               HELICES      :',I5,/,
	2	  ' SECND :               STRANDS      :',I5,/,
	3	  ' SECND :               TURNS        :',I5,/,
	4	  ' SECND :               UNCLASSIFIED :',I5,/,
	5	  ' SECND :               DOUBLE COUNTS:',I5)

	WRITE(2,5001)

	RETURN
	
	END	

C=======================================================================

C==================   SUBROUTINE LIBRARY   ===========================
C=======================   AXISCHK   =================================
C=================== (updated on 8 Mar 88) ===========================
C=====================================================================

	SUBROUTINE AXISCHK(DIST,DIST2,MAXLEN,JMAX,IPEND,
	1		IFIRST,ILAST,SYMBOL,CURV,ABC,INT,
	2		CHKHHD,CHKHSD,CHKSHD,XA,YA,ZA,
	3		HELIX,H310,STRAND,TRN,TRN5,
	4		HLXCOS,HLXPT,STRCOS,STRPT,IHSTRT,IHEND,ISSTRT,
	5		ISEND,HRAD,HDIST,SRAD,SDIST,ICP,NH,NS,
	6		HTHA,HTHMAX,STHA,STHMAX)

C****** 	FUNCTIONS AND SUBROUTINES CALLED:
C******		HAXIS   	SAXIS		VDOTV
C******		  DET3R		  
C******		  VXV		  
C******		  POINT

	DIMENSION DIST(MAXLEN,MAXLEN)
	DIMENSION IFIRST(100),ILAST(100),LENGTH(100),CURV(100)
	DIMENSION  XA(MAXLEN),YA(MAXLEN),ZA(MAXLEN)
	DIMENSION HLXCOS(ICP,3),HLXPT(ICP,3),STRCOS(ICP,3),STRPT(ICP,3)
	DIMENSION IHSTRT(ICP),IHEND(ICP),ISSTRT(ICP),ISEND(ICP)
	DIMENSION HRAD(ICP),HDIST(ICP),SRAD(ICP),SDIST(ICP)
	DIMENSION HP(3,50),ISEGTH(50,50),V1(3),V2(3)
	REAL*4    REJECT3(20),KEEP3(20), APT(3)

	DIMENSION HTHA(ICP),HTHMAX(ICP),STHA(ICP),STHMAX(ICP)

	INTEGER  HCNT,SCNT,TCNT,UCNT,DCNT,REJECT2(4,20),KEEP2(4,20)

	BYTE      DIST2(MAXLEN,MAXLEN),ABC(26),INT(12)
	BYTE	  TEST,HELIX,H310,STRAND,TRN,TRN5,SYMBOL(100)
	BYTE	  REJECT1(20),KEEP1(20)

	CHARACTER*3  FMID(30) 
	CHARACTER*4  FSTRT
	CHARACTER*6  FEND
	CHARACTER*13 FMT

C*******************************************************************

	DATA FMID/' 1X',' 3X',' 6X',' 9X','12X','15X','18X','21X',
	1	  '24X','27X','30X','33X','36X','39X','42X','45X',
	2	  '48X','51X','54X','57X','60X','63X','66X','69X',
	3	  '72X','75X','78X','81X','84X','87X'/

	DATA FSTRT/'(I3,'/
	DATA FEND /',25I3)'/

	FMT(1:4)=FSTRT
	FMT(8:13)=FEND

C*******************************************************************
	PI=ACOS(-1.0)
	RTOD=180./PI
C*******************************************************************
C******  Read secondary structure list. Pick out helices and strands
C******  and determine direction cosines of axes and N-terminal axis
C******  points.

	ISYMB=1
	NH=0		! NH=MAXH = # of helical segments in list
	NS=0		! NS=MAXS = # of strand  segments in list
	IR=0		! IR= number of pairs rejected for combination
	IA=0		! IA= number of pairs accepted for combination

100	DO N=ISYMB,100
		IF(IFIRST(N).EQ.0) GOTO	1000	! finished defining axes
		IF(SYMBOL(N).EQ.HELIX)  GOTO 200
		IF(SYMBOL(N).EQ.STRAND) GOTO 400	
	ENDDO

C-------------------- Start Helix Check ----------------------------

200	IS=IFIRST(N)
	IE=ILAST(N)
	ISYMB=N+1
	NN=N

C******  Helix must be at least 5 residues long.

	IEE=IE
	IF(SYMBOL(N+1).EQ.H310) THEN
		IEE=ILAST(N+1)
		ISYMB=N+2
		NN=N+1
	ENDIF
		
	IF((IEE-IS).LT.4) GOTO 100	!Pick next seg. in SS list

C******  Determine direction cosines of helix axis
C******  (CX,CY,CZ)-->HLXCOS(,)
C******   Determine N-terminal point on helix axis
C******   HP(I,1)-->HLXPT(,)

	CALL HAXIS(XA,YA,ZA,MAXLEN,IS,IEE,CX,CY,CZ,THA,THMAX,HP,D,R)

C******  Load arrays with direction cosines and coordinates of 
C******  points on helix axes.

	NH=NH+1
	HLXCOS(NH,1)=CX
	HLXCOS(NH,2)=CY
	HLXCOS(NH,3)=CZ
	HLXPT (NH,1)=HP(1,1)	! Only the first position in HP is
	HLXPT (NH,2)=HP(2,1)	! used. The other helix axis points
	HLXPT (NH,3)=HP(3,1)	! in HP are available for use before
	IHSTRT(NH)  =IS		! the next call to HAXIS or SAXIS
	IHEND(NH)   =IEE
	HRAD(NH)    =R
	HDIST(NH)   =D
  	HTHA(NH)    =THA	! replace with later value, see below
	HTHMAX(NH)  =THMAX	! replace with later value, see below

	L=IEE-IS
	DO K=1,3
		APT(K)=HP(K,1)+L*D*HLXCOS(NH,K)
	ENDDO
	WRITE(4,6000) SYMBOL(N),(HP(I,1),I=1,3),(APT(I),I=1,3),IS,IEE	
6000	FORMAT('ELEMNT',1X,A1,6F10.3,2I5)

C******   Calculate helix segment-angle matrix

 	NEXT=1
	IE3=IEE
	LEN1=IE3-IS+1

250	CONTINUE
	JE=IE3-IS-2
	ME=JE+1
		DO J=2,JE
			DO K=1,3
				V1(K)=HP(K,J+1)-HP(K,J)
			ENDDO
			DO M=J+1,ME
				DO K=1,3
					V2(K)=HP(K,M+1)-HP(K,M)
				ENDDO
				CALL VDOTV(V1,V2,RTOD,AB,TH)
				ISEGTH(J,M)=TH
			ENDDO	
		ENDDO
		TOT=0
		MCNT=0
		ISEGMX=0
		DO J=2,JE
			DO M=J+1,ME
				IF(ISEGMX.LT.ISEGTH(J,M)) ISEGMX=ISEGTH(J,M)
				TOT=TOT+ISEGTH(J,M)**2
				MCNT=MCNT+1
			ENDDO
		ENDDO
		IF(TOT.NE.0) THEN
			TOT=SQRT(TOT/MCNT)
		ENDIF

D		WRITE(2,5000) IS,IE3,(I,I=3,30)
D5000		FORMAT(/,' HELIX  ',I3,' - ',I3,/,6X,28I3,/)
D		IF(JE.LT.30) THEN
D			DO J=2,JE
D				FMT(5:7)=FMID(J)
D				WRITE(2,FMT) J,(ISEGTH(J,NJ),NJ=J+1,ME)
D			ENDDO
D		ENDIF
D		WRITE(2,5500) TOT
D5500		FORMAT(' RMS VALUE OF MATRIX ELEMENTS =  ',F5.0,/)
		IF(NEXT.EQ.1) CURVTMP1=TOT
		IF(NEXT.EQ.0) CURVTMP2=TOT
	IF(NEXT.EQ.0) THEN
	
		IF(TOT.LE.25.0) THEN
C******  If the rms angle deviation of the combined helices is less than
C******  25 degrees, then the helices should be combined and considered 
C******  curved. If helix 1 ends in a 310 segment, this segment is 
C******  considered to be alpha in the combined helix.
C******  Reset ILAST(N) to the end of the second alpha helix in the 
C******  new joint helix.
		
			IF(NN.EQ.N) THEN 
				ILAST(N)=ILAST(N+1)
				NSHFT=1
			ELSEIF(NN.GT.N) THEN
				ILAST(N)=ILAST(N+2)
				NSHFT=2
			ENDIF
			DO NC=N+1,98	! Contract stack beyond N 
				IFIRST(NC)=IFIRST(NC+NSHFT)
				ILAST (NC)=ILAST (NC+NSHFT)
				SYMBOL(NC)=SYMBOL(NC+NSHFT)
			ENDDO
D			WRITE(2,7000) IS,IE3,LEN1,LEN2
			GOTO 270
		ELSE
D			WRITE(2,7002)IS,IE,LEN1,LEN2,JS,ME
			IR=IR+1		! load rejected pairs
			REJECT1(IR)  =HELIX
			REJECT2(1,IR)=IS
			REJECT2(2,IR)=IEE
			REJECT2(3,IR)=IEE
			REJECT2(4,IR)=IE3
			REJECT3(IR)  =TOT
		ENDIF
		CURV(NN)=CURVTMP1
		GOTO 100

270		HTHA(NH)   =TOT
		HTHMAX(NH) =ISEGMX
		HLXCOS(NH,1)=CX
		HLXCOS(NH,2)=CY
		HLXCOS(NH,3)=CZ
		HLXPT (NH,1)=HP(1,1)	! Only the first position in HP is
		HLXPT (NH,2)=HP(2,1)	! used. The other helix axis points
		HLXPT (NH,3)=HP(3,1)	! in HP are available before next call
		IHSTRT(NH)  =IS		! start of combined helix
		IHEND(NH)   =IE3	! end of combined helix
		HRAD(NH)    =R		! mean radius of combined helix
		HDIST(NH)   =D		! mean repeat of combined helix

		L=IE3-IS
		DO K=1,3
			APT(K)=HP(K,1)+L*D*HLXCOS(NH,K)
		ENDDO
		BACKSPACE 4
		WRITE(4,6000) SYMBOL(N),(HP(I,1),I=1,3),(APT(I),I=1,3),IS,IE3	

		IA=IA+1			! load accepted pair into KEEP*
		KEEP1(IA)  =HELIX
		KEEP2(1,IA)=IS
		KEEP2(2,IA)=IEE
		KEEP2(3,IA)=IEE
		KEEP2(4,IA)=IE3
		KEEP3(IA)  =TOT
C******   Revise DIST2. Overwrite the original two helices with the new 
C******   combined helices, maintaining a final 310 helix if present.
		DO I=IS,IE2
			DO J=I+1,IE2
				DIST2(I,J)=HELIX
			ENDDO
		ENDDO
		IF(IE3.GT.IE2) THEN
			DO I=IS,IE2
				DO J=IE2+1,IE3
					DIST2(I,J)=H310
				ENDDO
			ENDDO 
			DO I=IE2+1,IE3
				DO J=I+1,IE3
					DIST2(I,J)=H310
				ENDDO
			ENDDO
		ENDIF
		ISYMB=N+1
		CURV(N)=CURVTMP2
		IF(IE3.GT.IE2) THEN
			ISYMB=N+2
			CURV(N)=0
			CURV(N+1)=CURVTMP2
		ENDIF
		GOTO100
	ENDIF

C******  Check following segment to see if it may be an extension of the 
C******  present helix as a single curved segment. The last atom of 
C******  present segment and first atom of next segment must be the same
C******  as well as the value of SYMBOL.
	
300	HTHA(NH)=TOT
	HTHMAX(NH)=ISEGMX
	IF(SYMBOL(NN+1).EQ.HELIX.AND.IFIRST(NN+1).EQ.ILAST(NN)) THEN
		IE2=ILAST(NN+1)
		IE3=IE2
		IF(SYMBOL(NN+2).EQ.H310) THEN
			IE3=ILAST(NN+2)
		ENDIF
		LEN2=IE3-IFIRST(NN+1)+1	
		CALL HAXIS(XA,YA,ZA,MAXLEN,IS,IE3,CX,CY,CZ,THA,THMAX,HP,D,R)
		NEXT=0
		GOTO 250
	ENDIF
	CURV(NN)=CURVTMP1
	GOTO 100
C--------- End of helix check section -------------------------------------
C--------- Start strand check section -------------------------------------

400	IS=IFIRST(N)
	IE=ILAST(N)
	ISYMB=N+1

C******  Determine direction cosines of strand axis
C******  (CX,CY,CZ)-->STRCOS(,)
C******  Determine N-terminal point on strand axis.
C******  HP(I,1)-->STRPT(,)

	IF((IE-IS).LT.3) GOTO 100

	CALL SAXIS(XA,YA,ZA,MAXLEN,IS,IE,CX,CY,CZ,THA,THMAX,HP,D,R)

	NS=NS+1
	STRCOS(NS,1)=CX
	STRCOS(NS,2)=CY
	STRCOS(NS,3)=CZ
	STRPT (NS,1)=HP(1,1)
	STRPT (NS,2)=HP(2,1)
	STRPT (NS,3)=HP(3,1)
	ISSTRT(NS)  =IS
	ISEND (NS)  =IE
	SRAD(NS)    =R
	SDIST(NS)   =D
	STHA(NS)    =THA	! replace with later numbers, see below
	STHMAX(NS)  =THMAX	! replace with later numbers, see below

	L=IE-IS
	DO K=1,3
		APT(K)=HP(K,1)+L*D*STRCOS(NS,K)
	ENDDO
	WRITE(4,6000) SYMBOL(N),(HP(I,1),I=1,3),(APT(I),I=1,3),IS,IE	

C******  Calculate strand segment-angle matrix

	IE3=IE
	NEXT=1

450	CONTINUE
	IF(NEXT.EQ.1) THEN 
		LEN1=IE3-IS+1
		IF(LEN1.LT.5) THEN
D			WRITE(2,7100) IS,IE3
D7100			FORMAT(' SHORT STRAND:',2I5)
			GOTO 500
		ENDIF
	ENDIF
	JE=IE3-IS-2
	ME=JE+1
		DO J=2,JE	! calc. angles and store as matrix elements
			DO K=1,3
				V1(K)=HP(K,J+1)-HP(K,J)
			ENDDO
			DO M=J+1,ME
				DO K=1,3
					V2(K)=HP(K,M+1)-HP(K,M)
				ENDDO
				CALL VDOTV(V1,V2,RTOD,AB,TH)
				ISEGTH(J,M)=TH
			ENDDO	
		ENDDO
		TOT=0
		MCNT=0
		ISEGMX=0
		DO J=2,JE	! sum squares of matrix elements
			DO M=J+1,ME
				IF(ISEGMX.LT.ISEGTH(J,M)) ISEGMX=ISEGTH(J,M)
				TOT=TOT+ISEGTH(J,M)**2
				MCNT=MCNT+1
			ENDDO
		ENDDO
		IF(TOT.NE.0) THEN 	! detm rms value of elements
			TOT=SQRT(TOT/MCNT)
		ENDIF

D		WRITE(2,5020) IS,IE3,(I,I=3,25)
D5020		FORMAT(/,' STRAND  ',I3,' - ',I3,/,6X,23I3,/)
D		DO J=2,JE
D			FMT(5:7)=FMID(J)
D			WRITE(2,FMT) J,(ISEGTH(J,NJ),NJ=J+1,ME)
D		ENDDO
D		WRITE(2,5500) TOT
		CURV(N)=TOT

	IF(NEXT.EQ.0) THEN
		JS=2
		MME=ME

		IF(TOT.LE.25.0) THEN
C******  If the rms angle deviation of the combined strands is less than
C******  25 degrees, then the strands should be combined and considered 
C******  curved. No new unassigned residues are produced.
C******  Reset ILAST(N) to the end of the joint strand
			ILAST(N)=ILAST(N+1)
			DO NC=N+1,99		! Contract stack beyond N
				IFIRST(NC)=IFIRST(NC+1)
				ILAST (NC)=ILAST (NC+1)
				SYMBOL(NC)=SYMBOL(NC+1)
			ENDDO

D			WRITE(2,7000) IS,IE3,LEN1,LEN2
D7000			FORMAT(' ACCEPTED PAIR: IS,IE3,LEN1,LEN2:'4I5)

			GOTO 470

		ELSEIF(TOT.GT.25.0.AND.LEN1.GT.4.AND.LEN2.GT.4) THEN
C******  If the rms angle deviation of the combined strands is greater than 
C******  25 degrees and if the length of each strand is greater than 4, then
C******  the strands should not be combined.
D			WRITE(2,7002)IS,IE,LEN1,LEN2,JS,ME

		ELSEIF(TOT.GT.25.0.AND.LEN1.EQ.4.AND.LEN2.GT.4) THEN
C******  If the rms angle deviation of the combined strands is greater than
C******  25 degrees and one strand is of length 4, then the short strand 
C******  should be peeled back from the end to see if an acceptable combined
C******  strand can be produced which is longer than the original long
C******  strand.
			DO JS=3,4
				TOT=0
				MCNT=0
				DO J=JS,JE
			  		DO M=J+1,ME
						TOT=TOT+ISEGTH(J,M)**2
						MCNT=MCNT+1
					ENDDO
				ENDDO
				IF(TOT.NE.0) TOT=SQRT(TOT/MCNT)
				IF(TOT.LE.25.) THEN 
	IWS=IS+JS-2
	IWE=IE3-ME+MME
D	WRITE(2,7001)IWS,IWE,LEN1,LEN2,JS,ME
D7001	FORMAT(' ACCEPTED PAIR: IS,IE3,LEN1,LEN2,JS,ME:',6I5)
C******  accept pair - goto 460 and reset LIST2ND tabulation
					GOTO 460
				ENDIF
			ENDDO
D	WRITE(2,7002)IS,IE3,LEN1,LEN2,JS,ME
D7002	FORMAT(' REJECTED PAIR: IS,IE3,LEN1,LEN2,JS,ME:',6I5)

		ELSEIF(TOT.GT.25.0.AND.LEN1.GT.4.AND.LEN2.EQ.4) THEN		
			DO MME=ME-1,ME-2,-1
				TOT=0
				MCNT=0
				DO J=2,MME-1
			  		DO M=J+1,MME
						TOT=TOT+ISEGTH(J,M)**2
						MCNT=MCNT+1
					ENDDO
				ENDDO
				IF(TOT.NE.0) THEN
					TOT=SQRT(TOT/MCNT)
				ENDIF
				IF(TOT.LE.25.) THEN
D	WRITE(2,7001)IS,IE3,LEN1,LEN2,JS,ME
C******  accept pair - goto 460 and reset LIST2ND tabulation
					GOTO 460
				ENDIF
			ENDDO
D	WRITE(2,7002)IS,IE,LEN1,LEN2,JS,ME
		ELSE
D	WRITE(2,7002)IS,IE,LEN1,LEN2,JS,ME
		ENDIF
			IR=IR+1		! load rejected pairs
			REJECT1(IR)  =STRAND
			REJECT2(1,IR)=IS
			REJECT2(2,IR)=IE
			REJECT2(3,IR)=IE
			REJECT2(4,IR)=IE3
			REJECT3(IR)  =TOT
		GOTO 100

C******  Reset IFIRST and ILAST depending on values of JS and MME
460		IFIRST(N)=IFIRST(N)+JS-2
		ILAST(N) =ILAST(N+1)-ME+MME
C******  If SYMBOL(N-1)=?, increment ILAST(N-1) for each shift in JS.
C******  If SYMBOL(N-1)=/?, insert new line N for SYMBOL(N)=?.
C******  If ILAST(N-1).NE.IFIRST(N), insert new line N but only if JS
C******  was incremented by more than 1. (If delta JS=1 there are no
C******  additional '?' entries to be listed).
		IF(JS.GT.2) THEN
			IF(SYMBOL(N-1).EQ.'?') ILAST(N-1)=ILAST(N-1)+JS-2
			IF(ILAST(N-1).EQ.IFIRST(N)-1) THEN
				DO NC=N+1,99	! Contract stack
					IFIRST(NC)=IFIRST(NC+1)
					ILAST (NC)=ILAST (NC+1)
					SYMBOL(NC)=SYMBOL(NC+1)
				ENDDO
			ELSE	! shift N to N+1 and insert new line N
				IFIRST(N+1)=IFIRST(N)
				ILAST (N+1)=ILAST (N)
				SYMBOL(N+1)=SYMBOL(N)
				CURV(N+1)  =CURV(N)
				IFIRST(N)  =ILAST(N-1)+1
				ILAST (N)  =IFIRST(N)+JS-2
				SYMBOL(N)  ='?'
				CURV(N)    =0
			ENDIF
		ELSEIF(MME.LT.ME) THEN
			IF(SYMBOL(N+2).EQ.'?') IFIRST(N+2)=IFIRST(N+2)-ME+MME
			IF(IFIRST(N+2).EQ.ILAST(N)+1) THEN
				DO NC=N+1,99	! Contract stack
					IFIRST(NC)=IFIRST(NC+1)
					ILAST (NC)=ILAST (NC+1)
					SYMBOL(NC)=SYMBOL(NC+1)
				ENDDO
			ELSE	! insert new line at N+1. no contraction.
				IFIRST(N+1)=ILAST(N)+1
				ILAST(N+1) =IFIRST(N+2)-1
				SYMBOL(N+1)='?'
				CURV(N+1)  = 0
			ENDIF
		ENDIF

470		STHA(NS)   =TOT
		STHMAX(NS) =ISEGMX
		STRCOS(NS,1)=CX
		STRCOS(NS,2)=CY
		STRCOS(NS,3)=CZ
		STRPT (NS,1)=HP(1,1)
		STRPT (NS,2)=HP(2,1)
		STRPT (NS,3)=HP(3,1)
		ISSTRT(NS)  =IS
		ISEND (NS)  =IE
		SRAD(NS)    =R
		SDIST(NS)   =D

		L=IE3-IS
		DO K=1,3
			APT(K)=HP(K,1)+L*D*STRCOS(NS,K)
		ENDDO
		BACKSPACE 4
		WRITE(4,6000) SYMBOL(N),(HP(I,1),I=1,3),(APT(I),I=1,3),IS,IE3	

		IA=IA+1		! load accepted pair
		KEEP1(IA)  =STRAND
		KEEP2(1,IA)=IS
		KEEP2(2,IA)=IE
		KEEP2(3,IA)=IE
		KEEP2(4,IA)=IE3
		KEEP3(IA)  =TOT
C******   Revise DIST2. Delete original strands and replace with 
C******   combined strands, shortened if required.
		DO I=IS,IE3
			DO J=I+1,IE3
				DIST2(I,J)=' '
			ENDDO
		ENDDO
		I2S=IS+JS-2
		I2E=IE3-ME+MME
		DO I=I2S,I2E
			DO J=I+1,I2E
				DIST2(I,J)=STRAND
			ENDDO
		ENDDO 
		ISYMB=N
		NS=NS-1
		GOTO100
	ENDIF

C******  Check following segment to see if it may be an extension of the 
C******  present strand as a single curved segment. Check SAXIS for the
C******  joint strands and check for fit to acceptable curve.

500	STHA(NS)=TOT
	STHMAX(NS)=ISEGMX
	IF(SYMBOL(N+1).EQ.STRAND.AND.IFIRST(N+1).EQ.ILAST(N)) THEN
		LEN2=ILAST(N+1)-IFIRST(N+1)+1
		IE3=ILAST(N+1)
		CALL SAXIS(XA,YA,ZA,MAXLEN,IS,IE3,CX,CY,CZ,THA,THMAX,HP,D,R)
		NEXT=0
		GOTO 450
	ENDIF

	GOTO 100

C-------------------- End of strand check section

1000	CONTINUE	! Finished defining segment axes and N-terminal points
	WRITE(2,4999)
4999	FORMAT(1H1)
	WRITE(2,5009)

2000	 WRITE(2,5040)
5040	FORMAT(/,' SECND :  SECOND PASS - CHECKED FOR CURVED SEGMENTS',/,
	1        ' SECND :   SEQUENCE  SEQUENCE SEGMENT SECONDARY MEAN_ANGLE',/,
	2        ' SECND :      START       END  LENGTH STRUCTURE WITH_AXIS',/)
	DO I=1,100
		IF(IFIRST(I).EQ.0) GOTO 3000
		LENGTH(I)=ILAST(I)-IFIRST(I)+1
		IF(SYMBOL(I).EQ.HELIX.OR.
	1	   SYMBOL(I).EQ.H310) THEN
		WRITE(2,5050) IFIRST(I),ILAST(I),LENGTH(I),SYMBOL(I),CURV(I)
5050			FORMAT(' SECND : ',2I10,I8,2X,A1,8X,F6.0)
		ELSEIF(SYMBOL(I).EQ.STRAND) THEN
		WRITE(2,5060) IFIRST(I),ILAST(I),LENGTH(I),SYMBOL(I),CURV(I)
5060			FORMAT(' SECND : ',2I10,I8,4X,A1,6X,F6.0)
		ELSEIF(SYMBOL(I).EQ.TRN.OR.
	1	       SYMBOL(I).EQ.TRN5) THEN
		WRITE(2,5070) IFIRST(I),ILAST(I),LENGTH(I),SYMBOL(I),CURV(I)
5070			FORMAT(' SECND : ',2I10,I8,6X,A1,4X,F6.0)
		ELSE
		WRITE(2,5080) IFIRST(I),ILAST(I),LENGTH(I),SYMBOL(I),CURV(I)
5080			FORMAT(' SECND : ',2I10,I8,8X,A1,2X,F6.0)
		ENDIF

	ENDDO

3000	DO I=1,100
		IF(IFIRST(I).EQ.0) GOTO 3100
		WRITE(4,5090) IFIRST(I),ILAST(I),LENGTH(I),SYMBOL(I),CURV(I)
5090		FORMAT('SECND ',3I5,A1,F5.0)
		IF(SYMBOL(I).EQ.HELIX.OR.
	1	   SYMBOL(I).EQ.H310) THEN
			HCNT=HCNT+ILAST(I)-IFIRST(I)+1
		ELSEIF(SYMBOL(I).EQ.STRAND) THEN
			SCNT=SCNT+ILAST(I)-IFIRST(I)+1
		ELSEIF(SYMBOL(I).EQ.TRN.OR.
	1	       SYMBOL(I).EQ.TRN5) THEN
			TCNT=TCNT+ILAST(I)-IFIRST(I)+1
		ELSEIF(SYMBOL(I).EQ.'?') THEN
			UCNT=UCNT+ILAST(I)-IFIRST(I)+1
		ENDIF
		
		IF(I.GT.1) THEN
			IF(ILAST(I-1).EQ.IFIRST(I)) DCNT=DCNT-1
		ENDIF
	ENDDO

3100	WRITE(2,5100) HCNT,SCNT,TCNT,UCNT,DCNT
5100	FORMAT(//,' SECND : TOTAL NUMBER OF RESIDUES IN:',/,
	1	  ' SECND :               HELICES      :',I5,/,
	2	  ' SECND :               STRANDS      :',I5,/,
	3	  ' SECND :               TURNS        :',I5,/,
	4	  ' SECND :               UNCLASSIFIED :',I5,/,
	5	  ' SECND :               DOUBLE COUNTS:',I5)

	WRITE(2,5200)
5200	FORMAT(//,' SECND : PAIRS ACCEPTED, ALL OR PART,',
	1	 ' AS SINGLE SEGMENTS',/,
	1	 ' SECND :                              RMS VALUE OF',/,
	1	 ' SECND :  TYPE  SEGMENT 1  SEGMENT 2  MATRIX ELEMENTS',/)
	IF(IA.EQ.0) WRITE(2,*) 'SECND :  NONE'
	DO I=1,IA
		WRITE(2,5110) KEEP1(I),(KEEP2(J,I),J=1,4),KEEP3(I)
5110		FORMAT(' SECND :',5X,A1,2X,I3,' - ',I3,3X,I3,' - ',I3,5X,F6.0)
	ENDDO
	WRITE(2,5210)
5210	FORMAT(//,' SECND : PAIRS REJECTED AS SINGLE SEGMENTS',/,
	1	 ' SECND :                              RMS VALUE OF',/,
	1	 ' SECND :  TYPE  SEGMENT 1  SEGMENT 2  MATRIX ELEMENTS',/)
	IF(IR.EQ.0) WRITE(2,*) 'SECND :  NONE'
	DO I=1,IR
		WRITE(2,5110) REJECT1(I),(REJECT2(J,I),J=1,4),REJECT3(I)
	ENDDO
	WRITE(2,5002)
5002	FORMAT(  /,' SECND : HELIX DIMENSIONS',/,
	1	   ' SECND :	   			       MEAN  SEGMENT',
	2	   '__ANGLES',/,
	3          ' SECND : NO.   HELIX   LENGTH  RADIUS  REPEAT   RMS  ',
	4	   '   MAX',/)
5003	FORMAT(' SECND :',  I4, I4,' -',I3,F8.1,2F8.3,2F8.1)   
	DO I=1,NH
		ALEN=FLOAT(IHEND(I)-IHSTRT(I))*HDIST(I)
		WRITE(2,5003) I,IHSTRT(I),IHEND(I),ALEN,HRAD(I),
	1			HDIST(I),HTHA(I),HTHMAX(I)
		WRITE(4,5103) I,IHSTRT(I),IHEND(I),ALEN,HRAD(I),
	1			HDIST(I),HTHA(I),HTHMAX(I)
5103		FORMAT('DIMENH',3I5,F6.1,2F6.3,2F6.1)
	ENDDO
	WRITE(2,5004)
5004	FORMAT(//, ' SECND : STRAND DIMENSIONS',/,
	1	   ' SECND :				       MEAN  SEGMENT',
	2	   '__ANGLES',/,
	3          ' SECND : NO.  STRAND   LENGTH  RADIUS  REPEAT   RMS  ',
	4	   '   MAX',/)




	DO I=1,NS
		ALEN=FLOAT(ISEND(I)-ISSTRT(I))*SDIST(I)
		WRITE(2,5003) I,ISSTRT(I),ISEND(I),ALEN,SRAD(I),
	1			SDIST(I),STHA(I),STHMAX(I)
		WRITE(4,5104) I,ISSTRT(I),ISEND(I),ALEN,SRAD(I),
	1			SDIST(I),STHA(I),STHMAX(I)
5104		FORMAT('DIMENS',3I5,F6.1,2F6.3,2F6.1)
	ENDDO
	WRITE(2,5009)
5009	FORMAT(/,' ---------------------------------------------------')
	DO L=1,100
		IF(IFIRST(L).NE.0) THEN
			LENGTH(L)=ILAST(L)-IFIRST(L)+1  
		ENDIF
	ENDDO

	WRITE(2,5009)

	RETURN

	END

C********************************************************************

	 SUBROUTINE HAXIS(XA,YA,ZA,MAXLEN,IS,IE,AX,AY,AZ,THA,THMAX,HP,D,R)

C		Subroutine calculates the mean axis of a helix from
C		alpha carbon positions and the point on the axis opposite
C		the N terminal C alpha.

C		SUBROUTINES CALLED: DET3R,VXV,POINT

C		XA,YA,ZA = alpha carbon coordinate list of dimension MAXLEN
C		IS,IE    = starting and ending indices of the helix atoms
C		AX,AY,AZ = direction cosines of the mean helix axis
C		THA      = mean angle of the helix axis segments with the axis
C		THMAX    = max angle of the segments with the axis
C		HP       = point on axis opposite the N terminal C alpha 
C               D = average distance between helix axis points
C               R = average radius of helix


	DIMENSION XA(MAXLEN),YA(MAXLEN),ZA(MAXLEN),A(4,3)
	DIMENSION HCOS(3,50),HP(3,50),P(3)
	DIMENSION V1(3),V2(3),V3(3),V4(3),V5(3)
	DIMENSION PLANES(4,3),SUM(3)

	DATA A(4,1)/1.0/,A(4,2)/1.0/,A(4,3)/1.0/

	AX=0
	AY=0
	AZ=0
	N=1
	DO I=IS+1,IE-2
		N=N+1
		DO J=1,3
			A(1,J)=XA(I+J-2)	! load points i-1,i,i+1
			A(2,J)=YA(I+J-2)	! defining plane 1
			A(3,J)=ZA(I+J-2)
		ENDDO

		CALL DET3R(A,B1,B2,B3,B4)	! Determine plane 

		PLANES(1,1)= B1			! sign on B4 is correct
		PLANES(2,1)=-B2			! for plane with equation
		PLANES(3,1)= B3			! Ax+By+Cz+D=0
		PLANES(4,1)=-B4

		V1(1)=XA(I-1)+XA(I+1)-2.*XA(I)	! vector perpendicular
		V1(2)=YA(I-1)+YA(I+1)-2.*YA(I)	! to helix axis from i
		V1(3)=ZA(I-1)+ZA(I+1)-2.*ZA(I)

		V2(1)=XA(I)+XA(I+2)-2.*XA(I+1)	! vector perpendicular
		V2(2)=YA(I)+YA(I+2)-2.*YA(I+1)	! to helix axis from i+1
		V2(3)=ZA(I)+ZA(I+2)-2.*ZA(I+1)

		CALL VXV(V1,V2,V3)		! V3 = unit vector parallel
		VL= SQRT(V3(1)**2+V3(2)**2+V3(3)**2)
		DO J=1,3			! to helix axis as defined
			HCOS(J,N)=V3(J)/VL	! by V1 and V2 i.e.
		ENDDO				! by points i-1,i,i+1,i+2

		CALL VXV(V1,V3,V4)		! V4 norm. to plane of V1,V3
		CALL VXV(V2,V3,V5)		! V4 norm. to plane of V2,V3
		
		PLANES(1,2)= V4(1)
		PLANES(2,2)= V4(2)
		PLANES(3,2)= V4(3)

		PLANES(1,3)= V5(1)
		PLANES(2,3)= V5(2)
		PLANES(3,3)= V5(3)

		PLANES(4,2)=-XA(I)  *V4(1)-YA(I)  *V4(2)-ZA(I)  *V4(3)
		PLANES(4,3)=-XA(I+1)*V5(1)-YA(I+1)*V5(2)-ZA(I+1)*V5(3)

C		All entries in PLANES and DET3R are for planes in the form
C			Ax+By+Cz+D=0

		CALL POINT(PLANES,P,ICHK)	! P = point on helix axis 
		IF(ICHK.EQ.-1) THEN		! near CAi
			TYPE 6000,I
6000			FORMAT(' HAXIS: POINT; ICHK=-1 ; I=',I4)
			STOP
		ENDIF
		DO J=1,3
			HP(J,N) = P(J)
		ENDDO
	ENDDO

C******   Put in last calculable axis point at IE-1
C******   V3,V4 and V5 are already loaded for proper planes 2 and 3
C******   Put in new plane 1.

	N=N+1
	DO J=1,3
		JJ=IE-1+J-2
		A(1,J)=XA(JJ)
		A(2,J)=YA(JJ)
		A(3,J)=ZA(JJ)
	ENDDO
	CALL DET3R(A,B1,B2,B3,B4)
	PLANES(1,1)=B1
	PLANES(2,1)=-B2
	PLANES(3,1)=B3
	PLANES(4,1)=-B4
	CALL POINT(PLANES,P,ICHK)
	IF(ICHK.EQ.-1) THEN		
		TYPE 6001
6001		FORMAT(' HAXIS: POINT; ICHK=-1 ; I= IE-1')
		STOP
	ENDIF
	DO J=1,3
		HP(J,N) = P(J)
	ENDDO
C******   HP and HCOS loading complete. Calc averages.
	SL=0
	DO J=1,3
		SUM(J) = 0
	ENDDO
	DO J=1,3
		DO I=2,N-1
			SUM(J)=SUM(J)+HCOS(J,I)
		ENDDO
		SUM(J)=SUM(J)/(N-2)	
		SL=SL+SUM(J)**2
	ENDDO
	SL=SQRT(SL)
	AX=SUM(1)/SL		! AX,AY,AZ = directions cosines of
	AY=SUM(2)/SL		! mean helix axis
	AZ=SUM(3)/SL
	THMAX=0
	THA=0

	DO I=2,N-1
		TH=ACOS(HCOS(1,I)*AX+HCOS(2,I)*AY+HCOS(3,I)*AZ)
		IF(THMAX.LT.TH) THMAX=TH
		THA=THA+TH
	ENDDO			! THMAX= maximum value of the angle
C				! of any individual axes with the 
	RAD=180./ACOS(-1.)	! mean helix axis.
C				! THA = average value of these angles
	THA=THA*RAD/(N-2)
	THMAX=THMAX*RAD
	DSQ=0
  	DO I=2,N-1 	!*pjf*
		DO J=1,3
			DSQ=DSQ+(HP(J,I+1)-HP(J,I))**2
		ENDDO
	ENDDO
 	D=SQRT(DSQ/(N-2))	!*pjf*
	RSQ=0
	DO I=2,N
		RSQ=RSQ+ (XA(IS+I-1)-HP(1,I))**2+
	1		 (YA(IS+I-1)-HP(2,I))**2+
	2	         (ZA(IS+I-1)-HP(3,I))**2
	ENDDO
	R=SQRT(RSQ/(N-1))
C		D = average distance between helix axis points
C		R = average radius of helix

	HP(1,1)=HP(1,2)-D*AX
	HP(2,1)=HP(2,2)-D*AY
	HP(3,1)=HP(3,2)-D*AZ

	RETURN

	END

C====================================================================

	 SUBROUTINE SAXIS(XA,YA,ZA,MAXLEN,IS,IE,AX,AY,AZ,THA,THMAX,HP,D,R)

C		Subroutine calculates the mean axis of a strand from

	DIMENSION XA(MAXLEN),YA(MAXLEN),ZA(MAXLEN),A(4,3)
	DIMENSION AAX(50),AAY(50),AAZ(50),HP(3,50)

	DATA A(4,1)/1.0/,A(4,2)/1.0/,A(4,3)/1.0/

	RTOD=180./ACOS(-1.0)
	N=0	
	AX=0
	AY=0
	AZ=0

C******   In a puckered straight chain the vector i-->i+2 is // to 
C******   strand axis

	DO I=IS,IE-2
		N=N+1		! N=number of (i,i+2) atom pairs
		DX=XA(I+2)-XA(I)
		DY=YA(I+2)-YA(I)
		DZ=ZA(I+2)-ZA(I)
		DL=SQRT(DX**2+DY**2+DZ**2)
		AAX(N)=DX/DL	! AAX,AAY,AAZ = direction cosines
		AAY(N)=DY/DL	! of axis segments
		AAZ(N)=DZ/DL
		AX=AX+AAX(N)
		AY=AY+AAY(N)
		AZ=AZ+AAZ(N)
	ENDDO

	EAL=SQRT(AX**2+AY**2+AZ**2)
	AX=AX/EAL	! AX,AY,AZ = mean direction cosines of
	AY=AY/EAL	! strand axis
	AZ=AZ/EAL
	
	THMAX=0
	THA=0
	DO I=1,N
		TH=ACOS(AAX(I)*AX+AAY(I)*AY+AAZ(I)*AZ)
		IF(THMAX.LT.TH) THMAX=TH
		THA=THA+TH
	ENDDO			! THMAX= maximum value of the angle
	THA=(THA/N)*RTOD	! of any individual axes with the 
	THMAX=THMAX*RTOD	! mean strand axis.
C				! THA= average value of these angles

C******  Establish helix axis points

	N=1		! N=number of helix axis points
	DO I=IS+1,IE-1
		N=N+1
		HP(1,N)=(XA(I-1)+XA(I+1)-2.*XA(I))/4.+XA(I)
		HP(2,N)=(YA(I-1)+YA(I+1)-2.*YA(I))/4.+YA(I)
		HP(3,N)=(ZA(I-1)+ZA(I+1)-2.*ZA(I))/4.+ZA(I)
	ENDDO
C******  Determine rms "radius" of strand
	RSQ=0
	DO I=2,N
		RSQ=RSQ+ (XA(IS+I-1)-HP(1,I))**2+
	1		 (YA(IS+I-1)-HP(2,I))**2+
	2		 (ZA(IS+I-1)-HP(3,I))**2
	ENDDO
	R=SQRT(RSQ/(N-1))
C******  Determine rms distance between strand axis points
	DSQ=0
	DO I=2,N-1
		DO J=1,3
			DSQ=DSQ+(HP(J,I+1)-HP(J,I))**2
		ENDDO
	ENDDO
	D=SQRT(DSQ/(N-2))
C******  Set up N terminal strand axis point
	HP(1,1)=HP(1,2)-D*AX
	HP(2,1)=HP(2,2)-D*AY
	HP(3,1)=HP(3,2)-D*AZ

	RETURN

	END
C=====================================================================

      SUBROUTINE DET3R(A,B1,B2,B3,B4)
C
C          A IS A 4X3 ARRAY. B1 IS THE VALUE OF THE DETERMINANT
C          COMPOSED OF ROWS 2,3,4; B2 OF ROWS 1,3,4; B3 OF ROWS 1,2,4;
C          B4 OF ROWS 1,2,3.
C
C          FUNCTIONS OR SUBROUTINES CALLED:  NONE
C
      DIMENSION A(4,3)
C
      A12=A(1,1)*A(2,2)-A(2,1)*A(1,2)
      A13=A(1,1)*A(3,2)-A(3,1)*A(1,2)
      A14=A(1,1)*A(4,2)-A(4,1)*A(1,2)
      A23=A(2,1)*A(3,2)-A(3,1)*A(2,2)
      A24=A(2,1)*A(4,2)-A(4,1)*A(2,2)
      A34=A(3,1)*A(4,2)-A(4,1)*A(3,2)
      B1=A23*A(4,3)-A24*A(3,3)+A34*A(2,3)
      B2=A13*A(4,3)-A14*A(3,3)+A34*A(1,3)
      B3=A12*A(4,3)-A14*A(2,3)+A24*A(1,3)
      B4=A12*A(3,3)-A13*A(2,3)+A23*A(1,3)
      RETURN
      END

C======================================================================

	SUBROUTINE POINT(PLANES,P,ICHK)

C******   Cramer's rule for point of intersection of 3 planes, each 
C******	  with the equation Ax+By+Cz+D=0
C******			| A1 A2 A3 |
C******			| B1 B2 B3 |
C******			| C1 C2 C3 |
C******			| D1 D2 D3 |
C******   (Note that signs are altered if equation for plane is taken
C******   as Ax+By+Cz=D )	
	
	DIMENSION PLANES(4,3),P(3)

	CALL DET3R(PLANES,B1,B2,B3,B4)

	ICHK=0

	IF(B4.NE.0) GOTO 100
		ICHK=-1
		RETURN

100	P(1)=-B1/B4
	P(2)= B2/B4
	P(3)=-B3/B4	

	RETURN

	END

C====================================================================

	SUBROUTINE VXV(V1,V2,V3)

C	Cross product of two vectors: V3 = V1 x V2

	DIMENSION V1(3),V2(3),V3(3)

	V3(1)= (V1(2)*V2(3) - V2(2)*V1(3))
	V3(2)=-(V1(1)*V2(3) - V2(1)*V1(3))
	V3(3)= (V1(1)*V2(2) - V2(1)*V1(2))

	RETURN

	END
C=====================================================================

	SUBROUTINE VDOTV(V1,V2,RTOD,AB,TH)

C******  First Vector
C******  Second Vector
C******  AB = Product of magnitudes
C******  TH = angle between vectors in degrees

	DIMENSION V1(3),V2(3)

	AB=0
	TH=0
	VDV = V1(1)*V2(1) + V1(2)*V2(2) +V1(3)*V2(3)
	A2  = V1(1)**2    + V1(2)**2    +V1(3)**2
	B2  = V2(1)**2    + V2(2)**2    +V2(3)**2
	AB  = SQRT(A2*B2)
	IF(AB.EQ.0)  RETURN

	TH  = ACOS(VDV/AB)*RTOD

	RETURN

	END

C=================================================================

C=================   SUBROUTINE LIBRARY   ============================
C=======================   SUPERSS   =================================
C================= (Updated on 29 Dec 87) ============================
C=====================================================================
C		BSHEET		RELATE
C		LISTBSH		MAPBOX
C		SUPERSS		CODE
C=====================================================================

C=======================   BSHEET   ==================================

C******  Check for beta strand pairs either parallel or antiparallel
C******  Start with search for neighboring atom which is not HELIX or TRN
C******  and which is closer than CHKDST. Label it SHORT, then check for
C******  a strand in both directions for match, as specified by CHKSHT,
C******  to the mask BB. Label all contiguous matching positions as SHEET.

	SUBROUTINE BSHEET (DIST,DIST2,MAXLEN,JMAX,STRAND,SHORT,
	1			SHEET,BB,HELIX,H310,TRN,TRN5,CHKSHT,CHKSDT)

	DIMENSION DIST(MAXLEN,MAXLEN),BB(13)

	BYTE	DIST2(MAXLEN,MAXLEN),STRAND,SHORT,SHEET,HELIX,H310
	BYTE	TEST,TRN,TRN5
	
	DO 1500 I=1,JMAX
		IF(DIST2(I,I+1).NE.STRAND) GOTO 1500
		DO 1400 J=I+5,JMAX
			IF(DIST2(J-1,J).NE.STRAND)   GOTO 1400
			IF(DIST(I,J).LT.CHKSDT) THEN
				DIST2(I,J) = SHORT
C******	Check up to 12 positions to left/right of SHORT for match
C****** to BB mask i.e. positions in strand identified by sequence 
C****** positions J+-12.
				LE=J-I-5
				IF(LE.GT.12) LE=12
				DO 1300 L=1,LE
				  IF(DIST2(J-L,J-L+1).NE.STRAND) GOTO 1320
				  IF(DIST2(I,J-L).GT.' ') GOTO 1300
					DD=ABS(DIST(I,J-L)-BB(L+1))
C******* ( Mismatch allowed to extent of CHKSHT)
					IF(DD.LE.CHKSHT) THEN
						DIST2(I,J-L)= SHEET
						GOTO 1300
					ELSE
						GOTO 1320
					ENDIF
1300				ENDDO
1320				LE=J+12
				IF(LE.GT.JMAX) LE=JMAX
				DO 1350 L=1,LE
				  IF(DIST2(J+L-1,J+L).NE.STRAND) GOTO 1400
				  IF(DIST2(I,J+L).GT.' ' ) GOTO 1350
					DD=ABS(DIST(I,J+L)-BB(L+1))
					IF(DD.LE.CHKSHT) THEN
						DIST2(I,J+L)= SHEET
						GOTO 1350
					ELSE
						GOTO 1400
					ENDIF
1350				ENDDO
			ENDIF
1400		CONTINUE
1500	CONTINUE

C	Eliminate accidental short distances and single line entries

	DO I=1,JMAX
		DO J=I+5,JMAX
			IF(DIST2(I,J).EQ.SHORT) THEN
				N=0
				IF(I.EQ.1) GOTO 2000
				IF(DIST2(I-1,J-1).EQ.SHORT) N=N+1
				IF(DIST2(I-1,J  ).EQ.SHORT) N=N+1
				IF(J.EQ.JMAX) GOTO 2000
				IF(DIST2(I-1,J+1).EQ.SHORT) N=N+1
2000				IF(I.EQ.JMAX) GOTO 2100
				IF(DIST2(I+1,J-1).EQ.SHORT) N=N+1
				IF(DIST2(I+1,J  ).EQ.SHORT) N=N+1
				IF(J.EQ.JMAX) GOTO 2100
				IF(DIST2(I+1,J+1).EQ.SHORT) N=N+1
2100				IF(N.EQ.0) THEN
					DIST2(I,J)=' '
					DO L=J-1,J-12,-1
						IF(DIST2(I,L).EQ.SHEET) THEN
							DIST2(I,L)=' '
						ELSE
							GOTO 2200
						ENDIF
					ENDDO
2200					DO L=J+1,J+12
						IF(DIST2(I,L).EQ.SHEET) THEN
							DIST2(I,L)=' '
						ELSE
							GOTO 2300
						ENDIF
					ENDDO
				ENDIF
			ENDIF
2300		ENDDO
	ENDDO	

	RETURN

	END

C=========================   LISTBSH   ===============================

C******  This subroutine lists the nearest neighbor atoms in the 
C******  identified sheets.

	SUBROUTINE LISTBSH (DIST,DIST2,MAXLEN,IPEND,JMAX,
	1		IFIRST,ILAST,SYMBOL,PAIR,HELIX,H310,TRN,TRN5,
	2		STRAND,LOOP,SHORT,SHEET)

	DIMENSION DIST(MAXLEN,MAXLEN)
	DIMENSION IFIRST(100),ILAST(100)
	
	INTEGER	    PAIR(100,6)
	BYTE        DIST2(MAXLEN,MAXLEN),SYMBOL(100),ALAST,TEST
	BYTE	    HELIX,H310,TRN,TRN5,STRAND,LOOP,SHORT,SHEET

C******  Look for off-diagonal atoms flagged as SHORT and set
C******  up neighbor list PAIR.

	L = 1
	DO 3100 I=1,IPEND
		PAIR(L,1) = I
		N=2
		DO 3050 J=I+5,IPEND
			IF(DIST2(I,J).EQ.SHORT) THEN
				PAIR(L,N) = J
				N=N+1
			ENDIF
3050		CONTINUE
		IF(N.GT.2) L=L+1
3100	CONTINUE
	PAIR(L,1)=0

	WRITE(2,5109)
5109	FORMAT(1H1)
	WRITE(2,5110)
5110	FORMAT(/,' SSS   : ADJACENT PAIRS IN SHEETS',/,
	1        ' SSS   :   FIRST      ADJACENT',/,
	2        ' SSS   : RESIDUE      RESIDUES',/)
	DO I=1,100
		IF(PAIR(I,1).EQ.0) GOTO 860
		WRITE(2,5120) (PAIR(I,N),N=1,6)
5120		FORMAT(' SSS   :',6I8)
		WRITE(4,5125) (PAIR(I,N),N=1,6)
5125		FORMAT('PAIRS ',6I5)
		IF(PAIR(I+1,1).NE.PAIR(I,1)+1) THEN
			WRITE(2,5130)
5130			FORMAT(' ------------------------------',
	1	'--------------------------')
		ENDIF
	ENDDO
860	CONTINUE

3200	RETURN

	END

C=======================   SUPERSS   =================================

	SUBROUTINE SUPERSS(DIST,DIST2,MAXLEN,JMAX,IPEND,
	1		IFIRST,ILAST,SYMBOL,ABC,INT,
	2		CHKHHD,CHKHSD,CHKSHD,XA,YA,ZA,
	3		LINS,IX,IY,BIT,COLR,NCNT,HELIX,H310,STRAND,
	4		SEPMAX,SEPMIN,AXDMAX,AXDMIN,
	5		PT1MAX,PT1MIN,PT2MAX,PT2MIN,
	6		HHPRNT,HSPRNT,SSPRNT,
	7		HLXCOS,HLXPT,STRCOS,STRPT,IHSTRT,IHEND,ISSTRT,
	8		ISEND,HRAD,HDIST,SRAD,SDIST,ICP,MAXH,MAXS)

C****** 	FUNCTIONS AND SUBROUTINES CALLED:
C******		RELATE	  MAPBOX    CODE	

	DIMENSION DIST(MAXLEN,MAXLEN),IFIRST(100),ILAST(100)
	DIMENSION  XA(MAXLEN),YA(MAXLEN),ZA(MAXLEN)
	DIMENSION HLXCOS(ICP,3),HLXPT(ICP,3),STRCOS(ICP,3),STRPT(ICP,3)
	DIMENSION IHSTRT(ICP),IHEND(ICP),ISSTRT(ICP),ISEND(ICP)
	DIMENSION HP(3,50),HRAD(ICP),HDIST(ICP),SRAD(ICP),SDIST(ICP)

	INTEGER*4 IX(3,LINS),IY(3,LINS),BIT(3,LINS),NCNT(3)

	BYTE      DIST2(MAXLEN,MAXLEN),ABC(26),INT(12)
	BYTE	  TEST,HELIX,H310,STRAND,SYMBOL(100)
	BYTE	  HHPRNT,HSPRNT,SSPRNT

	CHARACTER*7 COLR(3,LINS),COLBOX,COLN1,COLC1,COLN2,COLC2

C*******************************************************************
	PI=ACOS(-1.0)
C*******************************************************************

C**************************  PART I  **************************

C******  Determine angle and distance relations between segment
C******  axes.
C******  Start with helices.
C***************************

5001	FORMAT()
	WRITE(2,5007)
5007	FORMAT(///)
	WRITE(2,5010) 
5010	FORMAT(' SSS   : HELIX-HELIX ORIENTATIONS')

	WRITE(2,5020)
5020	FORMAT( ' SSS   : ',55X,'FRACTIONAL DISTANCE',/,
	1       ' SSS   : ',55X,'FROM N TERM. POINT ',/,
	2	' SSS   :                   MIN     ATOMS AT',
	3	'   H1-H2      AXIAL  TO MIN SEPN. POINT',/, 
	4	' SSS   : HELIX 1  HELIX 2  DIST    MIN DIST',
	5	'   ANGLE      SEPN.  HELIX 1  HELIX 2',/)

C*************************

	DO I=1,MAXH-1
		IS=IHSTRT(I)
		IE=IHEND(I)
		LI=IE-IS
		NW=0	
		DO J=I+1,MAXH
			JS=IHSTRT(J)
			JE=IHEND(J)
			LJ=JE-JS
			CALL RELATE(DIST,MAXLEN,I,J,IS,IE,JS,JE,
	1			HLXCOS,HLXCOS,HLXPT,HLXPT,ICP,ICHK,
	2			CHKHHD,IWRITE,DMIN,IM,JM,ANGL,D6,D1,D2)

			IF(IWRITE.EQ.0) GOTO 1500

			NW=NW+1
			D1=D1/(LI*HDIST(I))
			D2=D2/(LJ*HDIST(J))
			IF(D1.GT. 9.9) D1= 9.9
			IF(D1.LT.-9.9) D1=-9.9
			IF(D2.GT. 9.9) D2= 9.9
			IF(D2.LT.-9.9) D2=-9.9
			WRITE(2,5030) IS,IE,JS,JE,DMIN,IM,JM,ANGL,D6,D1,D2
5030			FORMAT(' SSS   : ',I3,'-',I3,2X,I3,'-',I3,F6.1,4X,
	1			I3,2X,I3,3X,F5.0,3(3X,F6.1))
			WRITE(4,5031) IS,IE,JS,JE,DMIN,IM,JM,ANGL,D6,D1,D2
5031			FORMAT('HLXHLX',4I5,F6.1,2I5,F6.0,3F6.1)
C******   	Load data in  DIST2 for line printer output
			ISIZE=2
			IF(LI.LT.4.OR.LJ.LT.4) ISIZE=1
			CALL CODE(DIST2,MAXLEN,INT,IS,IE,JS,JE,
	1				ANGL,D6,D1,D2,ISIZE)
			IF(HHPRNT.EQ.'N') GOTO 1500
                        IF(HHPRNT.EQ.'n') GOTO 1500

C******		Load data for MXMAP
			IF(DMIN.LT.SEPMIN.OR.DMIN.GT.SEPMAX) GOTO 1500
			IF(D6  .LT.AXDMIN.OR.D6  .GT.AXDMAX) GOTO 1500
			IF(D1  .LT.PT1MIN.OR.D1  .GT.PT1MAX) GOTO 1500
			IF(D2  .LT.PT2MIN.OR.D2  .GT.PT2MAX) GOTO 1500
			COLBOX='GREEN  '
			COLN1 ='BLUE   '
			COLC1 ='RED    '
			COLN2 ='BLUE   '
			COLC2 ='RED    '			
			CALL MAPBOX(JMAX,IS,IE,JS,JE,ANGL,
	1			LINS,IX,IY,BIT,COLR,NCNT,
	2			COLBOX,COLN1,COLC1,COLN2,COLC2)
1500		ENDDO
		IF(NW.NE.0) WRITE(2,5001)
	ENDDO

C**************************  PART II  *************************************

C******   Set up strand-helix relations

	WRITE(2,5007)
	WRITE(2,5011) 
5011	FORMAT(' SSS   : STRAND-HELIX  ORIENTATIONS')

	WRITE(2,5021)
5021	FORMAT(' SSS   : ',35X,'                    FRACTIONAL DISTANCE',/,
	1      ' SSS   : ',35X,'                    FROM N TERM. POINT ',/,
	2	   ' SSS   :                   MIN     ATOMS AT',
	3	   '   S1-H2      AXIAL  TO MIN SEPN. POINT',/, 
	4	   ' SSS   : STRND 1  HELIX 2  DIST    MIN DIST',
	5	   '   ANGLE      SEPN.  STRND 1  HELIX 2',/)

C**************************

	DO I=1,MAXS
		IS=ISSTRT(I)
		IE=ISEND(I)
		LI=IE-IS
		NW=0
		DO J=1,MAXH
			JS=IHSTRT(J)
			JE=IHEND(J)
			LJ=JE-JS
			IF(JS.LT.IE) GOTO 1600
			CALL RELATE(DIST,MAXLEN,I,J,IS,IE,JS,JE,
	1			STRCOS,HLXCOS,STRPT,HLXPT,ICP,ICHK,
	2			CHKSHD,IWRITE,DMIN,IM,JM,ANGL,D6,D1,D2)

			IF(IWRITE.EQ.0) GOTO 1600
			NW=NW+1
			
			D1=D1/(LI*SDIST(I))
			D2=D2/(LJ*HDIST(J))
			IF(D1.GT. 9.9) D1= 9.9
			IF(D1.LT.-9.9) D1=-9.9
			IF(D2.GT. 9.9) D2= 9.9
			IF(D2.LT.-9.9) D2=-9.9
			
			WRITE(2,5030) IS,IE,JS,JE,DMIN,IM,JM,ANGL,D6,D1,D2
			WRITE(4,5032) IS,IE,JS,JE,DMIN,IM,JM,ANGL,D6,D1,D2
5032			FORMAT('STRHLX',4I5,F6.1,2I5,F6.0,3F6.1)
C******   	Load data in DIST2 for line printer output
			ISIZE=2
			IF(LI.LT.4.OR.LJ.LT.4) ISIZE=1
			CALL CODE(DIST2,MAXLEN,INT,IS,IE,JS,JE,
	1				ANGL,D6,D1,D2,ISIZE)
			IF(HSPRNT.EQ.'N') GOTO 1600
                        IF(HSPRNT.EQ.'n') GOTO 1600

C******		Load data for MXMAP
			IF(DMIN.LT.SEPMIN.OR.DMIN.GT.SEPMAX) GOTO 1600
			IF(D6  .LT.AXDMIN.OR.D6  .GT.AXDMAX) GOTO 1600
			IF(D1  .LT.PT1MIN.OR.D1  .GT.PT1MAX) GOTO 1600
			IF(D2  .LT.PT2MIN.OR.D2  .GT.PT2MAX) GOTO 1600
			COLBOX='WHITE  '
			COLN1 ='GREEN  '
			COLC1 ='WHITE  '
			COLN2 ='BLUE   '
			COLC2 ='RED    '			
			CALL MAPBOX(JMAX,IS,IE,JS,JE,ANGL,
	1			LINS,IX,IY,BIT,COLR,NCNT,
	2			COLBOX,COLN1,COLC1,COLN2,COLC2)
1600		ENDDO
		IF(NW.NE.0) WRITE(2,5001)
	ENDDO

C**************************  PART III  ********************************

C******   Set up helix-strand relations

	WRITE(2,5007)

	WRITE(2,5012) 
5012	FORMAT(' SSS   : HELIX-STRAND  ORIENTATIONS')

	WRITE(2,5022)
5022	FORMAT(' SSS   : ',35X,'                    FRACTIONAL DISTANCE',/,
	1      ' SSS   : ',35X,'                    FROM N TERM. POINT ',/,
	2	   ' SSS   :                   MIN     ATOMS AT',
	3	   '   H1-S2      AXIAL  TO MIN SEPN. POINT',/, 
	4	   ' SSS   : HELIX 1  STRND 2  DIST    MIN DIST',
	5	   '   ANGLE      SEPN.  HELIX 1  STRND 2',/)

C*****************************

	DO I=1,MAXH
		IS=IHSTRT(I)
		IE=IHEND(I)
		LI=IE-IS
		NW=0
		DO J=1,MAXS
			JS=ISSTRT(J)
			JE=ISEND(J)
			LJ=JE-JS
			IF(JS.LT.IE) GOTO 1700
			CALL RELATE(DIST,MAXLEN,I,J,IS,IE,JS,JE,
	1			HLXCOS,STRCOS,HLXPT,STRPT,ICP,ICHK,
	2			CHKHSD,IWRITE,DMIN,IM,JM,ANGL,D6,D1,D2)

			IF(IWRITE.EQ.0) GOTO 1700
			NW=NW+1

			D1=D1/(LI*HDIST(I))
			D2=D2/(LJ*SDIST(J))
			IF(D1.GT. 9.9) D1= 9.9
			IF(D1.LT.-9.9) D1=-9.9
			IF(D2.GT. 9.9) D2= 9.9
			IF(D2.LT.-9.9) D2=-9.9

			WRITE(2,5030) IS,IE,JS,JE,DMIN,IM,JM,ANGL,D6,D1,D2
			WRITE(4,5033) IS,IE,JS,JE,DMIN,IM,JM,ANGL,D6,D1,D2
5033			FORMAT('HLXSTR',4I5,F6.1,2I5,F6.0,3F6.1)
C******   	Load data in DIST2 for line printer output
			ISIZE=2
			IF(LI.LT.4.OR.LJ.LT.4) ISIZE=1
			CALL CODE(DIST2,MAXLEN,INT,IS,IE,JS,JE,
	1				ANGL,D6,D1,D2,ISIZE)
			IF(HSPRNT.EQ.'N') GOTO 1700
                        IF(HSPRNT.EQ.'n') GOTO 1700

C******		Load data for MXMAP
			IF(DMIN.LT.SEPMIN.OR.DMIN.GT.SEPMAX) GOTO 1700
			IF(D6  .LT.AXDMIN.OR.D6  .GT.AXDMAX) GOTO 1700
			IF(D1  .LT.PT1MIN.OR.D1  .GT.PT1MAX) GOTO 1700
			IF(D2  .LT.PT2MIN.OR.D2  .GT.PT2MAX) GOTO 1700
			COLBOX='WHITE  '
			COLN1 ='BLUE   '
			COLC1 ='RED    '
			COLN2 ='GREEN  '
			COLC2 ='WHITE  '			
			CALL MAPBOX(JMAX,IS,IE,JS,JE,ANGL,
	1			LINS,IX,IY,BIT,COLR,NCNT,
	2			COLBOX,COLN1,COLC1,COLN2,COLC2)
1700		ENDDO
		IF(NW.NE.0) WRITE(2,5001)
	ENDDO


C***************************  PART IV *****************************

C******   Check angles between adjacent strand segments 
C******    (evidence for long curves?)

	WRITE(2, 5007)

	WRITE(2,5013) 
5013	FORMAT(' SSS   : ANGLES FOR STRAND-STRAND PAIRS')

	WRITE(2,5023)
5023	FORMAT(' SSS   : ',35X,'                    FRACTIONAL DISTANCE',/,
	1      ' SSS   : ',35X,'                    FROM N TERM. POINT ',/,
	2	   ' SSS   :                   MIN     ATOMS AT',
	3	   '   A1-A2      AXIAL  TO MIN SEPN. POINT',/, 
	4	   ' SSS   : STRND 1  STRND 2  DIST    MIN DIST',
	5	   '   ANGLE      SEPN.  STRND 1  STRND 2',/)

C*****************************

	CHK= CHKHHD
	DO I=1,MAXS-1
		IS=ISSTRT(I)
		IE=ISEND(I)
		LI=IE-IS
		NW=0
		DO J=I+1,MAXS
			JS=ISSTRT(J)
			JE=ISEND (J)
			LJ=JE-JS
			
			CALL RELATE(DIST,MAXLEN,I,J,IS,IE,JS,JE,
	1			STRCOS,STRCOS,STRPT,STRPT,ICP,ICHK,
	2			CHK,IWRITE,DMIN,IM,JM,ANGL,D6,D1,D2)

			IF(IWRITE.EQ.0) GOTO 1800
			NW=NW+1

			D1=D1/(LI*SDIST(I))
			D2=D2/(LJ*SDIST(I+1))
			IF(D1.GT. 9.9) D1= 9.9
			IF(D1.LT.-9.9) D1=-9.9
			IF(D2.GT. 9.9) D2= 9.9
			IF(D2.LT.-9.9) D2=-9.9

			WRITE(2,5030) IS,IE,JS,JE,DMIN,IM,JM,ANGL,D6,D1,D2
			WRITE(4,5034) IS,IE,JS,JE,DMIN,IM,JM,ANGL,D6,D1,D2
5034			FORMAT('STRSTR',4I5,F6.1,2I5,F6.0,3F6.1)
			NC=0
			DO K=IS,IE-1
				DO L=JS+1,JE
					IF(DIST2(K,L).NE.' ') NC=NC+1
				ENDDO
			ENDDO
			IF(NC.NE.0) GOTO 1800
C******   	Load data in DIST2 for line printer output
			ISIZE=2
			IF(LI.LT.4.OR.LJ.LT.4) ISIZE=1
			CALL CODE(DIST2,MAXLEN,INT,IS,IE,JS,JE,
	1				ANGL,D6,D1,D2,ISIZE)
			IF(SSPRNT.EQ.'N') GOTO 1800
                        IF(SSPRNT.EQ.'n') GOTO 1800

C******		Load data for MXMAP
			IF(DMIN.LT.SEPMIN.OR.DMIN.GT.SEPMAX) GOTO 1800
			IF(D6  .LT.AXDMIN.OR.D6  .GT.AXDMAX) GOTO 1800
			IF(D1  .LT.PT1MIN.OR.D1  .GT.PT1MAX) GOTO 1800
			IF(D2  .LT.PT2MIN.OR.D2  .GT.PT2MAX) GOTO 1800
			DO K=IS,IE
				DO L=JS,JE
					IF(DIST2(K,J).NE.' ')
	1					GOTO 1800
				ENDDO
			ENDDO

			COLBOX='BLUE '
			COLN1 ='GREEN  '
			COLC1 ='WHITE  '
			COLN2 ='GREEN  '
			COLC2 ='WHITE  '			
			CALL MAPBOX(JMAX,IS,IE,JS,JE,ANGL,
	1			LINS,IX,IY,BIT,COLR,NCNT,
	2			COLBOX,COLN1,COLC1,COLN2,COLC2)
1800			ENDDO
		IF(NW.NE.0) WRITE(2,5001)
	ENDDO


	RETURN

	END

C========================   RELATE   ================================

	SUBROUTINE RELATE(DIST,MAXLEN,I,J,IS,IE,JS,JE,
	1		FCOS,GCOS,APT,BPT,ICP,ICHK,
	2		CHK,IWRITE,DMIN,IM,JM,ANGL,D6,D1,D2)

C******  DIST      = distance matrix of dimension MAXLEN**2
C******  I,J       = helix numbers in secondary structure list
C******  IS,IE     = start and end residue numbers for helix I
C******  JS,JE     = start and end residue numbers for helix J
C******  FCOS,GCOS = direction cosine arrays
C******  APT,BPT   = N terminal axis point arrays
C******  ICP       = Size for above 4 arrays
C******  CHK       = calculations omitted for min distance between
C******                segments greater than this value.

	DIMENSION  DIST(MAXLEN,MAXLEN)
	DIMENSION  FCOS(ICP,3),GCOS(ICP,3),APT(ICP,3),BPT(ICP,3)
	DIMENSION  V0(3),V1(3),V2(3),VP(3),P1(3),P2(3),PLANES(4,3)

C******	 Find minimum distance in rectangle IS,IE,JS,JE
	
			IWRITE=1
			DMIN=10000.
			DO K=IS,IE
				DO L=JS,JE
					IF(DMIN.GT.DIST(K,L)) THEN
						DMIN=DIST(K,L)
						IM=K
						JM=L
					ENDIF
				ENDDO
			ENDDO

			IF(DMIN.GT.CHK) THEN
				IWRITE=0
				RETURN
			ENDIF

C******  V0=vector product = FCOS(I,) x GCOS(J,)

			V0(1)= (FCOS(I,2)*GCOS(J,3)-
	1			FCOS(I,3)*GCOS(J,2))
			V0(2)= (FCOS(I,3)*GCOS(J,1)-
	1			FCOS(I,1)*GCOS(J,3))
			V0(3)= (FCOS(I,1)*GCOS(J,2)-
	1			FCOS(I,2)*GCOS(J,1))

C******  V1=vector product = V0 x FCOS(I,)

			V1(1)= V0(2)*FCOS(I,3)-V0(3)*FCOS(I,2)
			V1(2)= V0(3)*FCOS(I,1)-V0(1)*FCOS(I,3)
			V1(3)= V0(1)*FCOS(I,2)-V0(2)*FCOS(I,1)

C******  V2=vector product = V0 x GCOS(J,)

			V2(1)= V0(2)*GCOS(J,3)-V0(3)*GCOS(J,2)
			V2(2)= V0(3)*GCOS(J,1)-V0(1)*GCOS(J,3)
			V2(3)= V0(1)*GCOS(J,2)-V0(2)*GCOS(J,1)

			DO K=1,3
				PLANES(K,1)=V0(K)
				PLANES(K,2)=V1(K)
				PLANES(K,3)=V2(K)
			ENDDO
			PLANES(4,1)=-(V0(1)*APT(I,1)+
	1			    V0(2)*APT(I,2)+
	2			    V0(3)*APT(I,3))
			PLANES(4,2)=-(V1(1)*APT(I,1)+
	1			    V1(2)*APT(I,2)+
	2			    V1(3)*APT(I,3))
			PLANES(4,3)=-(V2(1)*BPT(J,1)+
	1			    V2(2)*BPT(J,2)+
	2			    V2(3)*BPT(J,3))

			D1=0
			D2=0
			D6=0
			ANGL=0
			SUM1=0		! P1 = intersection of min distance
			SUM2=0		!      line with FCOS(I)
			SUM3=0		! P2 = intersection of min distance
			SUM4=0		!      line with GCOS(J)	
			SUM5=0
			SUM6=0
			SUM7=0

			CALL POINT(PLANES,P1,ICHK)
			IF(ICHK.EQ.-1) THEN
				TYPE 6000,I
6000				FORMAT(' RELATE: POINT(1); ICHK=-1 ; I=',I4)
				RETURN
			ENDIF

			PLANES(4,1)=-(V0(1)*BPT(J,1)+
	1			    V0(2)*BPT(J,2)+
	2			    V0(3)*BPT(J,3))

			CALL POINT(PLANES,P2,ICHK)
			IF(ICHK.EQ.-1) THEN
				TYPE 6001,I
6001				FORMAT(' RELATE: POINT(2); ICHK=-1 ; I=',I4)
				RETURN
			ENDIF

			DO K=1,3
				VP(K)=P2(K)-P1(K)
			ENDDO


			DO K=1,3
				SUM1=SUM1+(P1(K)-APT(I,K))**2 ! SQD,APT-P1
				SUM2=SUM2+(P2(K)-BPT(J,K))**2 ! SQD,BPT-P2
				SUM3=SUM3+FCOS(I,K)*(P1(K)-APT(I,K)) ! dot prod
			 	SUM4=SUM4+GCOS(J,K)*(P2(K)-BPT(J,K)) ! dot prod
				SUM5=SUM5+FCOS(I,K)*GCOS(J,K) ! dot prod.
				SUM6=SUM6+VP(K)**2            ! SQD,P2-P1
				SUM7=SUM7+V0(K)*VP(K)
			ENDDO

			D1=SQRT(SUM1)		! Dist from APT to P1
			D2=SQRT(SUM2)		! Dist from BPT to P2
			D6=SQRT(SUM6)		! Dist from P1  to P2
			D1=D1*SUM3/ABS(SUM3)	! +=same dir. as helix axis
			D2=D2*SUM4/ABS(SUM4)	! -=opp. dir. to helix axis
			ANGL=ACOS(SUM5)*180./ACOS(-1.0)
			ANGL=ANGL*SUM7/ABS(SUM7)
C******  ANGL is -(0-180) for ccw rotation around VP looking in 
C******  positive direction. ANGL is +(0-180) for CW rotation to move
C******  V1 into V2.

D	TYPE *,I,J,P1(1),P1(2),P1(3)
D	TYPE *,'                    ',APT(I,1),APT(I,2),APT(I,3)
	RETURN

	END

C========================   MAPBOX   ======================================

	SUBROUTINE MAPBOX(JMAX,IS,IE,JS,JE,ANGL,
	1		LINS,IX,IY,BIT,COLR,NCNT,
	2		COLBOX,COLN1,COLC1,COLN2,COLC2)

C**************************************************************************

C	ANGL = Angle between secondary segment axes.
C	JMAX = maximum number of CA atoms in input list
c	IS,IE,JS,JE = indices identifying the corners of the box to be outlined
C		      The indices are multiplied by 10 to increase the 
C		      resolution of the output drawing by UGSYS.
C	IX,IY,BIT,COLR arrays containing lines to be passed to MXMAP
C	COLN1 etc are the widths of the lines to be drawn. The COLOR parameter
C	is used by G.Furnish' version of UGSYS to control line width.
C	The first index in IX,IY,BIT,COLR refers to the panel to be drawn
C	identified by NQ.
C			For NQ=1: I=   1,260	J=   1,260
C			For NQ=2: I=   1,260	J= 241,500
C			For NQ=3: I= 241,500	J= 241,500
C***************************************************************************
	DIMENSION  IX(3,LINS),IY(3,LINS),BIT(3,LINS),NCNT(3)

	CHARACTER*7  COLR(3,LINS),COLN1,COLC1,COLN2,COLC2,COLBOX	


	IS10=IS*10
	IE10=IE*10	
	JS10=JS*10
	JE10=JE*10
	II10=(IS10+IE10)/2	! identify center of box
	JJ10=(JS10+JE10)/2	! 
	I2=(IE10-IS10)/2	! half edge lengths of box
	J2=(JE10-JS10)/2	!
	PI=ACOS(-1.0)
	RANGL=ANGL*PI/180.

	IF(IS.GE.241.AND.JS.GE.241) NQ=3
	IF(IE.LE.260.AND.JS.GE.241) NQ=2
	IF(IE.LE.260.AND.JE.LE.260) NQ=1
	IF(IS.LT.241.AND.IE.GT.260) RETURN
	IF(JS.LT.241.AND.JE.GT.260) RETURN

C******   Outline box for MXMAP output.

	NCNT(NQ)=NCNT(NQ)+1	
	IX(NQ,NCNT(NQ))  =IS10
	IY(NQ,NCNT(NQ))  =JS10
	BIT(NQ,NCNT(NQ)) =0
	COLR(NQ,NCNT(NQ))='       '
		NCNT(NQ)=NCNT(NQ)+1
	IX(NQ,NCNT(NQ))  =IS10
	IY(NQ,NCNT(NQ))  =JE10
	BIT(NQ,NCNT(NQ)) =1
	COLR(NQ,NCNT(NQ))=COLBOX
		NCNT(NQ)=NCNT(NQ)+1
	IX(NQ,NCNT(NQ))  =IE10
	IY(NQ,NCNT(NQ))  =JE10
	BIT(NQ,NCNT(NQ)) =1
	COLR(NQ,NCNT(NQ))=COLBOX
		NCNT(NQ)=NCNT(NQ)+1
	IX(NQ,NCNT(NQ))  =IE10
	IY(NQ,NCNT(NQ))  =JS10
	BIT(NQ,NCNT(NQ)) =1
	COLR(NQ,NCNT(NQ))=COLBOX
		NCNT(NQ)=NCNT(NQ)+1
	IX(NQ,NCNT(NQ))  =IS10
	IY(NQ,NCNT(NQ))  =JS10
	BIT(NQ,NCNT(NQ)) =1
	COLR(NQ,NCNT(NQ))=COLBOX

C******   Put in first segment axis for MXMAP as vertical line 
C******   with N terminus at top

		NCNT(NQ)=NCNT(NQ)+1
	IX(NQ,NCNT(NQ))  =IS10
	IY(NQ,NCNT(NQ))  =JJ10
	BIT(NQ,NCNT(NQ)) =0
	COLR(NQ,NCNT(NQ))='       '
		NCNT(NQ)=NCNT(NQ)+1
	IX(NQ,NCNT(NQ))  =II10
	IY(NQ,NCNT(NQ))  =JJ10
	BIT(NQ,NCNT(NQ)) =1
	COLR(NQ,NCNT(NQ))=COLN1
		NCNT(NQ)=NCNT(NQ)+1
	IX(NQ,NCNT(NQ))  =IE10
	IY(NQ,NCNT(NQ))  =JJ10
	BIT(NQ,NCNT(NQ)) =1
	COLR(NQ,NCNT(NQ))=COLC1

C******   Determine terminal points on box edges for second segment
C******   axis.

	KN=I2*TAN(RANGL)
	KNA=ABS(KN)
	LN=J2*TAN(PI/2.-RANGL)
	LNA=ABS(LN)
	IF(ANGL.GE.0) THEN
		IF(KNA.LE.J2) THEN	! start on line IS10 or IE10
			IN10=IS10	! first half
			IF(KN.LT.0) IN10=IE10
			JN10=JJ10-KNA
			IC10=IE10
			IF(KN.LT.0) IC10=IS10
			JC10=JJ10+KNA
		ELSE			! start on line JS10
			IN10=II10-LN
			JN10=JS10
			IC10=II10+LN
			JC10=JE10
		ENDIF
	ELSE			! ANGL.LT.0
		IF(KNA.LE.J2) THEN	! start on line IS10 or IE10
			IN10=IS10	! second half
			IF(KN.GT.0) IN10=IE10
			JN10=JJ10+KNA
			IC10=IE10
			IF(KN.GT.0) IC10=IS10
			JC10=JJ10-KNA
		ELSE			! start on line JE10
			IN10=II10+LN
			JN10=JE10
			IC10=II10-LN
			JC10=JS10
		ENDIF
	ENDIF

C******  Load second segment

		NCNT(NQ) =NCNT(NQ)+1
	IX  (NQ,NCNT(NQ))=IN10
	IY  (NQ,NCNT(NQ))=JN10
	BIT (NQ,NCNT(NQ))=0
	COLR(NQ,NCNT(NQ))='       '
		NCNT(NQ) =NCNT(NQ)+1
	IX  (NQ,NCNT(NQ))=II10
	IY  (NQ,NCNT(NQ))=JJ10
	BIT (NQ,NCNT(NQ))=1
	COLR(NQ,NCNT(NQ))=COLN2
		NCNT(NQ) =NCNT(NQ)+1
	IX  (NQ,NCNT(NQ))=IC10
	IY  (NQ,NCNT(NQ))=JC10
	BIT (NQ,NCNT(NQ))=1
	COLR(NQ,NCNT(NQ))=COLC2

	RETURN

	END

C==========================   CODE   ===================================

	SUBROUTINE CODE(DIST2,MAXLEN,INT,IS,IE,JS,JE,ANGL,D6,
	1			D1,D2,ISIZE)

	BYTE DIST2(MAXLEN,MAXLEN)
	BYTE INT(12)

	IF(ISIZE.EQ.1) GOTO 100
	IF(ISIZE.EQ.2) GOTO 200

C******  Set up ANGL in units of 10 degrees: 0-->+-18

100	IF(ANGL.LT.0) DIST2(IS,JE-2)='-'
	LANGL=(ABS(ANGL)+0.5)/10.0
	L1ANGL=LANGL/10
	L2ANGL=LANGL-10*L1ANGL
	DIST2(IS,JE-1)=INT(L1ANGL+1)
	DIST2(IS,JE)  =INT(L2ANGL+1)

C******  D6 truncate to the nearest Angstrom unit.
C******  Set up D6 under ANGL
	
	LD6=(ABS(D6)+0.5)
	L1D6=LD6/100
	L2D6=(LD6-100*L1D6)/10
	L3D6=LD6-100*L1D6-10*L2D6
	DIST2(IS+1,JE-2) = INT(L1D6+1)
	DIST2(IS+1,JE-1) = INT(L2D6+1)
	DIST2(IS+1,JE  ) = INT(L3D6+1)

	RETURN

200	IF(ANGL.LT.0) DIST2(IS,JE-3)='-'
	LANGL=(ABS(ANGL)+0.5)
	L1ANGL=LANGL/100
	L2ANGL=(LANGL-100*L1ANGL)/10
	L3ANGL=LANGL-100*L1ANGL-10*L2ANGL
	DIST2(IS,JE-2)=INT(L1ANGL+1)
	DIST2(IS,JE-1)=INT(L2ANGL+1)	 
	DIST2(IS,JE  )=INT(L3ANGL+1)
C******  set up D6 under ANGL
	LD6=(ABS(D6)+0.05)*10.
	L1D6=LD6/100
	L2D6=(LD6-100*L1D6)/10
	L3D6=LD6-100*L1D6-10*L2D6
	DIST2(IS+1,JE-3)=INT(L1D6+1)
	DIST2(IS+1,JE-2)=INT(L2D6+1)
	DIST2(IS+1,JE-1)='.'
	DIST2(IS+1,JE  )=INT(L3D6+1)
C******  set up D1 under D6
	IF(D1.LT.0) DIST2(IS+2,JE-3)='-'
	LD1=(ABS(D1)+0.05)*10.
	L1D1=LD1/10
	L2D1=LD1-L1D1*10
	DIST2(IS+2,JE-2)=INT(L1D1+1)
	DIST2(IS+2,JE-1)='.'
	DIST2(IS+2,JE  )=INT(L2D1+1)
C******  set up D2 under D1
	IF(D2.LT.0) DIST2(IS+3,JE-3)='-'
	LD2=(ABS(D2)+0.05)*10.
	L1D2=LD2/10
	L2D2=LD2-L1D2*10
	DIST2(IS+3,JE-2)=INT(L1D2+1)
	DIST2(IS+3,JE-1)='.'
	DIST2(IS+3,JE  )=INT(L2D2+1)

	RETURN

	END
C====================================================================
C=================   SUBROUTINE MXMAP   =====================
C================== (updated on 8 Mar 88) ===================
C============================================================
C    CALLS to:	DATE		UGLINE		
C		TIME		UGTEXT
C		UGFONT		UGMARK
C		UGPICT		UGWRIT
C		UGDSPC		UGXERR
C		UGINIT
C============================================================

 	SUBROUTINE MXMAP(DIST2,MAXLEN,JMAX,NQ,FNAME,HELIX,H310,
	1	STRAND,TRN,TRN5,LOOP,SHEET,SHORT,
	2	LINES,IX,IY,BIT,COLR,NCNT,DEVICE,
	3	SEPMAX,SEPMIN,AXDMAX,AXDMIN,PT1MAX,PT1MIN,
	4	PT2MAX,PT2MIN)

C******  Subroutine to print out a line map as an approximation
C******  of the DIST2 character matrix.

	IMPLICIT NONE

	COMMON/SEGMENT/GRF

	INTEGER	     LINES,IS,JS,IE,JE,IYINCR,JXINCR,NPASS,ICOLOR
	INTEGER      GRF(500),I,J,N,JMAX,NQ,MAXLEN,JJ,ISEP
	INTEGER*4    IX(3,LINES),IY(3,LINES),BIT(3,LINES),NCNT(3)
	INTEGER*4    BBIT
	INTEGER	     NI,NJ,JBEGIN,IYSTRT,JXEND

	REAL*4       SEPMAX,SEPMIN,AXDMAX,AXDMIN,PT1MAX,PT1MIN
	REAL*4       PT2MAX,PT2MIN
	REAL*4	     X,Y,XMAX,YMAX,XSHFT,YSHFT,MATSIZ, XMIN, YMIN
	REAL*4	     XEND,YEND

	CHARACTER*3  B
	CHARACTER*5  BB
	CHARACTER*8  FNAME
	CHARACTER*8  PRTIME
	CHARACTER*9  PRDATE
	CHARACTER*7  COLOR,COLR(3,LINES),DEVICE

	BYTE	DIST2(MAXLEN,MAXLEN),TEST
	BYTE	HELIX,H310,STRAND,TRN,TRN5,LOOP,SHEET,SHORT

	IF(JMAX.GT.500) THEN 
		WRITE(2,5000)
5000		FORMAT(' JMAX.GT.500 - STOP')
		STOP
	ENDIF

	CALL DATE(PRDATE)
	CALL TIME(PRTIME)	

	XMAX=2870.	! XMAX,YMAX = maximum size of drawing space
	YMAX=2870.	! in drawing units
	XSHFT=200.	! XSHFT,YSHFT = origin of matrix in drawing
	YSHFT=50.	! space ( lower left corner )
	MATSIZ=2600.	! MATSIZ =  square matrix size in drawing units.

C******  Indices IS,IE,JS,JE refer to positions in the matrix DIST2(I,J)
C******  NQ = 1,2,3 identifies the panel to be drawn
C******  NI,NJ,IYINCR,JXINCR are given in drawing space units.

	IE=JMAX
	JE=JMAX
	IF(NQ.EQ.1.AND.JMAX.GT.260) THEN
		IE=260
		JE=260
	ENDIF
100	CONTINUE	! Set up indices for different NQ values
	IF(NQ.EQ.1) THEN
		IS=1
		JS=1
		NI=0
		NJ=0
		IYINCR=0
		JXINCR=0
	ELSEIF(NQ.EQ.2) THEN
		IS=1
		JS=241
		NI=0
		NJ=400
		IYINCR=0
		JXINCR=2400
	ELSEIF(NQ.EQ.3) THEN
		IS=241
		JS=241
		NI=400
		NJ=400
		IYINCR=2400
		JXINCR=2400
	ENDIF

150	CONTINUE	! Initialize drawing program

c*	CALL UGFONT('SIMPLEX')
c*	CALL UGPICT('CLEAR',0)
c*	CALL UGDSPC ('PUT',XMAX,YMAX,1.0)
c*	CALL UGINIT('CLEAR',GRF,500)

        write (8,1111)
	CALL UGLINE(1,  0.,  0.,0,0,GRF)	!
	CALL UGLINE(1,XMAX,  0.,1,0,GRF)	! Outline
	CALL UGLINE(1,XMAX,YMAX,1,0,GRF)	! drawing space
	CALL UGLINE(1,  0.,YMAX,1,0,GRF)	!
	CALL UGLINE(1,  0.,  0.,1,1,GRF)	!
	write (8,1114)
	write (8,1120)

1110	FORMAT('1.0 setlinewidth')
1111	FORMAT('newpath')
1114	FORMAT('closepath')
1115    FORMAT('1.5 setlinewidth')
1116    FORMAT('3.0 setlinewidth')
1120	FORMAT('stroke')
1121	FORMAT('showpage')

	XEND=MATSIZ+XSHFT
	YEND=MATSIZ+YSHFT

200     write (8,1111)
	CALL UGLINE(1, XSHFT, YEND,0,0,GRF)	!
	CALL UGLINE(1, XEND, YEND,1,0,GRF)	! Outline matrix 
	CALL UGLINE(1, XEND,YSHFT,1,0,GRF)	! 
	CALL UGLINE(1,XSHFT,YSHFT,1,0,GRF)	!
		Y=YEND-150.			!
		IF(NQ.EQ.2) Y=YEND
	CALL UGLINE(1,XSHFT,    Y,1,1,GRF)	!

	IF(NQ.NE.2) THEN			
		CALL UGLINE(1,XSHFT, YEND,0,0,GRF)	! Draw diagonal
		CALL UGLINE(1, XEND,YSHFT,1,1,GRF)	!
	ENDIF

300	CONTINUE 	! Put in tick marks on axes and diagonal
	IF(NQ.NE.2) THEN
		DO I=500-NI,MATSIZ,500			! Put in horizontal
			X=XSHFT				! tick marks for
			Y=YEND-FLOAT(I)			! blocks with
			CALL UGLINE(1,X,Y,0,0,GRF)	! diagonal triangle
			X=XSHFT+FLOAT(I)-50.
			CALL UGLINE(1,X,Y,1,1,GRF)
			X=XEND-50.
			CALL UGLINE(1,X,Y,0,0,GRF)
			CALL UGLINE(1,XEND,Y,1,1,GRF)
 		ENDDO
		DO J=500-NJ,MATSIZ,500			! Put in vertical
			Y=YSHFT				! tick marks for
			X=XSHFT+FLOAT(J)		! blocks with
			CALL UGLINE(1,X,Y,0,0,GRF)	! diagonal triangle
			Y=YEND-FLOAT(J)-120.
			CALL UGLINE(1,X,Y,1,1,GRF)
			Y=YEND-50.
			CALL UGLINE(1,X,   Y,0,0,GRF)
			CALL UGLINE(1,X,YEND,1,1,GRF)
		ENDDO

	ELSEIF(NQ.EQ.2) THEN
		NI=0
		NJ=400
		DO I=500-NI,MATSIZ,500			! Put in horizontal
			X=XSHFT				! tick marks for off
			Y=YEND-FLOAT(I)			! diagonal block
			CALL UGLINE(1,X,Y,0,0,GRF)
			X=XSHFT+50.
			CALL UGLINE(1,   X,Y,1,1,GRF)
			X=XEND-50.
			CALL UGLINE(1,   X,Y,0,0,GRF)
			CALL UGLINE(1,XEND,Y,1,1,GRF)
		ENDDO
		DO J=500-NJ,MATSIZ,500			! Put in vertical
			Y=YSHFT				! tick marks for off
			X=XSHFT+FLOAT(J)		! diagonal block
			CALL UGLINE(1,X,Y,0,0,GRF)
			Y=YSHFT+50.
			CALL UGLINE(1,X,Y,1,1,GRF)
			Y=YEND-50.
			CALL UGLINE(1,X,Y,0,0,GRF)
			CALL UGLINE(1,X,YEND,1,1,GRF)
		ENDDO
	ENDIF

	IF(NQ.NE.2) THEN
		DO I=1,259	! Put tick marks on diagonal
			X=XSHFT+FLOAT(10*I)-50.
			Y=YEND-FLOAT(10*I)
				CALL UGLINE(1,X,Y,0,0,GRF)
			X=X-50.
			IF(JMOD(I,10).EQ.0) X=X-50.
				CALL UGLINE(1,X,Y,1,1,GRF)
		ENDDO
	ENDIF
			write (8,1114)
		        write (8,1120)
			write (8,1111)

500	CONTINUE	! Axis labels
		X=XSHFT-50.
		IYSTRT=500
		IF(NQ.EQ.3) IYSTRT=1000
 	DO I=IYSTRT-NI,MATSIZ,500		! Put in labels on
		Y=YEND-FLOAT(I)			! Y axis
		WRITE(B,FMT='(I3)') (I+IYINCR)/10
C*		CALL UGTEXT('DSIZE=0.015,RIGHT',X,Y,B,GRF)
		X=XSHFT-75.
		CALL UGTEXT(0.015,X,Y,B,GRF)
	ENDDO
c		Y=YEND+50.
		Y=YEND+25.
		JXEND=0
		IF(NQ.EQ.1) JXEND=500
	DO J=500-NJ,MATSIZ-JXEND,500		! Put in labels on
c		X=XSHFT+FLOAT(J)		! X axis
		X=XSHFT+FLOAT(J)-25.
		WRITE(B,FMT='(I3)') (J+JXINCR)/10
C*		CALL UGTEXT('DSIZE=0.015,CENTER',X,Y,B,GRF)
		CALL UGTEXT(0.015,X,Y,B,GRF)
	ENDDO	

600	CONTINUE	! Put in header
	Y=YMAX-70.		
	X=XSHFT-150.
c*		CALL UGTEXT('DSIZE=0.025,LEFT',X,Y,FNAME,GRF)
		CALL UGTEXT(0.025,X,Y,FNAME,GRF)
	X=XSHFT+450.
	Y=YMAX-40.
c*		CALL UGTEXT('DSIZE=0.012,LEFT',X,Y,PRDATE,GRF)
		CALL UGTEXT(0.012,X,Y,PRDATE,GRF)
	Y=YMAX-100.
c*		CALL UGTEXT('DSIZE=0.012,LEFT',X,Y,PRTIME,GRF)
		CALL UGTEXT(0.012,X,Y,PRTIME,GRF)
	X=XSHFT+900.
	Y=YMAX-40.
c*		CALL UGTEXT('DSIZE=0.012,LEFT',X,Y,'MAXATOMS',GRF)
		CALL UGTEXT(0.012,X,Y,'MAXATOMS',GRF)
	X=XSHFT+970.
	Y=YMAX-100.
	WRITE(B,FMT='(I3)') JMAX
c*		CALL UGTEXT('DSIZE=0.012,LEFT',X,Y,     B,GRF)
		CALL UGTEXT(0.012,X,Y,B,GRF)
	X=XSHFT+1200.
	Y=YEND+120.
c*		CALL UGTEXT('DSIZE=0.015,LEFT',X,Y,'ALPHA CARBON NUMBER',GRF)
		CALL UGTEXT(0.015,X,Y,'ALPHA CARBON NUMBER',GRF)
	IF(NQ.GT.1) GOTO 800	! Leave out parameters in panels after 
C				! first

	X=XSHFT+2125.
	Y=YMAX-34
c		CALL UGTEXT('DSIZE=0.008,LEFT',X,Y,'AT.SEPN.',GRF)
		CALL UGTEXT(0.008,X,Y,'AT.SEPN.',GRF)
	Y=YMAX-85
c		CALL UGTEXT('DSIZE=0.008,LEFT',X,Y,'AX.DIST.',GRF)
		CALL UGTEXT(0.008,X,Y,'AX.DIST.',GRF)
	Y=YMAX-136
c		CALL UGTEXT('DSIZE=0.008,LEFT',X,Y,' POINT 1',GRF)
		CALL UGTEXT(0.008,X,Y,' POINT 1',GRF)
	Y=YMAX-187
c		CALL UGTEXT('DSIZE=0.008,LEFT',X,Y,' POINT 2',GRF)
		CALL UGTEXT(0.008,X,Y,' POINT 2',GRF)
	X=XSHFT+2350
	Y=YMAX-34
		WRITE(BB,FMT='(F5.1)') SEPMAX
c		CALL UGTEXT('DSIZE=0.008,LEFT',X,Y,BB,GRF)		
		CALL UGTEXT(0.008,X,Y,BB,GRF)
	X=XSHFT+2500
		WRITE(BB,FMT='(F5.1)') SEPMIN
c		CALL UGTEXT('DSIZE=0.008,LEFT',X,Y,BB,GRF)		
		CALL UGTEXT(0.008,X,Y,BB,GRF)
	X=XSHFT+2350
	Y=YMAX-85
		WRITE(BB,FMT='(F5.1)') AXDMAX
c		CALL UGTEXT('DSIZE=0.008,LEFT',X,Y,BB,GRF)		
		CALL UGTEXT(0.008,X,Y,BB,GRF)
	X=XSHFT+2500
		WRITE(BB,FMT='(F5.1)') AXDMIN
c		CALL UGTEXT('DSIZE=0.008,LEFT',X,Y,BB,GRF)		
		CALL UGTEXT(0.008,X,Y,BB,GRF)
	X=XSHFT+2350
	Y=YMAX-136
		WRITE(BB,FMT='(F5.1)') PT1MAX
c		CALL UGTEXT('DSIZE=0.008,LEFT',X,Y,BB,GRF)		
		CALL UGTEXT(0.008,X,Y,BB,GRF)
	X=XSHFT+2500
		WRITE(BB,FMT='(F5.1)') PT1MIN
c		CALL UGTEXT('DSIZE=0.008,LEFT',X,Y,BB,GRF)		
		CALL UGTEXT(0.008,X,Y,BB,GRF)
	X=XSHFT+2350
	Y=YMAX-187
		WRITE(BB,FMT='(F5.1)') PT2MAX
c		CALL UGTEXT('DSIZE=0.008,LEFT',X,Y,BB,GRF)		
		CALL UGTEXT(0.008,X,Y,BB,GRF)
	X=XSHFT+2500
		WRITE(BB,FMT='(F5.1)') PT2MIN
c		CALL UGTEXT('DSIZE=0.008,LEFT',X,Y,BB,GRF)		
		CALL UGTEXT(0.008,X,Y,BB,GRF)

	write (8,1114)
	write (8,1120)

	write (8,1115)
	write (8,1111)
	
C**********************  ENTER DATA FROM DIST2  ***********************

800	IF(NQ.EQ.2) GOTO 1140

        CALL UGCOMMENT ('starting for putting cis peptides')

	COLOR='WHITE  '		! List cis peptides on the diagonal
	DO I=IS,IE
		IF(DIST2(I,I).EQ.'C') THEN
			X=XSHFT+FLOAT(10*I-JXINCR)
			Y=YEND -FLOAT(10*I-IYINCR)
c			CALL UGMARK('MARK=3,DSIZE=0.006,COLOR',X,Y,GRF)
			CALL UGMARK(3,0.006,X,Y,GRF)
c			type *,'mark=3'
		ENDIF
	ENDDO

	TEST = HELIX		! Set to draw helix

	CALL UGCOMMENT ('starting for putting HELIX ')

	do i=is,ie		!Put helix on diagonal
                DO J=I+1,JE+1   ! J-1 will = end of symbol line
			IF(DIST2(I,J).EQ.TEST) THEN
                                X=XSHFT+FLOAT(10*I-JXINCR)
                                Y=YEND - FLOAT(10*I-IYINCR)
                                CALL UGLINE(5,X,Y,0,0,GRF)
                                X=XSHFT+FLOAT(10*(J-1)-JXINCR)
                                CALL UGLINE(5,X,Y,1,0,GRF)
                        ENDIF

                ENDDO
	enddo
        call ugline(5,x,y,0,1,grf)

	call ugcomment ('putting in h310 on diag')
        TEST = H310            ! Set to draw h310

        do i=is,ie              !Put h31 on diagonal
                DO JJ=J,JE      ! Put in H310 line where required
                        IF(DIST2(I,JJ).EQ.H310) THEN    !JJ will = end
                        ELSE                            !of H310 line
                                X=XSHFT+FLOAT(10*(JJ-1)-JXINCR)
                                CALL UGLINE(4,X,Y,1,1,GRF)
                        ENDIF
                ENDDO
	enddo
	call ugline(4,x,y,0,1,grf)

        call ugcomment ('putting in strand on diag')
        TEST = STRAND            ! Set to draw strand

        do i=is,ie              !Put strand on diagonal
                DO J=I+1,JE+1   ! J-1 will = end of symbol line
                        IF(DIST2(I,J).EQ.TEST) THEN
                                X=XSHFT+FLOAT(10*I-JXINCR)
                                Y=YEND - FLOAT(10*I-IYINCR)
                                CALL UGLINE(3,X,Y,0,0,GRF)
                                X=XSHFT+FLOAT(10*(J-1)-JXINCR)
                                CALL UGLINE(3,X,Y,1,0,GRF)
                        ENDIF
                ENDDO
        enddo
	call ugline(3,x,y,0,1,grf)


	call ugcomment ('Putting TRN triangle on diagonal')
	COLOR='RED    '
	icolor=2
	DO I=IS,IE		! Put triangle on diagonal for TRN
		J=I+3
		IF(DIST2(I,J).EQ.TRN) THEN
			X=XSHFT+FLOAT(10*I-JXINCR)
			Y=YEND -FLOAT(10*I-IYINCR)
			CALL UGLINE(1,X,Y,0,0,GRF)
			X=X+30.
			CALL UGLINE(icolor,X,Y,1,0,GRF)
			Y=Y-30.
			CALL UGLINE(icolor,X,Y,1,1,GRF)
		ENDIF
	ENDDO

	call ugcomment ('putting TRN5 triangle on diagonal')
	DO I=IS,IE		! Put triangle on diagonal for TRN5
		J=I+4
		IF(DIST2(I,J).EQ.TRN5) THEN
			X=XSHFT+FLOAT(10*I-JXINCR)
			Y=YEND -FLOAT(10*I-IYINCR)
			CALL UGLINE(1,X,Y,0,0,GRF)
			X=X+40.
			CALL UGLINE(icolor,X,Y,1,0,GRF)
			Y=Y-40.
			CALL UGLINE(icolor,X,Y,1,1,GRF)
		ENDIF
	ENDDO

c       write (8,1114)                                                         
c       write (8,1120)

c	write (8,1111)

C******  Start here for off-diagonal elements
	call ugcomment ('starting off-diagonal elements')
1140	COLOR='WHITE  '
	icolor=2
	call ugcomment ('put in lines and marks for BSHEET')
	DO 1300 I=IS,IE		! Put in lines and marks for BSHEET
		J=I+5
		IF(NQ.EQ.2) J=JS
1150		IF(J.GT.JE) GOTO 1300
		IF(DIST2(I,J).NE.SHEET.AND.DIST2(I,J)
	1		.NE.SHORT) THEN	! Set J=start of SHEETor SHORT
			J=J+1
			GOTO 1150
		ELSE
			DO JJ=J,JE+1
				IF(DIST2(I,JJ).EQ.SHEET.OR.DIST2(I,JJ)
	1				.EQ.SHORT) THEN	!Set JJ= other 
C						! end of SHEET or SHORT
					IF(JJ.GE.JE+1) GOTO 1152
				ELSE
1152					X=XSHFT + FLOAT(10*J-JXINCR)
					Y=YEND  - FLOAT(10*I-IYINCR)
					CALL UGLINE(1,X,Y,0,0,GRF)
					X=XSHFT+FLOAT(10*(JJ-1)-JXINCR)
					CALL UGLINE(icolor,X,Y,1,1,GRF)
					GOTO 1200
				ENDIF
			ENDDO
1200			J=JJ
			GOTO 1150
		ENDIF
1300	CONTINUE

	call ugcomment (' Put in mardkers for short sheet distances')
	DO I=IS,IE	!Put in markers for short distances in sheets
		JBEGIN=I
		IF(NQ.EQ.2) JBEGIN=JS
		DO J=JBEGIN,JE
			IF(DIST2(I,J).EQ.SHORT) THEN
				X=XSHFT+FLOAT(10*J-JXINCR)
				Y=YEND-FLOAT(10*I-IYINCR)
c				CALL UGMARK('MARK=1,DSIZE=0.004,
c	1				COLOR',X,Y,GRF)
				CALL UGMARK(1,0.004,x,y,grf)
 			ELSEIF(DIST2(I,J).EQ.'<') THEN	!Entry for LOOPS
 				X=XSHFT+FLOAT(10*J-JXINCR)	
 				Y=YEND- FLOAT(10*I-IYINCR)
c				CALL UGMARK('MARK=2,DSIZE=0.004,
c	1				COLOR',X,Y,GRF)
 				CALL UGMARK(2,0.004,X,Y,GRF)
			ENDIF

		ENDDO
	ENDDO

	call ugcomment ('Putting in supersecondary boxes')
	DO I=1,NCNT(NQ)		! Put in supersecondary boxes
		X=XSHFT+FLOAT(IY(NQ,I)-JXINCR)
		Y=YEND-FLOAT(IX(NQ,I)-IYINCR)
		BBIT=BIT(NQ,I)
		COLOR=COLR(NQ,I)
		IF (COLOR. EQ. ' ') ICOLOR=1
		IF (COLOR. EQ. 'RED  ') ICOLOR=1
		IF (COLOR. EQ. 'WHITE  ') ICOLOR=2
		IF (COLOR. EQ. 'BLUE  ') ICOLOR=3
		IF (COLOR. EQ. 'GREEN  ') ICOLOR=4
c		write(8,1112) x,y
c1112    FORMAT(F5.0,1X,F5.0,' moveto')
		CALL UGLINE(icolor,X,Y,BBIT,1,GRF)
	ENDDO

C****************************************************************

c	CALL UGWRIT(' ',0,GRF)
	
c	IF(DEVICE.EQ. 'TEK4010 ') ACCEPT 10,B
10	FORMAT(A)

	write (8,1114)
	write (8,1120)
	WRITE (8,1121)
        GOTO 10000

10000	END

C========================   UGXERR   ==================================

D	SUBROUTINE UGXERR(LEVEL,SNAME,INDEX)
C******  This subroutine is the error processor for UGSYS.
C******  level indicates the level number of the error received.
C******  sname name of the subrutine which encountered the error
C******  index index for the error. See UGSYS manual.
	
D	IMPLICIT NONE
D	INTEGER*4 LEVEL, INDEX
D	CHARACTER*8 SNAME
D	COMMON/SEGMENT/GRF
D	INTEGER*4 GRF(500)

D	IF(INDEX.EQ.11) THEN
D		CALL UGWRIT(' ',0,GRF)
D		CALL UGINIT('CONTINUE',GRF,500)
D		LEVEL=0
D	ENDIF

D	RETURN

D	END

C=======================================================================


	SUBROUTINE UGCOMMENT (string)
	character*(*) string

	WRITE (8,100) string
100	FORMAT('%',A)
	return
	end
	
C====================================================================
C=================   SUBROUTINE UGLINE   =====================
C================== (created 2/10/94)    =====================
C=============================================================
	SUBROUTINE UGLINE(ICOLOR,X,Y,II,IJ,GRF)
c**	produce line in the PS file
c**	II	define "move to" or " line to" (0=moveto, 1=lineto)
c**	IJ	define  finished or unfinished line(0=unfinished, 1=finished)
c**	COLOR	define line width
	
        INTEGER      GRF(500),I,J,N,JMAX,NQ,MAXLEN,JJ,ISEP
        INTEGER      II,ICOLOR,IJ
        REAL*4       X,Y
	real*4       colors (5)
        data         colors / 0.5, 1.0, 4.0, 6.0, 9.0/


	if (II .EQ. 0) then
	        write (8,1112) X,Y
        else
		write (8,1113) X,Y
		write (8,1115) colors(icolor)      !set linewidth 
	endif	
	if (IJ .EQ. 1) then
	   	write (8,1120)                    !stroke
 	   	write (8,1130)                    !newpath
 	   	write (8,1112) X,Y                !moveto here
	endif
	return

1112    FORMAT(F5.0,1X,F5.0,' moveto')
1113    FORMAT(F5.0,1X,F5.0,' lineto')
1115    FORMAT(F4.2, ' setlinewidth')
1120    FORMAT('stroke')
1130	format('newpath')
	end

C====================================================================
C=================   SUBROUTINE UGMARK   =====================
C================== (created 2/23/94)    =====================
C=============================================================
        SUBROUTINE UGMARK(II,SIZE,X,Y,GRF)
c**	subroutine to produce different marks in the PS file

        INTEGER      GRF(500),I,J,N,JMAX,NQ,MAXLEN,JJ,ISEP
        INTEGER      II                                                         
        REAL*4       X,Y,DSIZE,SIZE

        xmax=2870.
        dsize=size*xmax*0.5

c**	selected digit is:
c**	1	diagonal cross
c**	2	diamond
c**	3	square

        if (II .EQ. 1) then
		write (8,1111)
                write (8,1112) x-dsize,y+dsize
		write (8,1113) x+dsize,y-dsize
		write (8,1112) x-dsize,y-dsize
		write (8,1113) x+dsize,y+dsize
		write (8,1114)
		write (8,1115)
        else
         if (II .EQ. 2) then
                write (8,1111)
                write (8,1112) x-dsize,y
		write (8,1113) x,y+dsize
		write (8,1113) x+dsize,y
		write (8,1113) x,y-dsize
		write (8,1113) x-dsize,y
                write (8,1114)
                write (8,1115)
	 else
	   if (II .EQ. 3) then
                write (8,1111)
                write (8,1112) x-dsize,y+dsize
                write (8,1113) x+dsize,y+dsize
                write (8,1113) x+dsize,y-dsize
                write (8,1113) x-dsize,y-dsize
                write (8,1113) x-dsize,y+dsize                
                write (8,1114)
                write (8,1115)
           endif
	 endif
	endif

1111	format('newpath')
1112    FORMAT(F5.0,1X,F5.0,' moveto')
1113    FORMAT(F5.0,1X,F5.0,' lineto')
1114	format('2 setlinewidth')
1115    format('stroke')

1200    end

C====================================================================
C=================   SUBROUTINE UGTEXT   =====================                  
C================== (created 2/15/94)    =====================
C=============================================================
        SUBROUTINE UGTEXT(SIZE,X,Y,TEXT,GRF)
c**	subroutine to produce text in the PS file
c**	size: relative size of text

        INTEGER      GRF(500),I,J,N,JMAX,NQ,MAXLEN,JJ,ISEP                      
        REAL*4       XMAX,YMAX,X,Y,SIZE,DSIZE
	CHARACTER*(*)  TEXT

	xmax=2870.
	dsize=size*xmax
	write (8,1122)
	write (8,1123) dsize
	write (8,1124)
	write (8,1112) x,y
	write (8,1125) text

1112    FORMAT(F5.0,1X,F5.0,' moveto')
1122	FORMAT ('/Times-Roman findfont')
1123	FORMAT (F5.0,' scalefont')
1124	FORMAT ('setfont')
1125	FORMAT ('(',A,')',' show')
	end
