########################################################################## ##PROCEDURE(help,notest) plots/surfintplot ##AUTHOR Robert Israel ##DATE April 2006 ## ##CALLINGSEQ ##- surfintplot(surface1, surface2, options) ##- surfintplot(expr1, expr2, x=a..b, y=c..d, options) ##- surfintplot(proc1, proc2, a..b, c..d, options) ##- surfintplot(expr1, expr2, x=a..b, y=c..d, z=e..f, options) ##- surfintplot(proc1, proc2, a..b, c..d, e..f, options) ## ##PARAMETERS ##- 'surface1','surface2': surface specifications ##- 'expr1','expr2': expressions representing implicit surfaces ##- 'proc1','proc2': procedures representing implicit surfaces ##- 'x','y','z': plotting variables (axes) ##- 'a','b','c','d','e','f': real values for ranges ##- 'options': sequence of options to apply to plot ## ##DESCRIPTION ##- The ~surfintplot~ command plots the 3 dimensional intersections between ## two 3 dimensional surfaces. ## The command is quite flexible, as the surfaces can be specified as ## functions of ~x,y~, as parametric surfaces, or as implicit surfaces. ## ##- Valid surface specifications include the following: ## ##- ~f(x,y), x=a..b, y=c..d~, where ~f(x,y)~ is an expression depending ## only on ~x~,~y~, and ~a~,~b~,~c~,~d~ are real constants. ## This represents the surface ~z=f(x,y)~ for plotting coordinates ## ~(x,y,z)~. ## ##- ~[x(s,t),y(s,t),z(s,t)],s=a..b,t=c..d~, where ~x(s,t),y(s,t),z(s,t)~ are ## expressions depending only on ~s~,~t~, and ~a~,~b~,~c~,~d~ are real ## constants. ## This represents a parametric surface where the three coordinates are ## expressed in terms of the parameters ~s~ and ~t~. ## A vector or Vector may be used in place of the list. ## ##- ~F(x,y,z)=G(x,y,z),x=a..b,y=c..d,z=e..f~, where ~F(x,y,z)~ and ~G(x,y,z)~ ## are expressions depending only on ~x~,~y~,~z~ and ~a~,~b~,~c~,~d~,~e~,~f~ ## are real constants. ## This represents a surface given implicitly by the equation ## ~F(x,y,z)=G(x,y,z)~. ## The equation F(x,y,z) = 0 may be replaced by the expression F(x,y,z). ## ##- In addition, each of these surface specifications can be replaced by ## a procedure form, where the expressions are replaced by procedures, ## and the provided ranges are unnamed: ## ## ~f,a..b,c..d~ ## ## ~[x,y,z],a..b,c..d~ ## ## ~F=G,a..b,c..d,e..f~ or ~F,a..b,c..d,e..f~ ## ##- If both surfaces are parametric or of the form ~z=f(x,y)~ then the ranges ## need only be specified once, and the 2nd and 3rd calling sequences above can be ## used. ## ##- Similarly, if both surfaces are implicit, then the ranges need only ## be specified once, and the 4rth and 5th calling sequences above can be ## used. ## ##- If one surface is parametric and the other is implicit, ## ~-infinity..infinity~ may be used in the ranges for the implicit surface. ## Otherwise all ranges must be positive and nonempty, with finite constant ## bounds. ## ##EXAMPLES(noexecute) ## ##> with(plots,surfintplot); ## ##-(lead=false) Two implicit surfaces ##> surfintplot(x^2+y^2+z^2=1, x^2-y^2=z*(z-1), x=-1..1, y=-1..1, ## z=-1..1, axes=box, thickness=2, orientation=[70,40]); ## ##-(lead=false) Plotting along with the sphere ##> plots[display](%,plottools[sphere]([0,0,0],1)); ## ##-(lead=false) Same surface intersections using the parametric form for the sphere ##> surfintplot( ## [sin(s)*cos(t),sin(s)*sin(t),cos(s)], s=0..Pi, t=0..2*Pi, ## x^2-y^2=z*(z-1), x=-1..1, y=-1..1, z=-1..1, ## axes=box, thickness=2, orientation=[70,40]); ## ##-(lead=false) Same surface intersections with both parametric ##> surfintplot( ## [sin(s)*cos(t),sin(s)*sin(t),cos(s)], s=0..Pi, t=0..2*Pi, ## [sinh(v)/2,cosh(v)*cos(u)/2,(cosh(v)*sin(u)+1)/2], u=0..2*Pi, v=-3..3, ## axes=box, thickness=2, orientation=[70,40]); ## ##-(lead=false) A difficult example, as the surfaces are nearly tangent along a curve ##> P1 := surfintplot([cos(u+v),sin(u-v),cos(u)+cos(v)],u=0..2*Pi,v=0..2*Pi, ## [sin(s)*cos(t),sin(s)*sin(t),cos(s)],s=0..Pi,t=-Pi..Pi, ## orientation=[145,90]): ##> P1; ##> P2 := plot3d([cos(u+v),sin(u-v),cos(u)+cos(v)],u=0..2*Pi,v=0..2*Pi): ##> P3 := plottools[sphere]([0,0,0],1,colour=blue): ##> plots[display]([P1,P2,P3],orientation=[145,90]); ## # # This code was written by Robert Israel under SFU's MITACS project. # Dated: April 2006. # surfintplot := module() local ModuleApply, parpar, parimp, gridtris, qtri, qtri2, tridiv, triinterp, triint, qpair, pair, qtet1, qtet2, tetinterp, tetradiv, impimp, cubetets, gridtets; ModuleApply := proc() local surfs, bds, rangeargs, i,j, astart, functype, axname, fproc, f12proc, xyzproc, surftype, badargs, lobd, hibd, parnames, R, badarg, z; R := 1: R := 'R'; # Mint warnings # Classification uses arguments that are ranges or name = range rangeargs := map(type,[args],{`..`, 'And'('name','Not'('identical'('view')))=`..`}); # cases allowed: # ((false, true, true) or (false,true,true,true)) twice # or false, false, true, true, true : Handled by the code below # or false, false, true, true : Transformed for handling by later code if nops(rangeargs)>=5 and rangeargs[1..5]=[false,false,true,true,true] then # case of two implicit surfaces with same ranges bds := [args[3..5]]; if type(bds,'list'(`=`)) then axname := map(lhs,bds); bds := map(rhs,bds); if not type(args[1],{'algebraic',`=`}) or not type(args[2],{'algebraic',`=`}) or indets([args[1..2]],'name') minus {constants,op(axname)}<>{} then error "input surface expressions must be algebraic, and only depend upon the plotting variables"; end if; functype := 'expression'; elif type(bds, 'list'(`..`)) then axname := ['none','none','none']; if not type(args[1],'procedure') or not type(args[2],'procedure') then error "input using unnamed ranges must be specified using procedures"; end if; functype := 'procedure'; else error "invalid input: coordinates must be all named or all unnamed" end if; for i to 2 do surfs[i] := `if`(type(args[i],`=`),(lhs-rhs)(args[i]),args[i]); if functype = 'expression' then fproc[i] := unapply(subs(seq(axname[j]=R[j],j=1..3),'evalf'(surfs[i])), R); else fproc[i] := surfs[i]; end if; end do; f12proc := subs(A=fproc[1],B=fproc[2], proc(R) option remember; [A(R), B(R)] end ); return impimp(f12proc, op(bds), axname, [args[6..-1]]); elif nops(rangeargs)>=4 and rangeargs[1..4]=[false,false,true,true] then if type(args[1],{list,array,Array}) or type(args[2],{list,array,Array}) then # case of two parametric surfaces with same ranges: # call recursively with range args repeated return procname(args[1],args[3..4],args[2],args[3..-1]) else # case of two z=f(x,y) surfaces with same ranges: # do one as implicit, other parametric bds := [args[3..4]]; if type(bds,'list'(`=`)) then axname := map(lhs,bds); if not type(args[1],'algebraic') or not type(args[2],'algebraic') or indets([args[1..2]],'name') minus {constants,op(axname)}<>{} then error "input surface expressions must be algebraic or parametric, and only depend upon the plotting variables"; end if; return procname(args[1]-z, args[3],args[4],z=-infinity..infinity, [lhs(args[3]),lhs(args[4]),args[2]], args[3..-1]); elif type(bds,'list'(`..`)) then if not type(args[1],'procedure') or not type(args[2],'procedure') then error "input using unnamed ranges must be specified using procedures"; end if; return procname(subs(Q=args[1],(x,y,z)->Q(x,y)-z), args[3], args[4], -infinity..infinity, [(x,y)->x, (x,y)->y, args[2]], args[3..-1]); else error "invalid input: coordinates must be all named or all unnamed"; end if; end if; end if; # We should now have only # rangeargs = ((false, true, true) or (false,true,true,true)) twice astart := 1; for i to 2 do if nops(rangeargs)>=astart+3 and rangeargs[astart..astart+3]=[false,true,true,true] then # this one is implicit surftype[i] := 'implicit'; surfs[i] := `if`(type(args[astart],`=`),(lhs-rhs)(args[astart]),args[astart]); bds[i] := args[astart+1..astart+3]; astart := astart+4; elif nops(rangeargs)>=astart+2 and rangeargs[astart..astart+2]=[false,true,true] then # this one is parametric surftype[i] := 'parametric'; surfs[i] := args[astart]; bds[i] := args[astart+1..astart+2]; astart := astart + 3; else # can't figure it out badargs := rangeargs[astart..astart+2]; badargs[1] := not badargs[1]; member(false,badargs,'badarg'); error "invalid input: argument %1 is of the wrong type", args[badarg+astart-1]; end if; # check all ranges named or all unnamed # if named, surface is an expression, otherwise a procedure if type([bds[i]], 'list'(`=`)) then axname := map(lhs,[bds[i]]); if not type(surfs[i],{'algebraic','list'('algebraic')}) or indets(surfs[i],'name') minus {constants,op(axname)}<>{} then error "input surface expressions must be algebraic, and only depend upon the plotting variables"; end if; functype[i] := 'expression'; elif type([bds[i]], 'list'(`..`)) then functype[i] := 'procedure' else error "invalid input: parameters for each surface must be all named or all unnamed"; end if; end do; # At this point, we have two surfaces, 1, 2, and we now handle the cases # where these have different forms (i.e. implicit or parametric) if surftype[1] = 'implicit' and surftype[2] = 'implicit' then # two implicit surfaces with separately specified ranges: intersect them axname := [``,``,``]; for i to 3 do lobd[i] := -infinity; hibd[i] := infinity; for j to 2 do if functype[j] = 'expression' then axname[i] := op(1,bds[j][i]); lobd[i] := max(lobd[i], op([2,1],bds[j][i])); hibd[i] := min(hibd[i], op([2,2],bds[j][i])); else lobd[i] := max(lobd[i], op(1,bds[j][i])); hibd[i] := min(hibd[i], op(2,bds[j][i])); end if; end do; if lobd[i]>=hibd[i] then error "ranges must be non-empty"; end if; end do; # the axis names, if any, must be the same if functype[1] = 'expression' and functype[2] = 'expression' and map(lhs,[bds[1]]) <> map(lhs,[bds[2]]) then error "conflicting names for axes: %1 and %2", map(lhs,[bds[1]]), map(lhs,[bds[2]]); end if; # get the procedure for [f1,f2] for i to 2 do if functype[i] = 'expression' then fproc[i] := unapply(subs(seq(axname[j]=R[j],j=1..3),'evalf'(surfs[i])), R); else fproc[i] := surfs[i]; end if; end do; f12proc := subs(A=fproc[1],B=fproc[2], proc(R) option remember; [A(R), B(R)]; end proc); impimp(f12proc, seq(lobd[i]..hibd[i],i=1..3), axname, [args[astart..-1]]); elif surftype[1] = 'implicit' and surftype[2] = 'parametric' then # implicit and parametric # get function and axes for the implicit surface if functype[1] = 'expression' then axname := map(lhs,[bds[1]]); bds[1] := op(map(rhs,[bds[1]])); fproc := unapply(subs(seq(axname[j]=R[j],j=1..3),'evalf'(surfs[1])), R); else axname := [``,``,``]; fproc := surfs[1]; end if; # parametric surface if functype[2] = 'expression' then # expression case parnames := map(lhs,[bds[2]]); bds[2] := op(map(rhs,[bds[2]])); # check if parametric is z as function of x,y if type(surfs[2],{'list','array','Array'}) then xyzproc := [seq(unapply('evalf'(surfs[2][i]),parnames[1], parnames[2]),i=1..3)]; else xyzproc := [unapply('x','x','y'), unapply('y','x','y'), unapply('evalf'(surfs[2]),parnames[1],parnames[2])]; end if; else # procedure case for parametric surface if type(surfs[2],{'list','array','Array'}) then xyzproc := surfs[2]; if type(xyzproc,{'array','Array'}) then xyzproc := convert(xyzproc,'list'); end if; else xyzproc := [unapply('x','x','y'), unapply('y','x','y'), surfs[2]]; end if; end if; parimp(xyzproc, bds[2], fproc, bds[1], axname, [args[astart..-1]]) elif surftype[1] = 'parametric' and surftype[2] = 'implicit' then # parametric and implicit # get function and axes for the implicit surface if functype[2] = 'expression' then axname := map(lhs,[bds[2]]); bds[2] := op(map(rhs,[bds[2]])); fproc := unapply(subs(seq(axname[j]=R[j],j=1..3),'evalf'(surfs[2])), R); else axname := ['none','none','none']; fproc := surfs[2]; end if; # functype[2] = 'expression' # parametric surface if functype[1] = 'expression' then # expression case parnames := map(lhs,[bds[1]]); bds[1] := op(map(rhs,[bds[1]])); # check if parametric is z as function of x,y if type(surfs[1],{'list','array','Array'}) then xyzproc := [seq(unapply('evalf'(surfs[1][i]),parnames[1], parnames[2]),i=1..3)]; else xyzproc := [unapply('x','x','y'), unapply('y','x','y'), unapply('evalf'(surfs[1]),parnames[1],parnames[2])]; end if; else # procedure case for parametric surface if type(surfs[1],{'list','array','Array'}) then xyzproc := surfs[1]; if type(xyzproc,{'array','Array'}) then xyzproc := convert(xyzproc,'list'); end if; else xyzproc := [unapply('x','x','y'), unapply('y','x','y'), surfs[1]]; end if; end if; parimp(xyzproc, bds[1], fproc, bds[2], axname, [args[astart..-1]]); elif surftype[1] = 'parametric' and surftype[2] = 'parametric' then # parametric and parametric for i to 2 do if functype[i] = 'expression' then # expression case parnames[i] := map(lhs,[bds[i]]); bds[i] := op(map(rhs,[bds[i]])); # check if parametric is z as function of x,y if type(surfs[i],{'list','array','Array'}) then xyzproc[i] := subsop(3=remember, unapply(([seq('evalf'( surfs[i][j]),j=1..3), parnames[i][1],parnames[i][2]]), parnames[i][1], parnames[i][2])); else xyzproc[i] := subsop(3=remember,unapply(([parnames[i][1], parnames[i][2],'evalf'(surfs[i]), parnames[i][1], parnames[i][2]]), parnames[i][1], parnames[i][2])); end if; else # procedure case for parametric surface if type(surfs[i],{'list','array','Array'}) then xyzproc[i] := subs(Q=surfs[i], proc(x,y) option remember; ([Q[1](x,y),Q[2](x,y),Q[3](x,y),x,y]) end); else xyzproc[i] := subsop(3=remember,unapply((['evalf'(surfs[i][1]), 'evalf'(surfs[i][2]),'evalf'(surfs[i][3]),parnames[i][1], parnames[i][2]]),parnames[i][1], parnames[i][2])); end if; end if; end do; parpar(xyzproc[1], bds[1], xyzproc[2], bds[2], [args[astart..-1]]); else error "should not happen"; end if; end proc; triint := proc(a,b, c,d, e,f) local ndef, nab, nac, ab, ac, de, df, dede, dfdf, dedf, u, v, t, h, g, lo, hi, ude, vdf, t2, t1, hd, p1, p2, q1, q2; ab := b-a; ac := c-a; de := e-d; df := f-d; # vector differences ndef := [de[2]*df[3]-de[3]*df[2],de[3]*df[1]-de[1]*df[3], de[1]*df[2]-de[2]*df[1]]; # normal to plane through d,e,f if abs(ndef[1])+abs(ndef[2])+abs(ndef[3]) <= 10^(-ceil(Digits/2))* (abs(de[1])+abs(de[2])+abs(de[3]))*(abs(df[1])+abs(df[2])+abs(df[3])) then return; # *** treat this special case as a line? *** end if; nab := ndef[1]*ab[1] + ndef[2]*ab[2] + ndef[3]*ab[3]; nac := ndef[1]*ac[1] + ndef[2]*ac[2] + ndef[3]*ac[3]; # dot products of ndef with b-a and c-a if abs(nab) < abs(nac) then # interchange b,c if needed so |nab| >= |nac| t := nab; nab := nac; nac := t; t := ab; ab := ac; ac := t; end if; if abs(nab) <= 10^(-ceil(Digits/2))*(abs(ab[1])+abs(ab[2])+abs(ab[3]))* (abs(ndef[1])+abs(ndef[2])+abs(ndef[3])) then WARNING("planes nearly parallel"); return; # *** maybe return a polygon *** end if; dede := de[1]^2 + de[2]^2 + de[3]^2; dfdf := df[1]^2 + df[2]^2 + df[3]^2; dedf := de[1]*df[1] + de[2]*df[2] + de[3]*df[3]; u := de - dedf/dfdf * df; ude := dede - (dedf)^2/dfdf; v := df - dedf/dede * de; vdf := dfdf - (dedf)^2/dede; # vectors u,v in plane of de, df so u.df = 0, v.de = 0 # ude = u.de, vdf = v.df t2 := nac/nab; g := ac - t2*ab; t1 := (ndef[1]*(d[1]-a[1]) + ndef[2]*(d[2]-a[2]) + ndef[3]*(d[3]-a[3]))/nab; h := a + t1 *ab; hd := h - d; # now x = h + s*g is intersection of planes. # where x = a + (t1 - s t2) ab + s ac # Get upper and lower bounds on s lo := 0; hi := 1; # need t1 - s * t2 >= 0 if t2 > 0 then hi := min(hi, t1/t2); elif t2 < 0 then lo := max(lo, t1/t2); elif t1 < 0 then return; # if t2 = 0, need t1 >= 0 end if; # need t1 + s * (1 - t2) <= 1 if t2 < 1 then hi := min(hi, (1-t1)/(1-t2)); elif t2 > 1 then lo := max(lo, (1-t1)/(1-t2)); elif t1 > 1 then return; # if t2 = 1, need t1 <= 1 end if; p1 := (u[1]*hd[1] + u[2]*hd[2] + u[3]*hd[3])/ude; p2 := (u[1]*g[1] + u[2]*g[2] + u[3]*g[3])/ude; q1 := (v[1]*hd[1] + v[2]*hd[2] + v[3]*hd[3])/vdf; q2 := (v[1]*g[1] + v[2]*g[2] + v[3]*g[3])/vdf; # x = d + (p1 + s*p2)*de + (q1 + s*q2)*df # need p1 + s * p2 >= 0 if p2 > 0 then lo := max(lo,-p1/p2); elif p2 < 0 then hi := min(hi,-p1/p2); elif p1 < 0 then return; # if p2 = 0, need p1 >= 0 end if; # need q1 + s * q2 >= 0 if q2 > 0 then lo := max(lo,-q1/q2); elif q2 < 0 then hi := min(hi,-q1/q2); elif q1 < 0 then return; end if; # need (p1 + q1) + s*(p2 + q2) <= 1 if p2 + q2 > 0 then hi := min(hi, (1-p1-q1)/(p2+q2)); elif p2 + q2 < 0 then lo := max(lo, (1-p1-q1)/(p2+q2)); elif p1 + q1 > 1 then return; end if; if lo < hi then # OK, return the line segment [h + lo*g, h + hi*g] end if; end proc; parpar := proc( xyz1, S1::range(realcons), T1::range(realcons), xyz2, S2::range(realcons), T2::range(realcons), opts::list ) local Defaults, Options, others, a, b, c, d, f, S, M, N, A, NBIN, NSUB, i, ds, dt, j, k, cmin, cmax, dbin, L, tris, ntris, tri, tmin, tmax, bmin, bmax, b1, b2, b3, bin, Mark, r, nstris, stris, B, P1, P2, P3, Res, i1, i2, pp, v, features; Defaults := table(['OptionParms' = table([ 'grid1' = [[10,10], {posint, [posint,posint]}, 'makelist'], 'grid2' = [[10,10], {posint, [posint,posint]}, 'makelist'], 'bins' = [5,posint], 'subdivisions' = [5,posint], 'coords' = ['cartesian',{name, function}]]), 'PositionalParms' = table()]); others := [ProcessOptions(0, opts, Defaults, Options)]; if hasoption(others, {'color','colour'}, 'v', 'others') then Options['colour'] := v; else Options['colour'] := 'red'; end if; features := `plot3d/options3d`(op(others), 'colour'=Options['colour']); M := Options['grid1']; if nops(M) = 1 then M := [M[1],M[1]]; end if; N := Options['grid2']; if nops(N) = 1 then N := [N[1],N[1]]; end if; NBIN := Options['bins']; NSUB := Options['subdivisions']; a[1] := evalf(lhs(S1)); b[1] := evalf(rhs(S1)); c[1] := evalf(lhs(T1)); d[1] := evalf(rhs(T1)); a[2] := evalf(lhs(S2)); b[2] := evalf(rhs(S2)); c[2] := evalf(lhs(T2)); d[2] := evalf(rhs(T2)); A[1] := Array(1..M[1]+1,1..N[1]+1); A[2] := Array(1..M[2]+1,1..N[2]+1); # Calculate grid points. Note that the vertices must have 5 entries: # their x,y,z coordinates and also their parameter s,t values (for subdividing) # need to fully evaluate the functions xyz1 and xyz2. f[1] := eval(xyz1); f[2] := eval(xyz2); for i to 2 do ds[i] := (b[i]-a[i])/M[i]; dt[i] := (d[i]-c[i])/N[i]; for j from 0 to M[i] do for k from 0 to N[i] do A[i][j+1,k+1] := evalf(f[i](a[i]+j*ds[i],c[i]+k*dt[i])); end do; end do; end do; # get mins and maxes, and set up bins for k to 3 do S := seq(seq(seq(A[L][i,j][k],i=1..M[L]+1),j=1..N[L]+1),L=1..2); cmin[k] := min(S); cmax[k] := max(S); dbin[k] := (cmax[k]-cmin[k])/NBIN; if dbin[k] = 0 then dbin[k] := 0.1; end if; end do; # get triangles for L to 2 do tris[L] := [ seq(seq([A[L][i,j],A[L][i+1,j],A[L][i,j+1]],i=1..M[L]),j=1..N[L]), seq(seq([A[L][i+1,j],A[L][i,j+1],A[L][i+1,j+1]],i=1..M[L]),j=1..N[L]) ]; ntris[L] := nops(tris[L]); end do; # enter triangles in bins for L to 2 do for i to ntris[L] do tri := tris[L][i]; for d to 3 do tmin := min(tri[1][d],tri[2][d],tri[3][d]); tmax := max(tri[1][d],tri[2][d],tri[3][d]); bmin[d] := ceil((tmin-cmin[d])/dbin[d]); bmax[d] := min(NBIN,ceil((tmax-cmin[d])/dbin[d])); end do; for b1 from bmin[1] to bmax[1] do for b2 from bmin[2] to bmax[2] do for b3 from bmin[3] to bmax[3] do if not assigned(bin[L][b1,b2,b3]) then bin[L][b1,b2,b3] := i; else bin[L][b1,b2,b3] := i,bin[L][b1,b2,b3]; end if; end do; end do; end do; end do; end do; # mark triangles for subdivision if they are in a bin that has triangles from # both surfaces Mark[1] := Array[1..ntris[1]]; Mark[2] := Array[1..ntris[2]]; for r in {indices(bin[1])} intersect {indices(bin[2])} do for L to 2 do for i in [bin[L][op(r)]] do Mark[L][i] := 1; end do; end do; end do; # get second level triangles in stris[1] and stris[2] cmin := [infinity,infinity,infinity]; cmax := [-infinity,-infinity,-infinity]; for L to 2 do nstris[L] := 0; stris[L] := TABLE(); for i to ntris[L] do if Mark[L][i] = 1 then B := Array(0..NSUB,0..NSUB); P1 := tris[L][i][1][4..5]; P2 := tris[L][i][2][4..5]; P3 := tris[L][i][3][4..5]; for j from 0 to NSUB do for k from 0 to NSUB - j do pp := P1+j*(P2-P1)/NSUB+k*(P3-P1)/NSUB; B[j,k] := evalf(f[L](pp[1],pp[2])[1..3]); for d to 3 do cmin[d] := min(cmin[d],B[j,k][d]); cmax[d] := max(cmax[d],B[j,k][d]); end do; end do; end do; for j from 0 to NSUB-1 do for k from 0 to NSUB-1-j do stris[L][nstris[L]+1] := [B[j,k],B[j+1,k],B[j,k+1]]; if j+k+2<=NSUB then nstris[L] := nstris[L]+2; stris[L][nstris[L]] := [B[j+1,k],B[j+1,k+1],B[j,k+1]]; else nstris[L] := nstris[L]+1; end if; end do; end do; end if; end do; end do; # now put in second level bins bin := 'bin'; NBIN := NBIN*5; for d to 3 do dbin[d] := (cmax[d]-cmin[d])/NBIN; end do; for L to 2 do for i to nstris[L] do tri := stris[L][i]; for d to 3 do tmin := min(tri[1][d],tri[2][d],tri[3][d]); tmax := max(tri[1][d],tri[2][d],tri[3][d]); bmin[d] := ceil((tmin-cmin[d])/dbin[d]); bmax[d] := ceil((tmax-cmin[d])/dbin[d]); end do; for b1 from bmin[1] to bmax[1] do for b2 from bmin[2] to bmax[2] do for b3 from bmin[3] to bmax[3] do if not assigned(bin[L][b1,b2,b3]) then bin[L][b1,b2,b3] := i; else bin[L][b1,b2,b3] := i,bin[L][b1,b2,b3]; end if; end do; end do; end do; end do; end do; # search bins for the second-level pairs that might intersect # storing results in Res # remember those pairs since they might arise in another bin pair := 'pair'; Res := PLOT3D(CURVES(seq(seq(seq(qpair(i1,i2,stris[1][i1],stris[2][i2]), i2 = bin[2][op(r)]),i1 = bin[1][op(r)]), r = {indices(bin[1])} intersect {indices(bin[2])})), features); if Options['coords'] <> 'cartesian' then Res := plots[changecoords](Res, Options['coords']); end if; Res; end proc: qpair := proc(i1,i2,st1,st2) if not assigned(pair[i1,i2]) then pair[i1,i2] := 1; triint(op(st1),op(st2)); end if; end proc; # intersection of parametric and implicit parimp := proc( xyzproc::list, S1::range, S2::range, Fproc, xbd::range, ybd::range, zbd::range, axname::list, opts::list ) local a,b,c,d, Defaults, Options, features, others, M, MAXLEV, ds, dt, fproc, A, agenda, lobds, hibds, inbds, Res, i, j, tri, v; a := lhs(S1); b := rhs(S1); c := lhs(S2); d := rhs(S2); Defaults := table(['OptionParms' = table([ 'grid' = [[25,25], {posint, [posint,posint]}, 'makelist'], 'maxlev' = [5,posint], 'coords' = ['cartesian',{name, function}], 'view' = [[xbd, ybd, zbd], [range, range, range]]]), 'PositionalParms' = table()]); others := [ProcessOptions(0, opts, Defaults, Options)]; if hasoption(others,{'color','colour'}, 'v', 'others') then Options['colour'] := v; else Options['colour'] := 'red'; end if; if has(Options['view'],-infinity..infinity) then Options['view']:= subs(-infinity..infinity = DEFAULT, Options['view']) fi; features := `plot3d/options3d`(op(others), 'colour'=Options['colour'], 'view'=Options['view']); M := Options['grid']; if nops(M) = 1 then M := [M[1],M[1]]; end if; MAXLEV := Options['maxlev']; _Envsignum0 := 0; ds := (b-a)/M[1]; dt := (d-c)/M[2]; fproc := proc(s,t) option remember; Fproc(xyzproc(s,t)); end; A := Array(0..M[1],0..M[2],(i,j) -> fproc(a+i*ds,c+j*dt)); agenda := [seq](seq(gridtris(i,j,A,ds,dt,a,c),i=0..M[1]-1),j=0..M[2]-1); # filter triangles in agenda for bounds lobds := map(lhs,[xbd, ybd, zbd]); hibds := map(rhs,[xbd, ybd, zbd]); inbds := proc(tri) local P, i, j; P := xyzproc(tri[1][1],tri[1][2]); for i to 3 do if P[i] < lobds[i] then for j from 2 to 3 do if xyzproc[i](tri[j][1],tri[j][2]) < lobds[i] then return false; end if; end do; elif P[i] > hibds[i] then for j from 2 to 3 do if xyzproc[i](tri[j][1],tri[j][2]) > hibds[i] then return false; end if; end do; end if; end do; true; end proc; agenda := select(inbds, agenda); for i to MAXLEV while nops(agenda) < 50000 do agenda := [seq](tridiv(op(tri),fproc), tri = agenda); end do; Res := PLOT3D(CURVES(seq(triinterp(tri,xyzproc),tri=agenda)),features); if Options['coords'] <> 'cartesian' then Res := plots[changecoords](Res, Options['coords']); end if; Res; end proc; gridtris := proc(i,j,A,ds,dt,a,c) # A is array of values for f(x(s,t),y(s,t),z(s,t)) # for points [s,t] in grid # return candidate grid triangles [s,t,A value] in rectangle # [a+i*ds .. a+(i+1)*ds, c+j*dt .. c+(j+1)*dt] local a1,a2,c1,c2,sigs; a1 := a+i*ds; a2 := a1+ds; c1 := c+j*dt; c2 := c1 + dt; sigs := [signum(A[i,j]), signum(A[i+1,j]), signum(A[i,j+1]), signum(A[i+1,j+1])]; qtri([a1,c1,A[i,j]],sigs[1],[a2,c1,A[i+1,j]],sigs[2], [a2,c2,A[i+1,j+1]],sigs[4]), qtri([a1,c1,A[i,j]],sigs[1],[a1,c2,A[i,j+1]],sigs[3], [a2,c2,A[i+1,j+1]],sigs[4]); end proc; qtri := proc(v1,s1,v2,s2,v3,s3) local sigs; sigs := {s1,s2,s3}; if sigs = {1} or sigs = {-1} then NULL; else [v1,v2,v3]; end if; end; qtri2 := proc(v1,s1,v2,s2,v3,s3) local sigs; sigs := {s1,s2,s3}; if nops(sigs) = 1 then NULL; else [v1,v2,v3]; end if; end; tridiv := proc(A,B,C, f) # A,B,C consist of points with f values for triangle vertices. # subdivide and return up to 4 sub-triangles local a,b,c,sa,sb,sc,sab,sbc,sac,ab,ac,bc,vab,vac,vbc; a := A[1..2]; b := B[1..2]; c := C[1..2]; sa := signum(A[3]); sb := signum(B[3]); sc := signum(C[3]); ab := (a+b)/2; ac := (a+c)/2; bc := (b+c)/2; vab := [ab[1],ab[2],f(ab[1],ab[2])]; sab := signum(vab[3]); vac := [ac[1],ac[2],f(ac[1],ac[2])]; sac := signum(vac[3]); vbc := [bc[1],bc[2],f(bc[1],bc[2])]; sbc := signum(vbc[3]); qtri2(A,sa,vab,sab,vac,sac), qtri2(vac,sac,vab,sab,vbc,sbc), qtri2(C,sc,vac,sac,vbc,sbc), qtri2(B,sb,vab,sab,vbc,sbc); end; triinterp := proc(T,xyzproc) # T a triangle consisting of [s,t,f] points # interpolate curve and map to [x,y,z] coordinates. local T1,T2,T3,s1,s2,s3,r1,p11,p12,r2,p21,p22; T1 := T[1]; s1 := signum(T1[3]); T2 := T[2]; s2 := signum(T2[3]); T3 := T[3]; s3 := signum(T3[3]); if s1 = s2 then # curve from T1,T3 line to T2,T3 line if s1 = s3 then if s1 = 0 then return [xyzproc(T1[1],T1[2]), xyzproc(T2[1],T2[2]), xyzproc(T3[1],T3[2]), xyzproc(T1[1],T1[2])]; else # shouldn't happen that all + or all - error "all same sign in %a", T end if; end if; r1 := T1[3]/(T1[3] - T3[3]); p11 := T1[1]+r1*(T3[1]-T1[1]); p12 := T1[2]+r1*(T3[2]-T1[2]); r2 := T2[3]/(T2[3] - T3[3]); p21 := T2[1]+r2*(T3[1]-T2[1]); p22 := T2[2]+r2*(T3[2]-T2[2]); else # s1 <> s2; one of these is nonzero # T1,T2 line is there r1 := T1[3]/(T1[3] - T2[3]); p11 := T1[1]+r1*(T2[1]-T1[1]); p12 := T1[2]+r1*(T2[2]-T1[2]); if (s1 = s3) or (s1 = 0 and s2 <> s3) then # go to T2,T3 line r2 := T2[3]/(T2[3] - T3[3]); p21 := T2[1]+r2*(T3[1]-T2[1]); p22 := T2[2]+r2*(T3[2]-T2[2]); else # must be T1,T3 line r2 := T1[3]/(T1[3] - T3[3]); p21 := T1[1]+r2*(T3[1]-T1[1]); p22 := T1[2]+r2*(T3[2]-T1[2]); end if; end if; [xyzproc(p11,p12), xyzproc(p21,p22)]; end proc; qtet1 := proc(v1,v2,v3,v4,V1,V2,V3,V4) # check that values of neither function are all positive or all negative # and if so return [point,values] pairs for tetrahedron local sig1,sig2; sig1 := map(t -> signum(t[1]),{V1,V2,V3,V4}); sig2 := map(t -> signum(t[2]),{V1,V2,V3,V4}); if sig1 <> {-1} and sig1 <> {1} and sig2 <> {-1} and sig2 <> {1} then [[v1,V1],[v2,V2],[v3,V3],[v4,V4]]; else NULL; end if; end proc; qtet2 := proc(v1,v2,v3,v4,V1,V2,V3,V4) # check that values of neither function are all same sign # and if so return [point,values] pairs for tetrahedron # this differs from qtet1 in not accepting signs all 0 local sig1, sig2; sig1 := map(t -> signum(t[1]),{V1,V2,V3,V4}); sig2 := map(t -> signum(t[2]),{V1,V2,V3,V4}); if nops(sig1) > 1 and nops(sig2) > 1 then [[v1,V1],[v2,V2],[v3,V3],[v4,V4]]; else NULL; end if; end proc; tetinterp := proc(T) # a tetrahedron consisting of [vertex, fvalue] pairs local A1, A2, N1, N2, n1n1, n1n2, n2n2, det, N3, a, b, X0, lo, hi, s, plo, phi; A1 := T[1][2][1]; A2 := T[1][2][2]; N1 := [T[2][2][1]-A1, T[3][2][1]-A1, T[4][2][1]-A1]; N2 := [T[2][2][2]-A2, T[3][2][2]-A2, T[4][2][2]-A2]; n1n1 := N1[1]^2 + N1[2]^2 + N1[3]^2; n1n2 := N1[1]*N2[1]+N1[2]*N2[2]+N1[3]*N2[3]; n2n2 := N2[1]^2 + N2[2]^2 + N2[3]^2; det := n1n1*n2n2-n1n2^2; if det = 0 then return; end if; N3 := [N1[2]*N2[3]-N1[3]*N2[2], N1[3]*N2[1]-N1[1]*N2[3], N1[1]*N2[2]-N1[2]*N2[1]]; a := (n1n2*A2 - n2n2*A1)/det; b := (n1n2*A1-n1n1*A2)/det; X0 := a*N1 + b*N2; lo := -infinity; hi := infinity; if N3[1] > 0 then lo := - X0[1]/N3[1]; elif N3[1] < 0 then hi := -X0[1]/N3[1]; elif X0[1] < 0 then return; end if; if N3[2] > 0 then lo := max(lo, -X0[2]/N3[2]); elif N3[2] < 0 then hi := min(hi, -X0[2]/N3[2]); elif X0[2] < 0 then return; end if; if N3[3] > 0 then lo := max(lo, -X0[3]/N3[3]); elif N3[3] < 0 then hi := min(hi, -X0[3]/N3[3]); elif X0[3] < 0 then return; end if; s := N3[1]+N3[2]+N3[3]; if s > 0 then hi := min(hi,(1 - X0[1]-X0[2]-X0[3])/s); elif s < 0 then lo := max(lo, (1 - X0[1]-X0[2]-X0[3])/s); elif X0[1] + X0[2] + X0[3] > 1 then return; end if; if lo >= hi then return; end if; plo := X0 + lo*N3; phi := X0 + hi*N3; plo := T[1][1] + plo[1]*(T[2][1]-T[1][1]) + plo[2]*(T[3][1]-T[1][1]) + plo[3]*(T[4][1]-T[1][1]); phi := T[1][1] + phi[1]*(T[2][1]-T[1][1]) + phi[2]*(T[3][1]-T[1][1]) + phi[3]*(T[4][1]-T[1][1]); [plo,phi]; end proc; tetradiv := proc(A,B,C,D, f) # A,B,C,D pairs consisting of the vertices of a tetrahedron with f values # return sequence of up to 8 tetrahedra, each a list of 4 vertices with f1 and f2 values # which are either vertices of the original or bisectors of its edges. local a,b,c,d,ab,ac,ad,bc,bd,cd,AB,AC,AD,BC,BD,CD; a := A[1]; b := B[1]; c := C[1]; d := D[1]; ab := (a+b)/2; ac := (a+c)/2; ad := (a+d)/2; AB := f(ab); AC := f(ac); AD := f(ad); bc := (b+c)/2; bd := (b+d)/2; cd := (c+d)/2; BC := f(bc); BD := f(bd); CD := f(cd); qtet2(a,ab,ac,ad,A[2],AB,AC,AD), qtet2(b,ab,bc,bd,B[2],AB,BC,BD), qtet2(c,ac,bc,cd,C[2],AC,BC,CD), qtet2(d,ad,bd,cd,D[2],AD,BD,CD), qtet2(ab,bc,ad,bd,AB,BC,AD,BD), qtet2(ad,cd,bc,bd,AD,CD,BC,BD), qtet2(ab,bc,ad,ac,AB,BC,AD,AC), qtet2(ad,cd,bc,ac,AD,CD,BC,AC); end proc; impimp := proc( f12proc, R1::range, R2::range, R3::range, axname::list, opts::list ) # intersection of implicit surfaces # f1 and f2 expressions in coordinates x,y,z, r1, r2 and r3 give ranges # for those coordinates local Defaults, others, Options, features, GRID1, GRID2, GRID3, MAXLEV, Res, a1, a2, a3, b1, b2, b3, agenda, dx, dy, dz, A, i, j, k, v, tet; a1 := lhs(R1); a2 := lhs(R2); a3 := lhs(R3); b1 := rhs(R1); b2 := rhs(R2); b3 := rhs(R3); Defaults := table(['OptionParms' = table([ 'grid' = [[10,10,10], {posint, [posint,posint,posint]}, 'makelist'], 'maxlev' = [3,posint], 'coords' = ['cartesian',{name, function}], 'view' = [[R1, R2, R3], [range, range, range]]]), 'PositionalParms' = table()]); others := [ProcessOptions(0, opts, Defaults, Options)]; if hasoption(others,{'color','colour'}, 'v', 'others') then Options['colour'] := v; else Options['colour'] := 'red'; end if; features := `plot3d/options3d`(op(others), 'colour'=Options['colour'], 'view'=Options['view']); if nops(Options['grid']) = 1 then Options['grid'] := [Options['grid'][1]$3]; end if; GRID1, GRID2, GRID3 := op(Options['grid']); MAXLEV := Options['maxlev']; _Envsignum0 := 0; dx := (b1-a1)/GRID1; dy := (b2-a2)/GRID2; dz := (b3-a3)/GRID3; # tetrahedra list for unit cube A := Array(0..GRID1,0..GRID2,0..GRID3,1..2, (j,k,L,i) -> f12proc([a1+j*dx,a2+k*dy,a3+L*dz])[i]); agenda := [seq](seq(seq(gridtets(i,j,k,A,dx,dy,dz,a1,a2,a3), i=0..GRID1-1),j=0..GRID2-1),k=0..GRID3-1); for i to MAXLEV do # protection against too many segments if nops(agenda) > 50000 then break; end if; agenda := [seq](tetradiv(op(tet),f12proc), tet = agenda); end do; Res := PLOT3D(CURVES(seq(tetinterp(tet), tet = agenda)),features); if Options['coords'] <> 'cartesian' then Res := plots[changecoords](Res, Options['coords']); end if; Res; end proc; # tetrahedra list for unit cube cubetets := [{[0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1]}, {[1, 1, 0], [1, 0, 0], [0, 1, 0], [1, 1, 1]}, {[0, 1, 1], [0, 1, 0], [0, 0, 1], [1, 1, 1]}, {[1, 0, 1], [1, 0, 0], [0, 0, 1], [1, 1, 1]}, {[1, 0, 0], [0, 1, 0], [0, 0, 1], [1, 1, 1]}]; gridtets := proc(i,j,k,A,dx,dy,dz,a1,a2,a3) # return sequence of candidate tetrahedra from a grid cell local L, sigs, t, v; for L to 2 do sigs[L] := map(signum, {A[i,j,k,L],A[i+1,j,k,L],A[i+1,j+1,k,L], A[i,j+1,k,L],A[i,j,k+1,L],A[i+1,j,k+1,L],A[i+1,j+1,k+1,L], A[i,j+1,k+1,L]}); end do; if sigs[1]={1} or sigs[1]={-1} or sigs[2]={1} or sigs[2]={-1} then return; end if; seq(qtet1(seq([a1+(i+t[v][1])*dx,a2+(j+t[v][2])*dy,a3+(k+t[v][3])*dz], v=1..4),seq([seq](A[i+t[v][1],j+t[v][2],k+t[v][3],L],L=1..2),v=1..4)), t = cubetets); end proc; end module: #savelib('surfintplot');