Premier chargement des fichiers
[mp-geo.git] / doc / mp-geotrans.mp
diff --git a/doc/mp-geotrans.mp b/doc/mp-geotrans.mp
new file mode 100644 (file)
index 0000000..11ba98d
--- /dev/null
@@ -0,0 +1,248 @@
+prologues:=2;
+
+input geometriesyr16;
+
+color tan,payscolor,cielfonce;
+tan=(0.824,0.705,0.55);
+payscolor=tan;
+cielfonce=(0.25,0.75,1);
+
+vardef Projgeo(expr X)=
+  pair $;
+  numeric Xobs,Yobs,Zobs,XProj,YProj;
+  Xobs = -redpart(X)*Aux1 + greenpart(X)*Aux3;
+  Yobs = -redpart(X)*Aux5 - greenpart(X)*Aux6 + bluepart(X)*Aux4;
+  Zobs = -redpart(X)*Aux7 - greenpart(X)*Aux8 - bluepart(X)*Aux2 + Rho;
+  XProj = DE*Xobs/Zobs;
+  YProj = DE*Yobs/Zobs;
+  $=(XProj,YProj);
+  $
+enddef;
+
+vardef FX(expr t,v)=rayon*cosd(c*t)*cosd(c*v)
+enddef;
+
+vardef FY(expr t,v)=rayon*cosd(c*t)*sind(c*v)
+enddef;
+
+vardef FZ(expr t,v)=rayon*sind(c*t)
+enddef;
+
+vardef arcsind(expr x)=%Définition mathématique en degré ici :)
+  angle((sqrt(1-x**2),x))
+enddef;
+
+numeric nbpts,nblec,nbcapitales;
+
+string arborescence;
+arborescence:="../DATA/";
+
+path pays[];path paysnonvus[];
+
+vardef lecture(expr nomfichier,fond)=
+  color Coord[],Pays[],Paysnonvus[];
+  numeric ll,lnv;
+  ll:=0;lnv:=0;
+  nbpts:=scantokens readfrom nomfichier;
+  for k=1 upto (nbpts div 3):
+    pair latlon;
+    latlon=((scantokens readfrom nomfichier)+(scantokens readfrom nomfichier)+(scantokens readfrom nomfichier))/3;
+    Coord[k]=rayon*(cosd(xpart(latlon/60))*cosd(ypart(latlon/60)),cosd(xpart(latlon/60))*sind(ypart(latlon/60)),sind(xpart(latlon/60)));
+    if ProduitScalaire(Coord[k]-pte3,Oeil-pte3)>0:
+      if ((xpart(latlon/60)>phim) and (xpart(latlon/60)<phip)):
+       ll:=ll+1;
+       Pays[k]=Coord[k];
+       Paysnonvus[k]=2*Coord[k];
+      else:
+       Pays[k]=2*Coord[k];
+       Paysnonvus[k]=2*Coord[k];
+      fi;
+    else:
+      Pays[k]=2*Coord[k];
+      lnv:=lnv+1;
+      Paysnonvus[k]=Coord[k];
+    fi;
+  endfor;
+  closefrom nomfichier;
+  if ll>0:
+    nbvus:=nbvus+1;
+    pays[nbvus]=Projgeo(Pays[1])
+    for l=2 upto (nbpts div 3):
+      --Projgeo(Pays[l])
+    endfor;
+  else:
+    nbnonvus:=nbnonvus+1;
+    %show nbnonvus;
+    paysnonvus[nbnonvus]=Projgeo(Paysnonvus[1])
+    for l=2 upto (nbpts div 3):
+      --Projgeo(Paysnonvus[l])
+    endfor;
+  fi;
+  if lnv>0:
+    nbnonvus:=nbnonvus+1;
+    %show nbnonvus;
+    paysnonvus[nbnonvus]=Projgeo(Paysnonvus[1])
+    for l=2 upto (nbpts div 3):
+      --Projgeo(Paysnonvus[l])
+    endfor;
+  fi;
+enddef;
+
+string NomFichier;
+
+vardef Lecture(expr Nomfichier)=
+  NomFichier:=arborescence&Nomfichier;
+  nblec:=scantokens readfrom NomFichier;
+  for w=1 upto nblec:
+    lecture(scantokens readfrom NomFichier);
+  endfor;
+  closefrom NomFichier;
+enddef;
+
+string nomfichiermul;
+
+vardef Lectureiles=
+  nomfichiermul:=arborescence&"iles.dat";
+  nblec:=scantokens readfrom nomfichiermul;
+  for p=1 upto nblec:
+    color Coord[],fond,Pays[],Paysnonvus[];
+    numeric ll,lnv;
+    ll:=0;lnv:=0;
+    nbpts:=scantokens readfrom nomfichiermul;
+    fond=scantokens readfrom nomfichiermul;
+    for k=1 upto nbpts:
+      pair latlon;
+      latlon=scantokens readfrom nomfichiermul;
+      Coord[k]=rayon*(cosd(xpart(latlon/60))*cosd(ypart(latlon/60)),cosd(xpart(latlon/60))*sind(ypart(latlon/60)),sind(xpart(latlon/60)));
+      if ProduitScalaire(Coord[k]-pte3,Oeil-pte3)>0:
+       if ((xpart(latlon/60)>phim) and (xpart(latlon/60)<phip)):
+         ll:=ll+1;
+         Pays[k]=Coord[k];
+         Paysnonvus[k]=2*Coord[k];
+       else:
+         Pays[k]=2*Coord[k];
+         Paysnonvus[k]=2*Coord[k];
+       fi;
+      else:
+       Pays[k]=2*Coord[k];
+       lnv:=lnv+1;
+       Paysnonvus[k]=Coord[k];
+      fi;
+    endfor;
+    if ll>0:
+      nbvus:=nbvus+1;
+      pays[nbvus]=Projgeo(Pays[1])
+      for l=2 upto nbpts:
+       --Projgeo(Pays[l])
+      endfor;
+    else:
+      nbnonvus:=nbnonvus+1;
+      %show nbnonvus;
+      paysnonvus[nbnonvus]=Projgeo(Paysnonvus[1])
+      for l=2 upto nbpts:
+       --Projgeo(Paysnonvus[l])
+      endfor;
+    fi;
+    if lnv>0:
+      nbnonvus:=nbnonvus+1;
+      %show nbnonvus;
+      paysnonvus[nbnonvus]=Projgeo(Paysnonvus[1])
+      for l=2 upto nbpts:
+       --Projgeo(Paysnonvus[l])
+      endfor;
+    fi;
+  endfor;
+  closefrom nomfichiermul;
+enddef;
+
+vardef Recap=
+  %message("nombre de pays non vus (partiellement ou total)");
+  %show nbnonvus;
+  %message("nombre de pays vus (partiellement ou total)");
+  %show nbvus;
+  for k=1 upto nbnonvus:
+    remplis paysnonvus[k]--cycle withcolor tan;
+    trace paysnonvus[k];
+  endfor;
+  fillcolor:=ciel;
+  transparence cercles(pte3,pte1,pte3,pte1,pte4);
+  for k=1 upto nbvus:
+    remplis pays[k]--cycle withcolor tan;
+    trace pays[k];
+  endfor;
+  clip currentpicture to cercles(pte3,pte1,pte3,pte1,pte4);
+enddef;
+
+nbvus:=0;
+nbnonvus:=0;
+
+rayon:=2;
+
+vardef MaillageS=
+  path Maillage[];
+  color CentMaillage[];
+  for k=Udebut step 0.5*pasU until (Ufin):
+    for l=Vdebut step pasV until (Vfin):
+      Maillage[100*k+l]=Projgeo((FX(k,l),FY(k,l),FZ(k,l)))--Projgeo((FX(k,l+pasV),FY(k,l+pasV),FZ(k,l+pasV)))--Projgeo((FX(k+pasU,l+pasV),FY(k+pasU,l+pasV),FZ(k+pasU,l+pasV)))--Projgeo((FX(k+pasU,l),FY(k+pasU,l),FZ(k+pasU,l)))--cycle;
+      if ProduitScalaire(1/2[(FX(k,l),FY(k,l),FZ(k,l)),(FX(k+pasU,l+pasV),FY(k+pasU,l+pasV),FZ(k+pasU,l+pasV))]-pte3,Oeil-pte3)>0:
+       draw Maillage[100*k+l];
+      fi;
+    endfor;
+  endfor;
+enddef;
+
+vardef InitialiseMaillage(expr ud,uf,up,vd,vf,vp)=
+  Udebut:=ud;
+  Ufin:=uf;
+  pasU:=up;
+  Vdebut:=vd;
+  Vfin:=vf;
+  pasV:=vp;
+enddef;
+
+vardef MaillageSphere=
+  InitialiseMaillage(((phim div 10)+1)*pi/18,((phip div 10)-1)*pi/18,pi/18,-pi,pi,pi/36);
+  MaillageS;
+enddef;
+
+vardef Mappemonde(expr longobs,latobs)=
+  figureespace(-100u,-100u,100u,100u);
+  Initialisation(5,longobs,latobs,750);
+  numeric phim,phip,phii;%phi moins -- phi plus - phi intermédiaire
+  phim=Phi+arcsind(rayon/Rho)-90;
+  phip=Phi+90-arcsind(rayon/Rho);
+  color pte[];
+  pte1=rayon*(cosd(phim)*cosd(Theta),cosd(phim)*sind(Theta),sind(phim));
+  pte2=rayon*(cosd(phip)*cosd(Theta),cosd(phip)*sind(Theta),sind(phip));
+  pte3=iso(pte1,pte2);
+  pte4-pte3=Normal((0,0,0),pte1,pte2);
+  if (Phi>90):
+    phip:=180-phip;
+    phii:=180-phim;
+    phim:=phip;
+    phip:=phii;
+  fi;
+  if (Phi<-90):
+    phip:=-180-phip;
+    phii:=-180-phim;
+    phim:=phip;
+    phip:=phii;
+  fi;
+  Lecture("Cameriquesud.dat");
+  Lecture("Cameriquecentrale.dat");
+  Lecture("Ccaraibes.dat");
+  Lecture("Cameriquenord.dat");
+  Lecture("Casie.dat");
+  Lecture("Ceurope.dat");
+  Lecture("Cafrique.dat");
+  %Lecturelacs;
+  Lectureiles;
+  %Lecturerivieres;
+  Recap;
+  trace cercles(pte3,pte1,pte3,pte1,pte4);
+  finespace;
+enddef;
+
+Mappemonde(-35,50);
+
+end
\ No newline at end of file

Licence Creative Commons Les fichiers de Syracuse sont mis à disposition (sauf mention contraire) selon les termes de la
Licence Creative Commons Attribution - Pas d’Utilisation Commerciale - Partage dans les Mêmes Conditions 4.0 International.