Program Artilleur;


 { Epreuve informatique de l'Ecole Polytechnique                       94.232

 
 
 
1/	Soit  v  la vitesse absolue d'un projectile en mouvement dans l'air
et w celle du vent. On admet que le projectile est soumis à deux forces, la
pesanteur P=mg , et la résistance F de l'air, orientée comme w-v  et de
module égal à  |F|=k|v-w|² .

Ecrire un programme Pascal qui calcule pas à pas le mouvement du projectile
dans un repère fixe Oxyz, en supposant que z est l'altitude.
Montrer qu'en l'absence d'obstacle, la vitesse v tend vers une valeur
 asymptotique v infini  .


2 	Dans toute la suite, on supposera que le vecteur vent w  est
horizontal et constant, égal à w0  et que  le projectile est lancé à
l'instant t0  avec une vitesse initiale   du point M0 d'altitude
z0.  
 
Ecrire un programme qui calcule le point d'impact au sol, d'équation   et le temps mis pour
l'atteindre. 
 
 
3 	Pour M0  fixé, on suppose maintenant que   a un module constant, mais que son 
orientation est libre : tracer le contour du domaine représentatif de tous les points du sol qui peuvent 
être atteints. Comment varie ce contour en fonction du module de    ? Illustrer graphiquement. 
 
 
4 	On cherche maintenant à atteindre une cible animée d'un mouvement uniforme 
  dans l'espace. 
 
Le point M0  restant fixé, calculer une valeur de    permettant d'atteindre la cible. 
 
 
5 	Montrer que la solution à la question précédente n'est pas unique et écrire le 
programme Pascal permettant de trouver celle qui minimise le temps d'interception. 
 
Tracer graphiquement la trajectoire correspondante. 
 
       
 
NB : On pourra tester avec les valeurs métriques suivantes : 
 
	               
 
Imprimer tous les résultats             

en indiquant chaque fois à quoi ils correspondent  

 
-=-=-=-  }

Uses crt, modubase ;

Type


  {  Myreal = Double;  }
        Myreal = real;

    Temps  = Myreal;
    Distances = Myreal;
    Coordonnees = record
           x,y,z : Distances;
    end;

    Angle         = MyReal;
    Forces        = Coordonnees;
    Positions     = Coordonnees;
    Vitesses      = Coordonnees  ;
    Accelerations = Coordonnees;
    Vents         = Vitesses    ;
    Projectiles = record
               Position     : Positions;
               Vitesse      : Vitesses  ;
    end;
    Orientation = record
               Azimut,
               Site         : Angle;
    end;

var
   Mortier:Boolean;
   t,dt:Temps;
      xmin,xmax,
      ymin,ymax,
      zmin,zmax:real;
    vent : Vents;
    P,P1,P2:Projectiles;
    Initial   : projectiles;
    ksurm : Myreal;
    v0    : MyReal ;
   Ch:char;
   c:integer;
   pas : MyReal;
   xini : MyReal;
   xinipas : MyReal;
   choisir: Boolean;

   sizeY:MyReal;
Const
    pesanteur : Accelerations = (x:0;y:0;z:-9.81);
    M0:Positions=(x:0;y:0;z:5);
    C0:Positions=(x:10;y:0;z:10);
    P0:Projectiles=(position:(x:0;y:0;z:5);
                    vitesse:(x:0;y:10;z:10));
    vent0:Vents =(x:4;y:5;z:0);

Function vect(P,Q:Positions):real;
begin
     vect:=P.x*Q.y-P.y*Q.x;
end;

Function scalaire(P,Q:Positions):real;
begin
     scalaire:=P.x*Q.x+P.y*Q.y+P.z*Q.z;
end;

Function distance(P,Q:Positions):distances;
begin
     distance:=sqrt(sqr(P.x-Q.x)+sqr(P.y-Q.y)+sqr(P.z-Q.z));
end;

Function module(P:Positions):real;
begin
     module:=sqrt(scalaire(P,P));
end;


Function Afficher(P:Coordonnees):string;
var
   sx,sy,sz:string;
