;:ts=8

	BLANKS	ON

	far	code
	far	data

	include	"rays.i"

	cseg

	public	_ipattern

_ipattern:
	movem.l	d4-d7/a2,-(sp)
	movea.l	24(sp),a0		;address of patch
	movea.l	28(sp),a2		;address of object

	move.l	(a0)+,d0		;hit point
	move.l	(a0)+,d1
	move.l	(a0)+,d2

	lea	object+so_posi(a2),a0	;address of position of object
	sub.l	(a0)+,d0		;subtract from hit point
	sub.l	(a0)+,d1
	sub.l	(a0)+,d2

	movem.l	d0-d2,-(sp)

	lea	object+so_cvct(a2),a1	;transform to local coordinates
	move.l	sp,a0
	bsr	_ndot
	move.l	d0,d4			;Z component
	sub.l	op_para+24(a2),d4

	lea	object+so_avct(a2),a1
	move.l	sp,a0
	bsr	_ndot
	move.l	d0,d3			;X component
	sub.l	op_para+20(a2),d3

	lea	12(sp),sp		;fix stack

	move.l	op_para+12(a2),d7	;spacing
	beq	chk10

	move.w	op_para+32(a2),d0	;random number seed

	bsr	rnd
	move.w	d0,-(sp)
	move.l	d0,d1
	add.l	d1,d1
	move.l	d3,d0
	jsr	_muln32
	move.l	d7,d1
	jsr	_divs32
	jsr	_sin32
	move.l	d0,d1
	move.l	d7,d0
	jsr	_muln32
	asr.l	#2,d0
	move.l	d0,d5

	move.w	(sp)+,d0
	bsr	rnd
	move.w	d0,-(sp)
	move.l	d0,d1
	add.l	d1,d1
	move.l	d3,d0
	jsr	_muln32
	move.l	d7,d1
	jsr	_divs32
	jsr	_sin32
	move.l	d0,d1
	move.l	d7,d0
	jsr	_muln32
	asr.l	#2,d0
	move.l	d0,d6

	move.w	(sp)+,d0
	bsr	rnd
	move.w	d0,-(sp)
	move.l	d0,d1
	add.l	d1,d1
	move.l	d4,d0
	jsr	_muln32
	move.l	d7,d1
	jsr	_divs32
	jsr	_sin32
	move.l	d0,d1
	move.l	d7,d0
	jsr	_muln32
	asr.l	#2,d0
	add.l	d0,d5

	move.w	(sp)+,d0
	bsr	rnd
	move.w	d0,-(sp)
	move.l	d0,d1
	add.l	d1,d1
	move.l	d4,d0
	jsr	_muln32
	move.l	d7,d1
	jsr	_divs32
	jsr	_sin32
	move.l	d0,d1
	move.l	d7,d0
	jsr	_muln32
	asr.l	#2,d0
	add.l	d0,d6

	move.w	(sp)+,d0
	bsr	rnd
	move.w	d0,-(sp)
	move.l	d0,d1
	add.l	d1,d1
	move.l	d3,d0
	jsr	_muln32
	move.l	d7,d1
	jsr	_divs32
	jsr	_sin32
	move.l	d0,d1
	move.l	d7,d0
	jsr	_muln32
	asr.l	#2,d0
	add.l	d0,d5

	move.w	(sp)+,d0
	bsr	rnd
	move.l	d0,d1
	add.l	d1,d1
	move.l	d4,d0
	jsr	_muln32
	move.l	d7,d1
	jsr	_divs32
	jsr	_sin32
	move.l	d0,d1
	move.l	d7,d0
	jsr	_muln32
	asr.l	#2,d0
	add.l	d0,d6

	move.l	op_para+28(a2),d0
	move.l	d5,d1
	jsr	_muls32
	add.l	d0,d3

	move.l	op_para+28(a2),d0
	move.l	d6,d1
	jsr	_muls32
	add.l	d0,d4

	move.l	d3,d0
	move.l	d3,d1
	jsr	_mulr64
	move.l	d3,-(sp)
	move.l	d2,-(sp)
	move.l	d4,d0
	move.l	d4,d1
	jsr	_mulr64
	move.l	(sp)+,d0
	add.l	(sp)+,d3
	addx.l	d0,d2
	move.l	d3,-(sp)
	move.l	d2,-(sp)
	move.l	sp,a1
	jsr	_sqrt64
	addq.w	#8,sp

	move.l	d7,d1
	jsr	_divs32
	move.l	#411775,d1		;2*PI
	jsr	_muls32
	jsr	_sin32
	add.l	#$10000,d0
	asr.l	#1,d0
	move.w 	op_para+16(a2),d3
	subq.w	#1,d3
	bmi	lb2

	move.l	d0,d2
