to pretzel3
	;Pretzel by Mike Sandy
	;NOTE DRAW USES 'THROW' TO TRAP A PARTICULAR ERROR

	;reset
	setsc 0
	perspective
	cs
	ht

	catch "stopping
	[	comment 
		[	rs AND rl SHOULD BE SMALL, RELATIVELY PRIME INTEGERS. DO NOT PUT rs = rl
			fh (0 - 1) SIZE OF h RELATIVE TO rs
			size CONTROLS SIZE OF PLOT
			n NUMBER OF SEGMENTS/CYCLE FOR WHOLE CURVE
			lor (values 0 or 1) DETERMINES WHETHER KNOT IS LEFT OR RIGHT HANDED
		]	

		;PRETZEL
		draw [0 360*rs 1]
			[size1*( dr*(cos u)+h*cos (dr/rs*u) )]
			[size1*( dr*(sin u)-h*sin (dr/rs*u) )]
			[size1*sin (fr*u+180*lor)]
			[
				rs 2 rl_ 3 n 45 size 150
				lor 1  h 0.2*rs dr rl_-rs size1 size/dr
				fr rl_/rs
			]
			[offset [0 0 0] eyepos [0 0 700]]
			[bradius 0.25*size b_incr 15]
	
		(rotateScene 1)
		cs
	
		;AMMONITE!
		draw[-9 2*pi*n 0.05]
			[size*c*(exp a*u)*((exp b)+1)*radsin u]
			[size*c*(exp a*u)*((exp b)+1)*radcos u]
			[0]
			[n 2.2 size 204 c 0.05 a 0.1 b a*2*pi]
			[offset [-12 0 0] eyepos [0 -800 1000]]
			[bradius size*c*(exp a*u)*((exp b)-1) b_incr 15]
	
		rotateScene
		cs
	
		;CONICAL SPIRAL
		draw[0 2*pi*n 0.05]
			[size*a*u/k*radsin u]
			[size*a*u/k*radcos u]
			[size*u/k]
			[n 2.9 size 57 a 0.4 c 1 k 2*pi*c]
			[offset [0 0 -95] eyepos [0 1000 100]]
			[bradius 1.5*u b_incr 15]
	
		rotateScene
		cs
	
		;7-KNOT FROM LISSAJOUS CURVE
		draw[-0.2 2*pi*k 0.01]
			[size*radsin (a*u+b*pi)]
			[size*radsin u]
			[size*0.5*radcos u*7/3]
			[size 200 k 3 a 2/k b 3/k ]
			[offset [0 0 0] eyepos [0 0 -1000]]
			[bradius 0.15*size b_incr 15]
	
		(rotateScene 3)
	]
end

