(*
           C'est un "petit" problème qui m'a assez longtemps
           agacé pour que je lui consacre un peu de temps.
           Evidemment pour l'expert ce genre de travail est
           monnaie courante, voire un sujet de pure dérision.

           Mais quand on est un néophite il n'est pas désagréable
           de tomber sur ce genre d'explication et d'utilitaire.

           L'objet de ce programme:

            Relier DEUX spheres dans l'espace
            par un cone "tangent" commun . Pour
            passer les coordonnées résultantes
            à MNM2.1.

            Ou éventuellement à Povray.

            ATTENTION!!!

            Notez que la version 2 envoie des données directement
            à povray et lance ce dernier pour visualiser le résultat.
            Ce type de chose fait planter misérablement TP7 deux fois
            sur deux et demi, quand à bp7 il refuse purement et
            simplement de lancer une application 32 bits.
            A noter aussi que même un truc aussi  "évolué" que tp7
            continue à utiliser des fcbs (!!) que l'on pouvait
            raisonablement penser comme remplacés par "files".

            PAR CONTRE ... ça ne pose aucun problème à partir du dos.

   A noter pour un éventuel retour sur ce problème
   dans les mois à venir:

     1/ la solution passe par une mise à plat du problème.
        Un plan 2D.

     2/ Dans ce cas la tangente à deux cercles forme avec les
        deux rayons qui lui sont perpendiculaires un début
        de quadrilatère qu'il faut compléter.

   Organiser les données permet de faciliter le traitement.

   PAr exemple certaines opérations de déplacement se font
   par rapport à la dimension du rayon de la sphère. Il est
   préférable d'avoir le plus grand en premier. Mais on ne
   va pas s'amuser à refuser les données qui ne respecteraient
   pas cet ordre. Voir la procedure hierarchisation_sphere.

   Il doit y avoir translation et rotation pour positionner
   correctement le cone. Les rotation se font oujours par
   rapport au centre de la plus grande sphère. Et dans povray
   elles se font par rapport à l'origine du monde 3d en  coor
   données [0,0,0]. Un transformation des données pour que le
   centre de la sphère corresponde à l'origine s'impose. Ca
   permettra de faciliter certains calculs pour les rotations.
   Et  cela fixe d'entrée les valeurs de translation.
   Voir la procedure Translation_.
*)


{$A+,B-,D+,E+,F-,G+,I+,L+,N+,O-,P-,Q-,R-,S+,T-,V+,X+,Y+}
{$M $2000,0,0}
program liensph1;
uses crt,dos;

 type   reel=double;
        sphere3d=record x,y,z,
                         r:reel; end;
        cone=record
                     gr,pr,
                     dim,
                     angzy,angzx,
                     tx,ty,tz: reel;
             end;

   const

    Povpath   :string[16]  =  'c:\graphic\pov\';
    povfile   :string[12]  =  'lientst.pov';
    povini    :string[12]  =  'tst.cmd';
    tgafile   :string[12]  =  'lientst.tga';

    deb_povfile:array[1..6]of string[24]=(
                 (**)    '// JM Soler Opérateur de',
                 (**)    '// Liens de Sphères v0.3',
                 (**)    '// pour       POVray 2.2',
                 (**)    '//           Avril 1996 ','','');

    errdon: array[1..2] of string[22]=
                 (**)          ('  Erreur de donnée: ',
                 (**)           '  Erreur de valeur: ');
    errparam: array[1..4] of string[40]=
                 (**)         ('  Pas de paramètre en ligne de commande!',
                 (**)          '  Nombre de paramétre insuffisant!',
                 (**)          '  Un des paramétres n''est pas un réel!',
                 (**)          ''
                 (**)          );
    aidelogic: array[1..6] of string[50]=
                 (**)         ('  ---------------------------------------------',
                 (**)          '  | Liaison de Sphères :                      |',
                 (**)          '  | L''objectif de ce logiciel est de relier  |',
                 (**)          '  | deux sphères localisées dans l''espace 3d |',
                 (**)          '  | à l''aide d''un cone tangent.             |',
                 (**)          '  ---------------------------------------------'
                 (**)          );

    aideparam: array[1..9] of string[50]=
                 (**)         ('  ---------------------------------------------',
                 (**)          '  | Aide.                                     |',
                 (**)          '  | Les données doivent être fournies sur la  |',
                 (**)          '  | ligne de commande  dans l''ordre suivant: |',
                 (**)          '  |    Rayon1, coord x1, coord y1, coord z1,  |',
                 (**)          '  |    Rayon2, coord x2, coord y2, coord z2 . | ',
                 (**)          '  -------------------------------------------- |',
                 (**)          '    ',
                 (**)          ' Jm Soler 1996');

    chiffre:string=('0123456789.-+');

    string_doc: array[1..21] of string[45]=
    (*  1 *)(' Distance des deux spheres sur le plan XZ:',
    (*  2 *) ' Distance réelle des deux spheres:',
    (*  3 *) ' Hiérarchisation  des cercles:',
    (*  4 *) ' Ouverture d''angle entre le rayon sphere et',
    (*  5 *) ' le rayon cone:',
    (*  6 *) ' Plus grand rayon du cone:',
    (*  7 *) ' Plus petit rayon du cone:',
    (*  8 *) ' Décalage du cone à partir du grand cercle :',
    (*  9 *) ' Décalage du cone à partir du petit cercle :',
    (* 10 *) ' Dimension du cone:',
    (* 11 *) '',
    (* 12 *) ' n2 est le petit coté du triangle rectangle ',
    (* 13 *) ' dont d12 est l''hypothénuse.',
    (* 14 *) '',
    (* 15 *) ' n3 est le grand coté du triangle rectangle ',
    (* 16 *) ' dont d12 est l''hypothénuse.',
    (* 17 *) '',
    (* 18 *) ' Rotation sur le plan ZY :',
    (* 19 *) ' Rotation sur le plan ZX :',
    (* 20 *) '',
    (* 21 *) ' Translation réelle du cone :'
    (*  *)      );

