Logo Search packages:      
Sourcecode: raster3d version File versions  Download package

balls.f

 
      PROGRAM BALLS
*
*         Program to set up input for RENDER
*
* This program takes a Brookhaven PDB-format file of atomic coordinates
* concatenated with a "colours" file that specifies the colouring and spherical
* radius of atoms according to atom type, residue type, or any other criterion,
* and produces a file that can be used directly as input to RENDER.
*
* Usage:    cat colors.pdb protein.pdb | balls [-h] | render
*
* auxiliary files:      setup.matrix      optional file containing 3x3 view matrix
*                 setup.angles      optional file defining view matrix
*                             in terms of 3 rotation angles
*
*     I/O units for colour/co-ordinate input, specs output, user output
*     
      INCLUDE 'VERSION.incl'
*
      INTEGER INPUT, OUTPUT, NOISE
      PARAMETER (INPUT=5, OUTPUT=6, NOISE=0)
      PARAMETER (MAXCOL=5000, MAXATM=300000)
      REAL RGB(3,MAXCOL),RADIUS(MAXCOL)
      REAL SPAM(7,MAXATM)
      CHARACTER*24 MASK(MAXCOL),TEST
      CHARACTER*80 ATOM(MAXATM),CARD
      LOGICAL MATCH
      logical           hflag
      character*80      flags
c
C     Ethan Merritt     Oct 1988
C     Modified to read in 3x3 view matrix (e.g. from CCP FRODO view command)
C     from file setup.matrix.  Matrix is applied before
C     finding translation, center, and scale.  Afterwards the input matrix
C     to RENDER is therefore the identity matrix.
C                                                     
      common /matrix/ matrix, postrn, coords
      real*4            matrix(3,3), postrn(3), coords(3)
C
C     Default to CPK colors and VDW radii
      character*60 defcol(7)
      data defcol /
     & 'COLOUR#######C################   0.625   0.625   0.625  1.70',
     & 'COLOUR#######N################   0.125   0.125   1.000  1.60',
     & 'COLOUR#######O################   0.750   0.050   0.050  1.50',
     & 'COLOUR#######S################   1.000   1.000   0.025  1.85',
     & 'COLOUR#######H################   1.000   1.000   1.000  1.20',
     & 'COLOUR#######P################   0.400   1.000   0.400  1.80',
     & 'COLOUR########################   1.000   0.000   1.000  2.00'
     &            /
c
    3 format(a,a)
c
      do i=1,3
      do j=1,3
          matrix(i,j)=0.
      enddo
      matrix(i,i)=1.
      enddo
      call view_matrix
c
c     Ethan Merritt Apr 1992
c     -h option suppresses header records in output
c
      hflag = .false.
      narg = iargc()
      do i = 1, narg
          call getarg( i, flags )
          if (flags(1:2) .eq. '-h') hflag = .true.
      end do
c
      if (.not. hflag) then
      WRITE(OUTPUT,'(A,A)') 'balls - Raster3D ', VERSION
      WRITE(OUTPUT,'(A)') '80  64    tiles in x,y'
      WRITE(OUTPUT,'(A)') ' 8   8    pixels (x,y) per tile'
      WRITE(OUTPUT,'(A)') '4         anti-aliasing 3x3 -> 2x2 pixels'
      WRITE(OUTPUT,'(A)') '0 0 0     black background'
      WRITE(OUTPUT,'(A)') 'T         yes, I LIKE shadows!'
      WRITE(OUTPUT,'(A)') '25        Phong power'
      WRITE(OUTPUT,'(A)') '0.15      secondary light contribution'
      WRITE(OUTPUT,'(A)') '0.05      ambient light contribution'
      WRITE(OUTPUT,'(A)') '0.25      specular reflection component'
      WRITE(OUTPUT,'(A)') '4.0       eye position'
      WRITE(OUTPUT,'(A)') '1 1 1     main light source position'
      end if
c
      ASPECT = 1280./1024.
      NCOL = 0
      NATM = 0
