'***************************************************************** '* This program seeks the equation of a conical passing through * '* 5 given points. * '* ------------------------------------------------------------- * '* SAMPLE RUN: * '* * '* SEEK A CONICAL PASSING THROUGH 5 POINTS * '* * '* X1 Y1 = 0 2 * '* X2 Y2 = 1 5 * '* X3 Y3 = 5 77 * '* X4 Y4 = -4 50 * '* X5 Y5 = 2 14 * '* * '* * '* The coefficients of the conical ax2+bxy+cy2+dx+ey+f = 0 are: * '* * '* a = 1.500000 * '* b = 0.000000 * '* c = -0.000000 * '* d = -0.000000 * '* e = -0.500000 * '* f = 1.000000 * '* * '* Basic version by J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '***************************************************************** ' Explanations: ' ------------ ' We seek the cartesian equation ax^2+bxy+cy^2+dx+ey+f=0 of the ' conical passing through 5 points: (x1,y1),(x2,y2)... to (x5,y5). ' The unknowns are a,b,c,d,e,f. Since the equation can be multiplied ' by any real number, the actual number of independant unknowns is 5. ' That is the reason why only five different points are necessary to ' determine the conical. The linear system to solve is a system of ' five homogeneous equations with 6 unknowns. ' The resolution of the system will be made by successive steps: ' Let us suppose that f<>0 and let us fix f=1, then the other ' parameters a,b,c,d,e are solutions of the linear system: ' x1^2 a + x1y1 b + y1^2 c + x1 d + y1 e = -1 ' x2^2 a + x2y2 b + y2^2 c + x2 d + y2 e = -1 ' ------------------------------------------- ' x5^2 a + x5y5 b + y5^2 c + x5 d + y5 e = -1 ' This linear system of size 5 is solved by a classical Gauss-Jordan ' method found in procedure SolveSystem (see also program MATRICES/sysmat.bas). ' if the system has a solution (det <> 0) the problem is finished, else ' the hypothesis f<>0 was wrong and we then try e<>0 (that is e=1). ' The linear system to solve is then: ' x1^2 a + x1y1 b + y1^2 c + x1 d + f = -y1 ' x2^2 a + x2y2 b + y2^2 c + x2 d + f = -y2 ' ------------------------------------------- ' x5^2 a + x5y5 b + y5^2 c + x5 d + f = -y5 ' if the system has a solution (det <> 0) the problem is finished, else ' the hypothesis e<>0 was again wrong and we successively explore d<>0, ' c<>0... to a<>0. ' So in the most unfavourable case, that is if only a degenerated conical ' answers the problem, we can have up to six linear systems of size 5 ' to solve. ' Usually (f<>0 or e<>0), we only solve one or two linear systems. '---------------------------------------------------------------------------------- defdbl a-h,o-z defint i-n 'integer SeekConical% (return value of Sub 2000: 1=success 0=failure) DIM AA(5,5) 'left hand matrix DIM BB(5,1) 'right hand matrix (here number of columns=1) EPSMACH=2E-16 'machine dependant smallest real number 'input data section cls print print " SEEK A CONICAL PASSING THROUGH 5 POINTS" print INPUT " X1 Y1 = ", x1, y1 INPUT " X2 Y2 = ", x2, y2 INPUT " X3 Y3 = ", x3, y3 INPUT " X4 Y4 = ", x4, y4 INPUT " X5 Y5 = ", x5, y5 print print gosub 2000 'call SeekConical 'write results section F\$="#####.######" if SeekConical%=1 then print " The coefficients of the conical ax2+bxy+cy2+dx+ey+f = 0 are:" print print " a = "; : print using F\$; a print " b = "; : print using F\$; b print " c = "; : print using F\$; c print " d = "; : print using F\$; d print " e = "; : print using F\$; e print " f = "; : print using F\$; f else print " No conical found." end if print print END 'of main program '******************************************* '* SOLVING A LINEAR MATRIX SYSTEM AX = B * '* with Gauss-Jordan method using full * '* pivoting at each step. During the pro- * '* cess, original AA and BB matrices are * '* destroyed to spare storage location. * '* --------------------------------------- * '* INPUTS: AA MATRIX N*N * '* BB MATRIX N*M * '* --------------------------------------- * '* OUTPUTS: AA INVERSE OF AA N*N * '* DET DETERMINANT OF AA * '* BB SOLUTION MATRIX N*M * '* --------------------------------------- * '* NOTA - If M=0 inversion of AA matrix * '* only. * '* * '* Basic version by J-P Moreau * '* (here N=5, M=1). * '******************************************* 'Subroutine SolveSystem(N, M, AA, BB, DET) 'N,M,AA,BB,DET must be defined by calling program 1000 DIM PC(5), PL(5), CS(5) 'auxiliary vectors 'LABELS: 100, 105, 110, 150 ' Initializations DET = 1 for I=1 to N PC(I) = 0# PL(I) = 0# CS(I) = 0# next ' Main loop for K=1 to N ' Searching greatest pivot PV=AA(K,K) IK=K JK=K PAV=abs(PV) for I=K to N for J=K to N temp = abs(AA(I,J)) if (temp > PAV) then PV=AA(I,J) PAV=abs(PV) IK=I JK=J end if next J next I ' Search terminated, the pivot is in location I=IK, J=JK. ' Memorizing pivot location: PC(K)=JK PL(K)=IK ' DETERMINANT DET is actualised ' If DET=0, RETURN (main program takes care) ' Machine dependant EPSMACH equals here 2E-16 if IK <> K then DET=-DET if JK.<> K then DET=-DET DET=DET*PV temp= abs(DET) if (temp < EPSMACH) then goto 150 ' POSITIONNING PIVOT IN K,K: if IK <> K then for I=1 to N ' EXCHANGE LINES IK and K of matrix AA: TT=AA(IK,I) AA(IK,I)=AA(K,I) AA(K,I)=TT next end if if M <> 0 then for I=1 to M TT=BB(IK,I) BB(IK,I)=BB(K,I) BB(K,I)=TT next end if ' OT is at correct line if JK <> K then for I=1 to N ' EXCHANGE COLUMNS JK and K of matrix AA: TT=AA(I,JK) AA(I,JK)=AA(I,K) AA(I,K)=TT next end if ' The PIVOT is at correct column. ' and is located in K,K. ' Column K of matrix AA is stored in CS vector ' then column K is set to zero. for I=1 to N CS(I)=AA(I,K) AA(I,K)= 0.d0 next CS(K)= 0# AA(K,K)= 1# ' Line K of matrix AA is modified: temp = abs(PV) 'if pivot is too small, return with message (not very frequent) if temp < EPSMACH then print " Sorry, pivot too small, program fails!" goto 150 end if for I=1 to N AA(K,I)=AA(K,I)/PV next if M <> 0 then for I=1 to M BB(K,I)=BB(K,I)/PV next end if ' Other lines of matrix AA are modified: for J=1 to N if J=K then goto 100 for I=1 to N ' Line J of matrix AA is modified: AA(J,I)=AA(J,I)-CS(J)*AA(K,I) next if M <> 0 then for I=1 to M BB(J,I)=BB(J,I)-CS(J)*BB(K,I) next end if 100 next J ' Line K is ready. next K ' end of K loop ' MATRIX AA INVERSION IS forNE - REARRANGEMENT OF MATRIX AA ' EXCHANGE LINES for I=N to 1 step -1 IK=int(PC(I)) if IK=I then goto 105 ' EXCHANGE LINES I and PC(I) of matrix AA: for J=1 to N TT=AA(I,J) AA(I,J)=AA(IK,J) AA(IK,J)=TT next if M <> 0 then for J=1 to M TT=BB(I,J) BB(I,J)=BB(IK,J) BB(IK,J)=TT next end if ' NO MORE EXCHANGE is NECESSARY ' Go to next line. 105 next I ' end of loop for I=N to 1 step -1 ' EXCHANGE COLUMNS for J=N to 1 step -1 JK=int(PL(J)) if JK=J then goto 110 ' EXCHANGE COLUMNS J ET PL(J) of matrix AA: for I=1 to N TT=AA(I,J) AA(I,J)=AA(I,JK) AA(I,JK)=TT next ' NO MORE EXCHANGE is NECESSARY ' Go to next column. 110 next J ' REARRANGEMENT TERMINATED. 150 RETURN ' END of Sub 1000 '************************************************** '* This subroutine calls Subroutine 1000 to solve * '* the linear system AA*X=BB and rearranges the * '* solution BB according to ii case selected. * '************************************************** 'subroutine used by Subroutine 2000 (SeekConical) 'Subroutine Solve(ii, AA, BB, a, b, c, d, e, f,det) 1900 n=5: m=1: gosub 1000 ' call SolveSystem(5,1,AA,BB,det) if det<>0 then if ii=1 then ' case f<>0 a=BB(1,1): b=BB(2,1): c=BB(3,1): d=BB(4,1): e=BB(5,1): f=1# end if if ii=2 then ' case e<>0 a=BB(1,1): b=BB(2,1): c=BB(3,1): d=BB(4,1): e=1#: f=BB(5,1) end if if ii=3 then ' case d<>0 a=BB(1,1): b=BB(2,1): c=BB(3,1): d=1#: e=BB(4,1): f=BB(5,1) end if if ii=4 then ' case c<>0 a=BB(1,1): b=BB(2,1): c=1#: d=BB(3,1): e=BB(4,1): f=BB(5,1) end if if ii=5 then ' case b<>0 a=BB(1,1): b=1#: c=BB(2,1): d=BB(3,1): e=BB(4,1): f=BB(5,1) end if if ii=6 then ' case a<>0 a=1#: b=BB(1,1): c=BB(2,1): d=BB(3,1): e=BB(4,1): f=BB(5,1) end if end if return '************************************************************ '* This subroutine seeks a conical passing through 5 points * '* -------------------------------------------------------- * '* INPUT: x1,y1,x2,y2... x5,y5 (the five given points) * '* OUTPUT: a,b,c,d,e,f the 6 coefficients of conical * '* equation: ax2+bxy+cy2+dx+ey+f = 0 * '* -------------------------------------------------------- * '* USES: Subroutine 1900 calling Subroutine 1000. * '************************************************************ 'Function SeekConical(x1,y1,x2,y2,x3,y3,x4,y4,x5,y5,a,b,c,d,e,f) 'Label: 2100 2000 DIM x(5), y(5) 'auxiliary vectors SeekConical%=0 x(1)=x1: y(1)=y1 x(2)=x2: y(2)=y2 x(3)=x3: y(3)=y3 x(4)=x4: y(4)=y4 x(5)=x5: y(5)=y5 ' verify that the 5 given points are different for i=2 to 5 for j=1 to i-1 if x(i)=x(j) and y(i)=y(j) then return next j next i ' try f<>0 for i=1 to 5 AA(i,1)=x(i)*x(i): AA(i,2)=x(i)*y(i): AA(i,3)=y(i)*y(i) AA(i,4)=x(i): AA(i,5)=y(i): BB(i,1)=-1# next 'call Solve(1,AA,BB,a,b,c,d,e,f,det) ii=1: gosub 1900 if det <> 0 then goto 2100 ' try e<>0 for i=1 to 5 AA(i,1)=x(i)*x(i): AA(i,2)=x(i)*y(i): AA(i,3)=y(i)*y(i) AA(i,4)=x(i): AA(i,5)=1#: BB(i,1)=-y(i) next 'call Solve(2,AA,BB,a,b,c,d,e,f,det) ii=2: gosub 1900 if det <> 0 then goto 2100 ' try d<>0 for i=1 to 5 AA(i,1)=x(i)*x(i): AA(i,2)=x(i)*y(i): AA(i,3)=y(i)*y(i) AA(i,4)=y(i): AA(i,5)=1#: BB(i,1)=-x(i) next 'call Solve(3,AA,BB,a,b,c,d,e,f,det) ii=3: gosub 1900 if det <> 0 then goto 2100 ' try c<>0 for i=1 to 5 AA(i,1)=x(i)*x(i): AA(i,2)=x(i)*y(i): AA(i,3)=x(i) AA(i,4)=y(i): AA(i,5)=1#: BB(i,1)=-y(i)*y(i) next 'call Solve(4,AA,BB,a,b,c,d,e,f,det) ii=4: gosub 1900 if det <> 0 then goto 2100 ' try b<>0 for i=1 to 5 AA(i,1)=x(i)*x(i): AA(i,2)=y(i)*y(i): AA(i,3)=x(i) AA(i,4)=y(i): AA(i,5)=1#: BB(i,1)=-x(i)*y(i) next 'call Solve(5,AA,BB,a,b,c,d,e,f,det) ii=5: gosub 1900 if det <> 0 then goto 2100 ' last try a<>0 for i=1 to 5 AA(i,1)=x(i)*y(i): AA(i,2)=y(i)*y(i): AA(i,3)=x(i) AA(i,4)=y(i): AA(i,5)=1#: BB(i,1)=-x(i)*x(i) next 'call Solve(6,AA,BB,a,b,c,d,e,f,det) ii=6: gosub 1900 if det <> 0 then goto 2100 else return end if 2100 SeekConical%=1 return ' det<>0 means success ' end of file conical1.bas