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