Browse Source

[nucleus] use 4 points to compute curvature instead of 3

Gaspard Jankowiak 11 months ago
parent
commit
b7146d1253
1 changed files with 56 additions and 9 deletions
  1. 56
    9
      src/Nucleus.jl

+ 56
- 9
src/Nucleus.jl View File

@@ -70,6 +70,61 @@ function g_p(x::Vector{Float64}, α::Float64)
70 70
     return -(2α*min.(α*x.-1, 0.0).*log.(α*x) .+ min.(α*x.-1, 0.0).^2 ./x)
71 71
 end
72 72
 
73
+function compute_curvature(c::NucleusCoords)
74
+    A = zeros(3,3)
75
+    b = zeros(3)
76
+
77
+    Nnuc = size(c.Y, 1)
78
+
79
+    for i in 1:Nnuc
80
+        A[1,1] = 8*( c.Y[c.circ_idx.m2[i],1]^2
81
+                    +c.Y[c.circ_idx.m1[i],1]^2
82
+                    +c.Y[i,1]^2
83
+                    +c.Y[c.circ_idx.p1[i],1]^2)
84
+        A[1,2] = 8*( c.Y[c.circ_idx.m2[i],1]*c.Y[c.circ_idx.m2[i],2]
85
+                     +c.Y[c.circ_idx.m1[i],1]*c.Y[c.circ_idx.m1[i],2]
86
+                     +c.Y[i,1]*c.Y[i,2]
87
+                     +c.Y[c.circ_idx.p1[i],1]*c.Y[c.circ_idx.p1[i],2])
88
+        A[2,1] = A[1,2]
89
+        A[1,3] = -4*( c.Y[c.circ_idx.m2[i],1]
90
+                    +c.Y[c.circ_idx.m1[i],1]
91
+                    +c.Y[i,1]
92
+                    +c.Y[c.circ_idx.p1[i],1])
93
+        A[3,1] = A[1,3]
94
+        A[2,2] = 8*( c.Y[c.circ_idx.m2[i],2]^2
95
+                    +c.Y[c.circ_idx.m1[i],2]^2
96
+                    +c.Y[i,2]^2
97
+                    +c.Y[c.circ_idx.p1[i],2]^2)
98
+        A[2,3] = -4*( c.Y[c.circ_idx.m2[i],2]
99
+                    +c.Y[c.circ_idx.m1[i],2]
100
+                    +c.Y[i,2]
101
+                    +c.Y[c.circ_idx.p1[i],2])
102
+        A[3,2] = A[2,3]
103
+        A[3,3] = 8
104
+
105
+        b[1] = 4*( c.Y[c.circ_idx.m2[i],1]^3
106
+                    +c.Y[c.circ_idx.m1[i],1]^3
107
+                    +c.Y[i,1]^3
108
+                    +c.Y[c.circ_idx.p1[i],1]^3
109
+                    +c.Y[c.circ_idx.m2[i],1]*c.Y[c.circ_idx.m2[i],2]^2
110
+                    +c.Y[c.circ_idx.m1[i],1]*c.Y[c.circ_idx.m1[i],2]^2
111
+                    +c.Y[i,1]*c.Y[i,2]^2
112
+                    +c.Y[c.circ_idx.p1[i],1]*c.Y[c.circ_idx.p1[i],2]^2)
113
+        b[2] = 4*( c.Y[c.circ_idx.m2[i],2]^3
114
+                    +c.Y[c.circ_idx.m1[i],2]^3
115
+                    +c.Y[i,2]^3
116
+                    +c.Y[c.circ_idx.p1[i],2]^3
117
+                    +c.Y[c.circ_idx.m2[i],1]^2*c.Y[c.circ_idx.m2[i],2]
118
+                    +c.Y[c.circ_idx.m1[i],1]^2*c.Y[c.circ_idx.m1[i],2]
119
+                    +c.Y[i,1]^2*c.Y[i,2]
120
+                    +c.Y[c.circ_idx.p1[i],1]^2*c.Y[c.circ_idx.p1[i],2])
121
+        b[3] = -0.25*(A[1,1] + A[2,2])
122
+
123
+        x = A\b
124
+        c.k[i] = 1/sqrt(x[1]^2 + x[2]^2 - x[3])
125
+    end
126
+end
127
+
73 128
 function recompute_nucleus_coords(c::NucleusCoords)
74 129
     delta_Y = c.Y - c.Y[c.circ_idx.m1,:]
75 130
     c.r[:] = sqrt.(sum(abs2, delta_Y; dims=2))
@@ -79,15 +134,7 @@ function recompute_nucleus_coords(c::NucleusCoords)
79 134
     delta_θ = rem.(c.θ - c.θ[c.circ_idx.m1], 2π, RoundNearest)
80 135
     c.θ[:] = cumsum([c.θ[1]; delta_θ[1:end-1]])
81 136
 
82
-    # http://cral-labo.univ-lyon1.fr/labo/fc/Ateliers_archives/ateliers_2005-06/cercle_3pts.pdf
83
-    xc = (0.5*(c.Y[c.circ_idx.p1,1].^2 - c.Y[:,1].^2 + c.Y[c.circ_idx.p1,2].^2 - c.Y[:,2].^2) ./ (c.Y[c.circ_idx.p1,2]-c.Y[:,2])
84
-          - 0.5*(c.Y[:,1].^2 - c.Y[c.circ_idx.m1,1].^2 + c.Y[:,2].^2 - c.Y[c.circ_idx.m1,2].^2) ./ (c.Y[:,2]-c.Y[c.circ_idx.m1,2])) ./ ((c.Y[:,1] - c.Y[c.circ_idx.m1,1])./(c.Y[:,2] - c.Y[c.circ_idx.m1,2]) - (c.Y[c.circ_idx.p1,1] - c.Y[:,1])./(c.Y[c.circ_idx.p1,2] - c.Y[:,2]))
85
-
86
-    yc = -(c.Y[:,1] - c.Y[c.circ_idx.m1,1])./(c.Y[:,2] - c.Y[c.circ_idx.m1,2]) .* xc + 0.5*(c.Y[:,1].^2 - c.Y[c.circ_idx.m1,1].^2 + c.Y[:,2].^2 - c.Y[c.circ_idx.m1,2].^2) ./ (c.Y[:,2]-c.Y[c.circ_idx.m1,2])
87
-    radius = sqrt.((xc - c.Y[:,1]).^2 + (yc - c.Y[:,2]).^2)
88
-    c.k[:] = 1 ./radius
89
-    c.k[:] = 0.5*(c.k + c.k[c.circ_idx.m1])
90
-
137
+    compute_curvature(c)
91 138
     c.η[:] = log.(c.r)
92 139
     c.L = sum(c.r)
93 140
 end

Loading…
Cancel
Save