                            program HELANAL
C****************************************************************************
C                                                                           *
C PROGRAM TO CHARACTERISE THE GEOMETRIES OF ALPHA HELICES IN PROTEINS.      *
C                                                                           *
C AUTHORS: Sandeep Kumar and Prof. Manju Bansal,                            *
C          MBU, Indian Institute of Science, Bangalore 560012, India        *
C                                                                           *
C  e-mail address: mb@mbu.iisc.ernet.in                                     *
C                                                                           *
C Method outlined in the following paper:                                   *
C                                                                           *
C Kumar, S. and Bansal, M. (1996). Structural and Sequence Characteristics  * 
C of Long Alpha Helices in Globular Proteins. Biophysical J.,71, 1574-1586. *
C                                                                           *
C If any Bug is found, please report it to the authors.                     *
C                                                                           *
C SUMMARY OF THE ALGORITHM                                                  *
C                                                                           *
C Geometry of an alpha helix is characterised in terms of the angles between*
C local helix axes and the path traced by the local helix origins,          *
C calculated using the procedure of Sugeta and Miyazawa (1967), for every   *
C set of four contiguous C-alpha atoms, and sliding this window over the    *
C length of the helix in steps of one C-alpha atom.                         *
C Matrix M(I, J) contains the bending angles between local helix axes I & J *
C which are used to characterize the overall geometry of the  helix.        * 
C The local helix origins trace out the path described by the helix in 3-D  *
C space. These origins are reoriented in X-Y plane and the reoriented       *
C points are used to fit a circle and a line by least squares method.       *
C Unit twist and unit height of the alpha helix are also calculated.        *
C                                                                           *
C A maximum of 5000 helices, each with 100 or less number of residues can   *
C be analysed at one time.                                                  * 
C                                                                           *
C In order to analyse a larger number of helices, increase the value of 'i' *
C in the  following statement:                                              *
C                                                                           *
C      parameter (i=5000) in the main program.                              *
C                                                                           *
C In order to analyse helices with lengths greater than 100 residues,       *
C increase the dimensions of the appropriate variables in the dimension     *
C statements of the main program and subroutines.                           *
C                                                                           *
C This program has several options for input and is fully interactive.      *
C                                                                           *
C INPUT FILES :                                                             *
C                                                                           *
C Input files to this program can be different depending upon the options   *
C chosen. This program can directly read HELIX records in one or more PDB   *
C files. In order to analyse helices found in helix records of a PDB file,  *
C give the PDB file name as input to the program. In order to analyse the   *
C helices found in the HELIX records of more than one PDB  files, give the  *
C PDB file names sequentially when prompted by the program or input the     *
C name of the file containing the PDB files names in format (5x,a11) as the *
C input to the program.                                                     *
C                                                                           *
C HELIX records can be read from PDB files or can be read from a file       *
C      containing information about the helix start/end residues written in *
C      the same format as PDB HELIX records or in a different format, in    *
C      which case the format has to be keyed in when running the program.   *
C In the last case, files other than the PDB files or files containing      *
C      C-alpha coordinates in the PDB format, can also be given as input to *
C      the program, upon specifying their format.                           * 
C                                                                           *
C                                                                           *
C All PDB files and other input files (if any) should be in the same        *
C directory.                                                                *
C                                                                           *
C OUTPUT FILES :                                                            *
C                                                                           *
C For each  run of HELANAL on file(s) containing alpha helices with length  *
C greater than or equal to 9 residues, the following output files are       *
C created.                                                                  *
C                                                                           *
C RUN.ANS contains the questions and their answers during a run of HELANAL. *
C                                                                           *
C HELINFO.OUT file created only when the HELIX records in the PDB files are *
C used for information on helix start/end residues, contains these helix    *
C records.                                                                  *
C                                                                           *
C HELCA.OUT contains Coordinates of the C-alpha atoms constituting the      *
C helices.                                                                  *
C                                                                           *
C AXES.OUT contains the local helix axes fitted to 4 consecutive C-alpha    *
C atoms along with a matrix  M(I, J) whose elements are the angles between  *
C local helix axes I and J.                                                 *
C                                                                           *
C ANGLE.OUT contains the angle between successive local helix axes. It also *
C lists mean bending angle and maximum bending angle for each helix.        *
C                                                                           *
C ORIGIN.OUT contains the local helix origins for the helix along with the  *
C statistics obtained by fitting least square plane, circle and line to the *
C local helix origins.                                                      *
C                                                                           *
C NH.OUT contains unit height and unit twist for every turn of the helix as *
C well as average unit height and average unit twist for the whole helix.   *  
C                                                                           *
C ***.PRM contains summary of various parameters obtained by HELANAL.       *
C                                                                           *
C ***.TAB contains the parameters from ***.PRM file in a tabular form along *
C the overall geometry assignment.                                          *
C                                                                           *
C****************************************************************************
c
           parameter (i=5000)
           character*80 line, prmfile,forma,formacor,tabfile
           character*40 inpfile
           character*11 fnm,pdbfile
           character*3 rb,re,helnm,res
           character*4 atom
           character*1 cb,ce,ans1,ans2,ans3,ans4,ans5,chain
           character*1 geom
c          
           dimension rb(i),cb(i),nrb(i)
           dimension re(i),ce(i),nre(i),fnm(i)
           dimension helnm(i),pdbfile(i)
           dimension xca(100),yca(100),zca(100)
           dimension ox(100),oy(100),oz(100)
           dimension twist(100),rnou(100),angl(100)
           dimension vec12(3),vec23(3),vec34(3),dv13(3),dv24(3)
           dimension uloc1(100),uloc2(100),uloc3(100),kount(100)
           dimension rad(3),ibendangl(100,100),height(100)