lb1	move.l	d2,d1
	jsr	_muln32
	dbra	d3,lb1

lb2	move.l	d0,d2
	move.l	#$10000,d3
	sub.l	d2,d3

	movea.l	24(sp),a0
	lea	ptc_col(a0),a0		;address to store color

	lea	op_para(a2),a1		;address of parameters 0-2 (color)

	move.l	(a1)+,d0
	asr.l	#8,d0
	move.l	d2,d1
	jsr	_muln32
	move.l	d0,d4
	move.l	(a0),d0
	move.l	d3,d1
	jsr	_muln32
	add.l	d4,d0
	move.l	d0,(a0)+		;red part

	move.l	(a1)+,d0
	asr.l	#8,d0
	move.l	d2,d1
	jsr	_muln32
	move.l	d0,d4
	move.l	(a0),d0
	move.l	d3,d1
	jsr	_muln32
	add.l	d4,d0
	move.l	d0,(a0)+		;green part

	move.l	(a1)+,d0
	asr.l	#8,d0
	move.l	d2,d1
	jsr	_muln32
	move.l	d0,d4
	move.l	(a0),d0
	move.l	d3,d1
	jsr	_muln32
	add.l	d4,d0
	move.l	d0,(a0)+		;blue part

chk10	movem.l	(sp)+,d4-d7/a2
	rts

;
;	multiply vector by unit normal (a0).(a1) = d0   (dot product)
;	- unit normal in (a1)
;

_ndot:
	move.l	(a0)+,d0
	move.l	(a1)+,d1
	bsr	_muln32
	move.l	d0,-(sp)

	move.l	(a0)+,d0
	move.l	(a1)+,d1
	bsr	_muln32
	add.l	d0,(sp)

	move.l	(a0),d0
	move.l	(a1),d1
	bsr	_muln32
	add.l	(sp)+,d0

	rts

_mulr64:
	tst.l	d1
	bpl.s	r6401

	neg.l	d1

	tst.l	d0
	bpl.s	r6402

	neg.l	d0

r6400	move.l	d0,d2
	swap	d2
	move.w	d0,d3
	mulu	d1,d3		;d0.L * d1.L
	move.w	d1,-(sp)
	swap	d1
	mulu	d1,d0		;d0.L * d1.H
	mulu	d2,d1		;d0.H * d1.H
	mulu	(sp)+,d2	;d1.L * d0.H
	add.l	d0,d2
	move.w	d2,d0
	clr.w	d2
	swap	d2
	swap	d0
	clr.w	d0
	add.l	d0,d3
	addx.l	d1,d2
	rts

r6401	tst.l	d0
	bpl	r6400

	neg.l	d0

r6402	move.l	d0,d2
	swap	d2
	move.w	d0,d3
	mulu	d1,d3		;d0.L * d1.L
	move.w	d1,-(sp)
	swap	d1
	mulu	d1,d0		;d0.L * d1.H
	mulu	d2,d1		;d0.H * d1.H
	mulu	(sp)+,d2	;d1.L * d0.H
	add.l	d0,d2
	move.w	d2,d0
	clr.w	d2
	swap	d2
	swap	d0
	clr.w	d0
	add.l	d0,d3
	addx.l	d1,d2
	neg.l	d3
	negx.l	d2
	rts

;
;	multiply d0 by "unit number" in d1, leave result is d0
;

_muln32:
	tst.l	d1
	bpl.s	n3203

	neg.l	d1

	btst.l	#16,d1
	beq	n3201

	neg.l	d0

n3200	rts

n3201	tst.l	d0
	bpl.s	n3204

	neg.l	d0