begin
     With P do
           begin
                Str(x:5:2,sx);
                Str(y:5:2,sy);
                Str(z:5:2,sz);
           end;
     Afficher:='(x:'+sx+', y:'+sy+', z:'+sz+')';
end;

Function Affiche_Coordonnees(P:Coordonnees):string;
var
   sx,sy,sz:string;
begin
     With P do
           begin
                Str(x:5:2,sx);
                Str(y:5:2,sy);
                Str(z:5:2,sz);
           end;
     Affiche_coordonnees:='('+sx+','+sy+','+sz+')';
end;

Function Affiche_Angle(a:angle):string;
var
   sa:string;
begin
                Str(round(a*180/pi):3,sa);
     Affiche_angle:=sa+'° ';
end;

Procedure Ajouter(VAR u:coordonnees;v,w:coordonnees);
begin
     u.x:=v.x+w.x;
     u.y:=v.y+w.y;
     u.z:=v.z+w.z;
end;

Procedure Homothetie(VAR u:coordonnees;k:Myreal;v:coordonnees);
begin
     u.x:=k*v.x;
     u.y:=k*v.y;
     u.z:=k*v.z;
end;

Procedure Soustraire(VAR u:coordonnees;v,w:coordonnees);
begin
     Homothetie(w,-1,w);
     Ajouter(u,v,w);
end;


Procedure Move(var p :Projectiles;dt:temps);
var
   du:positions;
   vitesse_relative,dv:vitesses ;
   Frottement : Forces;
   Acceleration : Accelerations;
begin
  With P do
       begin
            Soustraire(vitesse_relative,vent,vitesse);
            Homothetie(Frottement,ksurm*module(vitesse_relative),vitesse_relative);
            Ajouter(Acceleration,Pesanteur,Frottement);
            Homothetie(du,dt,P.vitesse);
            Ajouter(P.position,P.position,du);
            Homothetie(dv,dt,acceleration);
            Ajouter(P.vitesse,P.vitesse,dv);
       end;
end;
                                 
Procedure Pointer(VAR P:projectiles;O:orientation);
begin
     with O do
          With P.vitesse do
               begin
                    x:=v0*cos(site)*cos(azimut);
                    y:=v0*cos(site)*sin(azimut);
                    z:=v0*sin(site);
               end;
end;

Procedure Profil;
var
 P,P1,oldp:projectiles;
 O:orientation;

begin
     dt:=0.01;
     t:=0;
     P1.position:=M0;
     xmin:=-50;
     xmax:= 50;
     zmin:=-5;
     ModeGraphique;
     IsoFenetre(xmin,xmax,zmin);
     x_axe(0,0,10);
     y_axe(0,0,10);
     Couleur(Brillant);
     Deplace(xmin,-zmin*10);
     Ecris('Point de chute du projectile');
     Deplace(xmin,-zmin*9.5);
     Ecris('k/m=');Ecrisreel(ksurm);
     Deplace(xmin,-zmin*9.0);
     Ecris('v0=');Ecrisreel(v0);
     Deplace(xmin,-zmin*8.5);
     Ecris('M0  ='+Affiche_Coordonnees(P1.position));
     Deplace(xmin,-zmin*8.0);
     Ecris('vent='+Affiche_Coordonnees(vent));
     with o do
      begin
        azimut:=0;
        site:=pi/2;
        Repeat
           site:=site-0.1*sqrt(2);
           if site<-(pi/2) then
              begin
                   site:=site+pi;
                   azimut:=azimut+pi;
              end;
           azimut:=azimut+0.1;
           if azimut>(2*pi) then
              azimut:=azimut-2*pi;
           Couleur(-(1+Trunc(8*Random)));
           Deplace(xmin,-zmin*7.5);
           Ecris('Site ='+affiche_angle(Site));
           Deplace(xmin,-zmin*7.0);
           Ecris('Azimut ='+affiche_angle(Azimut));
             P:=P1;
           Pointer(P,O);
           With P.position do
            begin
                t:=0;
                repeat
                      t:=t+dt;
                      OldP:=P;
                      Move(P,dt);
                      with OldP.position do Deplace (x,z);
                      Trace(x,z);
                      with OldP.position do Deplace (x,40+y);
                      Trace(x,40+y);
                until P.position.z<=0;
                Deplace(xmin,-zmin*6.5);
                Ecris('t=');EcrisEntier(round(t));

            end; {  with P.position }
     until KeyPressed;
     Ch:=readKey;
     Pause;
   end; {with}