var o1,o2,
    a,b,c,d,
    trans_          :sphere3d;
    l0,l1,l2        :string;
    pc,n            :word;
    r3d             :array[1..8] of reel;
    D12, Nxz        :reel;
    n0,  n1,
    n2,  n3,
    long_c,
    PGray_c, Ppray_c,
    DPcone,  Dgcone,
    Dimcone, ang1,
    Ang_zx,ang_zy   :reel;
    k1,k2           :integer;


(*
--------------------------------+-------+--------+------+-------+--------+
+-------------------------------┤ o1,o2 |        |      |       |        |
|          \   variables---->   |       | Pgray_c|      |       |        |
|            \                  |       | ppray_c|      |       |        |
| procedures   \                |       |        |dgcone|       |        |
|    et          \______________|       |        |dpcone|       |        |
| fonctions       | u=utilisées |       |        |      |dimcone|        |
|     |           | p=produites |       |        |      |       |ang_zy  |
|    \|/          |             |       |        |      |       |ang_zx  |
+-------------------------------+-------+--------+------+-------+--------┤
|proc hierarchisation_cercle;   | *  * u|        |      |       |        |
|proc translation_;             |    * u|        |      |       |        |
|proc translation_2;            | *  * u|        |      |       |        |
|proc distance_sphere;          |    * u|        |      |       |        |
|proc PGC_triangle;             | *  * u|        |      |       |        |
|proc pgcrayon_cone;            |       |  *  p  |  * p |  *  p |        |
|func rotation_(x1,y1:reel):reel|       |        |      |  *  p |  * * p |
|proc rotation;                 |       |        |      |       |  * * p |
|proc translation_relle;        |       |        |      |       |        |
|                               |       |        |      |       |        |
+-------------------------------+-------+--------+------+-------+--------+

  Structure de o1 et o2

                      x, |
                      y, |-> coordonnées 3d
                      z, |

                      r  : rayon

  Ces données peuvent être produites ou lues n'importe oŚ et n'importe
  comment...

  Dans l'idéal une procedure Prod_cone devrait pouvoir recevoir
  o1 et o2 en entrée et fournir le cone en sortie.

*)

procedure hierarchisation_cercle;
begin
     (*-------------------------------*
      Hiérarchisation  des cercles.

       Quel que soit l'ordre dans lequel sont
       fournies les données on range la plus
       grande sphère en premier.
       L'ordre dans lequel elles apparaissent
       ne joue pas sur le résultat mais facilite
       le traitement.

       Si la seconde est plus grande on permute.

       Si elles sont identique on  ne fait rien.
       Le traitement se fait toujours sur la
       première série de données
      *-------------------------------*)
     if o1.r0 then ang1:=arctan(n2/n3)
           else
            if n2>0 then ang1:=pi/2
            else
            if n2=0 then ang1:=0
          ;
     writeln(string_doc[4]);
     writeln(string_doc[5]);

     writeln(' Angle = ',ang1/pi*180:4:6);
end;


procedure pgcrayon_cone;
begin
   Pgray_c:=cos(ang1)*n0;
   Ppray_c:=cos(ang1)*n1;
       writeln(string_doc[6]);
         writeln('  ',Pgray_c:4:6);
       writeln(string_doc[7]);
         writeln('  ',Ppray_c:4:6);

   DGcone:=sin(ang1)*n0;
   Dpcone:=sin(ang1)*n1;
      writeln(string_doc[8]);
        writeln('  ',DGcone:4:6);
      writeln(string_doc[9]);
        writeln('  ',Dpcone:4:6);

   Dimcone:=(D12+Dpcone)-dgcone;
        writeln(string_doc[10]);
        writeln('  ',Dimcone:4:6);
