C
C	PLT1   (OMNIGRAPH VERSION)
C
C	MODIFIED BY R.S. ZHOU AT BOSTON COLLEGE, 09/14/89.
C	TO LINK TO GRAFIT'S PLOTTING LIBRARY.
C
C***********************************************************************
C
	COMMON /STRO/ ISTRO, STROX, STROY, SSTROX, SSTROY, MSEC, MSEC2, 
	1              IINZ, NZMX
C
	COMMON /DAT/ TI(10), XMIN(3), XMAX(3), DX(3), XCHAR(3), ACHAR(3), 
	1  IDIM(2), JDIM(2), NEXT, XSHX, XSHY, XSHZ, YSHY, YSHZ, ZSHZ, 
	2  XMPMIN, XMPMAX, YMPMIN, YMPMAX, ISCALE, NSDONE, NIP, ZMAP, A, 
     	3  RMIN, RMAX, RINC, NSEC, ANG, RC, RMINN, RMAXN, RINCN
C
	COMMON /IOSET/ IO2, IO3, IO4, IO5, IO6, IO7, IO8
C
	COMMON /MOD/ AM(3, 3, 36), AT(3, 36), XMIN3(3), XMAX3(3), 
	1 DZ3, DZ3N, DZ2, DZ2N, DMX(4), CHAR(4, 20), NSYM, PWIDE(4), ICHR(4), 
	2 HT(4), HT2(4), SCALE, RAD(4, 20), RDEL(4), IMASK(36, 4), 
	3 NSET(4, 20), ATMFIL(80, 4), VRAD(4, 20), VTYP(4, 20), IFVDW(4), 
	4 IFVDW0, IRESOL, IPEN(20,4)
	COMMON/OFFSET/XOFF,XFAC,YOFF,YFAC,OFFX,OFFY
C
	DIMENSION ID(3), A(3), ANG(3)
	DIMENSION XMIN2(3)
	INTEGER PWIDE, CHAR, VTYP
	LOGICAL*1 ATMFIL
C
     	INTEGER   CRT
     	LOGICAL   YES
     	CHARACTER PLTFIL*32, DESFIL*32
     	CHARACTER COMM*16, PLOTID*32
     	DIMENSION HPAR(20)
     	DATA HPAR/ 4.,1.,1.,1.,0,0,1.,1.,1.,0,5.,1.,1.,0,0,0,0,0,0,0 /
     	DATA IO2, IO3, IO4, IO6 / 2, 3, 4, 16 /
C
C  PARSE THE ACTIVATING DCL COMMAND LINE...
C
C     	CRT = IF_DCL_PUT( 'PLOT1Z' )
C     	IF( CRT.LT.0 ) GOTO 888
     	CALL IFINIT_SESSION
     	IF( .NOT.IF_DCL_HARDCOPY() ) THEN
     	    CRT = IFINIT_ITM(0)
     	    CALL IFANSI_CLEAR(0)
     	ENDIF
     	XORG = 0.0
     	YORG = 0.0
     	OFFX = 0.0
     	OFFY = 0.0
     	IRESOL = 2
     	PLOTID = ' '
C
C  GET THE PLOT FILE...
C
     	CALL IF_DCL_FILE( 1, PLTFIL )
 4	IF( PLTFIL.NE.' ' ) THEN
     	    CALL IF_DCL_LAYOUT( DESFIL )
     	    IF( DESFIL.EQ.' ' ) THEN
     		DESFIL = PLTFIL
     		CALL IFMISC_NAMEXT( DESFIL, '.DES' )
     	    ENDIF
     	ELSE
 5     	    WRITE(6,5000)
5000	    FORMAT(' Enter plot file [.PLT] ',$)
     	    READ(5,5001,END=999) PLTFIL
5001	    FORMAT(A)
     	    IF( PLTFIL.EQ.' ' ) GOTO 5
     	    CALL IFMISC_DEFEXT( PLTFIL, '.PLT' )
     	    INQUIRE( FILE=PLTFIL, EXIST=YES )
     	    IF( .NOT.YES ) THEN
     		WRITE(6,6)
 6		FORMAT(' oops: can not locate the plot file.')
     		GOTO 5
     	    ENDIF
 8	    DESFIL = PLTFIL
     	    CALL IFMISC_NAMEXT( DESFIL, '.DES' )
5101	    OFFX = 0.0
     	    OFFY = 0.0
     	    WRITE(6,5100)
5100	    FORMAT(' Plot offsets on screen (x,y) [0.0,0.0] ',$)
     	    READ(5,5102,ERR=5101,END=999) STROX,STROY
     	    OFFX = OFFX + STROX
     	    OFFY = OFFY + STROY
5102	    FORMAT(2F10.0)
5201	    WRITE(6,5200)
5200	    FORMAT(' Space-filling line resolution',
     +' (Low/Medium/High) [L] ',$)
    	    READ(5,5001,ERR=5201,END=999) YES
     	    IRESOL = 0
     	    IF( YES.EQ.'M' .OR. YES.EQ.'m' ) IRESOL = 1
     	    IF( YES.EQ.'H' .OR. YES.EQ.'h' ) IRESOL = 2
     	ENDIF
C
C  READ THE PLOT FILE....
C
	CALL TRSTUP( PLTFIL,
     +  XMIN2, MDIR, NDZ, RSCX, RSCY, XFRAM, YFRAM, IFNFLD, XFNFLD)

	ZMAPS  = XMIN(3)
	SC     = 1.0 / SCALE
	SCX    = 1.0 / ( SCALE * RSCX)
	SCY    = 1.0 / ( SCALE * RSCY)
	SSTROX = DX( 3) * STROX
	SSTROY = DX( 3) * STROY
	STML0  = YMPMIN
	STMT0  = XMPMIN - 8.0 / SCALE
	IF( ISTRO.NE.0 ) STML0 = STML0 - FLOAT( MSEC - 1) * 0.5 * SSTROY
	IF( ISTRO.NE.0 ) STMT0 = STMT0 - FLOAT( MSEC - 1) * 0.5 * SSTROX
	ASTROX = STROX
	ASTROY = STROY
C
C  START THE FRAME LOOP...
C
     	III = 0
 1234	III = III + 1
     	IF( III.LT.1 .OR. III.GT.NSEC ) THEN
     	    IF( IF_DCL_HARDCOPY() ) GOTO 4321
     	    CALL IFANSI_VIDEO(70)
     	    WRITE(6,1235)
 1235	    FORMAT(' oops: Frame number out of bound.')
     	    IF( III.LT.1 ) III = 1
     	    IF( III.GT.NSEC ) III = NSEC
     	    GOTO 444
     	ENDIF
     	IF( IF_DCL_HARDCOPY() ) WRITE(6,1236) III
 1236	FORMAT(' Drawing frame #',I2,'...',$)
	ZMAP = ZMAPS + MSEC2*DX(3)*(III-1)
	ZMAPO=ZMAP*ZSHZ
	ZMAPF=ZMAP/A(3)
	WRITE(IO6, 230) XCHAR(3), ZMAPO, ZMAPF
 230	FORMAT('0BEGINNING SECTION ', A1, ' = ', F7.3,'A (',F6.3,')')
     	STROSN = 0
	DO 260 JSTRO=1,2
	IF(ISTRO.EQ.0.AND.JSTRO.EQ.2)GO TO 260
C
C  INITIALIZE THE PLOTTER...
C
     	IF( JSTRO.EQ.1 ) THEN
     	    CALL IFINIT_LINE(3)
     	    CALL IFINIT_PLOT( XORG, YORG, PLTFIL )
     	    IF( IFCTRL_C() ) GOTO 333
C
C  OPEN A GRAPHICS SEGMENT TO SAVE THE PLOT...
C
     	    IF( CRT.EQ.2 ) THEN
     		CALL IFSEG_PIVOT( 0.0, 0.0 )
     		CALL IFSEG_OPEN( 15 )
     	    ENDIF
     	    CALL IFINIT_DIF(0)
     	    CALL IFSPEN(7)
     	    COMM = 'CURRENT_FRAME:'
     	    ENCODE( 2, 240, COMM(15:16) ) III
  240	    FORMAT(I2)
     	    CALL IFTEXT( 0.1-XORG, 7.0-YORG, COMM, 0.1, 0.1, 0.0 )
     	    CALL IFINIT_DIF(1)
     	    CALL IFSPEN(1)
     	ENDIF
     	IF( ISTRO.NE.0 ) STROSN = ( -1.0 )**( JSTRO + 1 )
     	STROX = ASTROX * STROSN
     	STROY = ASTROY * STROSN
     	IF( IFCTRL_C() ) GOTO 333
     	CALL IFSPEN(1)
C     	CALL FRAME( SC, RSCX, RSCY, XFRAM, YFRAM, IFNFLD, XFNFLD)
C     	IF( IFCTRL_C() ) GOTO 333
C     	CALL PLOT1Z_OFFSET( STMT0, SCX, STML0, SCY)
C     	CALL PLOT1Z_FACE( ISTRO, JSTRO)
     	IF( JSTRO.EQ.1 ) CALL PLOT1Z_TITLE( XORG, YORG, TI )  !R.S. ZHOU
     	IF( IFCTRL_C() ) GOTO 333
     	CALL PLOT1Z_OFFSET( STMT0, SCX, STML0, SCY)
     	CALL MODEL( STROX, STROY)
     	IF( IFCTRL_C() ) GOTO 333
     	IF( IFVDW0.EQ.3.OR.IFVDW0.EQ.5 ) THEN   !R.S. ZHOU
     	    CALL SRFDRW( SSTROX, SSTROY, STML0, STMT0, MSEC, 
     +                   STROSN, SCX, SCY, STROX, STROY)
     	END IF
     	IF( IFVDW0.EQ.2.OR.IFVDW0.EQ.5 ) THEN   !R.S. ZHOU
     	    CALL EDGDRW( SSTROX, SSTROY, STML0, STMT0, MSEC, 
     +                   STROSN, SCX, SCY, STROX, STROY)
     	END IF
	NSDONE = NSDONE + 1
C	IF(IFNFLD.NE.0) THEN
C	   XXXX = XFNFLD * RSCX
C	   CALL PLOT1Z_QPLOT( XXXX, 0.0, -3)
C	END IF
	IF(ISTRO.GT.0.AND.JSTRO.EQ.1) CALL IFPLOT_VECT(4.5,0.0,-3)
 260	CONTINUE
	WRITE(IO6, 231)
 231	FORMAT('0SECTION COMPLETED.')
C
C  CLOSE THE FRAME...
C
     	IF( CRT.EQ.2 ) CALL IFSEG_OPEN( 0 ) 	!CLOSE THE SEGMENT
    	IF( DESFIL.NE.' ' ) CALL LAYOUT_AUX(0,XORG,YORG,DESFIL)
  333	CALL IFINIT_PLOT_END
     	IF( IF_DCL_HARDCOPY() ) GOTO 666
     	IF( CRT.EQ.1 ) THEN
     	    CALL IFANSI_CURSOR(1,1)
     	    CALL IFANSI_CLEAR(6)
     	ENDIF
C
C  PAUSE FOR INTERACTIVE INPUT....
C
  444	CALL IFANSI_VIDEO(47)
     	WRITE(6,1000)
 1000	FORMAT(' PLOT1Z> Please press a desired function key:')
  445	WRITE(6,1001)
 1001	FORMAT(' PLOT1Z> ?>Help 0>Exit 1>Hardcopy 2>Layout',
     +' 3>Newfile 4>Down_ 5>Up_frame ...')
     	CALL IFANSI_CURSOR(0,8)
     	CALL IFQIOW_KEY( KEY )

     	IF( KEY.EQ.26 .OR. KEY.EQ.48 .OR. KEY.EQ.69 ) GOTO 4321
     	IF( KEY.EQ.49 .OR. KEY.EQ.72 ) THEN		!hardcopy
     	    CALL IFCOPY( HPAR(1), PLOTID )
     	    CALL IFINIT_LINE( -3 )
     	ELSE IF( KEY.EQ.-80 ) THEN			!hardocpy
     	    CALL IFANSI_CURSOR(0,14)
     	    CALL IFANSI_LEDS(1)
     	    CALL IFFEED( HPAR(1), PLOTID, ' ' )
     	    CALL IFANSI_LEDS(0)
     	ELSE IF( KEY.EQ.50 .OR. KEY.EQ.76 ) THEN	!layout
     	    IF( DESFIL.EQ.' ' ) THEN
     		DESFIL = PLTFIL
     		CALL IFMISC_NAMEXT( DESFIL, '.DES' )
     	    ENDIF
     	    CALL LAYOUT_AUX( 2, XORG, YORG, DESFIL )
     	ELSE IF( KEY.EQ.51 .OR. KEY.EQ.-82 ) THEN	!newfile
     	    PLTFIL = ' '
    	    GOTO 4
     	ELSE IF( KEY.EQ.52 .OR. KEY.EQ.68 ) THEN	!down_frame
    	    GOTO 1234
     	ELSE IF( KEY.EQ.53 .OR. KEY.EQ.85 ) THEN	!up_frame
     	    III = III-2
     	    GOTO 1234
     	ELSE IF( KEY.EQ.54 ) THEN		!redraw the same frame
     	    III = III-1
     	    GOTO 1234
     	ELSE IF( KEY.EQ.13 .OR. KEY.EQ.32
     +  .OR. KEY.EQ.63 .OR. KEY.EQ.-81 ) THEN !help
     	    CALL IFMISC_HLP( 'PLOT1Z', '$ MENU', -2, 0 )
     	ELSE IF( KEY.EQ.36 .OR. KEY.EQ.-83 ) THEN	!dcl prompt
     	    CALL IFMISC_SPAWN( '$', 0 )
     	    CALL IFINIT_LINE( -3 )
     	ELSE IF( KEY.EQ.64 ) THEN
554	    WRITE(6,555)
555	    FORMAT(' Which frame: ',$)
     	    READ(5,556,ERR=554,END=444) KEY
556	    FORMAT(I)
     	    IF( KEY.LT.1 .OR. KEY.GT.NSEC ) THEN
     		CALL IFANSI_VIDEO(70)
     		WRITE(6,557)
557		FORMAT('_oops: non-existing frame.')
     		GOTO 554
     	    ENDIF
     	    III = KEY - 1
     	    GOTO 1234
     	ELSE
     	    CALL IFMISC_OOPS( 5, 0, 0, ' ' )
    	ENDIF
     	GOTO 445
C
C  HARDCOPY THE CURRENT FRAME...
C
 666	CALL IFCOPY( HPAR(1), PLOTID )
     	GOTO 1234			!GET THE NEXT FRAME
4321	CLOSE(UNIT=IO6)
999	CALL IFINIT_EXIT
888	CALL EXIT
	END

C
C==========================================================================
C
	SUBROUTINE QUIT( NSDONE)
	COMMON /IOSET/ IO2, IO3, IO4, IO5, IO6, IO7, IO8
	WRITE( IO6, 100) NSDONE
 100	FORMAT( '0HALTING AFTER COMPLETING ', I2, ' SECTIONS.')
	IF( NSDONE .EQ. 0) STOP 'ABNORMAL TERMINATION TO PLOT1Z.'
	CLOSE( UNIT = IO6)
	STOP 'PLOT1Z NORMAL TERMINATION.'
	END

C
C==========================================================================
C
	SUBROUTINE FRAME( SC, RSCX, RSCY, XFRAM, YFRAM, IFNFLD, XFNFLD)
C
	CALL IFSPEN( 1)
	CALL PLOT1Z_QPLOT( 0.0, 8.0, -3)
	CALL PLOT1Z_QPLOT( 0.0, 0.0,  +3)
	XX  = XFRAM * RSCX
	XXX = XFNFLD * RSCX
	YY  = YFRAM * RSCY
	CALL PLOT1Z_QPLOT( 0.0, YY,  +2)
	CALL PLOT1Z_QPLOT( XX, YY,  +2)
	CALL PLOT1Z_QPLOT( XX, 0.0,  +2)
	CALL PLOT1Z_QPLOT( 0.0, 0.0,  +2)
	IF( IFNFLD.NE.0.0) CALL PLOT1Z_QPLOT( XXX, 0.0,  +2)
	CALL IFSPEN( 2)
	Y1 = -8.0
	Y2 = YY
	X1 = 0.0
	NX = ( XX + 20.0 ) / 2.0
	IF( NX.GT.XXX ) NX = XXX / 2.0
	NY = 24
	IF( Y2 + 8.0.GT.260.0 ) GO TO 8
	NY = ( 260.1 - Y2 - 8.0 ) / 2.0
	IF( NY.GT.24) NY = 8
 8	CONTINUE
	NX = ( XXX - XX ) / 2.0
	IF( NX.GT.20 ) NX = 20
	NY = YY / 4.0
	X1 = XX
	CALL IFSPEN( 1)
	RETURN
	END

     	SUBROUTINE PLOT1Z_TITLE( XORG, YORG, TI )
C
C  DRAW PLOT TITLE...
C
     	DIMENSION TI(10), TII(10)
     	CHARACTER TITLE*40
     	EQUIVALENCE (TITLE(1:1),TII(1))
     	DO I=1, 10
     	   TII(I) = TI(I)
     	ENDDO
     	IF( TITLE.NE.' ' ) CALL IFTEXT_CENT( -XORG, 10.0-XORG,
     +  6.8-YORG, 6.8-YORG, TITLE, 0.12, 0.12 )
     	RETURN
     	END

C
C==========================================================================
C
	SUBROUTINE PLOT1Z_FACE(ISTRO,JSTRO)
C
C      Labels the section, etc.
C
	COMMON /DAT/ TI(10), XMIN(3), XMAX(3), DX(3), XCHAR(3), ACHAR(3), 
	1  IDIM(2), JDIM(2), NEXT, XSHX, XSHY, XSHZ, YSHY, YSHZ, ZSHZ, 
	2  XMPMIN, XMPMAX, YMPMIN, YMPMAX, ISCALE, NSDONE, NIP, ZMAP, A, 
     	3  RMIN, RMAX, RINC, NSEC, ANG, RC, RMINN, RMAXN, RINCN
	COMMON /IOSET/ IO2,IO3,IO4,IO5,IO6,IO7,IO8
	REAL*4 IX1,IY1,IX2,IY2,IX3,IY3,IX4,IY4
	REAL*4 IX5,IY5,IX6,IY6,IX7,IY7,IX8,IY8
     	CHARACTER TITLE*40, SECT*40
	LOGICAL*1 STRCHR(5),STROC
	DATA STRCHR/' ','T','B','L','R'/
     	EQUIVALENCE (TITLE(1:1),TI(1))
C
C      Title the section
C
     	IF( TITLE.NE.' ' ) THEN
     	    CALL IFSPEN(3)
     	    CALL IFTEXT( 0.0, -0.6, TITLE, 0.12, 0.12, 0.0 )
     	ENDIF
	K=1+JSTRO+2*(ISTRO-1)
	IF(ISTRO.EQ.0)K=1
	STROC=STRCHR(K)
	ISC=IABS(ISCALE)
     	ENCODE(40,1,SECT) ISC,XCHAR(3),ZMAP,ACHAR(2),ACHAR(1),STROC
  1	FORMAT('G=',I2,'A ',A1,'=',F6.2,'A, ',A1,
     +  '-HORI,',A1,'-VERT, S.P.-',A1,'.')
     	CALL IFSPEN(1)
	CALL IFTEXT( 0.0, -0.2, SECT, 0.1, 0.1, 0.0 )
C
C      outline the plotting area
C
	IX1 = XMPMIN
	IX2 = XMPMAX
	IY1 = YMPMIN
	IY2 = YMPMAX
	CALL PLOT1Z_LINEV(IX1, IY1, IX1, IY2)
	CALL PLOT1Z_LINEV(IX2, IY1, IX2, IY2)
	CALL PLOT1Z_LINEV(IX1, IY1, IX2, IY1)
	CALL PLOT1Z_LINEV(IX1, IY2, IX2, IY2)
C
C      place scaled grid
C
	IF(ISCALE.EQ.0)GO TO 6
	SC=IABS(ISCALE)
	SC1=SC/10.0
	NTKXI=IFIX(XMPMIN/SC + 1000.0)- 1000
	NTKYI=IFIX(YMPMIN/SC + 1000.0)- 1000
	NTKX=(XMPMAX-SC*NTKXI)/SC
	NTKY=(YMPMAX-SC*NTKYI)/SC
	IF(ISCALE.LT.0)GO TO 11
	NTKXI1=IFIX(XMPMIN/SC1 + 1000.0)- 1000
	NTKYI1=IFIX(YMPMIN/SC1 + 1000.0)- 1000
	NTKX1=(XMPMAX-SC1*NTKXI1)/SC1
	NTKY1=(YMPMAX-SC1*NTKYI1)/SC1
	DO 4 I = 1, NTKX
	   X = I + NTKXI
	   X = SC*X
	   DO 8 J=1,NTKY1
	      Y=J+NTKYI1
	      Y=SC1*Y
	      CALL PLOT1Z_QPLOT(X,Y,+13)
	      CALL PLOT1Z_QPLOT(X,Y,+12)
 8	   CONTINUE
 4	CONTINUE
	DO 5 I = 1, NTKY
	   Y = I + NTKYI
	   Y = SC*Y
	   DO 10 J=1,NTKX1
	      X=J+NTKXI1
	      X=X*SC1
	      CALL PLOT1Z_QPLOT(X,Y,+13)
	      CALL PLOT1Z_QPLOT(X,Y,+12)
 10	   CONTINUE
 5	CONTINUE
	GO TO 6
 11	CONTINUE
	DO 12 I=1,NTKX
	   X=I+NTKXI
	   X=X*SC
	   DO 12 J=1,NTKY
	      Y=SC*(J+NTKYI)
	      CALL PLOT1Z_MARK(X,Y,0.08,9)
 12	   CONTINUE
 6	CONTINUE
	IX5 = X
	RETURN
	END

C
C==========================================================================
C
	SUBROUTINE PLOT1Z_LINEV( X1, Y1, X2, Y2)
C
C	SUBROUTINE TO DRAW A LINE FROM X1, Y1 TO X2, Y2.
C
	CALL PLOT1Z_QPLOT( X1, Y1,  + 13)
	CALL PLOT1Z_QPLOT( X2, Y2,  + 12)
	RETURN
	END

C
C==========================================================================
C
	SUBROUTINE MODEL1( NDIR, XYZO, MR)
C
	COMMON /DAT/ TI(10), XMIN(3), XMAX(3), DX(3), XCHAR(3), ACHAR(3), 
	1  IDIM(2), JDIM(2), NEXT, XSHX, XSHY, XSHZ, YSHY, YSHZ, ZSHZ, 
	2  XMPMIN, XMPMAX, YMPMIN, YMPMAX, ISCALE, NSDONE, NIP, ZMAP, ABC, 
	3  RMIN, RMAX, RINC, NSEC, ANG, RC
C
	COMMON /IOSET/ IO2, IO3, IO4, IO5, IO6, IO7, IO8
C
	COMMON /STRO/ ISTRO, STROX, STROY, SSTROX, SSTROY, MSEC, MSEC2, IINZ