end;



Procedure Impact(VAR I:Positions;P:Projectiles;Cible:Positions);
begin
     t:=0;
     if mortier then
      begin
         dt:=0.01;
         repeat
           t:=t+dt;
           Move(P,dt);
         until (P.vitesse.z<0) and (P.position.z<Cible.z);
         I:=P.position;
         I.z:=Cible.z;
      end
     else
      begin
       dt:=0.01;
       I:=P.position;
       repeat
           t:=t+dt;
           Move(P,dt);
       until (P.position.x>Cible.x) or (P.position.z<0);
       If P.Position.x>Cible.x then
       begin
        I:=P.position;
        I.x:=Cible.x;
       end;
      end;
end;


Function Calculer_Impact(VAR I:positions;P:projectiles;O:orientation;Cible:positions):distances;
begin
     Pointer(P,O);
     Impact(I,P,Cible);
     Calculer_impact:=distance(I,Cible);
end;

Procedure Trajectoire(P0:Projectiles);
var
   w:Myreal;
   O:orientation;
   Q1,Q2:Positions;
   I,OldP:Positions;
   Cible:Positions;
   d:distances;

begin
     t:=0;
     dt:=0.005;
     Cible:=C0;
     xmin:=-50;
     xmax:=-xmin;
     ymin:=xmin*0.6;
     ModeGraphique;
     IsoFenetre(xmin,xmax,ymin);
     x_axe(0,0,10);
     y_axe(0,0,10);
     Couleur(Brillant);
     Deplace(xmin,-ymin);
     Ecris('Point de chute du projectile');
     Deplace(xmin,-ymin*0.95);
     Ecris('k/m=');Ecrisreel(ksurm);
     Deplace(xmin,-ymin*0.90);
     Ecris('v0=');Ecrisreel(v0);
     Deplace(xmin,-ymin*0.85);
     Ecris('M0  ='+Affiche_Coordonnees(P0.position));
     Deplace(xmin,-ymin*0.80);
     Ecris('vent='+Affiche_Coordonnees(vent));

    with Q1 do
      begin
           x:=xmin/2;
           y:=-ymin/2;
           z:=0;
      end;
      Ajouter(Q2,Q1,vent);
      With Q1 do deplace(x,y);
      With Q2 do trace (x,y);
      Ecris('Vent');
    with O do
     begin
      site:=-Pi/2;
      Repeat
           site:=site+0.1;
           Couleur(-(1+Trunc(8*Random)));
           Deplace(xmin,-ymin*0.75);
           Ecris('Site ='+Affiche_angle(Site));
           azimut:=0;
           repeat
           azimut:=azimut+0.02;
             P:=P0;
            Pointer(P,O);
            With P.position do
             begin
                t:=0;
                repeat
                      t:=t+dt;
                      Move(P,dt);
                until P.position.z<=0;
              With OldP do
                Deplace(x,y);
              If azimut>0.03 then
                 With P do
                Trace(x,y);
              OldP:=P.position;
           
           Deplace(xmin,-ymin*0.70);
           Ecris('t=');Ecrisentier(Round(t));

            end; {  with P.position }
           until azimut>2*pi;
          until KeyPressed;
         Ch:=readKey;
      Pause;
     end; { with O }
end;



Function Calculer_tir(VAR O:orientation;Var d:Distances;v0:real;M0,cible:coordonnees):Temps;
Var
   n:integer;
   DeltaO:orientation;

Function Rapprocher(Var O,deltaO:orientation):distances;

