/*************************************************************************************************************************************/ /* */ /* D E T E R M I N A T I O N D E L ' E Q U A T I O N D ' U N E E L L I P S E : */ /* */ /* */ /* Author of '$xtc/EquationEllipse.02$vv$c' : */ /* */ /* Jean-Francois COLONNA (LACTAMME, 20221011183131). */ /* */ /*************************************************************************************************************************************/ #include "INCLUDES.01.I" extern double sqrt(); static int GenererImage=FAUX; /* Afin de pouvoir verifier qu'il s'agit bien d'une ellipse... */ #define ZOOM \ 2 #define FACTEUR \ 1 #define REDUC(z) \ ((z)/FACTEUR) #define CARRE(z) \ ((z)*(z)) #define CHECK(x,y) \ (CARRE(x)/CARRE(a))+(CARRE(y)/CARRE(b)) #define INFINI \ 1.0e100 #define NOMBRE \ 41 extern void *malloc(); static int dimX=0; #define Xmin 0 #define Xmax (Xmin + (dimX-1)) /* Definition des abscisses. */ static int dimY=0; #define Ymin 0 #define Ymax (Ymin + (dimY-1)) /* Definition des ordonnees. */ #define IMAGE(x,y) \ (*(image + (((y-Ymin)*dimX) + (x-Xmin)))) \ /* Acces a un point de l'image. */ #define store(n,x,y) \ { \ IMAGE(x,y) = n; \ } \ /* Rangement d'un point valide d'une image. */ #define ELLIPSE(Lx,Ly) \ { \ int cX=(((Lx[index]-minimum_2_X)/(maximum_2_X-minimum_2_X))*(Xmax-Xmin))+Xmin; \ int cY=(((Ly[index]-minimum_2_Y)/(maximum_2_Y-minimum_2_Y))*(Ymax-Ymin))+Ymin; \ \ cX = (((cX-Xmin)-((Xmax-Xmin)/2))/ZOOM)+(dimX/2)+Xmin; \ cY = (((cY-Ymin)-((Ymax-Ymin)/2))/ZOOM)+(dimY/2)+Ymin; \ \ store(BLANC,cX,cY); \ } \ /* Trace de l'ellipse... */ void main() { int index; double Lx1[NOMBRE]={-800000000000 ,-799166000064.454468 ,-796663511901.509888 ,-792491283632.184937 ,-786647243443.405396 ,-779128521231.340820 ,-769931489097.738647 ,-759051804715.996460 ,-746484479038.770020 ,-732223953778.206543 ,-716264213954.533447 ,-698598919432.062622 ,-679221575849.836548 ,-658125753714.453125 ,-635305359835.672241 ,-610754984887.974731 ,-584470346608.554565 ,-556448846633.977905 ,-526690275857.960571 ,-495197717087.342590 ,-461978695533.236389 ,-427046653020.436279 ,-390422834068.982605 ,-352138716161.846497 ,-312239139907.355164 ,-270786351349.391632 ,-227865229686.255005 ,-183590042709.856995 ,-138113159832.700897 ,-91636235028.870025 ,-44424416024.877166 ,+3175900065.465378 ,+50715503185.226310 ,+97617427472.923400 ,+143150528185.577606 ,+186406172868.677582 ,+226286945182.115662 ,+261521714723.657501 ,+290725629051.373169 ,+312520841585.724304 ,+325716668536.934021 }; double Ly1[NOMBRE]={+1000 ,+27989568947.760650 ,+55921375246.897300 ,+83735923255.243881 ,+111373744767.720016 ,+138774353113.066803 ,+165875883044.374756 ,+192614711478.117554 ,+218925054332.794769 ,+244738533244.007599 ,+269983706706.834137 ,+294585557945.559204 ,+318464929193.731323 ,+341537896320.715393 ,+363715071637.104004 ,+384900821220.337158 ,+404992390454.005859 ,+423878913452.105347 ,+441440297714.340515 ,+457545961486.339233 ,+472053411434.564392 ,+484806639343.538879 ,+495634329810.908386 ,+504347874400.330322 ,+510739199645.671082 ,+514578458131.797607 ,+515611666175.225647 ,+513558468956.451233 ,+508110328739.113464 ,+498929666769.876892 ,+485650807662.503113 ,+467884075341.501648 ,+445225103997.960083 ,+417272362288.405457 ,+383656913569.580444 ,+344089133453.590576 ,+298426404636.323669 ,+246761858210.234467 ,+189524738428.388153 ,+127567417275.454803 ,+62198223006.654739 }; double Lx2[NOMBRE]={+329515043551.962524 ,+323674532868.917603 ,+308564150133.509460 ,+285081892772.509705 ,+254475501947.702972 ,+218137658787.323669 ,+177436906761.184479 ,+133610597143.845917 ,+87716052194.003341 ,+40622578393.841774 ,-6973667863.918034 ,-54525136275.164947 ,-101606497534.667068 ,-147890330250.330505 ,-193126365338.586426 ,-237124520192.660980 ,-279741380978.299683 ,-320869592701.490845 ,-360429601866.614075 ,-398363259766.155457 ,-434628874995.758789 ,-469197382259.000793 ,-502049380352.338501 ,-533172829268.901184 ,-562561259638.667480 ,-590212378315.997192 ,-616126969700.465210 ,-640308042395.343994 ,-662760159839.397949 ,-683488909924.952026 ,-702500482089.756958 ,-719801340182.548706 ,-735397962263.524780 ,-749296632401.443970 ,-761503277891.523193 ,-772023340802.756470 ,-780861679188.498535 ,-788022488122.065308 ,-793509238390.982666 ,-797324630363.287231 ,-799470563564.701050 }; double Ly2[NOMBRE]={-4882226246.836296 ,-71755092708.735641 ,-136528805597.224564 ,-197575237592.599060 ,-253681601590.071869 ,-304090926985.294800 ,-348453819930.926392 ,-386735626516.374207 ,-419117050514.681030 ,-445909525249.911804 ,-467492196848.978088 ,-484269211220.109741 ,-496642941958.444824 ,-504998521330.766052 ,-509695854117.710144 ,-511066339505.118835 ,-509412401795.952881 ,-505008611894.748901 ,-498103619840.776733 ,-488922441210.968323 ,-477668816546.808472 ,-464527498768.808533 ,-449666385627.008789 ,-433238466668.945557 ,-415383573394.559753 ,-396229942333.709900 ,-375895603252.519348 ,-354489608547.058411 ,-332113120234.844666 ,-308860374532.823303 ,-284819536244.112122 ,-260073461792.595642 ,-234700379460.007568 ,-208774502202.248444 ,-182366580523.182190 ,-155544404924.997314 ,-128373267044.422150 ,-100916386548.792847 ,-73235308198.392166 ,-45390273860.971466 ,-17440583527.914199 }; /* Ces listes {x1,y1} et {x2,Y2} sont des points "diametralement opposes" d'une ellipse */ /* obtenue par integration des equations de Newton ('v $xiirk/NCOS_81' et aussi */ /* 'v $xiirk/$Fnota Debut_listG_NCOS_81_coordonnees_XY')... */ double minimum_1_X=+INFINI,maximum_1_X=-INFINI; double minimum_1_Y=+INFINI,maximum_1_Y=-INFINI; double Translation_X; double Translation_Y; for (index=0 ; index<NOMBRE ; index++) { minimum_1_X = MIN2(Lx1[index],minimum_1_X); minimum_1_X = MIN2(Lx2[index],minimum_1_X); maximum_1_X = MAX2(Lx1[index],maximum_1_X); maximum_1_X = MAX2(Lx2[index],maximum_1_X); minimum_1_Y = MIN2(Ly1[index],minimum_1_Y); minimum_1_Y = MIN2(Ly2[index],minimum_1_Y); maximum_1_Y = MAX2(Ly1[index],maximum_1_Y); maximum_1_Y = MAX2(Ly2[index],maximum_1_Y); /* Premier calcul des extrema des coordonnees afin de centrer l'ellipse... */ } Translation_X = (maximum_1_X+minimum_1_X)/2; Translation_Y = (maximum_1_Y+minimum_1_Y)/2; for (index=0 ; index<NOMBRE ; index++) { Lx1[index]=Lx1[index]-Translation_X; Ly1[index]=Ly1[index]-Translation_Y; Lx2[index]=Lx2[index]-Translation_X; Ly2[index]=Ly2[index]-Translation_Y; /* Operation de centrage de l'ellipse afin que son equation soit "canonique"... */ } if (GenererImage == FAUX) { double moyenne_a=0; double moyenne_b=0; int nombre=0; for (index=0 ; index<NOMBRE ; index++) { Lx1[index]=REDUC(Lx1[index]); Ly1[index]=REDUC(Ly1[index]); Lx2[index]=REDUC(Lx2[index]); Ly2[index]=REDUC(Ly2[index]); /* Operation eventuelle de reduction des coordonnees au cas ou elles seraient vraiment */ /* trop grandes... */ } for (index=0 ; index<NOMBRE ; index++) { double x1=Lx1[index]; double y1=Ly1[index]; double x2=Lx2[index]; double y2=Ly2[index]; /* Les points {x1,y1} et {x2,Y2} sont deux points quasiment "diametralement opposes"... */ double X1=CARRE(x1); double Y1=CARRE(y1); double X2=CARRE(x2); double Y2=CARRE(y2); if ((X1 != X2) && (Y1 != Y2)) /* Ce test est necessaire car, en effet, les coordonnees sont recentrees par rapport a */ /* leurs extrema. Il peut donc y avoir deux points "opposees" dont les coordonnees sont */ /* opposees et donc leurs carres sont egaux. Cela donnerait donc des 'iA's et/ou des 'iB's */ /* nuls... */ { double a,A,iA; double b,B,iB; double determinant=((Y2*X1)-(Y1*X2)); /* Les coordonnees ayant ete "centrees", on peut utiliser l'equation "canonique" de */ /* l'ellipse... */ /* */ /* On va travailler simplement sur les carres en passant de : */ /* */ /* 2 2 2 2 */ /* x1 y1 x2 y2 */ /* ----- + ----- = 1 ----- + ----- = 1 */ /* 2 2 2 2 */ /* a b a b */ /* */ /* a : */ /* */ /* X1 Y1 X2 Y2 */ /* ---- + ---- = 1 ---- + ---- = 1 */ /* A B A B */ /* */ /* puis a : */ /* */ /* iA.X1 + iB.Y1 = 1 iA.X2 + iB.Y2 = 1 */ /* */ /* {iA,iB} etant les inverses de {A,B}. On est donc en presence d'un systeme lineaire de */ /* deux equations a deux inconnues trivial a resoudre... */ iB = (X1-X2)/determinant; iA = (Y2-Y1)/determinant; A = 1/iA; B = 1/iB; a = sqrt(A); b = sqrt(B); moyenne_a = moyenne_a+a; moyenne_b = moyenne_b+b; printf("{{%+f,%+f},{%+f,%+f}} a=%+f b=%+f",x1,y1,x2,y2,a,b); printf(" Check=(%+f,%+f)\n" ,CHECK(x1,y1) ,CHECK(x2,y2) ); nombre++; } else { } } moyenne_a = moyenne_a/nombre; moyenne_b = moyenne_b/nombre; printf("\n"); printf("Moyenne : a=%+f b=%+f sur %d couples de points.\n",moyenne_a,moyenne_b,nombre); } else { int x,y; /* Definition des coordonnees des points de l'image. */ unsigned char *image; /* Definition de l'image a generer... */ double minimum_2_X=+INFINI,maximum_2_X=-INFINI; double minimum_2_Y=+INFINI,maximum_2_Y=-INFINI; Get(dimX,"dimX"); Get(dimY,"dimY"); /* Recuperation des dimensions en 'X' et en 'Y' de l'image a generer. */ image=(unsigned char*)malloc(dimY*dimX); /* Definition de l'image a generer... */ for (y=Ymin ; y<=Ymax ; y++) { for (x=Xmin ; x<=Xmax ; x++) { store(NOIR,x,y); /* Initialisation de l'image finale... */ } } for (index=0 ; index<NOMBRE ; index++) { minimum_2_X = MIN2(Lx1[index],minimum_2_X); minimum_2_X = MIN2(Lx2[index],minimum_2_X); maximum_2_X = MAX2(Lx1[index],maximum_2_X); maximum_2_X = MAX2(Lx2[index],maximum_2_X); minimum_2_Y = MIN2(Ly1[index],minimum_2_Y); minimum_2_Y = MIN2(Ly2[index],minimum_2_Y); maximum_2_Y = MAX2(Ly1[index],maximum_2_Y); maximum_2_Y = MAX2(Ly2[index],maximum_2_Y); /* Deuxieme calcul des extrema des coordonnees afin de centrer l'image... */ } for (index=0 ; index<NOMBRE ; index++) { ELLIPSE(Lx1,Ly1); ELLIPSE(Lx2,Ly2); } write(1,image,dimX*dimY); /* Sortie de l'image afin de verifier qu'il s'agit bien d'une ellipse... */ } }