c
c            FORMATS USED
c
    1    format (/1x,'Do You wish to use HELIX records in PDB files?')
    2   format (/1x,'Key in name of the INPUT file containing PDB file n
     -ames in (5x,a11) format',/1x,'e.g. ***.fil :')
    3     format (/1x,'Key in name of the INPUT file containing helix st
     -art/end information',/1x,'e.g. ***.inp :')
    4      format (1x,'C-ALPHA ATOMS',/)
    5      format (30x,3f8.3) 
    6      format (1x,'LOCAL HELIX AXES FOR CA ATOMS (I - (I+3))',3x,
     -     'L',7x,'M',7x,'N',/)
    7      format (5x,2(i5,5x),15x,3f8.3)
    8      format (5x,3(f8.3))
    9      format (1x,'LOCAL HELIX ORIGINS FOR CA ATOMS (I - (I+3))'
     -     ,2x,'X',7x,'Y',7x,'Z',/)
   10   format(/1x,'Is the information in the same format as HELIX recor
     -ds in PDB files?')
   11      format (5x,'O ',3(i5,2x),15x,3f8.3)
   12      format (1x,'LOCAL BENDING ANGLES CALCULATED FROM LOCAL HELIX 
     -AXES FITTED TO CA ATOMS'/1x,'((I-3) - I)  AND (I - (I+3))'/)
   13      format (4(i5,2x),5x,f8.3)
   14      format (1x,'LOCAL HELIX PARAMETERS CALCULATED FOR CA ATOMS (I
     - - (I+3))')
   16      format (1x,'MATRIX M(I,J) FOR ANGLES BETWEEN LOCAL HELIX AXES
     - I AND J')
   17      format (i3,1x,40i3)
   18      format (4x,40i3)
   19      format(1x,'AVERAGE TWIST (DEG.)',6x,f8.2,2x,'S.D.',f8.2,
     -     5x,'ADV',f8.2) 
   20      format(1x,'AVERAGE N',17x,f8.2,2x,'S.D.',f8.2,5x,'ADV'
     -     ,f8.2) 
   21      format(1x,'AVERAGE UNIT HEIGHT (ANG.)',f8.2,2x,'S.D.',f8.2,
     -      5x,'ADV',f8.2) 
   22 format(/1x,'If not, key in the input data format:'/1x,'(Specify fi 
     -elds for coordinate file name, helix name, name, chain and number'
     -/1x,'for first and last residues. Use three letter code for residu
     -e names.)'/)
   23      format(a80)
   24     format (/1x,'Do you wish to type in the PDB file names?')
   25      format(5x,'N------>C'/)
   26      format(1x,'MAXIMUM BENDING ANGLE',5x,f8.3)
   27      format(1x,'MEAN BENDING ANGLE',8x,f8.3,2x,'S.D.',f8.3,
     -     5x,'ADV',f8.3) 
   28      format(1x,'SUMMARY OF PARAMETERS OBTAINED BY HELANAL'/) 
   29      format(1x,'MEAN BENDING ANGLE (DEG.)',15x,f8.3,2x,'S.D.',
     -      f8.3,2x,'ADV',f8.3) 
   30      format(1x,'MAXIMUM BENDING ANGLE (DEG.)',12x,f8.3)
   31      format (42x,'TWIST',6x, 'N',4x,'HEIGHT'/)
   32      format(a80)
   33      format(13x,a3,1x,a3,1x,a1,1x,i3,4x,3f8.3)
   34      format(11x,a3,1x,2(a3,1x,a1,1x,i4,2x),33x,a4)
   35      format(/1x,a11,2x,'HELIX',2x,a3,1x,2(a3,1x,a1,1x,i3,4x))
   36      format(1x,'AVERAGE N',31x,f8.3,2x,'S.D.',f8.3,2x,'ADV'
     -     ,f8.3) 
   37      format(1x,'AVERAGE UNIT HEIGHT (ANG.)',14x,f8.3,2x,'S.D.',
     -     f8.3,2x,'ADV',f8.3) 
   38      format(1x,'LEAST SQUARES PLANE, CIRCLE AND LINE FITTED TO LOC
     -AL HELIX ORIGINS')     
   39      format(5x,a11)
   40   format(/1x,'Key in a PDB file name :')
   42   format(/1x,'QUESTIONS AND THEIR ANSWERS'/)
   43   format(/1x,'All questions and their answers are written into the 
     - file "RUN.ANS". '/)
   44   format(/1x,'Is',1x,a11,1x,'a PDB file?')
   45   format(/1x,'If not, key in the input data format:',/1x,'(Specify 
     - fields for atom name (4 characters), res. name, chain, res. numbe
     -r',/1x,'and XYZ coordinates. Use three letter code for residue nam
     -es.)'/)
   46   format(a40)
   47   format(1x,'HELICES TAKEN FROM :',1x,a40)
   48   format(/1x,'WRITING THE RESULTS TO :',1x,a40,/,1x,'AND IN TABULA
     -R FORM TO :',1x,a40)
   49   format(1x,a40)
   50   format(1x,'HELIX IGNORED : IT CONTAINS LESS THAN 9 RESIDUES.')
   51   format(/1x,'Do you wish to use another PDB file?')
   52   format(/1x,a11,' DOES NOT CONTAIN HELIX RECORDS')
   53  format(3x,'File Name',6x,'Helix',5x,'n',4x,'h',3x,'Aver',2x,
     - 'Max',2x,'Radius',2x,'rmsC',2x,'rmsL',2x,'r2',1x,'Geometry',
     -/38x,'BA',3x,'BA',/)
   54    format(1x, 'Helix     : Chain identifier, number and name of fi
     -rst and last residues'/13x, 'of the helix.')
   55    format (1x,'n         : Average number of residues per turn of 
     -the helix.')
   56 format(1x,'h         : Average unit height in the helix (Angstroms
     -).')
   57    format(1x,'Aver BA   : Average Bending angle between successive 
     - local helix axes (Deg.).')
   58  format(1x,'Max BA    : Maximum Bending angle between succesive lo
     -cal helix axes (Deg.).')
   59  format(1x,'Radius    : Radius of least squares circle fitted to t
     -he local helix '/13x,'origins (Angstroms).')
   60 format(1x,'rmsC      : Root Mean Square Deviation for least squares
     - circle fitted to'/13x,'the local helix origins (Angstroms).')
   61 format(1x,'rmsL      : Root Mean Square Deviation for least squares 
     - line fitted to the'/13x,'local helix origins (Angstroms).
     -')
   62   format(1x,'r2        : Square of linear correlation coefficient  
     -for least squares line'/13x,'fittted to the local helix 1origins.
     -')
   63   format(1x,'Geometry  : Overall geometry of the helix, Linear (L)
     -, Curved (C),'/13x,'Kinked (K) or unassigned (*).')
   64   format(1x,'Std. Dev. : Standard Deviations of the average parame 
     -ters for the helix'/13x,'in the previous line.')
   66   format(/1x,'DESCRIPTION OF THE TABLE'/)
   67   format(/'*******************************************************
     -************************')
   68    format(/,5x,'NL=',i5,5x,'NC=',i5,5x,'NK=',i5,5x,'NA=',i5
     -,5x,'NH=',i5)
   69 format(1x,'NL, NC, NK: Number of helices with overall linear, curv
     -ed, kinked,',/1x,'NA, NH',6x,'unassigned geometry and total number
     - of helices.')

               nlinear=0
               ncurved=0
               nkinked=0
               nogeom=0
c
C OPENING INPUT AND OUTPUT FILES
c
           open (unit=1,file='helca.out')
           open (unit=2,file='axes.out')
           open (unit=3,file='origin.out')
           open (unit=7,file='RUN.ANS')
           open (unit=8,file='angle.out') 
           open (unit=9,file='nh.out')
              write(1,4)
              write(2,6)
              write(3,9)
              write(7,42)
              write(8,12)
              write(9,14)
              write(9,31)
              write(6,43)

c  ASKING THE FIRST QUESTION : USE HELIX RECORDS IN PDB FILE?

  120        write (6,1)
             write (7,1)
             read(5,'(a1)')ans1
             write(7,'(1x,a1)')ans1
             if ((ans1.ne.'n').and.(ans1.ne.'N').and.(ans1.ne.'y')
     -   .and.(ans1.ne.'Y')) then
                 goto 120
               endif
c
c READING HELIX RECORDS FROM PDB FILES
c
             if ((ans1.eq.'y').or.(ans1.eq.'Y')) then

c ASKING SECOND QUESTION : TYPE IN THE PDB FILE NAMES ?

 119           write(6,24)
               write(7,24)
             read(5,'(a1)') ans2
             write(7,'(1x,a1)') ans2
             if ((ans2.ne.'n').and.(ans2.ne.'N').and.(ans2.ne.'y')
     -   .and.(ans2.ne.'Y')) then
                 goto 119
               endif

            if((ans2.eq.'n').or.(ans2.eq.'N')) then
  121          write(6,2)
               write(7,2)
              read(5,46) inpfile
              write(7,49) inpfile
    
              kkk = index(inpfile(1:40),':')
             if (kkk.eq.0) then
                 kkk=1
               else
                kkk=kkk+1
              endif
              kkkk=index(inpfile(kkk:40),'.')
              kkkk=kkkk+kkk
             prmfile=inpfile(kkk:kkkk-1)//'prm'
             tabfile=inpfile(kkk:kkkk-1)//'tab'
             
      if((inpfile(kkkk:kkkk+2).eq.'ent').or.(inpfile(kkkk:kkkk+2)
     -.eq.'pdb').or.(inpfile(1:3).eq.'pdb')) goto 121            

           open (unit=12,file=inpfile,status='OLD',err=121)
           open (unit=11,file='helinfo.out')
           open (unit=10,file=prmfile)
           open (unit=13,file=tabfile)

            write(13,53)
                nfiles=0
              do 501 j=1,i
            read(12,39,end=129,err=121) pdbfile(j)
                nfiles=nfiles+1
  501          continue
  129      continue 
              
             do  502 l=1,nfiles
            open (unit=15,file=pdbfile(l),status='OLD',err=121)
                do 503 j=1,i
            read(15,32,end=131) line
            if (line(1:7).eq.'HELIX  ') then
            write(11,32) line 
              endif
503         continue
131      continue
  
             close (15)
502      continue
              rewind(11)
           write(1,47) inpfile
           write(2,47) inpfile
           write(3,47) inpfile
           write(8,47) inpfile
           write(9,47) inpfile
              write(10,28)
           write(10,47) inpfile
           write(6,48) prmfile,tabfile
           write(7,48) prmfile,tabfile
               goto 100

             endif
             endif

 155          if((ans1.eq.'y').or.(ans1.eq.'Y')) then
              if((ans2.eq.'y').or.(ans2.eq.'Y')) then
 133          write(6,40)
              write(7,40)
              read(5,46) inpfile
              write(7,49) inpfile


              kkk = index(inpfile(1:40),':')
             if (kkk.eq.0) then
                 kkk=1
               else
                kkk=kkk+1
              endif
              kkkk=index(inpfile(kkk:40),'.')
              kkkk=kkkk+kkk

      if((inpfile(kkkk:kkkk+2).ne.'ent').and.(inpfile(kkkk:kkkk+2)
     -.ne.'pdb')) goto 133            
             
             prmfile=inpfile(kkk:kkkk-1)//'prm'
             tabfile=inpfile(kkk:kkkk-1)//'tab'


           open (unit=12,file=inpfile,status='OLD',err=133)
           open (unit=11,file='helinfo.out')
           open (unit=10,file=prmfile)
           open (unit=13,file=tabfile)

            write(13,53)
                      nlines = 0
                   do 497 j=1,i
            read(12,32,end=132) line
            if (line(1:7).eq.'HELIX  ') then
            write(11,32) line 
                 nlines=nlines+1
              endif
497         continue
132      continue

             if(nlines.eq.0) then
           write(6,52) inpfile
           write(7,52) inpfile
               goto 156
               endif

              rewind(11)
            close (12)
           write(1,47) inpfile
           write(2,47) inpfile
           write(3,47) inpfile
           write(8,47) inpfile
           write(9,47) inpfile
              write(10,28)
           write(10,47) inpfile
           write(6,48) prmfile,tabfile
           write(7,48) prmfile,tabfile

            go to 100

                endif
                endif


c   READING HELIX INFORMATION FROM A SEPARATE FILE
 
             if ((ans1.eq.'n').or.(ans1.eq.'N')) then
  122        write (6,3)
             write (7,3)
             read(5,46) inpfile
             write(7,49) inpfile

              kkk = index(inpfile(1:40),':')
             if (kkk.eq.0) then
                 kkk=1
               else
                kkk=kkk+1
              endif
              kkkk=index(inpfile(kkk:40),'.')
              kkkk=kkkk+kkk
             prmfile=inpfile(kkk:kkkk-1)//'prm'
             tabfile=inpfile(kkk:kkkk-1)//'tab'
      if((inpfile(kkkk:kkkk+2).eq.'ent').or.(inpfile(kkkk:kkkk+2)
     -.eq.'pdb').or.(inpfile(1:3).eq.'pdb')) goto 122            


           open (unit=11,file=inpfile,status='OLD',err=122)
           open (unit=10,file=prmfile)
           open (unit=13,file=tabfile)

            write(13,53)
           write(1,47) inpfile
           write(2,47) inpfile
           write(3,47) inpfile
           write(8,47) inpfile
           write(9,47) inpfile
              write(10,28)
           write(10,47) inpfile

C ASKING THE THIRD QUESTION : IS THE FORMAT SAME AS HELIX RECORDS IN PDB FILES?
c
 123          write(6,10)
              write(7,10)
              read(5,'(a1)') ans3
             write(7,'(1x,a1)') ans3
             if ((ans3.ne.'n').and.(ans3.ne.'N').and.(ans3.ne.'y')
     -   .and.(ans3.ne.'Y')) then
                 goto 123
               endif


                if ((ans3.eq.'N').or.(ans3.eq.'n')) then
 127            rewind(11)
                write(6,22)
                write(7,22)
                read(5,'(a80)') forma 
                write(7,'(1x,a79)') forma 


                   indexhel=0
                 do 498 j=1,i
          read(11,forma,end=126,err=127) fnm(j),helnm(j),rb(j),cb(j),
     1          nrb(j),re(j),ce(j),nre(j)
                    indexhel=indexhel+1
 499          continue     
 498          continue     
 126          continue     

               go to 125
               endif
c
                if((ans3.eq.'Y').or.(ans3.eq.'y')) then

              goto 100
            endif
             endif
c
c READING HELIX INFORMATION
c
  100              indexhel=0
             do 504 j=1,i
             read(11,34,end=99) helnm(j),rb(j),cb(j),nrb(j),
     1          re(j),ce(j),nre(j),fnm(j)(4:7)
                  do 505 k=5,7
               if(fnm(j)(k:k).eq.'A') fnm(j)(k:k)='a'
               if(fnm(j)(k:k).eq.'B') fnm(j)(k:k)='b'
               if(fnm(j)(k:k).eq.'C') fnm(j)(k:k)='c'
               If(fnm(j)(k:k).eq.'D') fnm(j)(k:k)='d'
               if(fnm(j)(k:k).eq.'E') fnm(j)(k:k)='e'
               if(fnm(j)(k:k).eq.'F') fnm(j)(k:k)='f'
               if(fnm(j)(k:k).eq.'G') fnm(j)(k:k)='g'
               if(fnm(j)(k:k).eq.'H') fnm(j)(k:k)='h'
               if(fnm(j)(k:k).eq.'I') fnm(j)(k:k)='i'
               if(fnm(j)(k:k).eq.'J') fnm(j)(k:k)='j'
               if(fnm(j)(k:k).eq.'K') fnm(j)(k:k)='k'
               if(fnm(j)(k:k).eq.'L') fnm(j)(k:k)='l'
               if(fnm(j)(k:k).eq.'M') fnm(j)(k:k)='m'
               if(fnm(j)(k:k).eq.'N') fnm(j)(k:k)='n'
               if(fnm(j)(k:k).eq.'O') fnm(j)(k:k)='o'
               if(fnm(j)(k:k).eq.'P') fnm(j)(k:k)='p'
               if(fnm(j)(k:k).eq.'Q') fnm(j)(k:k)='q'
               if(fnm(j)(k:k).eq.'R') fnm(j)(k:k)='r'
               if(fnm(j)(k:k).eq.'S') fnm(j)(k:k)='s'
               if(fnm(j)(k:k).eq.'T') fnm(j)(k:k)='t'
               if(fnm(j)(k:k).eq.'U') fnm(j)(k:k)='u'
               if(fnm(j)(k:k).eq.'V') fnm(j)(k:k)='v'
               if(fnm(j)(k:k).eq.'W') fnm(j)(k:k)='w'
               if(fnm(j)(k:k).eq.'X') fnm(j)(k:k)='x'
               if(fnm(j)(k:k).eq.'Y') fnm(j)(k:k)='y'
               if(fnm(j)(k:k).eq.'Z') fnm(j)(k:k)='z'
 505          continue     
              fnm(j)(1:3)='pdb'
              fnm(j)(8:11)='.ent'
                indexhel=indexhel+1
 504          continue
  99      continue
               close (11)
                
c  WRITING HELIX INFORMATION ON OUTPUT FILES
c 
 125          do  506 k1=1,indexhel

             len = nre(k1)-nrb(k1)+1
              if (len.le.8) then
       write(6,35) fnm(k1),helnm(k1),rb(k1),cb(k1),nrb(k1),
     1 re(k1),ce(k1),nre(k1) 
       write(7,35) fnm(k1),helnm(k1),rb(k1),cb(k1),nrb(k1),
     1 re(k1),ce(k1),nre(k1) 
             write(6,50) 
             write(7,50) 
              else
          open(unit=4,file=fnm(k1))
       write(1,35) fnm(k1),helnm(k1),rb(k1),cb(k1),nrb(k1),
     1 re(k1),ce(k1),nre(k1) 
       write(2,35) fnm(k1),helnm(k1),rb(k1),cb(k1),nrb(k1),
     1 re(k1),ce(k1),nre(k1) 
       write(3,35) fnm(k1),helnm(k1),rb(k1),cb(k1),nrb(k1),
     1 re(k1),ce(k1),nre(k1) 
       write(8,35) fnm(k1),helnm(k1),rb(k1),cb(k1),nrb(k1),
     1 re(k1),ce(k1),nre(k1) 
       write(9,35) fnm(k1),helnm(k1),rb(k1),cb(k1),nrb(k1),
     1 re(k1),ce(k1),nre(k1) 
         if (fnm(k1).ne.fnm(k1-1)) write(10,67)
       write(10,35) fnm(k1),helnm(k1),rb(k1),cb(k1),nrb(k1),
     1 re(k1),ce(k1),nre(k1) 
c
c IS IT A PDB FILE? : FOURTH QUESTION
c

           if ((ans1.eq.'n').or.(ans1.eq.'N')) then
             if(fnm(k1).ne.fnm(k1-1)) then
 134         write(6,44) fnm(k1)
             write(7,44) fnm(k1)
            read(5,'(a1)') ans4
            write(7,'(1x,a1)') ans4
             if ((ans4.ne.'n').and.(ans4.ne.'N').and.(ans4.ne.'y')
     1   .and.(ans4.ne.'Y')) then
                 go to 134
               endif

             if((ans4.eq.'y').or.(ans4.eq.'y')) then
         write(6,48) prmfile,tabfile
         write(7,48) prmfile,tabfile
                     go to 135
                      endif

                if ((ans4.eq.'N').or.(ans4.eq.'n')) then
  139           rewind(4)
                write(6,45)
                write(7,45)
                read(5,'(a80)') formacor 
                write(7,'(1x,a79)') formacor 
              
                   j=1
 137             continue
            read(4,32,end=136) line
      read(line,formacor,err=139) atom,res,chain,nores,xx,yy,zz
          if(chain.eq.cb(k1)) then
         if ((atom.eq.'CA  ').or.(atom.eq.'  CA').or.
     -    (atom.eq.' CA ')) then
        if((nores.ge.nrb(k1)).and.(nores.le.nre(k1))) then
	      kk=nores-nrb(k1)+1
	      xca(kk)=xx
	      yca(kk)=yy
	      zca(kk)=zz
	      write(1,32) line
              endif
              endif
              endif
                   j=j+1
              go to 137
  136     continue
           write(6,48) prmfile,tabfile
           write(7,48) prmfile,tabfile
              go to 138 
            endif
            endif
            endif

c
C EXTRACTING C-ALPHA  ATOMS FROM PDB FILES
c
 135        read(4,32,end=999) line
       if ((line(1:4).eq.'ATOM').and.(line(14:15).eq.'CA')) then    
        read(line,33) atom,res,chain,nores,xx,yy,zz
         if(chain.eq.cb(k1)) then
        if((nores.ge.nrb(k1)).and.(nores.le.nre(k1))) then
	      kk=nores-nrb(k1)+1
	      xca(kk)=xx
	      yca(kk)=yy
	      zca(kk)=zz
	      write(1,32) line
              endif
              endif
              endif
               go to 135
  999     continue
c
c TAKING 4 CONSECUTIVE CA ATOMS
c
  138          kk3=kk-3
              pi=180.0/acos(-1.0)  
              do 508 l=1,kk3
              kount(l)=l
c
c VECTORS JOINING CA ATOMS 
c
              vec12(1)=xca(l+1)-xca(l)
              vec12(2)=yca(l+1)-yca(l)
              vec12(3)=zca(l+1)-zca(l)
c
              vec23(1)=xca(l+2)-xca(l+1)
              vec23(2)=yca(l+2)-yca(l+1)
              vec23(3)=zca(l+2)-zca(l+1)
c
              vec34(1)=xca(l+3)-xca(l+2)
              vec34(2)=yca(l+3)-yca(l+2)
              vec34(3)=zca(l+3)-zca(l+2)
c
c DIFFERENCE OF VECTORS JOINING CA ATOMS
c
              dv13(1)=vec12(1)-vec23(1)
              dv13(2)=vec12(2)-vec23(2)
              dv13(3)=vec12(3)-vec23(3)
c           
              dv24(1)=vec23(1)-vec34(1)  
              dv24(2)=vec23(2)-vec34(2)  
              dv24(3)=vec23(3)-vec34(3)  
c
c DIRECTION OF THE LOCAL HELIX AXIS
c            
      call cross(dv13,dv24,el,em,en)
              write(2,7) l,l+3, el,em,en
                uloc1(l)=el
                uloc2(l)=em
                uloc3(l)=en
c
                dmag=dot(dv13,dv13)
                dmag=sqrtf(dmag)
                emag=dot(dv24,dv24)
                emag=sqrtf(emag)
                costheta=dot(dv13,dv24)/(dmag*emag)
c
c    TWIST OF THE LOCAL HELIX, and thus number of  units per turn
c
                 twist(l)=acos(costheta)*pi
                 rnou(l)=360.0/twist(l)
c
c    RADIUS OF LOCAL HELIX CYLINDER
c
                   costheta1=1.0-costheta       
              radmag=sqrtf(dmag*emag)/(2.0*costheta1)
c
c    UNIT HEIGHT  OF THE LOCAL HELIX
c
      height(l)=vec23(1)*uloc1(l)+vec23(2)*uloc2(l)+vec23(3)*uloc3(l)    
           write(9,7)  l,l+3, twist(l),rnou(l),height(l)
c
c      LOCAL HELIX ORIGIN
c         
                   dv13(1)=dv13(1)/dmag
                   dv13(2)=dv13(2)/dmag
                   dv13(3)=dv13(3)/dmag
c                                    
                   rad(1)=radmag*dv13(1)
                   rad(2)=radmag*dv13(2)
                   rad(3)=radmag*dv13(3)
c          
                   ox(l)=xca(l+1)-rad(1)
                   oy(l)=yca(l+1)-rad(2)
                   oz(l)=zca(l+1)-rad(3)
c                   
                     dv24(1)=dv24(1)/emag
                     dv24(2)=dv24(2)/emag
                     dv24(3)=dv24(3)/emag
c
                     rad(1)=radmag*dv24(1)
                     rad(2)=radmag*dv24(2)
                     rad(3)=radmag*dv24(3)
c
                    ox(l+1)=xca(l+2)-rad(1)
                    oy(l+1)=yca(l+2)-rad(2)
                    oz(l+1)=zca(l+2)-rad(3)
 508                    continue
c
                        do  509 k=1,kk-2 
             if(abs(ox(k)).le.1.0E-04) ox(k) = 0.0000
             if(abs(oy(k)).le.1.0E-04) oy(k) = 0.0000
             if(abs(oz(k)).le.1.0E-04) oz(k) = 0.0000
              write(3,11) k,k+1,k+2,ox(k),oy(k),oz(k) 
 509                continue
                  write(3,38)
c
c LOCAL (k--k+3, k+3--k+6) BENDING ANGLES
c
                      kk6=kk3-3
                     do  510 k=1,kk6
         angle=uloc1(k)*uloc1(k+3)+uloc2(k)*uloc2(k+3)+
     -   uloc3(k)*uloc3(k+3)
               if(abs(angle-1.0).le.1.0E-06) angle = 1.0
                  angl(k)=acos(angle)*pi
                 write(8,13) k,k+3,K+3,k+6,angl(k)
 510                continue
c
c MEAN BENDING ANGLE FOR THE WHOLE HELIX
c
             call average(angl,kk6,avang,sdang,Advang)
                   write(8,27) avang,sdang,Advang
                   write(10,29) avang,sdang,Advang
c
c  MAXIMUM BENDING ANGLE IN THE WHOLE HELIX
c
                    call max(angl,kk6,anglmax)
                    write(8,26) anglmax
                    write(10,30) anglmax
c
c LOCAL BENDING ANGLES MATRIX
c
                 do  511 j=1,kk3
                 do  512 k=1,kk3
      proj=uloc1(j)*uloc1(k)+uloc2(j)*uloc2(k)+uloc3(j)*uloc3(k)
              if(abs(proj-1).le.1.0E-06) proj=1.0
              bendangl=acos(proj)*pi
              ibend=int(bendangl)
              bend2=bendangl-float(ibend)
c
               if(bend2.lt.0.5) then
              ibendangl(j,k)=int(bendangl)
                    else
              ibendangl(j,k)=int(bendangl+1.0)
                     endif
c
 512             continue   
 511             continue   
              write(2,16)
              write(2,25) 
               write(2,18) (kount(j),j=1,kk3)
               do  513 j=1,kk3
           write(2,17) j,(ibendangl(j,k),k=1,kk3)
 513            continue
c
c   AVERAGE HELICAL PARAMETERS
c
c    AVERAGE TWIST
c
             call average(twist,kk3,av,sd,Adv)
               write(9,19) av,sd,Adv
c
c    AVERAGE N
c
             call average(rnou,kk3,avn,sdn,Advn)
                write(9,20) avn,sdn,Advn
                write(10,36) avn,sdn,Advn
c
c    AVERAGE UNIT HEIGHT
c
             call average(height,kk3,avh,sdh,Advh)
               write(9,21) avh,sdh,Advh
               write(10,37) avh,sdh,Advh
c
c  FITTING LEAST SQUARES PLANE, LINE AND CIRCLE TO THE LOCAL HELIX ORIGINS
c
             call fit(ox,oy,oz,kk-2,radc,rmsdc,rmsdl,r2)
c
c  WRITING RESULTS IN TABULAR FORM
c
       call table(fnm(k1),cb(k1),rb(k1),nrb(k1),re(k1),nre(k1),
     - avn,sdn,avh,sdh,avang,sdang,anglmax,radc,rmsdc,rmsdl,r2,
     - geom)
               if(geom.eq.'L')  nlinear = nlinear+1
               if(geom.eq.'C')  ncurved = ncurved+1
               if(geom.eq.'K')  nkinked = nkinked+1
               if(geom.eq.'*')  nogeom = nogeom+1

               close (4)
                endif
 506           continue

                nhelices = nlinear+ncurved+nkinked+nogeom
          write(13,68) nlinear,ncurved,nkinked,nogeom,nhelices
c
C WRITING SOME TEXT TO TABLULAR FORM OF THE RESULTS. 
C
                 write (13,66)
                 write (13,54)
                 write (13,55)
                 write (13,56)
                 write (13,57)
                 write (13,58)
                 write (13,59)
                 write (13,60)
                 write (13,61)
                 write (13,62)
                 write (13,63)
                 write (13,64)
                 write(13,69) 
C
                 write(10,67) 
c

C FIFTH QUESTION : DO YOU WISH TO USE ANOTHER PDB FILE?
c
 156             if((ans1.eq.'y').or.(ans1.eq.'Y')) then
                 if((ans2.eq.'y').or.(ans2.eq.'Y')) then
 153           write(6,51)
               write(7,51)
                 read(5,'(a1)') ans5
                 write(7,'(1x,a1)') ans5
             if ((ans5.ne.'n').and.(ans5.ne.'N').and.(ans5.ne.'y')
     -   .and.(ans5.ne.'Y')) then
                 goto 153
               endif
             
                if ((ans5.eq.'y').or.(ans5.eq.'Y')) then
                  go to 155 
                 endif

                if ((ans5.eq.'n').or.(ans5.eq.'N')) then
                  go to 154 
                 endif
                 endif
                 endif
      
 154            continue
                stop 
                end

c
      function dot (x,y)                                                
c
c COMPUTES AND RETURNS THE DOT PRODUCT OF  THE VECTORS X AND Y                        
c
      dimension x(3),y(3)                                           
      dot = 0.0                                                   
      do 1000 i=1,3                                             
      dot = dot + x(i) * y(i)                                  
 1000 continue                                               
      return                                                
      end                                                  
c
c    SUBROUTINE TO CALCULATE CROSS PRODUCTS OF VECTORS A and B.
c
       subroutine cross(a,b,el,em,en)
         dimension a(3), b(3), c(3)
       c(1) = a(2) * b(3) - a(3) * b(2)
       c(2) = a(3) * b(1) - a(1) * b(3)
       c(3) = a(1) * b(2) - a(2) * b(1)
       fac1 = c(1) * c(1) + c(2) * c(2) + c(3) * c(3)
       fac1 = sqrtf(fac1)
       if(fac1.le.0.00001) go to 100
       fac1 = 1.0/fac1
 100    continue
       el = c(1) * fac1
       em = c(2) * fac1
       en = c(3) * fac1
      return
      end

c
c    SUBROUTINE TO CALCULATE AVERAGES
c
                 subroutine average(b,npts,av,sd,Adv)                                  
c CALCULATES AVERAGE, STANDARD DEVIATION AND MEAN ABSOLUTE DEVIATION
                dimension b(100)
                sum = 0.0 
                do 514 j = 1, npts
                sum = sum + b(j)
 514            continue
                av= sum/float(npts)
                sum1=0.0
                sum2=0.0
                 do 515 j=1,npts
                sum1=sum1+(b(j)-av)*(b(j)-av)
                sum2=sum2+abs(b(j)-av)
 515               continue
                 sum1=sum1/float(npts-1)
                 sd=sqrtf(sum1)
                 Adv=sum2/float(npts)
                  return
                  end
c
c   SUBROUTINE TO FIT PLANE,CIRCLE AND LINE TO LOCAL HELIX ORIGINS

              subroutine fit(ox,oy,oz,np,radc,rmsdc,rmsdl,r2)
c
              dimension ox(100),oy(100),oz(100)
              dimension rx(100),ry(100),rz(100)
              dimension rmp(100),rmc(100),rml(100)
              dimension matp(3,3),matc(3,3),matl(2,2),rotmt(3,3)
              dimension pmat(3,3),cmat(3,3),lmat(2,2)
              dimension v(3),ap(3),ac(3),al(2)
c             real matp,matc,matl,lmat
       double precision matp,pmat,matc,cmat,matl,lmat
c
   22  format(1x,'R.M.S. DEVIATION FROM BEST PLANE',8x,f8.3,2x,'ANG.')
   23  format(1x,'L, M AND N OF THE BEST PLANE',11x,3(1x,f8.3))
   24  format(1x,'R.M.S. DEVIATION FROM BEST CIRCLE',7x,f8.3,2x,'ANG.')
   25  format(1x,'RADIUS OF THE BEST CIRCLE',13x,f10.3,2x,'ANG.')
   26  format(1x,'CENTRE OF THE BEST CIRCLE',13x,2(f10.3,1x))
   27  format(1x,'R.M.S. DEVIATION FROM BEST LINE',9x,f8.3,2x,'ANG.')
   28  format(1x,'SLOPE OF THE BEST LINE',18x,f8.3)
   29  format(1x,'INTERCEPT OF THE BEST LINE',14x,f8.3)
   30  format(1x,'SQUARE OF LINEAR CORRELATION COEFFICIENT',f8.3)
   31  format(1x,'FIT TO THE CIRCLE IS NOT GOOD')
   32  format(1x,'FIT TO THE LINE IS NOT GOOD')
   33  format(1x,'FIT TO THE PLANE IS NOT GOOD')
   34  format(1x,'PERPENDICULAR DISTANCE OF BEST PLANE',4x,f8.3,
     1 2x,'ANG.')
c
c FITTING LEAST SQUARES PLANE TO LOCAL HELIX ORIGINS
c
                 x2=0.0000000 
                 y2=0.0000000
                 z2=0.0000000   
                 xy=0.0000000 
                 xz=0.0000000 
                 yz=0.0000000 
                  x=0.0000000  
                  y=0.0000000 
                  z=0.0000000 
c
                do 516 j=1,np
                x2=ox(j)*ox(j)+x2 
                y2=oy(j)*oy(j)+y2 
                z2=oz(j)*oz(j)+z2 
                xy=ox(j)*oy(j)+xy 
                xz=ox(j)*oz(j)+xz 
                yz=oy(j)*oz(j)+yz 
                x=ox(j)+x 
                y=oy(j)+y 
                z=oz(j)+z 
 516            continue
c
                   do 100 l=1,3
                   do 101 m=1,3
                  matp(l,m)=0.000000
 101            continue
 100            continue
                matp(1,1)=x2+matp(1,1)
                matp(1,2)=xy+matp(1,2)
                matp(1,3)=xz+matp(1,3)
                matp(2,1)=xy+matp(2,1)
                matp(2,2)=y2+matp(2,2)
                matp(2,3)=yz+matp(2,3)
                matp(3,1)=xz+matp(3,1)
                matp(3,2)=yz+matp(3,2)
                matp(3,3)=z2+matp(3,3)
 
               call matinv(matp,pmat)
c
                bp1=x 
                bp2=y 
                bp3=z 
c
                 do 517 j=1,3
             ap(j)=pmat(j,1)*bp1+pmat(j,2)*bp2+pmat(j,3)*bp3
 517               continue
c
                sump=0.0 
                do  518 j=1,np
            rmp(j)=ap(1)*ox(j)+ap(2)*oy(j)+ap(3)*oz(j)-1 
                sump=sump+rmp(j)*rmp(j) 
 518            continue
c
                if(sump.lt.0.0) then
                write(3,33)
                else
                rmsdp=(sqrtf(sump/np))
                write(3,22) rmsdp
                endif
c
c   CONVERTING TO THE NORMAL FORM OF THE PLANE
c deno is denominator sqrt(ap(1)*ap(1)+ap(2)*ap(2)+ap(3)*ap(3))
c pp is the distance of plane from XY plane
c
           deno=sqrtf(ap(1)*ap(1)+ap(2)*ap(2)+ap(3)*ap(3)) 
                     pp=1.0/deno 
                     pl=ap(1)*pp  
                     pm=ap(2)*pp 
                     pn=ap(3)*pp 
             write(3,23) pl,pm,pn
             write(3,34) pp
c
c TO REORIENT THE POINTS SO THAT THE BEST FIT PLANE COINCIDES WITH 
c (X-Y) PLANE
c The l,m,n of the normal to (X-Y) plane are (0,0,1)
c Angle between the normal to best fit plane & the (X-Y) plane.
c
c     
c VECTOR NORMAL TO THE L,M,N OF THE BEST FIT & (X-Y) PLANE
c
                  do 102 j=1,3
                  v(j) = 0.0000000
 102            continue
                  v(1)=pm+v(1) 
                  v(2)=-pl+v(2) 
                  v(3)=0.0+v(3) 
             vmag=sqrtf(v(1)*v(1)+v(2)*v(2)+v(3)*v(3)) 
                  v(1)=v(1)/vmag 
                  v(2)=v(2)/vmag 
                  v(3)=v(3)/vmag 
c  THE ROTATION MATRIX REQUIRED TO MAKE l,m,n as 0,0,1
c 
                         costheta=pn 
                	 costheta1=1.0-pn 
                         theta=acos(pn) 
	                 sintheta=sin(theta) 
	                 a1=v(1)*sintheta 
	                 a2=v(2)*sintheta 
	                 a3=v(3)*sintheta 
	                 b1=v(2)*v(3)*costheta1 
	                 b2=v(3)*v(1)*costheta1 
	                 b3=v(1)*v(2)*costheta1 
                 do 103 j=1,3
                 do 104 k=1,3
               rotmt(j,k)=0.0000000
 104           continue
 103           continue
            rotmt(1,1)=costheta+v(1)*v(1)*costheta1+rotmt(1,1) 
            rotmt(2,2)=costheta+v(2)*v(2)*costheta1+rotmt(2,2) 
	    rotmt(3,3)=costheta+v(3)*v(3)*costheta1+rotmt(3,3)           
     	    rotmt(1,2)=b3-a3+rotmt(1,2) 
	    rotmt(2,1)=b3+a3+rotmt(2,1) 
	    rotmt(3,1)=b2-a2+rotmt(3,1) 
	    rotmt(1,3)=b2+a2+rotmt(1,3) 
	    rotmt(2,3)=b1-a1+rotmt(2,3) 
	    rotmt(3,2)=b1+a1+rotmt(3,2) 
c
c ROTATING THE POINTS SO THAT THEY LIE IN (X-Y) PLANE
c
                 do 105 j=1,np
                rx(j) = 0.00000000
                ry(j) = 0.00000000
                rz(j) = 0.00000000
 105             continue
                do  519 j=1,np
      rx(j)=rotmt(1,1)*ox(j)+rotmt(1,2)*oy(j)+rotmt(1,3)*oz(j)+rx(j)  
      ry(j)=rotmt(2,1)*ox(j)+rotmt(2,2)*oy(j)+rotmt(2,3)*oz(j)+ry(j) 
      rz(j)=rotmt(3,1)*ox(j)+rotmt(3,2)*oy(j)+rotmt(3,3)*oz(j)+rz(j) 
              if(abs(rx(j)).le.1.0E-04) rx(j)=0.0000
              if(abs(ry(j)).le.1.0E-04) ry(j)=0.0000
              if(abs(rz(j)).le.1.0E-04) rz(j)=0.0000
  519           continue
c
c TRANSLATING THE POINTS SO THAT THEY LIE IN X-Y PLANE OF THE FORM
c ax+by=0
c
                do 520 j=1,np
                rx(j)=rx(j)  
                ry(j)=ry(j) 
                rz(j)=rz(j)-pp 
 520            continue
c
c FITTING CIRCLE TO THESE REORIENTED POINTS
c
c  Equation used is (x+a)^2+(y+b)^2=r^2
c center of circle (-a,-b) ,c=r^2-a^2-b^2
c delta=x^2+y^2+2ax+2by-c
c circle finds out radius of curvature and coordinates of center
c from a set of given points using LEAST SQUARE FIT method.
c
                x2=0.0000000 
                y2=0.0000000 
                x3=0.0000000 
                y3=0.0000000 
                xy2=0.0000000 
                x2y=0.0000000  
                xy=0.0000000 
                x=0.0000000 
                y=0.0000000 
c
                do 521 j=1,np
                x3=rx(j)*rx(j)*rx(j)+x3 
                y3=ry(j)*ry(j)*ry(j)+y3 
                x2y=rx(j)*rx(j)*ry(j)+x2y 
                xy2=rx(j)*ry(j)*ry(j)+xy2 
                x2=rx(j)*rx(j)+x2 
                y2=ry(j)*ry(j)+y2 
                xy=rx(j)*ry(j)+xy 
                x=rx(j)+x 
                y=ry(j)+y 
 521            continue
c
                  do 106 j=1,3
                  do 107 k=1,3
                matc(j,k)=0.00000000
 107             continue
 106             continue
                matc(1,1)=x2+matc(1,1)
                matc(1,2)=xy+matc(1,2) 
                matc(1,3)=-x+matc(1,3) 
                matc(2,1)=xy+matc(2,1)
                matc(2,2)=y2+matc(2,2) 
                matc(2,3)=-y+matc(2,3) 
                matc(3,1)=x+matc(3,1) 
                matc(3,2)=y+matc(3,2) 
                matc(3,3)=-np+matc(3,3) 
                  
c
                 call matinv(matc,cmat)
c
                bc1=-x3-xy2 
                bc2=-x2y-y3 
                bc3=-x2-y2 
        
                 do 108 l=1,3
                 ac(l) = 0.00000000
 108             continue
                 do 522 j=1,3
                ac(j)=cmat(j,1)*bc1+cmat(j,2)*bc2+cmat(j,3)*bc3+ac(j)
 522             continue

                ac(1)=ac(1)/2.0
                ac(2)=ac(2)/2.0
                ac(3)=ac(3)
                radcsq=ac(1)*ac(1)+ac(2)*ac(2)+ac(3)
		if (radcsq.gt.0.0) then 
                radc=sqrtf(radcsq)
                write(3,25) radc
                write(10,25) radc
                write(3,26) -ac(1), -ac(2)
                   else
		   radc = 0.0
                    write(10,31)
		    endif
c
                  sumc=0.00000000 
c
c dev=sqrt((x+a)^2+(y+b)^2)-r
c rms dev=sqrt(dev^2/n)
c
                   do 523 j=1,np
      rem=(rx(j)+ac(1))*(rx(j)+ac(1))+(ry(j)+ac(2))*(ry(j)+ac(2)) 
                 rmc(j)=sqrtf(rem)-radc 
                sumc=sumc+rmc(j)*rmc(j) 
 523               continue
                  if(sumc.lt.0.0) then
                    write(3,31)
                    write(10,31)
                    else
                    rmsdc=sqrtf(sumc/np)
                    write(3,24) rmsdc
                    write(10,24) rmsdc
                   endif
c
c FITTING LINE TO THE REORIENTED  POINTS
c EQUATION OF LINE Y=MX+C
c                 
c
                     do 109 l=1,2
                     do 110 m=1,2
                    matl(l,m) = 0.00000000
 110              continue
 109              continue
                 
                    matl(1,1)=x2+matl(1,1) 
                    matl(1,2)=x+matl(1,2) 
                    matl(2,1)=x+matl(2,1) 
                    matl(2,2)=np+matl(2,2) 
 
                   call matinv2(matl,lmat)
c
                    bl1=xy 
                    bl2=y 
c
                    do 111 k=1,2
                    al(k) = 0.0000000
 111                continue
                    do  524 j=1,2
                    al(j)=lmat(j,1)*bl1+lmat(j,2)*bl2+al(j)
  524                continue
c
                     write(3,28) al(1)
                     write(3,29) al(2)
c
                    suml=0.00000000 
c
                    do  525 j=1,np
                    rml(j)=ry(j)-al(1)*rx(j)-al(2) 
                    suml=suml+rml(j)*rml(j) 
 525               continue
                    if(suml.lt.0.0) then
                    write(3,32)
                    write(10,32)
                    else
                    rmsdl=(sqrtf(suml/np))
                    write(3,27) rmsdl
                    write(10,27) rmsdl
                    endif
c
c LINEAR CORRELATION COEFFICIENT
c r=(n*(xy)-x*y)/sqrt((n*x2-x*x)*(n*y2-y*y))
c
                    rnum=np*xy-x*y 
                    rdeno=(np*x2-x*x)*(np*y2-y*y) 
                    r=rnum/sqrtf(rdeno) 
                    r2=r*r
                    write(3,30) r2
                    write(10,30) r2
                    return
                    end

C
       FUNCTION SQRTF(X)
  1    FORMAT(2X,'SQUARE ROOT OF A NEGATIVE QUANTITY?',E16.8)
       IF(X.LT.0.0) GO TO 101
          SQRTF = SQRT(X)
          RETURN
  101   CONTINUE
       WRITE(*,1) X
          SQRTF = 0.0
          RETURN
          END
c
c SUBROUTINE TO COMPUTE INVERSE OF A 3 X 3 MATRIX
c
                subroutine matinv(h,s)
c CALCULATES INVERSE OF A 3x3 MATRIX
           dimension h(3,3),s(3,3)
       double precision c11,c12,c13,c21,c22,c23,c31,c32,c33
       double precision deth1,deth2,deth3,deth,s,h
                do 112 j = 1,3
                do 113 k = 1,3
                s(j,k) = 0.0000000
 113      continue
 112      continue
              c11=h(2,2)*h(3,3)-h(3,2)*h(2,3)
              c12=h(2,1)*h(3,3)-h(3,1)*h(2,3)
              c13=h(2,1)*h(3,2)-h(3,1)*h(2,2)
              c21=h(1,2)*h(3,3)-h(1,3)*h(3,2)
              c22=h(1,1)*h(3,3)-h(1,3)*h(3,1)
              c23=h(1,1)*h(3,2)-h(3,1)*h(1,2)
              c31=h(1,2)*h(2,3)-h(1,3)*h(2,2)
              c32=h(1,1)*h(2,3)-h(1,3)*h(2,1)
              c33=h(1,1)*h(2,2)-h(1,2)*h(2,1) 
            deth1=h(1,1)*c11                            
            deth2=-h(1,2)*c12                          
            deth3=h(1,3)*c13                          
            deth=deth1+deth2+deth3
	       if(deth.ne.0.0) then
            s(1,1)=c11/deth+s(1,1)
            s(1,2)=-c21/deth+s(1,2)
            s(1,3)=c31/deth+s(1,3)
            s(2,1)=-c12/deth+s(2,1)
            s(2,2)=c22/deth+s(2,2)
            s(2,3)=-c32/deth+s(2,3)
            s(3,1)=c13/deth+s(3,1)
            s(3,2)=-c23/deth+s(3,2)
            s(3,3)=c33/deth+s(3,3)
	      else
	      write (*,*) 'MATRIX IS SINGULAR'
	      endif
            return
            end  
c
c SUBROUTINE TO COMPUTE INVERSE OF A 2 X 2 MATRIX
c
             subroutine matinv2(p,q)
C  INVERSE OF A 2x2 MATRIX
           dimension p(2,2),q(2,2)
       double precision c11,c12,c21,c22
       double precision detp1,detp2,detp,p,q
               do 114 j = 1,2
               do 115 k = 1,2
             q(j,k) = 0.00000000
 115         continue
 114         continue
              c11=p(2,2)
              c12=p(2,1)
              c21=p(1,2)
              c22=p(1,1)
            detp1=p(1,1)*c11                            
            detp2=-p(1,2)*c12                          
            detp=detp1+detp2
             if(detp.ne.0.0) then
            q(1,1)=c11/detp+q(1,1)
            q(1,2)=-c21/detp+q(1,2)
            q(2,1)=-c12/detp+q(2,1)
            q(2,2)=c22/detp+q(2,2)
                else
            write(*,*) 'MATRIX IS SINGULAR'
               endif
            return
            end  
c
c SUBROUTINE TO FIND LARGEST NUMBER

               subroutine max(x,no,xmax)
c              FINDS MAXIMUM NO.
                  dimension x(100)
                     xmax=x(1)
                    do 525 j=1,no
                   if (x(j).gt.xmax) xmax=x(j)
 525                continue
                     return
                     end

      subroutine table(code,cbg,rbg,nrbg,red,nred,turnn,sdn,
     - height,sdh,bav,sdb,bam,rad,rmsc,rmsl,rsq,geom)

c  SUBROUTINE TO WRITE HELIX ANALYSIS RESULTS IN TABULAR FORM.
c

		       character*11 code
                       character*3 rbg,red
		       character*1 cbg,rb1,re1
		       character*1 geom


    9  format(1x,a11,1x,a1,1x,i3,a1,' - ',i3,a1,1x,f4.2,1x,
     -     f4.2,1x,f4.1,1x,f4.1,2x,i4,2x,2(f5.2,1x),f4.2,3x,a1)
   10    format(16x,'Std. Dev.',2x,f4.2,1x,f4.1,1x,f4.1)

               if(rbg.eq.'ALA') rb1='A'
               if(rbg.eq.'CYS') rb1='C'
               if(rbg.eq.'ASP') rb1='D'
               if(rbg.eq.'GLU') rb1='E'
               if(rbg.eq.'PHE') rb1='F'
               if(rbg.eq.'GLY') rb1='G'
               if(rbg.eq.'HIS') rb1='H'
               if(rbg.eq.'ILE') rb1='I'
               if(rbg.eq.'LYS') rb1='K'
               if(rbg.eq.'LEU') rb1='L'
               if(rbg.eq.'MET') rb1='M'
               if(rbg.eq.'ASN') rb1='N'
               if(rbg.eq.'PRO') rb1='P'
               if(rbg.eq.'GLN') rb1='Q'
               if(rbg.eq.'ARG') rb1='R'
               if(rbg.eq.'SER') rb1='S'
               if(rbg.eq.'THR') rb1='T'
               if(rbg.eq.'VAL') rb1='V'
               if(rbg.eq.'TRP') rb1='W'
               if(rbg.eq.'TYR') rb1='Y'

               if(red.eq.'ALA') re1='A'
               if(red.eq.'CYS') re1='C'
               if(red.eq.'ASP') re1='D'
               if(red.eq.'GLU') re1='E'
               if(red.eq.'PHE') re1='F'
               if(red.eq.'GLY') re1='G'
               if(red.eq.'HIS') re1='H'
               if(red.eq.'ILE') re1='I'
               if(red.eq.'LYS') re1='K'
               if(red.eq.'LEU') re1='L'
               if(red.eq.'MET') re1='M'
               if(red.eq.'ASN') re1='N'
               if(red.eq.'PRO') re1='P'
               if(red.eq.'GLN') re1='Q'
               if(red.eq.'ARG') re1='R'
               if(red.eq.'SER') re1='S'
               if(red.eq.'THR') re1='T'
               if(red.eq.'VAL') re1='V'
               if(red.eq.'TRP') re1='W'
               if(red.eq.'TYR') re1='Y'

	       if(bam.ge.20.0) geom='K'

	       if (bam.lt.20.0) then
                if ((rmsc.gt.1.0).and.(rmsl.gt.1.0)) then
                         geom = '*'
                          else 
	             ratio=rmsl/rmsc
	       if (ratio.gt.1.0) then 
                      geom='C' 
                       endif
	       if (ratio.le.0.7) then
               if (rsq.ge.0.8) then 
                     geom='L'
                      else
                      geom='*'
                     endif
                     endif

               if ((ratio.gt.0.7).and.(ratio.le.1.0)) then
                if(rsq.le.0.5) then 
                   geom='C'
                       endif
                if(rsq.ge.0.8) then
                     geom='L'
                        endif
                if ((rsq.gt.0.5).and.(rsq.lt.0.8)) then 
                    geom = '*'
                    endif
	       endif
                          endif
                          endif
                if ((rad-int(rad)).lt.0.5) then
                   irad=int(rad)
                   else
                   irad=int(rad+0.5)
                   endif
            write(13,9) code,cbg,nrbg,rb1,nred,re1,turnn,height,bav,
     -      bam,irad,rmsc, rmsl,rsq,geom
            write(13,10)  sdn, sdh, sdb

		return
		end

