From 5ad855f5358bec2bda470d39e0fc23b3582da601 Mon Sep 17 00:00:00 2001 From: Jean-Denis Vauguet Date: Wed, 27 May 2009 16:08:08 +0200 Subject: [PATCH] plop --- .../IPGP Saint Maur/MATLAB/vectorial/modesT_vect.m | 13 +- .../MATLAB/vectorial/modesT_vect_wind_M.m | 537 --------------------- 2 files changed, 5 insertions(+), 545 deletions(-) delete mode 100644 docs/fac/stage/IPGP Saint Maur/MATLAB/vectorial/modesT_vect_wind_M.m diff --git a/docs/fac/stage/IPGP Saint Maur/MATLAB/vectorial/modesT_vect.m b/docs/fac/stage/IPGP Saint Maur/MATLAB/vectorial/modesT_vect.m index f6e7575..8d1d33f 100644 --- a/docs/fac/stage/IPGP Saint Maur/MATLAB/vectorial/modesT_vect.m +++ b/docs/fac/stage/IPGP Saint Maur/MATLAB/vectorial/modesT_vect.m @@ -47,12 +47,11 @@ dz(n_mod) = dz(n_mod-1); % FIXME ajouté pour la version vectorisée, b dz = dz'; % conditions limites sur la vitesse horizontale : pour le calcul du gradient de vent -u0b = -50; % vitesse en bas (bottom) -u0t = -52; % vitesse en haut (top) +u0b = 0; % vitesse en bas (bottom) +u0t = 0; % vitesse en haut (top) % coefficient d'un gradient linéaire des vents égal à d(u0)/dz alpha = (u0t-u0b)/((Z(n_mod)-Z(1))*1.e3); alpha -pause % profil des vitesses moyennes u0 = Z.*(1.e3*alpha) + u0b; @@ -188,8 +187,6 @@ for ik = 1:nlon % (nlon-1)/2 + 2:nlon p1 = i/(k*k)*(alpha*k - sqrt(kz2)*(omega - k*u0(1)) - 0.5*(omega - k*u0(1))*dlgrho(1))*w1 end - pause - % calcul des coefficients pour les conditions initiales au step 1 %b1 = (i.*dz(1).*k.*alpha)./BigOmega(1) - 0.5.*dz(1).*dlgrho(1) - 1 %b2 = (i*k*k)./BigOmega(1) @@ -352,11 +349,11 @@ for ik = 1:nlon % (nlon-1)/2 + 2:nlon figure(100); F = getframe(gcf); [X, Map] = frame2im(F); - file = ['lineaire_L' num2str(round(1/k*1.e-3)) '_T_' num2str(2*pi/omega/60) '.tif']; + file = ['vent0_L' num2str(round(1/k*1.e-3)) '_T_' num2str(2*pi/omega/60) '.tif']; imwrite(X, file, 'tif', 'Compression', 'none'); disp('PRESS ENTER TO CONTINUE...'); - pause + %pause close figure(300); @@ -381,7 +378,7 @@ for ik = 1:nlon % (nlon-1)/2 + 2:nlon hold off; xlabel('lambda Z (km)', 'fontsize', 12); - pause + %pause close figure(2*ik); diff --git a/docs/fac/stage/IPGP Saint Maur/MATLAB/vectorial/modesT_vect_wind_M.m b/docs/fac/stage/IPGP Saint Maur/MATLAB/vectorial/modesT_vect_wind_M.m deleted file mode 100644 index 628af69..0000000 --- a/docs/fac/stage/IPGP Saint Maur/MATLAB/vectorial/modesT_vect_wind_M.m +++ /dev/null @@ -1,537 +0,0 @@ -% TEST GRAVITY WAVES % - -% FIXME -% linéarisation et taille des vecteurs de données (cf. plus bas) -% ensuite, changer le profil des vents, et regarder comment est modifiée l'onde selon la direction, pour une longueur d'onde donnée et un mode d'excitation donné… - -% version de jd, vectorisée - -% pour le moment, il n'y a pas distinction entre les deux axes horizontaux, -% c'est du 3D planaire disons - -clear all; -close all; - -% modele1: -% altitude (0 - 500km) -% densité -% dérivée de la densité -% dérivée du logarithme de la densité -load modele1 -load wind_M - - -%Z = RHOout(:,1); -%Z = modele1(:,1); -%nnn = max(size(Z)); % 9999 (1D) -%rho = modele1(:,2); -%drho = modele1(:,3); -%dlgrho = modele1(:,4); - -% tests à résolution variable -res = 1 % 10, 100, etc. -binf = 1 -bsup = 9981 % jusqu'à 499.00 km -Z = modele1(binf:res:bsup,1); -rho = modele1(binf:res:bsup,2); -drho = modele1(binf:res:bsup,3); -dlgrho = modele1(binf:res:bsup,4); - -xx = size(Z); -n_mod = xx(1) - -I(1:n_mod) = 1; - -% calcul de la dérivée spatiale verticale -% TODO à vectoriser ^^ -for j = 1:n_mod-1 - dz(j) = (Z(j+1)-Z(j))*1.e3; % pas spatial en mètres, utilisé dans les -end % dérivées en dz j'imagine -dz(n_mod) = dz(n_mod-1); % FIXME ajouté pour la version vectorisée, bof -dz = dz'; - -%% conditions limites sur la vitesse horizontale : pour le calcul du gradient de vent -%u0b = -50; % vitesse en bas (bottom) -%u0t = -52; % vitesse en haut (top) -%% coefficient d'un gradient linéaire des vents égal à d(u0)/dz -%alpha = (u0t-u0b)/((Z(n_mod)-Z(1))*1.e3); -%alpha = alpha.*ones(n_mod,1); -%% profil des vitesses moyennes -%u0 = Z.*(1.e3.*alpha) + u0b; - -u0 = wind_M(:,2); - -% interpolation linéaire au pas de modele1 -%size(u0) - -res = []; -for j = 2:size(u0) - tmp = linspace(u0(j-1), u0(j), 21); - res = [res tmp(1:end-1)]; % end-1 pour ne pas avoir de doublons sur la borne sup -end -res = [res tmp(end)]; -u0 = res'; -%u0(1:3) -%u0(end-3:end) - -%size(u0) -%size(dz) -%u0diff = diff(u0); -%size(u0diff) -%alpha = u0diff./dz; -%alpha = gradient(u0, dz); -%size(alpha) -alpha = []; -for j = 2:n_mod-1 - alpha = [alpha (u0(j+1) - u0(j-1))/dz(j)]; -end -alpha = [(u0(2)-u0(1))/dz(2) alpha (u0(end)-u0(end-1))/dz(end)]'; -alpha(1:5) -alpha(end-5:end) -size(alpha) - -pause - -M = 5.9768e+24; % masse de la Terre -G = 6.67e-11; % constante gravitationnelle -R = 6378.e3; % rayon de la Terre - -%---------------------------------- -% omega, kx & ky é la main -%---------------------------------- - - % si j'ai bien compris, - % on étudie la réponse à une altitude donnée pour -T = [50 100 150].*60 % trois périodes de temps : 3000, 6000, 9000 secondes -wn = 2*pi./T; % les trois fréquences associées... - -h = 2500; % hauteur de ? -g = 9.8; -c = sqrt(g*h); % vitesse de phase des ondes acoustiques simples - % elle est imposée par continuité à l'interface, c'est - % une condition initiale/limite -kx = wn./c; % nombres d'ondes horizontaux -ky = kx; % . - -%---------------------------------- -% itérateurs spatio-temporels - -nlat = max(size(ky)); % en fait, on a fait en sorte que nlon = nlat... -nlon = max(size(kx)); -nt = max(size(wn)); - -%---------------------------------- -% Calculs principaux -% L'objectif est de calculer les modes pour les IGW, ie. la valeur des paramètres -% des équations du système à chaque altitude (« mode ») pour laquelle est fixé -% un kz approprié, sous la contrainte des conditions limites. La condition limite -% bottom fait office de condition initiale, couplée à une condition sur w, une -% condition sur P (continuité à l'interface). -% -% On considère trois périodes de temps, donc trois fréquences wn, donc trois -% nombre d'ondes k à travers la vitesse de phase imposée comme condition au limite -% implicite par continuité à l'interface. Il y a donc trois cas étudiés (ik). -% Pour chaque cas, on analyse les trois fréquences (iw). La boucle sur les iw -% résout le propagateur, sous les conditions initiales/limites. -%---------------------------------- - -load BestView2 - -%---------------------------------- - -% boucle sur les nombres d'onde (sur les longeurs d'onde), ici trois. -% une figure par mode ik -for ik = 1:nlon % (nlon-1)/2 + 2:nlon - ik - drawnow; % DRAWNOW causes figure windows and their children to update and - % flushes the system event queue. Any callbacks generated by incoming - % events - e.g. mouse or key events - will be dispatched before - % DRAWNOW returns. - - k = kx(ik) % nombre d'onde horizontal courant - - figure(2*ik); % FIGURE(H) makes H the current figure, forces it to become visible, - % and raises it above all other figures on the screen. If Figure H - % does not exist, and H is an integer, a new figure is created with - % handle H. - set(2*ik, 'position', [44 250 500 750]); - set(gcf, 'Color', 'w'); - set(gca, 'LineWidth', [2]) - set(gca, 'fontsize', 18); - view(VV); - title(['Hwater = ', num2str(h), 'm and kx = ', num2str(k), ' rad/m'], 'fontsize', 18); - % xlabel('omega w (rad/s)', 'fontsize', 18); - xlabel('period T (min)', 'fontsize', 18); - ylabel('Vertical V_r_e_a_l (m/s)', 'fontsize', 18); - zlabel('altitude (km)', 'fontsize', 18); - grid on; - hold on; - - figure(2*ik+1); - set(2*ik+1, 'position', [44 250 500 750]); - set(gcf, 'Color', 'w'); - set(gca, 'LineWidth', [2]) - set(gca, 'fontsize', 18); - view(VV); - title(['Hwater = ', num2str(h), 'm and kx = ', num2str(k), ' rad/m'], 'fontsize', 18); - % xlabel('omega (rad/s)', 'fontsize', 18); - xlabel('period T (min)', 'fontsize', 18); - ylabel('Vertical V_i_m_a_g (m/s)', 'fontsize', 18); - zlabel('altitude (km)', 'fontsize', 18); - grid on; - hold on; - - % FIXME relou ces figures :) - close(2*ik); - close(2*ik+1); - - % boucle sur les trois fréquences par mode, sachant qu'un mode est une altitude ici - % avec nt = max(size(wn)) == 3 ici - for iw = 1:nt %(nt-1)/2 + 2:nt - - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - % vectorisation en K*U = 0 - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - omega = wn(iw) - BigOmega = omega - k.*u0; - - % avec le propagateur en différences finies - - % calcul des coefficients pour les blocs diagonaux de la matrice creuse - %c1 = (i.*2.*dz.*k.*alpha)./BigOmega - dz.*dlgrho; - %c2 = (i*k*k)./BigOmega; - %c3 = dz.*dlgrho; - %c4 = -(2*i).*dz.*(BigOmega + (g./BigOmega).*dlgrho); - - % tiré de la version de ninto, au signe près - c1 = -(-k.*alpha.*2.*dz./BigOmega + dlgrho.*dz); - c2 = -i.*k.*k.*2.*dz./BigOmega; - c3 = -(-dlgrho.*dz); - c4 = -(i.*2.*dz.*(BigOmega + dlgrho.*g./BigOmega)); - - % condition initiale/limite sur P et w - % calcul de la pression initiale sous contrainte de continuité - w1 = complex(1,0); % FIXME justification ? - kz2 = k*k*(-g/omega/omega*dlgrho(1) - 1) - 0.25/6.4e7; - N2 = -g*dlgrho(1); - if N2 >= omega*omega - % cas d'une onde propagative - disp 'onde propagative vent profil' - p1 = i/(k*k)*(alpha(1)*k - i*sqrt(kz2)*(omega - k*u0(1)) - 0.5*(omega - k*u0(1))*dlgrho(1))*w1 - else - % cas d'une onde non propagative - p1 = i/(k*k)*(alpha(1)*k - sqrt(kz2)*(omega - k*u0(1)) - 0.5*(omega - k*u0(1))*dlgrho(1))*w1 - end - - pause - - % calcul des coefficients pour les conditions initiales au step 1 - %b1 = (i.*dz(1).*k.*alpha)./BigOmega(1) - 0.5.*dz(1).*dlgrho(1) - 1 - %b2 = (i*k*k)./BigOmega(1) - %b3 = 0.5.*dz(1).*dlgrho(1) - 1 - %b4 = -(2*i).*dz(1).*(BigOmega(1) + (g./BigOmega(1)).*dlgrho(1)) - - % tiré de la version de ninto, au signe près - b1 = -(-i*k*dz(1) * (1 / BigOmega(1) * (-i.*alpha(1))) + 0.5*dlgrho(1)*dz(1) + 1) - b2 = -(-i*k*dz(1) * (1 / BigOmega(1) * (k))) - b3 = -(-0.5*dlgrho(1)*dz(1) + 1) - b4 = -(i*BigOmega(1)*dz(1) + i/BigOmega(1)*dlgrho(1)*dz(1)*g) - - % vecteurs colonnes à 0 alternés pour spdiags() - % attention, permutation circulaire pour les diagonales d'indices impairs* - C13 = [c1'; c3']; - C2 = [zeros(1, size(c2,1)); c2']; % * - C4 = [c4'; zeros(1, size(c4,1))]; - - diagC13 = C13(:); - %size(diagC13) - diagC2 = C2(:); - %size(diagC2) - % FIXME pb de signe ici ! en attendant, je mets un moins pour retrouver - % des parties imaginaires positives... - diagC4 = -C4(:); - %size(diagC4) - - % matrice des conditions limites - cl0 = [1 0; 0 1]; - cl1 = [b1 b2 1 0; b4 b3 0 1]; - %cl = [cl0 zeros(size(cl0,1), n_mod-size(cl0,1)) ; cl1 zeros(size(cl1,1), n_mod-size(cl1,1))] - n_mod2 = 2*n_mod; - cl = [cl0 zeros(2, n_mod2-2) ; cl1 zeros(2, n_mod2-4)]; - - % construction de la matrice creuse des coefficients du système linéaire pour - % l'itération en différences finies - % partie utile de la matrice diagonale à laquelle sont ajoutées les 2*2 - % conditions limites/initiales - - plop = ones(n_mod2, 1); - preK = spdiags([-plop diagC4 diagC13 diagC2 plop], 0:4, n_mod2, n_mod2); - preK = preK(1:n_mod2 - 4, :); - %size(cl) - %size(preK) - K = [cl; preK]; - - % calcul du condition number, histoire de… - %kn = svd(K)/svd(1/K); - %disp('kn = ') - %kn - - %Kf = full(K); - %Kf(1:8,1:8) - %diagC2(1) - %diagC2(2) - %diagC2(3) - %diagC2(4) - %Kf(n_mod2-7:end,n_mod2-7:end) - %size(Kf) - - % résolution du système linéaire K*U = 0 ------------------------------- - % attention : faux pour la condition initiale qui n'est pas exprimée en - % différences finies ! - % on construit le vecteur reste R nul partout sauf pour la CL - R = [w1; p1; zeros(n_mod2-2, 1)]; - U = K\R; - size(U) - U(1:10) - U(n_mod2-9:end) - p1 - - % récupération des résultats en w et p - w = U(1:2:end); - %size(w) - p = U(2:2:end); - %size(p) - - % calcul du profil des vitesses intégrées explicites - u = (-i.*alpha.*w + k.*p)./BigOmega; - - % calcul des kz explicites - % FIXME coefficient à remplacer par la formule analytique, non ? - kz2 = (k*k).*((-g/omega/omega).*dlgrho-1) - 0.25/6.4e7; - - % calcul des rho1 explicites - rho1 = (-i.*dlgrho.*w)./BigOmega; - - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - %------------------------------------------------------------ - % da usare con kx e omega é la main - % - wmax(ik, iw) = max(abs(real(w))); - umax(ik, iw) = max(abs(real(u))); - pmax(ik, iw) = max(abs(real(p))); - rhomax(ik, iw) = max(abs(real(rho1))); - KZ(ik, iw, :) = sqrt(kz2); - - if max(abs(real(w))) >= 2. - ctrw(ik, iw) = NaN; - else - ctrw(ik, iw) = 1.; - end - - if sqrt(kz2(j)) >= 2*pi/50.e3 - A(ik, iw) = 1; - else - A(ik, iw) = NaN; - end - - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - figure(100); - set(100, 'position', [10 1 750 600]) - set(gcf, 'Color', 'w'); - fnt = 20; - fnt2 = 18; - Zmax = max(Z); - Zmin = min(Z); - subplot(1, 4, 1); - plot(real(w), Z, 'r', 'linewidth', 2); - hold on; - plot(imag(w), Z, 'g', 'linewidth', 2); - hold off; - %axis([-2 2 0 max(Z)]) - v = axis; - axis([v(1) v(2) Zmin Zmax]) - xlabel('Vertical V (m/s)', 'fontsize', 12); - ylabel('Altitude (km)', 'fontsize', 16); - title(['Period = ', num2str(2*pi/omega/60), ' mn' ' (omega = ', num2str(omega), ' rad/s)'], 'fontsize', 12); - subplot(1, 4, 2); - plot(real(u), Z, 'r', 'linewidth', 2); - hold on; - plot(imag(u), Z, 'g', 'linewidth', 2); - hold off; - %axis([-20 20 0 max(Z)]) - v = axis; - axis([v(1) v(2) Zmin Zmax]); - xlabel('Horizontal V (m/s)', 'fontsize', 12); - subplot(1, 4, 3); - plot(real(p), Z, 'r', 'linewidth', 2); - hold on; - plot(imag(p), Z, 'g', 'linewidth', 2); - hold off; - %axis([-3000 3000 0 max(Z)]) - v = axis; - axis([v(1) v(2) Zmin Zmax]); - xlabel('Pressure (Pa)', 'fontsize', 12); - subplot(1, 4, 4); - plot(real(rho1), Z, 'r', 'linewidth', 2); - hold on; - plot(imag(rho1), Z, 'g', 'linewidth', 2); - hold off; - %axis([-0.1 0.1 0 max(Z)]) - v = axis; - axis([v(1) v(2) Zmin Zmax]); - xlabel('Density (kg/m3)', 'fontsize', 12); - title(['Lambda = ', num2str(round(1/k*1e-3)), ' km' ' (kx = ', num2str(k), ' rad/m)'], 'fontsize', 12); - - figure(100); - F = getframe(gcf); - [X, Map] = frame2im(F); - file = ['Earth_Mode_L' num2str(round(1/k*1.e-3)) '_T_' num2str(2*pi/omega/60) '.tif']; - file = ['profil_L' num2str(round(1/k*1.e-3)) '_T' num2str(2*pi/omega/60) '.tif']; - imwrite(X, file, 'tif', 'Compression', 'none'); - - disp('PRESS ENTER TO CONTINUE...'); - pause - close - - figure(300); - subplot(1, 3, 1); - plot(real(sqrt(kz2)), Z, 'r'); - hold on; - plot(imag(sqrt(kz2)), Z, 'g'); - xlabel('kz (rad/m)', 'fontsize', 12); - ylabel('Altitude (km)', 'fontsize', 16); - hold off; - subplot(1, 3, 2); - plot(sqrt(N2), Z, 'r'); - hold on; - v = axis; - plot([omega omega], [v(3) v(4)], 'g'); - hold off; - xlabel('N and omega (rad/s)', 'fontsize', 12); - subplot(1, 3, 3); - plot(2*pi./real(sqrt(kz2))*1.e-3, Z, 'r'); - hold on; - plot(2*pi./imag(sqrt(kz2))*1.e-3, Z, 'g'); - hold off; - xlabel('lambda Z (km)', 'fontsize', 12); - - pause - close - - figure(2*ik); - plot3(I*2*pi/omega/60, real(w), Z, 'k'); - plot3(I*2*pi/omega/60, imag(w), Z, 'r'); - - figure(2*ik+1); - plot3(I*2*pi/omega/60, imag(w), Z, 'k'); - - clear w u p rho1 kz2 - end -end - -Nsol = sqrt(N2(1)); -Nmax = max(sqrt(N2)); -Nmin = min(sqrt(N2)); - -njump = 0; -njump2 = size(wmax); - -n1 = 1; -n2 = 1; - -h = figure; -set(h, 'position', [2 700 800 600]); -set(gcf, 'Color', 'w'); -font_size = 22; -font_sizeN = 18; - -flag_plot = input('-Plot omega-k or T-lambda ? (0 or 1) '); - -subplot(2, 2, 1); -if flag_plot == 0 - pcolor(wn(n1:nt), kx(n2:nlon), log10(wmax(:,1+njump:njump2(2)))); -else - pcolor(2*pi./wn(n1:nt)/60, 2*pi./kx(n2:nlon)*1.e-3, log10(wmax(:,1+njump:njump2(2)))); -end -shading flat; -v = axis; -hold on; -plot([Nmax Nmax], [v(3) v(4)], 'k'); -plot([Nsol Nsol], [v(3) v(4)], 'k'); -plot([Nmin Nmin], [v(3) v(4)], 'k'); -colorbar; -title('W_r_e_a_l', 'fontsize', font_size); -if flag_plot == 0 - ylabel('kx (rad/m)', 'fontsize', font_size); -else - ylabel('lambda (km)', 'fontsize', font_size); -end -set(gca, 'LineWidth', [2]) -set(gca, 'fontsize', font_sizeN); - -subplot(2, 2, 2); -if flag_plot == 0 - pcolor(wn(n1:nt), kx(n2:nlon), log10(umax(:,1+njump:njump2(2)))); -else - pcolor(2*pi./wn(n1:nt)/60, 2*pi./kx(n2:nlon)*1.e-3, log10(umax(:,1+njump:njump2(2)))); -end -shading flat; -v = axis; -hold on; -plot([Nmax Nmax], [v(3) v(4)], 'k'); -plot([Nsol Nsol], [v(3) v(4)], 'k'); -plot([Nmin Nmin], [v(3) v(4)], 'k'); -colorbar; -title('U_r_e_a_l', 'fontsize', font_size); -set(gca, 'LineWidth', [2]) -set(gca, 'fontsize', font_sizeN); - -subplot(2, 2, 3); -if flag_plot == 0 - pcolor(wn(n1:nt), kx(n2:nlon), log10(pmax(:,1+njump:njump2(2)))); -else - pcolor(2*pi./wn(n1:nt)/60, 2*pi./kx(n2:nlon)*1.e-3, log10(pmax(:,1+njump:njump2(2)))); -end -shading flat; -v = axis; -hold on; -plot([Nmax Nmax], [v(3) v(4)], 'k'); -plot([Nsol Nsol], [v(3) v(4)], 'k'); -plot([Nmin Nmin], [v(3) v(4)], 'k'); -colorbar; -if flag_plot == 0 - ylabel('kx (rad/m)', 'fontsize', font_size); - xlabel('omega (rad/s)', 'fontsize', font_size); -else - ylabel('lambda (km)', 'fontsize', font_size); - xlabel('period T (min)', 'fontsize', font_size); -end -title('P_r_e_a_l', 'fontsize', font_size); -set(gca, 'LineWidth', [2]) -set(gca, 'fontsize', font_sizeN); - -subplot(2, 2, 4); -if flag_plot == 0 - pcolor(wn(n1:nt), kx(n2:nlon), log10(rhomax(:,1+njump:njump2(2)))); -else - pcolor(2*pi./wn(n1:nt)/60, 2*pi./kx(n2:nlon)*1.e-3, log10(rhomax(:,1+njump:njump2(2)))); -end -shading flat; -v = axis; -hold on; -plot([Nmax Nmax], [v(3) v(4)], 'k'); -plot([Nsol Nsol], [v(3) v(4)], 'k'); -plot([Nmin Nmin], [v(3) v(4)], 'k'); -colorbar; -if flag_plot == 0 - xlabel('omega (rad/s)', 'fontsize', font_size); -else - xlabel('period T (min)', 'fontsize', font_size); -end -title('Rho_r_e_a_l', 'fontsize', font_size); -set(gca, 'LineWidth', [2]) -set(gca, 'fontsize', font_sizeN); - -- 2.11.4.GIT