{  On va chercher O=(azimut,site) par approximations successives

Soient y et z les coordonnées du projectiles lorsqu'il passe dans le
plan  x=Cible.x

On a l'équation différentielle :

       dy         dy
  dy= ---- da +  ---- ds
       da         ds

       dz         dz
  dz= ---- da +  ---- ds
       da         ds

       En inversant le système, on obtient les paramètres d'ajustement

                 dy     dz      dy     dz
  determinant = ---- * ----  - ---- * ----
                 da     ds      ds     da

        dz         dy
  da=[ ---- dy -  ----- dz ]/determinant
        ds         ds

         dz         dy
  ds=[- ---- dy +  ----- dz ]/determinant
         da         da
            }

var
   Oa,Os,Oi,dOi:orientation;
   I,Ia,Is,Ii:positions;
   dx_da,dx_ds,
   dy_da,dy_ds,
   dz_da,dz_ds,
   determinant : real;
   Ecart_x,
   Ecart_y,
   Ecart_z,
   dist,dista,disti,dists,
   d : distances;
begin    {  rapprocher }
     Os:=O; with Os do site  :=site  +deltaO.site  ;
     Oa:=O; with Oa do azimut:=azimut+deltaO.azimut;
   { with deltaO do WriteLN('ds=',site:7:4,' da=',azimut:7:4);}
     dist:=Calculer_Impact(I ,P0,O ,cible);
   { WriteLN('Impact_I ='+affiche_coordonnees(I)+' dist=',dist:8:3);}
     dists:=Calculer_Impact(Is,P0,Os,cible);
   { WriteLN('Impact_Is ='+affiche_coordonnees(Is)+' dist=',dists:8:3);}
     dista:=Calculer_Impact(Ia,P0,Oa,cible);
   { WriteLN('Impact_Ia ='+affiche_coordonnees(Ia)+' dist=',dista:8:3); }
     d:=dist;
     if dists<dist then
        begin
             dist:=dists;
             O:=Os
        end;
     if dista<dist then
        begin
             dist:=dista;
             O:=Oa;
        end;
    If d<=dist then
        begin
              if mortier then
                 begin
                      Ecart_x:=Cible.x-I.x;
                      Ecart_y:=Cible.y-I.y;
                      with deltaO do
                           begin
                                dx_da:=(Ia.x-I.x)/azimut;
                                dx_ds:=(Is.x-I.x)/site;
                                dy_da:=(Ia.y-I.y)/azimut;
                                dy_ds:=(Is.y-I.y)/site;
                           end;
                      determinant:=dx_da*dy_ds-dx_ds*dy_da;
                   {    WriteLN('determinant=',determinant:9:4); }
                      dOi:=deltaO;
                      with dOi do
                           if determinant<>0 then
                              begin
                                   azimut:=(dy_ds*Ecart_x-dx_ds*Ecart_y)/Determinant;
                                   site  :=(dx_da*Ecart_y-dy_da*Ecart_x)/Determinant;
                                   Oi.azimut:=O.azimut+azimut;
                                   Oi.site  :=O.site  +site;
                              end;
                  end  {if mortier}
              else
                  begin
                       Ecart_y:=Cible.y-I.y;
                       Ecart_z:=Cible.z-I.z;
                       with deltaO do
                            begin
                                 dy_da:=(Ia.y-I.y)/azimut;
                                 dy_ds:=(Is.y-I.y)/site;
                                 dz_da:=(Ia.z-I.z)/azimut;
                                 dz_ds:=(Is.z-I.z)/site;
                            end;
                       determinant:=dy_da*dz_ds-dy_ds*dz_da;
                       dOi:=deltaO;
                       with dOi do
                            if determinant<>0 then
                               begin
                                    azimut:=(dz_ds*Ecart_y-dy_ds*Ecart_z)/Determinant;
                                    site  :=(dy_da*Ecart_z-dz_da*Ecart_y)/Determinant;
                                    Oi.azimut:=O.azimut+azimut;
                                    Oi.site  :=O.site  +site;
                               end;
                  end; {else}
              disti:=Calculer_Impact(Ii,P0,Oi,cible);
              If disti<dist then
                 begin
                   dist:=disti;
                   O:=Oi;
                   deltaO:=dOi
                 end
              else
                  with deltaO do
                       begin
                            azimut:=-azimut/2;
                            site:=-site/2;
                       end;
        end;   {else}
     Rapprocher:=dist;
