diff --git a/Documentation_texfol/Graphiques/.gitignore b/Documentation_texfol/Graphiques/.gitignore index 12f9b6c9a28702b21c3aada099e1160924442634..5015fed72c6c632394ccbf0e8c4a6d36caf5ee2d 100644 --- a/Documentation_texfol/Graphiques/.gitignore +++ b/Documentation_texfol/Graphiques/.gitignore @@ -7,3 +7,4 @@ regions.pdf window.pdf periodicity.pdf copy.pdf +set_all_outerm.pdf diff --git a/Documentation_texfol/Graphiques/GNUmakefile b/Documentation_texfol/Graphiques/GNUmakefile index 4d57024ff306fc4c36d49c58fabe9c2f604a0fc9..fd3e40cccbd0c76de71b1586807e670c3c6e025d 100644 --- a/Documentation_texfol/Graphiques/GNUmakefile +++ b/Documentation_texfol/Graphiques/GNUmakefile @@ -1,4 +1,4 @@ -sources = window processes m3 call_graph input_output degeneracy periodicity copy +sources = window processes m3 call_graph input_output degeneracy periodicity copy set_all_outerm objects := $(addsuffix .pdf, ${sources}) diff --git a/Documentation_texfol/Graphiques/set_all_outerm.odg b/Documentation_texfol/Graphiques/set_all_outerm.odg new file mode 100644 index 0000000000000000000000000000000000000000..821cd1c1a7f49f3bc92e6c4109fd221f302bd866 Binary files /dev/null and b/Documentation_texfol/Graphiques/set_all_outerm.odg differ diff --git a/Documentation_texfol/documentation.tex b/Documentation_texfol/documentation.tex index 3fc1d30842a53fcbd5fb52fe119bfe9ac89216a6..421054a5cbe8f9e9f1997700a7b476c82cd41ee3 100644 --- a/Documentation_texfol/documentation.tex +++ b/Documentation_texfol/documentation.tex @@ -519,7 +519,7 @@ Les procédures géométriques utilisées : travaillent dans le plan et non sur la sphère. Ces procédures sont appelées par : \begin{itemize} -\item spherical\_polyline\_area +\item spher\_polyline\_area \item inside\_4 \item good\_contour \end{itemize} @@ -544,7 +544,7 @@ procédures qui appellent les procédures géométriques doivent se charger d'ajouter le multiple de $2 \pi$ pour ramener les longitudes dans un intervalle de longueur strictement inférieure à $\pi$. La contrainte de domaine de longitudes s'applique aux procédures -\verb+get_1_outerm+ et \verb+set_max_speed+. +\verb+get_1_outerm+ et \verb+set_max_speed+. Cf. § \ref{sec:set_all_outerm}. \subsection{Sous-algorithmes} @@ -656,6 +656,7 @@ Cas de périodicité : field(i + m - 2, :) = field(i, :). Il manque des tests pour lesquels direction vaut 1 et 4. \subsubsection{set\_all\_outerm} +\label{sec:set_all_outerm} Notons ici : \begin{align*} @@ -666,6 +667,15 @@ Notons ici : & w = \mathtt{corner\_window(1)} \\ & m = \mathtt{max\_radius(1)} \end{align*} +Cf. figure (\ref{fig:set_all_outerm}). +\begin{figure}[htbp] + \centering + \includegraphics{set_all_outerm} + \caption{Contrainte sur les longitudes dans set\_all\_outerm. Les + valeurs en vert sont dans l'espace des indices, les valeurs en + noir dans l'espace des longitudes.} + \label{fig:set_all_outerm} +\end{figure} Pour que la contrainte de domaine de longitudes dans \verb+get_1_outerm+ soit satisfaite, il faut que, dans set\_all\_outerm, au moment de l'appel à \verb+get_1_outerm+ : @@ -734,7 +744,7 @@ weight à u et v, c'est-à -dire qu'il faudrait conserver u, v à plusieurs dates en dehors de \verb+get_snpashot+. Cela remettrait en question tout l'algorithme. -\subsubsection{spherical\_polyline\_area} +\subsubsection{spher\_polyline\_area} La procédure de calcul des contours donne une suite de sommets. Les courbes reliant ces sommets ne sont en fait pas imposés @@ -971,7 +981,7 @@ Les procédures géométriques utilisées : \end{itemize} Ces procédures sont appelées par : \begin{itemize} -\item spherical\_polyline\_area +\item spher\_polyline\_area \item \verb+successive_overlap+ \end{itemize} Cf. discussion autour de l'équation~(\ref{eq:length_pi}). @@ -1512,6 +1522,25 @@ dans polygon\_1 et polygon\_2, il faudrait faire des copies de gagnerait rien en quantité de copie, et on perdrait en mémoire vive occupée. +Contrainte géométrique : les extremums i1 et i2 doivent être ramenés, +via une modulation de $2 \pi$, à une distance inférieure à $pi$ en +longitude. Si $\lambda_1$ désigne la longitude de l'extremum 1, et +$\lambda$ une longitude quelconque, on cherche $m \in \mathbb{Z}$ tel +que : +\begin{equation*} + |\lambda + 2 m \pi - \lambda_1| < \pi +\end{equation*} +Si $m$ existe (ce qui est le cas pour $\lambda$ correspondant à i2) +alors : +\begin{equation*} + m = E\left(\frac{\lambda_1 - \lambda}{2 \pi} + \frac{1}{2} \right) +\end{equation*} +Dans \verb+set_all_outerm+, le multiple de $2 \pi$ à ajouter aux +longitudes de \verb+outside_points+ peut varier selon le point. Ici, +le multiple de $2 \pi$ à ajouter aux longitudes doit être le même pour +tous les points du contour, et égal au multiple à ajouter à la +longitude de l'extremum i2. + \subsubsection{weight} D'autant plus proche de 0 que les tourbillons sont diff --git a/extraction_eddies.f90 b/extraction_eddies.f90 index 53d7c896a6424143d591488734a4d6bce8b23d52..bcf2744974d45c6285a2ae86858853e548773892 100644 --- a/extraction_eddies.f90 +++ b/extraction_eddies.f90 @@ -112,8 +112,9 @@ program extraction_eddies ! "periodic" must be consistent with the values of step(1) and nlon: periodic = nint(twopi / step(1)) == nlon print *, "periodic = ", periodic - if (periodic) call assert(2 * max_radius(1) * step(1) < pi, & - "extraction_eddies max_radius") + + call assert(2 * max_radius(1) * step(1) < pi, "extraction_eddies max_radius") + ! (Let us require this even if not periodic. This is clearer.) copy = merge(max_radius(1), 0, periodic) allocate(ssh(1 - copy:nlon + copy, nlat), u(1 - copy:nlon + copy, nlat), & diff --git a/successive_overlap.f90 b/successive_overlap.f90 index 6a9c64bd5fcb2dbc5b60c605acae9cd0a6816794..e6088a9992258391f90d30f24ff5fcb2cd41a6a9 100644 --- a/successive_overlap.f90 +++ b/successive_overlap.f90 @@ -17,6 +17,7 @@ contains ! Libraries: use contour_531, only: polyline use gpc_f, only: gpc_polygon_clip_f, GPC_INT, polygon + use nr_util, only: twopi use candidate_overlap_m, only: candidate_overlap use derived_types, only: snapshot @@ -32,7 +33,7 @@ contains ! Local: - integer i1, i2, l, n_select + integer i1, i2, l, n_select, m type(polyline) polyline_1, polyline_2 type(polygon) res_pol @@ -74,6 +75,13 @@ contains polyline_2 = merge(flow(j)%list_vis(i2)%speed_cont%polyline, & flow(j)%list_vis(i2)%out_cont%polyline, & flow(j)%list_vis(i2)%speed_cont%n_points /= 0) + + ! Shift the longitudes of polyline_2 to values close to the + ! longitude of extremum i1: + m = floor((flow(j - 1)%list_vis(i1)%coord_extr(1) & + - flow(j)%list_vis(i2)%coord_extr(1)) / twopi + 0.5) + polyline_2%points(1, :) = polyline_2%points(1, :) + m * twopi + call gpc_polygon_clip_f(GPC_INT, polygon(nparts = 1, & part = [polyline_1], hole = [.false.]), & polygon(nparts = 1, part = [polyline_2], hole = [.false.]), &