end;

function rotx(x,y,r:reel):reel;
var z:reel;
begin
     z:=cos(r)*x-sin(r)*y;
     rotx:=z;
end;

function roty(x,y,r:reel):reel;
var z:reel;
begin
     z:=sin(r)*x+cos(r)*y;
     roty:=z;
end;

function rotation_(x1,y1:reel):reel;
var ang0:reel;
    begin
    if x1<>0 then
           ang0:=(pi/2-arctan(abs(y1/x1)))
                          else
                          begin
                             if y1>0 then ang0:=0
                                  else
                                  if y1<0 then ang0:=pi;
                          end;

             (* si nxz est égal à zéro cela signifie que
                la sphère 2 est sur l'axe y.
                Dans ce cas deux possibilités se présentent:
                soit les centres sont confondus, soit le
                second est légèrement plus loin sur cet axe,
                de toute façon les deux spheres ne pouvant pas
                être l'une dans l'autre, la première possibilité
                est exclue dès le départ; dans le second cas, le
                cone existe mais un problème d'orientation se pose
                car tous les calculs sont en positifs l'orientation
                pointe en Y  donc soit  y est positif et il n'y
                a rien à  faire soit y est négatif et il faut
                tourner de 180 degrés.
             *)


       if x1>0  then
                  begin
                     if y1>=0 then ang0:=ang0 else
                        if y1<0 then ang0:=pi-ang0;
                  end;

       if x1<0  then
                  begin
                     if y1>=0 then ang0:=2*pi-ang0 else
                        if y1<0 then ang0:=pi+ang0  ;
                  end;


   rotation_:=ang0;
end;

procedure rotation;
begin
(*
  En deux partie:

  1/ On se sert des coordonnées de o2 puisqu'elles sont
     correctement liées à l'origine.
     Arctan(o2.z/o2.x) donne l'angle de rotation sur le plan ZX

  2/ Nxz=sqrt(sqr(o2.z)+sqr(o2.x)) donne l'élément manquant pour obtenir
     Arctan(o2.y/Nxz) la rotation en ZY.

  Le cone étant orienté en [0,1,0]
*)

             ang_zy:=rotation_(nxz,o2.y);
             writeln(string_doc[18]);
             writeln('  ',ang_zy*180/pi:4:6, '(',o2.y:2:4,'/',nxz:2:4,')');

             ang_zx:=rotation_(o2.x,o2.z);
             writeln(string_doc[19]);
             writeln('  ', ang_zx*180/pi:4:6, '(',o2.z:2:4,'/',o2.x:2:4,')');

end;

var  x,y,z,
     x1,y1,z1:reel;
     t,t1:integer;