end;

Procedure Tracer_trajectoire(P:projectiles);
Const
     dt=0.01;
VAr
        OldP:Projectiles;
begin
     With P.position do
       begin
          Deplace(x,y);
          repeat
          OldP:=P;
           Move(P,dt);
           With OldP.position do deplace(y,z);
           Trace(y,z);
           until (z<0) or (P.position.x>Cible.x);
        end; {with}

end;

begin                {  Calculer_tir }
   P0.position:=M0;
   with O do
   begin
     n:=0;
    with deltaO do
      begin
           if mortier then
              begin
                   site:=- 0.5;
                   azimut:=1;
              end
           else
             begin
               site:=0.2;
               azimut:=1;
             end;
      end;
   repeat
     Inc(n);
{     Write('.');}
     d:=Rapprocher(O,deltaO);
     with O do
      begin
       if site>pi/2 then
        begin
            site:=pi-site;
            azimut:=azimut+pi;
        end;
       while azimut>2*pi do
             azimut:=azimut-2*pi;
       while azimut<0 do
             azimut:=azimut+2*pi;
      end; {  with }
{      if mortier then
         Write('Tir mortier n°')
         else
         Write('Tir tendu n°');
      WriteLN(n,' site=',Affiche_angle(site),
      ' azimut=',Affiche_angle(azimut),
      ' dist=',d:6:3,' durée=',t:5:2); }
   until (n>30) or (d<0.1);
   if d>=0.1 then    { si le tir a échoué }
         t:=0;
    Calculer_tir:=t;
   
 end ; { with o}
end;


Procedure Presentation;
begin;
      Modetexte;
      Efface;
      WriteLN(' ====== Projectile d''artillerie  ============');
      WriteLN('                                        ');
      WriteLN('                                        ');
      WriteLN('   (1) Profil de la trajectoire   ');
      WriteLN('   (2) Domaine utile en plan   ');
      WriteLN('   (3) Cible fixe   ');
      WriteLN('   (4) Cible mobile  ');
      WriteLN('                              ');
      Write('   Tapez votre Choix  :  ');
end;





procedure Question1;
begin
     Profil;
end;

Procedure Question2;
begin
  Trajectoire(P0);
end;

procedure Question3;
Var
   M0,Cible,Impact:Positions;
   i,j:integer;
   O:orientation;
   duree:temps;

procedure Resoudre(s,a:angle;m:Boolean);
var
   d:distances;
begin
   With O do
        begin
             site:=s;
             azimut:=a;
        end;
   mortier:=m;
   duree:=Calculer_tir(O,d,v0,M0,Cible);
   if mortier then
      Write('Tir mortier : ')
   else
       Write('Tir tendu : ');
   With O do
        WriteLN(' site='  ,Affiche_angle(site),
                ' azimut=',Affiche_angle(azimut),
                ' dist=',d:8:3,' durée=',duree:5:2);
        
end;

begin
     Efface;
     WriteLN('Pointage sur une cible fixe');
     WriteLN('k/m=',ksurm:8:3);
     WriteLN('v0=',v0:5:0);
     M0:=P0.position;
     WriteLN('M0  ='+Affiche_Coordonnees(M0));
     WriteLN('vent='+Affiche_Coordonnees(vent));
     Cible:=C0;
     WriteLN('Cible='+Affiche_coordonnees(Cible));
     with cible do
       begin
          Resoudre(pi/2-0.2,arctan((y-M0.y)/(x-M0.x)),true);
          Resoudre(0,arctan((y-M0.y)/(x-M0.x)),false);
       end;
     Pause;
end;


procedure Question4;
Type
    Triplet = record
           t0,tc,duree: Temps;
    end;
Var
   Cible,Impact:Positions;
   i,j:integer;
   O:orientation;
   tr,tr_old:triplet;


procedure Resoudre(t:temps);
var
   d:distances;