to draw  urange_l  xexpr_l yexpr_l zexpr_l  fnpars_l plotpars_l  bradius_l

	comment
	[	urange_l - RANGE LIST FOR THE VARIABLE u - [START END INCR]
		u CAN ONLY APPEAR IN THE PARAMETRIC EQUATIONS AND BRADIUS VALUE
		xexpr_l etc. - PARAMETRIC EQUNS
		fnpars - AS LIST [var1 val1 var2 val2.. ]
		A VARIABLE MAY BE DEFINED IN TERMS OF PREVIOUSLY DEFINED VARIABLES/PARAMETERS
		BUT u IS NOT INCLUDED ANYWHERE IN THIS LIST
		plotpars_l - [offset [  ] eyepos [  ]]
		eyepos TURTLE -1 POSITION
		bradius - RADIUS OF BAND. u CAN BE INCLUDED IN ITS VALUE
		bradius AND b_incr MUST BE DECLARED
		b_incr band incremental angle in degreeS

		EXCEPT IN THE CASE OF A SIMPLE CURVE PLOT
		IF BRADIUS IS TOO SMALL COMPARED TO U THE PLOT IS STOPPED
	]
	
	clearShadows
	enableShadows
	local 
	[	eyepos offset 
		b_incr br_l bradius 
		u_start u u_incr u_end 
		iposn iposn0 
		edge2_l edge1_l 
		fl
		timestart
	]
	timestart=timefine
	eyepos=[0 0 1000]		;SETS DEFAULT VALUE
	offset=[0 0 0]			;SETS DEFAULT VALUE

	b_incr=run (list last bradius_l)		;ASSIGNS BAND INCREMENT
	br_l=butlast butlast bf bradius_l           ;BR_L bradius FORMULA
	assign_val fnpars_l assign_val plotpars_l             ;ASSIGN VARIABLE VALUES
	u_start=run (list first urange_l)
	u=u_start                                   ;SET U TO START
	urange_l=bf urange_l                             ;REMOVE 
	u_incr=run (list last urange_l)            ;SET UP U INCREMENT
	u_end=run butlast urange_l                  ;SET UP U_END
	iposn=offset+(list x y z)            ;STORES FIRST CURVE POINT
	iposn0=[]
	edge2_l=[] 
	edge1_l=[]             ;STORES FOR BAND EDGE COORD
	fl=0                                         ;FLAG FOR START
	setLightSpecular "white
	setMaterialSpecular "white
	setMaterialShininess 10
	pd
	SurfaceStart
	repeat (round (u_end-u_start)/u_incr)+1 
	[	u=u+u_incr
		setPC hsb 360*u/(u_end-u_start) 1 1
		bradius=run br_l                     ;ASSIGN BRADIUS FROM FORMULA
		local [tcoord]
		tcoord=offset+(list x y z) ;FIND NEXT CURVE POINT
		ifelse bradius==0 
		[
			setposxyz iposn 
			pd 
			setposxyz tcoord 
			pu
		]
		[
			;	CURVE ONLY PLOT
			if (and (abs u)>8*u_incr
				(sqrt sumsq iposn-tcoord)>1.5*bradius)
			[
				(throw "ERROR [BRADIUS TOO SMALL, OR U INCR TOO LARGE])
			]
			edge2 bradius                     ;GENERATES EDGE VALUES
			SurfaceColumn
;			if fl==1 [band]
;			edge1_l=edge2_l 
;			edge2_l=[]
;			fl=1
;			iposn0=iposn
		]
		iposn=tcoord
		if keyP [throw "stopping]
	]
	SurfaceEnd
	castShadows
	(pr [Drawn in] timefine-timestart [seconds])
end

to assign_val val_list                                ;SETS UP PARAMETERS
	if empty? val_list [stop]
	make  (first val_list)  run (list first bf val_list)
	assign_val bf bf val_list
end

to edge2 r
	;MAKES LIST OF BAND EDGE COORDS
	local [norm_ s theta phi ang vtx vtz]
	norm_=iposn-tcoord                 ;SPHERICAL COORD OF NORMAL TO BAND
	s=Norm norm_
	theta=(atan first norm_ first bf norm_)
	phi=arccos (last norm_)/s
	ang=0

	repeat (round (360/b_incr))+1 
	[	vtx=tx r ang
		vtz=tz r ang 
		setPosXYZ (iposn+(list xc ang yc ang zc ang))
		ang=ang+b_incr
	]
end

to edge1 r
	;AS edge2 BUT START PI OUT OF PHASE, CORRECTS ANOMALY
	local [norm_ s theta phi ang vtx vtz]
	norm_=iposn0-iposn
	s=Norm norm_
	theta=(atan first norm_ first bf norm_)
	phi=arccos (last norm_)/s      ;CHANGE THE OVERLAP BETWEEN BANDS HERE
	ang=180
	edge1_l=[]

	repeat (round (360/b_incr))+1
	[	vtx=tx r ang
		vtz=tz r ang 
		edge1_l=fput (iposn0+(list xc ang yc ang zc ang)) edge1_l
		ang=ang+b_incr
	]
	edge1_l=fput last edge1_l edge1_l
end

to band
	local [p1 p2 q1 q2 dist]
	p1=first edge1_l 
	q1=first edge2_l
	dist=(sqrt sumsq p1-q1)
	
	if dist>2*bradius 
	[                             ;ALLOWS FOR OUT OF PHASE CONDITION
		edge1 bradius
		p1=first edge1_l 
		q1=first edge2_l
	]
	hband bf edge1_l bf edge2_l
end

to hband l1 l2
	if empty? l1 [stop]
	p2=first l1
	q2=first l2
	rect
	p1=p2 q1=q2
	hband bf l1 bf l2
end

to sumsq vec
	op 0+vec*vec
;	if empty? vec [op 0 stop]
;	local [1stel]
;	1stel=first vec
;	op 1stel*1stel+sumsq bf vec
end

to rect
;	pu
;	setposxyz p1
	pd
	polystart
	_setposxyz p2
	_setposxyz q2
	_setposxyz q1
	_setposxyz p1
	polyend
	pu
end

to tx r ang	;X FOR BAND AFTER Z ROTN
	op r*(cos ang)*cos phi
end

to tz r ang	;X FOR BAND AFTER Z ROTN
	op r*(sin ang)
end

to xc ang		;X BAND AFTER Y ROTN
	op vtx*(cos theta)-vtz*sin theta
end

to yc ang
	op vtx*(sin theta)+vtz*cos theta
end

to zc ang
	op  -r*(cos ang)*sin phi
end

to swap
	p1=p0
	p3=p2
end

to x                                                    ;X FOR CURVE
	op run xexpr_l
end

to y
	op run yexpr_l
end

to z
	op run zexpr_l
end

to atan x_ y_
	if (and x_==0 y_==0) [op   0 stop]
	if (and x_==0 y_>0 ) [op  90 stop]
	if (and x_==0 y_<0)  [op -90 stop]
	op (arctan x_ y_)
end