n3202	move.w	d0,-(sp)
	swap	d0
	mulu	d1,d0
	mulu	(sp)+,d1
	add.w	d1,d1
	clr.w	d1
	swap	d1
	addx.l	d1,d0
	rts

n3203	btst.l	#16,d1
	bne.s	n3200

	tst.l	d0
	bpl	n3202

	neg.l	d0

n3204	move.w	d0,-(sp)
	swap	d0
	mulu	d1,d0
	mulu	(sp)+,d1
	add.w	d1,d1
	clr.w	d1
	swap	d1
	addx.l	d1,d0
	neg.l	d0
	rts

;
;	muls32 - multiply d0 by d1, leave result is d0
;

_muls32:
	tst.l	d1
	bpl.s	m3201

	neg.l	d1

	tst.l	d0
	bpl.s	m3202

	neg.l	d0

m3200	move.l	d0,-(sp)
	move.l	d1,-(sp)
	mulu.w	d0,d1		;d0.L * d1.L
	swap	d0
	mulu.w	(sp),d0		;d0.H * d1.H
	swap	d0
	clr.w	d0
	add.w	d1,d1		;roundoff bit
	swap	d1
	addx.w	d1,d0
	move.w	(sp)+,d1	;d1.H
	mulu.w	4(sp),d1	;d0.L * d1.H
	add.l	d1,d0
	move.w	(sp)+,d1	;d1.L
	mulu.w	(sp),d1		;d0.H * d1.L
	add.l	d1,d0
	addq.w	#4,sp
	rts

m3201	tst.l	d0
	bpl	m3200

	neg.l	d0

m3202	move.l	d0,-(sp)
	move.l	d1,-(sp)
	mulu.w	d0,d1		;d0.L * d1.L
	swap	d0
	mulu.w	(sp),d0		;d0.H * d1.H
	swap	d0
	clr.w	d0
	add.w	d1,d1		;roundoff bit
	swap	d1
	addx.w	d1,d0
	move.w	(sp)+,d1	;d1.H
	mulu.w	4(sp),d1	;d0.L * d1.H
	add.l	d1,d0
	move.w	(sp)+,d1	;d1.L
	mulu.w	(sp),d1		;d0.H * d1.L
	add.l	d1,d0
	neg.l	d0
	addq.w	#4,sp
	rts

;
;	_muli	- multiply integer in d0 by FRACT in d1, result in d0
;

_muli:
	tst.w	d0
	bpl	mli2

	neg.w	d0
	move.w	d1,-(sp)
	swap	d1
	mulu	d0,d1
	mulu	(sp)+,d0
	swap	d1
	clr.w	d1
	add.l	d1,d0
	neg.l	d0
	rts

mli2	move.w	d1,-(sp)
	swap	d1
	mulu	d0,d1
	mulu	(sp)+,d0
	swap	d1
	clr.w	d1
	add.l	d1,d0
	rts

;
;	divs32 - divides two 32 fractional integers d0 / d1 = d0
;
;	note: dividing by very small numbers can result in an infinite
;	result. -oo = 8000.0000, oo = 7fff.ffff
;

_divs32:
	movem.l	d2-d3,-(sp)

	tst.l	d1			;test d1
	beq	d3231			;divide by zero?
	bpl	d3201			;is it positive?

	neg.l	d1			;negate d1
	neg.l	d0			;and d0

d3201	move.l	d0,d3			;test d0 - save sign in d3
	bpl	d3202

	neg.l	d0			;negate d1

d3202	cmp.l	d1,d0			;branch if d0 >= d1
	bcc.s	d3203

	moveq	#-1,d2			;init quotient to zero
	move.w	#15,d3			;16 bits to do (since d1 > d0)