10    CONTINUE
        READ(INPUT,'(A80)',END=50) CARD
        IF (CARD(1:4).EQ.'COLO') THEN
          NCOL = NCOL + 1
          IF (NCOL.GT.MAXCOL) THEN
            WRITE(NOISE,*) 'Colour table overflow.  Increase ',
     &                     'MAXCOL and recompile.'
            STOP 10
          ENDIF
          READ(CARD,'(6X,A24,3F8.3,F6.2)') MASK(NCOL),
     &          (RGB(I,NCOL),I=1,3),RADIUS(NCOL)
        ELSEIF (CARD(1:4).EQ.'ATOM'.OR.CARD(1:4).EQ.'HETA') THEN
          NATM = NATM + 1
          IF (NATM.GT.MAXATM) THEN
            WRITE(NOISE,*) 'Atom array overflow.  Increase ',
     &                     'MAXATM and recompile.'
            STOP 20
          ENDIF
          ATOM(NATM) = CARD
        ELSEIF (CARD(1:3).EQ.'END') THEN
          GO TO 50
        ENDIF
        GO TO 10
*     Come here when EOF or 'END' record is reached
50    CONTINUE
      IF (NATM.EQ.0) THEN
        WRITE(NOISE,*) 'No atoms in input.'
        STOP 30
      ENDIF
*     Load default colors after any that were read in
      IF (NCOL.LT.MAXCOL-8) THEN
        DO i = 1,7
          NCOL = NCOL + 1
          READ(defcol(i),'(6X,A24,3F8.3,F6.2)') MASK(NCOL),
     &          (RGB(J,NCOL),J=1,3), RADIUS(NCOL)
        ENDDO
      ENDIF
*
      IF (NCOL.EQ.0) THEN
        WRITE(NOISE,*) 'No colours in input.'
        STOP 40
      ENDIF
      XMAX = -1E20
      XMIN =  1E20
      YMAX = -1E20
      YMIN =  1E20
      ZMAX = -1E20
      ZMIN =  1E20
      DO 100 IATM=1,NATM
        CARD = ATOM(IATM)
        TEST = CARD(7:30)
        DO 80 ICOL=1,NCOL
          IF (MATCH(TEST,MASK(ICOL))) THEN
c           READ(CARD,'(30X,3F8.3)') X,Y,Z
c EAM Oct88
            READ(CARD,'(30X,3F8.3)') coords
            x = coords(1)*matrix(1,1) + coords(2)*matrix(2,1) 
     1        + coords(3)*matrix(3,1)
            y = coords(1)*matrix(1,2) + coords(2)*matrix(2,2) 
     1        + coords(3)*matrix(3,2)
            z = coords(1)*matrix(1,3) + coords(2)*matrix(2,3)
     1        + coords(3)*matrix(3,3)
c EAM Oct88
            RAD = RADIUS(ICOL)
            SPAM(1,IATM) = X
            SPAM(2,IATM) = Y
            SPAM(3,IATM) = Z
            SPAM(4,IATM) = RAD
            SPAM(5,IATM) = RGB(1,ICOL)
            SPAM(6,IATM) = RGB(2,ICOL)
            SPAM(7,IATM) = RGB(3,ICOL)
            XMAX = MAX(XMAX,X+RAD)
            XMIN = MIN(XMIN,X-RAD)
            YMAX = MAX(YMAX,Y+RAD)
            YMIN = MIN(YMIN,Y-RAD)
            ZMAX = MAX(ZMAX,Z+RAD)
            ZMIN = MIN(ZMIN,Z-RAD)
            GO TO 100
          ENDIF
80      CONTINUE
        WRITE(NOISE,*) 'No colour table mask matches this atom:'
        WRITE(NOISE,*) ATOM(IATM)
        STOP 90
100   CONTINUE
      XMID = (XMAX+XMIN)/2.
      YMID = (YMAX+YMIN)/2.
      ZMID = (ZMAX+ZMIN)/2.
      TX = -XMID
      TY = -YMID
      TZ = -ZMID
      IF (ASPECT.GE.1.) THEN
*       The X direction is wider than the Y
        XROOM = ASPECT
        YROOM = 1.
        ZROOM = 2.
      ELSE
        XROOM = 1.
        YROOM = ASPECT
        ZROOM = 2.
      ENDIF
      XSPAN = XMAX-XMIN
      YSPAN = YMAX-YMIN
      ZSPAN = ZMAX-ZMIN
      SCALE = MAX(XSPAN/XROOM,YSPAN/YROOM,ZSPAN/ZROOM)
*     Leave a little extra room as a border:
      SCALE = SCALE / 0.90
c
      if (.not. hflag) then
      WRITE(OUTPUT,120) TX,TY,TZ,SCALE
120   FORMAT('1 0 0 0   input co-ordinate + radius transformation'/
     &       '0 1 0 0'/
     &       '0 0 1 0'/
     &       4F10.3)
      WRITE(OUTPUT,'(A)') '3         mixed object types'
      WRITE(OUTPUT,'(A)') '*'
      WRITE(OUTPUT,'(A)') '*'
      WRITE(OUTPUT,'(A)') '*'
      end if
