; docformat = 'rst' ; ; NAME: ; interLineEllipse ; PURPOSE: ; Evaluate the coordinates of the intersection(s) between an ellipse ; and a segment or a line ; ;+ ; :Description: ; Evaluate the coordinates of the intersection(s) between an ellipse ; and a segment or a line ; ; :Categories: ; Maths, Intersection ; ; :Params: ; A: in, required, type=fltarr(2) ; Cartesian coordinates of the first point of segment [A-B] ; inB: in, required, type=fltarr(2) ; Cartesian coordinates of the second point of segment [A-B] ; PARAM: in, required, type=float ; The parameters of the ellipse:: ; P[0] = X center ; P[1] = Y center ; P[2] = Semi-major axis ; P[3] = Semi-minor axis ; P[4] = Angular tilt (rad) ; ; :Returns: The cartesian coordinates of the intersection ; ; :Keywords: ; FLAG: out, optional, type=integer ; Provides a summary of the intersection result:: ; 0: No intersection ; 1: One Intersection ; 2: Two Intersections ; 3: Null Segment or Degenerate Circle ; 4: Bug ; LINE: in, optional, type=boolean, default=0 ; Treat segment A-B as a line ; DIST: out, optional, type=float ; Distance of the intersection from A ; ; :Author: ; B.Carry (IMCCE) ; ; :History: ; Change History:: ; Original Version written in June 2013, B. Carry (IMCCE) ; 2013 Sep. - B.Carry (IMCCE) - Header cleaned ;- function interLineEllipse, A, inB, param, flag=flag, line=line, dist=dist ;-------------------------------------------------------------------------------- ;--I-- Initialization And Input Verification ----------------------------------- ;--I.1-- Check Presence of Arguments if N_params() LT 3 then begin print,'Syntax - Y = interLineEllipse( A, B, P )' return, -1 endif ;--I.2-- Flag Initialization flag=4 ;--I.3-- Null Segment or Degenerate Circle if (A(0) eq inB(0) and A(1) eq inB(1)) or param(0) eq 0 or param(1) eq 0 then begin flag=3 return, [0,0] endif ;-------------------------------------------------------------------------------- ;--II-- Shift and Rotate the System -------------------------------------------- ;--II.1-- Translate the System so A is at the Origin B = inB - float(A) C = param(0:1) - float(A) ;--II.2-- Length of Segment A-B distAB = sqrt( B(0)*B(0) + B(1)*B(1)); ;--II.3-- Rotate the System so B is on the positive X axis. theCos = B(0)/distAB theSin = B(1)/distAB newX = C(0)*theCos + C(1)*theSin C(1) = -C(0)*theSin + C(1)*theCos C(0) = newX tilt = param(4) - atan( B(1)/B(0) ) ;--II.4-- Compute Discrimant for Function Roots ;--II.4.1-- Rewrite Variables x0 = float(C(0)) ;-Center X y0 = float(C(1)) ;-Center Y rx = float(param(2)) ;-Radius ry = float(param(3)) ;-Radius c = cos(tilt) ;-Rotation s = sin(tilt) ;-Rotation ;--II.4.2-- Compute Discriminant Parts ;--II.4.2/A-- Squared c2 = c*c s2 = s*s rx2= rx*rx ry2= ry*ry ;--II.4.2/B-- Discriminant "A,B,C" dA = ry2*c2 + rx2*s2 dB = 2* ( rx2 * s * (-x0*s + y0*c ) - $ ry2 * c * ( x0*c + y0*s ) ) dC = ry2 * ( ( x0*c + y0*s )^2. ) + $ rx2 * ( (-x0*s + y0*c )^2. ) - $ rx2 *ry2 ;-------------------------------------------------------------------------------- ;--III-- Search for Intersections ----------------------------------------------- delta = dB*dB - 4*dA*dC ;------------------------------------------------------------ ;--III.1-- No intersection if delta lt 0 then begin flag=0 return, [0,0] endif else begin ;------------------------------------------------------------ ;--III.2-- Single intersection if delta eq 0 then begin ;--III.2.1-- Position in the A-B coordinate system crossPos = -dB / (2*dA) dist=crossPos ;--III.2.2-- Position in the original coordinate system xOut = A(0) + crossPos*theCos yOut = A(1) + crossPos*theSin ;--III.2.3-- Intersection in the segment or line-circle case if ( crossPos ge 0 and crossPos le distAB ) OR keyword_set(line) then begin flag=1 return, [xOut,yOut] ;--III.2.4-- Intersection outside the segment AND segment-circle case endif else begin flag=0 return, [0,0] endelse ;------------------------------------------------------------ ;--III.3-- Two intersections endif else begin ;--III.3.1-- Positions in the A-B coordinate system crossPos1 = ( -dB + sqrt(delta) ) / (2*dA) crossPos2 = ( -dB - sqrt(delta) ) / (2*dA) dist=[crossPos1, crossPos2] ;--III.3.2-- Position in the original coordinate system xOut1 = A(0) + crossPos1*theCos yOut1 = A(1) + crossPos1*theSin xOut2 = A(0) + crossPos2*theCos yOut2 = A(1) + crossPos2*theSin ;--III.3.3-- Double intersection in the segment or line-circle case if keyword_set(line) or $ ( ( crossPos1 ge 0 and crossPos1 le distAB ) AND $ ( crossPos2 ge 0 and crossPos2 le distAB ) ) then begin flag=2 return, [ [xOut1,yOut1], [xOut2,yOut2] ] endif else begin ;--III.3.4-- Both intersections outside the segment AND segment-circle case if ( ( crossPos1 lt 0 and crossPos1 gt distAB ) AND $ ( crossPos2 lt 0 and crossPos2 gt distAB ) ) then begin flag=0 return, [0,0] ;--III.3.5-- One intersections inside the segment AND segment-circle case endif else begin ;--III.3.5.1-- Intersection corresponds to the 1st point if ( crossPos1 ge 0 and crossPos1 le distAB ) then begin flag=1 return, [xOut1,yOut1] ;--III.3.5.2-- Intersection corresponds to the 2nd point endif else begin if ( crossPos2 ge 0 and crossPos2 le distAB ) then begin flag=1 return, [xOut2,yOut2] ;--III.3.5.3-- Deadend: no intersection in the segment endif else begin flag=0 return, [0,0] endelse endelse endelse endelse endelse endelse end