1
2*	PROGRAM RIBBON
3*
4*       Program to set up input for RENDER (RASTER3D package)
5*	to draw ribbon diagram.  The RIBBON routine itself is simply
6*	extracted from CCP FRODO.  The original invoked a bspline feature
7*	of the ps300;  I have replaced this with a spline equation gotten
8*	from Larry Andrews.  Conversion from ribbon edges to solid rendering
9*	is my own hacking.
10*				Ethan Merritt	-  8-Nov-1988
11*	Slightly modified code to guarantee output of triangles with
12*	vertices in correct order for triangular mesh algorithms EAM Sep 90
13*
14*	Usage:	ribbon [-h] [-dn] pdbfile > setup.r3d
15*		ribbon [-h] -dn -         to take pdb records from stdin
16*
17*	Input:	pdbfile
18*     			Brookhaven PDB-format file of atomic co-ordinates
19*			only C-alpha and O atoms are needed
20*		setup.matrix or setup.angles
21*			rotation matrix or angles applied to PDB coords
22*			(see writeup for SETUP/RENDER).
23*	Output:	stdout (new for DS5000 version)
24*			file suitable for input to RENDER
25*
26*	Interactive parameters:
27*		WIDTH of ribbon in Angstroms
28*		NUMBER of interpolated coordinates between successive C-alphas.
29*		COLOR scheme for ribbon
30*		    0 or 1: solid color (RGB values from color1 below)
31*		    2:	    shade from color1 at 1st res to color2 at last res
32*		    3:      front of ribbon is color1, back of ribbon is color2
33*		    4:      shade front as in scheme 2, back is color 3
34*                   5:      each chain is new color (from successive input
35*                               (COLOR cards at start of input file)
36*		    6:      use prefixed COLOR cards (as in SETUP/RENDER)
37*				(implemented 4-Aug-1997 EAM)
38*		COLOR1,COLOR2,COLOR3	RGB components (9f8.0)
39*
40      INCLUDE 'VERSION.incl'
41c
42      INTEGER INPUT, OUTPUT, NOISE
43      PARAMETER (OUTPUT=6, NOISE=0)
44      PARAMETER (MAXCOL=5000, MAXATM=10000)
45C     REAL RGB(3,MAXCOL)
46      REAL RADIUS(MAXCOL)
47      REAL RAD
48      CHARACTER*24 MASK(MAXCOL),TEST
49      CHARACTER*80 ATOM(MAXATM),CARD
50      LOGICAL SMATCH
51c
52C	Ethan Merritt	Oct 1988
53C	Modified to read in 3x3 view matrix (e.g. from CCP FRODO view command)
54C	from file.  Matrix is applied before
55C	finding translation, center, and scale.  Afterwards the input matrix
56C	to RENDER is therefore the identity matrix.
57C EAM Aug 1997 - Honor COLOUR requests
58C EAM Nov 1999 - remove all (q) formats
59C EAM Jan 2010 - declare and initialize RAD
60C
61c
62	common /COLORS/ ischeme, cindex, COLOR1(3), COLOR2(3), COLOR3(3)
63     &  		,RGB(3,MAXCOL)
64 	integer		cindex
65	common /SPAM/   natm, SPAM(4,MAXATM), SCAM(MAXATM)
66	integer		SCAM
67	common /FLAGS/  mflag, hflag, dflag
68	logical		mflag, hflag, dflag
69c
70	character*64	in_file
71	character*32	flags
72	character*80	line
73	common /matrix/ matrix, coords
74	real		matrix(3,3), coords(3)
75	data		matrix / 1.,0.,0.,0.,1.,0.,0.,0.,1. /
76c
77c	-h causes the header records not to be printed
78c	-m [now obsolete because always in force] uses format
79c	   mixed object types in output file
80c	-d suppresses interactive input
81c
82	hflag = .FALSE.
83	dflag = .FALSE.
84	mflag = .TRUE.
85c
86	narg  = iargc()
87	do i = 1,narg
88	    call getarg( i, flags )
89	    if (flags(1:2) .eq. '-h') then
90	    	hflag = .TRUE.
91	    else if (flags(1:2) .eq. '-d') then
92	    	dflag = .TRUE.
93		read (flags(3:4),'(I1)') ischeme
94	    end if
95	end do
96c
97	call getarg( narg, in_file )
98	if (in_file(1:1) .eq. '-') then
99	    INPUT = 5
100	else
101	    INPUT = 1
102	    open( unit=INPUT, file=in_file, status='OLD' )
103	end if
104c
105    3	format(a,a)
106c
107	call view_matrix
108c
109	NCOL = 0
110	NATM = 0
111	ASPECT = 1280./1024.
112c
113	if (hflag) goto 10
114c
115      WRITE(OUTPUT,'(A,A)') 'C-alpha ribbon - Raster3D ',VERSION
116      WRITE(OUTPUT,'(A)') '80 64     tiles in x,y'
117      WRITE(OUTPUT,'(A)') ' 8  8     pixels (x,y) per tile'
118      WRITE(OUTPUT,'(A)') '4         anti-aliasing 3x3 into 2x2 pixels'
119      WRITE(OUTPUT,'(A)') '0 0 0     black background'
120      WRITE(OUTPUT,'(A)') 'F         no, ribbons cast funny shadows'
121      WRITE(OUTPUT,'(A)') '25        Phong power'
122      WRITE(OUTPUT,'(A)') '0.15      secondary light contribution'
123      WRITE(OUTPUT,'(A)') '0.05      ambient light contribution'
124      WRITE(OUTPUT,'(A)') '0.25      specular reflection component'
125      WRITE(OUTPUT,'(A)') '4.0       eye position'
126      WRITE(OUTPUT,'(A)') '1 1 1     main light source position'
127c
12810    CONTINUE
129        READ(INPUT,'(A80)',END=50) CARD
130        IF (CARD(1:4).EQ.'COLO') THEN
131          NCOL = NCOL + 1
132          IF (NCOL.GT.MAXCOL) THEN
133            WRITE(NOISE,*) 'Colour table overflow.  Increase ',
134     &                     'MAXCOL and recompile.'
135            STOP 10
136          ENDIF
137          READ(CARD,'(6X,A24,3F8.3,F6.2)') MASK(NCOL),
138     &          (RGB(I,NCOL),I=1,3),RADIUS(NCOL)
139        ELSEIF ((CARD(1:4).EQ.'ATOM') .AND.
140     &         ( CARD(14:16).EQ.'CA ' .OR. CARD(14:16).EQ.'O  ')) THEN
141          NATM = NATM + 1
142          IF (NATM.GT.MAXATM) THEN
143            WRITE(NOISE,*) 'Atom array overflow.  Increase ',
144     &                     'MAXATM and recompile.'
145            STOP 20
146          ENDIF
147          ATOM(NATM) = CARD
148        ELSEIF (CARD(1:3).EQ.'END') THEN
149          GO TO 50
150        ENDIF
151        GO TO 10
152*     Come here when EOF or 'END' record is reached
15350    CONTINUE
154      IF (NATM.EQ.0) THEN
155        WRITE(NOISE,*) 'No atoms in input.'
156        STOP 30
157      ELSE
158	WRITE(NOISE,*) NATM,' atoms accepted from input.'
159      ENDIF
160      IF (NCOL.EQ.0) THEN
161        WRITE(NOISE,*) 'No colours in input.'
162c       STOP 40
163      ENDIF
164C
165      XMAX = -1E20
166      XMIN =  1E20
167      YMAX = -1E20
168      YMIN =  1E20
169      ZMAX = -1E20
170      ZMIN =  1E20
171      RAD  =  1.7
172
173      DO 100 IATM=1,NATM
174        CARD = ATOM(IATM)
175        TEST = CARD(7:30)
176            READ(CARD,82) coords
177   82	    format(30x,3f8.3)
178		x = coords(1)*matrix(1,1) + coords(2)*matrix(2,1)
179     1  	  + coords(3)*matrix(3,1)
180		y = coords(1)*matrix(1,2) + coords(2)*matrix(2,2)
181     1  	  + coords(3)*matrix(3,2)
182		z = coords(1)*matrix(1,3) + coords(2)*matrix(2,3)
183     1  	  + coords(3)*matrix(3,3)
184            SPAM(1,IATM) = X
185            SPAM(2,IATM) = Y
186            SPAM(3,IATM) = Z
187c           SPAM(4,IATM) = RAD
188C
189C	    EAM Aug 1997 - finally get around to honoring atom colors
190	    DO 84 ICOL = 1, NCOL
191		IF (SMATCH(TEST,MASK(ICOL))) THEN
192		    SCAM(IATM) = ICOL
193		    RAD = RADIUS(ICOL)
194		    SPAM(4,IATM) = RAD
195		    GOTO 86
196		ENDIF
197   84	    CONTINUE
198   86	    CONTINUE
199C
200            XMAX = MAX(XMAX,X+RAD)
201            XMIN = MIN(XMIN,X-RAD)
202            YMAX = MAX(YMAX,Y+RAD)
203            YMIN = MIN(YMIN,Y-RAD)
204            ZMAX = MAX(ZMAX,Z+RAD)
205            ZMIN = MIN(ZMIN,Z-RAD)
206100   CONTINUE
207      XMID = (XMAX+XMIN)/2.
208      YMID = (YMAX+YMIN)/2.
209      ZMID = (ZMAX+ZMIN)/2.
210      TX = -XMID
211      TY = -YMID
212      TZ = -ZMID
213      IF (ASPECT.GE.1.) THEN
214*       The X direction is wider than the Y
215        XROOM = ASPECT
216        YROOM = 1.
217        ZROOM = 2.
218      ELSE
219        XROOM = 1.
220        YROOM = ASPECT
221        ZROOM = 2.
222      ENDIF
223      XSPAN = XMAX-XMIN
224      YSPAN = YMAX-YMIN
225      ZSPAN = ZMAX-ZMIN
226      SCALE = MAX(XSPAN/XROOM,YSPAN/YROOM,ZSPAN/ZROOM)
227*     Leave a little extra room as a border:
228      SCALE = SCALE / 0.90
229      if (hflag) goto 129
230      WRITE(OUTPUT,120) TX,TY,TZ,SCALE
231120   FORMAT('1 0 0 0   input co-ordinate + radius transformation'/
232     &       '0 1 0 0'/
233     &       '0 0 1 0'/
234     &       4F10.3)
235      if (mflag) then
236	WRITE (OUTPUT,'(A)') '3         mixed object types'
237	WRITE (OUTPUT,'(A)') '(9F8.3,2x,3f5.2)'
238	WRITE (OUTPUT,'(A)') '(11F8.3)'
239	WRITE (OUTPUT,'(A)') '(11F8.3)'
240      else
241	WRITE (OUTPUT,'(A)') '1         all objects are triangles'
242	WRITE (OUTPUT,'(A)') '(9F8.3,2x,3f5.2)'
243      end if
244  129	continue
245	write (noise,'(/)')
246	write (noise,153) 'X  min max:', XMIN, XMAX
247	write (noise,153) 'Y  min max:', YMIN, YMAX
248	write (noise,153) 'Z  min max:', ZMIN, ZMAX
249	write (noise,153) '     scale:', SCALE
250  153	format(1x,a,3f8.2)
251c
252c
253	if (dflag) then
254	    width = 1.5
255	    offset = 1.2
256	    nchord = 5
257	    if (ischeme .le. 0 .or. ischeme .gt. 6) ischeme = 2
258	     call vload( color1, 0.0, 0.0, 0.4 )
259	     call vload( color2, 0.5, 0.0, 0.0 )
260	     call vload( color3, 0.6, 0.6, 0.6 )
261	else
262	width = 0
263	write (noise,3) 'Width of ribbon (default 1.5A): '
264	read  (5,'(A80)') line
265	read  (line,*,end=154,err=154) width
266  154	continue
267	if (width.le.0) width = 1.5
268c	Original RIBBON used bspline smoothing, which requires "offset"
269c	because smoothed curve doesn't go through guide points.
270	write (noise,3) 'Offset from CA position (default 1.2A): '
271	read  (5,'(A80)') line
272	read  (line,*,end=156,err=156) offset
273  156	continue
274	if (offset.le.0) offset = 1.2
275	write (noise,3) 'Chords per residue (default = 10): '
276	read  (5,'(A80)') line
277	read  (line,*,end=158,err=158) nchord
278  158	continue
279	if (nchord.le.1) nchord = 10
280
281	write (noise,160)
282  160	format(' Coloring schemes available:',
283     1	/,' 0 or 1: solid color (RGB values from color1 below)',
284     2	/,'      2: shade from color1 at 1st res to color2 at last res',
285     3	/,'      3: front of ribbon is color1, back of ribbon is color2',
286     4	/,'      4: shade front as in scheme 2, back is color 3',
287     5	/,'      5: new color for each chain (requires COLOUR cards)')
288	write (noise,3) 'Coloring scheme: '
289	read  (5,'(A80)') line
290	read  (line,*,end=161,err=161) ischeme
291  161	continue
292	if (ischeme.le.0 .or. ischeme.gt.6) ischeme = 1
293	if (ischeme .eq. 1) write (noise,3)
294     1      'COLOR1 (RGB values, 3f8.0): '
295	if (ischeme .eq. 2) write (noise,3)
296     1      'COLOR1, COLOR2 (RGB values, 6f8.0): '
297	if (ischeme .eq. 3) write (noise,3)
298     1      'COLOR1, COLOR2 (RGB values, 6f8.0): '
299	if (ischeme .eq. 4) write (noise,3)
300     1      'COLOR1, COLOR2, COLOR3 (RGB values, 9f8.0): '
301	if (ischeme .lt. 5) then
302	    read  (5,'(A80)') line
303	    if (line.eq.' ') goto 163
304	    read  (line,*,end=163,err=163) color1,color2,color3
305        endif
306	goto 164
307  163	continue
308	     call vload( color1, 0.0, 0.0, 0.4 )
309	     call vload( color2, 0.5, 0.0, 0.0 )
310	     call vload( color3, 0.6, 0.6, 0.6 )
311  164	continue
312	if (ischeme .eq. 3) then
313		color3(1) = color2(1)
314		color3(2) = color2(2)
315		color3(3) = color2(3)
316	end if
317c	end of -d suppression
318	end if
319	write (noise,169) ischeme,color1,color2,color3
320  169	format(' color scheme',i3,/,3(3x,3f6.3))
321	cindex = 1
322c
323	call ribbon( 2, width, nchord, offset, natm )
324c
325      END
326      LOGICAL FUNCTION SMATCH (SUBJ, MASK)
327      CHARACTER*24 SUBJ,MASK
328      SMATCH = .FALSE.
329      DO 10 I = 1, 24
330        IF (SUBJ(I:I).NE.MASK(I:I) .AND. MASK(I:I).NE.'#') RETURN
33110    CONTINUE
332      SMATCH = .TRUE.
333      RETURN
334      END
335
336	subroutine view_matrix
337c
338	common /matrix/ matrix, coords
339	real		matrix(3,3), coords(3)
340c
341	real		phiX, phiY, phiZ
342	parameter (noise = 0)
343	parameter (R2D = 180./3.1415927)
344
345	open (unit=3, file='setup.matrix', status='OLD', err=100)
346		write (noise,3) ' View Matrix from file '
347		read (3,*) ((matrix(i,j),i=1,3),j=1,3)
348		write (noise,'(1x,3f9.5)') ((matrix(i,j),i=1,3),j=1,3)
349		close (3)
350
351		det = matrix(1,1) * matrix(2,2) * matrix(3,3)
352     1  	    + matrix(1,2) * matrix(2,3) * matrix(3,1)
353     2  	    + matrix(2,1) * matrix(3,2) * matrix(1,3)
354     3  	    - matrix(1,3) * matrix(2,2) * matrix(3,1)
355     4  	    - matrix(1,2) * matrix(2,1) * matrix(3,3)
356     5  	    - matrix(1,1) * matrix(2,3) * matrix(3,2)
357		write (noise,'(''       determinant ='',f8.3)') det
358
359		phiX = atan2( -matrix(3,2), matrix(3,3) )
360		phiY = atan2(  matrix(3,1), matrix(3,3) / cos(phiX) )
361		phiZ = atan2( -matrix(2,1), matrix(1,1) )
362		write (noise,3) ' View Angles from matrix',' '
363		write (noise,2) phiZ*R2D, phiY*R2D, phiX*R2D
364		return
365  100	continue
366
367	open (unit=3, file='setup.angles', status='OLD', err=200)
368		write (noise,3) ' View Angles from file '
369		read (3,*) phiZ, phiY, phiX
370		close (3)
371		write (noise,2) phiZ, phiY, phiX
372		cx = cos(phiX/R2D)
373		sx = sin(phiX/R2D)
374		cy = cos(phiY/R2D)
375		sy = sin(phiY/R2D)
376		cz = cos(phiZ/R2D)
377		sz = sin(phiZ/R2D)
378		matrix(1,1) = cz*cy
379		matrix(1,2) = sz*cx + cz*sy*sx
380		matrix(1,3) = sz*sx - cz*sy*cx
381		matrix(2,1) = -sz*cy
382		matrix(2,2) = cz*cx - sx*sy*sz
383		matrix(2,3) = cz*sx + sz*sy*cx
384		matrix(3,1) = sy
385		matrix(3,2) = -sx*cy
386		matrix(3,3) = cx*cy
387		write (noise,3) ' View Matrix from angles',' '
388		write (noise,'(1x,3f9.5)') ((matrix(i,j),i=1,3),j=1,3)
389		return
390  200	continue
391
392    2 	format(1x,'   phiZ =',f8.2,'   phiY =',f8.2,'   phiX =',f8.2)
393    3	format(/a,a)
394
395    	write (noise,*) ' No view matrix or angles provided'
396	return
397	end
398
399
400C
401	SUBROUTINE RIBDRW(GUIDE,NRIB,MAXRES,NPT,NCHORD)
402	integer	npt			! number of guide points
403	real	guide(4,MAXRES,NRIB)    ! 4 dim because E&S wanted it that way
404	integer	nchord			! how many interpolations per guide pt
405	parameter (MAXCOL = 5000)
406	integer	OUTPUT
407	parameter (OUTPUT = 6)
408C
409C	splining from Larry Andrews 7-Nov-1988
410C
411	parameter (nspln = 5000)	! maximum of (npt*nchord)
412	parameter (ndata =  500)	! maximum # guidepoints
413	parameter (ndata1 = 501)
414c
415	common /COLORS/ ischeme, cindex, COLOR1(3), COLOR2(3), COLOR3(3)
416     &			,RGB(3,MAXCOL)
417 	integer		cindex
418	common /FLAGS/  mflag, hflag, dflag
419	logical		mflag, hflag, dflag
420c
421c	real	s(ndata1)
422c	REAL	XP(4,NDATA)
423	real	smooth(4,nspln,2)	! npt*nchord points on splined curve
424	real	color(3)
425c
426	if (npt .gt. ndata) stop 'spline - TOO MANY GUIDE POINTS'
427	if (npt*nchord .gt. nspln) stop 'spline - NPT*NCHORD > 5000'
428c
429c	fill 4th coord with fraction of chain traced
430c
431	color_inc = 1.000 / float(npt)
432	fraction  = 0.0
433	if (ischeme.le.5) then
434	    do i = 1, npt
435		guide(4,i,1) = fraction
436		guide(4,i,2) = fraction
437		fraction = fraction + color_inc
438	    end do
439	endif
440c
441c	calculate spline segments
442c
443	tinc = 1./float(nchord)
444	do 1000 irib = 1, 2
445	iout = 1
446	do  900  ipt = 2, npt-1
447	t = 0.0
448	do i = 1, nchord
449	    iout = iout + 1
450	    call bspline( guide(1,ipt-1,irib), guide(1,ipt,irib),
451     1     		  guide(1,ipt+1,irib), t, smooth(1,iout,irib) )
452	    t = t + tinc
453	end do
454  900	continue
455 1000	continue
456c
457c	Add end segments (splines go midpoint-to-midpoint)
458c
459	iout = iout + 1
460	do 1100 irib = 1, 2
461	    do i = 1, 4
462		smooth(i, 1,    irib) = guide(i, 1,   irib )
463	        smooth(i, iout, irib) = guide(i, npt, irib )
464	    end do
465 1100	continue
466C
467    2	format(9f8.3,2x,3f5.2)
468    3	format('1',/,9f8.3,2x,3f5.2)
469c
470c	Start loop over spline segments
471c
472	ires = 1
473	jres = 1
474	kres = 2
475 2000	continue
476c	do 2100 ires = 1, iout-1
477		fraction =   smooth(4,ires,  1)
478c
479c	Make sure the two sides of the ribbon stay in register
480c
481		inext = ires + 1
482   55		dist0 = dist(smooth(1,inext,1),smooth(1,kres,2))
483		dist1 = dist(smooth(1,inext,1),smooth(1,kres+1,2))
484		if ((dist1 .lt. dist0) .and. (kres .lt. iout)) then
485			kres = kres + 1
486			goto 55
487		end if
488   56		dist0 = dist(smooth(1,inext,1),smooth(1,kres,2))
489		dist1 = dist(smooth(1,inext+1,1),smooth(1,kres,2))
490		if ((dist1 .lt. dist0) .and. (inext .lt. iout)) then
491			inext = inext + 1
492			goto 56
493		end if
494c
495		call colorit( color, fraction,
496     1  	  smooth(1,ires,1), smooth(1,jres,2), smooth(1,inext,1))
497c
498		if (mflag) then
499			write (output,3) (smooth(i,ires,  1),i=1,3),
500     1  		          (smooth(i,jres,  2),i=1,3),
501     2  		          (smooth(i,inext,1),i=1,3),
502     3  		          color
503 		else
504			write (output,2) (smooth(i,ires,  1),i=1,3),
505     1  		          (smooth(i,jres,  2),i=1,3),
506     2  		          (smooth(i,inext,1),i=1,3),
507     3  		          color
508 		endif
509c
510		if (jres .eq. kres) goto 2100
511		call colorit( color, fraction,
512     1  	  smooth(1,kres,2), smooth(1,inext,1), smooth(1,jres,2))
513		if (mflag) then
514			write (output,3) (smooth(i,jres,  2),i=1,3),
515     1  		          (smooth(i,inext,1),i=1,3),
516     2  		          (smooth(i,kres,  2),i=1,3),
517     3  		          color
518		else
519			write (output,2) (smooth(i,jres,  2),i=1,3),
520     1  		          (smooth(i,inext,1),i=1,3),
521     2  		          (smooth(i,kres,  2),i=1,3),
522     3  		          color
523		end if
524		jres = kres
525		if (kres .lt. iout) kres = kres + 1
526 2100	continue
527	ires = inext
528	if (ires .lt. iout) goto 2000
529c
530c	End loop over spline segments
531c
532	cindex = cindex + 1
533	return
534	end
535
536
537 	function dist(v1, v2)
538	real diff(3)
539	call vdif(diff,v1,v2)
540	dist = dot(diff,diff)
541	return
542	end
543
544	subroutine vload( v, s1, s2, s3 )
545	real v(3)
546	v(1) = s1
547	v(2) = s2
548	v(3) = s3
549	return
550	end
551
552	subroutine colorit( color, fraction, point1, point2, point3 )
553	real	color(3), point1(3), point2(3), point3(3)
554c
555c	   scheme 1	solid color (COLOR1)
556c	   scheme 2	shade from COLOR1 at 1st residue to COLOR2 at last
557c	   scheme 3	COLOR1 on front, COLOR3 (=COLOR2) on back
558c	   scheme 4	combination of 2 and 3 above
559c	   scheme 5	color each new chain a new color from RGB
560c
561        PARAMETER (MAXCOL=5000, MAXATM=10000)
562	common /COLORS/ ischeme, cindex, COLOR1(3), COLOR2(3), COLOR3(3)
563     &			,RGB(3,MAXCOL)
564 	integer		cindex
565	common /SPAM/ NATM, SPAM(4,MAXATM), SCAM(MAXATM)
566	integer SCAM
567	real	vec1(3), vec2(3), vec3(3)
568c
569	if ((ischeme .eq. 3) .or. (ischeme .eq. 4)) then
570		call vdif( vec1, point2, point1 )
571		call vdif( vec2, point3, point1 )
572		call cross( vec1, vec2, vec3 )
573		if (vec3(3) .lt. 0) then
574			color(1) = color3(1)
575			color(2) = color3(2)
576			color(3) = color3(3)
577		else if (ischeme .eq. 4) then
578			color(1) = fraction*color2(1)
579     &				 + (1.-fraction)*color1(1)
580	 		color(2) = fraction*color2(2)
581     &				 + (1.-fraction)*color1(2)
582			color(3) = fraction*color2(3)
583     &				 + (1.-fraction)*color1(3)
584		else
585			color(1) = color1(1)
586			color(2) = color1(2)
587			color(3) = color1(3)
588		end if
589	else if (ischeme .eq. 2) then
590		color(1) = fraction*color2(1) + (1.-fraction)*color1(1)
591		color(2) = fraction*color2(2) + (1.-fraction)*color1(2)
592		color(3) = fraction*color2(3) + (1.-fraction)*color1(3)
593	else if (ischeme .eq. 5) then
594		call vload( color,
595     &			    RGB(1,cindex), RGB(2,cindex), RGB(3,cindex))
596c	else if (ischeme .eq. 6) then
597c		ICOL = SCAM(fraction)
598c		color(1) = RGB(1,icol)
599c		color(2) = RGB(2,icol)
600c		color(3) = RGB(3,icol)
601	else
602		call vload( color, color1(1), color1(2), color1(3) )
603	end if
604	return
605	end
606
607	subroutine bspline( v1, v2, v3, t, v4 )
608	real v1(4), v2(4), v3(4)
609	real t
610	real v4(4)
611c
612	frac3 = 0.5 * t*t
613	frac1 = 0.5 * (1.-t) * (1.-t)
614	frac2 = 1. - (frac1 + frac3)
615	do i = 1, 4
616		v4(i) = frac1 * v1(i) + frac2 * v2(i) + frac3 * v3(i)
617	end do
618	return
619	end
620