| 2201 | |
| 2202 | trackParabolaCore( Cy, Cx, a, b, Zy, Vy, Vx ); |
| 2203 | } |
| 2204 | |
| 2205 | return; |
| 2206 | } |
| 2207 | |
| 2208 | void trackEllipse( std::vector< std::vector<double> > &Cxs, std::vector< std::vector<double> > &Cys, |
| 2209 | const double *const PHI, const double *const Vx, const double *const Vy ) |
| 2210 | { |
| 2211 | const double &lambda1 = PHI[0]; const double &lambda2 = PHI[1]; |
| 2212 | const double &D = PHI[6]; const double &E = PHI[7]; const double &F = PHI[8]; |
| 2213 | |
| 2214 | assert( lambda1*lambda2 > 0 ); |
| 2215 | |
| 2216 | // lambda1*(X')^2 + lambda2*(Y')^2 + F = 0 |
| 2217 | // Y' = 0 => X'= sqrt( -F/lambda1 ), it means that -F/lambda1 > 0 |
| 2218 | if( -F/lambda1 <= 0 ) |
| 2219 | return; |
| 2220 | |
| 2221 | // lambda1 * X*X + lambda2 * Y*Y + F = 0 |
| 2222 | // <=> (-lambda1/F) * X*X + (-lambda2/F) * Y*Y = 1 |
| 2223 | // <=> (X/a)^2 + (Y/b)^2 = 1, 1/a = sqrt(-lambda1/F), 1/b = sqrt(-lambda2/F) |
| 2224 | const double a = sqrt( -F / lambda1 ); |
| 2225 | const double b = sqrt( -F / lambda2 ); |
| 2226 | |
| 2227 | // Ellipse is tangent to an edge of the triangle element, and localtes outside of the triangle |
| 2228 | // Examine both opposite sides (a,0) and (-a,0) are belong to inside of the triangle element |
| 2229 | if( !isInsideTriangle( -a, 0, Vx, Vy ) && !isInsideTriangle( -a, 0, Vx, Vy ) ) |
| 2230 | return; |
| 2231 | |
| 2232 | const double PI = atan(static_cast<double>(1)) * 4; |
| 2233 | const double c = 35*(32/(PI*PI*PI) - 96/(PI*PI*PI*PI)) - static_cast<double>(13)/12; |
| 2234 | const double p1 = c*b; |
| 2235 | const double p2 = c*a; |
| 2236 | |
| 2237 | // X' = X+D/(2*lambda1), Y' = Y+E/(2*lambda2) <=> X = X'-D/(2*lambda1), Y = Y'-E/(2*lambda2) |
| 2238 | // [x;y] = P[X;Y] |
| 2239 | |
| 2240 | std::vector<double> Cx, Cy; |
| 2241 | |
| 2242 | // starting point |
| 2243 | Cx.push_back( a ); |
| 2244 | Cy.push_back( 0 ); |
| 2245 | |
| 2246 | // quater 1 |
| 2247 | Cx.push_back( a ); Cy.push_back( p1 ); |
| 2248 | Cx.push_back( p2 ); Cy.push_back( b ); |
| 2249 | Cx.push_back( 0 ); Cy.push_back( b ); |
| 2250 | |
| 2251 | // quater 2 |
| 2252 | Cx.push_back( -p2 ); Cy.push_back( b ); |
| 2253 | Cx.push_back( -a ); Cy.push_back( p1 ); |
| 2254 | Cx.push_back( -a ); Cy.push_back( 0 ); |
| 2255 | |
| 2256 | // quater 3 |
| 2257 | Cx.push_back( -a ); Cy.push_back( -p1 ); |
| 2258 | Cx.push_back( -p2 ); Cy.push_back( -b ); |
| 2259 | Cx.push_back( 0 ); Cy.push_back( -b ); |
| 2260 | |