C
	COMMON /MOD/ AM(3, 3, 36), AT(3, 36), XMIN3(3), XMAX3(3), 
	1 DZ3, DZ3N, DZ2, DZ2N, DMX(4), CHAR(4, 20), NSYM, PWIDE(4), ICHR(4), 
	2 HT(4), HT2(4), SCALE, RAD(4, 20), RDEL(4), IMASK(36, 4), 
	3 NSET(4, 20), ATMFIL(80, 4), VRAD(4, 20), VTYP(4, 20), IFVDW(4), 
	4 IFVDW0, IRESOL, IPEN(20,4)
C
	REAL*4    MR(3, 3), XYZO(3), A(3, 3), T(3), B(3, 3), ABC(3)
	INTEGER   PWIDE, CHAR, VTYP
	LOGICAL*1 ATMFIL
C
	K1 = 1
	K2 = 2
	K3 = 3
C	CALL IROT( K1, K2, K3, NDIR)
	READ( IO2, 1)NSYM
 1	FORMAT( I5)
	IF( NSYM.EQ.0)RETURN
	DO 2 I = 1, 4
	   DMX( I)   = 0.0
	   PWIDE( I) = 0
	   ICHR( I)  = 0
	   HT( I)    = 0.0
	   READ( IO2, 3) IQ, ( ATMFIL( K, I), K = 1, IQ)
	   IF( IQ.LE.3)GO TO 2
 3	   FORMAT( Q, 80A1)
 	   READ( IO2, 11) DMX( I), PWIDE( I), ICHR( I), HT( I), 
	1                 ( NSET( I, K), K = 1, 20)
 11	   FORMAT( F5.2, 2I5, F5.2, 20I5)
	   IF( NSET( I, 2).EQ.0)NSET( I, 2) = 32000
	   WRITE( IO6, 111)I, ( ATMFIL( K, I), K = 1, IQ)
	   WRITE( IO6, 113) DMX( I), PWIDE( I), ICHR( I), 
	1    HT( I), ( NSET( I, K), K = 1, 20)
 111	   FORMAT( '1*****   INPUT PARAMETERS FOR MODEL PORTION', I2, /, 
	1    '  MODEL FROM FILE ', 80A1)
 113	   FORMAT( '   INITIAL MAX BOND DIST. = ', F6.2, /, 
	2     '   PEN WIDTH  = ', I3, /, 
	3     '   CHARACTER TYPE  = ', I3, /, 
	4     '   CHARACTER HEIGHT  = ', F4.2, 'MM.', 
	5    10( /, '    ACCEPT ATOMS FROM #', I5, ' TO #', I5))
	   HT2( I) = HT( I) / 2.0
	   IF( ICHR( I).EQ.1)GO TO 32
	   HT2( I) = 0.0
	   READ( IO2, 33) ( CHAR( I, J), J = 1, 20)
 33	   FORMAT( 20I4)
	   IF( ICHR( I).NE.0 ) ICHR( I) = -1
	   GO TO 34
 32	   READ( IO2, 35) ( CHAR( I, J), J = 1, 20)
 35	   FORMAT( 20( 3X, A1))
 34	   CONTINUE
	   READ( IO2, 37) RDEL( I), ( RAD( I, J), J = 1, 20)
 37	   FORMAT( 21F5.2)
	   READ( IO2, 13) IFVDW( I), 
	1  ( VTYP( I, J), VRAD( I, J), J = 1, 20), (IPEN(J,I),J=1,20)
 13	   FORMAT( I5, 10( I3, F5.1), /, 5X, 10( I3, F5.1),/,5X,20I3)
	   DO 12 J = 1, 20
	      IF( VRAD( I, J).LT.0.0.OR.VRAD( I, J).GT.5.0 )
	1              VRAD( I, J) = 0.0
C	      IF( VTYP( I, J).LT.0.OR.VTYP( I, J).GT.99) VTYP( I, J) = 0
	      IF( VTYP( I, J).LT.1.OR.VTYP( I, J).GT.9) VTYP( I, J) = 1
 12	   CONTINUE
	   IF( ICHR( I))42, 46, 46
 42	   WRITE( IO6, 43) RDEL( I), ( J, CHAR( I, J), RAD( I, J), 
 	1          VRAD( I, J), VTYP( I, J), J = 1, 20)
 43	   FORMAT( '0 INCREMENT IN BOND SEARCH  = ', F5.3, /, 
	1   '0 TYPE  CHARACTER  BOND RADIUS VDW RADIUS TYPE', /, 
	2    ( ' ', I5, 2X, I4, 8X, F6.3, F10.3, I3))
	   GO TO 48
 46	   WRITE( IO6, 47) RDEL( I), ( J, CHAR( I, J), RAD( I, J), 
 	1        VRAD( I, J), VTYP( I, J), J = 1, 20)
 47	   FORMAT( '0 INCREMENT IN BOND SEARCH  = ', F5.3, /, 
	1   '0 TYPE  CHARACTER  BOND RADIUS VDW RADIUS TYPE', /, 
	2    ( ' ', I5, 5X, A1, 8X, F6.3, F10.3, I3))
 48	   CONTINUE
	   READ( IO2, 39)( IMASK( J, I), J = 1, NSYM)
 39	   FORMAT( 36I2)
	   WRITE( IO6, 49)( IMASK( J, I), J = 1, NSYM)
 49	   FORMAT( '0  SYMMETRY OPERATORS TO BE USED', 36I2)
 2	CONTINUE
	IFVDW0 = 0
	DO I = 1, 4
	   IF( IFVDW( I).EQ.2 )THEN
	      IF( IFVDW0.LT.2 )THEN
	         IFVDW0 = 2
	      ELSE IF( IFVDW0.EQ.3 )THEN
	         IFVDW0 = 5
	      END IF
C	   ELSE IF( IFVDW( I).EQ.2 )THEN
	   ELSE IF( IFVDW( I).EQ.3 )THEN	!R.S. ZHOU
	      IF( IFVDW0.LT.2 )THEN
	         IFVDW0 = 3
	      ELSE IF( IFVDW0.EQ.2 )THEN
	         IFVDW0 = 5
	      ELSE IF( IFVDW0.EQ.3 )THEN
	         IFVDW0 = 3
	      END IF
	   ELSE IF( IFVDW( I).EQ.5 )THEN
	      IFVDW0 = 5
	   END IF
	END DO
	DMXMX = 3.5
	IF( 1.5 * ABS( DMX( 1)).GT.DMXMX ) DMXMX = ABS( DMX( 1)) * 1.5
	IF( 1.5 * ABS( DMX( 2)).GT.DMXMX ) DMXMX = ABS( DMX( 2)) * 1.5
	IF( 1.5 * ABS( DMX( 3)).GT.DMXMX ) DMXMX = ABS( DMX( 3)) * 1.5
	IF( 1.5 * ABS( DMX( 4)).GT.DMXMX ) DMXMX = ABS( DMX( 4)) * 1.5
	XMIN3( 1) = XMPMIN - DMXMX
	XMAX3( 1) = XMPMAX + DMXMX
	XMIN3( 2) = YMPMIN - DMXMX
	XMAX3( 2) = YMPMAX + DMXMX
	DZ3  = DX( 3) * 0.6666667 * FLOAT( MSEC)
	DZ3N = -DZ3
	DZ2  = DX( 3) * 0.5 * FLOAT( MSEC)
	DZ2N = -DZ2
	IF( NSYM.LE.36)GO TO 22
	WRITE( IO6, 213)NSYM
 213	FORMAT( ' **** LIMIT OF 36 SYMMETRY OPERATORS EXCEEDED, NSYM = ',
 	1        I3)
	STOP
 22	WRITE( IO6, 21) NSYM
 21	FORMAT( '1', I5, ' SYMMETRY EQUIVALENT ATOMS TO BE GENERATED')
 	DO 10 I = 1, NSYM
	   READ( IO2, 17)A( K1, 1), A( K1, 2), A( K1, 3), A( K2, 1), 
	1       A( K2, 2), A( K2, 3), A( K3, 1), A( K3, 2), A( K3, 3), 
	2       T( K1), T( K2), T( K3)
	   WRITE( IO6, 23)( A( K, 1), A( K, 2), A( K, 3), T( K), K = 1, 3)
 23	   FORMAT( '0', 3( /, '  ( ', 3F10.5, ' )  + ', F10.5))
 17	   FORMAT( 12F6.3)
	   DO I1 = 1, 3
	      AT( I1, I) = 0.0
	      DO I2 = 1, 3
	         AM( I1, I2, I) = 0.0
	         DO K = 1, 3
	            AM( I1, I2, I) = AM( I1, I2, I) 
	1                          + MR( I1, K) * A( K, I2)
		 END DO
	      END DO
	      DO K = 1, 3
	         AT( I1, I) = AT( I1, I) 
	1                   + MR( I1, K) * ( T( K) - XYZO( K) )
	      END DO
	   END DO
 10	CONTINUE
	WRITE( IO6, 27)
 27	FORMAT( '1')
	RETURN
	END

C
C==========================================================================
C
	SUBROUTINE TRSTUP( PLTFIL,
     +  XMIN2, MDIR, NDZ, RSCX, RSCY, XFRAM, YFRAM, IFNFLD, XFNFLD)
C
	COMMON /DAT/ TI(10), XMIN(3), XMAX(3), DX(3), XCHAR(3), ACHAR(3), 
	1  IDIM(2), JDIM(2), NEXT, XSHX, XSHY, XSHZ, YSHY, YSHZ, ZSHZ, 
	2  XMPMIN, XMPMAX, YMPMIN, YMPMAX, ISCALE, NSDONE, NIP, ZMAP, ABC, 
	3  RMIN, RMAX, RINC, NSEC, ANG, RC, RMINN, RMAXN, RINCN
C
	COMMON /STRO/ ISTRO, STROX, STROY, SSTROX, SSTROY, MSEC, MSEC2, 
	1  IINZ, NZMX
C
	COMMON /IOSET/ IO2, IO3, IO4, IO5, IO6, IO7, IO8
C
	COMMON /MOD/ AM(3, 3, 36), AT(3, 36), XMIN3(3), XMAX3(3), 
	1 DZ3, DZ3N, DZ2, DZ2N, DMX(4), CHR(4, 20), NSYM, PWIDE(4), ICHR(4), 
	2 HT(4), HT2(4), SCALE, RAD(4, 20), RDEL(4), IMASK(36, 4), 
	3 NSET(4, 20), ATMFIL(80, 4), VRAD(4, 20), VTYP(4, 20), IFVDW(4), 
	4 IFVDW0, IRESOL, IPEN(20,4)
C
	DIMENSION ID(3), ABC(3), ANG(3)
	DIMENSION XMIN2(3), XMAX2(3)
	DIMENSION XMIN1(3), XMAX1(3)
	REAL*4  ML, MT
	REAL*4  XYZO(3), MO(3, 3)
	INTEGER PWIDE, CHR, VTYP
	DIMENSION CHAR(6)
	LOGICAL*1 ATMFIL
     	CHARACTER LOGFIL*32, PLTFIL*(*)
	DATA CHAR/'X', 'Y', 'Z', 'A', 'B', 'C'/
C
	LEXT = 0
	DGTRD = 0.017453293
	DO I = 1, 3
	   J = I + 3
	   XCHAR( I) = CHAR( I)
	   ACHAR( I) = CHAR( J)
	END DO
	NIP = 12
     	LOGFIL = PLTFIL
     	CALL IFMISC_NAMEXT( LOGFIL, '.LOG' )
     	IF( IF_DCL_HARDCOPY() ) WRITE(6,1) LOGFIL
 1	FORMAT(/' *** Run log file: ',A32,/)
	OPEN(UNIT=IO6,FILE=LOGFIL,FORM='FORMATTED',STATUS='UNKNOWN')
	OPEN(UNIT=IO2,FILE=PLTFIL,TYPE='OLD',FORM='FORMATTED')
	RC = 3.141592653 / 180.0
	READ ( IO2, 100) TI
 100	FORMAT( 10A4)
	WRITE( IO6, 101) TI
 101	FORMAT( ' ', 10A4)
	READ( IO2, 102) ABC, ANG
 102	FORMAT( 3F10.4, 3F10.6)
	WRITE( IO6, 103) ABC, ANG
 103	FORMAT( '0A = ', F10.4, 10X, 'B = ', F10.4, 10X, 'C = ', F10.4/
	1  '0ALPHA = ', F10.6, 10X, 'BETA = ', F10.6, 10X, 'GAMMA = ', 
	2  F10.6)
	READ( IO2, 1103)( XMIN1( I), XMAX1( I), I = 1, 3), DZ, NSEC
 1103	FORMAT( 6F10.4, F10.4, I5)
	DX( 1) = 0.0
	DX( 2) = 0.0
	DX( 3) = DZ
	WRITE( IO6, 1104)( XMIN1( I), XMAX1( I), I = 1, 3), DZ, NSEC
 1104	FORMAT( '0 VIEW LIMIT AND SECTION THICKNESS', 
	1 /, '   ', 6F10.3, F10.3, /, '   MAX NUMBER OF SECTIONS  = ', I5)
	SN1 = SIN( ANG( 1) * RC)
	CS1 = COS( ANG( 1) * RC)
	SN2 = SIN( ANG( 2) * RC)
	CS2 = COS( ANG( 2) * RC)
	SN3 = SIN( ANG( 3) * RC)
	CS3 = COS( ANG( 3) * RC)
	XSHX = 1.0
	XSHY = CS3
	XSHZ = CS2
	YSHY = SN3
	YSHZ = ( CS1 - CS2 * CS3) / SN2
	ZSHZ = ( 1.0 - CS1**2 - CS2**2 - CS3**2 + 2.0 * CS1 * CS2 * CS3) 
	ZSHZ = SQRT( ZSHZ / SN3**2 )
	WRITE( IO6, 202) XSHX, XSHY, XSHZ, YSHY, YSHZ, ZSHZ
 202	FORMAT( '0SHIFT FACTORS:  XSHX = ', F6.3, ' XSHY = ', F6.3, 
	1  ' XSHZ = ', F6.3, ' YSHY = ', F6.3, ' YSHZ = ', F6.3, 
	2  ' ZSHZ = 'F6.3)
C	GET MAXIMUM AND MINUMUM MAP COORDINATES
	DO 4 I = 1, 3
	DO 6 J = 1, 3
 6	MO( I, J) = 0.0
 4	CONTINUE
	MO( 1, 1) = ABC( 1) * XSHX
	MO( 2, 2) = ABC( 2) * YSHY
	MO( 3, 3) = ABC( 3) * ZSHZ
	MO( 1, 2) = ABC( 2) * XSHY
	MO( 1, 3) = ABC( 3) * XSHZ
	MO( 2, 3) = ABC( 3) * YSHZ
	CALL ORGDIR( XYZO, MO)
	MDIR = 1
