+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