doc/.svn/text-base/mp-geotrans.mp.svn-base

   1 prologues:=2;
   2
   3 input geometriesyr16;
   4
   5 color tan,payscolor,cielfonce;
   6 tan=(0.824,0.705,0.55);
   7 payscolor=tan;
   8 cielfonce=(0.25,0.75,1);
   9
  10 vardef Projgeo(expr X)=
  11   pair $;
  12   numeric Xobs,Yobs,Zobs,XProj,YProj;
  13   Xobs = -redpart(X)*Aux1 + greenpart(X)*Aux3;
  14   Yobs = -redpart(X)*Aux5 - greenpart(X)*Aux6 + bluepart(X)*Aux4;
  15   Zobs = -redpart(X)*Aux7 - greenpart(X)*Aux8 - bluepart(X)*Aux2 + Rho;
  16   XProj = DE*Xobs/Zobs;
  17   YProj = DE*Yobs/Zobs;
  18   $=(XProj,YProj);
  19   $
  20 enddef;
  21
  22 vardef FX(expr t,v)=rayon*cosd(c*t)*cosd(c*v)
  23 enddef;
  24
  25 vardef FY(expr t,v)=rayon*cosd(c*t)*sind(c*v)
  26 enddef;
  27
  28 vardef FZ(expr t,v)=rayon*sind(c*t)
  29 enddef;
  30
  31 vardef arcsind(expr x)=%Définition mathématique en degré ici :)
  32   angle((sqrt(1-x**2),x))
  33 enddef;
  34
  35 numeric nbpts,nblec,nbcapitales;
  36
  37 string arborescence;
  38 arborescence:="../DATA/";
  39
  40 path pays[];path paysnonvus[];
  41
  42 vardef lecture(expr nomfichier,fond)=
  43   color Coord[],Pays[],Paysnonvus[];
  44   numeric ll,lnv;
  45   ll:=0;lnv:=0;
  46   nbpts:=scantokens readfrom nomfichier;
  47   for k=1 upto (nbpts div 3):
  48     pair latlon;
  49     latlon=((scantokens readfrom nomfichier)+(scantokens readfrom nomfichier)+(scantokens readfrom nomfichier))/3;
  50     Coord[k]=rayon*(cosd(xpart(latlon/60))*cosd(ypart(latlon/60)),cosd(xpart(latlon/60))*sind(ypart(latlon/60)),sind(xpart(latlon/60)));
  51     if ProduitScalaire(Coord[k]-pte3,Oeil-pte3)>0:
  52       if ((xpart(latlon/60)>phim) and (xpart(latlon/60)<phip)):
  53         ll:=ll+1;
  54         Pays[k]=Coord[k];
  55         Paysnonvus[k]=2*Coord[k];
  56       else:
  57         Pays[k]=2*Coord[k];
  58         Paysnonvus[k]=2*Coord[k];
  59       fi;
  60     else:
  61       Pays[k]=2*Coord[k];
  62       lnv:=lnv+1;
  63       Paysnonvus[k]=Coord[k];
  64     fi;
  65   endfor;
  66   closefrom nomfichier;
  67   if ll>0:
  68     nbvus:=nbvus+1;
  69     pays[nbvus]=Projgeo(Pays[1])
  70     for l=2 upto (nbpts div 3):
  71       --Projgeo(Pays[l])
  72     endfor;
  73   else:
  74     nbnonvus:=nbnonvus+1;
  75     %show nbnonvus;
  76     paysnonvus[nbnonvus]=Projgeo(Paysnonvus[1])
  77     for l=2 upto (nbpts div 3):
  78       --Projgeo(Paysnonvus[l])
  79     endfor;
  80   fi;
  81   if lnv>0:
  82     nbnonvus:=nbnonvus+1;
  83     %show nbnonvus;
  84     paysnonvus[nbnonvus]=Projgeo(Paysnonvus[1])
  85     for l=2 upto (nbpts div 3):
  86       --Projgeo(Paysnonvus[l])
  87     endfor;
  88   fi;
  89 enddef;
  90
  91 string NomFichier;
  92
  93 vardef Lecture(expr Nomfichier)=
  94   NomFichier:=arborescence&Nomfichier;
  95   nblec:=scantokens readfrom NomFichier;
  96   for w=1 upto nblec:
  97     lecture(scantokens readfrom NomFichier);
  98   endfor;
  99   closefrom NomFichier;
 100 enddef;
 101
 102 string nomfichiermul;
 103
 104 vardef Lectureiles=
 105   nomfichiermul:=arborescence&"iles.dat";
 106   nblec:=scantokens readfrom nomfichiermul;
 107   for p=1 upto nblec:
 108     color Coord[],fond,Pays[],Paysnonvus[];
 109     numeric ll,lnv;
 110     ll:=0;lnv:=0;
 111     nbpts:=scantokens readfrom nomfichiermul;
 112     fond=scantokens readfrom nomfichiermul;
 113     for k=1 upto nbpts:
 114       pair latlon;
 115       latlon=scantokens readfrom nomfichiermul;
 116       Coord[k]=rayon*(cosd(xpart(latlon/60))*cosd(ypart(latlon/60)),cosd(xpart(latlon/60))*sind(ypart(latlon/60)),sind(xpart(latlon/60)));
 117       if ProduitScalaire(Coord[k]-pte3,Oeil-pte3)>0:
 118         if ((xpart(latlon/60)>phim) and (xpart(latlon/60)<phip)):
 119           ll:=ll+1;
 120           Pays[k]=Coord[k];
 121           Paysnonvus[k]=2*Coord[k];
 122         else:
 123           Pays[k]=2*Coord[k];
 124           Paysnonvus[k]=2*Coord[k];
 125         fi;
 126       else:
 127         Pays[k]=2*Coord[k];
 128         lnv:=lnv+1;
 129         Paysnonvus[k]=Coord[k];
 130       fi;
 131     endfor;
 132     if ll>0:
 133       nbvus:=nbvus+1;
 134       pays[nbvus]=Projgeo(Pays[1])
 135       for l=2 upto nbpts:
 136         --Projgeo(Pays[l])
 137       endfor;
 138     else:
 139       nbnonvus:=nbnonvus+1;
 140       %show nbnonvus;
 141       paysnonvus[nbnonvus]=Projgeo(Paysnonvus[1])
 142       for l=2 upto nbpts:
 143         --Projgeo(Paysnonvus[l])
 144       endfor;
 145     fi;
 146     if lnv>0:
 147       nbnonvus:=nbnonvus+1;
 148       %show nbnonvus;
 149       paysnonvus[nbnonvus]=Projgeo(Paysnonvus[1])
 150       for l=2 upto nbpts:
 151         --Projgeo(Paysnonvus[l])
 152       endfor;
 153     fi;
 154   endfor;
 155   closefrom nomfichiermul;
 156 enddef;
 157
 158 vardef Recap=
 159   %message("nombre de pays non vus (partiellement ou total)");
 160   %show nbnonvus;
 161   %message("nombre de pays vus (partiellement ou total)");
 162   %show nbvus;
 163   for k=1 upto nbnonvus:
 164     remplis paysnonvus[k]--cycle withcolor tan;
 165     trace paysnonvus[k];
 166   endfor;
 167   fillcolor:=ciel;
 168   transparence cercles(pte3,pte1,pte3,pte1,pte4);
 169   for k=1 upto nbvus:
 170     remplis pays[k]--cycle withcolor tan;
 171     trace pays[k];
 172   endfor;
 173   clip currentpicture to cercles(pte3,pte1,pte3,pte1,pte4);
 174 enddef;
 175
 176 nbvus:=0;
 177 nbnonvus:=0;
 178
 179 rayon:=2;
 180
 181 vardef MaillageS=
 182   path Maillage[];
 183   color CentMaillage[];
 184   for k=Udebut step 0.5*pasU until (Ufin):
 185     for l=Vdebut step pasV until (Vfin):
 186       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;
 187       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:
 188         draw Maillage[100*k+l];
 189       fi;
 190     endfor;
 191   endfor;
 192 enddef;
 193
 194 vardef InitialiseMaillage(expr ud,uf,up,vd,vf,vp)=
 195   Udebut:=ud;
 196   Ufin:=uf;
 197   pasU:=up;
 198   Vdebut:=vd;
 199   Vfin:=vf;
 200   pasV:=vp;
 201 enddef;
 202
 203 vardef MaillageSphere=
 204   InitialiseMaillage(((phim div 10)+1)*pi/18,((phip div 10)-1)*pi/18,pi/18,-pi,pi,pi/36);
 205   MaillageS;
 206 enddef;
 207
 208 vardef Mappemonde(expr longobs,latobs)=
 209   figureespace(-100u,-100u,100u,100u);
 210   Initialisation(5,longobs,latobs,750);
 211   numeric phim,phip,phii;%phi moins -- phi plus - phi intermédiaire
 212   phim=Phi+arcsind(rayon/Rho)-90;
 213   phip=Phi+90-arcsind(rayon/Rho);
 214   color pte[];
 215   pte1=rayon*(cosd(phim)*cosd(Theta),cosd(phim)*sind(Theta),sind(phim));
 216   pte2=rayon*(cosd(phip)*cosd(Theta),cosd(phip)*sind(Theta),sind(phip));
 217   pte3=iso(pte1,pte2);
 218   pte4-pte3=Normal((0,0,0),pte1,pte2);
 219   if (Phi>90):
 220     phip:=180-phip;
 221     phii:=180-phim;
 222     phim:=phip;
 223     phip:=phii;
 224   fi;
 225   if (Phi<-90):
 226     phip:=-180-phip;
 227     phii:=-180-phim;
 228     phim:=phip;
 229     phip:=phii;
 230   fi;
 231   Lecture("Cameriquesud.dat");
 232   Lecture("Cameriquecentrale.dat");
 233   Lecture("Ccaraibes.dat");
 234   Lecture("Cameriquenord.dat");
 235   Lecture("Casie.dat");
 236   Lecture("Ceurope.dat");
 237   Lecture("Cafrique.dat");
 238   %Lecturelacs;
 239   Lectureiles;
 240   %Lecturerivieres;
 241   Recap;
 242   trace cercles(pte3,pte1,pte3,pte1,pte4);
 243   finespace;
 244 enddef;
 245
 246 Mappemonde(-35,50);
 247
 248 end