Overall issue: I am finding the intersection of two planes. I have the values of both planes at ungridded points, so I subtract their values (call the difference z). Then I have x,y,z values and I find the contour for z=0. I then plan to take this contour and map it back to my original surfaces. There are several ways to accomplish that last step.
Specific issue: This is the code I have so far:
# random data to simulate z
f(x, y) = sin(x + 10sin(x)) + cos(y)
xx = rand(10000);
yy = rand(10000);
d = hcat(hcat(xx,yy),f.(xx,yy));
# create Delauny triangularization
using GMT
tri = GMT.triangulate(d) #the output are index numbers of the original points based on line
# create contour at z=0
cont = GMT.contour(d,index=tri, dump=true, contour=[0]) #this gives x,y,z values
Given a specific point in cont, say cont[1][1,1:2], is there an easy way to know which triangle/edge it is in/on?
I believe that I can code a search through all the triangles, but I feel like there should be an easier way. I can also take the x,y points in the contour and plug them back into my functions (not necessarily what i want to do though).
Most of the documentation I’m reading has to do with plotting or output of tri and cont, but I don’t know how to effectively use the structure of tri in anything but plotting.
Do you want to plot a trisurf, i.e. a triangulated surface of equation z=f(x,y), (x,y)\in [a,b]\times[c,d]?
If the answer is positive, then I don’t think that GMT provides a function to plot it. Instead you can use DelaunayTriangulations.jl to triangulate the rectangle of definition.
Namely, you should define a mesgrid over that rectangle
and triangulate the set of meshgrid points, (x,y), stored in a vector pts.
Then lifting the vertices of the Delaunay triangulation to the surface, i.e. to
3d points, (x,y, f(x,y)), and keeping the triangles resulted from the planar points triangulation, we get a surface triangulation that can be plotted as Plots.mesh3d:
Example:
using DelaunayTriangulation
using Plots, ColorSchemes
plotlyjs();
f(x, y) = sin(x) + cos(y)
m = 75
n = 75 #n can be different from m
xl = range(0, 2pi, n) # or xl= 2pi*rand(n)
yl=range(0, 2pi, m) # yl= 2pi*rand(m)
x = vec(ones(m)*xl')
y = vec(yl *ones(n)')
z = f.(x, y)
pts = [[xe, ye] for (xe, ye) in zip(x, y)]
tri = triangulate(pts)
cm = stack(tri.triangles, dims=1) #cm is a Matrix of size (N, 3)
plt = plot(seriestype=:mesh3d, x, y, z,
connections=tuple([(c .-1) for c in eachcol(cm)]...), #mesh3d works with tri inds starting from 0
legend=false, intensity=z,
seriescolor = reverse(cgrad(ColorSchemes.matter.colors)), size=(400, 400))
