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