d3221	add.l	d0,d0			;shift d0 left until it's >= d1
	cmp.l	d1,d0			;(zero bits in quotient's LSB word)
	dbcc	d3,d3221		;decr bit count all but last time

	bcc.s	d3205			;do rest of loop if d0 >= d1
					;d3 is number of bits to do minus one

	moveq	#0,d0			;underflow if d3 counted out
	bra	d3207

d3203	moveq	#14,d2			;start with 15 leading zero bits
					; in quotient MSB word (+ sign bit)

d3204	add.l	d1,d1			;double d1 until d1 >= d0
	cmp.l	d1,d0
	dbcs	d2,d3204		;decr # of zeros all but last time

	bcs	d3232			;branch if d1 >= d0 now

d3231	move.l	#$7fffffff,d0		;overflow if # of zeros hits -1
	bra	d3207

d3232	lsr.l	#1,d1			;shift d1 to the right so d0 >= d1
	move.w	#30,d3			;total number of bits to do is
	sub.w	d2,d3			; 31 minus # of leading zero bits
	moveq	#-1,d2			;init quotient to zero

d3205	sub.l	d1,d0			;is d1 <= "remainder"
	bcc.s	d3206
	add.l	d1,d0			;no - fix up d0 again
d3206	addx.l	d2,d2			;shift in quotient bit (inverted)
	add.l	d0,d0			;double remainder
	dbeq	d3,d3205		;loop while bits left to do and
					; remainder is non-zero

	bne	d3261			;branch if didn't get zero remainder

	rol.l	d3,d2			;else shift in 'd3' zero bits

d3261	not.l	d2			;(un)invert quotient bits
	move.l	d2,d0

d3207	tst.l	d3			;was this a positive quotient
	bpl	d3208

	neg.l	d0			;negate quotient

d3208	movem.l (sp)+,d2-d3
	rts

;
;	sqrt64 - finds the square root of a 64 bit fractional integer
;		 sqrt (a1) = d0
;

_sqrt64:
	move.l	d4,-(sp)

	move.l	(a1)+,d2
	move.l	(a1),d3

	move.l	#$40000000,d0
	move.l	#$c0000000,d1

	moveq	#30,d4

sq_1	cmp.l	d0,d2
	bcs.s	sq_2

	sub.l	d0,d2
	or.l	d1,d0

sq_2	add.l	d3,d3
	addx.l	d2,d2
	lsr.l	#1,d1
	eor.l	d1,d0
	dbra	d4,sq_1

	sub.l	#$80000000,d3
	sub.l	d0,d2
	bcs.s	sq_3

	or.l	d1,d0

sq_3	movem.l	(sp)+,d4

	rts

rnd:	mulu	#$41c6,d0
	add.w	#$3039,d0
	and.l	#$7fff,d0
	rts

;
;	sin32	- takes the sin of a FRACT (radians) angle in d0
;

pi2	equ	411775
pi	equ	205887
pid2	equ	102944

_sin32	move.l	d0,-(sp)
	move.l	#pi2,d1
	jsr	_divs32
	clr.w	d0
	swap	d0
	move.l	#pi2,d1
	jsr	_muli
	move.l	(sp)+,d1
	sub.l	d0,d1		;0 <= d1 < 2*PI

	moveq	#0,d0
	cmp.l	#pi,d1
	ble	sn3201
	sub.l	#pi,d1
	moveq	#-1,d0

sn3201	cmp.l	#pid2,d1
	ble	sn3202
	sub.l	#pi,d1
	neg.l	d1

sn3202	move.w	d0,-(sp)	;sign
	move.l	d1,-(sp)	;x

	move.l	#$10000,-(sp)	;sin(x)/x

	move.l	d1,d0
	jsr	_muls32
	move.l	d0,-(sp)	;x*x

	move.l	#10923,d1	;1/(2*3)
	jsr	_muln32
	sub.l	d0,4(sp)

	move.l	d0,d1		;x*x/2/3
	move.l	(sp),d0		;x*x
	jsr	_muln32
	move.l	#3277,d1	;1/(4*5)
	jsr	_muln32
	add.l	d0,4(sp)

	move.l	d0,d1		;x*x*x*x/2/3/4/5
	move.l	(sp),d0		;x**2
	jsr	_muln32
	move.l	#1560,d1	;1/(6*7)
	jsr	_muln32
	sub.l	d0,4(sp)

	move.l	d0,d1		;x*x*x*x*x*x/2/3/4/5/6/7
	move.l	(sp)+,d0	;x**2
	jsr	_muln32
	move.l	#910,d1		;1/(8*9)
	jsr	_muln32
	add.l	(sp)+,d0

	move.l	d0,d1
	move.l	(sp)+,d0
	jsr	_muln32

	move.w	(sp)+,d1
	bpl	sn3203

	neg.l	d0

sn3203	rts