1210	CONTINUE
C      ROTATE PARAMETER LABELS TO STANDARD FORM

 1208	NDIR = MDIR
	NSDONE = 0
 25	CONTINUE
 106	FORMAT( 3F10.0, I5, I5)
	READ( IO2, 126)SCALE, ISCALE, RSCX, RSCY, XFRAM, YFRAM, 
	1 IFNFLD, XFNFLD, ISTRO, STRANG, MSEC
 126	FORMAT( F8.3, I5, 4F8.3, I5, F8.4, I5, F8.4, I5)
	WRITE( IO6, 371)SCALE, RSCX, RSCY, XFRAM, YFRAM, IFNFLD, 
	1 XFNFLD, ISTRO, STRANG, MSEC
 371	FORMAT( '  SCALE IS', F6.2, 'MM/A.  ', /, 
	1 '   DISTORTION CORRECTION IN VERT. AND HORIZ. = ', 2F10.3, /, 
	2 '   VERT. AND HORIZ. FRAME  = ', 2F10.3, ' MM.', /, 
	3 '   FAN FOLD PAPER?', I2, '  PAGE LENGTH  = ', F8.1, 'MM.', /, 
	4 '   STERO ?', I2, ' STERO ANGLE  = ', F5.2, I5, 
	5 ' DENSITY SECTIONS/SECTION', //)
 388	CONTINUE
	MSEC2 = MSEC / 2
	STRANG = DGTRD * STRANG * 0.5
	STROX = 0.00
	STROY = 0.00
	SSTROX = 0.00
	SSTROY = 0.0
	IF( ISTRO.EQ.1)STROX = TAN( STRANG)
	IF( ISTRO.EQ.2)STROY = TAN( STRANG)
	NSEC1 = ( XMAX1( 3) - XMIN1( 3)) / DX( 3) + 1.01
	IF( NSEC.GT.NSEC1.OR.NSEC.EQ.0)NSEC = NSEC1
	NZMX = NSEC1
 40	CONTINUE
C
	NDZ = 0
 510	CONTINUE
C	IINZ = ( XMIN1( 3) - XMIN( 3)) / DX( 3) + 0.1
	IINZ = IINZ - MSEC2 - 1
 500	CONTINUE
	DO 530 I = 1, 3
	XMAX( I) = XMAX1( I)
	XMIN( I) = XMIN1( I)
 530	CONTINUE
 520	CONTINUE
	MEXT = 1
	IF( NEXT .EQ. 1) NEXT = 3
	IF( ( NEXT .EQ. 0) .AND. ( LEXT .GT. 0)) NEXT = LEXT
C      SCALING SECTION
C      GET THE SHIFT FACTORS
C	CONVERT TO ANGSTROM SCALE
	DO 414 I = 1, 3
	XMIN( I) = XMIN1( I)
	XMAX( I) = XMAX1( I)
	XMIN2( I) = XMIN1( I)
	XMAX2( I) = XMAX1( I)
 414	CONTINUE
C	GET MAXIMUM AND MINUMUM MAP COORDINATES
	XMPMIN = XMIN( 1)
	XMPMAX = XMAX( 1)
	YMPMIN = XMIN( 2)
	YMPMAX = XMAX( 2)
	WRITE( IO6, 115) XMPMIN, XMPMAX, YMPMIN, YMPMAX
 115	FORMAT( '0A CARTESIAN COORDINATE SYSTEM FOR THE MAPS HAS ',
 	1 'BEEN ESTABLISHED'/' X RUNS FROM ', F7.2, ' TO ', F7.2/
	2 ' Y RUNS FROM ', F7.2, ' TO ', F7.2/
	3 ' THE COORDINATES OF EACH GRID POINT WILL BE TRANSFORMED ',
	4 'TO THIS SYSTEM BEFORE CONTOURING.')
 26	CONTINUE
C	GET SCALE FACTOR
	WRITE( IO6, 117) ISCALE
 117	FORMAT( '0GRID SPACING IS ', I3, 'A.')
	CALL MODEL1( NDIR, XYZO, MO)
	CLOSE( UNIT = IO2)
	RETURN
	END

C
C==========================================================================
C
	SUBROUTINE ORGDIR(  XYZO, MO)
C********************************
C	SUBROUTINE ORGDIR
C		ESTABLISHES ORIGIN AND DIRECTION OF DIAGRAM
C		XYZO(3) IS ORIGIN IN A.
C		ROTATION IS GIVEN BY
C			(XN)   (      )   ((X)   (XO))
C			(YN) = (  MR  ) * ((Y) - (YO))
C			(ZN)   (      )   ((Z)   (ZO))
C
C			WHERE XN, YN, ZN ARE NEW COORDS. IN A.
C			X, Y, Z ARE OLD COORDS. IN A.
C			AND XO, YO, ZO ARE ORIGIN COORDS. IN A.
C********************************************************
	COMMON /IOSET/ IO2, IO3, IO4, IO5, IO6, IO7, IO8
C
	REAL*8    NAME, ORIGIN
	REAL*4    MR(3, 3), MO(3, 3), X(3)
	REAL*4    VX1(3), VX(3), VY(3), VZ(3), MM(3, 3), XT(3), ANG(10)
	REAL*4    XYZY1(3), XYZX1(3), XYZX2(3)
	REAL*4    XYZO(3), XYZP(3, 50), XYZY2(3)
	INTEGER   NP(50), IAX(10), NP1(16)
	LOGICAL*1 ORGFIL(80)
	DATA ORIGIN/'ORIGIN'/
C
	READ( IO2, 13) IQ, ( ORGFIL( K), K = 1, IQ)
 13	FORMAT( Q, 80A1)
	OPEN( UNIT = IO4, NAME = ORGFIL( 1), TYPE = 'OLD', 
	1 FORM = 'FORMATTED', READONLY, SHARED)
	CALL ORTHIN( MO)
C	WRITE( IO6, 203) MO
	CALL ATMGET( 0, NN, NAME, IK, X( 1), X( 2), X( 3), BT, 
	1             IO4, IFOK)
	READ( IO2, 11) NO, NY1, NY2, NX1, NX2
	IF( NX1.LE.0 .AND. NX2.LE.0 )THEN
	   NX2 = NY2
	   NX1 = NO
	   NY2 = NY1
	   NY1 = NO
	ELSE IF( NX2.LE.0 )THEN
	   NX2 = NX1
	   NX1 = NY1
	END IF
	IF( NO.LE.0 )THEN
	   XYZO( 1) = 0.0
	   XYZO( 2) = 0.0
	   XYZO( 3) = 0.0
	   NN = NO
	   NAME = ORIGIN
	ELSE
	   CALL ATMGET( NO, NN, NAME, IK, X( 1), X( 2), X( 3), 
	1               BT, IO4, IFOK )
	   XYZO( 1) = X( 1)
	   XYZO( 2) = X( 2)
	   XYZO( 3) = X( 3)
	END IF
	WRITE( IO6, 21)
 21	FORMAT( '0 ORIGIN DEFINING ATOM  -  NO.  NAME')
	WRITE( IO6, 23) NN, NAME
 23	FORMAT( '0', 24X, I5, 2X, A8)
C
	WRITE( IO6, 25)
 25	FORMAT( '0 Y-AXIS DEFINING ATOMS -  NO.  NAME')
	CALL ATMGET( NY1, NN, NAME, IK, X( 1), X( 2), X( 3), BT, 
	1             IO4, IFOK )
	WRITE( IO6, 23)NN, NAME
	CALL ORTHO( X( 1), X( 2), X( 3), XYZY1( 1), XYZY1( 2), 
	1           XYZY1( 3) )
	CALL ATMGET( NY2, NN, NAME, IK, X(1), X(2), X(3), BT, 
	1            IO4, IFOK )
	WRITE(IO6, 23)NN, NAME
	CALL ORTHO( X(1), X(2), X(3), XYZY2(1), XYZY2(2), XYZY2(3) )
C
	WRITE( IO6, 27)
 27	FORMAT( '0 X-AXIS DEFINING ATOMS -  NO.  NAME')
	CALL ATMGET( NX1, NN, NAME, IK, X(1), X(2), X(3), BT, IO4, IFOK )
	WRITE( IO6, 23)NN, NAME
	CALL ORTHO( X(1), X(2), X(3), XYZX1(1), XYZX1(2), XYZX1(3) )
	CALL ATMGET( NX2, NN, NAME, IK, X(1), X(2), X(3), BT, IO4, IFOK )
	WRITE( IO6, 23)NN, NAME
	CALL ORTHO( X(1), X(2), X(3), XYZX2(1), XYZX2(2), XYZX2(3) )
C
	READ( IO2, 11) ( NP1( I), I = 1, 16 )
 11	FORMAT( 16I5)
	NPC = 1
C
C	THE LIST OF NUMBERS ON LINE FIVE OF VEWCC.DAT
C	ARE PLACED INTO A BUFFER ARRAY CALLED NP1, BEFORE
C	ACTUALLY GOING TO THE REST OF THE PROGRAM.
C	THIS IS SO THAT MORE NUMBERS CAN BE USED
C	IN DETERMINING THE LEAST SQUARES PLANE.
C	THE DO LOOP LOOKS FOR HYPHENS BY THINKING OF
C	THEM AS MINUS SIGNS.  NPC IS A COUNTER OF THE
C	NUMBER OF ACTUAL ATOMS USED TO CALCULATE THE
C	PLANE.
C
	DO I = 1, 16
	   IF( NP1(1).NE.0 )THEN
	      IF( NP1(I).LT.0 ) CALL FILL( NP, NP1, I, NPC )
	      NP( NPC) = ABS( NP1( I) )
	      NPC = NPC + 1
	   END IF
	END DO
	K = NPC - 1
C	SORT NP( I)
	DO I1 = 1, K
	   I11 = I1 + 1
	   DO 20 I2 = I11, NPC
	      IF( NP( I2).GT.NP( I1) )GO TO 20
	      NPT    = NP( I2)
	      NP( I2) = NP( I1)
	      NP( I1) = NPT
 20	   CONTINUE
	END DO
	NF = 0
	WRITE( IO6, 28)
 28	FORMAT( '0 PLANE  DEFINING ATOMS -  NO.  NAME')
	DO 40 N = 1, NPC
	   IF( NP( N).EQ.0)GO TO 40
	   CALL ATMGET( NP(N), NN, NAME, IK, X(1), X(2), X(3), BT, 
	1               IO4, IFOK )
	   IF( IFOK.LE.0 )THEN
	      WRITE( IO6, 31)NP( N)
 31	      FORMAT( ' ***** ERROR-REQUESTED ATOM NOT IN FILE-ATOM NO.',
 	1                 I5)
	      STOP
	   END IF
	   NF = NF + 1
	   WRITE( IO6, 23)NN, NAME
	   X(1) = X(1) - XYZO(1)
	   X(2) = X(2) - XYZO(2)
	   X(3) = X(3) - XYZO(3)
	   CALL ORTHO( X(1), X(2), X(3), XYZP(1, NF), XYZP(2, NF), 
	1               XYZP(3, NF))
 40	CONTINUE
	CLOSE( UNIT = IO4)
	DO K1 = 1, 3
	   DO K2 = 1, 3
	      MM( K1, K2) = 0.0
	   END DO
	   MM( K1, K1) = 1.0
	END DO
	IF( NF.EQ.0)GO TO 112
	CALL PLANE( NF, XYZP, VZ, DP, MM)
C	WRITE( IO6, 203) VZ, DP
	DVZ = SQRT( VZ(1)**2 + VZ(2)**2 + VZ(3)**2 )
	DVY = 0.0
	DVX = 0.0
	DO I = 1, 3
	   VZ(I) = VZ(I) / DVZ
	   VY(I) = XYZY2(I) - XYZY1(I)
	   VX(I) = XYZX2(I) - XYZX1(I)
	   DVX   = DVX + VX(I)**2
	   DVY   = DVY + VY(I)**2
	END DO
	DVY = SQRT( DVY)
	DVX = SQRT( DVX)
	DO K = 1, 3
	   VY(K)  = VY(K) / DVY
	   VX1(K) = VX(K) / DVX
	END DO
	CALL V1XV2( VY, VZ, VX)
	DVXVX1 = VX(1) * VX1(1) + VX(2) * VX1(2) + VX(3) * VX1(3)
	DVX    = SQRT( VX(1)**2 + VX(2)**2 + VX(3)**2)
	A1 = 1.0
	IF( DVXVX1.LT.0.0 ) A1 = -1.0
	DO K = 1, 3
	   VX(K) = A1 * VX(K) / DVX
	   VZ(K) = A1 * VZ(K)
	END DO
	CALL V1XV2( VZ, VX, VY)
	DO K = 1, 3
	   MR(1, K) = VX(K)
	   MR(2, K) = VY(K)
	   MR(3, K) = VZ(K)
	END DO
C	WRITE( IO6, 203) MR
C 203	FORMAT( 3F10.4)
	DO K1 = 1, 3
	   DO K2 = 1, 3
	      MM( K1, K2) = 0.0
	      DO K3 = 1, 3
	         MM( K1, K2) = MM( K1, K2) + MR( K1, K3) * MO( K3, K2)
	      END DO
	   END DO
	END DO
 112	CONTINUE
	DO K1 = 1, 3
	   DO K2 = 1, 3
	      MO( K1, K2) = MM( K1, K2)
	   END DO
	END DO
C	WRITE( IO6, 203) MO
	READ( IO2, 151) ( IAX(K), ANG(K), K = 1, 10)
 151	FORMAT( 10( I3, F5.0))
	DGTRD = 0.017453293
	DO 160 K = 1, 10
	   IF( IAX(K).EQ.0)GO TO 160
	   WRITE( IO6, 161) K, IAX(K), ANG(K)
 161	   FORMAT( '  ROTATION', I3, ', AXIS  = ', I3, ', ', F7.2, 'DEG')
	   IF( IAX(K).LT.0)ANG(K) = -ANG(K)
	   AA = ANG(K) * DGTRD
	   CSANG = COS( AA)
	   SNANG = SIN( AA)
	   DO K1 = 1, 3
	      DO K2 = 1, 3
	         MR( K1, K2) = 0.0
	      END DO
	      MR( K1, K1) = 1.0
	   END DO
	   GO TO ( 210, 220, 230) IABS( IAX(K))
 210	      MR( 2, 2) = CSANG
	      MR( 2, 3) = -SNANG
	      MR( 3, 2) = SNANG
	      MR( 3, 3) = CSANG
	      GO TO 240
 220	      MR( 1, 1) = CSANG
	      MR( 1, 3) = SNANG
	      MR( 3, 1) = -SNANG
	      MR( 3, 3) = CSANG
	      GO TO 240
 230	      MR( 1, 1) = CSANG
	      MR( 1, 2) = -SNANG
	      MR( 2, 1) = SNANG
	      MR( 2, 2) = CSANG
 240	      CONTINUE
	   DO K1 = 1, 3
	      DO K2 = 1, 3
	         MM( K1, K2) = 0.0
	         DO K3 = 1, 3
	            MM( K1, K2) = MM( K1, K2) + MR( K1, K3) * MO( K3, K2)
	         END DO
	      END DO
	   END DO
	   DO K1 = 1, 3
	      DO K2 = 1, 3
	         MO( K1, K2) = MM( K1, K2)
	      END DO
	   END DO
C	   WRITE( IO6, 203) MO
 160	CONTINUE
	RETURN
	END

C
C==========================================================================
C
	SUBROUTINE V1XV2( V1, V2, V3)
	DIMENSION V1(3), V2(3), V3(3)
	V3(1) = V1(2) * V2(3) - V1(3) * V2(2)
	V3(2) = V1(3) * V2(1) - V1(1) * V2(3)
	V3(3) = V1(1) * V2(2) - V1(2) * V2(1)
	RETURN
	END

C
C==========================================================================
C
	SUBROUTINE SRFDRW( SSTROX, SSTROY, STML0, STMT0, MSEC, STROSN, 
	1           SCX, SCY, STROX, STROY)
C
	COMMON /IOSET/ IO2, IO3, IO4, IO5, IO6, IO7, IO8, MAXTYP
C
	COMMON /MOD/ AM(3, 3, 36), AT(3, 36), XMIN3(3), XMAX3(3), 
	1 DZ3, DZ3N, DZ2, DZ2N, DMX(4), CHAR(4, 20), NSYM, PWIDE(4), 
	2 ICHR(4), HT(4), HT2(4), SCALE, RAD(4, 20), RDEL(4), 
	3 IMASK(36, 4), NSET(4, 20), ATMFIL(80, 4), VRAD(4, 20), 
	4 VTYP(4, 20), IFVDW(4), IFVDW0, IRESOL, IPEN(20,4)
C
	COMMON /DAT/ TI(10), XMIN(3), XMAX(3), DX(3), XCHAR(3), ACHAR(3), 
	1  IDIM(2), JDIM(2), NEXT, XSHX, XSHY, XSHZ, YSHY, YSHZ, ZSHZ, 
	3 XMPMIN, XMPMAX, YMPMIN, YMPMAX, ISCALE, NSDONE, NIP, ZMAP, A, 
	3  RMIN, RMAX, RINC, NSEC, ANG, RC
C
	COMMON /ATMLST/ XI(3, 4500), KI(4500), J
C
	DIMENSION X(3), X1(3)
	DIMENSION YIN(20, 2), YTR(20, 2)
	LOGICAL*1 KI, ATMFIL
	INTEGER   PWIDE, CHAR, VTYP, IFCOLR(10)
C
	IIN = 15
C
	IF( NSYM.EQ.0)RETURN
	CALL CIRCLE( IRESOL )
	CALL IFSPEN( 1)
	DO 100 II = 1, 4
     	IF( IFCTRL_C() ) RETURN
	IF( IFVDW(II).NE.3 .AND. IFVDW(II).NE.5 )GO TO 100
	IF( DMX(II).EQ.0.0.AND.HT(II).EQ.0.0)GO TO 100
	DMXX = 0.0
	DO I1 = 1, 20
	   DO I2 = 1, 20
	      IF( VRAD( II, I1) + VRAD( II, I2).GT.DMXX)
	1	DMXX = VRAD( II, I1) + VRAD( II, I2)
	   END DO
	END DO
	DMXX     = DMXX + 2.0
	DMX2     = DMXX**2
	XMIN3(3) = -ABS( DMX(II)) - DZ2
	XMAX3(3) = ABS( DMX(II)) + DZ2
	ZMAPO    = ZMAP
	J        = 0
C
     	IF( IFCTRL_C() ) RETURN
	CALL ATMLST( ZMAPO, II)
C	
	JN = J
C
C	JN ATOMS TO BE INCLUDED HAVE BEEN COLLECTED.
C
	IF( JN.EQ.0)GO TO 100
C
C	CYCLE OVER ATOM TYPES...   !R.S. ZHOU
C
	DO 99 JTYPE = 1, MAXTYP
     	    ICOLOR = 0
	DO 1000 ISEC = 1, MSEC
     	IF( IFCTRL_C() ) RETURN
	DZ = 0.5 * FLOAT( 2 * ISEC - 1 - MSEC)
	STML = STML0 - DZ * SSTROY * STROSN
	STMT = STMT0 - DZ * SSTROX * STROSN
	DZ = DZ * DX(3)
	CALL PLOT1Z_OFFSET( STMT, SCX, STML, SCY)
	DO 40 J = 1, JN
     	IF( IFCTRL_C() ) RETURN
     	IF( KI(J).NE.JTYPE ) GOTO 40
     	JUNK = VTYP(II,KI(J))   !CIRCLE TYPE
     	IF( JUNK.LT.1 .OR. JUNK.GT.9 .OR. JUNK.EQ.4
     + .OR. JUNK.EQ.5 .OR. JUNK.EQ.7 ) GOTO 40
     	IF( ICOLOR.EQ.0 ) THEN
     	    ICOLOR = 1
     	    CALL IFSPEN( IPEN(JTYPE,II) )
     	ENDIF
C
C	ESTABLISH LIMITS FOR ATOM J.
C	NOTE:	CIRCLE IS DEVIDED INTO TO HALVES TOP AND BOTTOM
C		FOR TOP   - YIN IS LEFT  AND YTR IS RIGHT
C		FOR BOTTOM- YIN IS RIGHT AND YTR IS LEFT
C
	IF( ( XI( 3, J)).GT.DZ)GO TO 40
	Z0 = XI( 3, J) - DZ
	X0 = XI( 1, J)
	Y0 = XI( 2, J)
C	X0 = X0 + Z0 * STROX
C	Y0 = Y0 + Z0 * STROY
	IK = KI(J)
	IF( MOD( VTYP( II, IK), 10).EQ.0)GO TO 40
	R0 = VRAD( II, IK)**2 - ( XI( 3, J) - DZ)**2
	IF( R0.LE.0.0)GO TO 40
	DO 42 K = 1, IIN
	YIN( K, 1) = 0.0
	YTR( K, 1) = 0.0
	YIN( K, 2) = 0.0
	YTR( K, 2) = 0.0
 42	CONTINUE
	R0 = SQRT( R0)
	XIN = X0 - R0
	YIN( 1, 1) = Y0 - R0
	XTR = X0 + R0
	YTR( 1, 1) = Y0 + R0
	YIN( 1, 2) = Y0 + R0
	YTR( 1, 2) = Y0 - R0
 60	IF( YIN( 1, 1).GT.YMPMIN)GO TO 62
	IF( YTR( 1, 1).LE.YMPMIN)GO TO 40
	YIN( 1, 1) = YMPMIN
	YTR( 1, 2) = YMPMIN
	IF( YMPMIN.LT.Y0)GO TO 62
	DDX = SQRT( R0**2 - ( YIN( 1, 1) - Y0)**2)
	XIN = X0 - DDX
	XTR = X0 + DDX
 62	IF( YTR( 1, 1).LT.YMPMAX)GO TO 64
	IF( YIN( 1, 1).GT.YMPMAX)GO TO 40
	YTR( 1, 1) = YMPMAX
	YIN( 1, 2) = YMPMAX
	IF( Y0.LE.YMPMAX)GO TO 64
	DDX = SQRT( R0**2 - ( YTR( 1, 1) - Y0)**2)
	XIN = X0 - DDX
	XTR = X0 + DDX
 64	IF( XIN.GT.XMPMIN)GO TO 66
	IF( XTR.LT.XMPMIN)GO TO 40
	XIN = XMPMIN
	IF( X0.LT.XMPMIN)GO TO 65
	YTR( 2, 1) = YTR( 1, 1)
	DY = SQRT( R0**2 - ( X0 - XMPMIN)**2)
	YTR( 1, 1) = Y0 - DY
	YIN( 2, 1) = Y0 + DY
	GO TO 66
 65	YIN( 1, 1) = 0.0
	YTR( 1, 1) = 0.0
	DY = SQRT( R0**2 - ( X0 - XMPMIN)**2)
	YTR( 1, 2) = Y0 - DY
	YIN( 1, 2) = Y0 + DY
 66	IF( XTR.LT.XMPMAX)GO TO 68
	IF( XIN.GT.XMPMAX)GO TO 40
	XTR = XMPMAX
	IF( X0.GT.XMPMAX)GO TO 67
	YTR( 2, 2) = YTR( 1, 2)
	DY = SQRT( R0**2 - ( XMPMAX - X0)**2)
	YIN( 2, 2) = Y0 - DY
	YTR( 1, 2) = Y0 + DY
	GO TO 68
 67	YTR( 1, 2) = 0.0
	YIN( 1, 2) = 0.0
	DY = SQRT( R0**2 - ( XMPMAX - X0)**2)
	YIN( 1, 1) = Y0 - DY
	YTR( 1, 1) = Y0 + DY
 68	CONTINUE
C
C	LIMITS ESTABLISHED.
C	DETERMINE CONTACTING ATOMS.
C
	DO  160 J2 = 1, JN
	IF( J.EQ.J2)GO TO 160
	Z2 = XI( 3, J2) - DZ
	X2 = XI( 1, J2) + Z2 * STROX
	D  = ( X0 - X2)**2
	IF( D.GT.DMX2)GO TO 160
	Y2 = XI( 2, J2) + Z2 * STROY
	D  = D + ( Y0 - Y2)**2
	IF( D.GT.DMX2 )GO TO 160
 162	CONTINUE
	IF( XI( 3, J).GT.XI( 3, J2)) GO TO 410
C
C	ATOM-J2 IN FRONT OF ATOM-J
C
	D = SQRT( D)
	KJ1 = KI(J)
	KJ2 = KI(J2)
	IF( D.GT.R0 + VRAD( II, KJ2))GO TO 160
C
C	ATOMS OVERLAP.
C
	   R2 = VRAD( II, KJ2)
	   D = D**2 + ( XI( 3, J) - XI( 3, J2))**2
	IF( D.GT.DMX2)GO TO 420
	   D = SQRT( D)
	IF( D.GT.( R0 + VRAD( II, KJ2)))GO TO 420
	IF( XI( 3, J2).GT.DZ)GO TO 420
	GO TO 430
 410	CONTINUE
C
C	ATOM-J IN FRONT OF ATOM-J2
C
	D = D + ( XI( 3, J) - XI( 3, J2))**2
	IF( D.GT.DMX2)GO TO 160
	D = SQRT( D)
	KJ1 = KI( J)
	KJ2 = KI( J2)
	IF( D.GT.R0 + VRAD( II, KJ2))GO TO 160
C
C	ATOMS IN CONTACT.
C
 430	R2 = VRAD( II, KJ2)**2 - ( XI( 3, J2) - DZ)**2
	IF( R2.LE.0.0)GO TO 160
	R2 = SQRT( R2)
 420	RR = ( X2 - X0)**2 + ( Y2 - Y0)**2
	RR = SQRT( RR)
C
C	IF ATOM-0 TOTALLY INSIDE ATOM-2 GO TO NEXT ATOM.
C
C	IF( ( R0 + RR).LE.R2 .AND. Z0.LT.Z2 )GO TO 40
	IF( ( R0 + RR).LE.R2)GO TO 40
C
C	FIND INTERSECTION OF CIRCLES. NOTE X3 < OR = X4.
C
	CALL INTRSC( X0, Y0, R0, X2, Y2, R2, X3, Y3, X4, Y4)
	IF( X3.EQ.0.0.AND.X4.EQ.0.0.AND.Y3.EQ.0.0.AND.Y4.EQ.0.0)GO TO 160
C
C	IF ATOM-0 INSIDE CIRCLE OF ATOM-2 GO TO 240.
C
	IF( RR.LT.R2)GO TO 240
C
C	ATOM-0 OUTSIDE CIRCLE OF ATOM-2.
C	IF ATOM-2 RIGHT OF ATOM-0 GO TO 170.
	IF( Y2.GT.Y0)GO TO 170
	IF( X4.LT.X0)GO TO 180
	IF( X3.GT.X0)GO TO 190
C
C	ATOM-2 LEFT OF ATOM-0.
C
	DO 202 K = 1, IIN
	IF( YIN( K, 1).EQ.0.0.AND.YTR( K, 1).EQ.0.0)GO TO 204
	IF( YIN( K, 1).GE.Y3)GO TO 204
	YIN( K, 1) = Y3
	IF( YIN( K, 1).LT.YTR( K, 1))GO TO 204
C
C	YIN AND YTR HAVE CROSSED, THIS PORTION BLANKED OUT.
C
	YIN( K, 1) = 0.0
	YTR( K, 1) = 0.0
 204	IF( YIN( K, 2).EQ.0.0.AND.YTR( K, 2).EQ.0.0)GO TO 202
	IF( YTR( K, 2).GE.Y4)GO TO 202
	YTR( K, 2) = Y4
	IF( YTR( K, 2).LT.YIN( K, 2))GO TO 202
C
C	YIN AND YTR HAVE CROSSED OVER, THIS PORTION OBLITERATED.
C
	YIN( K, 2) = 0.0
	YTR( K, 2) = 0.0
 202	CONTINUE
	GO TO 164
 180	CONTINUE
C
C	ATOM-2 ABOVE ATOM-0.
C	PLACE Y3 ON LEFT AND Y4 ON RIGHT.
C
	Y = Y3
	Y3 = Y4
	Y4 = Y
	Y = X3
	X3 = X4
	X4 = Y
	GO TO 210
 190	CONTINUE
C
C	ATOM-2 BELOW ATOM-0 (>X).
C	PLACE Y4 ON LEFT AND Y3 ON RIGHT.
C
	Y = Y3
	Y3 = Y4
	Y4 = Y
	Y = X3
	X3 = X4
	X4 = Y
	GO TO 220
 170	CONTINUE
C
C	ATOM-2 RIGHT OF ATOM-0.
C
	IF( X4.LT.X0)GO TO 210
	IF( X3.GT.X0)GO TO 220
	DO 232 K = 1, IIN
	IF( YIN( K, 1).EQ.0.0.AND.YTR( K, 1).EQ.0.0)GO TO 234
	IF( YTR( K, 1).LE.Y3)GO TO 234
	YTR( K, 1) = Y3
	IF( YTR( K, 1).GT.YIN( K, 1))GO TO 234
C
C	YIN AND YTR CROSSED OVER, PORTION OBLITERATED.
C
	YIN( K, 1) = 0.0
	YTR( K, 1) = 0.0
 234	IF( YIN( K, 2).EQ.0.0.AND.YTR( K, 2).EQ.0.0)GO TO 232
	IF( YIN( K, 2).LE.Y4)GO TO 232
	YIN( K, 2) = Y4
	IF( YTR( K, 2).LT.YIN( K, 2))GO TO 232
C
C	YIN AND YTR CROSSED OVER, PORTION OBLITERATED.
C
	YIN( K, 2) = 0.0
	YTR( K, 2) = 0.0
 232	CONTINUE
	GO TO 164
 210	CONTINUE
C
C	ATOM-2 ABOVE ATOM-0 (X2&X3&X4 < X0)
C	ALL IN TOP HALF - Y3 LEFT, Y4 RIGHT
	DO 212 K = 1, IIN
	IF( YIN( K, 1).EQ.0.0.AND.YTR( K, 1).EQ.0.0)GO TO 212
	IF( YIN( K, 1).LT.Y3.AND.YTR( K, 1).GT.Y4)GO TO 214
	IF( Y3.LT.YTR( K, 1).AND.Y4.GT.YTR( K, 1))YTR( K, 1) = Y3
	IF( Y3.LT.YIN( K, 1).AND.Y4.GT.YIN( K, 1))YIN( K, 1) = Y4
	IF( YIN( K, 1).LT.YTR( K, 1))GO TO 212
C
C	YIN AND YTR CROSSED OVER, PORTION OBLITERATED.
C
	YIN( K, 1) = 0.0
	YTR( K, 1) = 0.0
	GO TO 212
C
C	SEGMENT MUST BE DEVIDED.
C
 214	K1 = K
 216	K1 = K1 + 1
	IF( K1.GT.10)GO TO 212
	IF( YIN( K1, 1).NE.0.0.OR.YTR( K1, 1).NE.0.0)GO TO 216
	YTR( K1, 1) = YTR( K, 1)
	YTR( K, 1) = Y3
	YIN( K1, 1) = Y4
 212	CONTINUE
	GO TO 164
 220	CONTINUE
C
C	ALL IN BOTTOM HALF - Y4 LEFT, Y3 RIGHT
C
	DO 222 K = 1, IIN
	IF( YIN(K, 2).EQ.0.0.AND.YTR(K, 2).EQ.0.0)GO TO 222
	IF( YIN(K, 2).GT.Y3.AND.YTR(K, 2).LT.Y4)GO TO 224
	IF( Y3.GT.YIN(K, 2).AND.Y4.LT.YIN(K, 2))YIN(K, 2) = Y4
	IF( Y3.GT.YTR(K, 2).AND.Y4.LT.YTR(K, 2))YTR(K, 2) = Y3
	IF( YTR(K, 2).LT.YIN(K, 2))GO TO 222
C
C	YIN AND YTR CROSSED OVER, PORTION OBLITERATED.
C
	YTR(K, 2) = 0.0
	YIN(K, 2) = 0.0
	GO TO 222
C
C	SEGMENT MUST BE DEVIDED.
C
 224	K1 = K
 226	K1 = K1 + 1
	IF( K1.GT.10)GO TO 222
	IF( YIN(K1, 2).NE.0.0.OR.YTR(K1, 2).NE.0.0)GO TO 226
	YTR(K1, 2) = YTR(K, 2)
	YTR(K, 2) = Y3
	YIN(K1, 2) = Y4
 222	CONTINUE
	GO TO 164
C
C	ATOM-0 INSIDE CIRCLE OF ATOM-2.
C
 240	IF( Y2.LT.Y0)GO TO 270
C
C	ATOM-2 RIGHT OF ATOM-0.
C
	IF( X4.LT.X0)GO TO 280
	IF( X3.GT.X0)GO TO 290
C
C	ATOM-2 RIGHT OF ATOM-0.
C
	DO 302 K = 1, IIN
	IF( YIN(K, 1).EQ.0.0.AND.YTR(K, 1).EQ.0.0)GO TO 304
	IF( YTR(K, 1).LE.Y3)GO TO 304
	YTR(K, 1) = Y3
	IF( YIN(K, 1).LT.YTR(K, 1))GO TO 304
C
C	YIN AND YTR CROSSED OVER, PORTION OBLITERATED.
C
	YIN(K, 1) = 0.0
	YTR(K, 1) = 0.0
 304	IF( YIN(K, 2).EQ.0.0.AND.YTR(K, 2).EQ.0.0)GO TO 302
	IF( YIN(K, 2).LE.Y4)GO TO 302
	YIN(K, 2) = Y4
	IF( YTR(K, 2).LT.YIN(K, 2))GO TO 302
C
C	YIN AND YTR CROSSED OVER, PORTION OBLITERATED.
C
	YIN(K, 2) = 0.0
	YTR(K, 2) = 0.0
 302	CONTINUE
	GO TO 164
 280	CONTINUE
C
C	ATOM-2 ABOVE ATOM-0 (X2&X3&X4 < X0)
C	PLACE Y3 ON LEFT AND Y4 ON RIGHT.
C
	Y = Y3
	Y3 = Y4
	Y4 = Y
	Y = X3
	X3 = X4
	X4 = Y
	GO TO 310
 290	CONTINUE
C
C	INTERSECTION BELOW ATOM-0 ( X3&X4 > X0 )
C	PLACE Y4 ON LEFT AND Y3 ON RIGHT.
C
	Y = Y3
	Y3 = Y4
	Y4 = Y
	Y = X3
	X3 = X4
	X4 = Y
	GO TO 320
 270	CONTINUE
C
C	ATOM-2 LEFT OF ATOM-0.
C
	IF( X4.LT.X0)GO TO 310
	IF( X3.GT.X0)GO TO 320
C
C	ATOM-2 LEFT OF ATOM-0.
C
	DO 332 K = 1, IIN
	IF( YIN(K, 1).EQ.0.0.AND.YTR(K, 1).EQ.0.0)GO TO 334
	IF( YIN(K, 1).GE.Y3)GO TO 334
	YIN(K, 1) = Y3
	IF( YTR(K, 1).GT.YIN(K, 1))GO TO 334
C
C	YIN AND YTR CROSSED OVER, PORTION OBLITERATED.
C
	YIN(K, 1) = 0.0
	YTR(K, 1) = 0.0
 334	IF( YIN(K, 2).EQ.0.0.AND.YTR(K, 2).EQ.0.0)GO TO 332
	IF( YTR(K, 2).GE.Y4)GO TO 332
	YTR(K, 2) = Y4
	IF( YTR(K, 2).LT.YIN(K, 2))GO TO 332
C
C	YIN AND YTR CROSSED OVER, PORTION OBLITERATED.
C
	YIN(K, 2) = 0.0
	YTR(K, 2) = 0.0
 332	CONTINUE
	GO TO 164
 310	CONTINUE
C
C	ALL IN TOP HALF - Y3 LEFT, Y4 RIGHT
C
	IF( X2.LT.X0)GO TO 180
C
C	BOTTOM HALF OBLITERATED.
C
	DO 312 K = 1, IIN
	YIN(K, 2) = 0.0
	YTR(K, 2) = 0.0
	IF( YIN(K, 1).EQ.0.0.AND.YTR(K, 1).EQ.0.0)GO TO 312
	IF( YIN(K, 1).GT.Y3.AND.YTR(K, 1).LT.Y4)GO TO 312
	IF( Y3.GT.YIN(K, 1) )YIN(K, 1) = Y3
	IF( Y4.LT.YTR(K, 1) ) YTR(K, 1) = Y4
	IF( YIN(K, 1).LT.YTR(K, 1))GO TO 312
C
C	YIN AND YTR CROSSED OVER, PORTION OBLITERATED.
C
	YIN(K, 1) = 0.0
	YTR(K, 1) = 0.0
 312	CONTINUE
	GO TO 164
 320	CONTINUE
C
C	ALL IN BOTTOM HALF - Y3 LEFT, Y4 RIGHT
C
	IF( X2.GT.X0)GO TO 190
	DO 322 K = 1, IIN
C
C	TOP HALF OBLITERATED.
C
	YIN(K, 1) = 0.0
	YTR(K, 1) = 0.0
	IF( YIN(K, 2).EQ.0.0.AND.YTR(K, 2).EQ.0.0)GO TO 322
	IF( YIN(K, 2).LT.Y4.AND.YTR(K, 2).GT.Y3)GO TO 322
	IF( Y3.LT.YIN(K, 2) ) YIN(K, 2) = Y3
	IF( Y4.GT.YTR(K, 2) ) YTR(K, 2) = Y4
	IF( YTR(K, 2).LT.YIN(K, 2))GO TO 322
C
C	YIN AND YTR CROSSED OVER, PORTION OBLITERATED.
C
	YTR(K, 2) = 0.0
	YIN(K, 2) = 0.0
 322	CONTINUE
 164	CONTINUE
	DO 350 K = 1, IIN - 1
	IF( YTR(K, 1).GT.YIN(K, 1))GO TO 354
	YIN(K, 1) = YIN(K + 1, 1)
	YIN(K + 1, 1) = 0.0
	YTR(K, 1) = YTR(K + 1, 1)
	YTR(K + 1, 1) = 0.0
 354	IF( YIN(K, 2).GT.YTR(K, 2))GO TO 350
	YIN(K, 2) = YIN(K + 1, 2)
	YIN(K + 1, 2) = 0.0
	YTR(K, 2) = YTR(K + 1, 2)
	YTR(K + 1, 2) = 0.0
 350	CONTINUE
 160	CONTINUE
	DO 90 K = 1, IIN
	Y1 = YIN(K, 1)
	Y2 = YTR(K, 1)
	Y3 = YIN(K, 2)
	Y4 = YTR(K, 2)
	IF( Y1.LT.YMPMIN) Y1 = YMPMIN
	IF( Y2.GT.YMPMAX) Y2 = YMPMAX
	IF( Y3.GT.YMPMAX) Y3 = YMPMAX
	IF( Y4.LT.YMPMIN) Y4 = YMPMIN
	IF( YIN(K, 1).NE.0.0)
	1     CALL CIRDRW(X0, Y0, R0, Y1, Y2, MOD(VTYP(II, IK), 10))
	IF( YIN(K, 2).NE.0.0)
	1     CALL CIRDRW(X0, Y0, R0, Y3, Y4, MOD(VTYP(II, IK), 10))
	YIN(K, 1) = 0.0
	YIN(K, 2) = 0.0
	YTR(K, 1) = 0.0
	YTR(K, 2) = 0.0
 90	CONTINUE
 40	CONTINUE
 1000	CONTINUE
  99    CONTINUE
 100	CONTINUE
	CALL IFSPEN( 1)
	RETURN
	END

C
C==========================================================================
C
	SUBROUTINE EDGDRW(SSTROX, SSTROY, STML0, STMT0, MSEC, STROSN, 
	1         SCX, SCY, STROX, STROY)
C
	COMMON /IOSET/ IO2, IO3, IO4, IO5, IO6, IO7, IO8, MAXTYP
C
	COMMON /MOD/ AM(3, 3, 36), AT(3, 36), XMIN3(3), XMAX3(3), 
	1 DZ3, DZ3N, DZ2, DZ2N, DMX(4), CHAR(4, 20), NSYM, PWIDE(4), 
	2 ICHR(4), HT(4), HT2(4), SCALE, RAD(4, 20), RDEL(4), 
	3 IMASK(36, 4), NSET(4, 20), ATMFIL(80, 4), VRAD(4, 20), 
	4 VTYP(4, 20), IFVDW(4), IFVDW0, IRESOL, IPEN(20,4)
C
	COMMON /DAT/ TI(10), XMIN(3), XMAX(3), DX(3), XCHAR(3), ACHAR(3), 
	1  IDIM(2), JDIM(2), NEXT, XSHX, XSHY, XSHZ, YSHY, YSHZ, ZSHZ, 
	3 XMPMIN, XMPMAX, YMPMIN, YMPMAX, ISCALE, NSDONE, NIP, ZMAP, A, 
	3  RMIN, RMAX, RINC, NSEC, ANG, RC
C
	COMMON /ATMLST/ XI(3, 4500), KI(4500), J
C
	COMMON /SEAM/ IQ, U0(3, 50), U1(3, 50), U2(3, 50), RU0(50), 
	1 PHII(50), PHIT(50)
C
	COMMON /EDGE/THI(50), THT(50), NA(50), NT(50), 
	1  XN(50), YN(50), ZN(50), RN(50), NN
C
	COMMON /ANGL/PI, TWOPI, DT, YUN(126)
C
	DIMENSION X(3), X1(3)
	LOGICAL*1 KI, ATMFIL
	INTEGER   PWIDE, CHAR, VTYP
	DATA PI/3.141592654/
C
	IF( NSYM.EQ.0 ) RETURN
	TWOPI = PI * 2.0
	THREPI = PI * 3.0
     	IF( IFCTRL_C() ) RETURN
	CALL ANGTBL
	DO 100 II = 1, 4
     	IF( IFCTRL_C() ) RETURN
	IF( IFVDW(II).NE.2 .AND. IFVDW(II).NE.5 )GO TO 100
	IF( DMX(II).EQ.0.0.AND.HT(II).EQ.0.0 )GO TO 100
	DMXX = 0.0
	DO I1 = 1, 20
	   DO I2 = 1, 20
	      IF( VRAD(II, I1) + VRAD(II, I2).GT.DMXX )
	1	DMXX = VRAD(II, I1) + VRAD(II, I2)
 	   END DO
	END DO
	DMXX     = DMXX + 3.0
	DMX2     = DMXX**2
	XMIN3(3)  = -ABS(DMX(II)) - DZ2
	XMAX3(3) = ABS(DMX(II)) + DZ2
	ZMAPO    = ZMAP
	J        = 0
C
	CALL ATMLST( ZMAPO, II)
C	
	JN = J
C
C	JN ATOMS TO BE INCLUDED HAVE BEEN COLLECTED.
C
	IF( JN.EQ.0 )GO TO 100
C
C  CYCLE OVER ATOM TYPE...   !R.S. ZHOU
C
     	DO JTYPE = 1, MAXTYP
     	   ICOLOR = 0
	DO 40 J = 1, JN
     	IF( IFCTRL_C() ) RETURN
     	IF( KI(J).NE.JTYPE ) GOTO 40
     	IF( ICOLOR.EQ.0 ) THEN
     	    ICOLOR = 1
     	    CALL IFSPEN( IPEN(JTYPE,II) )
     	ENDIF
C
C	ESTABLISH LIMITS FOR ATOM J.
C
	X0   = XI(1, J)
	Y0   = XI(2, J)
	DZ   = XI(3, J)
	Z0   = 0.0
	STML = STML0 - DZ * STROY
	STMT = STMT0 - DZ * STROX
	CALL PLOT1Z_OFFSET(STMT, SCX, STML, SCY)
	IK = KI(J)
	IF( VTYP(II, IK).EQ.0 )GO TO 40
	R0 = VRAD(II, IK)
	NN = 0
	DO K = 1, 50
	   THI(K) = 0.0
	   THT(K) = 0.0
	   NT(K)  = 0
	END DO
	XIN = X0 - R0
	YIN = Y0 - R0
	XTR = X0 + R0
	YTR = Y0 + R0
	IQ  = 0
 60	IF( YIN.GT.YMPMIN )GO TO 62
	IF( YTR.LE.YMPMIN )GO TO 40
	YIN = YMPMIN
	Y3  = YIN - Y0
	DDX = SQRT( R0**2 - (Y3)**2 )
	X3  =  DDX
	TH1 = ATAN2(X3, Y3)
	IF( TH1.GT.PI )TH1 = TH1 - TWOPI
	TH2 = TWOPI - TH1
	CALL THSET(TH1, TH2, 2, 1, 0.0, YIN, 0.0, 9999.9)
 62	IF( YTR.LT.YMPMAX )GO TO 64
	IF( YIN.GT.YMPMAX )GO TO 40
	YTR = YMPMAX
	Y3  = YTR - Y0
	DDX = SQRT( R0**2 - (Y3)**2 )
	X3  =  -DDX
	TH1 = ATAN2(X3, Y3)
	TH2 =  -TH1
	CALL THSET(TH1, TH2, 2, 2, 0.0, YTR, 0.0, 9999.9)
 64	IF( XIN.GT.XMPMIN )GO TO 66
	IF( XTR.LT.XMPMIN )GO TO 40
	XIN = XMPMIN
	X3  = XIN - X0
	DY  = SQRT( R0**2 - X3**2 )
	Y3  = -DY
	TH1 = ATAN2(X3, Y3)
	IF( TH1.LT.0.0 )TH2 = -PI - TH1
	IF( TH1.GE.0.0 )TH2 = THREPI - TH1
	CALL THSET(TH1, TH2, 2, 3, XIN, 0.0, 0.0, 9999.9)
 66	IF( XTR.LT.XMPMAX )GO TO 68
	IF( XIN.GT.XMPMAX )GO TO 40
	XTR = XMPMAX
	X3  = XTR - X0
	DY  = SQRT(R0**2 - X3**2 )
	Y3  = DY
	TH1 = ATAN2(X3, Y3)
	TH2 = PI - TH1
	CALL THSET(TH1, TH2, 2, 4, XTR, 0.0, 0.0, 9999.9)
 68	CONTINUE
C
C	LIMITS ESTABLISHED.
C	DETERMINE CONTACTING ATOMS.
C
	DO  160 J2 = 1, JN
     	IF( IFCTRL_C() ) RETURN
	IFRNT = 0
	NNN = 0
	IF( J.EQ.J2 )GO TO 160
	D = 0.0
	DO K = 1, 2
	   D = D + (XI(K, J) - XI(K, J2))**2
	   IF( D.GT.DMX2)GO TO 160
	END DO
	X2  = XI(1, J2)
	Y2  = XI(2, J2)
	Z2  = XI(3, J2) - DZ
	X20 = X2
	Y20 = Y2
	KJ2 = KI(J2)
	R2  = VRAD(II, KJ2)
	X2  = X2 + (Z2) * STROX
	Y2  = Y2 + (Z2) * STROY
	D   = (X2 - X0)**2 + (Y2 - Y0)**2
	IF( D.GT.DMX2 )GO TO 160
	IF( XI(3, J).GT.XI(3, J2) )GO TO 410
C
C	ATOM-J2 IN FRONT OF ATOM-J
C
	D = SQRT(D)
	KJ1 = KI(J)
	KJ2 = KI(J2)
	IF( D.GT.VRAD(II, KJ1) + VRAD(II, KJ2) )GO TO 160
C
C	ATOMS OVERLAP.
C
	R2    = VRAD(II, KJ2)
	R21   = R2
	D     = D**2 + Z2**2
	NNN   = 4
	IFRNT = 1
	IF( D.GT.DMX2 )GO TO 420
	D     = SQRT(D)
	IF( D.GT.(VRAD(II, KJ1) + VRAD(II, KJ2)) )GO TO 420
C	IF( (VRAD(II, KJ1)**2 + VRAD(II, KJ2)**2 - Z2**2).LT.
C	1    ((X0 - X2)**2 + (Y0 - Y2)**2) )GO TO 420
	NNN   = 3
	GO TO 420
 410	CONTINUE
C
C	ATOM-J IN FRONT OF ATOM-J2
C
	NNN = 1
	D   = D + Z2**2
	IF( D.GT.DMX2 )GO TO 160
	D   = SQRT(D)
	KJ1 = KI(J)
	KJ2 = KI(J2)
	IF( D.GT.VRAD(II, KJ1) + VRAD(II, KJ2) )GO TO 160
C
C	ATOMS IN CONTACT.
C
 430	R21 = VRAD(II, KJ2)**2 - Z2**2
	IF( R21.LE.0.0 )GO TO 160
	R21 = SQRT(R21)
 420	RR  = (X2 - X0)**2 + (Y2 - Y0)**2
	RR  = SQRT(RR)
	IF( (RR + R21) .GT. R0 )GO TO 422
C
C	ATOM-J2 TOTALLY INSIDE ATOM-J
C
	TH1 = 0.0
	TH2 = 0.0
	NNN = 5
	GO TO 470
 422	CONTINUE
C
C	IF ATOM-J TOTALLY INSIDE ATOM-J2 GO TO NEXT ATOM.
C
	IF( ((R0 + RR).LE.R21).AND.(Z2.GT.0.0) )GO TO 40
C
C	FIND INTERSECTION OF CIRCLES. NOTE X3 < OR = X4.
C
	CALL INTRSC(X0, Y0, R0, X2, Y2, R21, X3, Y3, X4, Y4)
	TH1 = -PI * 0.5
	TH2 = PI * 1.5
	TH3 = 0.0
	IF( ((R0 + RR).LE.R21).AND.(Z2.LE.0.0) )GO TO 470
	IF( X3.EQ.0.0.AND.X4.EQ.0.0.AND.Y3.EQ.0.0.AND.Y4.EQ.0.0 )GO TO 160
	TH1 = ATAN2((X3 - X0), (Y3 - Y0))
	TH2 = ATAN2((X4 - X0), (Y4 - Y0))
	TH3 = ATAN2((X2 - X0), (Y2 - Y0))
	IF( (TH1 - TH3).GT.PI )TH1 = TH1 - TWOPI
	IF( (TH3 - TH1).GT.PI )TH1 = TH1 + TWOPI
	IF( (TH2 - TH3).GT.PI )TH2 = TH2 - TWOPI
	IF( (TH3 - TH2).GT.PI )TH2 = TH2 + TWOPI
	IF( TH2.GT.TH3 )GO TO 440
	TH3 = TH2
	TH2 = TH1
	TH1 = TH3
 440	CONTINUE
	IF( TH1.GT.-PI ) GO TO 460
	TH1 = TH1 + TWOPI
	TH2 = TH2 + TWOPI
 460	IF( TH1.LT.PI )GO TO 470
	TH1 = TH1 - TWOPI
	TH2 = TH2 - TWOPI
 470	CONTINUE
C	CALL THSET(TH1, TH2, NNN, J2, X20, Y20, Z2, R2)
	CALL THSET(TH1, TH2, NNN, J2, X2, Y2, Z2, R2)
C
 160	CONTINUE
	CALL DRWEDG(X0, Y0, Z0, R0, STROX, STROY)
 40	CONTINUE
     	ENDDO
 100	CONTINUE
	RETURN
	END

C
C==========================================================================
C
	SUBROUTINE CIRCLE( IR )
	DIMENSION XUN(73), YUN(73)
	PI = 3.141592654
     	IF( IR.EQ.0 ) IRR = 18
     	IF( IR.EQ.1 ) IRR = 36
     	IF( IR.EQ.2 ) IRR = 72
     	DT = PI / FLOAT(IRR)
	DO I = 1, IRR+1
	   YUN(I) = -COS(DT * FLOAT(I - 1))
	   XUN(I) = SQRT(1.0 - YUN(I)**2)
	END DO
	RETURN
C
C==========================================================================
C
	ENTRY CIRDRW(X0, Y0, R0, Y1, Y2, ITYP)
	IF( (ITYP.LT.1) .OR. (ITYP.GT.9) ) ITYP = 7
	IF( Y2.GE.Y1 ) THEN
	   T1 = ( (Y1 - Y0) / R0)
	   T2 = ( (Y2 - Y0) / R0)
	   IF( T1.LT.-1.0) T1 = -1.0
	   IF( T1.GT.1.0)  T1 = 1.0
	   IF( T2.LT.-1.0) T2 = -1.0
	   IF( T2.GT.1.0)  T2 = 1.0
	   T1 = ACOS(-T1)
	   T2 = ACOS(-T2)
	   I1 = (T1 / DT) + 2.0
	   I2 = (T2 / DT) + 1.0
	   I3 = (3.0 / R0)
	   IF( ITYP.GT.7 ) I3 = ( 1.5 / R0 )
	   IF( I3.LT.1)I3 = 1
	   IF( I3.GT.8)I3 = 8
	   Y = Y1
	   X = (R0**2 - (Y - Y0)**2)
	   IF( X.LT.0.0)X = 0.0
	   X = X0 - SQRT(X)
	   GO TO (12, 112, 212, 312, 412, 512, 612, 112, 212, 612)ITYP
 12	   CONTINUE
C
C	   SOLID CIRCLES
C
	   CALL PLOT1Z_QPLOT(X, Y, 13)
	   IF( I2.LE.I1)GO TO 40
	   DO 30 I = I1, I2, I3
	   X = X0 - R0 * XUN(I)
	   Y = Y0 + R0 * YUN(I)
	   CALL PLOT1Z_QPLOT(X, Y, 12)
 30	   CONTINUE
 40	   CONTINUE
	   Y = Y2
	   X = (R0**2 - (Y - Y0)**2)
	   IF( X.LT.0.0)X = 0.0
	   X = X0 - SQRT(X)
	   CALL PLOT1Z_QPLOT(X, Y, 12)
	   RETURN
C
C
 112	   CONTINUE
C
C	   DOTTED CIRCLE
C
	   IF( I2.LE.I1)GO TO 140
	   DO I = I1, I2, I3
	      X = X0 - R0 * XUN(I)
	      Y = Y0 + R0 * YUN(I)
	      CALL PLOT1Z_QPLOT(X, Y, 13)
	      CALL PLOT1Z_QPLOT(X, Y, 12)
	   END DO
 140	   CONTINUE
	   Y = Y2
	   X = (R0**2 - (Y - Y0)**2)
	   IF( X.LT.0.0)X = 0.0
	   X = X0 - SQRT(X)
	   CALL PLOT1Z_QPLOT(X, Y, 13)
	   CALL PLOT1Z_QPLOT(X, Y, 12)
	   RETURN
C
C
 212	   CONTINUE
C
C	   DASHED CIRCLE
C
	   CALL PLOT1Z_QPLOT(X, Y, 13)
	   I4 = 13
	   I5 = 12
	   IF( I2.LE.I1)GO TO 240
	   DO 230 I = I1, I2, I3
	   X = X0 - R0 * XUN(I)
	   Y = Y0 + R0 * YUN(I)
	   CALL PLOT1Z_QPLOT(X, Y, I4)
	   I6 = I4
	   I4 = I5
	   I5 = I6
 230	   CONTINUE
 240	   CONTINUE
	   Y = Y2
	   X = (R0**2 - (Y - Y0)**2)
	   IF( X.LT.0.0)X = 0.0
	   X = X0 - SQRT(X)
	   CALL PLOT1Z_QPLOT(X, Y, I4)
	   RETURN
C
C
 312	   CONTINUE
C
C	   SAW TOOTH
C
	   CALL PLOT1Z_QPLOT(X, Y, 13)
	   K = -1
	   IF( I2.LE.I1)GO TO 340
	   DO 330 I = I1, I2, I3
	   K = -1 * K
	   X = X0 - (R0 + K * 0.05) * XUN(I)
	   Y = Y0 + (R0 + K * 0.05) * YUN(I)
	   CALL PLOT1Z_QPLOT(X, Y, 12)
 330	   CONTINUE
 340	   CONTINUE
	   Y = Y2
	   X = (R0**2 - (Y - Y0)**2)
	   IF( X.LT.0.0)X = 0.0
	   X = X0 - SQRT(X)
	   CALL PLOT1Z_QPLOT(X, Y, 12)
	   RETURN
C
C
 412 	   CONTINUE
C
C	   0.1*RO SPIKED CIRCLES
C
	   CALL PLOT1Z_QPLOT(X, Y, 13)
	   IF( I2.LE.I1)GO TO 440
	   R3 = R0 * 0.9
	   DO 430 I = I1, I2, I3
	   X = X0 - R0 * XUN(I)
	   Y = Y0 + R0 * YUN(I)
	   CALL PLOT1Z_QPLOT(X, Y, 12)
	   X3 = X0 - R3 * XUN(I)
	   Y3 = Y0 + R3 * YUN(I)
	   CALL PLOT1Z_QPLOT(X3, Y3, 12)
	   CALL PLOT1Z_QPLOT(X, Y, 13)
 430 	   CONTINUE
 440 	   CONTINUE
	   Y = Y2
	   X = (R0**2 - (Y - Y0)**2)
	   IF( X.LT.0.0)X = 0.0
	   X = X0 - SQRT(X)
	   CALL PLOT1Z_QPLOT(X, Y, 12)
	   RETURN
C
C
 512 	   CONTINUE
C
C	   0.15A SPIKED CIRCLES
C
	   CALL PLOT1Z_QPLOT(X, Y, 13)
	   IF( I2.LE.I1)GO TO 540
	   R4 = 0.15
	   DO 530 I = I1, I2, I3
	   X = X0 - R0 * XUN(I)
	   Y = Y0 + R0 * YUN(I)
	   CALL PLOT1Z_QPLOT(X, Y, 12)
	   X4 = X + R4 * XUN(I)
	   Y4 = Y - R4 * YUN(I)
	   CALL PLOT1Z_QPLOT(X4, Y4, 12)
	   CALL PLOT1Z_QPLOT(X, Y, 13)
 530 	   CONTINUE
 540 	   CONTINUE
	   Y = Y2
	   X = (R0**2 - (Y - Y0)**2)
	   IF( X.LT.0.0)X = 0.0
	   X = X0 - SQRT(X)
	   CALL PLOT1Z_QPLOT(X, Y, 12)
	   RETURN
C
C
 612 	   CONTINUE
C
C	   BLANK
C
	   RETURN
C
C
C
	ELSE
	   T2 = ((Y1 - Y0) / R0)
	   T1 = ((Y2 - Y0) / R0)
	   IF( T1.LT. - 1.0)T1 = -1.0
	   IF( T1.GT.1.0)T1 = 1.0
	   IF( T2.LT. - 1.0)T2 = -1.0
	   IF( T2.GT.1.0)T2 = 1.0
	   T1 = ACOS(-T1)
	   T2 = ACOS(-T2)
	   I1 = (T1 / DT) + 2.0
	   I2 = (T2 / DT) + 1.0
	   I3 = (3.0 / R0)
	   IF( ITYP.GT.7 ) I3 = ( 1.5 / R0 )
	   IF( I3.LT.1)I3 = 1
	   IF( I3.GT.8)I3 = 8
	   Y = Y2
	   X = (R0**2 - (Y - Y0)**2)
	   IF( X.LT.0.0)X = 0.0
	   X = X0 + SQRT(X)
	   GO TO (14, 114, 214, 314, 414, 514, 614, 114, 214, 614)ITYP
 14	   CONTINUE
	   CALL PLOT1Z_QPLOT(X, Y, 13)
	   IF( I2.LE.I1)GO TO 60
	   DO 50 I = I1, I2, I3
	   X = X0 + R0 * XUN(I)
	   Y = Y0 + R0 * YUN(I)
	   CALL PLOT1Z_QPLOT(X, Y, 12)
 50	   CONTINUE
 60	   CONTINUE
	   Y = Y1
	   X = (R0**2 - (Y - Y0)**2)
	   IF( X.LT.0.0)X = 0.0
	   X = X0 + SQRT(X)
	   CALL PLOT1Z_QPLOT(X, Y, 12)
	   RETURN
C
C
 114	   CONTINUE
	   CALL PLOT1Z_QPLOT(X, Y, 13)
	   IF( I2.LE.I1)GO TO 160
	   DO 150 I = I1, I2, I3
	   CALL PLOT1Z_QPLOT(X, Y, 12)
	   X = X0 + R0 * XUN(I)
	   Y = Y0 + R0 * YUN(I)
	   CALL PLOT1Z_QPLOT(X, Y, 13)
 150	   CONTINUE
 160	   CONTINUE
	   Y = Y1
	   X = (R0**2 - (Y - Y0)**2)
	   IF( X.LT.0.0)X = 0.0
	   X = X0 + SQRT(X)
	   CALL PLOT1Z_QPLOT(X, Y, 13)
	   CALL PLOT1Z_QPLOT(X, Y, 12)
	   RETURN
C
C
 214	   CONTINUE
	   CALL PLOT1Z_QPLOT(X, Y, 13)
	   I4 = 12
	   I5 = 13
	   IF( I2.LE.I1)GO TO 260
	   DO 250 I = I1, I2, I3
	   X = X0 + R0 * XUN(I)
	   Y = Y0 + R0 * YUN(I)
	   CALL PLOT1Z_QPLOT(X, Y, I4)
	   I6 = I4
	   I4 = I5
	   I5 = I6
 250	   CONTINUE
 260	   CONTINUE
	   Y = Y1
	   X = (R0**2 - (Y - Y0)**2)
	   IF( X.LT.0.0)X = 0.0
	   X = X0 + SQRT(X)
	   CALL PLOT1Z_QPLOT(X, Y, 12)
	   RETURN
C
C
 314	   CONTINUE
	   CALL PLOT1Z_QPLOT(X, Y, 13)
	   K = -1
	   IF( I2.LE.I1)GO TO 360
	   DO 350 I = I1, I2, I3
	   K = -1 * K
	   X = X0 + (R0 + K * 0.05) * XUN(I)
	   Y = Y0 + (R0 + K * 0.05) * YUN(I)
	   CALL PLOT1Z_QPLOT(X, Y, 12)
 350	   CONTINUE
 360	   CONTINUE
	   Y = Y1
	   X = (R0**2 - (Y - Y0)**2)
	   IF( X.LT.0.0)X = 0.0
	   X = X0 + SQRT(X)
	   CALL PLOT1Z_QPLOT(X, Y, 12)
	   RETURN
C
C
414	   CONTINUE
	   CALL PLOT1Z_QPLOT(X, Y, 13)
	   IF( I2.LE.I1) GO TO 460
	   R3 = R0 * 0.9
	   DO 450 I = I1, I2, I3
	   X = X0 + R0 * XUN(I)
	   Y = Y0 + R0 * YUN(I)
	   CALL PLOT1Z_QPLOT(X, Y, 12)
	   X3 = X0 + R3 * XUN(I)
	   Y3 = Y0 + R3 * YUN(I)
	   CALL PLOT1Z_QPLOT(X3, Y3, 12)
	   CALL PLOT1Z_QPLOT(X, Y, 13)
 450 	   CONTINUE
 460 	   CONTINUE
	   Y = Y1
	   X = (R0**2 - (Y - Y0)**2)
	   IF( X.LT.0.0)X = 0.0
	   X = X0 + SQRT(X)
	   CALL PLOT1Z_QPLOT(X, Y, 12)
	   RETURN
C
C
 514 	   CONTINUE
	   CALL PLOT1Z_QPLOT(X, Y, 13)
	   IF( I2.LE.I1)GO TO 560
	   R4 = 0.15
	   DO 550 I = I1, I2, I3
	   X = X0 + R0 * XUN(I)
	   Y = Y0 + R0 * YUN(I)
	   CALL PLOT1Z_QPLOT(X, Y, 12)
	   X4 = X - R4 * XUN(I)
	   Y4 = Y - R4 * YUN(I)
	   CALL PLOT1Z_QPLOT(X4, Y4, 12)
	   CALL PLOT1Z_QPLOT(X, Y, 13)
 550 	   CONTINUE
 560 	   CONTINUE
	   Y = Y1
	   X = (R0**2 - (Y - Y0)**2)
	   IF( X.LT.0.0)X = 0.0
	   X = X0 + SQRT(X)
	   CALL PLOT1Z_QPLOT(X, Y, 12)
	   RETURN
C
C
614	   CONTINUE
	END IF
	RETURN
C
C
	END

C
C==========================================================================
C
	SUBROUTINE CIRDR1(X0, Y0, R0, ITYP)
	COMMON /ANGL/ PI, TWOPI, DT, XUN(126)
	COMMON /EDGE/TH1(50), TH2(50), NA(50), NT(50), 
	1  XN(50), YN(50), ZN(50), RN(50), NN
	DIMENSION YUN(101)
	EQUIVALENCE (YUN(1), XUN(26))
C	WRITE(5, 11)NN, X0, Y0, Z0, R0, (I, TH1(I), TH2(I), NA(I), NT(I), 
C	1    XN(I), YN(I), ZN(I), RN(I), I = 1, NN)
C   11	FORMAT(/, I4, 4F6.1, /, (4X, I4, 2F6.2, 2I4, 4F6.1))
	TMX = -PI
	T1 = -PI
	T2 = PI
	NN1 = NN + 1
	DO 10 J = 1, NN1
	   IF( J.GT.NN)GO TO 12
	   T1 = -PI
	   T2 = TH1(J)
	   IF( TH2(J).GT.TMX)TMX = TH2(J)
	   DO 20 I2 = 1, NN
C	      IF( I2.EQ.J)GO TO 20
	      T3 = TH2(I2)
	      T4 = TH1(I2)
	      IF( I2.GE.J)T3 = T3 - TWOPI
	      IF( T3.GT.T1)T1 = T3
 20	      CONTINUE
 12	   IF( T1.GT.T2)GO TO 14
	   I1 = IFIX(T1 / DT + 52.0)
	   I2 = IFIX(T2 / DT + 50.0)
	   I3 = (0.5 + 3.5 / R0)
	   IF( I3.LT.1)I3 = 1
	   IF( I3.GT.5)I3 = 5
	   X = X0 + R0 * SIN(T1)
	   Y = Y0 + R0 * COS(T1)
	   CALL PLOT1Z_QPLOT(X, Y,  + 13)
C	WRITE(5, 13)J, T1, X, Y
C  13	FORMAT(' ', I5, 3F7.2)
	   IF( I2.LE.I1)GO TO 40
	   DO 30 I = I1, I2, I3
	      X = X0 + R0 * XUN(I)
	      Y = Y0 + R0 * YUN(I)
	      CALL PLOT1Z_QPLOT(X, Y, 12)
 30	      CONTINUE
 40	   CONTINUE
	   Y = Y0 + R0 * COS(T2)
	   X = X0 + R0 * SIN(T2)
	   CALL PLOT1Z_QPLOT(X, Y, 12)
C	WRITE(5, 13)J, T2, X, Y
C
 14	   T1 = TMX
	   T2 = PI
 10	   CONTINUE
	RETURN
	END

C
C==========================================================================
C
	SUBROUTINE DRWEDG(X0, Y0, Z0, R0, STROX, STROY)
	COMMON /SEAM/ IQ, U0(3, 50), U1(3, 50), U2(3, 50), RU0(50), 
	1 PHII(50), PHIT(50)
	COMMON /ANGL/PI, TWOPI, DT, YUN(126)
	COMMON /EDGE/ TH1(50), TH2(50), NA(50), NT(50), 
	1  XN(50), YN(50), ZN(50), RN(50), NN
C
C	ELIMINATE COMPLETELY HIDDEN SEGMENTS
C
	DO 10 I = 1, 50
	   IF( NT(I).EQ.5)GO TO 10
	   DO 20 J = 1, 50
	      IF( I.EQ.J)GO TO 20
	      IF( NT(J).EQ.0)GO TO 20
	      IF( NT(J).EQ.5)GO TO 20
	      IF( TH1(J).LT.TH1(I))GO TO 22
	      IF( TH2(J).GT.TH2(I))GO TO 22
C
C	J INSIDE OF I
C
	      IF( NT(J).LT.NT(I))NT(J) = 0
	      GO TO 20
 22	      IF( TH2(I).LT.PI)GO TO 20
	      IF( TH1(J) + TWOPI.LT.TH1(I))GO TO 20
	      IF( TH2(J) + TWOPI.GT.TH2(I))GO TO 20
C
C	J INSIDE OF I
C
	      IF( NT(J).LT.NT(I))NT(J) = 0
 20	      CONTINUE
 10	   CONTINUE
C
C	CONDENCE
C
	DO 30 I = 1, 50
 32	   IF( I.GT.NN)GO TO 30
	   IF( NT(I).NE.0)GO TO 30
	   I1 = I + 1
	   NN = NN - 1
	   IF( I1.EQ.50)GO TO 30
	   DO 40 J = I1, 50
	      J1      = J - 1
	      TH1(J1) = TH1(J)
	      TH2(J1) = TH2(J)
	      NT(J1)  = NT(J)
	      NA(J1)  = NA(J)
	      XN(J1)  = XN(J)
	      YN(J1)  = YN(J)
	      ZN(J1)  = ZN(J)
	      RN(J1)  = RN(J)
 40	      CONTINUE
	   GO TO 32
 30	   CONTINUE
	IF( NN.EQ.1)GO TO 100
C
C	SORT
C
	NN1 = NN - 1
	DO 50 I = 1, NN1
	   I1 = I + 1
	   DO 60 J = I1, NN
	      IF( TH1(J).GT.TH1(I))GO TO 60
	      TT      = TH1(I)
	      TH1(I)  = TH1(J)
	      TH1(J)  = TT
	      TT      = TH2(I)
	      TH2(I)  = TH2(J)
	      TH2(J)  = TT
	      IT      = NT(I)
	      NT(I)   = NT(J)
	      NT(J)   = IT
	      IT      = NA(I)
	      NA(I)   = NA(J)
	      NA(J)   = IT
	      TT      = XN(I)
	      XN(I)   = XN(J)
	      XN(J)   = TT
	      TT      = YN(I)
	      YN(I)   = YN(J)
	      YN(J)   = TT
	      TT      = ZN(I)
	      ZN(I)   = ZN(J)
	      ZN(J)   = TT
	      TT      = RN(I)
	      RN(I)   = RN(J)
	      RN(J)   = TT
 60	      CONTINUE
 50	   CONTINUE
 100	CONTINUE
C
C	DRAW OUTER CIRCLE SEGMENTS
C
	ITYP = 1
	CALL CIRDR1(X0, Y0, R0, ITYP)
C
C	ESTABLISH LIMITS OF SEGMENT AND DRAW IT
C
C	
C	SET SEGMENT PARAMETERS
C
	IN = 0
	DO I = 1, NN
	   IF( NT(I).EQ.1)IN = IN + 1
	END DO
	IF( IN.EQ.0)RETURN
	DO 120 I = 1, 50
	   IF( I.GT.NN)GO TO 120
	   CALL INCIRC( X0, Y0, Z0, R0, XN(I), YN(I), ZN(I), RN(I), 
	1                 I, NT(I))
 120	   CONTINUE
	IQ = NN
	DO 140 I = 1, 50
	   IF( I.GT.IQ)GO TO 140
	   IF( NT(I).NE.1)GO TO 140
	   DO 160 J = 1, 50
	      IF( J.GT.IQ)GO TO 160
	      IF( J.EQ.I)GO TO 160
	      IF( NA(J).EQ.NA(I))GO TO 160
	      IF( NT(J).GE.4)GO TO 164
	      IF( TH2(J).LT.TH1(I).OR.TH1(J).GT.TH2(I))GO TO 162
	      CALL INTR( I, J, NT(J), XN(J), YN(J), ZN(J), RN(J), 
	1                STROX, STROY)
	      GO TO 160
 162	      CONTINUE
	      IF( (TH1(I) .GT. (TH2(J) - TWOPI) ) .AND.
	1         (TH2(I) .LT. (TH1(J) + TWOPI) ) ) GO TO 160
 164	      CALL INTR( I, J, NT(J), XN(J), YN(J), ZN(J), RN(J), 
 	1                 STROX, STROY)
 160	      CONTINUE
	   CALL DRAWIT(STROX, STROY, I)
 140	   CONTINUE
	RETURN
	END

C
C==========================================================================
C
	SUBROUTINE ANGTBL
	COMMON /ANGL/PI, TWOPI, DT, XUN(126)
	DT = PI / 50.
	DO I = 1, 126
	   XUN(I) = SIN(DT * FLOAT(I - 51))
	END DO
	RETURN
	END

C
C==========================================================================
C
	SUBROUTINE THSET(THI, THT, NTYP, NATM, X, Y, Z, R)
	COMMON /EDGE/TH1(50), TH2(50), NA(50), NT(50), 
	1  XN(50), YN(50), ZN(50), RN(50), NN
 50	CONTINUE
	IF( NN.GE.50)GO TO 100
	NN      = NN + 1
	TH1(NN) = THI
	TH2(NN) = THT
	NT(NN)  = NTYP
	NA(NN)  = NATM
	XN(NN)  = X
	YN(NN)  = Y
	ZN(NN)  = Z
	RN(NN)  = R
	RETURN
C
C	CONDENCE
C
 100	CONTINUE
	DO 10 I = 1, 50
	   DO 20 J = 1, 50
	      IF( I.EQ.J)GO TO 20
	      IF( NT(J).EQ.0)GO TO 20
	      IF( TH1(J).LT.TH1(I))GO TO 20
	      IF( TH2(J).GT.TH2(I))GO TO 20
C
C	J INSIDE OF I
C
	      IF( NT(J).LE.NT(I))NT(J) = 0
 20	      CONTINUE
 10	   CONTINUE
	DO 30 I = 1, 50
	   IF( I.GT.NN)GO TO 30
	   IF( NT(I).NE.0)GO TO 30
	   I1 = I + 1
	   NN = NN - 1
	   IF( I1.EQ.50)GO TO 30
	   DO 40 J = I1, 50
	      J1 = J - 1
	      TH1(J1) = TH1(J)
	      TH2(J1) = TH2(J)
	      NT(J1) = NT(J)
	      NA(J1) = NA(J)
	      XN(J1) = XN(J)
	      YN(J1) = YN(J)
	      ZN(J1) = ZN(J)
	      RN(J1) = RN(J)
 40	      CONTINUE
 30	   CONTINUE
	IF( NN.EQ.50)RETURN
	GO TO 50
	END

C
C==========================================================================
C
	SUBROUTINE INTR(II, IJ, ITYP, X, Y, Z, R, STROX, STROY)
	COMMON /SEAM/ IQ, U0(3, 50), U1(3, 50), U2(3, 50), RU0(50), 
	1 PHII(50), PHIT(50)
	COMMON /EDGE/ THI(50), THT(50), NA(50), NT(50), 
	1  XM(50), YM(50), ZM(50), RM(50), NN
	REAL*4 R1, R2, RI, RT, PHIN, PHI1, PHI2, R, RN, SN, CS
	PI = 3.141592654
	TWOPI = 2.0 * PI
	I = II
	IF( ITYP.EQ.1)GO TO 300
	IF( ITYP.EQ.2)GO TO 200
	IF( ITYP.EQ.3)GO TO 300
C	IF( ITYP.EQ.4)GO TO 500
C	IF( ITYP.EQ.5)GO TO 500
 500	CONTINUE
	R2 = (U0(1, I) - X)**2 + (U0(2, I) - Y)**2
	R2 = SQRT(R2)
D	WRITE(5, 17)I, IJ, PHII(I), PHIT(I), RU0(I), R, R2
	IF( (RU0(I) + R).LE.R2)GO TO 100
C	IF( (R2 + RU0(I)).GT.R)GO TO 10
C	GO TO 90
 10	CONTINUE
	SN  = SIN(PHII(I))
	CS  = COS(PHII(I))
	ZPI = U0(3, I) + RU0(I) * (CS * U1(3, I) + SN * U2(3, I))
	YPI = U0(2, I) + RU0(I) * (CS * U1(2, I) + SN * U2(2, I))
	YPI = YPI + ZPI * STROY
	XPI = U0(1, I) + RU0(I) * (CS * U1(1, I) + SN * U2(1, I))
	XPI = XPI + ZPI * STROX
	SN  = SIN(PHIT(I))
	CS  = COS(PHIT(I))
	ZPT = U0(3, I) + RU0(I) * (CS * U1(3, I) + SN * U2(3, I))
	YPT = U0(2, I) + RU0(I) * (CS * U1(2, I) + SN * U2(2, I))
	YPT = YPT + ZPT * STROY
	XPT = U0(1, I) + RU0(I) * (CS * U1(1, I) + SN * U2(1, I))
	XPT = XPT + ZPT * STROX
	RI = SQRT((XPI - X)**2 + (YPI - Y)**2)
	RT = SQRT((XPT - X)**2 + (YPT - Y)**2)
	IF( RI.GT.R.OR.RT.GT.R)GO TO 30
	GO TO 90
 30	CONTINUE
D	WRITE(5, 17) II, IJ, PHII(II), PHIT(II), XPI, YPI, ZPI, RI, 
D	1            XPT, YPT, ZPT, RT
D  17	FORMAT( ' INT1 ', 2I5, 10F7.2)
	IF( (RI.LE.R).OR.(RT.LE.R))GO TO 40
	ACS = ((X - U0(1, I)) * U1(1, I) + (Y - U0(2, I)) * U1(2, I))
	ASN = ((X - U0(1, I)) * U2(1, I) + (Y - U0(2, I)) * U2(2, I))
	A1 = SQRT(ASN**2 + ACS**2)
	IF( ASN.LT.0.0)A1 = -A1
	CS = ACS / A1
	SN = ASN / A1
	ZN = U0(3, I) + RU0(I) * (U1(3, I) * CS + U2(3, I) * SN)
	XN = U0(1, I) + RU0(I) * (U1(1, I) * CS + U2(1, I) * SN)
	XN = XN + ZN * STROX
	YN = U0(2, I) + RU0(I) * (U1(2, I) * CS + U2(2, I) * SN)
	YN = YN + ZN * STROY
	RN = SQRT((X - XN)**2 + (Y - YN)**2)
	IF( RN.LT.R)GO TO 34
	PHIM = ATAN2(SN, CS)
	IF( PHIM.LE.PHII(I) .OR. PHIM.GT.PHIT(I) )GO TO 100
	PHIN = ( PHIM + PHII(I) ) * 0.5
	CS   = COS(PHIN)
	SN   = SIN(PHIN)
	ZN   = U0(3, I) + RU0(I) * (U1(3, I) * CS + U2(3, I) * SN)
	XN   = U0(1, I) + RU0(I) * (U1(1, I) * CS + U2(1, I) * SN)
	XN   = XN + ZN * STROX
	YN   = U0(2, I) + RU0(I) * (U1(2, I) * CS + U2(2, I) * SN)
	YN   = YN + ZN * STROY
	RN   = SQRT((X - XN)**2 + (Y - YN)**2)
	IF( RN.LT.R)GO TO 34
	PHIN = ( PHIM + PHIT(I) ) * 0.5
	CS   = COS(PHIN)
	SN   = SIN(PHIN)
	ZN   = U0(3, I) + RU0(I) * (U1(3, I) * CS + U2(3, I) * SN)
	XN   = U0(1, I) + RU0(I) * (U1(1, I) * CS + U2(1, I) * SN)
	XN   = XN + ZN * STROX
	YN   = U0(2, I) + RU0(I) * (U1(2, I) * CS + U2(2, I) * SN)
	YN   = YN + ZN * STROY
	RN   = SQRT((X - XN)**2 + (Y - YN)**2)
	IF( RN.GE.R)GO TO 100
 34	IF( IQ.GE.50)GO TO 40
	IF( IQ.GE.50)GO TO 40
	PHIN = ATAN2(SN, CS)
	IF( PHIN.LE.PHII(I).OR.PHIN.GE.PHIT(I))GO TO 100
	IQ      = IQ + 1
	NT(IQ)  = NT(I)
	RU0(IQ) = RU0(I)
	THI(IQ) = THI(I)
	THT(IQ) = THT(I)
	XM(IQ)  = XM(I)
	YM(IQ)  = YM(I)
	ZM(IQ)  = ZM(I)
	RM(IQ)  = RM(I)
	NA(IQ)  = NA(I)
	DO K = 1, 3
	   U0(K, IQ) = U0(K, I)
	   U1(K, IQ) = U1(K, I)
	   U2(K, IQ) = U2(K, I)
	END DO
	PHIT(IQ) = PHIT(I)
	PHIT(I)  = PHIN - 0.005
	RT       = RN
	PHII(IQ) = PHIN + 0.005
D	WRITE(5, 17)II, IJ, PHII(I), PHIT(II)
D	WRITE(5, 17)II, IQ, PHII(IQ), PHIT(IQ)
 40	CONTINUE
	R1   = RI
	R2   = RT
	PHI1 = PHII(I)
	PHI2 = PHIT(I)
	IF( R1.LT.R2)GO TO 46
	RRR  = R1
	R1   = R2
	R2   = RRR
	RRR  = PHI1
	PHI1 = PHI2
	PHI2 = RRR
 46	CONTINUE
C	WRITE(5, 11)I, R1, R2, R, PHI1, PHI2
C  11	FORMAT(I5, 6F7.3)
	RN   = 0.0
	IK   = 0
C   44	PHIN = PHI1 + (PHI2 - PHI1) * (R - R1) / (R2 - R1)
 44	PHIN = (PHI1 + PHI2) * 0.5
	IF( IK.GT.8)GO TO 50
C	WRITE(5, 11)I, PHIN, PHI1, PHI2, RN, R1, R2
	SN = SIN(PHIN)
	CS = COS(PHIN)
	ZN = U0(3, I) + RU0(I) * (U1(3, I) * CS + U2(3, I) * SN)
	XN = U0(1, I) + RU0(I) * (U1(1, I) * CS + U2(1, I) * SN)
	XN = XN + ZN * STROX
	YN = U0(2, I) + RU0(I) * (U1(2, I) * CS + U2(2, I) * SN)
	YN = YN  + ZN * STROY
	RN = SQRT((X - XN)**2 + (Y - YN)**2)
	IF( RN.GT.R .AND. RN - R.LT.0.001)GO TO 50
	IK = IK + 1
	IF( RN.GT.R)GO TO 42
	R1 = RN
	PHI1 = PHIN
	GO TO 44
 42	R2 = RN
	PHI2 = PHIN
	GO TO 44
 50	CONTINUE
	IF( RT.LT.R)PHIT(I) = PHIN
	IF( RI.LT.R)PHII(I) = PHIN
D	WRITE(5, 17)I, IJ, PHII(I), PHIT(I), XN, YN, ZN, RN
C	WRITE(5, 11)I, PHII(I), PHIT(I), PHI1, PHI2
	GO TO 100
 90	CONTINUE
	RU0(I) = 0.0
	RETURN
 100	CONTINUE
	RETURN
 200	CONTINUE
C
C	EDGE SEGMENT INTERSECTION
C
	IF( RU0(I).EQ.0.0)RETURN
	KK  = 1
	IF( NA(IJ).LE.2)KK = 2
	IF( U1(KK, I).EQ.0.0)GO TO 250
	A   = -U2(KK, I) / U1(KK, I)
	IF( KK.EQ.1)EDG = XM(IJ)
	IF( KK.EQ.2)EDG = YM(IJ)
	B   = (EDG - U0(KK, I)) / (RU0(I) * U1(KK, I))
	SN  =  (A * B)**2 - (1.0 + A**2)*(B**2 - 1.0)
	IF( SN.LT.0.0) RETURN
	SN   = -A * B + SQRT(SN)
	SN   = SN / (1.0 + A**2)
	CS   = A * SN + B
	PHI  = ATAN2(SN, CS)
	END1 = U0(KK, I)  + 
	1  RU0(I) * (COS(PHII(I)) * U1(KK, I) + SIN(PHII(I)) * U2(KK, I))
	END2 = U0(KK, I)  + 
	1  RU0(I) * (COS(PHIT(I)) * U1(KK, I) + SIN(PHIT(I)) * U2(KK, I))
	IF( (NA(IJ).NE.2.AND.NA(IJ).NE.4))GO TO 210
	IF( EDG.LT.END1) PHII(I) = PHI
	IF( EDG.LT.END2) PHIT(I) = PHI
	RETURN
 210	CONTINUE
	IF( EDG.GT.END1) PHII(I) = PHI
	IF( EDG.GT.END2) PHIT(I) = PHI
	RETURN
 250	CONTINUE
	RETURN
 300	CONTINUE
C	
C	SEGMENT - SEGMENT INTERSECTION
C	
	R2 = (U0(1, I) - X)**2 + (U0(2, I) - Y)**2
	R2 = SQRT(R2)
D	WRITE(5, 19)I, IJ, PHII(I), PHIT(I), RU0(I), R, R2
	IF( (RU0(I) + R).LE.R2)GO TO 400
 310	CONTINUE
	SN = SIN( PHII(I))
	CS = COS( PHII(I))
	ZPI = U0(3, I) + RU0(I) * (CS * U1(3, I) + SN * U2(3, I))
	YPI = U0(2, I) + RU0(I) * (CS * U1(2, I) + SN * U2(2, I))
	YPI = YPI + ZPI * STROY
	XPI = U0(1, I) + RU0(I) * (CS * U1(1, I) + SN * U2(1, I))
	XPI = XPI + ZPI * STROX
	SN = SIN( PHIT(I))
	CS = COS( PHIT(I))
	ZPT = U0(3, I) + RU0(I) * (CS * U1(3, I) + SN * U2(3, I))
	YPT = U0(2, I) + RU0(I) * (CS * U1(2, I) + SN * U2(2, I))
	YPT = YPT + ZPT * STROY
	XPT = U0(1, I) + RU0(I) * (CS * U1(1, I) + SN * U2(1, I))
	XPT = XPT + ZPT * STROX
	RI = SQRT( (XPI - X)**2 + (YPI - Y)**2 + (ZPI - Z)**2)
	RT = SQRT( (XPT - X)**2 + (YPT - Y)**2 + (ZPT - Z)**2)
	IF( RI.GT.R.OR.RT.GT.R)GO TO 330
	PHIN = (PHII(I) + PHIT(I)) * 0.50
	SN = SIN(PHIN)
	CS = COS(PHIN)
	ZN = U0(3, I) + RU0(I) * (CS * U1(3, I) + SN * U2(3, I))
	YN = U0(2, I) + RU0(I) * (CS * U1(2, I) + SN * U2(2, I))
	YN = YN + ZN * STROY
	XN = U0(1, I) + RU0(I) * (CS * U1(1, I) + SN * U2(1, I))
	XN = XN + ZN * STROX
	RN = SQRT((X - XN)**2 + (Y - YN)**2 + (Z - ZN)**2)
	IF( RN.GT.R)GO TO 316
	PHIN = PHII(I) + (PHIT(I) - PHII(I)) * 0.25
	SN = SIN(PHIN)
	CS = COS(PHIN)
	ZN = U0(3, I) + RU0(I) * (CS * U1(3, I) + SN * U2(3, I))
	YN = U0(2, I) + RU0(I) * (CS * U1(2, I) + SN * U2(2, I))
	YN = YN + ZN * STROY
	XN = U0(1, I) + RU0(I) * (CS * U1(1, I) + SN * U2(1, I))
	XN = XN + ZN * STROX
	RN = SQRT((X - XN)**2 + (Y - YN)**2 + (Z - ZN)**2)
	IF( RN.GT.R)GO TO 316
	PHIN = PHII(I) +  (PHIT(I) - PHII(I)) * 0.75
	SN = SIN(PHIN)
	CS = COS(PHIN)
	ZN = U0(3, I) + RU0(I) * (CS * U1(3, I) + SN * U2(3, I))
	YN = U0(2, I) + RU0(I) * (CS * U1(2, I) + SN * U2(2, I))
	YN = YN + ZN * STROY
	XN = U0(1, I) + RU0(I) * (CS * U1(1, I) + SN * U2(1, I))
	XN = XN + ZN * STROX
	RN = SQRT((X - XN)**2 + (Y - YN)**2 + (Z - ZN)**2)
	IF( RN.LE.R)GO TO 390
 316	III = 3
	R3 = RN
	PHI3 = PHIN
	R1 = RI
	R2 = R3
	PHI1 = PHII(I)
	PHI2 = PHI3
	GO TO 346
 330	CONTINUE
D	WRITE(5, 19)II, IJ, PHII(II), PHIT(II), XPI, YPI, ZPI, RI, XPT, YPT, ZPT, RT
D  19	FORMAT(' INT2 ', 2I5, 10F7.2)
	IF( (RI.LE.R).OR.(RT.LE.R))GO TO 340
C	ACS = ((X - U0(1, I)) * U1(1, I) + (Y - U0(2, I)) * U1(2, I))
C	ASN = ((X - U0(1, I)) * U2(1, I) + (Y - U0(2, I)) * U2(2, I))
C	A1 = SQRT(ASN**2 + ACS**2)
C	IF( ASN.LT.0.0)A1 = -A1
C	CS = ACS / A1
C	SN = ASN / A1
	PHIM = (PHII(I) + PHIT(I)) * 0.50
	CS   = COS(PHIM)
	SN   = SIN(PHIM)
	ZN   = U0(3, I) + RU0(I) * (U1(3, I) * CS + U2(3, I) * SN)
	XN   = U0(1, I) + RU0(I) * (U1(1, I) * CS + U2(1, I) * SN)
	XN   = XN + ZN * STROX
	YN   = U0(2, I) + RU0(I) * (U1(2, I) * CS + U2(2, I) * SN)
	YN   = YN + ZN * STROY
	RN   = SQRT((X - XN)**2 + (Y - YN)**2 + (Z - ZN)**2)
	IF( ZN.LT.Z)GO TO 500
	IF( RN.LT.R)GO TO 320
C	PHIM = ATAN2(SN, CS)
C	IF( PHIM.LE.PHII(I) .OR. PHIM.GT.PHIT(I) )GO TO 400
	PHIN = ( PHIM + PHII(I) ) * 0.5
	CS = COS(PHIN)
	SN = SIN(PHIN)
	ZN = U0(3, I) + RU0(I) * (U1(3, I) * CS + U2(3, I) * SN)
	XN = U0(1, I) + RU0(I) * (U1(1, I) * CS + U2(1, I) * SN)
	XN = XN + ZN * STROX
	YN = U0(2, I) + RU0(I) * (U1(2, I) * CS + U2(2, I) * SN)
	YN = YN + ZN * STROY
	RN = SQRT((X - XN)**2 + (Y - YN)**2 + (Z - ZN)**2)
	IF( RN.LT.R)GO TO 320
	IF( ZN.LT.Z)GO TO 500
	PHIN = ( PHIM + PHIT(I) ) * 0.5
	CS = COS(PHIN)
	SN = SIN(PHIN)
	ZN = U0(3, I) + RU0(I) * (U1(3, I) * CS + U2(3, I) * SN)
	XN = U0(1, I) + RU0(I) * (U1(1, I) * CS + U2(1, I) * SN)
	XN = XN + ZN * STROX
	YN = U0(2, I) + RU0(I) * (U1(2, I) * CS + U2(2, I) * SN)
	YN = YN + ZN * STROY
	RN = SQRT((X - XN)**2 + (Y - YN)**2 + (Z - ZN)**2)
	IF( RN.LT.R)GO TO 320
	IF( ZN.LT.Z)GO TO 500
	GO TO 400
 320	CONTINUE
	IF( IQ.GE.50)GO TO 340
	PHIN = ATAN2(SN, CS)
	IF( PHIN.LE.PHII(I).OR.PHIN.GE.PHIT(I))GO TO 400
	IQ = IQ + 1
	NT(IQ) = NT(I)
	RU0(IQ) = RU0(I)
	THI(IQ) = THI(I)
	THT(IQ) = THT(I)
	XM(IQ) = XM(I)
	YM(IQ) = YM(I)
	ZM(IQ) = ZM(I)
	RM(IQ) = RM(I)
	NA(IQ) = NA(I)
	DO K = 1, 3
	   U0(K, IQ) = U0(K, I)
	   U1(K, IQ) = U1(K, I)
	   U2(K, IQ) = U2(K, I)
	END DO
	PHIT(IQ) = PHIT(I)
	PHIT(I) = PHIN - 0.005
	RT = RN
	PHII(IQ) = PHIN + 0.005
D	WRITE(5, 19)II, IJ, PHII(I), PHIT(II)
D	WRITE(5, 19)II, IQ, PHII(IQ), PHIT(IQ)
 340	CONTINUE
	R1 = RI
	R2 = RT
	PHI1 = PHII(I)
	PHI2 = PHIT(I)
	III = 2
	IF( R1.LT.R2)GO TO 346
	III = 1
	RRR = R1
	R1 = R2
	R2 = RRR
	RRR = PHI1
	PHI1 = PHI2
	PHI2 = RRR
 346	CONTINUE
C	WRITE(5, 11)I, R1, R2, R, PHI1, PHI2
C  11	FORMAT(I5, 6F7.3)
	RN = 0.0
	IK = 0
C   44	PHIN = PHI1 + (PHI2 - PHI1) * (R - R1) / (R2 - R1)
 344	PHIN = (PHI1 + PHI2) * 0.5
	IF( IK.GT.8)GO TO 350
C	WRITE(5, 11)I, PHIN, PHI1, PHI2, RN, R1, R2
	SN = SIN(PHIN)
	CS = COS(PHIN)
	ZN = U0(3, I) + RU0(I) * (U1(3, I) * CS + U2(3, I) * SN)
	XN = U0(1, I) + RU0(I) * (U1(1, I) * CS + U2(1, I) * SN)
	XN = XN + ZN * STROX
	YN = U0(2, I) + RU0(I) * (U1(2, I) * CS + U2(2, I) * SN)
	YN = YN + ZN * STROY
	RN = SQRT((X - XN)**2 + (Y - YN)**2 + (Z - ZN)**2)
	IF( ABS(RN - R).LT.0.001)GO TO 350
	IK = IK + 1
	IF( RN.GT.R)GO TO 342
	R1 = RN
	PHI1 = PHIN
	GO TO 344
 342	R2 = RN
	PHI2 = PHIN
	GO TO 344
 350	CONTINUE
	IF( III.EQ.3)GO TO 352
	IF( III.EQ.1) PHIT(I) = PHIN
	IF( III.EQ.2) PHII(I) = PHIN
D	WRITE(5, 19) I, IJ, PHII(I), PHIT(I), XN, YN, ZN, RN, Z
	IF( ZN.LT.Z)GO TO 500
C	WRITE(5, 11) I, PHII(I), PHIT(I), PHI1, PHI2
	GO TO 400
 352	CONTINUE
	PHII(I) = PHIN
	RI = RN
	III = 1
	R1 = RT
	R2 = R3
	PHI2 = PHI3
	PHI1 = PHIT(I)
D	WRITE(5, 19)I, IJ, PHII(I), PHIT(I), XN, YN, ZN, RN, Z
	GO TO 346
 390	CONTINUE
	RU0(I) = 0.0
	RETURN
 400	CONTINUE
	RETURN
	END

C
C==========================================================================
C
	SUBROUTINE BOUNDS(XMN, XMX, YMN, YMX, II)
	COMMON /SEAM/ IQ, IID(10), U0(3, 10), U1(3, 10), U2(3, 10), 
	1             RU0(10), PHII(10), PHIT(10)
	PI = 3.141592654
	TWOPI = PI * 2.0
	I = II
	IIQ = 0
	IF( (U0(1, I) - RU0(I).GE.XMN) .AND.
	1   (U0(1, I) + RU0(I).LE.XMX) )GO TO 40
	IF( U2(1, I).EQ.0.0)GO TO 2
	IT = 1
	XMM = XMN
	IF( U0(1, I) + RU0(I).GT.XMX)IT = 2
	IF( IT.EQ.2)XMM = XMX
	G1 = (XMM - U0(1, I)) / (RU0(I) * U2(1, I))
	G2 = -U1(1, I) / U2(1, I)
	A = 1.0 + G2**2
	B = 2.0 * G1 * G2
	C = G1**2 - 1.0
	G1 = B**2 - 4.0 * A * C
	IF( G1.LT.0.0)GO TO 40
	G1 = SQRT(G1)
	CS1 = (-B + G1) / (2.0 * A)
	PH1 = 999.
	IF( ABS(CS1).LE.1.0 ) PH1 = ACOS(CS1)
	IF( PH1.LT.0.0 ) PH1 = PH1 + TWOPI
 	CS1 = (-B - G1) / (2.0 * A)
	PH2 = 999.
	IF( ABS(CS1).LE.1.0)PH2 = ACOS(CS1)
	IF( PH2.LT.0.0)PH2 = PH2 + TWOPI
	IF( PH1.GT.100.0.AND.PH2.GT.100.0)GO TO 40
	IF( PH2.GT.PH1)GO TO 4
	P = PH1
	PH1 = PH2
	PH2 = P
 4	CONTINUE
	IF( (IT.EQ.1).AND.(U1(1, I).GT.0.0))GO TO 20
	IF( (IT.EQ.2).AND.(U1(1, I).LT.0.0))GO TO 20
	P = PH1
	PH1 = PH2
	PH2 = P
 20	CONTINUE
	IF( PHIT(I).LT.PH1.AND.PHII(I) + TWOPI.GT.PH2)GO TO 40
	IF( PHIT(I).GT.PH2.OR.PHII(I).LT.PH1)GO TO 6
	IID(I) = 0
	GO TO 40
 6	IF( PHII(I).LT.PH1.AND.PHIT(I).GT.PH2)GO TO 8
	IF( PHII(I).LT.PH1)PHIT(I) = PH2
	IF( PHIT(I).GT.PH2)PHII(I) = PH1
	GO TO 40
 8	IF( IQ.GE.10)GO TO 40
	IQ = IQ + 1
	DO 10 K = 1, 3
	U0(K, IQ) = U0(K, I)
	U1(K, IQ) = U1(K, I)
	U2(K, IQ) = U2(K, I)
 10	CONTINUE
	IID(IQ) = IID(I)
	RU0(IQ) = RU0(I)
	PHIT(IQ) = PHIT(I)
	PHII(IQ) = PH2
	PHIT(I) = PH1
	IIQ = IIQ + 1
	GO TO 40
 2	CONTINUE
	GO TO 40
 40	IF( (U0(2, I) - RU0(I).GT.YMN) .AND.
 	1   (U0(2, I) + RU0(I).LT.YMX) )GO TO 80
	IT = 1
	YMM = YMN
	IF( U0(2, I) + RU0(I).GT.YMX)IT = 2
	IF( IT.EQ.2)YMM = YMX
	IF( U2(2, I).EQ.0.0)GO TO 42
	G1 = (YMM - U0(2, I)) / (RU0(I) * U2(2, I))
	G2 = -U1(2, I) / U2(2, I)
	A = 1.0 + G2**2
	B = 2.0 * G1 * G2
	C = G1**2 - 1.0
	G1 = B**2 - 4.0 * A * C
	IF( G1.LT.0.0)GO TO 80
	G1 = SQRT(G1)
	CS1 = (-B + G1) / (2.0 * A)
	PH1 = 999.
	IF( ABS(CS1).LE.1.0)PH1 = ACOS(CS1)
	IF( PH1.LT.0.0)PH1 = PH1 + TWOPI
	PH2 = 999.
	CS1 = (-B - G1) / (2.0 * A)
	IF( ABS(CS1).LE.1.0)PH2 = ACOS(CS1)
	IF( PH2.LT.0.0)PH2 = PH2 + TWOPI
	IF( PH1.GT.100.0.AND.PH2.GT.100.0)GO TO 80
	IF( PH2.GT.PH1)GO TO 44
	P = PH1
	PH1 = PH2
	PH2 = P
 44	CONTINUE
	IF( (IT.EQ.1).AND.(U1(2, I).GT.0.0))GO TO 60
	IF( (IT.EQ.2).AND.(U1(2, I).LT.0.0))GO TO 60
	P = PH1
	PH1 = PH2
	PH2 = P
 60	CONTINUE
	IF( PHIT(I).LT.PH1.AND.PHII(I) + TWOPI.GT.PH2)GO TO 80
	IF( PHIT(I).GT.PH2.OR.PHII(I).LT.PH1)GO TO 46
	IID(I) = 0
	GO TO 80
 46	IF( PHII(I).LT.PH1.AND.PHIT(I).GT.PH2)GO TO 48
	IF( PHII(I).LT.PH1)PHIT(I) = PH2
	IF( PHIT(I).GT.PH2)PHII(I) = PH1
	GO TO 80
 48	IF( IQ.GE.10)GO TO 80
	IQ = IQ + 1
	DO 50 K = 1, 3
	U0(K, IQ) = U0(K, I)
	U1(K, IQ) = U1(K, I)
	U2(K, IQ) = U2(K, I)
 50	CONTINUE
	IID(IQ) = IID(I)
	RU0(IQ) = RU0(I)
	PHIT(IQ) = PHIT(I)
	PHII(IQ) = PH2
	PHIT(I) = PH1
	IIQ = IIQ + 1
	GO TO 80
 42	CONTINUE
 80	CONTINUE
	IF( IIQ.EQ.0)RETURN
	IIQ = IIQ - 1
	I = IQ - IIQ
	GO TO 40
	END

C
C==========================================================================
C
	SUBROUTINE DRAWIT(STROX, STROY, I)
	COMMON /SEAM/ IQ, U0(3, 50), U1(3, 50), U2(3, 50), RU0(50), 
	1 PHII(50), PHIT(50)
	COMMON /ANGL/PI, TWOPI, DPH, XUN(126)
	DIMENSION YUN(101)
	EQUIVALENCE (XUN(26), YUN(1))
	IF( RU0(I).EQ.0.0)RETURN
	IF( PHII(I).GE.PHIT(I))RETURN
D	WRITE(5, 11) I, RU0(I), (U0(K, I), K = 1, 3), PHII(I), PHIT(I)
D  11	FORMAT(5X, I5, 6F8.3)
	SN = SIN(PHII(I))
	CS = COS(PHII(I))
	ZP = U0(3, I) + RU0(I) * (CS * U1(3, I) + SN * U2(3, I))
	XP = U0(1, I) + RU0(I) * (CS * U1(1, I) + SN * U2(1, I))
	XP = XP + STROX * ZP
	YP = U0(2, I) + RU0(I) * (CS * U1(2, I) + SN * U2(2, I))
	YP = YP + STROY * ZP
	CALL PLOT1Z_QPLOT(XP, YP, 13)
D	WRITE(5, 13)I, PHII(I), XP, YP
D  13	FORMAT(' ', I5, 3F7.2)
	IPH1 = IFIX(PHII(I) / DPH + 52.0)
	IPH2 = IFIX(PHIT(I) / DPH + 50.0)
	IF( IPH2.LE.IPH1)GO TO 20
	DO 10 IIP = IPH1, IPH2
	IP = IIP
	IF( IP.GT.101)IP = IP - 100
	ZP = U0(3, I) + RU0(I) * ( YUN(IP) * U1(3, I) + XUN(IP) * U2(3, I) )
	YP = U0(2, I) + RU0(I)*( YUN(IP) * U1(2, I) + XUN(IP) * U2(2, I) )
	YP = YP + STROY * ZP
	XP = U0(1, I) + RU0(I) * ( XUN(IP) * U2(1, I) + YUN(IP) * U1(1, I) )
	XP = XP + STROX * ZP
	CALL PLOT1Z_QPLOT(XP, YP, 12)
 10	CONTINUE
 20	CONTINUE
	SN = SIN(PHIT(I))
	CS = COS(PHIT(I))
	ZP = U0(3, I) + RU0(I) * (CS * U1(3, I) + SN * U2(3, I))
	XP = U0(1, I) + RU0(I) * (CS * U1(1, I) + SN * U2(1, I))
	XP = XP + STROX * ZP
	YP = U0(2, I) + RU0(I) * (CS * U1(2, I) + SN * U2(2, I))
	YP = YP + STROY * ZP
	CALL PLOT1Z_QPLOT(XP, YP, 12)
D	WRITE(5, 13)I, PHIT(I), XP, YP
	RU0(I) = 0.0
	RETURN
	END

C
C==========================================================================
C
	SUBROUTINE INTRSC(X1, Y1, R1, X2, Y2, R2, X3, Y3, X4, Y4)
	IF( ABS(Y1 - Y2).EQ.0.00)GO TO 30
	A2 = 0.5 * (R1**2 - X1**2 - Y1**2 - R2**2 + X2**2 + Y2**2)
	A2 = A2 / (Y1 - Y2)
	A1 = -(X1 - X2) / (Y1 - Y2)
	A = 1 + A1**2
	B = -2.0 * (X1 + A1 * (A2 + Y1))
	C = X1**2 + (A2 + Y1)**2 - R1**2
	D = B**2 - 4.0 * A * C
	IF( D.LT.0.0)GO TO 10
	D = SQRT(D)
	X3 = (-B - D) / (2.0 * A)
	X4 = (-B + D) / (2.0 * A)
	Y3 = X3 * A1 - A2
	Y4 = X4 * A1 - A2
	RETURN
 10	X3 = 0.0
	Y3 = 0.0
	X4 = 0.0
	Y4 = 0.0
	RETURN
 30	CONTINUE
	X3 = R2**2 - R1**2 - X2**2 + X1**2
	X3 = 10000.0
	IF( (X1 - X2).NE.0.0 ) X3 = 0.5 * X3 / (X1 - X2)
	X4 = X3
	A = (R1**2 - (X3 - X1)**2)
	IF( A.LT.0.0)GO TO 10
	A = SQRT(A)
	IF( X2.LT.X0)GO TO 40
	Y3 = Y1 - A
	Y4 = Y1 + A
	RETURN
 40	Y3 = Y1 + A
	Y4 = Y1 - A
	RETURN
	END

C
C==========================================================================
C
	SUBROUTINE INCIRC(X1, Y1, Z1, R1, X2, Y2, Z2, R2, II, ITYP)
	COMMON /SEAM/ IQ, U0(3, 50), U1(3, 50), U2(3, 50), RU0(50), 
	1 PHII(50), PHIT(50)
	COMMON /EDGE/ THI(50), THT(50), NA(50), NT(50), 
	1  XN(50), YN(50), ZN(50), RN(50), NN
	DIMENSION V1(3), V2(3), V5(3)
C
	PI = 3.141592654
	TWOPI = PI * 2.0
	IF( ITYP.EQ.4)GO TO 100
	IF( ITYP.EQ.2)GO TO 300
D	WRITE(5, 11) X1, Y1, Z1, R1, X2, Y2, Z2, R2, II, ITYP
D 11	FORMAT(' ', 8F7.2, 2I4)
	V1(1) = X1
	V1(2) = Y1
	V1(3) = Z1
	V2(1) = X2
	V2(2) = Y2
	V2(3) = Z2
	R32 = ( V1(1) - V2(1))**2 + ( V1(2) - V2(2))**2 + 
	1     ( V1(3) - V2(3))**2 + 0.000001
	RATIO = (R2**2 - R1**2 + R32) / (2.0 * R32)
	R02 = R2**2 - R32 * RATIO**2
	RU0(II) = 0.0
D	WRITE(5, 11) V1, V2, R02, R32, II, ITYP
	IF( R02.LT.0.0)RETURN
	RU0(II) = SQRT(R02)
	R3 = SQRT(R32)
	DO 10 K = 1, 3
	   U0(K, II) = V2(K) + (V1(K) - V2(K)) * RATIO
	   V5(K) = (V2(K) - V1(K)) / R3
 10	CONTINUE
	IF( V5(1).NE.0.0.OR.V5(2).NE.0.0)GO TO 30
	U1(3, II) = 0.0
	U1(1, II) = 0.0
	U1(2, II) = 1.0
	GO TO 40
 30	CONTINUE
	U1(3, II) = 0.0
	U1(2, II) =  SQRT(V5(1)**2 / (V5(1)**2 + V5(2)**2))
	U1(1, II) = 1.0
	IF( V5(1).NE.0.0)U1(1, II) = -U1(2, II) * V5(2) / V5(1)
 40	CALL V1XV2(V5, U1(1, II), U2(1, II))
	IF( U2(3, II).GE.0.0)GO TO 50
	DO K = 1, 3
	   U1(K, II) = -U1(K, II)
	   U2(K, II) = -U2(K, II)
	END DO
 50	CONTINUE
	IF( (THT(II) - THI(II) + 0.001).LT.TWOPI)GO TO 52
	PHII(II) = -0.5 * PI
	PHIT(II) =  1.5 * PI
	GO TO 90
 52	T2 = 0.0
	T3 = 0.0
	T4 = 0.0
	T1 = 0.0
	V1(1) = X1 + R1 * SIN(THI(II))
	V1(2) = Y1 + R1 * COS(THI(II))
	V1(3) = Z1
	V2(1) = X1 + R1 * SIN(THT(II))
	V2(2) = Y1 + R1 * COS(THT(II))
	V2(3) = Z1
	DO K = 1, 3
	   T1 = T1 + U2(K, II) * (V1(K) - U0(K, II))
	   T2 = T2 + U1(K, II) * (V1(K) - U0(K, II))
	   T3 = T3 + U2(K, II) * (V2(K) - U0(K, II))
	   T4 = T4 + U1(K, II) * (V2(K) - U0(K, II))
	END DO
	PHII(II) = ATAN2(T1, T2)
	PHIT(II) = ATAN2(T3, T4)
	IF( PHIT(II).GT.PHII(II))GO TO 70
	P = PHIT(II)
	PHIT(II) = PHII(II)
	PHII(II) = P
 70	CONTINUE
C	IF( (ITYP.EQ.1) .AND. (PHIT(II) - PHII(II).LT.PI) ) GO TO 90
C	IF( (PHIT(II) - PHII(II).LT.PI))GO TO 90
	IF( (PHII(II).GT.-PI * 0.5).OR.(PHIT(II).LT.-PI * 0.5))GO TO 90
C	IF( (ITYP.EQ.3) .AND. (PHIT(II) - PHII(II).GE.PI) ) GO TO 90
	P = PHII(II)
	PHII(II) = PHIT(II)
	PHIT(II) = P + TWOPI
 90	CONTINUE
D	WRITE(5, 13)II, (U0(K, II), K = 1, 3), RU0(II), PHII(II), PHIT(II)
D 13	FORMAT('  **', I3, 3F7.1, 3F7.3)
	RETURN
 100	CONTINUE
C
C	SET PARAMETERS FOR CIRCLE OVERLAP
C
	RETURN
 200	CONTINUE
C	
C	SET PARAMETER FOR EDGE OVERLAP
C
	RETURN
 300	CONTINUE
C
C	SET PARAMETERS FOR EDGE OVERLAP
C
	RETURN
	END

C
C==========================================================================
C
	SUBROUTINE ATMGET(N, NN, NAME, ISF, X, Y, Z, B, IO, IFOK)
C	****	SUBROUTINE ATOMGET
C	GETS ATOM N FROM ATOMS FILE
C	***
	REAL*8 NAME
	IF( N.GT.0)GO TO 10
	   REWIND IO
	   N1   = 1
	   IFOK = 1
	   RETURN
 10	CONTINUE
	IF( N.GE.N1)GO TO 30
	   REWIND IO
	   N1 = 1
 30	CONTINUE
	IF( N1.EQ.N)GO TO 40
	   NI = N1
	   NF = N - 1
	   DO 50 N1 = NI, NF
	      READ(IO, 51, END = 99)
 51	      FORMAT(5X)
 50	   CONTINUE
 40	CONTINUE
	READ(IO, 11, END = 99)NN, NAME, ISF, X, Y, Z, B
 11	FORMAT(I5, 2X, A8, I3, 4F10.4)
	N1   = N + 1
	IFOK = 1
	RETURN
 99	IFOK = 0
	REWIND IO
	N1 = 1
	RETURN
	END

C
C==========================================================================
C
	SUBROUTINE ORTHIN(M)
	REAL*4 MO(3, 3), M(3, 3), XYZ(3), XYZO(3)
	DO 2 K1 = 1, 3
	DO 4 K2 = 1, 3
	      MO(K1, K2) = M(K1, K2)
 4	CONTINUE
 2	CONTINUE
	RETURN
	ENTRY ORTHO( X, Y, Z, XO, YO, ZO)
	XYZ(1) = X
	XYZ(2) = Y
	XYZ(3) = Z
	DO K1 = 1, 3
	   XYZO(K1) = 0.0
	   DO K2 = 1, 3
	      XYZO(K1) = XYZO(K1) + MO(K1, K2) * XYZ(K2)
	   END DO
	END DO
	XO = XYZO(1)
	YO = XYZO(2)
	ZO = XYZO(3)
	RETURN
	END

C
C==========================================================================
C
	SUBROUTINE FILL( NP, NP1, I, NPC)
C	
C	SUBROUTINE TO FILL IN BETWEEN NUMBERS ON BOTH SIDES OF THE
C	HYPHEN ON LINE FIVE OF THE DATA IN VEWCC.DAT.
C	
	INTEGER NP(50), NP1(16)
C
	L = NPC - 1
C	
C	THE FOLLOWING IS A LOOP TO FILL
C	IN NUMBERS BETWEEN NP(L) AND NP(I)
C
	DO K = L, 49
	   NP(K + 1) = NP(K) + 1
	   IF( NP(K + 1).EQ.ABS(NP1(I)))GO TO 20
C
C	IF NP(K + 1) EQUALS THE LAST NUMBER WHICH IS
C	NP1(I), CONTROL RETURNS BACK TO ORGDIR.
C
	   NPC = NPC + 1
	END DO
 20	RETURN
	END

C
C==========================================================================
C
	SUBROUTINE PLANE(N, X, VM, D, G)
C
C****FIT A LEAST-SQUARES PLANE TO A SET OF POINTS
C    PROCEDURE OF SCHOMAKER ET AL.(ACTA CRYST., 12, 600, 1959)
C
	DIMENSION G(3, 3), X(3, 20), XS(3), XXS(3, 3), ADJ(3, 3), 
	1         B(3, 3), A(3, 3), VM(3), VMI(3), BV(3)
	DATA ZIP/1.0E-5/
	DATA MM/10/
C
C----SET UP THE A MATRIX
	SN = N
	DO I = 1, 3
	   XS(I) = 0
	   DO K = 1, N
	      XS(I) = XS(I) + X(I, K)
	   END DO
	END DO
	DO I = 1, 3
	   DO J = 1, 3
	      XXS(I, J) = 0
	      DO K = 1, N
	         XXS(I, J) = XXS(I, J) + X(I, K) * X(J, K)
	      END DO
	      A(I, J) = XXS(I, J) - XS(I) * XS(J) / SN
	   END DO
	END DO
C
C----EVALUATE MATRIX B = ADJ(A) * G
 65	ADJ(1, 1) = A(2, 2) * A(3, 3) - A(2, 3) * A(3, 2)
	ADJ(2, 1) = A(3, 1) * A(2, 3) - A(2, 1) * A(3, 3)
	ADJ(3, 1) = A(2, 1) * A(3, 2) - A(3, 1) * A(2, 2)
	ADJ(1, 2) = A(3, 2) * A(1, 3) - A(1, 2) * A(3, 3)
	ADJ(2, 2) = A(1, 1) * A(3, 3) - A(3, 1) * A(1, 3)
	ADJ(3, 2) = A(3, 1) * A(1, 2) - A(1, 1) * A(3, 2)
	ADJ(1, 3) = A(1, 2) * A(2, 3) - A(1, 3) * A(2, 2)
	ADJ(2, 3) = A(2, 1) * A(1, 3) - A(1, 1) * A(2, 3)
	ADJ(3, 3) = A(1, 1) * A(2, 2) - A(2, 1) * A(1, 2)
	DO I = 1, 3
	   DO J = 1, 3
	      B(I, J) = 0
	      DO K = 1, 3
	         B(I, J) = B(I, J) + ADJ(I, K) * G(K, J)
	      END DO
	   END DO
	END DO
C
C----CHOOSE THE LARGEST COLUMN VECTOR OF B AS THE INITIAL SOLUTION
	BV(1) = B(1, 1)**2 + B(2, 1)**2 + B(3, 1)**2
	BV(2) = B(1, 2)**2 + B(2, 2)**2 + B(3, 2)**2
	BV(3) = B(1, 3)**2 + B(2, 3)**2 + B(3, 3)**2
	KK = 1
	IF( BV(2).GT.BV( 1)) KK = 2
	IF( BV(3).GT.BV(KK)) KK = 3
	DO I = 1, 3
	   VMI(I) = B(I, KK)
	END DO
C
C----SOLVE TO CONVERGENCE BY ITERATION OF M(I) = B * M(I-1)
	DO NNN = 1, MM
	   VM(1) = B(1, 1) * VMI(1) + B(1, 2) * VMI(2) + B(1, 3) * VMI(3)
	   VM(2) = B(2, 1) * VMI(1) + B(2, 2) * VMI(2) + B(2, 3) * VMI(3)
	   VM(3) = B(3, 1) * VMI(1) + B(3, 2) * VMI(2) + B(3, 3) * VMI(3)
	   IF( VMI(1).LT.1.0E-5 ) VMI(1) = 1.0E-5
	   RATIO1 = VM(1) / VMI(1)
	   IF( VMI(2).LT.1.0E-5 ) VMI(2) = 1.0E-5
	   RATIO2 = VM(2) / VMI(2)
	   IF( VMI(3).LT.1.0E-5 ) VMI(3) = 1.0E-5
	   RATIO3 = VM(3) / VMI(3)
	   IF( RATIO1.LT.1.0E-5 ) RATIO1 = 1.0E-5
	   RAT12 = ABS( RATIO2 / RATIO1 - 1.0 )
	   RAT13 = ABS( RATIO3 / RATIO1 - 1.0 )
	   IF( RAT12.LT.ZIP .AND. RAT13.LT.ZIP )GO TO 110
	   DO I = 1, 3
	      IF( VM(1).GT.1.0E-5 ) VMI(I) = VM(I) / VM(1)
	   END DO
	END DO
C
C----NORMALIZE THE SOLUTION VECTOR AND EVALUATE D(PLANE TO ORIGIN DISTAN
 110	ORM = 0
	DO I = 1, 3
	   DO J = 1, 3
	      ORM = ORM + VM(I) * VM(J) * G(I, J)
	   END DO
	END DO
	ORM = SQRT(ORM)
	DO I = 1, 3
	   VM(I) = VM(I) / ORM
	END DO
	D = (VM(1) * XS(1) + VM(2) * XS(2) + VM(3) * XS(3)) / SN
	RETURN
	END

C
C==========================================================================
C
	SUBROUTINE LINEV1( XA, YA, XB, YB, XMN, XMX, YMN, YMX)
	X1 = XA
	Y1 = YA
	X2 = XB
	Y2 = YB
	IF( X1.GE.XMN .AND. X1.LE.XMX .AND. Y1.GE.YMN .AND. Y1.LE.YMX .AND.
	1   X2.GE.XMN .AND. X2.LE.XMX .AND. Y2.GE.YMN .AND. Y2.LE.YMX ) THEN
	   CALL PLOT1Z_LINEV( X1, Y1, X2, Y2 )
	   RETURN
	ELSE IF( ( X1.LE.XMN .AND. X2.LE.XMN ) .OR.
	1        ( X1.GE.XMX .AND. X2.GE.XMX ) .OR.
	2        ( Y1.LE.YMN .AND. Y2.LE.YMN ) .OR. 	
	3        ( Y1.GE.YMN .AND. Y2.GE.YMN ) )THEN
	   RETURN
	ELSE
	   A1 = 0.0
	   IF( X2 - X1.NE.0.0 ) A1 = ( Y2 - Y1 ) / ( X2 - X1 )
	   B1 = Y1 - A1 * X1
	   A2 = 0.0
	   IF( Y2 - Y1.NE.0.0 ) A2 = ( X2 - X1 ) / ( Y2 - Y1 )
	   B2 = X1 - A2 * Y1
	   IF( X1.LT.XMN ) THEN
	      X1 = XMN
	      Y1 = B1 + A1 * X1
	   END IF
	   IF( X2.LT.XMN ) THEN
	      X2 = XMN
	      Y2 = B1 + A1 * X2
	   END IF
	   IF( X1.GT.XMX ) THEN
	      X1 = XMX
	      Y1 = B1 + A1 * X1
	   END IF
	   IF( X2.GT.XMX ) THEN
	      X2 = XMX
	      Y2 = B1 + A1 * X2
	   END IF
	   IF( Y1.LT.YMN ) THEN
	      Y1 = YMN
	      X1 = B2 + A2 * Y1
	   END IF
	   IF( Y2.LT.YMN ) THEN
	      Y2 = YMN
	      X2 = B2 + A2 * Y2
	   END IF
	   IF( Y1.GT.YMX ) THEN
	      Y1 = YMX
	      X1 = B2 + A2 * Y1
	   END IF
	   IF( Y2.GT.YMX ) THEN
	      Y2 = YMX
	      X2 = B2 + A2 * Y2
	   END IF
	   IF( X1.GE.XMN .AND. X1.LE.XMX .AND. Y1.GE.YMN .AND. Y1.LE.YMX .AND.
	1      X2.GE.XMN .AND. X2.LE.XMX .AND. Y2.GE.YMN .AND. Y2.LE.YMX ) THEN
	      CALL PLOT1Z_LINEV( X1, Y1, X2, Y2 )
	      RETURN
	   ELSE
	      RETURN
	   END IF
	END IF
	RETURN
	END

C
C============================================================================
	SUBROUTINE MODEL( STROX, STROY)
C
C	Subroutine to mark atom positions and draw bonds
C	This version has the option of generating closest contact "bonds"
C	excluding those involving 1st, 2nd and 3rd bonded neighbors.
C
	COMMON /IOSET/ IO2, IO3, IO4, IO5, IO6, IO7, IO8
C
	COMMON /MOD/ AM(3, 3, 36), AT(3, 36), XMIN3(3), XMAX3(3), 
	1 DZ3, DZ3N, DZ2, DZ2N, DMX(4), CHAR(4, 20), NSYM, PWIDE(4), 
	2 ICHR(4), HT(4), HT2(4), SCALE, RAD(4, 20), RDEL(4), 
	3 IMASK(36, 4), NSET(4, 20), ATMFIL(80, 4), VRAD(4, 20), 
	4 VTYP(4, 20), IFVDW(4), IFVDW0
C
	COMMON /DAT/ TI(10), XMIN(3), XMAX(3), DX(3), XCHAR(3), ACHAR(3), 
	1  IDIM(2), JDIM(2), NEXT, XSHX, XSHY, XSHZ, YSHY, YSHZ, ZSHZ, 
	2  XMPMIN, XMPMAX, YMPMIN, YMPMAX, ISCALE, NSDONE, NIP, ZMAP, A, 
	3  RMIN, RMAX, RINC, NSEC, ANG, RC
C
	COMMON /ATMLST/ XI(3, 4500), KI(4500), J
C
	REAL*4    XI, X(3), X1(3), X2(3), X3(3)
	LOGICAL*1 KI, ATMFIL
	INTEGER*4 PWIDE, CHAR, BNDPAR(50), BNDTYP(50), NBOND
C	INTEGER*4 BNDK1(50), BNDK2(50)
	REAL*4    DBOND(50), ADIS(4), FJUNK
C	REAL*4    DBND1(50), DBND2(50)
     	CHARACTER LABEL*4
     	EQUIVALENCE( LABEL, LABLL ), ( LABEL, FJUNK )
C
	IF( NSYM.EQ.0)RETURN
C
	DO 100 II = 1, 4
     	IF( IFCTRL_C() ) RETURN
C  Changed 3/4/87 from return to go to 100. MMT.
	IF( DMX(II).EQ.0.0 .AND. HT(II).EQ.0.0 ) GO TO 100
	DMX2  = DMX(II)**2
	IF( DMX(II).LT.0.0 )THEN
	   IBTYP = 1
	   DO K = 1, 20
	      IF( ( RDEL(II) + 2. * VRAD(II, K) )**2.GT.DMX2 ) 
	1         DMX2 = ( RDEL(II) + 2. * VRAD(II, K) )**2
	   END DO
	ELSE IF( DMX(II).GT.0.0 )THEN
	   IBTYP    = 0
	   DO K = 1, 20
	      IF( ( RDEL(II) + 2. * RAD(II, K) )**2.GT.DMX2 ) 
	1         DMX2 = ( RDEL(II) + 2. * RAD(II, K) )**2
	   END DO
	END IF
	XMIN3(3) = -SQRT( DMX2) - DZ2
	XMAX3(3) =  SQRT( DMX2) + DZ2
	ZMAPO    = ZMAP * ZSHZ
	J = 0
C
	CALL ATMLST( ZMAPO, II)
C
C	mark atom position.
C
	JN = J
	IF( JN.EQ.0)GO TO 100
C
	CALL IFSPEN(PWIDE(II)-1)
	XMPMIN1 = XMPMIN + 1.2 * HT2(II) / SCALE
	XMPMAX1 = XMPMAX - 1.2 * HT2(II) / SCALE
	YMPMIN1 = YMPMIN + HT2(II) / SCALE
	YMPMAX1 = YMPMAX - HT2(II) / SCALE
	DO J = 1, JN
     	IF( IFCTRL_C() ) RETURN
	   IF( XI(3, J).GE.DZ2N .AND. XI(3, J).LE.DZ2 .AND.
	1  XI(1, J).GE.XMPMIN1 .AND. XI(1, J).LE.XMPMAX1 .AND.
	2  XI(2, J).GE.YMPMIN1 .AND. XI(2, J).LE.YMPMAX1 .AND.
	3  ICHR(II).NE.0 ) THEN
	      IK=KI(J)
     	      IF( CHAR(II,IK).NE.0 ) THEN
	      XC=XI(1,J)+STROX*XI(3,J)
	      YC=XI(2,J)+STROY*XI(3,J)
     	      LABLL = CHAR(II,IK)
     	      CALL PLOT1Z_TEXT(XC,YC,LABEL,HT(II))
     	      ENDIF
	   END IF
	END DO
C
C	Find and draw bonds.
C	If PWIDE < 0, fill bonds.
C	If PWIDE = 0, single line.
C	If PWIDE > 0, open bonds.
C
	IF( DMX(II).EQ.0.0 )GO TO 100
	CALL IFSPEN(PWIDE(II))
	DWID  = 0.01
	NWID  = IFIX( 0.5 + 1.414**IABS( PWIDE(II)) )
	NWID1 = 1
	IF( PWIDE(II).EQ.0 ) NWID = 0
	IF( PWIDE(II).GT.0 ) THEN
	   DWID  = 0.01 * NWID
	   NWID  = 1
	   NWID1 = 2
	END IF
	NBOND = 0
	DO J = 1, JN - 1
     	IF( IFCTRL_C() ) RETURN
	   IF( NBOND.GT.0 )THEN
	      DO J2 = 1, NBOND
	         BNDPAR(J2) = 0
	         BNDTYP(J2) = 0
	      END DO
	      NBOND = 0
	   END IF
	   IF( IBTYP.EQ.0 ) J1 = J + 1
	   IF( IBTYP.EQ.1 ) J1 = 1
	   DO J2 = J1, JN
     	IF( IFCTRL_C() ) RETURN
	      D  = 0.0
	      DR = 0.0
	      DO K = 1, 3
	         D  = D  + ( XI(K, J) - XI(K, J2) )**2
	         IF( D.GT.DMX2 )GO TO 160
	      END DO
	      IF( D.LT.0.001 )GO TO 160
	      D   = SQRT( D)
	      KJ1 = KI( J)
	      KJ2 = KI(J2)
	      IF( ( IBTYP.EQ.0 )                             .AND.
	1         ( D.LE.RDEL(II) + RAD(II, KJ1) + RAD(II, KJ2) )  .AND.
	1         ( XI(3, J).GT.DZ3N .OR. XI(3, J2).GT.DZ3N )  .AND.
	1         ( XI(3, J).LT.DZ3  .OR. XI(3, J2).LT.DZ3  )  .AND.
	1         ( NBOND.LT.50 )    .AND. ( J.NE.J2 )    )THEN
	         NBOND = NBOND + 1
	         BNDPAR(NBOND) = J2
	         BNDTYP(NBOND) = 1    ! covalent bond
	         DBOND(NBOND)  = D
c
	      ELSE IF( ( IBTYP.EQ.1 )  .AND.
	1           ( NBOND.LT.50 )    .AND. ( J.NE.J2 )    )THEN
	         NBOND = NBOND + 1
	         IF( D.LE.(RDEL(II) + RAD(II, KJ1) + RAD(II, KJ2)))THEN
	            BNDPAR(NBOND) = -J2
	            BNDTYP(NBOND) =   1	    ! covalent bond
	            DBOND(NBOND)  =   D
	         ELSE IF( 
	1           ( D.LE.(RDEL(II) + VRAD(II, KJ1) + VRAD(II, KJ2)) )
	1              	  .AND. 
	1     	    ( XI(3, J).GT.DZ3N .OR. XI(3, J2).GT.DZ3N )  .AND.
	1           ( XI(3, J).LT.DZ3  .OR. XI(3, J2).LT.DZ3  )  )THEN
	            BNDPAR(NBOND) =  J2
	            BNDTYP(NBOND) =   4	    ! good
		    DBOND(NBOND)  =   D
	         ELSE IF( D.LE.
	1	       	  ( RDEL(II) + VRAD(II, KJ1) + VRAD(II, KJ2) ) )THEN
	            BNDPAR(NBOND) = -J2
	            BNDTYP(NBOND) =   4	    ! good but out of section.
	            DBOND(NBOND)  =   D
	         ELSE
		    BNDPAR(NBOND) = -J2
		    BNDTYP(NBOND) =   5     ! to long but could be 1-3 in a 1-4
		    DBOND(NBOND)  =   D
	         END IF
	         IF( J2.LT.J )THEN
	            BNDPAR(NBOND) = -J2	    ! will be picked up in other direction.
	         END IF
	      END IF
 160	      CONTINUE
	   END DO
c
c	   if(ibtyp.eq.1)then
c	      write( 27, 201) J, JJ2, NBOND, 
c	1       ( J2, BNDPAR(J2), BNDTYP(J2), BNDK1(J2), BNDK2(J2), 
c	1     	       DBND1(J2), DBND2(J2), DBOND(J2), J2=1,NBOND )
c	   end if
c	
	   IF( IBTYP.EQ.1 .AND. NBOND.GT.1 )THEN
	      DO J2 = 1, NBOND
     	IF( IFCTRL_C() ) RETURN
	         J22 = IABS( BNDPAR(J2))
		 IF( BNDTYP(J2).GT.1 )THEN
	            DO J21 = 1, NBOND
     	IF( IFCTRL_C() ) RETURN
C
C	      	       EXCLUDE 1-3
C
	               IF( BNDTYP(J21).EQ.1 )THEN
	                  J23 	 =  IABS( BNDPAR(J21))
	                  D  	 =  0.0
	                  DR 	 =  0.0
	                  DO K = 1, 3
	                     D  = D + ( XI(K, J22) - XI(K, J23) )**2
	                  END DO
	                  D  	 =  SQRT(D)
	                  KJ1 	 =  KI(J22)
	                  KJ2 	 =  KI(J23)
	                  IF( D.LE.( RDEL(II) + 
	1		     	     RAD(II, KJ1) + RAD(II, KJ2) ) )THEN
	                     IF( BNDTYP(J2).GT.2 )THEN
	        	        BNDPAR( J2 ) = -IABS( BNDPAR( J2 ) )
	        	        BNDTYP( J2 ) =  2
	         	     ELSE
	      	       	     	BNDPAR( J2 ) = -IABS( BNDPAR( J2 ) )
	         	     END IF
	         	  END IF
	               END IF
	            END DO
		 END IF
	      END DO
C
C	      	       EXCLUDE 1-4
C
	      DO J2 = 1, NBOND
     	IF( IFCTRL_C() ) RETURN
	         J22 = BNDPAR(J2)
	         IF( J22.GT.0 .AND. BNDTYP(J2).EQ.4 )THEN
	            DO J21 = 1, NBOND
     	IF( IFCTRL_C() ) RETURN
	               IF( BNDTYP(J21).EQ.2 )THEN
	                  J23 	 =  IABS( BNDPAR(J21))
	                  D  	 =  0.0
	                  DR 	 =  0.0
	                  DO K = 1, 3
	                     D  = D + ( XI(K, J22) - XI(K, J23) )**2
	                  END DO
	                  D  	 =  SQRT(D)
	                  KJ1 	 =  KI(J22)
	                  KJ2 	 =  KI(J23)
	                  IF( D.LE.( RDEL(II) + 
	1		     	     RAD(II, KJ1) + RAD(II, KJ2) ) )THEN
	                     IF( BNDTYP(J2).GT.3 )THEN
			     	BNDPAR( J2 ) = -IABS( BNDPAR( J2 ) )
			     	BNDTYP( J2 ) =  3
		       	     ELSE
			        BNDPAR( J2 ) = -IABS( BNDPAR( J2 ) )
		       	     END IF
			  END IF
	               END IF
	            END DO
	         END IF
	      END DO
C
	   END IF
c
c	   if( ibtyp.eq.1 )then
c	      write( 27, 201) J, JJ2, NBOND, 
c	1       ( J2, BNDPAR(J2), BNDTYP(J2), BNDK1(J2), BNDK2(J2), 
c	1     	       DBND1(J2), DBND2(J2), DBOND(J2), J2=1,NBOND )
c201	      format( /, ' ', i5, 2i5, i3, i5, i2, 2i3, 3f6.2, 
c	1     	      /, ( ' ',15x, i3, i5, i2, 2i3, 3f6.2) )
c	   end if
c	
	   DO JJ2 = 1, NBOND
     	IF( IFCTRL_C() ) RETURN
	      IDIS = 0
	      J2 = BNDPAR(JJ2)
	      IF( J2.LE.0 )GO TO 60
	      XIN = XI(1, J)  + XI(3, J)  * STROX
	      YIN = XI(2, J)  + XI(3, J)  * STROY
	      XTR = XI(1, J2) + XI(3, J2) * STROX
	      YTR = XI(2, J2) + XI(3, J2) * STROY
	      IF( (XIN.LT.XMPMIN).AND.(XTR.LT.XMPMIN))GO TO 60
	      IF( (XIN.GT.XMPMAX).AND.(XTR.GT.XMPMAX))GO TO 60
	      IF( (YIN.LT.YMPMIN).AND.(YTR.LT.YMPMIN))GO TO 60
	      IF( (YIN.GT.YMPMAX).AND.(YTR.GT.YMPMAX))GO TO 60
	      DR = DBOND(JJ2)
	      IF( XI(3, J) .LT.DZ3N .OR. XI(3, J) .GT.DZ3 .OR.
	1         XI(3, J2).LT.DZ3N .OR. XI(3, J2).GT.DZ3 ) THEN
	         A1 = (XTR - XIN) / (XI(3, J2) - XI(3, J))
	         A2 = (YTR - YIN) / (XI(3, J2) - XI(3, J))
	         B1 = XIN - A1 * XI(3, J)
	         B2 = YIN - A2 * XI(3, J)
	         IF( XI(3, J).LT.DZ3N ) THEN
	            XIN = A1 * DZ3N + B1
	            YIN = A2 * DZ3N + B2
	         ELSE IF( XI(3, J).GT.DZ3 ) THEN
	            XIN = A1 * DZ3 + B1
	            YIN = A2 * DZ3 + B2
	         END IF
	         IF( XI(3, J2).LT.DZ3N ) THEN
	            XTR = A1 * DZ3N + B1
	            YTR = A2 * DZ3N + B2
	         ELSE IF( XI(3, J2).GT.DZ3 ) THEN
	            XTR = A1 * DZ3 + B1
	            YTR = A2 * DZ3 + B2
	         END IF
	      ELSE
	      	 IDIS = 1
	      END IF
	      DDX = XTR - XIN 
	      DDY = YTR - YIN 
	      DDR = SQRT( DDX**2 + DDY**2 )
	      IF( DDR.EQ.0.0 )THEN
	         DO ID = -NWID, NWID, NWID1
	            DDX  = DWID * ID
	            XIN1 = XIN - DDY
	            XTR1 = XTR - DDY
	            YIN1 = YIN - DDX
	            YTR1 = YTR - DDX
	            XIN1 = XIN - DDX
	            XTR1 = XTR - DDX
	            XIN2 = XIN + DDX
	            XTR2 = XTR + DDX
	            YIN1 = YIN - DDX
	            YTR1 = YTR + DDX
	            YIN2 = YIN - DDX
	            YTR2 = YTR + DDX
     	CALL LINEV1( XIN1, YIN1, XTR1, YTR1, XMPMIN, XMPMAX, 
     +				YMPMIN, YMPMAX)
     	CALL LINEV1( XIN2, YIN2, XTR2, YTR2, XMPMIN, XMPMAX, 
     +				YMPMIN, YMPMAX)
     	CALL LINEV1( XIN1, YIN1, XIN2, YIN2, XMPMIN, XMPMAX, 
     +				YMPMIN, YMPMAX)
     	CALL LINEV1( XTR1, YTR1, XTR2, YTR2, XMPMIN, XMPMAX, 
     +				YMPMIN, YMPMAX)
	         END DO
c
c	      	 if( ibtyp.eq.1 )then
c		    write( 27, 202) j, j2, xin, yin, xtr, ytr
c		 end if
c		 
	      ELSE
	         DDX1  =  DWID * DDX / DDR
	         DDY1  =  DWID * DDY / DDR
	         DO ID = -NWID, NWID, NWID1
	            DDY  = DDY1 * ID
	            DDX  = DDX1 * ID
	            XIN1 = XIN - DDY
	            XTR1 = XTR - DDY
	            YIN1 = YIN + DDX
	            YTR1 = YTR + DDX
     	CALL LINEV1( XIN1, YIN1, XTR1, YTR1, XMPMIN, XMPMAX,  
     +				YMPMIN, YMPMAX)
	         END DO
c
		 IF( IFVDW(II).LT.0 .AND. IDIS.NE.0 )THEN
C
C	write in bond length?
C
	      	    XC = ( XIN + XTR ) * 0.5
	      	    YC = ( YIN + YTR ) * 0.5
		    IF( ( ( XC -      HT2(II) / SCALE ).GT.XMPMIN ) .AND.
	1	        ( ( XC +      HT2(II) / SCALE ).LT.XMPMAX ) .AND.
	2	        ( ( YC - 4. * HT2(II) / SCALE ).GT.YMPMIN ) .AND.
	3	        ( ( YC + 4. * HT2(II) / SCALE ).LT.YMPMAX )	)THEN
		       ENCODE( 4, 1, ADIS) DR
 1	      	       FORMAT( F4.2)		    
     		       FJUNK = ADIS(1)
     	      	       CALL PLOT1Z_TEXT(XC,YC,LABEL,HT(II))
	      	    END IF
		 END IF
c
	      END IF
 60	      CONTINUE
	   END DO
	END DO
C
 100	CONTINUE
	RETURN
	END

	SUBROUTINE ATMLST( ZMAPO, II)
C
C	Subroutine to fill atom list.  7/18/85 G.J.Q.
C
	COMMON /IOSET/ IO2, IO3, IO4, IO5, IO6, IO7, IO8, MAXTYP
C
	COMMON /MOD/ AM(3, 3, 36), AT(3, 36), XMIN3(3), XMAX3(3), 
     +  DZ3, DZ3N, DZ2, DZ2N, DMX(4), CHAR(4, 20), NSYM, PWIDE(4), 
     +  ICHR(4), HT(4), HT2(4), SCALE, RAD(4, 20), RDEL(4), 
     +  IMASK(36, 4), NSET(4, 20), ATMFIL(80, 4), VRAD(4, 20), 
     +  VTYP(4, 20), IFVDW(4), IFVDW0
C
	COMMON /DAT/ TI(10), XMIN(3), XMAX(3), DX(3), XCHAR(3), 
     +  ACHAR(3), IDIM(2), JDIM(2), NEXT, XSHX, XSHY, XSHZ, YSHY,
     +  YSHZ, ZSHZ, XMPMIN, XMPMAX, YMPMIN, YMPMAX, ISCALE, NSDONE, 
     +  NIP, ZMAP, A, RMIN, RMAX, RINC, NSEC, ANG, RC
C
	COMMON /ATMLST/ XI(3, 4500), KI(4500), J
C
	REAL*4    XI, X(3), X1(3), X2(3)
	LOGICAL*1 KI, ATMFIL
C
C	master loop over seperate atom sets.
C
	J = 0
	OPEN( UNIT = IO4, NAME = ATMFIL(1, II), TYPE = 'OLD', 
	1 FORM = 'FORMATTED', READONLY)
C
	IEND = 0
     	MAXTYP = 0
 10	CONTINUE
C	DO WHILE( IEND.EQ.0 .AND. J.LE.4500 )
	   READ(IO4, 11, END = 99)N, IK, X(1), X(2), X(3)
 11	   FORMAT(I5, 10X, I3, 3F10.4)
	   DO K = 1, 19, 2
	      IF( (N.GE.NSET(II, K)).AND.(N.LE.NSET(II, K + 1)))GO TO 120
	   END DO
	   GO TO 10
 120	   CONTINUE
	   IF( IK.LE.0)GO TO 10
	   DO 20 I = 1, NSYM
	      IF( IMASK(I, II).EQ.0)GO TO 20
	      DO I1 = 1, 3
	         X1(I1) = AT(I1, I)
	         DO I2 = 1, 3
	            X1(I1) = X1(I1) + AM(I1, I2, I) * X(I2)
	         END DO
	      END DO
c	      CALL XYZ2ZPR( X1, X2, FIBREP, RSPRTR)
c	      X2(3) = X2(3) - ZMAPO
	      IOK   = 1
	      DO I1 = 1, 3
	         IF( IOK.EQ.1 )THEN
c  Changed 3/4/87. x2->x1. MMT
	            IF( X1(I1).LT.XMIN3(I1) .OR.
	1               X1(I1).GT.XMAX3(I1) ) IOK = 0
	         END IF
	      END DO
	      IF( IOK.EQ.1 )THEN
	         J = J + 1
	         IF( J.GT.4500 )GO TO 38
	         DO I1 = 1, 3
	            XI(I1, J) = X1(I1)
c	            XR(I1, J) = X2(I1)
	         END DO
	         IF( IK.GT.20 )IK = 20
	         KI(J) = IK
     		 IF( IK.GT.MAXTYP ) MAXTYP = IK
	      END IF
 20	   CONTINUE
	GO TO 10
C
 38	WRITE(IO6, 39)
 39	FORMAT('0*** TOO MANY ATOMS ***')
	CALL QUIT(0)
 99	CONTINUE
	CLOSE(UNIT = IO4)
C
	RETURN
	END

     	SUBROUTINE PLOT1Z_QPLOT(X,Y,IC)
C
C  PLOT A POINT OR A LINE OR MOVE PEN
C
C	X,Y=  COORDINATES, IN INCHES
C	IB= PEN CONTROL. IC=2, PEN IS DOWN DURING MOVE,
C			 IC=3, PEN IS UP DURING MOVE
C			 IC=-2 OR -3, A NEW ORIGIN IS DEFINED
C
     	PARAMETER CIN= 25.4
	INTEGER   IC, IB, IP
	COMMON/OFFSET/XOFF,XFAC,YOFF,YFAC,OFFX,OFFY

	XN = X
	YN = Y
	IB = IC
     	IP = IABS(IB)
     	IF( IP.GT.10 ) THEN
     	    IB = ISIGN(IP-10,IB)
     	    XN = (XN-XOFF) / XFAC / CIN
     	    YN = (YN-YOFF) / YFAC / CIN
     	ENDIF
	CALL IFPLOT_VECT( XN+OFFX, YN+OFFY, IB )
	RETURN
	END

	SUBROUTINE PLOT1Z_OFFSET(XOFFX,XFACX,YOFFX,YFACX)
C
	COMMON/OFFSET/XOFF,XFAC,YOFF,YFAC,OFFX,OFFY
	DATA XOFF,XFAC,YOFF,YFAC /0.0,1.0,0.0,1.0/
C
	XOFF = XOFFX
	XFAC = XFACX
	YOFF = YOFFX
	YFAC = YFACX
	RETURN
	END

     	SUBROUTINE PLOT1Z_MARK(X,Y,SIZE,ISYM)
C
C	THIS ROUTINE PUTS ON GRAPHICS MARKERS.
C
     	PARAMETER CIN= 25.4
	COMMON/OFFSET/XOFF,XFAC,YOFF,YFAC,OFFX,OFFY

     	XN = (X-XOFF) / XFAC / CIN + OFFX
     	YN = (Y-YOFF) / YFAC / CIN + OFFY
     	CALL IFTEXT_MARK( XN, YN, SIZE, 0.0, ISYM )
     	RETURN
     	END

     	SUBROUTINE PLOT1Z_TEXT(X,Y,TEXT,SIZE)
C
C	THIS ROUTINE HANDLES GRAPHTEXT GENERATION.
C
     	PARAMETER CIN= 25.4
     	CHARACTER TEXT*(*)
	COMMON/OFFSET/XOFF,XFAC,YOFF,YFAC,OFFX,OFFY

     	SZ = SIZE/CIN
     	OF = 0.5*SZ
     	XN = (X-XOFF) / XFAC / CIN - OF + OFFX
     	YN = (Y-YOFF) / YFAC / CIN - OF + OFFY
     	CALL IFTEXT( XN, YN, TEXT, SZ, SZ, 0.0 )
     	RETURN
     	END