begin
   Couleur(Blanc);
   With Cible do Croix(y,10);
   with tr do
    begin
         tc:=t;
         {WriteLN('cible='+affiche_coordonnees(Cible));}
         duree:=Calculer_tir(O,d,v0,M0,Cible);
         Deplace(xmin,ymax*0.70);
         Ecris('cible='+affiche_coordonnees(Cible));
         Deplace(xmin,ymax*0.65);
         If mortier then
            Ecris('Tir mortier : ')
         else
             Ecris('Tir tendu : ');
         With O do
              Ecris(' site='  +Affiche_angle(site)+
                ' azimut='+Affiche_angle(azimut));
         t0:=tc-duree;
           Couleur(vert) ;Ecris('t0=');EcrisReel(t0);
           Couleur(blanc);Ecris(' tc ');
           Couleur(jaune);Ecris('durée=');EcrisReel(duree);
         if (t<>0) and (duree>0) then
          begin
           Couleur(vert) ;with tr_old do Deplace(tc,t0);Trace(tc,t0);
           Couleur(blanc);with tr_old do Deplace(tc,tc);Trace(tc,tc);
           Couleur(jaune);with tr_old do Deplace(tc,duree);Trace(tc,duree);
          end;
         tr_old:=tr;
   {    if mortier then
          Write('Tir mortier : ')
       else
           Write('Tir tendu : ');
   With O do
        WriteLN(' site='  ,Affiche_angle(site),
                ' azimut=',Affiche_angle(azimut),
                ' dist=',d:8:3,' durée=',duree:5:2);}

    end; { with tr }
end;

begin
     xmin:=-40;
     xmax:=-xmin;
     ymin:=xmin*0.6;
     ymax:=-xmin;
     ModeGraphique;
     Fenetre(xmin,xmax,ymin,ymax);
     x_axe(0,0,10);
     y_axe(0,0,10);
     Couleur(Brillant);
     Deplace(xmin,ymax*0.95);
     Ecris('Tir sur une cible mobile');
     Deplace(xmin,ymax*0.90);
     Ecris('k/m=');Ecrisreel(ksurm);
     Deplace(xmin,ymax*0.85);
     Ecris('v0=');Ecrisreel(v0);
     Deplace(xmin,ymax*0.80);
     Ecris('M0  ='+Affiche_Coordonnees(P0.position));
     Deplace(xmin,ymax*0.75);
     Ecris('vent='+Affiche_Coordonnees(vent));

   {  WriteLN('Pointage sur une cible mobile');
     WriteLN('k/m=',ksurm:8:3);
     WriteLN('v0=',v0:5:0);
     WriteLN('M0  ='+Affiche_Coordonnees(P0.position));
     WriteLN('vent='+Affiche_Coordonnees(vent));}
     with O do
      begin
           Cible:=C0;
           site:=0;
           azimut:=0;
           repeat
                 Mortier:=false;
                 Resoudre(Cible.y);
                 Cible.y:=Cible.y+1;
           until Cible.y>35;
           Cible:=C0;
           site:=pi/2-0.1;
           azimut:=0;
           repeat
                 Mortier:=true;
                 Resoudre(Cible.y);
                 Cible.y:=Cible.y+1;
           until Cible.y>35;
           Cible:=C0;
           site:=0;
           azimut:=0;
           repeat
                 Mortier:=false;
                 Resoudre(Cible.y);
                 Cible.y:=Cible.y-1;
           until Cible.y<-30;
           Cible:=C0;
           site:=pi/2-0.1;
           azimut:=0;
             repeat
                 Mortier:=true;
                 Resoudre(Cible.y);
                 Cible.y:=Cible.y-1;
           until Cible.y<-30;
           Pause;
      end;

end;



begin
  InitGraphique;
  Randomize;
     repeat
          Presentation;
              Read(Ch);
          ksurm:=0.05;
          v0:=30;
          vent:=vent0;
          case ch of
               '1':  Question1;
               '2':  Question2;
               '3' : Question3;
               '4' : Question4;
               '0' : Halt
          end;
    until false;



end.