1************************************************************************
2*              Support routines for quadric surfaces                   *
3************************************************************************
4* EAM Jun 1997	- initial version, supports version 2.4(alpha) of render
5* EAM May 1998	- additional error checking to go with Parvati/rastep
6* EAM Jan 1999	- version 2.4i
7* EAM Mar 2008	- Gfortran optimization breaks the object accounting.
8*		  No real solution yet.  Move qinp.f into separate file
9************************************************************************
10*
11* Quadric surfaces include spheres, cones, ellipsoids, paraboloids, and
12* hyperboloids.  The motivation for this code was to allow rendering
13* thermal ellipsoids for atoms, so the other shapes have not been
14* extensively tested.
15* A quadric surface is described by 10 parameters (A ... J).
16* For efficiency during rendering it is also useful to know the center and
17* a bounding sphere. So a QUADRIC descriptor to render has 17 parameters:
18* 14 (object type QUADRIC)
19*    X Y Z RADLIM RED GRN BLU
20*    A B C D E F G H I J
21*
22* The surface itself is the set of points for which Q(x,y,z) = 0
23* where
24*      Q(x,y,z)	=  A*x^2 +  B*y^2 +  C*z^2
25*         	+ 2D*x*y + 2E*y*z + 2F*z*x
26*         	+ 2G*x   + 2H*y   + 2I*z
27*         	+  J
28*
29* It is convenient to store this information in a matrix QQ
30*	| QA QD QF QG |		QA = A    QB = B    QC = C
31* QQ =	| QD QB QE QH |		QD = D    QE = E    QF = F
32*	| QF QE QC QI |		QG = G    QH = H    QI = I
33*	| QG QH QI QJ |		QJ = J
34*
35* Then Q(x,y,x) = XT*QQ*X	where X = (x,y,z,1)
36* The point of this is that a 4x4 homogeneous transformation T can be
37* applied to QQ by matrix multiplication:     QQ' = TinvT * QQ * Tinv
38*
39* The surface normal is easily found by taking the partial derivatives
40* of Q(x,y,z) at the point of interest.
41************************************************************************
42* TO DO:
43*	- can we distinguish an ellipsoid from other quadrics? Do we care?
44************************************************************************
45
46
47CCC	Calculate the matrices for quadric transformation
48CC	QQ' = TINVT * QQ * TINV      (TINV is inverse of transposed TMAT)
49C
50	subroutine qsetup
51*
52      COMMON /MATRICES/ XCENT, YCENT, SCALE, EYEPOS, SXCENT, SYCENT,
53     &                  TMAT, TINV, TINVT, SROT, SRTINV, SRTINVT
54     &                 ,RAFTER, TAFTER
55      REAL   XCENT, YCENT, SCALE, EYEPOS, SXCENT, SYCENT
56      REAL   TMAT(4,4), TINV(4,4),   TINVT(4,4)
57      REAL   SROT(4,4), SRTINV(4,4), SRTINVT(4,4)
58      REAL   RAFTER(4,4), TAFTER(3)
59*
60      COMMON /ASSCOM/ noise, verbose
61      integer         noise
62      logical         verbose
63*
64	real   TRAN(4,4), POST(4,4), det
65*
66*	TMAT is a post-multiplier matrix, but unfortunately the
67*	quadric surface math was worked out for pre-multipliers
68*	So the quadric math uses the transpose of TMAT
69	call trnsp4( TRAN, TMAT )
70*
71*	Remove translation component from TMAT.
72*	I had to make a choice whether the X,Y,Z "center" of the quadric
73*	is implicitly coded into the coefficients, or whether it is to be
74*	maintained explicitly.  Since all other object types do the latter,
75*	I have chosen to keep all translation components out of the quadric
76*	coefficients
77	tran(1,4) = 0.0
78	tran(2,4) = 0.0
79	tran(3,4) = 0.0
80*
81*	July 1999 - Allow post-hoc rotation matrix also
82	call tmul4( POST, RAFTER, TRAN )
83*
84	det = tinv4( TINV, POST ) / TMAT(4,4)
85	call trnsp4( TINVT, TINV )
86*
87*	While we're at it, check for legality of rotation matrix
88	if (abs(1. - abs(det)) .gt. 0.02) then
89	    write (noise,901) abs(det)
90  901	    format('>>> Warning: Determinant of rotation matrix =',
91     &           F7.3,' <<<')
92	endif
93	if (TMAT(1,4).ne.0.or.TMAT(2,4).ne.0.or.TMAT(3,4).ne.0) then
94	    write (noise,902)
95  902       format('>>> Warning: Non-zero cross terms in 4th column of',
96     &             ' TMAT <<<')
97	endif
98*
99*	Same thing again for shadow rotation matrix
100	call  trnsp4( TRAN, SROT )
101	det = tinv4( SRTINV, TRAN )
102	call  trnsp4( SRTINVT, SRTINV )
103	if (abs(1. - abs(det)) .gt. 0.02) then
104	    write (noise,903) abs(det)
105  903	    format('>>> Warning: Determinant of shadow matrix =',
106     &           F7.3,' <<<')
107	endif
108*
109	return
110	end
111
112
113************************************************************************
114* This routine does the exact test for pixel impingement of a quadric  *
115* in both pixel and shadow space.                                      *
116* Find Z coordinate of point on quadric surface with given X, Y        *
117* Also return surface normal at that point                             *
118************************************************************************
119CCC
120CC
121C
122	function qtest( center, coeffs, xp, yp, zp, qnorm, shadowspace,
123     &                  backside )
124	IMPLICIT NONE
125	logical  qtest, shadowspace
126	logical  backside
127	real    center(3), coeffs(10), xp, yp, zp, qnorm(3)
128*
129      COMMON /MATRICES/ XCENT, YCENT, SCALE, EYEPOS, SXCENT, SYCENT,
130     &                  TMAT, TINV, TINVT, SROT, SRTINV, SRTINVT
131     &                 ,RAFTER, TAFTER
132      REAL   XCENT, YCENT, SCALE, EYEPOS, SXCENT, SYCENT
133      REAL   TMAT(4,4), TINV(4,4),   TINVT(4,4)
134      REAL   SROT(4,4), SRTINV(4,4), SRTINVT(4,4)
135      REAL   RAFTER(4,4), TAFTER(3)
136*
137	real    QA,QB,QC,QD,QE,QF,QG,QH,QI,QJ
138	real    AA,BB,CC,DD,EE
139	real    x,y,z
140*
141	QA = coeffs(1)
142	QB = coeffs(2)
143	QC = coeffs(3)
144	QD = coeffs(4)
145	QE = coeffs(5)
146	QF = coeffs(6)
147	QG = coeffs(7)
148	QH = coeffs(8)
149	QI = coeffs(9)
150	QJ = coeffs(10)
151*
152*	xp and yp are in pixel coordinates, rather than in normalized ones
153*	First displace center of quadric surface from the implicit origin
154*	to the point specified by center(3), then convert to pixel coords
155*
156	x = (xp - center(1)) / scale
157	y = (yp - center(2)) / scale
158*
159	AA = QC
160	BB = 2.0 * (QF*x + QE*y + QI)
161	CC = QA*x*x + QB*y*y + 2.0*(QD*x*y + QG*x + QH*y) + QJ
162	DD = BB*BB - 4*AA*CC
163	if (DD.LT.0) then
164	    qtest = .false.
165	    return
166	else
167	    qtest = .true.
168	    EE = sqrt(DD)
169	endif
170	if (AA .eq. 0.) then
171CDEBUG	    Can this happen????
172	    z = 9999.
173	else if (backside) then
174	    if (AA .le. 0.) then
175	    	z = (-BB + EE) / (2*AA)
176	    else
177	    	z = (-BB - EE) / (2*AA)
178	    endif
179	else
180	    if (AA .gt. 0.) then
181	    	z = (-BB + EE) / (2*AA)
182	    else
183	    	z = (-BB - EE) / (2*AA)
184	    endif
185	endif
186	zp = z * scale
187	zp = zp + center(3)
188*
189*	Surface normal comes from partial derivatives at this point
190	if (.not. shadowspace) then
191	    qnorm(1) = QA*x + QD*y + QF*z + QG
192	    qnorm(2) = QB*y + QD*x + QE*z + QH
193	    qnorm(3) = QC*z + QF*x + QE*y + QI
194	endif
195*
196	return
197	end
198
199CCC	Convert ANISOU description of anisotropic displacement parameters
200CC	into quadric surface enclosing a given probability volume
201C	Returns -1 if the Uij are non-positive definite, 1 otherwise
202C
203	function anitoquad( ANISOU, PROB, QUADRIC, EIGENS, EVECS)
204	implicit NONE
205	integer  anitoquad
206	real     ANISOU(6), PROB, QUADRIC(10)
207	real     EIGENS(3), EVECS(4,4)
208c
209	COMMON /ASSCOM/ noise, verbose
210	integer         noise
211	logical         verbose
212c
213	real	UU(4,4), UINV(4,4)
214	integer i
215c
216	real     tinv4, det
217c
218c	Save FPU traps later by checking that ANISOU is non-zero
219	do i = 1,3
220	    EIGENS(i) = 0.0
221	end do
222	if (ANISOU(1).eq.0.and.ANISOU(2).eq.0.and.ANISOU(3).eq.0 .and.
223     &      ANISOU(4).eq.0.and.ANISOU(5).eq.0.and.ANISOU(6).eq.0) then
224     		anitoquad = -2
225		return
226	end if
227c
228c	Build matrix from Uij coefficients and invert it
229	UU(1,1) = ANISOU(1)
230	UU(2,2) = ANISOU(2)
231	UU(3,3) = ANISOU(3)
232 	UU(4,4) = -(1/PROB**2)
233	UU(1,2) = ANISOU(4)
234	UU(2,1) = ANISOU(4)
235	UU(2,3) = ANISOU(6)
236	UU(3,2) = ANISOU(6)
237	UU(1,3) = ANISOU(5)
238	UU(3,1) = ANISOU(5)
239	UU(1,4) = 0.0
240	UU(4,1) = 0.0
241	UU(2,4) = 0.0
242	UU(4,2) = 0.0
243	UU(3,4) = 0.0
244	UU(4,3) = 0.0
245	det = tinv4( UINV, UU )
246c
247c	and from that we extract the coefficients of the surface
248	QUADRIC(1)  = UINV(1,1)
249	QUADRIC(2)  = UINV(2,2)
250	QUADRIC(3)  = UINV(3,3)
251	QUADRIC(4)  = UINV(1,2)
252	QUADRIC(5)  = UINV(2,3)
253	QUADRIC(6)  = UINV(1,3)
254	QUADRIC(7)  = UINV(1,4)
255	QUADRIC(8)  = UINV(2,4)
256	QUADRIC(9)  = UINV(3,4)
257	QUADRIC(10) = UINV(4,4)
258c
259c	Find eigenvalues of the ellipsoid
260	call jacobi( UU,   3, 4, EIGENS, EVECS )
261c
262c	Units of this matrix are A**2; we want values in A
263	anitoquad = 1
264	do i = 1,3
265	    if (EIGENS(i).gt.0.) then
266		EIGENS(i) = sqrt(EIGENS(i))
267	    else
268C		write (noise,*) 'Non-positive definite ellipsoid!'
269		EIGENS(i) = 0.
270		anitoquad = -1
271	    endif
272	enddo
273c
274	return
275	end
276
277CCC	Matrix manipulation routines for 4x4 homogeneous transformations
278CC
279C
280	subroutine tmul4( C, A, B )
281	real   C(4,4), A(4,4), B(4,4)
282	integer i,j
283	do i = 1,4
284	do j = 1,4
285	    C(i,j) = A(i,1)*B(1,j) + A(i,2)*B(2,j)
286     &		   + A(i,3)*B(3,j) + A(i,4)*B(4,j)
287	enddo
288	enddo
289	return
290	end
291
292	subroutine trnsp4( B, A )
293	real   B(4,4), A(4,4)
294	integer i,j
295	do i=1,4
296	do j=1,4
297	    B(i,j) = A(j,i)
298	enddo
299	enddo
300	return
301	end
302
303CCC	Matrix inversion for a 4x4 matrix A
304CC
305C
306	function tinv4( AI, A )
307	real     tinv4
308	real     AI(4,4), A(4,4)
309c
310	real	TMP(4,4), D
311	integer index(4)
312c
313	do i=1,4
314	    do j=1,4
315		TMP(i,j) = A(i,j)
316		AI(i,j) = 0.
317	    enddo
318	    AI(i,i) = 1.
319	enddo
320	call ludcmp( TMP, 4, index, D )
321	tinv4 = D * TMP(1,1)*TMP(2,2)*TMP(3,3)*TMP(4,4)
322	do j=1,4
323	    call lubksb( TMP, 4, index, AI(1,j) )
324	enddo
325	return
326	end
327
328************************************************************************
329*              Matrix inversion via LU decomposition                   *
330*	adapted from Numerical Recipes in Fortran (1986)               *
331************************************************************************
332*
333CCC	input  NxN matrix A is replaced by its LU decomposition
334CC	output index(N) records row permutation due to pivoting
335C	output D is parity of row permutations
336c
337	subroutine ludcmp( A, n, index, D )
338	implicit  NONE
339	integer   n,index(n)
340	real      A(n,n)
341	real      D
342c
343	integer   i,imax,j,k
344	real      aamax,dum,sum
345	integer    NMAX
346	parameter (NMAX=10)
347	real      vv(NMAX)
348c
349	d = 1.
350	do i=1,n
351	    aamax = 0.
352	    do j=1,n
353	    	if (abs(A(i,j)).gt.aamax) aamax = abs(A(i,j))
354	    enddo
355	    call assert(aamax.ne.0.,'Singular matrix')
356	    vv(i) = 1. / aamax
357	enddo
358c
359	imax = 1
360	do j=1,n
361	    if (j.gt.1) then
362		do i=1,j-1
363		    sum = A(i,j)
364		    if (i.gt.1) then
365			do k=1,i-1
366			    sum = sum - A(i,k)*A(k,j)
367			enddo
368			A(i,j) = sum
369		    endif
370		enddo
371	    endif
372	    aamax = 0.
373	    do i=j,n
374		sum = A(i,j)
375		if (j.gt.1) then
376		    do k=1,j-1
377			sum = sum - A(i,k)*A(k,j)
378		    enddo
379		    A(i,j) = sum
380		endif
381		dum = vv(i) * abs(sum)
382		if (dum.ge.aamax) then
383		    imax = i
384		    aamax = dum
385		endif
386	    enddo
387	    if (j.ne.imax) then
388		do k=1,n
389		    dum = A(imax,k)
390		    A(imax,k) = A(j,k)
391		    A(j,k) = dum
392		enddo
393		d = -d
394		vv(imax) = vv(j)
395	    endif
396	    index(j) = imax
397	    call assert(A(j,j).ne.0.,'Singular matrix')
398	    if (j.ne.n) then
399		dum = 1. / A(j,j)
400		do i=j+1,n
401		    A(i,j) = A(i,j) * dum
402		enddo
403	    endif
404	enddo
405	return
406	end
407
408CCC	corresponding back-substitution routine
409CC
410C
411	subroutine lubksb( A, N, index, B )
412	implicit NONE
413	integer  n, index(n)
414	real     A(n,n), B(n)
415c
416	integer  i,ii,j,ll
417	real     sum
418c
419	ii = 0
420	do i=1,n
421	    ll = index(i)
422	    sum = B(ll)
423	    B(ll) = B(i)
424	    if (ii.ne.0) then
425		do j=ii,i-1
426		    sum = sum - A(i,j)*B(j)
427		enddo
428	    else if (sum.ne.0.) then
429		ii = i
430	    endif
431	    B(i) = sum
432	enddo
433c
434	do i=n,1,-1
435	    sum = B(i)
436	    if (i.lt.n) then
437		do j=i+1,n
438		    sum = sum - A(i,j)*B(j)
439		enddo
440	    endif
441	    B(i) = sum / A(i,i)
442	enddo
443	return
444	end
445
446************************************************************************
447* Find eigenvalues and eigenvectors of nxn symmetric matrix using      *
448* cyclic Jacobi method. NP is (Fortran) physical storage dimension.    *
449* On return A is destroyed, D contains eigenvalues,and each column of  *
450* V is a normalized eigenvector                                        *
451* Adapted from Numerical Recipes in Fortran (1986)                     *
452* We only need it for 3x3 symmetric matrices, so it's overkill.        *
453************************************************************************
454CCC
455CC
456C
457	subroutine jacobi( A, n, np, D, V )
458	PARAMETER (NMAX=4)
459c	Machine dependent! (converge when off-diagonal sum is less than this)
460	PARAMETER (TINY = 1.e-37)
461	real A(np,np), D(n), V(np,np)
462	real B(NMAX), Z(NMAX)
463c
464	call assert(n.le.NMAX,'Matrix too big for eigenvector routine')
465c
466c	Initialize V to identity matrix, B and D to diagonal of A
467	do ip = 1,n
468	do iq = 1,n
469	    V(ip,iq) = 0.0
470	enddo
471	V(ip,ip) = 1.0
472	B(ip) = A(ip,ip)
473	D(ip) = B(ip)
474	Z(ip) = 0.0
475	enddo
476	nrot = 0
477c
478c	50 sweeps is never expected to happen; Press et al (1986) claim
479c	6 to 10 are typical for moderate matrix sizes
480c	Empirical trials of rastep show 6 sweeps, 8-10 rotations
481	do i = 1,50
482	    sm = 0.0
483	    do ip = 1,n-1
484	    do iq = ip+1,n
485		sm = sm+abs(A(ip,iq))
486	    enddo
487	    enddo
488	    if (sm .lt. TINY) return
489c	    After 4 sweeps skip rotation if the off-diagonal is small
490	    if (i .lt. 4) then
491		thresh = 0.2*sm/(n*n)
492	    else
493		thresh = 0.0
494	    endif
495	    do ip = 1,n-1
496	    do iq = ip+1,n
497		g = 100. * abs(A(ip,iq))
498		if ((i.gt.4) .and.
499     &		    (abs(D(ip)) + g .eq. abs(D(ip))) .and.
500     &		    (abs(D(iq)) + g .eq. abs(D(iq)))) then
501		    A(ip,iq) = 0.0
502		else if (abs(A(ip,iq)).gt.thresh) then
503		    h = D(iq) - D(ip)
504		    if (abs(h) + g .eq. abs(h)) then
505			t = A(ip,iq) / h
506		    else
507			theta = 0.5 * h / A(ip,iq)
508			t = 1. / (abs(theta) + sqrt(1.+theta*theta))
509			if (theta.lt.0.) t = -t
510		    endif
511		    c = 1./sqrt(1.+t*t)
512		    s = t * c
513		    tau = s/(1.+c)
514		    h = t * A(ip,iq)
515		    Z(ip) = Z(ip) - h
516		    Z(iq) = Z(iq) + h
517		    D(ip) = D(ip) - h
518		    D(iq) = D(iq) + h
519		    A(ip,iq) = 0.0
520		    do j = 1,ip-1
521			g = A(j,ip)
522			h = A(j,iq)
523			A(j,ip) = g - s*(h+g*tau)
524			A(j,iq) = h + s*(g-h*tau)
525		    enddo
526		    do j = ip+1,iq-1
527			g = A(ip,j)
528			h = A(j, iq)
529			A(ip,j) = g - s*(h+g*tau)
530			A(j,iq) = h + s*(g-h*tau)
531		    enddo
532		    do j = iq+1,n
533			g = A(ip,j)
534			h = A(iq,j)
535			A(ip,j) = g - s*(h+g*tau)
536			A(iq,j) = h + s*(g-h*tau)
537		    enddo
538		    do j = 1,n
539			g = V(j,ip)
540			h = V(j,iq)
541			V(j,ip) = g - s*(h+g*tau)
542			V(j,iq) = h + s*(g-h*tau)
543		    enddo
544		    nrot = nrot + 1
545		endif
546	    enddo
547	    enddo
548c
549c	    Update D with sum and reinitialize Z
550	    do ip = 1,n
551		B(ip) = B(ip) + Z(ip)
552		D(ip) = B(ip)
553		Z(ip) = 0.0
554	    enddo
555	enddo
556c
557	call assert(.false.,'Failed to find eigenvectors')
558	return
559	end
560c
561
562
563C	This one really has nothing to do with quadrics per se,
564C	but since it's invoked by qinp, both render and rastep
565C	need to be able to see it.
566C	Should really be in separate file of support routines
567C
568	FUNCTION PERSP( Z )
569	REAL PERSP, Z
570	COMMON /MATRICES/ XCENT, YCENT, SCALE, EYEPOS, SXCENT, SYCENT,
571     &                  TMAT, TINV, TINVT, SROT, SRTINV, SRTINVT
572     &                 ,RAFTER, TAFTER
573	REAL   XCENT, YCENT, SCALE, EYEPOS, SXCENT, SYCENT
574	REAL   TMAT(4,4), TINV(4,4),   TINVT(4,4)
575	REAL   SROT(4,4), SRTINV(4,4), SRTINVT(4,4)
576	REAL   RAFTER(4,4), TAFTER(3)
577	IF (Z/EYEPOS .GT. 0.999) THEN
578	    PERSP = 1000.
579	ELSE
580	    PERSP = 1. / (1. - Z/EYEPOS)
581	ENDIF
582	RETURN
583	END
584
585