procedure translation_reelle;
begin
  x1:=0;
  y1:=dgcone;
     (*
        Remarque :   pour obtenir un positionnement correct
        avec Povray, ce dernier étant orienté "main gauche", et
        pas seulement sur l'axe de rotation x,i l faut inverser
        le sens de rotation trigonométrique !!
     *)
  z:=rotx(x1,y1,2*pi-ang_zy);
  y:=roty(x1,y1,2*pi-ang_zy);

  y1:=z;

  x:=rotx(x1,y1,2*pi-ang_zx);
  z:=roty(x1,y1,2*pi-ang_zx);

 (*
  writeln;
  gotoxy(40,0);
  writeln(string_doc[14],' x= ',x:4:4,
                       ';  y= ', y:4:4,
                       ' ;  z= ',z:4:4,';') ;
  *)
                  for t:=1 to  24 do
                         begin
                            gotoxy(46,t);
                            write(#179);
                          end;
  x:=x+trans_.x;
  y:=y+trans_.y;
  z:=z+trans_.z;

   gotoxy(47,1);
   write(string_doc[21]);
   gotoxy(49,2);
   write(' x= ',x:4:4);
   gotoxy(49,3);
   write(' y= ', y:4:4);
   gotoxy(49,4);
   write(' z= ',z:4:4,';') ;

end;

procedure calcul_;
var c:char;
begin
            hierarchisation_cercle;
            translation_;
            distance_sphere;
            pgc_triangle;
            pgcrayon_cone;
            rotation;
            translation_reelle;
            repeat
            if keypressed then c:=readkey;
            until c=#27;

end;

procedure doc_;
begin
     writeln;
     for pc:=1 to 6 do  writeln(aidelogic[pc]);
     for pc:=1 to 9 do  writeln(aideparam[pc]);
     writeln;
end;

procedure test_param;
begin
     if paramcount=0 then
        begin
             writeln(errparam[1],' Paramcount: ', paramcount);
             doc_;
             halt;
        end;

     if paramcount<8 then
        begin
             writeln(errparam[2],' , ',
                8-paramcount,'  paramètre(s) manque(nt) ');
             doc_;
             halt;
        end;
end;

procedure transforme_param(n:integer; var x:reel);
var y,errcode:integer;
begin
     l0:= paramstr(n);
     writeln('Analyse  des paramètres: ', n );
     writeln(paramstr(n));
     for y:=1 to length(l0) do
     begin
      if pos(l0[y],chiffre)=0 then
         begin
              writeln;
              writeln(errdon[1],errparam[3]);
              writeln('Paramètre en ligne numéro : ',n,'.' );
              writeln('Le caractère ',#39,l0[y],
              #39,' n''est pas autorisé !');
              halt(0);
         end
        else writeln(l0[y], ' OK ,');
     end;
    val(l0,x,errcode);
    writeln('Valeur réelle :', x:4:6);
    writeln;
end;

procedure transfert_donnees;
begin
     o1.r:=r3d[1];
     o1.x:=r3d[2];
     o1.y:=r3d[3];
     o1.z:=r3d[4];

     o2.r:=r3d[5];
     o2.x:=r3d[6];
     o2.y:=r3d[7];
     o2.z:=r3d[8];
end;

procedure doc_povini;
var f:text;
begin
     assign(f,povini);
     rewrite(f);
     writeln(f,' +I',povfile,' +O',tgafile, ' +L',povpath,'include\  +DGt ');
     writeln(f,' +H200 +W320 +B64 -v +x +P +a0.2 +r2 +q9 +j0.2');
     close(f);
end;


const light1:string[15]=('light_source{ <');
      light2:string[20]=(' color rgb<1,1,1>}');
      sphere:string[20]=(' sphere{< ');
      pigment:string[30]=('pigment{ color rgb<1,1,1>}}');


procedure  doc_povfile;
var f:text;
    t:integer;
    u,v,w:reel;
begin

     assign(f,povfile);
     rewrite(f);
     for  t:=1 to 6 do writeln(f,deb_povfile[t]);

     (* lumières *)
     writeln(f,light1,-trans_.x*7.5:4:4,
             ',',trans_.y*7.5:4:4,',',trans_.z*10.5:4:4,
             '> ',light2);
     writeln(f,light1,trans_.x*5:4:4,
             ',',trans_.y*7.5:4:4,',',-trans_.z*10.5:4:4,
             '> ',light2);
     writeln(f,light1, '100,',
             trans_.y*7.5:4:4,',',trans_.z*10.5:4:4,
             '> ',light2);

     writeln(f,' ');

     u:= ((o1.x+o2.x)/2)*6;
     v:= ((o1.y+o2.y)/2)*6;
     w:= ((o1.z+o2.z)/2)*8;

     x1:=u;
     z1:=w;

     u:=rotx(x1,z1,pi/2);
     w:=roty(x1,z1,pi/2);

     writeln(f,'camera{');
     writeln(f,'        location <',u:4:4,',',v:4:4,',',w:4:4,'>');

     u:= (o1.x+o2.x)/2;
     v:= (o1.y+o2.y)/2;
     w:= (o1.z+o2.z)/2;

     writeln(f,'//     direction      <0, 0, 1.2071>');
     writeln(f,'        up       <0, 1, 0>');
     writeln(f,'        right      <1.333, 0, 0>');
     writeln(f,'        look_at  <',u:4:4,',',v:4:4,',',w:4:4,'>');

     writeln(f,' }');

     writeln(f,' ');

     writeln(f,sphere, o1.x:4:6,',',o1.y:4:6,',',o1.z:4:6,'>,',o1.r:4:6,pigment);
     writeln(f,sphere, o2.x:4:6,',',o2.y:4:6,',',o2.z:4:6,'>,',o2.r:4:6,pigment);
     writeln(f,'cone{<0,0,0> , ',Pgray_c:4:6);
     writeln(f,        ' <0,',dimcone:4:6,',0> , ',ppray_c:4:6);
     writeln(f,'       rotate<',ang_zy*180/pi:4:6,',',ang_zx*180/pi:4:6,',0>');
     writeln(f,'       translate<',(x):4:6,',',(y):4:6,',',(z):4:6,'>');
     writeln(f,'     ',pigment);

     close(f);
end;

procedure tstpov;
begin
     translation_2;
     doc_povini;
     doc_povfile;
if paramstr(9)='POV' then
   begin
     swapvectors;
     exec(povpath+'povray.exe',povini);
     swapvectors;
   end;
end;

begin
     clrscr;
     test_param;
                 for n:=1 to 8 do
                     transforme_param(n,r3d[n]);

               (* De r3d vers o1 et o2 *)
                  transfert_donnees;

               (* A partir de o1 et o2, création
                  du cone *)
                  calcul_;

               (* transfers des données en texte
                  pour Povray *)
                  tstpov;
end.