diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 3d34edd..b730cd5 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -29,7 +29,7 @@ jobs: python -m pip install --upgrade pip pip install numpy matplotlib - name: Checkout repository - uses: actions/checkout@v1 + uses: actions/checkout@main with: repository: 'NaokiHori/EllipsesInFlows' ref: 'main' @@ -67,7 +67,7 @@ jobs: echo "Date :" $(date) >> artifacts/ci.txt echo "Hash :" $(git rev-parse HEAD) >> artifacts/ci.txt - name: Upload artifacts - uses: actions/upload-artifact@v3 + uses: actions/upload-artifact@main with: name: LateralMigration${{ matrix.resolution }} path: artifacts @@ -88,7 +88,7 @@ jobs: python -m pip install --upgrade pip pip install numpy matplotlib - name: Checkout repository - uses: actions/checkout@v1 + uses: actions/checkout@main with: repository: 'NaokiHori/EllipsesInFlows' ref: 'main' @@ -123,7 +123,7 @@ jobs: echo "Date :" $(date) >> artifacts/ci.txt echo "Hash :" $(git rev-parse HEAD) >> artifacts/ci.txt - name: Upload artifacts - uses: actions/upload-artifact@v3 + uses: actions/upload-artifact@main with: name: SegreSilberberg${{ matrix.resolution }} path: artifacts @@ -144,7 +144,7 @@ jobs: python -m pip install --upgrade pip pip install numpy matplotlib - name: Checkout repository - uses: actions/checkout@v1 + uses: actions/checkout@main with: repository: 'NaokiHori/EllipsesInFlows' ref: 'main' @@ -182,7 +182,7 @@ jobs: echo "Date :" $(date) >> artifacts/ci.txt echo "Hash :" $(git rev-parse HEAD) >> artifacts/ci.txt - name: Upload artifacts - uses: actions/upload-artifact@v3 + uses: actions/upload-artifact@main with: name: RotatingEllipse${{ matrix.resolution }} path: artifacts @@ -205,60 +205,65 @@ jobs: - name: Install dependencies run: | sudo apt-get update && \ - sudo apt-get -y install gnuplot ghostscript imagemagick + sudo apt-get -y install make gnuplot ghostscript imagemagick - name: Edit ghostscript config run: | sudo sed -i 's/rights="none"/rights="read|write"/g' /etc/ImageMagick-6/policy.xml - name: Checkout repository - uses: actions/checkout@v1 + uses: actions/checkout@main with: repository: 'NaokiHori/EllipsesInFlows' ref: 'main' - name: Download LateralMigration16 - uses: actions/download-artifact@v3 + uses: actions/download-artifact@main with: name: LateralMigration16 path: docs/source/examples/case1/data/16 - name: Download LateralMigration32 - uses: actions/download-artifact@v3 + uses: actions/download-artifact@main with: name: LateralMigration32 path: docs/source/examples/case1/data/32 - name: Download LateralMigration48 - uses: actions/download-artifact@v3 + uses: actions/download-artifact@main with: name: LateralMigration48 path: docs/source/examples/case1/data/48 - name: Download SegreSilberberg032 - uses: actions/download-artifact@v3 + uses: actions/download-artifact@main with: name: SegreSilberberg32 path: docs/source/examples/case2/data/32 - name: Download SegreSilberberg064 - uses: actions/download-artifact@v3 + uses: actions/download-artifact@main with: name: SegreSilberberg64 path: docs/source/examples/case2/data/64 - name: Download SegreSilberberg96 - uses: actions/download-artifact@v3 + uses: actions/download-artifact@main with: name: SegreSilberberg96 path: docs/source/examples/case2/data/96 - name: Download RotatingEllipse48 - uses: actions/download-artifact@v3 + uses: actions/download-artifact@main with: name: RotatingEllipse48 path: docs/source/examples/case3/data/48 - name: Download RotatingEllipse64 - uses: actions/download-artifact@v3 + uses: actions/download-artifact@main with: name: RotatingEllipse64 path: docs/source/examples/case3/data/64 - name: Download RotatingEllipse96 - uses: actions/download-artifact@v3 + uses: actions/download-artifact@main with: name: RotatingEllipse96 path: docs/source/examples/case3/data/96 + - name: Run collision scripts + run: | + cd docs/source/collision/data + make all + ./a.out - name: Create eps and tex files run: | cd docs @@ -298,9 +303,9 @@ jobs: sphinxdoc/sphinx:5.0.1 \ /bin/bash sphinx.sh - name: Setup GitHub Pages - uses: actions/configure-pages@v1 + uses: actions/configure-pages@main - name: Upload HTML - uses: actions/upload-pages-artifact@v1 + uses: actions/upload-pages-artifact@main with: path: docs/build/html - name: Deploy to GitHub Pages diff --git a/docs/source/_static/.gitignore b/docs/source/_static/.gitignore deleted file mode 100644 index e69de29..0000000 diff --git a/docs/source/collision/data/Makefile b/docs/source/collision/data/Makefile new file mode 100644 index 0000000..d8c8dc9 --- /dev/null +++ b/docs/source/collision/data/Makefile @@ -0,0 +1,35 @@ +CC := cc +CFLAGS := -O3 -std=c99 -Wall -Wextra +DEPEND := -MMD +LIBS := -lm +INCLUDES := -Iinclude +SRCSDIR := src +OBJSDIR := obj +SRCS := $(foreach dir, $(shell find $(SRCSDIR) -type d), $(wildcard $(dir)/*.c)) +OBJS := $(addprefix $(OBJSDIR)/, $(subst $(SRCSDIR)/,,$(SRCS:.c=.o))) +DEPS := $(addprefix $(OBJSDIR)/, $(subst $(SRCSDIR)/,,$(SRCS:.c=.d))) +TARGET := a.out + +help: + @echo "all : create \"$(TARGET)\"" + @echo "clean : remove \"$(TARGET)\" and object files \"$(OBJSDIR)/*.o\"" + @echo "help : show this help message" + +all: $(TARGET) + +$(TARGET): $(OBJS) + $(CC) $(CFLAGS) $(DEPEND) -o $@ $^ $(LIBS) + +$(OBJSDIR)/%.o: $(SRCSDIR)/%.c + @if [ ! -e `dirname $@` ]; then \ + mkdir -p `dirname $@`; \ + fi + $(CC) $(CFLAGS) $(DEPEND) $(INCLUDES) -c $< -o $@ + +clean: + $(RM) -r $(OBJSDIR) $(TARGET) + +-include $(DEPS) + +.PHONY : help all clean + diff --git a/docs/source/collision/data/example.py b/docs/source/collision/data/example.py new file mode 100644 index 0000000..d0e8a0b --- /dev/null +++ b/docs/source/collision/data/example.py @@ -0,0 +1,190 @@ +import numpy as np +from matplotlib import pyplot +from matplotlib import patches + + +def shift_and_rotate(disp, angle, p0): + # shift a point p0, then rotate it + # return the resulting point p1 + p1 = dict() + x = p0["x"]+disp["x"] + y = p0["y"]+disp["y"] + p1["x"] = np.cos(angle)*x - np.sin(angle)*y + p1["y"] = np.sin(angle)*x + np.cos(angle)*y + return p1 + + +def rotate_and_shift(disp, angle, p0): + # rotate a point (x0, y0), then shift it + # return the resulting point p1 + p1 = dict() + x = p0["x"] + y = p0["y"] + p1["x"] = disp["x"] + np.cos(angle)*x - np.sin(angle)*y + p1["y"] = disp["y"] + np.sin(angle)*x + np.cos(angle)*y + return p1 + + +def compute_evolute(e, t): + # compute center of evolute and radius of fitted circle + evolute = dict() + a = e["a"] + b = e["b"] + evolute["x"] = a*(1.-np.power(b/a, 2.))*np.power(np.cos(t), 3.) + evolute["y"] = b*(1.-np.power(a/b, 2.))*np.power(np.sin(t), 3.) + num = np.power(np.power(a*np.sin(t), 2.)+np.power(b*np.cos(t), 2.), 1.5) + den = a*b + evolute["r"] = num/den + return evolute + + +def find_normal_t(e, p): + # compute parameter t, with which a point (a*cos(t), b*sin(t)) + # and the target point p gives a normal vector to the given ellipse + # ref: https://blog.chatfield.io/simple-method-for-distance-to-ellipse/ + t = 0.25*np.pi + while True: + xe = e["a"]*np.cos(t) + ye = e["b"]*np.sin(t) + evolute = compute_evolute(e, t) + dxe = xe-evolute["x"] + dye = ye-evolute["y"] + dxp = abs(p["x"])-evolute["x"] + dyp = abs(p["y"])-evolute["y"] + norme = np.hypot(dxe, dye) + normp = np.hypot(dxp, dyp) + dc = norme*np.arcsin((dxe*dyp-dye*dxp)/(norme*normp)) + dt = dc/np.sqrt(e["a"]**2.+e["b"]**2.-xe**2.-ye**2.) + t += dt + # limit in the first quadrant + t = min(t, 0.5*np.pi) + t = max(t, 0.) + if abs(dt) < 1.e-8: + break + # apply result to the other quadrants + if p["x"] < 0.: + t = -t+np.pi + if p["y"] < 0.: + t = -t + return t + + +def fit_circles(e0, e1): + # core function, fitting two circles + # read the documentation + for n in range(10): + if n == 0: + # initialise temporary evolute + p0 = {"x": e0["x"], "y": e0["y"]} + p1 = {"x": e1["x"], "y": e1["y"]} + # forward coordinate transform + disp = {"x": -e0["x"], "y": -e0["y"]} + angle = -e0["theta"] + p1_ = shift_and_rotate(disp, angle, p1) + # forward coordinate transform + disp = {"x": -e1["x"], "y": -e1["y"]} + angle = -e1["theta"] + p0_ = shift_and_rotate(disp, angle, p0) + # find optimum parameter t for each ellipse + e0_t = find_normal_t(e0, p1_) + e1_t = find_normal_t(e1, p0_) + # compute evolutes and their radii + evolute0_ = compute_evolute(e0, e0_t) + evolute1_ = compute_evolute(e1, e1_t) + # backward coordinate transform + disp = {"x": +e0["x"], "y": +e0["y"]} + angle = +e0["theta"] + p0 = rotate_and_shift(disp, angle, evolute0_) + # backward coordinate transform + disp = {"x": +e1["x"], "y": +e1["y"]} + angle = +e1["theta"] + p1 = rotate_and_shift(disp, angle, evolute1_) + evolute0 = {"x": p0["x"], "y": p0["y"], "r": evolute0_["r"]} + evolute1 = {"x": p1["x"], "y": p1["y"], "r": evolute1_["r"]} + return evolute0, evolute1 + + +def plot_result(e0, e1, c0, c1): + # wrapper of matplotlib.patches to visualise + # resulting ellipses and fitted circles + def get_patch_of_ellipse(e): + x = e["x"] + y = e["y"] + width = 2.*e["a"] + height = 2.*e["b"] + angle = 180./np.pi*e["theta"] + keywords = { + "xy": (x, y), + "width": width, + "height": height, + "angle": angle, + "fc": "none", + "ec": "#FF0000", + } + return patches.Ellipse(**keywords) + + def get_patch_of_circle(c): + x = c["x"] + y = c["y"] + radius = c["r"] + keywords = { + "xy": (x, y), + "radius": radius, + "fc": "none", + "ec": "#0000FF", + } + return patches.Circle(**keywords) + + def compute_bounds(es, cs): + # determine plot ranges + # i.e., arguments of set_[xy]lim + xmin = +np.inf + xmax = -np.inf + ymin = +np.inf + ymax = -np.inf + for e in es: + xmin = min(xmin, e["x"]-max(e["a"], e["b"])) + xmax = max(xmax, e["x"]+max(e["a"], e["b"])) + ymin = min(ymin, e["y"]-max(e["a"], e["b"])) + ymax = max(ymax, e["y"]+max(e["a"], e["b"])) + for c in cs: + xmin = min(xmin, c["x"]-c["r"]) + xmax = max(xmax, c["x"]+c["r"]) + ymin = min(ymin, c["y"]-c["r"]) + ymax = max(ymax, c["y"]+c["r"]) + margin = 0.25 + xmin -= margin + xmax += margin + ymin -= margin + ymax += margin + return (xmin, xmax), (ymin, ymax) + # prepare matplotlib objects + fig = pyplot.figure() + ax = fig.add_subplot(111) + # draw ellipses and circles + ax.add_patch(get_patch_of_ellipse(e0)) + ax.add_patch(get_patch_of_ellipse(e1)) + ax.add_patch(get_patch_of_circle(c0)) + ax.add_patch(get_patch_of_circle(c1)) + # set auxiliary parameters + xlim, ylim = compute_bounds([e0, e1], [c0, c1]) + keywords = { + "xlim": xlim, + "ylim": ylim, + "aspect": "equal", + } + ax.set(**keywords) + # visualise + pyplot.show() + # clean-up + pyplot.close() + + +if __name__ == "__main__": + # initialise ellipses + e0 = {"a": 2.0, "b": 1.5, "x": 0.5, "y": 0.5, "theta": 0.2} + e1 = {"a": 1.5, "b": 1.0, "x": 2.0, "y": 2.5, "theta": 2.0} + # fit circles + c0, c1 = fit_circles(e0, e1) + # plot + plot_result(e0, e1, c0, c1) diff --git a/docs/source/collision/data/fit-circle.gp b/docs/source/collision/data/fit-circle.gp new file mode 100644 index 0000000..cf3e230 --- /dev/null +++ b/docs/source/collision/data/fit-circle.gp @@ -0,0 +1,75 @@ +reset + +load 'fit-circle.dat' + +min(x, y) = x > y ? y : x +max(x, y) = x < y ? y : x + +xe(a, t) = a*cos(t) +ye(b, t) = b*sin(t) + +x_t0 = atan2(-e0_b*sin(e0_theta), e0_a*cos(e0_theta)) +x_t1 = x_t0+pi +y_t0 = atan2(+e0_b*cos(e0_theta), e0_a*sin(e0_theta)) +y_t1 = y_t0+pi +x0 = e0_x+xe(e0_a, x_t0)*cos(e0_theta)-ye(e0_b, x_t0)*sin(e0_theta) +x1 = e0_x+xe(e0_a, x_t1)*cos(e0_theta)-ye(e0_b, x_t1)*sin(e0_theta) +y0 = e0_y+xe(e0_a, y_t0)*sin(e0_theta)+ye(e0_b, y_t0)*cos(e0_theta) +y1 = e0_y+xe(e0_a, y_t1)*sin(e0_theta)+ye(e0_b, y_t1)*cos(e0_theta) +xmin = min(e0_xp, min(x0, x1)) +xmax = max(e0_xp, max(x0, x1)) +ymin = min(e0_yp, min(y0, y1)) +ymax = max(e0_yp, max(y0, y1)) + +mrgn = 0.25 + +xmin = xmin-mrgn +xmax = xmax+mrgn +ymin = ymin-mrgn +ymax = ymax+mrgn + +lx = xmax-xmin +ly = ymax-ymin + +set terminal epslatex standalone color size lx,ly font ',17.28' +set output 'fit-circle.tex' + +unset border + +set lmargin 0 +set rmargin 0 +set bmargin 0 +set tmargin 0 + +unset xlabel +unset ylabel + +set xrange [xmin:xmax] +set yrange [ymin:ymax] + +unset xtics +unset ytics + +set size ratio -1 + +set style line 1 lc rgb '#000000' lw 7 +set style arrow 1 heads ls 1 +set style arrow 2 nohead ls 1 dt 3 +asize = 0.5 + +set object ellipse at first e0_x, first e0_y size first 2.*e0_a, first 2.*e0_b angle 180./pi*e0_theta fillstyle empty fc rgb '#FF0000' lw 3 +e0x = e0_x +e0y = e0_y +e1x = e0_x+e0_a*cos(e0_theta) +e1y = e0_y+e0_a*sin(e0_theta) +e2x = e0_x-e0_b*sin(e0_theta) +e2y = e0_y+e0_b*cos(e0_theta) + +set object circle at first e0_xc, first e0_yc size first e0_r fillstyle empty fc rgb '#0000FF' lw 3 + +set arrow from first e0_xp, first e0_yp to first e0_xc, first e0_yc as 1 +set label '$\left( x_c, y_c \right)$' right at first e0_xc, first e0_yc +set label '$\left( x_p, y_p \right)$\ \ ' right at first e0_xp, first e0_yp + +plot \ + NaN notitle diff --git a/docs/source/collision/data/fit-circles.gp b/docs/source/collision/data/fit-circles.gp new file mode 100644 index 0000000..5721b9b --- /dev/null +++ b/docs/source/collision/data/fit-circles.gp @@ -0,0 +1,88 @@ + +do for [case = 0:1:1] { + + reset + + load sprintf('fit-circles-%d.dat', case) + + array e_xs[2] = [e0_x, e1_x] + array e_ys[2] = [e0_y, e1_y] + array e_as[2] = [e0_a, e1_a] + array e_bs[2] = [e0_b, e1_b] + array e_ts[2] = [e0_theta, e1_theta] + array e_xcs[2] = [e0_xc, e1_xc] + array e_ycs[2] = [e0_yc, e1_yc] + array e_rs[2] = [e0_r, e1_r] + + min(x, y) = x > y ? y : x + max(x, y) = x < y ? y : x + + xe(a, t) = a*cos(t) + ye(b, t) = b*sin(t) + + xmin = +1.e+16 + xmax = -1.e+16 + ymin = +1.e+16 + ymax = -1.e+16 + + do for [n = 1:2:1] { + x_t0 = atan2(-e_bs[n]*sin(e_ts[n]), e_as[n]*cos(e_ts[n])) + x_t1 = x_t0+pi + y_t0 = atan2(+e_bs[n]*cos(e_ts[n]), e_as[n]*sin(e_ts[n])) + y_t1 = y_t0+pi + x0 = e_xs[n]+xe(e_as[n], x_t0)*cos(e_ts[n])-ye(e_bs[n], x_t0)*sin(e_ts[n]) + x1 = e_xs[n]+xe(e_as[n], x_t1)*cos(e_ts[n])-ye(e_bs[n], x_t1)*sin(e_ts[n]) + y0 = e_ys[n]+xe(e_as[n], y_t0)*sin(e_ts[n])+ye(e_bs[n], y_t0)*cos(e_ts[n]) + y1 = e_ys[n]+xe(e_as[n], y_t1)*sin(e_ts[n])+ye(e_bs[n], y_t1)*cos(e_ts[n]) + x2 = e_xcs[n]-e_rs[n] + xmin = min(xmin, min(e_xcs[n]-e_rs[n], min(x0, x1))) + xmax = max(xmax, max(e_xcs[n]+e_rs[n], max(x0, x1))) + ymin = min(ymin, min(e_ycs[n]-e_rs[n], min(y0, y1))) + ymax = max(ymax, max(e_ycs[n]+e_rs[n], max(y0, y1))) + } + + mrgn = 0.25 + + xmin = xmin-mrgn + xmax = xmax+mrgn + ymin = ymin-mrgn + ymax = ymax+mrgn + + lx = xmax-xmin + ly = ymax-ymin + + set terminal epslatex standalone color size lx,ly font ',17.28' + set output sprintf('fit-circles-%d.tex', case) + + unset border + + set lmargin 0 + set rmargin 0 + set bmargin 0 + set tmargin 0 + + unset xlabel + unset ylabel + + set xrange [xmin:xmax] + set yrange [ymin:ymax] + + unset xtics + unset ytics + + set size ratio -1 + + set style line 1 lc rgb '#000000' lw 7 + set style arrow 1 heads ls 1 + set style arrow 2 nohead ls 1 dt 3 + asize = 0.5 + + set object ellipse at first e0_x, first e0_y size first 2.*e0_a, first 2.*e0_b angle 180./pi*e0_theta fillstyle empty fc rgb '#FF0000' lw 3 + set object circle at first e0_xc, first e0_yc size first e0_r fillstyle empty fc rgb '#0000FF' lw 3 + set object ellipse at first e1_x, first e1_y size first 2.*e1_a, first 2.*e1_b angle 180./pi*e1_theta fillstyle empty fc rgb '#FF0000' lw 3 + set object circle at first e1_xc, first e1_yc size first e1_r fillstyle empty fc rgb '#0000FF' lw 3 + + plot \ + NaN notitle +} + diff --git a/docs/source/collision/data/include/ellipse.h b/docs/source/collision/data/include/ellipse.h new file mode 100644 index 0000000..4f962d7 --- /dev/null +++ b/docs/source/collision/data/include/ellipse.h @@ -0,0 +1,9 @@ +#if !defined(ELLIPSE_H) +#define ELLIPSE_H + +extern double ellipse_compute_xc(const double a, const double b, const double t); +extern double ellipse_compute_yc(const double a, const double b, const double t); +extern double ellipse_compute_curvature(const double a, const double b, const double t); +extern double ellipse_find_normal_t(const double a, const double b, const double xp, const double yp); + +#endif // ELLIPSE_H diff --git a/docs/source/collision/data/include/example.h b/docs/source/collision/data/include/example.h new file mode 100644 index 0000000..6a0c760 --- /dev/null +++ b/docs/source/collision/data/include/example.h @@ -0,0 +1,7 @@ +#if !defined(EXAMPLE_H) +#define EXAMPLE_H + +extern int fit_circle(void); +extern int fit_circles(void); + +#endif // EXAMPLE_H diff --git a/docs/source/collision/data/include/others.h b/docs/source/collision/data/include/others.h new file mode 100644 index 0000000..060b69f --- /dev/null +++ b/docs/source/collision/data/include/others.h @@ -0,0 +1,12 @@ +#if !defined(OTHERS_H) +#define OTHERS_H + +#if !defined(M_PI) +#define M_PI 3.1415926535897932 +#endif + +extern int shift_and_rotate(const double dx, const double dy, const double dt, const double x0, const double y0, double *x1, double *y1); +extern int rotate_and_shift(const double dx, const double dy, const double dt, const double x0, const double y0, double *x1, double *y1); +extern double get_random_value(const double min, const double max); + +#endif // OTHERS_H diff --git a/docs/source/collision/data/problem-setup.gp b/docs/source/collision/data/problem-setup.gp new file mode 100644 index 0000000..7d56506 --- /dev/null +++ b/docs/source/collision/data/problem-setup.gp @@ -0,0 +1,90 @@ +reset + +array e_xs[2] = [0.5, 2.5] +array e_ys[2] = [0.5, 3.0] +array e_as[2] = [2.0, 1.5] +array e_bs[2] = [1.5, 1.0] +array e_ts[2] = [0.5, 1.2] + +min(x, y) = x > y ? y : x +max(x, y) = x < y ? y : x + +xe(a, t) = a*cos(t) +ye(b, t) = b*sin(t) + +xmin = +1.e+16 +xmax = -1.e+16 +ymin = +1.e+16 +ymax = -1.e+16 + +do for [n = 1:2:1] { + x_t0 = atan2(-e_bs[n]*sin(e_ts[n]), e_as[n]*cos(e_ts[n])) + x_t1 = x_t0+pi + y_t0 = atan2(+e_bs[n]*cos(e_ts[n]), e_as[n]*sin(e_ts[n])) + y_t1 = y_t0+pi + x0 = e_xs[n]+xe(e_as[n], x_t0)*cos(e_ts[n])-ye(e_bs[n], x_t0)*sin(e_ts[n]) + x1 = e_xs[n]+xe(e_as[n], x_t1)*cos(e_ts[n])-ye(e_bs[n], x_t1)*sin(e_ts[n]) + y0 = e_ys[n]+xe(e_as[n], y_t0)*sin(e_ts[n])+ye(e_bs[n], y_t0)*cos(e_ts[n]) + y1 = e_ys[n]+xe(e_as[n], y_t1)*sin(e_ts[n])+ye(e_bs[n], y_t1)*cos(e_ts[n]) + xmin = min(xmin, min(x0, x1)) + xmax = max(xmax, max(x0, x1)) + ymin = min(ymin, min(y0, y1)) + ymax = max(ymax, max(y0, y1)) +} + +mrgn = 0.25 + +xmin = xmin-mrgn +xmax = xmax+mrgn +ymin = ymin-mrgn +ymax = ymax+mrgn + +lx = xmax-xmin +ly = ymax-ymin + +set terminal epslatex standalone color size lx,ly font ',17.28' +set output 'problem-setup.tex' + +unset border + +set lmargin 0 +set rmargin 0 +set bmargin 0 +set tmargin 0 + +unset xlabel +unset ylabel + +set xrange [xmin:xmax] +set yrange [ymin:ymax] + +unset xtics +unset ytics + +set size ratio -1 + +set style line 1 lc rgb '#AAAAAA' lw 5 +set style arrow 1 head ls 1 +set style arrow 2 nohead ls 1 dt 3 +asize = 0.5 + +do for [n = 1 : 2 : 1] { + set object ellipse at first e_xs[n], first e_ys[n] size first 2.*e_as[n], first 2.*e_bs[n] angle 180./pi*e_ts[n] fillstyle empty fc rgb '#FF0000' lw 3 + e0x = e_xs[n] + e0y = e_ys[n] + e1x = e_xs[n]+e_as[n]*cos(e_ts[n]) + e1y = e_ys[n]+e_as[n]*sin(e_ts[n]) + e2x = e_xs[n]-e_bs[n]*sin(e_ts[n]) + e2y = e_ys[n]+e_bs[n]*cos(e_ts[n]) + set arrow from first e0x, first e0y to first e1x, first e1y as 1 + set arrow from first e0x, first e0y to first e2x, first e2y as 1 + set arrow from graph 0., first e0y to graph 1., first e0y as 2 + set object circle at first e0x, first e0y size first 0.5 arc [0.:180./pi*e_ts[n]] fillstyle empty fc rgb '#AAAAAA' lw 3 + set label sprintf('$a_%d$', n-1) center at first 0.5*(e0x+e1x), first 0.5*(e0y+e1y) front + set label sprintf('$b_%d$', n-1) center at first 0.5*(e0x+e2x), first 0.5*(e0y+e2y) front + set label sprintf('$\theta_%d$', n-1) center at first e0x+0.5*cos(0.5*e_ts[n]), first e0y+0.5*sin(0.5*e_ts[n]) front + set label sprintf('$( x_%d, y_%d )$', n-1, n-1) center at first e0x, first e0y-0.25 front +} + +plot \ + NaN notitle diff --git a/docs/source/collision/data/src/ellipse.c b/docs/source/collision/data/src/ellipse.c new file mode 100644 index 0000000..6c87b64 --- /dev/null +++ b/docs/source/collision/data/src/ellipse.c @@ -0,0 +1,62 @@ +#include +#include "others.h" +#include "ellipse.h" + + +// compute x-coordinate of evolute +double ellipse_compute_xc(const double a, const double b, const double t){ + return a*(1.-pow(b/a, 2.))*pow(cos(t), 3.); +} + +// compute y-coordinate of evolute +double ellipse_compute_yc(const double a, const double b, const double t){ + return b*(1.-pow(a/b, 2.))*pow(sin(t), 3.); +} + +double ellipse_compute_curvature(const double a, const double b, const double t){ + const double num = a*b; + const double den = pow(pow(a*sin(t), 2.)+pow(b*cos(t), 2.), 1.5); + return num/den; +} + +// https://blog.chatfield.io/simple-method-for-distance-to-ellipse/ +double ellipse_find_normal_t(const double a, const double b, const double xp, const double yp){ + /* ! consider in the 1st quadrant ! 2 ! */ + double abs_xp = fabs(xp); + double abs_yp = fabs(yp); + /* ! initialise t ! 1 ! */ + double t = 0.25*M_PI; + // iterate until converged + while(1){ + /* ! compute point on the ellipse ! 2 ! */ + double xe = a*cos(t); + double ye = b*sin(t); + /* ! compute center of the fitted circle ! 2 ! */ + double xc = ellipse_compute_xc(a, b, t); + double yc = ellipse_compute_yc(a, b, t); + /* ! compute residual ! 8 ! */ + double dxe = xe-xc; + double dye = ye-yc; + double dxp = abs_xp-xc; + double dyp = abs_yp-yc; + double norme = hypot(dxe, dye); + double normp = hypot(dxp, dyp); + double dc = norme*asin((dxe*dyp-dye*dxp)/(norme*normp)); + double dt = dc/sqrt(a*a+b*b-xe*xe-ye*ye); + /* ! update t ! 4 ! */ + t += dt; + // limit to the 1st quadrant + t = fmin(0.5*M_PI, t); + t = fmax( 0., t); + /* ! terminate iteration when the residual is sufficiently small ! 3 ! */ + if(fabs(dt) < 1.e-8){ + break; + } + } + /* ! recover information of the quadrants ! 2 ! */ + if(xp < 0.) t = M_PI-t; + if(yp < 0.) t = -t; + /* ! return final t ! 1 ! */ + return t; +} + diff --git a/docs/source/collision/data/src/fit_circle.c b/docs/source/collision/data/src/fit_circle.c new file mode 100644 index 0000000..a81b645 --- /dev/null +++ b/docs/source/collision/data/src/fit_circle.c @@ -0,0 +1,44 @@ +#include +#include +#include "others.h" +#include "ellipse.h" +#include "example.h" + + +int fit_circle(void){ + /* ! initialise an ellipse ! 5 ! */ + const double e0_a = 2.0; + const double e0_b = 1.5; + const double e0_x = 0.5; + const double e0_y = 0.5; + const double e0_theta = 0.2; + /* ! set target point (xp, yp) ! 2 ! */ + const double e0_xp = 3.0; + const double e0_yp = 2.0; + /* ! transform coordinate, forward ! 2 ! */ + double e0_xp_, e0_yp_; + shift_and_rotate(-e0_x, -e0_y, -e0_theta, e0_xp, e0_yp, &e0_xp_, &e0_yp_); + /* ! find desired t ! 4 ! */ + double e0_t = ellipse_find_normal_t(e0_a, e0_b, e0_xp_, e0_yp_); + double e0_xc_ = ellipse_compute_xc(e0_a, e0_b, e0_t); + double e0_yc_ = ellipse_compute_yc(e0_a, e0_b, e0_t); + double e0_r = 1./ellipse_compute_curvature(e0_a, e0_b, e0_t); + /* ! transform coordinate, backward ! 2 ! */ + double e0_xc, e0_yc; + rotate_and_shift(+e0_x, +e0_y, +e0_theta, e0_xc_, e0_yc_, &e0_xc, &e0_yc); + // output + FILE *fp = fopen("fit-circle.dat", "w"); + fprintf(fp, "e0_a = % .15f\n", e0_a ); + fprintf(fp, "e0_b = % .15f\n", e0_b ); + fprintf(fp, "e0_x = % .15f\n", e0_x ); + fprintf(fp, "e0_y = % .15f\n", e0_y ); + fprintf(fp, "e0_theta = % .15f\n", e0_theta); + fprintf(fp, "e0_xp = % .15f\n", e0_xp ); + fprintf(fp, "e0_yp = % .15f\n", e0_yp ); + fprintf(fp, "e0_xc = % .15f\n", e0_xc ); + fprintf(fp, "e0_yc = % .15f\n", e0_yc ); + fprintf(fp, "e0_r = % .15f\n", e0_r ); + fclose(fp); + return 0; +} + diff --git a/docs/source/collision/data/src/fit_circles.c b/docs/source/collision/data/src/fit_circles.c new file mode 100644 index 0000000..7c47949 --- /dev/null +++ b/docs/source/collision/data/src/fit_circles.c @@ -0,0 +1,74 @@ +#include +#include +#include "others.h" +#include "ellipse.h" +#include "example.h" + + +int fit_circles(void){ + /* ! initialise geometry of ellipses ! 6 ! */ + const double e0_a = 2.0; + const double e0_b = 1.5; + const double e0_theta = 0.2; + const double e1_a = 1.5; + const double e1_b = 1.0; + const double e1_theta = 2.0; + for(int index = 0; index < 2; index++){ + /* ! initialise centers of ellipses ! 4 ! */ + const double e0_x = index == 0 ? 0.5 : 0.2; + const double e0_y = index == 0 ? 0.5 : 0.2; + const double e1_x = index == 0 ? 2.0 : 2.3; + const double e1_y = index == 0 ? 2.5 : 2.8; + /* ! initialise evolutes using the centers of ellipses ! 4 ! */ + double e0_xc = e0_x; + double e0_yc = e0_y; + double e1_xc = e1_x; + double e1_yc = e1_y; + double e0_r, e1_r; + for(int n = 0; n < 10; n++){ + /* ! transform coordinates, forward ! 4 ! */ + double e0_xc_, e0_yc_; + double e1_xc_, e1_yc_; + shift_and_rotate(-e0_x, -e0_y, -e0_theta, e1_xc, e1_yc, &e0_xc_, &e0_yc_); + shift_and_rotate(-e1_x, -e1_y, -e1_theta, e0_xc, e0_yc, &e1_xc_, &e1_yc_); + /* ! find desired t ! 2 ! */ + double e0_t = ellipse_find_normal_t(e0_a, e0_b, e0_xc_, e0_yc_); + double e1_t = ellipse_find_normal_t(e1_a, e1_b, e1_xc_, e1_yc_); + /* ! update center of fitted circles ! 4 ! */ + e0_xc_ = ellipse_compute_xc(e0_a, e0_b, e0_t); + e0_yc_ = ellipse_compute_yc(e0_a, e0_b, e0_t); + e1_xc_ = ellipse_compute_xc(e1_a, e1_b, e1_t); + e1_yc_ = ellipse_compute_yc(e1_a, e1_b, e1_t); + /* ! transform coordinates, backward ! 2 ! */ + rotate_and_shift(+e0_x, +e0_y, +e0_theta, e0_xc_, e0_yc_, &e0_xc, &e0_yc); + rotate_and_shift(+e1_x, +e1_y, +e1_theta, e1_xc_, e1_yc_, &e1_xc, &e1_yc); + // radii of the fitted circles + e0_r = 1./ellipse_compute_curvature(e0_a, e0_b, e0_t); + e1_r = 1./ellipse_compute_curvature(e1_a, e1_b, e1_t); + } + char fname[128]; + sprintf(fname, "fit-circles-%d.dat", index); + FILE *fp = fopen(fname, "w"); + fprintf(fp, "e0_a = % .15f\n", e0_a ); + fprintf(fp, "e0_b = % .15f\n", e0_b ); + fprintf(fp, "e0_x = % .15f\n", e0_x ); + fprintf(fp, "e0_y = % .15f\n", e0_y ); + fprintf(fp, "e0_theta = % .15f\n", e0_theta); + fprintf(fp, "e0_xc = % .15f\n", e0_xc ); + fprintf(fp, "e0_yc = % .15f\n", e0_yc ); + fprintf(fp, "e0_r = % .15f\n", e0_r ); + fprintf(fp, "e1_a = % .15f\n", e1_a ); + fprintf(fp, "e1_b = % .15f\n", e1_b ); + fprintf(fp, "e1_x = % .15f\n", e1_x ); + fprintf(fp, "e1_y = % .15f\n", e1_y ); + fprintf(fp, "e1_theta = % .15f\n", e1_theta); + fprintf(fp, "e1_xc = % .15f\n", e1_xc ); + fprintf(fp, "e1_yc = % .15f\n", e1_yc ); + fprintf(fp, "e1_r = % .15f\n", e1_r ); + fclose(fp); + } + return 0; +} + + + diff --git a/docs/source/collision/data/src/main.c b/docs/source/collision/data/src/main.c new file mode 100644 index 0000000..28049e4 --- /dev/null +++ b/docs/source/collision/data/src/main.c @@ -0,0 +1,9 @@ +#include "example.h" + + +int main(void){ + fit_circle(); + fit_circles(); + return 0; +} + diff --git a/docs/source/collision/data/src/others.c b/docs/source/collision/data/src/others.c new file mode 100644 index 0000000..702be59 --- /dev/null +++ b/docs/source/collision/data/src/others.c @@ -0,0 +1,22 @@ +#include +#include +#include "others.h" + + +// coordinate transformations +int shift_and_rotate(const double dx, const double dy, const double dt, const double x0, const double y0, double *x1, double *y1){ + *x1 = cos(dt)*(x0+dx) - sin(dt)*(y0+dy); + *y1 = sin(dt)*(x0+dx) + cos(dt)*(y0+dy); + return 0; +} + +int rotate_and_shift(const double dx, const double dy, const double dt, const double x0, const double y0, double *x1, double *y1){ + *x1 = dx + cos(dt)*x0 - sin(dt)*y0; + *y1 = dy + sin(dt)*x0 + cos(dt)*y0; + return 0; +} + +double get_random_value(const double min, const double max){ + return (max-min)*(1.*rand()/RAND_MAX)+min; +} + diff --git a/docs/source/collision/main.rst b/docs/source/collision/main.rst new file mode 100644 index 0000000..3fee3dc --- /dev/null +++ b/docs/source/collision/main.rst @@ -0,0 +1,223 @@ +##################### +Collision of Ellipses +##################### + +In this part I consider the collisions between two ellipses; namely, in this picture: + +.. image:: data/problem-setup.png + :width: 400 + +I aim at investigating + +* whether the two ellipses collide, +* if they do, how much the penetration depth :math:`\delta` is. + +I assume the ellipses are completely smooth and rigid (i.e. no deformation), and :math:`\delta` is sufficiently small compared to the size of the ellipses. + +.. note:: + + This is a rough solution to obtain something plausible just to avoid the over-penetration between objects. + Neither the collision detection nor the quantification of :math:`\delta` are rigorous. + +.. seealso:: + + `Documentation for Japanese readers `_. + +************* +Problem setup +************* + +An ellipse whose center is at the origin and the major / minor axes are :math:`a_0` and :math:`b_0` is given by + +.. math:: + + \frac{x_e^2}{a_0^2} + + + \frac{y_e^2}{b_0^2} + = + 1. + +By introducing a parameter :math:`t_0`, this equation reads + +.. math:: + + \begin{pmatrix} + x_e \\ + y_e + \end{pmatrix} + = + \begin{pmatrix} + a_0 \cos t_0 \\ + b_0 \sin t_0 + \end{pmatrix}. + +To take into account the rotation of :math:`\theta_0` and off-centred objects whose centre is at :math:`\left( x_0, y_0 \right)`, I modify the equation as + +.. math:: + + \begin{pmatrix} + x_e \\ + y_e + \end{pmatrix} + = + \begin{pmatrix} + x_0 \\ + y_0 + \end{pmatrix} + + + \begin{pmatrix} + \cos \theta_0 & -\sin \theta_0 \\ + \sin \theta_0 & \cos \theta_0 + \end{pmatrix} + \begin{pmatrix} + a_0 \cos t_0 \\ + b_0 \sin t_0 + \end{pmatrix}. + +Hereafter subscripts :math:`0,1` are used to distinguish the two ellipses. + +The collision / intersection of the two ellipses results in a quartic equation with respect to the intersecting point(s) :math:`x` (or :math:`y`). +The solution can be categorised as follows: + +#. Four real solutions + + Four intersecting points + +#. Two real and two imaginary solutions + + Two intersecting points + +#. Four imaginary solutions + + No Collision + +#. And more + +We are interested in distinguishing the second and the third cases. +Note that the other cases are out of focus, since we assume that the collision depth is small enough. + +Although we know a quartic formula, it is non-trivial to treat it numerically. +Also, let's say we obtain a theoretical solution; it is, however, still unclear how we define the penetration depth :math:`\delta`. + +For these reasons, instead of solving quartic equations, I consider a simpler and approximated solution in this project. + +************************************** +Collision check using circular fitting +************************************** + +Although ellipses can rotate and be off-centered, by transforming the coordinate system: + +.. math:: + \begin{pmatrix} + x_p \\ + y_p + \end{pmatrix} + \leftarrow + \begin{pmatrix} + \cos ( -\theta) & - \sin ( -\theta) \\ + \sin ( -\theta) & \cos ( -\theta) + \end{pmatrix} + \begin{pmatrix} + x_p-x_0 \\ + y_p-y_0 + \end{pmatrix}, + +we are able to say one object is at the origin and does not rotate (the major axis coincides with the Cartesian :math:`x` axis). + +====================== +Circular approximation +====================== + +.. note:: + + This part is largely inspired by `Simple Method for Distance to Ellipse `_. + +Recall that I consider situations where :math:`\delta` is much smaller than the sizes of the colliding objects. +This allows me to think that collisions between two circles fitting the ellipses in the vicinity of the colliding points well approximates the collisions between ellipses. + +Collisions between two circles are much simpler; the penetration depth of two circles (not a rigorous definition but one sound candidate) leads to + +.. math:: + + \delta \equiv r_0 + r_1 - d, + +where :math:`r_0` and :math:`r_1` are radii of the two circles, and :math:`d` is center-to-center distance. + +Now the main focus is how to approximate an ellipse using a circle in a systematic way. +Actually there is an analytical circle which approximates an ellipse locally: + +.. math:: + + \left( + x_c, + y_c + \right) + = + \left( + a \left( 1 - \frac{b^2}{a^2} \right) \cos^3 t, + b \left( 1 - \frac{a^2}{b^2} \right) \sin^3 t + \right), + +which is known as `evolute `_. +The local curvature is given by + +.. math:: + + \kappa + = + \frac{ + ab + }{ + \sqrt{\left( a^2 \sin^2 t + b^2 \cos^2 t \right)^3} + }, + +whose reciprocal is the radius of the circle fitting the original ellipse. + +The only unknown is :math:`t` which gives a circle nicely approximating the ellipse locally. + +Finding a nice :math:`t` is equal to find a normal vector from the point :math:`( x_p, y_p )` to the ellipse, and also is equal to find the minimum distance to the ellipse. + +Although these lemmas are not proved here, the following picture gives a good indication: + +.. image:: data/fit-circle.png + :width: 400 + +Here, + +* the black arrow connecting :math:`( x_p, y_p )` and :math:`( x_c, y_c )` gives a normal vector to the ellipse, + +* the fitted circle gives a good approximation of the ellipse locally. + +========================= +Collision of two ellipses +========================= + +I use the above method to quantify the penetration depth :math:`\delta`, which is very simple: for the ellipse :math:`i`, the center of the fitted circle of the ellipse :math:`j` :math:`( x_{c_j}, y_{c_j} )` is used as the target point :math:`( x_{p_i}, y_{p_i} )` to fit a circle. +This process is iterated until the locations of :math:`( x_{c_i}, y_{c_i} )` converge. + +When the two ellipses are colliding, the fitting circles lead + +.. image:: data/fit-circles-0.png + :width: 400 + +When the two ellipses are not colliding, the final state leads + +.. image:: data/fit-circles-1.png + :width: 400 + +Since we now know nice circles, the penetration depth is simply + +.. math:: + + \delta \equiv r_0 + r_1 - d, + +where, again, :math:`r_i` are radii of the fitted circles, while :math:`d` is the distance between the two centers of the fitting circles. +Obviously two circles collide when :math:`\delta > 0` and do not otherwise. + +************************** +Stand-alone implementation +************************** + +.. literalinclude:: data/example.py + :language: python + diff --git a/docs/source/conf.py b/docs/source/conf.py index 481908e..8a9753b 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -1,35 +1,12 @@ -# Configuration file for the Sphinx documentation builder. -# -# This file only contains a selection of the most common options. For a full -# list see the documentation: -# https://www.sphinx-doc.org/en/master/usage/configuration.html - -# -- Path setup -------------------------------------------------------------- - -# If extensions (or modules to document with autodoc) are in another directory, -# add these directories to sys.path here. If the directory is relative to the -# documentation root, use os.path.abspath to make it absolute, like shown here. -# import os import sys -# sys.path.insert(0, os.path.abspath(".")) - - -# -- Project information ----------------------------------------------------- project = "Ellises in Flows" -copyright = "2022, Naoki Hori" author = "Naoki Hori" +copyright = f"2022, {author}" -# The full version, including alpha/beta/rc tags release = "0.1" - -# -- General configuration --------------------------------------------------- - -# Add any Sphinx extension module names here, as strings. They can be -# extensions coming with Sphinx (named "sphinx.ext.*") or your custom -# ones. sys.path.append(os.path.abspath("./ext")) extensions = [ "myliteralinclude", @@ -37,33 +14,26 @@ "sphinx.ext.mathjax", ] -# Add any paths that contain templates here, relative to this directory. -templates_path = ["_templates"] - -# List of patterns, relative to source directory, that match files and -# directories to ignore when looking for source files. -# This pattern also affects html_static_path and html_extra_path. -exclude_patterns = [] - - -# -- Options for HTML output ------------------------------------------------- - -# The theme to use for HTML and HTML Help pages. See the documentation for -# a list of builtin themes. -# html_theme = "alabaster" html_theme_options = { - "fixed_sidebar": "true", - "github_user": "NaokiHori", + "fixed_sidebar": "false", + "github_banner": "false", + "github_button": "true", + "github_count": "true", "github_repo": "EllipsesInFlows", - "github_type": "true", + "github_type": "star", + "github_user": "NaokiHori", + "navigation_with_keys": "true", + "nosidebar": "false", + "page_width": "95vw", + "show_powered_by": "true", + "show_related": "false", + "show_relbars": "false", + "sidebar_collapse": "true", + "sidebar_includehidden": "false", + "sidebar_width": "20vw", } -# Add any paths that contain custom static files (such as style sheets) here, -# relative to this directory. They are copied after the builtin static files, -# so a file named "default.css" will overwrite the builtin "default.css". -html_static_path = ["_static"] - rst_prolog = """ .. role:: c-lang(code) :language: c diff --git a/docs/source/index.rst b/docs/source/index.rst index 1a8b4c8..7df74fa 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -2,10 +2,6 @@ Documentation ############# -.. note:: - - This library is still under development. - .. image:: snapshot.png :width: 800 @@ -15,6 +11,7 @@ Documentation introduction governing_equations/main numerical_method/main + collision/main implementation/main examples/index references