#Table structure: # VP_CIRCLE, VP_TREE, VP_BIPARTITE, VP_SPRING, VP_FIXED, VP_USER, VP_DEFAULT::[vpos, lpos], vpos::list, lpos::list([numeric, numeric]) # V_COLOR # V_RADIUS::numeric # E_THICKNESS::Matrix # E_COLOR::Matrix macro( VP_FIXED="draw-pos-fixed", VP_CIRCLE="draw-pos-circular", VP_TREE="draw-pos-tree", VP_BIPARTITE="draw-pos-bipartite", VP_SPRING="draw-pos-spring", VP_USER="draw-pos-user", VP_DEFAULT="draw-pos-default" ): getdir := proc(G) option inline; op(1,G); end proc: getwt := proc(G) option inline; op(2,G); end proc: vlist := proc(G) option inline; op(3,G); end proc: listn := proc(G) option inline; op(4,G); end proc: ginfo := proc(G) option inline; op(5,G); end proc: eweight := proc(G) option inline; op(6,G); end proc: getops := proc(G) option inline; op(1..6,G); end proc: listint := proc(v) option inline; [$1..v]; end proc: seqint := proc(v) option inline; $1..v; end proc: makeweights := proc(G,E) option inline; subsop(6=E,2=weighted,G) end: makevertices := proc(G,V) option inline; subsop(3=V,G) end: macro( VERTEXTYPE = {integer, symbol, string, indexed} ): macro( COLORTYPE = {symbol, specfunc(anything,COLOR)} ): macro( vectorplot3d = `GraphTheory/vectorplot3d` ): macro( vertex_r = 0.03, pagewidth = 1, pageheight = 1, #Mahdi arrow_h_l = 0.02771638598, arrow_v_l = 0.01148050298): #Mahdi GraphDrawing := module() export HighlightEdges, HighlightSubgraph, HighlightTrail, HighlightVertex, DrawGraph, DrawNetwork; option package; local IsBiPartite, GetConnectedGraphs, GetPlotObjs, GetTransVect, GetVtxRadius, SetBipartiteLayout, SetDefaultLayout, SetCircleLayout, SetSpringLayout, SetTreeLayout, PlotSpringConnected, SetVertexPos, ZUniform, vectorplot2d, makeG, Findetxtloc, AnalyseOptions, plotG, projection, findmax, BuildProc3, edgeweight, BuildProc2, findmax2, makebox, Getedges, makeg, cnvxhull, U01; #---------------------------------------------------------------------------------- DrawGraph := proc(H) local s, G, stl, rt, dim, anmt, arg, rarg, larg, layouts, n, i, objs, T, cols, res, rdraw, shwts, shlbs; G := H; stl := 'none'; rt := -1; dim := 2; anmt := false; rdraw := false; shwts := `if`(GraphTheory:-NumberOfEdges(G) < 46 , true, false); #AE shlbs := `if`(GraphTheory:-NumberOfVertices(G) < 100 , true, false); #MBM for i from 2 to nargs do arg := args[i]; if type(arg, equation) then larg, rarg := op(arg); if larg = 'style' then if not rarg in {'bipartite', 'circle', 'default', 'spring', 'tree'} then error "argument %1 (optional) expected to be of the form style=bipartite/circle/default/spring/tree", i; end if; stl := rarg; elif larg = 'showweights' then #AE if rarg in {true,false} then shwts:=rarg; else error "expected showweights = true/false but recieved %1",arg fi; elif larg = 'showlabels' then #MBM if rarg in {true,false} then shlbs:=rarg; else error "expected showlabels= true/false but recieved %1",arg fi; elif larg = 'root' then if not rarg in vlist(G) then error "'%1' is not a vertex of the input graph", rarg; end if; rt := GraphTheory:-GraphInfo:-LabelToInteger(G)[rarg]; elif larg = 'dimension' then if rarg in {2, 3} then dim := rarg; else error "dimension should be 2 or 3"; end if; elif larg = 'animate' then if rarg in {true, false} then anmt := rarg; else error "the option animate is either true or false"; end if; else error "argument %1 (optional) expected to be of the form 'style/rt/dimension/animate=...'", i; end if; else if arg = 'redraw' then rdraw := true; elif arg = 'showweights' then shwts := true; elif arg = 'showlabels' then shlbs := true; else error "argument %1 (optional) is wrong.",i; end if; end if; end do; if stl = 'none' then if GraphTheory:-GetVPos(G, VP_USER) <> [] then stl := 'user'; elif GraphTheory:-GetVPos(G, VP_FIXED) <> [] then stl := 'fixed'; else stl := 'default'; end if; end if; if stl = 'tree' and rt =-1 then rt := 1; elif stl <> 'tree' and rt <> -1 then error "set style to 'tree' or remove '%1'", arg; end if; if dim = 2 then if stl <> 'spring' and anmt = true then error "'animate' option can be set only when style='spring' or dimension=3";#AE end if; else if stl <> 'spring' and stl <> 'default' then error "only 'spring' and 'default' styles are allowed for 3D view of the graph"; end if; end if; if dim = 2 then if stl = 'user' then s := VP_USER; elif stl = 'fixed' then s := VP_FIXED; elif stl = 'default' then SetDefaultLayout(G); s := VP_DEFAULT; elif stl = 'tree' then if not GraphTheory:-IsTree(G) then error "input graph is not a tree"; end if; SetTreeLayout(G, rt); s := VP_TREE; elif stl = 'circle' then SetCircleLayout(G, 0.5, 0); s := VP_CIRCLE; elif stl = 'bipartite' then res := IsBiPartite(G); if res = false then error "input graph is not bipartite"; end if; SetBipartiteLayout(G); s := VP_BIPARTITE; elif stl = 'spring' then return SetSpringLayout(G, dim, anmt, rdraw,shwts); else error "unexpected style"; end if; elif dim = 3 then return SetSpringLayout(G, dim, anmt, rdraw,shwts); end if; PLOT(GetPlotObjs(G, s, 'showweights'=shwts, 'showlabels'=shlbs), AXESSTYLE(NONE)); #AE end: #------------------------------------------------------------------------------- SetVertexPos := proc(G::GRAPHLN, v::{integer, list}, pos::list) local L, vpos, i, n; vpos := GraphTheory:-GetVPos(G, VP_USER); n := nops(vlist(G)); if vpos = [] then vpos := [seq([-i, -i], i = 1..n)]; end if; vpos := Array(1..n,vpos); # MBM L := GraphTheory:-GraphInfo:-LabelToInteger(G); if type(v, 'integer') then if not member(v, vlist(G)) then error "vertex %1 not in graph", v; end if; vpos[L[v]] := pos; else if nops(pos) <> nops(v) then error "list of vertex positions should be of size %1", nops(v); end if; for i to nops(v) do if not member(v[i], vlist(G)) then error "vertex %1 not in graph", v[i]; end if; vpos[L[v[i]]] := pos[i]; end do; end if; GraphTheory:-SetVPos(G, VP_USER, vpos); return; end: #------------------------------------------------------------------------------- HighlightVertex := proc( G::GRAPHLN, v::{VERTEXTYPE, set(VERTEXTYPE), list(VERTEXTYPE)}, c::{COLORTYPE, list(COLORTYPE)}) local SetColor, color, vcolor, L, n, i; SetColor := proc(v::{set, list}, c); if not {op(v)} subset convert(vlist(G), set) then error "2nd argument expected to be a subset of the vertices of the graph, but received %1", v; end if; for i to nops(v) do vcolor[L[v[i]]] := `if`(type(c, list), c[i], c); end do; end; if nargs = 3 then if type(c, list) then color := [seq(`plot/color`(c[i]), i=1..nops(c))]; else color := `plot/color`(c); end if; else if type(v, list) then n := nops(v); color := ['red, blue, green, black, gold, pink, navy, coral, cyan, brown, gray, khaki, magenta, maroon, orange, plum, sienna, tan, turquoise, violet, wheat, white, yellow, aquamarine', seq(COLOR(HUE, i/(n-24)), i = 1..n-24)]; color := [seq(`plot/color`(color[i]), i = 1..n)]; else color := `plot/color`('green'); end if; end if; if type(color, list) then if not type(v, list) then error "only one color expected but received %1", c; elif nops(v) <> nops(color) then error "a list of %1 colors expected but received %2", nops(v), c; end if; else if type(v, list) then # assign them all the same color color := [color$nops(v)]; end if; end if; vcolor := GraphTheory:-GetVColor(G); vcolor := Array(vcolor); # MBM L := GraphTheory:-GraphInfo:-LabelToInteger(G); if type(v, list) or type(v, set) then SetColor(v, color); else SetColor({v}, color); end if; GraphTheory:-SetVColor(G, vcolor); return; end; #------------------------------------------------------------------------------- HighlightTrail := proc( G::GRAPHLN, T::{list(VERTEXTYPE),specfunc(VERTEXTYPE,Trail)} ) local E,i; E := { seq({op(i..i+1,T)}, i=1..nops(T)-1) }; HighlightEdges( G, E, args[3..nargs] ) end; #------------------------------------------------------------------------------- HighlightSubgraph := proc( G::GRAPHLN, H::{GRAPHLN,list(GRAPHLN),set(GRAPHLN)}, C1::COLORTYPE, C2::COLORTYPE ) if nargs=2 then return HighlightSubgraph(G, H, 'red', 'green'); fi; if nargs=3 then error "missing 4th argument: vertex color" fi; HighlightEdges( G, GraphTheory:-Edges(H), C1 ); HighlightVertex( G, GraphTheory:-Vertices(H), C2 ); end; #------------------------------------------------------------------------------- HighlightEdges := proc( G::GRAPHLN, s::{GRAPHLN, set(GRAPHLN), list(GRAPHLN), specfunc(VERTEXTYPE,Trail), set(VERTEXTYPE), set(set(VERTEXTYPE)), list(set(VERTEXTYPE)), list(VERTEXTYPE), set(list(VERTEXTYPE)), list(list(VERTEXTYPE))}, c::{COLORTYPE, list(COLORTYPE)}) local IsOfTypeList, IsOfTypeSet, SetColor, color, n, i, GEdges, cmtx, tmtx, L, E; if type(s,specfunc(VERTEXTYPE,'Trail')) then # just convert to a set of edges E := {seq({op(i,s),op(i+1,s)}, i=1..nops(s)-1)}; return HighlightEdges(G,E,args[3..nargs]); fi; IsOfTypeList := proc(s); return type(s, list) and not type(s, list({integer, symbol})); end; IsOfTypeSet := proc(s); return type(s, set) and not type(s, set({integer, symbol})); end; SetColor := proc(s::{set, list, GRAPHLN}, c) local e, edges; edges := `if`(type(s, GRAPHLN), GraphTheory:-Edges(s), s); if not {op(edges)} subset GEdges then error "2nd argument expected to be a subset of the edges of the graph or a subgraph of the graph, but received %1", s; end if; for i to nops(edges) do e := L[edges[i][1]], L[edges[i][2]]; cmtx[e] := `if`(type(c, list), c[i], c); tmtx[e] := 3; end do; end; if nargs = 3 then if type(c, list) then color := [seq(`plot/color`(c[i]), i=1..nops(c))]; else color := `plot/color`(c); end if; else if IsOfTypeList(s) then n := nops(s); color := ['red, blue, green, black, gold, pink, navy, coral, cyan, brown, gray, khaki, magenta, maroon, orange, plum, sienna, tan, turquoise, violet, wheat, white, yellow, aquamarine', seq(COLOR(HUE, i/(n-24)), i = 1..n-24)]; color := [seq(`plot/color`(color[i]), i = 1..n)]; else color := `plot/color`('red'); end if; end if; if type(color, list) then if not IsOfTypeList(s) then error "only one color expected but received %1", c; elif nops(s) <> nops(color) then error "a list of %1 colors expected but received %2", nops(s), c; end if; else if IsOfTypeList(s) then error "a list of colors expected but received %1", c; end if; end if; GEdges := GraphTheory:-Edges(G); cmtx := GraphTheory:-GetEdgesColor(G); #tmtx := GraphTheory:-GetEdgesThickness(G); L := GraphTheory:-GraphInfo:-LabelToInteger(G); if IsOfTypeList(s) then if type(s, list(GRAPHLN)) then for i to nops(s) do SetColor(s[i], color[i]); end do; else SetColor(s, color); end if; elif IsOfTypeSet(s) then if type(s, set(GRAPHLN)) then for i to nops(s) do SetColor(s[i], color); end do; else SetColor(s, color); end if; else if type(s, GRAPHLN) then SetColor(s, color); else SetColor({s}, color); end if; end if; GraphTheory:-SetEdgesColor(G, cmtx); return; end; #------------------------------------------------------------------------------- IsBiPartite := proc(G::GRAPHLN) local DFS, a, A, i, n, isBiPart; DFS := proc(G::GRAPHLN, v::integer) local u; for u in A[v] do if a[u] = 0 then a[u] := -a[v]; if not DFS(G, u) then return false; end if; elif a[u] = a[v] then return false; end if; end do; return true; end; n := GraphTheory:-NumberOfVertices(G); A := listn(G); a := Array(1..n); # a[v] is 1 if v is in the first set and -1 if in the second set. otherwsise v has not been visited yet and a[v] is 0. for i to n do if a[i] = 0 then a[i] := 1; if not DFS(G, i) then return false; end if; end if; end do; return a; end; #------------------------------------------------------------------------------- SetBipartiteLayout := proc(G::GRAPHLN) local GetHorizDist, res, n1, n2, k1, k2, n, h1, h2, y1, y2, i, d1, d2, r, labeled, max, vpos, lpos, pos, isBiPart, arr::Array; GetHorizDist := proc(n); if n <= 1 then return pagewidth/2, pagewidth/2; else return pagewidth/(n - 1), 0; end if; end proc; if GraphTheory:-GetVPos(G, VP_BIPARTITE) <> [] then return GraphTheory:-GetVPos(G, VP_BIPARTITE); end if; res := IsBiPartite(G); if res = false then return false; end if; arr := res; n := ArrayNumElems(arr); n1, n2 := 0, 0; for i to n do if arr[i] = 1 then n1 := n1 + 1; else n2 := n2 + 1; end if; end do; max := 15; r := `if`(n1>max or n2>max, 0, vertex_r); labeled := (n1 < 26) and (n2 < 26); h1, d1 := GetHorizDist(n1); h2, d2 := GetHorizDist(n2); y1 := `if`(n1 > n2, 0, 1); y2 := 1 - y1; vpos := Array(1..n,[seq([-i, -i], i = 1..n)]); # MBM k1, k2 := 1, 1; for i to n do if arr[i] = 1 then vpos[i] := [(k1-1)*h1 + d1, y1]; k1 := k1 + 1; else vpos[i] := [(k2-1)*h2 + d2, y2]; k2 := k2 + 1; end if; end do; #whats the point of this if we don't return it? AE lpos := []; if r = 0 then lpos := [seq([vpos[i][1], `if`(vpos[i][2] = 0.03, 0, 1)], i = 1..n)]; end if; # GraphTheory:-SetLPos(G, VP_BIPARTITE, lpos); GraphTheory:-SetVPos(G, VP_BIPARTITE, vpos); return vpos; end proc: #------------------------------------------------------------------------------- GetConnectedGraphs := proc(G::GRAPHLN) local n, a, singletons, comps, vset, i, Subgraph; Subgraph := proc(vset) local H; if (nops(vset) = 0) then return NULL; end if; H := GraphTheory:-InducedSubgraph(G, vset); return H; end: comps := GraphTheory:-ConnectedComponents(G); # put all the singletons (isolated vertices) in one componnent singletons := {seq( `if`(nops(vset)=1, vset[1], NULL), vset=comps)}; comps := [seq( `if`(nops(vset)=1, NULL, vset), vset=comps)]; n := nops(comps) + 1; a := Array(1..n); a[1] := Subgraph(singletons); for i from 2 to n do a[i] := Subgraph(comps[i-1]); end do; return a; end: #------------------------------------------------------------------------------- GetTransVect := proc(index, col, d) local x, y; x, y := (index mod col) * (1 + d), -trunc(index/col) * (1 + d); return [x, y]; end; #------------------------------------------------------------------------------- SetDefaultLayout := proc(G::GRAPHLN) local a, UB, H, i, k, vpos, x, y, m, n, p, j, cols, vp, res, L, Translate; Translate := proc(L::{list,Array}) local i, n; if L::list then n := nops(L) else n := ArrayNumElems(L); fi; [seq(L[i] + GetTransVect(k, cols, .2), i = 1..n)]; end; vp := GraphTheory:-GetVPos(G, VP_DEFAULT); if vp <> [] then return vp end if; a := GetConnectedGraphs(G); m := ArrayNumElems(a); UB := a[1]; n := GraphTheory:-NumberOfVertices(G); vpos := [seq([-i, -i], i = 1..n)]; vpos := Array(1..n,vpos); # MBM L := GraphTheory:-GraphInfo:-LabelToInteger(G); k := 0; if UB <> NULL then vp := SetCircleLayout(UB, 0.5, 0); for i to nops(vp) do vpos[L[vlist(UB)[i]]] := vp[i]; end do; k := 1; end if; if m=3 then cols := 3 else cols := ceil(sqrt(m)) fi; for i from 2 to m do H := a[i]; res := IsBiPartite(H); if not GraphTheory:-IsDirected(H) and GraphTheory:-IsTree(H) then vp := SetTreeLayout(H, 1); #may not return list elif res <> false then vp := SetBipartiteLayout(H); else vp := SetCircleLayout(H, 0.5, 0); end if; vp := Translate(vp); k := k + 1; for j to nops(vp) do vpos[L[vlist(H)[j]]] := vp[j]; end do; end do; GraphTheory:-SetVPos(G, VP_DEFAULT, vpos); return vpos; end: #------------------------------------------------------------------------------ PlotSpringConnected := proc(G, H, dim, anmt, redraw,shwts) #AE local a, b, i, j, n, out, SOL, directed, nfr, U; a:=map(proc(x) [op(x)] end, convert(listn(H),list)): U := GraphTheory:-UnderlyingGraph(H); b:=map(proc(x) [op(x)] end, convert(listn(U),list)): n:=nops(a); out := `if`(anmt=true, 'anim', 'plotsol'); if n = 1 then SOL := `if`(dim=2, [[0,0]], [[0,0,0]]); elif redraw = true or anmt = true then # SOL := makeg(1/1100.0, n, 700 + 40 * n, b, out, 12, dim, 50); #MAHDI nfr := `if`(anmt = true, 50 , 0); # Try Allan's code first SOL := traperror( `dsolve/numeric/SpringIntegrate`(b , 1 , dim , nfr , 20000) ); if SOL=lasterror then # it failed, so try Moe's code # MBM: This lprint is left in so that if Allan's code for Maple 11 # MBM: below is not found then one of the tests will not pass. lprint("Could not find Allan's Maple 11 code: executing Moe's Maple 10 code."); SOL := makeg(1/1100.0, n, 700 + 40 * n, b, out, 12, dim, 50); #this is the place where animate=true option fails. if dim=2 then # convert output from dsolve/numeric to list of 2D points if anmt = false then #AE SOL:=[seq([rhs(SOL[2*i]),rhs(SOL[2*n+2*i])],i=1..n)]; else SOL:=[seq([seq([SOL[j,2*i],SOL[j,2*n+2*i]],i=1..n)],j=1..50)]; fi; elif dim=3 then if anmt = false then #AE SOL:=[seq([rhs(SOL[2*i]),rhs(SOL[2*n+2*i]),rhs(SOL[4*n+2*i])],i=1..n)]; else SOL:=[seq([seq([SOL[j,2*i],SOL[j,2*n+2*i],SOL[j,4*n+2*i]],i=1..n)],j=1..50)]; fi; end if; else # convert Allan's Array(1..n,1..2) to a list of lists # lprint("Executing Allan's code"); SOL:=convert(SOL,listlist); fi; #----- end if; directed := `if`(GraphTheory:-IsDirected(H), 1, 0); plotG(G, SOL, n, dim, a, b, directed, 50, out, vlist(H), eweight(H), 0, 1, `plot/color`('black'), `plot/color`('yellow'), `plot/color`('green'), `plot/color`('black'), `plot/color`('blue'), redraw, shwts); #AE end: #------------------------------------------------------------------------------ SetSpringLayout := proc(G, dim, anmt, redraw,shwts) #AE local arr, UB, L, m, r, n, vl, vc, gc, gt, ec, et, i, j,vpos; if not GraphTheory:-IsConnected(G) then error "sorry, for style=spring and dimension=3, the graph must be connected."; fi; #AE: make sure we are in the right dimension and that GetVPos executes properly #AE: unfortunatly there is only one place to put these coordinates so when #AE: graph is drawn in 2 then 3 dimensions it is forced to redraw. vpos:=GraphTheory:-GraphInfo:-GetVPos(G,VP_SPRING); if vpos=[] or dim <> nops(vpos[1]) or redraw then r:= true; else r:=false; fi; return PlotSpringConnected(G, G, dim, anmt, r, shwts); #AE # arr := GetConnectedGraphs(G); # m := ArrayNumElems(arr); # UB := arr[1]; # # if UB <> NULL then # SetCircleLayout(UB, 0.5, 0); # L := GraphTheory:-GraphInfo:-LabelToInteger(G); # n := GraphTheory:-NumberOfVertices(UB); # vl := vlist(UB); # # vc := GraphTheory:-GetVColor(G); # vc := [seq(vc[L[vl[i]]], i = 1..n)]; # GraphTheory:-SetVColor(UB, vc); # # gc := GraphTheory:-GetEdgesColor(G); # gt := GraphTheory:-GetEdgesThickness(G); # ec := GraphTheory:-GetEdgesColor(UB); # et := GraphTheory:-GetEdgesThickness(UB); # for i to n do # for j to n do # ec[i, j] := gc[L[vl[i]], L[vl[j]]]; # et[i, j] := gt[L[vl[i]], L[vl[j]]]; # end do; # end do; # GraphTheory:-SetEdgesColor(UB, ec); # GraphTheory:-SetEdgesThickness(UB, et); # # return(PLOT(GetPlotObjs(UB, VP_CIRCLE), AXESSTYLE(NONE))); # end if; # # r := `if`(GraphTheory:-GetVPos(G, VP_SPRING) = [], true, false) or redraw; # return PlotSpringConnected(G, arr[2], dim, anmt, r); end: #------------------------------------------------------------------------------ SetTreeLayout := proc(G::GRAPHLN, root::integer) local maxwidth, width, height, i, r, nolabel, n, vpos, lpos, vertdist, horizdist, GetChildren, GetSubTreeLayout, GetWidthAndHeight; GetChildren := proc(node, parent); convert(listn(G)[node] minus {parent}, 'list'); end: GetWidthAndHeight := proc(root, parent) local children, width, maxheight, child, w, h; children := GetChildren(root, parent); if nops(children) = 0 then return 1, 1; end if; width, maxheight := 0, 0; for child in children do w, h := GetWidthAndHeight(child, root); width := width + w; if h > maxheight then maxheight := h; end if; end do; return width, maxheight + 1; end: GetSubTreeLayout := proc(root, parent, next_leaf_x, level) local y, children, i, n, x; x := next_leaf_x; y := pageheight+level * vertdist; children := GetChildren(root, parent); n := nops(children); if n = 0 then #root is a leaf vpos[root] := [x, y]; # lpos[root] := [x, `if`(r <> 0, y, y - vertex_r)]; return x + horizdist; end if; for i from 1 to n do x := GetSubTreeLayout(children[i], root, x, level + 1); end do; vpos[root] := [(vpos[children[1]][1] + vpos[children[n]][1])/2, y]; # lpos[root] := [`if`(r <> 0, vpos[root][1], vpos[root][1] - vertex_r), y]; return x; end: if GraphTheory:-GetVPos(G, VP_TREE) <> [] then return GraphTheory:-GetVPos(G, VP_TREE); end if; maxwidth := 17; width, height := GetWidthAndHeight(root, 0); horizdist := `if`(width = 1, 0, pagewidth /(width - 1)); vertdist := `if`(height = 1, 0, -pageheight/(height - 1)); r := `if`(width <= maxwidth, vertex_r, 0); nolabel := `if`(width <= 25, false, true); n := GraphTheory:-NumberOfVertices(G); vpos := [seq([-i, -i], i = 1..n)]; vpos := Array(1..n,vpos); # MBM # lpos := [seq([-i, -i], i = 1..n)]; GetSubTreeLayout(root, 0, 0, 0); # GraphTheory:-SetLPos(G, VP_TREE, lpos); GraphTheory:-SetVPos(G, VP_TREE, vpos); return vpos; end: #------------------------------------------------------------------------------- SetCircleLayout := proc(G::GRAPHLN, graph_r, phase) local n, nolabel, label_r, c, i, angle, unit_x, unit_y, lpos, vpos, l; vpos := GraphTheory:-GetVPos(G, VP_CIRCLE); if vpos <> [] then return vpos; end if; n := GraphTheory:-NumberOfVertices(G); nolabel := `if`(n <= 80, false, true); label_r := `if`(n > 35, 1.1, 1); c := Array(1..n); #an array storing centers of the vertesis # l := Array(1..n); for i to n do angle := evalf(phase + Pi/2 * (1 - 4*(i - 1)/n)); unit_x, unit_y := evalf(cos(angle)), evalf(sin(angle)); #(unit_x, unit_y) is a point on the unit circle c[i] := [graph_r*(unit_x + 1), graph_r*(unit_y + 1)]; #a point on a circle with center (graph_r, graph_r) and radius graph_r # l[i] := [graph_r*(unit_x*label_r + 1), graph_r*(unit_y*label_r - 1)]; end do; vpos := [seq(c[i], i = 1..n)]; # lpos := []; # if label_r <> 1 then lpos := [seq(l[i], i = 1..n)]; end if; GraphTheory:-SetVPos(G, VP_CIRCLE, vpos); # GraphTheory:-SetLPos(G, VP_CIRCLE, lpos); return vpos; end: #------------------------------------------------------------------------------- GetVtxRadius := proc(G::GRAPHLN)::numeric; local r; r := ginfo(G)[V_RADIUS]; if type(r, indexed) then return vertex_r; end if; return r; end; #------------------------------------------------------------------------------ GetPlotObjs := proc(G::GRAPHLN, style::string) local n, m, vpos, lpos, i, v, l, e, vobjs, lobjs, eobjs, GetVtx, GetEdge, ecolor, ethickness, u, directed, arrow, GetXY, GetXY2, DirectedEdge, NonStandardEdge, vcolor, polyStyle, interface, Line, GetResolution, x, y, scaling, res, wobjs, weighted, M, r, mp, label, d, f, vtxfnt, vtxwid, m1,m2,k,N,edgfnt, shwts, shlbs; #---------- # (x, y) is the coordinate of the vertex center with height h and width w GetVtx := proc(xy, h, w) local x, y; x, y := xy[1], xy[2]; return [[x + w, y + h], [x + w, y - h], [x - w, y - h], [x - w, y + h]]; end: #---------- GetResolution := proc(p::list) local x, y, dx, dy; x := seq(p[i][1], i = 1..nops(p)); dx := max(x) - min(x); y := seq(p[i][2], i = 1..nops(p)); dy := max(y) - min(y); if nops(p) = 1 then return 1; end if; return max(dx, dy); end; #---------- GetXY2 := proc(p, q, a) local b, c, d; d := evalf(sqrt((p[1] - q[1])^2 + (p[2] - q[2])^2)); b := a*(q[2] - p[2])/d; c := a*(q[1] - p[1])/d; return p + [b, -c], q + [b, -c]; end; #---------- DirectedEdge := proc(f, t, type, interface) local Arrow; Arrow := proc(p1, p2, arrType) local x, y, dir, orth, a1, a2, NORM, GetArrow, rat, l1, l2; NORM := proc(p) local l; l := evalf(sqrt(p[1]^2+p[2]^2)); if l = 0 then return p; fi; return p / l; end: rat := `if`(arrType = SINGLE, 4/5, .4); l1, l2 := arrow_h_l*res, arrow_v_l*res; x := evalf(rat*p1[1]+(1-rat)*p2[1]); y := evalf(rat*p1[2]+(1-rat)*p2[2]); dir := NORM(p2 - p1); orth := l2*[-dir[2],dir[1]]; dir := -l1*dir; a1 := dir+orth; a2 := dir-orth; return Line([x,y],[x,y]+a1, interface), Line([x,y],[x,y]+a2, interface); end: return Line(f, t, interface), Arrow(f, t, type); end; #---------- GetXY := proc(p, q, a) local b, c, d; d := evalf(sqrt((p[1] - q[1])^2 + (p[2] - q[2])^2)); b := a*(q[2] - p[2])/d; c := a*(q[1] - p[1])/d; return p + [b, - c], p - [b, -c]; end; #---------- #NonStandardEdge := proc(p1, p2) local a, q1, q2, q3, q4; # #a := 0.0015; # too small # q1, q2 := GetXY(p1, p2, a); # q3, q4 := GetXY(p2, p1, a); # return [q1, q2, q3, q4]; #end proc; # MBM: needs to be scaled to the resolution of the plot NonStandardEdge := proc(p,q) local a, u, v, NORM; NORM := proc(p) evalf(sqrt(p[1]^2+p[2]^2)); end: u := q-p; v := [-u[2],u[1]]; # rotate anti-clockwise 90 degrees a := 0.008*res/NORM(v)/evalf(n^(1/4)); # n = |V| and res is the plot resolution v := a*v; [p+v, q+v, q-v, p-v]; end: #---------- Line := proc(p1, p2, interface); return `if`(interface = 'standard', [p1, p2], NonStandardEdge(p1, p2)); end; #---------- GetEdge := proc(u, v, interface) local a, b, cmtx; if not directed then return Line(vpos[u], vpos[v], interface); end if; if member(u, e[v]) then cmtx := GraphTheory:-GetEdgesColor(G); if cmtx[u, v] <> cmtx[v, u] then a, b := GetXY2(vpos[u], vpos[v], 0.015); return DirectedEdge(a, b, SINGLE, interface); else return DirectedEdge(vpos[u], vpos[v], DOUBLE, interface); end if; else return DirectedEdge(vpos[u], vpos[v], SINGLE, interface); end if; end; #---------- start GetPlotObjs shwts := true; shlbs := true; for i from 3 to nargs do if type(args[i],identical('showweights')=truefalse) then shwts := rhs(args[i]); elif type(args[i],identical('showlabels')=truefalse) then shlbs := rhs(args[i]); else error "unknown option: %1", args[i]; fi; end: n := GraphTheory:-NumberOfVertices(G); N := GraphTheory:-NumberOfEdges(G); vpos := GraphTheory:-GetVPos(G, style); for i to n do if not type(vpos[i], [numeric,numeric]) then error "position should be specified for all vertices."; end if; end do; lpos := GraphTheory:-GetLPos(G, style); l := vlist(G); #lobjs := seq(TEXT(lpos[i], convert(l[i], 'string')), i =1..n); lobjs := [seq( convert(l[i],string), i=1..n )]; m := max( seq(length(l), l=lobjs) ); # maximum string length if m = 1 then m := 1 else m := 2/3+m/3.0; fi; if shlbs then if n>50 then f:=9 elif n>20 then f:=10 elif n>9 then f:=11 else f:=12 fi; vtxfnt := FONT(HELVETICA,BOLD,f); lobjs := seq( TEXT(lpos[i],lobjs[i],vtxfnt), i=1..n ); else lobjs := NULL; fi; res := GetResolution(vpos); vtxwid := vertex_r * res * (1.0 - min(0.5,n/100)); if shlbs=false then m := 1; vtxwid := 0.8*vtxwid; fi; # MBM shrink it if no label v := seq(GetVtx(vpos[i], vtxwid, m*vtxwid), i = 1..n); vcolor := GraphTheory:-GetVColor(G); if shlbs then vobjs := POLYGONS(v, COLOR(RGB, seq(op(vcolor[i])[2..4], i = 1..n)), STYLE(PATCHNOGRID)); else vobjs := POLYGONS(v, COLOR(RGB, seq(op(vcolor[i])[2..4], i = 1..n)), STYLE(PATCH)); fi; e := listn(G); directed := GraphTheory:-IsDirected(G); weighted:=GraphTheory:-IsWeighted(G); #AE if weighted and shwts=true then #MBM edgfnt := FONT(HELVETICA,`if`(N>20,`if`(N>45,9,10),11)); wobjs := table(); k := 0; M := GraphTheory:-WeightMatrix(G); r:=.60; for u to n do for v in e[u] do if directed or v>u then mp:= [evalf(((1-r)*vpos[u,1]+r*vpos[v,1])),evalf(((1-r)*vpos[u,2]+r*vpos[v,2]))]: label:=M[u,v]; if type(label,integer) then label:= sprintf("%d",label); else label:=sprintf("%0.3g",label); fi; m2:=evalf(vpos[u,1]-vpos[v,1]); if abs(m2) < .01*res then m:= infinity; mp:= mp + [1,0]*res*0.007; else m1:=vpos[u,2]-vpos[v,2]; if abs(m1) < 0.01*res then #put the label above the arrow if we have to if directed and u in e[v] and vpos[u][1]>vpos[v][1] then mp:= mp + [0,1]*res*0.015; else mp:= mp + [0,1]*res*0.005; fi; m:=0; else m:=evalf(m1/m2); mp := mp + res*`if`(m>0,0.009,0.002)*[abs(m1),abs(m2)*`if`(m>0,-1,1)]/sqrt(m1^2+m2^2); fi; fi; label:= TEXT(mp,label,ALIGNRIGHT,`if`(m>0,ALIGNBELOW,ALIGNABOVE),edgfnt); k:=k+1; wobjs[k]:=label; fi; od; od; wobjs := seq( wobjs[u], u=1..k ); else wobjs:=NULL fi; ecolor := GraphTheory:-GetEdgesColor(G); ethickness := GraphTheory:-GetEdgesThickness(G); eobjs := table(); M := 0; for u to n do for v in e[u] do if directed or u>v then if IsWorksheetInterface('Standard') then polyStyle := STYLE(LINE); interface := 'standard'; else polyStyle := STYLE(PATCHNOGRID); interface := 'nonstandard'; end if; M := M+1; eobjs[M] := POLYGONS(GetEdge(u, v, interface), ecolor[u, v], THICKNESS(ethickness[u, v]), polyStyle); fi; end do; end do; eobjs := seq( eobjs[u], u=1..M ); if n > 1 then scaling := SCALING(CONSTRAINED); elif n = 0 then return TEXT([0,0],"empty graph") else x, y := op(vpos[1]); scaling := VIEW(x-.5..x+.5, y-.5..y+.5); end if; return (vobjs, lobjs, eobjs, wobjs, scaling); end; DrawNetwork := proc(G :: GRAPHLN) local GD,S,T,levels,rlevels,l,i,R,k,x,y,n,rs,ta,ad,maxh,temp,xmax,ymax,hstep,isn,Tb,ma,vl,ps,LS,SS,TT,hor, vp, drawoptions; Tb := GraphTheory:-GraphInfo:-LabelToInteger(G); hor := false; drawoptions := NULL; for i from 2 to nargs do if args[i]='horizontal' then hor := true; elif args[i]='vertical' then hor := false; elif args[i]='showweights' then drawoptions := drawoptions,showweights elif args[i]='showlabels' then drawoptions := drawoptions,showlabels elif type(args[i],VERTEXTYPE) and type(Tb[args[i]],integer) then if assigned(S) then T := args[i] else S := args[i] fi; elif type(args[i],`=`) then drawoptions := drawoptions,args[i] else error "unrecognized argument: %1",args[i]; fi od; if assigned(S) and assigned(T) then isn := GraphTheory:-IsNetwork(G, S, T); if isn = false then error "graph is not a network with source %1 and sink %2.",S,T; fi; S := {Tb[S]}; T := {Tb[T]}; elif assigned(S) then error "sink not specified" else SS, TT := GraphTheory:-IsNetwork(G); if SS = {} or TT = {} then error "graph is not a network: please use DrawGraph."; fi; S := {}; T := {}; for i in SS do S := S union {Tb[i]}; od; for i in TT do T := T union {Tb[i]}; od; fi; Tb := ginfo(G); if not type(Tb, table) then assign(Tb, table()) fi; #M.G. GD := GraphTheory:-CopyGraph(G); if not type(Tb[VERTEX_POS], table) then n := nops(vlist(GD)): levels := Array(1..n): for i from 1 to n do levels[i] := infinity; od: for i in S do levels[i] := 1; od: for i in T do levels[i] := 0; od: rs := {seq(i,i=1..n)} minus S: ps := {seq(i,i=1..n)}; ad := listn(GD): for i from 1 to n do ps := ps minus ad[i]; od; ps := ps minus S; rs := rs minus ps; for i in ps do levels[i] := 2; od; k := 2: maxh := 0; ma := 0; LS := S union ps; while nops(rs) > 0 do R := {}; for i in LS do R := R union ad[i]; od; for i in R do if i in T then levels[i] := max(k , levels[i]); else levels[i] := min(k , levels[i]); ma := max(ma , levels[i]); fi; maxh := max(maxh , levels[i]); od; k := k + 1; LS := R; rs := rs minus R; od: maxh := max(ma + 1, maxh); for i in T do levels[i] := maxh; od: rlevels := Array(1..maxh): for i from 1 to n do temp := levels[i]; rlevels[temp] := rlevels[temp] + 1; od: ta := Array(1..maxh): ymax := 1.5: xmax := 1.5: hstep := - ymax / (maxh + 1): vl := vlist(GD); vp := Array(1..nops(vl)); for i from 1 to n do l := levels[i]; x := (ta[l]+1)* (xmax / (rlevels[l] + 1)); y := l*hstep; if not hor then vp[i] := [x,y]; else vp[i] := [-y,x]; fi; ta[l] := ta[l] + 1; od: fi; GraphTheory:-GraphInfo:-SetVPos(GD, VP_DEFAULT, convert(vp, list)); DrawGraph(GD,drawoptions); end proc; #Mohammad Ali Ebrahimi, Michael Monagan, Allan Wittkopf #-------------------------------------------------------------------------------- #makeG #AE: makeG doesn't seem to be called by any other routine in this document (GDv25.mpl) makeG:=proc(G::GRAPHLN) local i,n,SOL,a,T,M,dim,con,OUTPUT,stps,Arguments,dir,nodes,GG,b,nframes, edgelab,verlab,verlabcol,vercol,edlabcol,ED,ned,edgcol,L,arrcol; Arguments:=[seq(args[i],i = 2..nargs)]; GG:=GraphTheory['UnderlyingGraph'](G); ED:=GraphTheory['Edges'](GG); L:=GraphTheory:-GraphInfo:-LabelToInteger(G); ned:=nops(ED); a:=map(proc(x) [op(x)] end, convert(op(4,G),list)): b:=map(proc(x) [op(x)] end, convert(op(4,GG),list)): M:=op(6,G); nodes:=op(3,G): n:=nops(a); AnalyseOptions(Arguments,n,ned,L,a,'OUTPUT','dim','con','nframes','T','stps','dir','edgelab', 'verlab','verlabcol','vercol','edlabcol','arrcol','edgcol'); #finds the final position of the nodes SOL := makeg(con,n,T,b,OUTPUT,stps,dim,nframes); #Drawing Graph! plotG(G, SOL,n,dim,a,b,dir,nframes,OUTPUT,nodes,M,edgelab,verlab,verlabcol,vercol, edlabcol,arrcol,edgcol); #AE: missing shwts option but since makeG is never #AE: called it doesn't matter. end: #--------------------------------------------------------------------------------- #AE: AnalyseOptions is only called by makeG which is not called by anything. #AE: recomment removing both of these routines. AnalyseOptions:=proc(Arguments,n,ned,L,a,OUTPUT,dim,con,nframes,T,stps,dir,edgelab,verlab, verlabcol,vercol,edlabcol,arrcol,edgcol) local CON,DIM,output,l,r,opt,narg,i,N,Time,STPS,Dir,Edgelab,Verlab, Verlabcol,Vercol,k,j,Edlabcol,Arrcol,Edgcol,Col; output:='plotsol'; CON:=1/1100.0; DIM:=2; Dir:=0; N:=50; Edgelab:=0; Verlab:=1; Verlabcol:='COLOR(RGB,0,0,0)'; Edlabcol:='COLOR(RGB,0,1,0)'; Arrcol:='COLOR(RGB,0,0,0)'; Vercol:='COLOR(RGB,1,1,0)'; Edgcol:='COLOR(RGB,0,0,1)'; Time:=700+40*n; STPS:=12; narg := nops(Arguments); if narg > 0 then for i from 1 to narg do opt:=Arguments[i]; l := op(1,opt); r := op(2,opt); if l='dimension' then if r in {2,3} then DIM:=r; else error"dimension should be 2 or 3"; end if; elif l='animate' then if r=false then output:='plotsol'; elif r=true then output:='anim'; else error "the option animate is either true or false"; end if; elif l='springcon' then if type(r,numeric) and (r > 0.01 or r < 0) then error "the springcon cannot be less than 0.01 or negative"; elif type(r,numeric) and r <= 0.01 then CON:=r; else error "the springcon must be numeric"; end if; elif l='numframes' then if type(r,posint) then N:=r; else error"the number of frames must be a positive integer"; end if; elif l='time' then if type(r,numeric) and r>0 then Time:=r; else error"time must be a positive number"; end if; elif l='stepsize' then if type(r,numeric) and r>5 then STPS:=r; else error "the stepsize option should be a positive number grater than 5"; end if; elif l='directed' then if r=true then Dir:=1; elif r=false then Dir:=0; else error "bad argument for directed option"; end if; elif l='edgelable' then if r=true then Edgelab:=1; elif r=false then Edgelab:=0; else error "bad argument for the edgelable option"; end if; elif l='vertexlable' then if r=true then Verlab:=1; elif r=false then Verlab:=0; else error "bad argument for the vertexlable option"; end if; elif l='vertexlablecolor' then if type(r,name) and not(type(r,procedure)) then Verlabcol := `plot/color`(r); elif op(0,r) = 'COLOR' or op(0,r) = 'COLOUR' then Verlabcol := r; else error "wrong color info for the vertexlablecolor"; end if; elif l='edgelablecolor' then if type(r,name) and not(type(r,procedure)) then Edlabcol := `plot/color`(r); elif op(0,r) = 'COLOR' or op(0,r) = 'COLOUR' then Edlabcol := r; else error "wrong color info for the edgelablecolor"; end if; elif l='arrowscolor' then if type(r,name) and not(type(r,procedure)) then Arrcol := `plot/color`(r); elif op(0,r) = 'COLOR' or op(0,r) = 'COLOUR' then Arrcol := r; else error "wrong color info for the arrowscolor"; end if; elif l='vertexcolor' then if type(r,name) and not(type(r,procedure)) then Vercol := `plot/color`(r); elif op(0,r) = 'COLOR' or op(0,r) = 'COLOUR' then Vercol := r; elif type(r,set) then if nops(r)<>n then error "the number of colors should be the same as number of nodes in the vertexcolor option"; end if; for k from 1 to n do if not type(L[r[k][1]],integer) then error "the node you have entered is not valid, in the vertexcolor option",r[k][1]; end if; if type(r[k][2],name) and not(type(r[k][2],procedure)) then Col:=`plot/color`(r[k][2]); else Col:=r[k][2]; end if; Vercol[L[r[k][1]]]:=Col; end do; for k from 1 to n do if not assigned(Vercol[k]) then error "the must have entered one the nodes twice in the vertexcolor option"; end if; end do; else error "wrong color info"; end if; elif l='edgecolor' then if type(r,name) and not(type(r,procedure)) then Edgcol := `plot/color`(r); elif op(0,r) = 'COLOR' or op(0,r) = 'COLOUR' then Edgcol := r; elif type(r,set) then if nops(r)<>ned then error "the number of colors should be the same as number of edges, in the edgecolor option"; end if; for k from 1 to ned do if not type([L[r[k][1][1]],L[r[k][1][2]]],[integer,integer]) then error "the node you have entered is not valid",r[k][1]; end if; if type(r[k][2],name) and not(type(r[k][2],procedure)) then Col:=`plot/color`(r[k][2]); else Col:=r[k][2]; end if; Edgcol[L[r[k][1][1]],L[r[k][1][2]]]:=Col; Edgcol[L[r[k][1][2]],L[r[k][1][1]]]:=Col; end do; for k from 1 to n do for j in a[k] do if not assigned(Edgcol[k,j]) then error "the must have entered one the nodes twice in the edgecolor option"; end if; end do; end do; else error "wrong color info, in the edgecolor option"; end if; else error "wrong option"; end if; end do; end if; OUTPUT:=output; dim:=DIM; con:=CON; nframes:=N; T:=Time; stps:=STPS; dir:=Dir; edgelab:=Edgelab; verlab:=Verlab; verlabcol:=Verlabcol; vercol:=Vercol; edlabcol:=Edlabcol; arrcol:=Arrcol; edgcol:=Edgcol; end proc: #------------------------------------------------------------------------------- makeg:=proc(con,n,T,a,OUTPUT,stps,dim,N) local ps,SOL,des,d,i,sol,Sqr,j,ind; _Env_Use_Proc:=true; _Env_Use_RKF:=true; if dim=2 then Sqr:=ceil(sqrt(n)/2); else Sqr:=ceil(n^(1/3)/2); end if; ps:=evalf([seq(([seq(U01(Sqr),ind=1..dim)]),i=1..n)],4): if dim=2 then if _Env_Use_Proc=true then des := [BuildProc2(a,con)]: des := 'procedure'=op(1,des),'procvars'=op(2,des),'start'=0, 'initial'=Array(1..4*n,map(op,[seq([ps[i][1],0],i=1..n), seq([ps[i][2],0],i=1..n)])); else for i from 1 to n do d[i]:=nops(a[i])*diff('x'[i]('t'),'t','t')=-sqrt(n)*con*add('x'[i]('t')-'x'[j]('t'),j=a[i])-nops(a[i])^2/200*diff('x'[i]('t'),'t') +con*add(`if`(j<>i,('x'[i]('t')-'x'[j]('t'))/(('x'[i]('t')-'x'[j]('t'))^2+('y'[i]('t')-'y'[j]('t'))^2)^(2/2),0),j=1..n), nops(a[i])*diff('y'[i]('t'),'t','t')=-sqrt(n)*con*add('y'[i]('t')-'y'[j]('t'),j=a[i])-nops(a[i])^2/200*diff('y'[i]('t'),'t') +con*add(`if`(j<>i,('y'[i]('t')-'y'[j]('t'))/(('x'[i]('t')-'x'[j]('t'))^2+('y'[i]('t')-'y'[j]('t'))^2)^(2/2),0),j=1..n) ; end do; des:={seq(d[i],i=1..n), seq('x'[i](0)=ps[i][1],i=1..n), seq('y'[i](0)=ps[i][2],i=1..n),seq(D('y'[i])(0)=0,i=1..n), seq(D('x'[i])(0)=0,i=1..n)}; end if; elif dim=3 then if _Env_Use_Proc=true then des := [BuildProc3(a,con)]: des := 'procedure'=op(1,des),'procvars'=op(2,des),'start'=0, 'initial'=Array(1..6*n,map(op,[seq([ps[i][1],0],i=1..n), seq([ps[i][2],0],i=1..n), seq([ps[i][3],0],i=1..n)])); else for i from 1 to n do d[i]:=nops(a[i])*diff('x'[i]('t'),'t','t')=-sqrt(n)*con*add('x'[i]('t')-'x'[j]('t'),j=a[i])-nops(a[i])^2/100*diff('x'[i]('t'),'t') +con*add(`if`(j<>i,('x'[i]('t')-'x'[j]('t'))/(('x'[i]('t')-'x'[j]('t'))^2+('y'[i]('t')-'y'[j]('t'))^2 +('z'[i]('t')-'z'[j]('t'))^2),0),j=1..n), nops(a[i])*diff('y'[i]('t'),'t','t')=-sqrt(n)*con*add('y'[i]('t')-'y'[j]('t'),j=a[i])-nops(a[i])^2/100*diff('y'[i]('t'),'t') +con*add(`if`(j<>i,('y'[i]('t')-'y'[j]('t'))/(('x'[i]('t')-'x'[j]('t'))^2+('y'[i]('t')-'y'[j]('t'))^2 +('z'[i]('t')-'z'[j]('t'))^2),0),j=1..n), nops(a[i])*diff('z'[i]('t'),'t','t')=-sqrt(n)*con*add('z'[i]('t')-'z'[j]('t'),j=a[i])-nops(a[i])^2/100*diff('z'[i]('t'),'t') +con*add(`if`(j<>i,('z'[i]('t')-'z'[j]('t'))/(('x'[i]('t')-'x'[j]('t'))^2+('y'[i]('t')-'y'[j]('t'))^2 +('z'[i]('t')-'z'[j]('t'))^2),0),j=1..n); end do; des:={seq(d[i],i=1..n), seq('x'[i](0)=ps[i][1],i=1..n), seq('y'[i](0)=ps[i][2],i=1..n),seq(D('y'[i])(0)=0,i=1..n), seq(D('x'[i])(0)=0,i=1..n),seq('z'[i](0)=ps[i][3],i=1..n),seq(D('z'[i])(0)=0,i=1..n)}; if(false) then # changes for maping 4dim to 3dim ps:=evalf([seq(([seq(U01(Sqr),ind=1..4)]),i=1..n)],4): for i from 1 to n do d[i]:=nops(a[i])*diff('x'[i]('t'),'t','t')=-sqrt(n)*con*add('x'[i]('t')-'x'[j]('t'),j=a[i])-nops(a[i])^2/100*diff('x'[i]('t'),'t') +con*add(`if`(j<>i,('x'[i]('t')-'x'[j]('t'))/(('x'[i]('t')-'x'[j]('t'))^2+('y'[i]('t')-'y'[j]('t'))^2 +('z'[i]('t')-'z'[j]('t'))^2+(Q[i]('t')-Q[j]('t'))^2),0),j=1..n), nops(a[i])*diff('y'[i]('t'),'t','t')=-sqrt(n)*con*add('y'[i]('t')-'y'[j]('t'),j=a[i])-nops(a[i])^2/100*diff('y'[i]('t'),'t') +con*add(`if`(j<>i,('y'[i]('t')-'y'[j]('t'))/(('x'[i]('t')-'x'[j]('t'))^2+('y'[i]('t')-'y'[j]('t'))^2 +('z'[i]('t')-'z'[j]('t'))^2+(Q[i]('t')-Q[j]('t'))^2),0),j=1..n), nops(a[i])*diff('z'[i]('t'),'t','t')=-sqrt(n)*con*add('z'[i]('t')-z[j]('t'),j=a[i])-nops(a[i])^2/100*diff('z'[i]('t'),'t') +con*add(`if`(j<>i,('z'[i]('t')-'z'[j]('t'))/(('x'[i]('t')-'x'[j]('t'))^2+('y'[i]('t')-'y'[j]('t'))^2 +('z'[i]('t')-'z'[j]('t'))^2+(Q[i]('t')-Q[j]('t'))^2),0),j=1..n), nops(a[i])*diff(Q[i]('t'),'t','t')=-sqrt(n)*con*add(Q[i]('t')-Q[j]('t'),j=a[i])-nops(a[i])^2/100*diff(Q[i]('t'),'t') +con*add(`if`(j<>i,(Q[i]('t')-Q[j]('t'))/((x[i]('t')-x[j]('t'))^2+(y[i]('t')-y[j]('t'))^2 +(z[i]('t')-z[j]('t'))^2+(Q[i]('t')-Q[j]('t'))^2),0),j=1..n); end do; des:={seq(d[i],i=1..n), seq('x'[i](0)=ps[i][1],i=1..n), seq('y'[i](0)=ps[i][2],i=1..n),seq(D('y'[i])(0)=0,i=1..n), seq(D('x'[i])(0)=0,i=1..n),seq('z'[i](0)=ps[i][3],i=1..n),seq(D('z'[i])(0)=0,i=1..n), seq(Q[i](0)=ps[i][3],i=1..n),seq(D(Q[i])(0)=0,i=1..n)}; fi; end if; end if; if OUTPUT='plotsol' then if _Env_Use_RKF=true then sol:=dsolve(des,'numeric','relerr'=0.0001,'abserr'=0.0001,'initstep'=0.1); else sol:=dsolve(des,'numeric','relerr'=0.0001,'abserr'=0.0001,'initstep'=0.1); #sol:=dsolve(des,'numeric','method'='classical'['rk4'],'stepsize'=stps); end if; SOL:=sol(T); elif OUTPUT='anim' then if _Env_Use_RKF=true then sol:=dsolve(des,'numeric','relerr'=0.0001,'abserr'=0.0001,'initstep'=0.1,'output'=array([0,seq((T/N)*i,i=1..N)])); else sol:=dsolve(des,'numeric','method'='classical'['rk4'],'stepsize'=stps,'output'=array([0,seq(T/N*i,i=1..N)])); end if; SOL:=op(sol)[2,1]; end if; return(evalf(SOL,4)); end proc: #-------------------------------------------------------------------------------- #AE: if dim=2 then plotG should take advantage of GetPlotObjs for consistency plotG:=proc(G::GRAPHLN,SOL,n,dim,a,b,dir,N,OUTPUT,nodes,M,edgelab,verlab,verlabcol,vercol, edlabcol,arrcol,edgcol,redraw,shwts) #AE local i,points,edges,size,p,Text,weight,j,k,numw,ori, L, vpos, vl, e, directed, weighted, NN, edgfnt, wobjs, MM, r, u, v, mp, label, m2, m, m1, displayoptions;#,size1,edges1; if n<2 then error "number of vertices must be greater than 1" fi; #{----------------- L := GraphTheory:-GraphInfo:-LabelToInteger(G); #AE so that it will always execute. now we need to check that we are not trying to draw in 3d w/ the 2d vertex positions. vpos := GraphTheory:-GraphInfo:-GetVPos(G, VP_SPRING); vl := [seq(L[nodes[i]], i = 1..n)]; #}----------------- #AE this is absurd. if redraw is false we are still solving the system? if OUTPUT='plotsol' then #{----------------- if redraw = true then #}----------------- if dim=2 then #MAHDI # points:=evalf([seq([rhs(SOL[2*i]),rhs(SOL[2*n+2*i])],i=1..n)],4); points := [ seq( [SOL[i,1] , SOL[i,2]] , i = 1..n) ]; elif dim=3 then # points:=evalf([seq([rhs(SOL[2*i]),rhs(SOL[2*n+2*i]),rhs(SOL[4*n+2*i])],i=1..n)],4); points := [seq([ SOL[i,1] , SOL[i,2] , SOL[i,3] ], i = 1..n)] end if; #---- #{----------------- GraphTheory:-GraphInfo:-SetVPos(G, VP_SPRING, points, vl); #AE else points := [seq(vpos[vl[i]], i = 1..n)]; end if; #}----------------- size:=max(seq(max(seq(points[i][j],i=1..n))-min(seq(points[i][j],i=1..n)),j=1..dim)); size:=2*size/(n+5.0)^(1/4); edges :=Getedges(G, points,a,n,size,dim,dir,arrcol,edgcol,nodes); if dim=3 then Text:=Findetxtloc(points,b,size,n); end if; #AE looks like there was the intention of having weight labels. #AE garbage? if edgelab=1 then weight:=edgeweight(points,a,size,n,M,dim,Text,'numw'); end if; if dim=2 then #AE begin weight labels e := listn(G); directed := GraphTheory:-IsDirected(G); weighted:=GraphTheory:-IsWeighted(G); #AE if weighted and shwts=true and dim = 2 then #MBM NN:=GraphTheory:-NumberOfEdges(G); edgfnt := FONT(HELVETICA,`if`(NN>20,`if`(NN>45,9,10),11)); wobjs := table(); k := 0; MM := GraphTheory:-WeightMatrix(G); r:=.60; for u to n do for v in e[u] do if directed or v>u then mp:= [evalf(((1-r)*points[u,1]+r*points[v,1])),evalf(((1-r)*points[u,2]+r*points[v,2]))]: label:=MM[u,v]; if type(label,integer) then label:= sprintf("%d",label); else label:=sprintf("%0.3g",label); fi; m2:=evalf(points[u,1]-points[v,1]); if abs(m2) < .01*size then m:= infinity; mp:= mp + [1,0]*size*0.007; else m1:=points[u,2]-points[v,2]; if abs(m1) < 0.01*size then # put the label above the arrow if we have to if directed and u in e[v] and points[u][1]>points[v][1] then mp:= mp + [0,1]*size*0.015; else mp:= mp + [0,1]*size*0.005; fi; m:=0; else m:=evalf(m1/m2); mp := mp + size*`if`(m>0,0.009,0.002)*[abs(m1),abs(m2)*`if`(m>0,-1,1)]/sqrt(m1^2+m2^2); fi; fi; label:= TEXT(mp,label,ALIGNRIGHT,`if`(m>0,ALIGNBELOW,ALIGNABOVE),edgfnt); k:=k+1; wobjs[k]:=label; fi; od; od; wobjs := seq( wobjs[u], u=1..k ); else wobjs:=NULL fi; # AE end weight labels p:=PLOT( `if`(verlab=1,seq(TEXT(points[i],convert(nodes[i],string),verlabcol,FONT('HELVETICA',10)),i=1..n),NULL), #AE garbage? `if`(edgelab=1,seq(TEXT(weight[i][1],convert(weight[i][2],string),FONT('HELVETICA',10),edlabcol),i=1..numw),NULL), makebox(G, points,size,n,dim,vercol,nodes), POLYGONS(edges,STYLE(PATCHNOGRID)), wobjs, #AE added wobjs SCALING(CONSTRAINED), # MBM May/06 AXESSTYLE(NONE)); else ori:=findmax(points); p:=PLOT3D( `if`(verlab=1,seq(TEXT(Text[i],convert(nodes[i],string),verlabcol,FONT('HELVETICA',10)),i=1..n),NULL), `if`(edgelab=1,seq(TEXT(weight[i][1],convert(weight[i][2],string),FONT('HELVETICA',10),edlabcol),i=1..numw),NULL), makebox(G,points,size,n,dim,vercol,nodes), POLYGONS(edges,STYLE(PATCHNOGRID)), SCALING(CONSTRAINED), # MBM May/06 LIGHTMODEL(LIGHT_4), # so it looks like it's in 3D (thanks Ha) AXESSTYLE(NONE),ORIENTATION(op(ori))); end if; return p; elif OUTPUT='anim' then #{----------------- #MAHDI #AE if vpos = [] then #AE don't we want to execute this anyway? #vpos := [seq([seq(0, j = 1..GraphTheory:-NumberOfVertices(G))], i = 1..N+1)]; vpos := [seq([seq(0, j = 1..GraphTheory:-NumberOfVertices(G))], i = 1..N)]; #vpos := Array(1..N,1..GraphTheory:-NumberOfVertices); #AE we should use an array b/c we will assign to it. #AE end if; #----- #}----------------- #MAHDI # for i from 1 to N+1 do for i from 1 to N do #---- #{----------------- # if redraw = true then # this is wrong, it has to always be executed for an animation if true then #}----------------- #MAHDI if dim=2 then #points:=evalf([seq([SOL[i,2*k],SOL[i,2*n+2*k]],k=1..n)],4); points := [seq([ SOL[i,k,1] , SOL[i,k,2] ], k = 1..n)] elif dim=3 then #points:=evalf([seq([SOL[i,2*k],SOL[i,2*n+2*k],SOL[i,4*n+2*k]],k=1..n)],4); points := [seq([SOL[i,k,1] , SOL[i,k,2],SOL[i,k,3]], k = 1..n)] end if; #----- #{----------------- for j to nops(points) do vpos[i][vl[j]] := points[j]; #AE we are assigning to a huge list here!!! end do; else points := [seq(vpos[i][vl[j]], j = 1..nops(vl))]; end if; #}----------------- size:=max(seq(max(seq(points[k][dim],k=1..n))-min(seq(points[k][dim],k=1..n)),j=1..dim)); edges :=Getedges(G, points,a,n,size,dim,dir,arrcol,edgcol,nodes); if dim=3 then Text:=Findetxtloc(points,b,size,n); end if; #AE garbage code? if edgelab=1 then weight:=edgeweight(points,a,size,n,M,dim,Text,'numw'); end if; if dim=2 then p[i]:=PLOT( `if`(verlab=1,seq(TEXT(points[k],convert(nodes[k],string),verlabcol,FONT('HELVETICA',10)),k=1..n),NULL), #AE garbage code? `if`(edgelab=1,seq(TEXT(weight[k][1],convert(weight[k][2],string),FONT('HELVETICA',10),edlabcol),k=1..numw),NULL), makebox(G,points,size,n,dim,vercol,nodes), POLYGONS(edges,STYLE(PATCHNOGRID)), AXESSTYLE(NONE),SCALING(CONSTRAINED)); displayoptions := NULL; else p[i]:=PLOT3D( `if`(verlab=1,seq(TEXT(Text[k],convert(nodes[k],string),verlabcol,FONT('HELVETICA',10)),k=1..n),NULL), makebox(G,points,size,n,dim,vercol,nodes), POLYGONS(edges,STYLE(PATCHNOGRID)), AXESSTYLE(NONE),SCALING(CONSTRAINED)); displayoptions := 'lightmodel=light4'; # thanks Ha end if; end do; #MAHDI # return plots['display'](seq(p[i],i=1..N+1),'insequence'=true); return plots['display']([seq(p[i],i=1..N)],'insequence'=true,displayoptions); #----- end if; end proc: #-------------------------------------------------------------------------------- makebox:=proc(G::GRAPHLN,p,size,n,dim,vercol,nodes) local dis,PL,i,K,j, vcolor, L; #{----------------- L := GraphTheory:-GraphInfo:-LabelToInteger(G); #}----------------- vcolor := GraphTheory:-GetVColor(G); # MBM: fixed ref to GraphTheory if dim=2 then dis:=evalf(0.030*size); PL:=POLYGONS(seq([[p[i][1]+dis,p[i][2]-dis],[p[i][1]+dis,p[i][2]+dis],[p[i][1]-dis,p[i][2]+dis],[p[i][1]-dis,p[i][2]-dis]],i=1..n), COLOR(RGB, seq(seq(op(vcolor[L[nodes[i]]])[j],j=2..4),i=1..n)), #`if`(type(vercol,table),seq(seq(op(vercol[i])[j],j=2..4),i=1..n), #seq(seq(op(vercol)[j],j=2..4),i=1..n))), STYLE(PATCHNOGRID)); else dis:=evalf(0.036*size); PL:=POLYGONS(seq([[p[i][1]+dis,p[i][2]-dis,p[i][3]+dis],[p[i][1]+dis,p[i][2]-dis,p[i][3]-dis], [p[i][1]+dis,p[i][2]+dis,p[i][3]-dis],[p[i][1]+dis,p[i][2]+dis,p[i][3]+dis]],i=1..n), seq([[p[i][1]+dis,p[i][2]+dis,p[i][3]+dis],[p[i][1]-dis,p[i][2]+dis,p[i][3]+dis], [p[i][1]-dis,p[i][2]+dis,p[i][3]-dis],[p[i][1]+dis,p[i][2]+dis,p[i][3]-dis]],i=1..n), seq([[p[i][1]+dis,p[i][2]-dis,p[i][3]+dis],[p[i][1]+dis,p[i][2]+dis,p[i][3]+dis], [p[i][1]-dis,p[i][2]+dis,p[i][3]+dis],[p[i][1]-dis,p[i][2]-dis,p[i][3]+dis]],i=1..n), seq([[p[i][1]-dis,p[i][2]-dis,p[i][3]+dis],[p[i][1]-dis,p[i][2]+dis,p[i][3]+dis], [p[i][1]-dis,p[i][2]+dis,p[i][3]-dis],[p[i][1]-dis,p[i][2]-dis,p[i][3]-dis]],i=1..n), seq([[p[i][1]+dis,p[i][2]-dis,p[i][3]-dis],[p[i][1]+dis,p[i][2]+dis,p[i][3]-dis], [p[i][1]-dis,p[i][2]+dis,p[i][3]-dis],[p[i][1]-dis,p[i][2]-dis,p[i][3]-dis]],i=1..n), seq([[p[i][1]+dis,p[i][2]-dis,p[i][3]+dis],[p[i][1]-dis,p[i][2]-dis,p[i][3]+dis], [p[i][1]-dis,p[i][2]-dis,p[i][3]-dis],[p[i][1]+dis,p[i][2]-dis,p[i][3]-dis]],i=1..n), COLOR(RGB, #`if`(type(vercol,table),seq(seq(seq(op(vercol[i])[j],j=2..4),i=1..n),K=1..6), seq(seq(seq(op(vcolor[L[nodes[i]]])[j],j=2..4),i=1..n),K=1..6)), #), STYLE(PATCHNOGRID)); end if; return PL; end proc: #-------------------------------------------------------------------------------- Getedges:=proc(G::GRAPHLN,points2,a,n,size,dim,dir,arrcol,edgcol,nodes) local ED,xp1,xp2,yp1,yp2,i,j,zp1,zp2,clr,k,v,ecolor, L; #{----------------- L := GraphTheory:-GraphInfo:-LabelToInteger(G); #}----------------- ED:=[]; clr:=[]; ecolor := GraphTheory:-GraphInfo:-GetEdgesColor(G); if dim=2 then for i from 1 to n do for j in a[i] do xp1:=points2[i][1]; xp2:=points2[j][1]; yp1:=points2[i][2]; yp2:=points2[j][2]; if dir=0 then ED:=[op(ED),evalf(vectorplot2d(xp1,yp1,xp2,yp2,size,"nodir"),4)]; clr := [op(clr), seq(op(ecolor[L[nodes[i]], L[nodes[j]]])[k], k=2..4)]; #clr:=[op(clr),`if`(not type(edgcol,table),seq(op(edgcol)[k],k=2..4),seq(op(edgcol[i,j])[k],k=2..4))]; else if member(i,a[j]) and ii then if assigned(vtbl[i,j]) then vrs := op(vtbl[i,j]); elif assigned(vtbl[j,i]) then vrs := [-vtbl[j,i][1],-vtbl[j,i][2],vtbl[j,i][3]]; else vrs := [tmp[ntmp+1],tmp[ntmp+2],tmp[ntmp+3]]; ntmp := ntmp+3: vtbl[i,j] := vrs; xj,yj := 2*(j-1)+1,2*nops(a)+2*(j-1)+1; nstmt := nstmt+1; stmt[nstmt] := vrs[1]=Y[x]-Y[xj], vrs[2]=Y[y]-Y[yj], vrs[3]=1/(vrs[1]^2+vrs[2]^2); end if; if member(j,a[i]) then acc1x := acc1x+vrs[1]; acc1y := acc1y+vrs[2]; end if; acc2x := acc2x+vrs[1]*vrs[3]; acc2y := acc2y+vrs[2]*vrs[3]; end if; end do; nstmt := nstmt+1; stmt[nstmt] := YP[xp]=kd*(-100*Y[xp]-r*acc1x+acc2x), YP[yp]=kd*(-100*Y[yp]-r*acc1y+acc2y): end do: stmt := codegen['makeproc']([tmp=wazzup(1..ntmp),seq(stmt[i],i=1..nstmt)], 'parameters'=[N,t,Y,YP],'locals'=[kd,r,tmp]): stmt := subs('wazzup'=hfarray,eval(stmt)); vrs := map(op,[seq([x[i](t),diff(x[i](t),t)],i=1..nops(a)), seq([y[i](t),diff(y[i](t),t)],i=1..nops(a))]); return eval(stmt),vrs; end proc: #--------------------------------------------------------------------------------- BuildProc3 := proc(a,k) local locs,tmp,ntmp,stmt,nstmt,vtbl,vrs,acc1x,acc1y,acc1z,acc2x,acc2y,acc2z,wazzup, x,y,z,xp,yp,zp,xj,yj,zj,N,t,Y,YP,i,j,kd,r; ntmp := 0: # Initial spring setting stmt := table([1=(r=5.0)]): nstmt := 1: vtbl := table(): for i to nops(a) do x,xp := 2*(i-1)+1,2*(i-1)+2; y,yp := 2*nops(a)+2*(i-1)+1,2*nops(a)+2*(i-1)+2; z,zp := 4*nops(a)+2*(i-1)+1,4*nops(a)+2*(i-1)+2; nstmt := nstmt+1; stmt[nstmt] := kd=k/nops(a[i]),YP[x]=Y[xp],YP[y]=Y[yp],YP[z]=Y[zp]; acc1x,acc1y,acc1z := 0,0,0; acc2x,acc2y,acc2z := 0,0,0; for j to nops(a) do if j<>i then if assigned(vtbl[i,j]) then vrs := op(vtbl[i,j]); elif assigned(vtbl[j,i]) then vrs := [-vtbl[j,i][1],-vtbl[j,i][2],-vtbl[j,i][3],vtbl[j,i][4]]; else vrs := [tmp[ntmp+1],tmp[ntmp+2],tmp[ntmp+3],tmp[ntmp+4]]; ntmp := ntmp+4: vtbl[i,j] := vrs; xj,yj,zj := 2*(j-1)+1,2*nops(a)+2*(j-1)+1,4*nops(a)+2*(j-1)+1; nstmt := nstmt+1; stmt[nstmt] := vrs[1]=Y[x]-Y[xj],vrs[2]=Y[y]-Y[yj], vrs[3]=Y[z]-Y[zj], vrs[4]=1/(vrs[1]^2+vrs[2]^2+vrs[3]^2); end if; if member(j,a[i]) then acc1x := acc1x+vrs[1]; acc1y := acc1y+vrs[2]; acc1z := acc1z+vrs[3]; end if; acc2x := acc2x+vrs[1]*vrs[4]; acc2y := acc2y+vrs[2]*vrs[4]; acc2z := acc2z+vrs[3]*vrs[4]; end if; end do; nstmt := nstmt+1; stmt[nstmt] := YP[xp]=kd*(-100*Y[xp]-r*acc1x+acc2x), YP[yp]=kd*(-100*Y[yp]-r*acc1y+acc2y), YP[zp]=kd*(-100*Y[zp]-r*acc1z+acc2z): end do: stmt := codegen['makeproc']([tmp=wazzup(1..ntmp),seq(stmt[i],i=1..nstmt)], 'parameters'=[N,t,Y,YP],'locals'=[kd,r,tmp]): stmt := subs('wazzup'=hfarray,eval(stmt)); vrs := map(op,[seq([x[i](t),diff(x[i](t),t)],i=1..nops(a)), seq([y[i](t),diff(y[i](t),t)],i=1..nops(a)), seq([z[i](t),diff(z[i](t),t)],i=1..nops(a))]); return eval(stmt),vrs; end proc: #--------------------------------------------------------------------------------- projection:=proc(points,orientation) local i,c,d,a,b,M,image,proj,CA,CB,SA,SB,A,B,C; c:=orientation[1]: d:=orientation[2]: a:=evalf(c*Pi/180):b:=evalf(d*Pi/180): SA:=sin(a); SB:=sin(b); CA:=cos(a); CB:=cos(b); A:=CA*SB; B:=SA*SB; C:=CB; M:=convert(points,Matrix): #Rot:=Matrix([[CB, CA*SB, SA*SB], [-SB, CA*CB, SA*CB], [0, -SA, CA]]): #project:=Matrix([[B^2+C^2,-A*B,-A*C],[-B*A,A^2+C^2,-B*C],[-C*A,-C*B,A^2+B^2]]): #proj:=LinearAlgebra:-MatrixMatrixMultiply(M,project): #proj:=LinearAlgebra:-MatrixMatrixMultiply(proj,Rot): image:=Matrix([[(SA^2*SB^2+CB^2)*CB+CA*SB^3*SA, (SA^2*SB^2+CB^2)*CA*SB-CA^2*SB^2*SA*CB+CA*SB*CB*SA, (SA^2*SB^2+CB^2)*SA*SB-CA*SB^2*SA^2*CB-CA^2*SB*CB], [-CA*SB^2*SA*CB-(CA^2*SB^2+CB^2)*SB, -CA^2*SB^3*SA+(CA^2*SB^2+CB^2)*CA*CB+SA^2*SB*CB, -CA*SB^3*SA^2+(CA^2*SB^2+CB^2)*SA*CB-CA*SB*CB*SA], [-CA*SB*CB^2+SA*SB^2*CB, -CA^2*SB^2*CB-SA*SB*CB^2*CA-(CA^2*SB^2+SA^2*SB^2)*SA, -CA*SB^2*SA*CB-SA^2*SB*CB^2+(CA^2*SB^2+SA^2*SB^2)*CA]]); proj:=LinearAlgebra:-MatrixMatrixMultiply(M,image); return [evalf(seq([proj[i,1],proj[i,2]],i=1..nops(points)),4)]; end proc: #-------------------------------------------------------------------------------- findmax:=proc(points) local ps,MAX,maxpos,temp,i,j; ps:=projection(points,[0,0]): MAX:=cnvxhull( ps, 'output'='area' ): maxpos:=[0,0]: for i from -90 to 90 by 45 do for j from -90 to 90 by 45 do ps:=projection(points,[i,j]): temp:=cnvxhull( ps, 'output'=area ): if MAX simplex[convexhull]( ps, output=area ) outputs the area of theconvex hull #--> simplex[convexhull]( ps, output=plot ) outputs a plot of the convexhull and the points #--> simplex[convexhull]( ps, output=hull ) outputs the list of pointson the convex hull # #macro( yellow=COLOR(RGB,1,1,0) ): cnvxhull := proc( ps::{list,set}([numeric,numeric]) ) local i, inds, is, j, n, sols, T, x, y, hh, M, l, ps1, k , H, P1, P2, P3, det, A; if not type(ps,set) and not type(ps,list) then error "invalid arguments" fi; if nargs=2 and not member( args[2], {'output=area', 'output=hull', 'output=plot'} ) then error "invalid output option: %1", args[2]; fi; n := nops(ps); if n < 3 then H := ps; if nargs=2 and args[2]='output=area' then return 0; fi; if nargs=2 and args[2]='output=plot' then return PLOT( POINTS(op(ps)), CURVES(H,yellow) ); fi; return H; fi; for l to n do M[l] := [op(evalf(ps[l])), op(ps[l])]; od; ps1 := convert(M, list); # Find one point on the hull (leftmost), and start there x := min( op( map( (p -> op(1,p)), ps1 ))); y := max( op( map( proc(p,x) if op(1,p)=x then op(2,p) fi end, ps1, x ))); # Form table T of tangents map( proc(p,x,y,T) local t; if p[2]<>y or p[1]<>x then if x=p[1] then t := (p[2]-y)*1000000 else t := (p[2]-y)/(p[1]-x) fi; if not type(t,numeric) then error "non-numeric data" fi; if assigned(T[t]) then if p[1] > T[t][1] then T[t] := p fi else T[t] := p fi fi end, ps1, x, y, T); # Sort all tangents inds := sort( map( op, [indices(T)]) ); # First two points in the hull are: for hh in ps1 do if hh[1] = x and hh[2] = y then sols[1] := [hh[3], hh[4]] fi; od; sols[2] := [T[inds[1]][3],T[inds[1]][4]]; if nops(inds) < 2 then H := [sols[1], sols[2]]; if nargs=2 and args[2]='output=area' then return 0; fi; if nargs=2 and args[2]='output=plot' then return PLOT( POINTS(op(ps)), CURVES(H,yellow) ); fi; return H; fi; sols[3] := T[inds[2]]; is := 3; # Main loop to find all other points for i from 3 to nops(inds) do is := is+1; sols[is] := T[inds[i]]; for is from is by -1 to 4 while (sols[is][2]-sols[is-2][2])*(sols[is-1][1]-sols[is-2][1]) - (sols[is][1]-sols[is-2][1])*(sols[is-1][2]-sols[is-2][2]) <= 0 do sols[is-1] := sols[is] od; od; for k from 3 to is do sols[k] := [sols[k][3],sols[k][4]] od; H := [ sols[j] $ j=1..is ]; # MBM: generate the required output if nargs = 2 and args[2] = 'output=area' then A := 0; P1 := H[1]; det := proc(u,v) u[1]*v[2]-u[2]*v[1] end: # area of the triangle for i from 2 to nops(H)-1 do P2,P3 := H[i],H[i+1]; A := A + abs(det(P2-P1,P3-P2)/2); # area of the triangle P1, P2, P3 od: H := A; fi; if nargs = 2 and args[2] = 'output=plot' then H := PLOT( POINTS(op(ps)), POLYGONS(H,STYLE(PATCHNOGRID),yellow) ); fi; H; end: end module: #savelib('GraphDrawing');