%j'ai du modifier main_memory
%sous Debian Etch
%modifier /etc/texmf/texmf.d/95NonPath.cnf
%main_memory:1500000->40000000 !!!!
%puis update-texmf
%puis fmtutil-sys --refresh
% ->not ok :( "independent variables" pose problème.
% solution: réduction du nombre de points placés.
% pb suivant : lorsque l'on a une position de l'observateur avec un
% phi>90, rien ne s'affiche->réglé

% autre solution envisagée ici : lecture des données dans un fichier
% extérieur au format dat par exemple

% il faut changer aussi le path size :) >2000 pour une version avec ..
% au lieu de -- mais cette version marche pas :(
% modifier /etc/texmf/texmf.d/95NonPath.cnf
% path_size.mpost=10000

% readfrom files=30. Comment changer ? Pas la peine : inclusion d'un
% closefrom (présent dans metapost depuis 0.64) pour clotûrer vraiment
% chaque fichier en fin de lecture. quand readfrom lit une ligne d'un
% fichier, il ne passe pas de suite à la prochaine ligne, il attend
% simplement le prochain ordre "readfrom".

% ajout de LectureMultiple pour les iles (utilisable pour les lacs :
% dépassement metapost ?)

prologues:=2;

input geometriesyr16;

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)=ray*cosd(c*t)*cosd(c*v)
enddef;

vardef FY(expr t,v)=ray*cosd(c*t)*sind(c*v)
enddef;

vardef FZ(expr t,v)=ray*sind(c*t)
enddef;

vardef arcsind(expr x)=%Définition mathématique en degré ici :)
  angle((sqrt(1-x**2),x))
enddef;

numeric nbpts;

vardef lecture(expr nomfichier,fond)=
  color Coord[];
  numeric ll;
  ll:=0;
  color Pays[];
  nbpts:=scantokens readfrom nomfichier;
  for k=1 upto nbpts:
    pair latlon;
    latlon=scantokens readfrom nomfichier;
    Coord[k]=ray*(cosd(xpart(latlon/60))*cosd(ypart(latlon/60)),
	cosd(xpart(latlon/60))*sind(ypart(latlon/60)),sind(xpart(latlon/60)));
    if ((xpart(latlon/60)>phim) and (xpart(latlon/60)<phip)):
      if ProduitScalaire(Coord[k]-pte3,Oeil-pte3)>0:
	ll:=ll+1;
	Pays[k]=Coord[k];
      else:
	Pays[k]=2*Coord[k];
      fi;
    else:
      Pays[k]=2*Coord[k];
    fi;
  endfor;
  closefrom nomfichier;
  path pays;
  if ll>0:
    pays=Projgeo(Pays[1])
    for l=2 upto nbpts:
      --Projgeo(Pays[l])
    endfor;
    remplis pays--cycle withcolor fond;
    trace pays;
    clip currentpicture to cercles(pte3,pte1,pte3,pte1,pte4);
  fi;
enddef;

vardef lecturemanuel(expr nomfichier,fond)=
  color Coord[];
  numeric ll;
  ll:=0;
  color Pays[];
  nbpts:=scantokens readfrom nomfichier;
  for k=1 upto nbpts:
    pair latlon;
    latlon=scantokens readfrom nomfichier;
    Coord[k]=ray*(cosd(xpart(latlon/60))*cosd(ypart(latlon/60)),
	cosd(xpart(latlon/60))*sind(ypart(latlon/60)),sind(xpart(latlon/60)));
    if ((xpart(latlon/60)>phim) and (xpart(latlon/60)<phip)):
      if ProduitScalaire(Coord[k]-pte3,Oeil-pte3)>0:
	ll:=ll+1;
	Pays[k]=Coord[k];
      else:
	Pays[k]=2*Coord[k];
      fi;
    else:
      Pays[k]=2*Coord[k];
    fi;
  endfor;
  closefrom nomfichier;
  path pays;
  if ll>0:
    pays=Projgeo(Pays[1])
    for l=2 upto nbpts:
      --Projgeo(Pays[l])
    endfor;
    remplis pays--cycle withcolor fond;
    trace pays;
    %clip currentpicture to cercles(pte3,pte1,pte3,pte1,pte4);
  fi;
enddef;

vardef lectureexemple(expr nomfichier,fond)=
  color Coord[];
  numeric ll;
  ll:=0;
  color Pays[];
  nbpts:=scantokens readfrom nomfichier;
  for k=1 upto nbpts:
    pair latlon;
    latlon=scantokens readfrom nomfichier;
    Coord[k]=ray*(cosd(xpart(latlon/60))*cosd(ypart(latlon/60)),
	cosd(xpart(latlon/60))*sind(ypart(latlon/60)),sind(xpart(latlon/60)));
  endfor;
  closefrom nomfichier;
  path pays;
  pays=Projgeo(Coord[1])
    for l=2 upto nbpts:
      --Projgeo(Coord[l])
    endfor;
    trace pays dashed evenly;
    %clip currentpicture to cercles(pte3,pte1,pte3,pte1,pte4);
enddef;


ray:=2;

figureespace(-100u,-100u,100u,100u);
Initialisation(5,-115,-10,750);
numeric phim,phip,phii;%phi moins -- phi plus - phi intermédiaire
phim=Phi+arcsind(ray/Rho)-90;
phip=Phi+90-arcsind(ray/Rho);
color pte[];
pte1=ray*(cosd(phim)*cosd(Theta),cosd(phim)*sind(Theta),sind(phim));
pte2=ray*(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;
remplis cercles(pte3,pte1,pte3,pte1,pte4) withcolor ciel;
%%amerique du sud :) Manque l'arctique.
lecturemanuel("../DATA/bresil.dat",vert);
trace cercles(pte3,pte1,pte3,pte1,pte4);
finespace;

figureespace(-100u,-100u,100u,100u);
Initialisation(5,-115,-10,750);
numeric phim,phip,phii;%phi moins -- phi plus - phi intermédiaire
phim=Phi+arcsind(ray/Rho)-90;
phip=Phi+90-arcsind(ray/Rho);
color pte[];
pte1=ray*(cosd(phim)*cosd(Theta),cosd(phim)*sind(Theta),sind(phim));
pte2=ray*(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;
remplis cercles(pte3,pte1,pte3,pte1,pte4) withcolor ciel;
%%amerique du sud :) Manque l'arctique.
lecture("../DATA/bresil.dat",vert);
lectureexemple("../DATA/bresil.dat",noir);
trace cercles(pte3,pte1,pte3,pte1,pte4);
finespace;
end