c
      DO 150 IATM=1,NATM
      SPAM(1,IATM) = SPAM(1,IATM) + postrn(1)
      SPAM(2,IATM) = SPAM(2,IATM) + postrn(2)
      SPAM(3,IATM) = SPAM(3,IATM) + postrn(3)
        WRITE(OUTPUT,140) (SPAM(I,IATM),I=1,7)
140   FORMAT(1H2,/,7(1X,F8.3))
150   CONTINUE
      write (noise,'(/)')
      write (noise,156) 'X  min max:', XMIN, XMAX
      write (noise,156) 'Y  min max:', YMIN, YMAX
      write (noise,156) 'Z  min max:', ZMIN, ZMAX
      write (noise,156) '     scale:', SCALE
  156 format(1x,a,3f8.2)
      END
      LOGICAL FUNCTION MATCH (SUBJ, MASK)
      CHARACTER*24 SUBJ,MASK
      MATCH = .FALSE.
      DO 10 I = 1, 24
        IF (SUBJ(I:I).NE.MASK(I:I) .AND. MASK(I:I).NE.'#') RETURN
10    CONTINUE
      MATCH = .TRUE.
      RETURN
      END

      subroutine view_matrix
c
      common /matrix/ matrix, postrn, coords
      real*4            matrix(3,3), postrn(3), coords(3)
c
      real*4            phiX, phiY, phiZ
      integer    noise
      parameter (noise = 0)
      parameter (R2D = 180./3.1415927)

      postrn(1) = 0.
      postrn(2) = 0.
      postrn(3) = 0.

      open (unit=3, file='setup.matrix', status='OLD', err=100)
            read (3,*) ((matrix(i,j),i=1,3),j=1,3)
            write (noise,'(1x,3f9.5)') ((matrix(i,j),i=1,3),j=1,3)
            close (3)

            det = matrix(1,1) * matrix(2,2) * matrix(3,3)
     1          + matrix(1,2) * matrix(2,3) * matrix(3,1)
     2          + matrix(2,1) * matrix(3,2) * matrix(1,3)
     3          - matrix(1,3) * matrix(2,2) * matrix(3,1)
     4          - matrix(1,2) * matrix(2,1) * matrix(3,3)
     5          - matrix(1,1) * matrix(2,3) * matrix(3,2)
            write (noise,'(''       determinant ='',f8.3)') det

            phiX = atan2( -matrix(3,2), matrix(3,3) )
            phiY = atan2(  matrix(3,1), matrix(3,3) / cos(phiX) )
            phiZ = atan2( -matrix(2,1), matrix(1,1) )
            write (noise,3) ' View Angles from matrix',' '
            write (noise,2) phiZ*R2D, phiY*R2D, phiX*R2D
            return
  100 continue

      open (unit=3, file='setup.angles', status='OLD', err=200)
            read (3,*) phiZ, phiY, phiX
            read (3,*,err=101,end=101) postrn(1), postrn(2), postrn(3)
  101       continue
            close (3)
            write (noise,2) phiZ, phiY, phiX
            write (noise,4) postrn(1), postrn(2), postrn(3)
            cx = cos(phiX/R2D)
            sx = sin(phiX/R2D)
            cy = cos(phiY/R2D)
            sy = sin(phiY/R2D)
            cz = cos(phiZ/R2D)
            sz = sin(phiZ/R2D)
            matrix(1,1) = cz*cy
            matrix(1,2) = sz*cx + cz*sy*sx
            matrix(1,3) = sz*sx - cz*sy*cx
            matrix(2,1) = -sz*cy
            matrix(2,2) = cz*cx - sx*sy*sz
            matrix(2,3) = cz*sx + sz*sy*cx
            matrix(3,1) = sy
            matrix(3,2) = -sx*cy
            matrix(3,3) = cx*cy
            write (noise,3) ' View Matrix from angles',' '
            write (noise,'(1x,3f9.5)') ((matrix(i,j),i=1,3),j=1,3)
            return
  200 continue

    2       format(1x,'   phiZ =',f8.2,'   phiY =',f8.2,'   phiX =',f8.2)
    3 format(/a,a)
    4       format(1x,'  tranX =',f8.2,'  tranY =',f8.2,'  tranZ =',f8.2)

      write (noise,*) ' No view matrix or angles provided'
      return
      end

Generated by  Doxygen 1.6.0   Back to index