le programme snell et la creation prealable d'une bathymetrie sont neccessaire au lancement de refract.m

% programme de réfraction de la houle
global omega profondeur

nb_point_vague=10;
nb_crete=15;
omega=0.5;
n=5;
nb_iteration=nb_crete*n

%creation de A
A(1,1:10)=[0,0,0,0,0,0,0,0,0,0];
%A(2,1:10)=[500,700,900,1100,1300,1500,1700,1900,2100,2300];
A(2,1:10)=[175:100:1075];
%depart en grands fonds
k=omega^2/9.81;
A(3,1:10)=[k,k,k,k,k,k,k,k,k,k];
A(4,1:10)=[0,0,0,0,0,0,0,0,0,0];
A(5,1:10)=[0,0,0,0,0,0,0,0,0,0];

A;

%boucle sur les points de la vague
for i=1:nb_point_vague
   i

%itération sur un point de la vague
for l=1:nb_crete

%itération sur un point de la vague
for j=1:n

   %deplacement x et y
   A(1,i)=A(1,i)+2*pi/(A(3,i)*n)*cos(A(4,i));
   A(2,i)=A(2,i)+2*pi/(A(3,i)*n)*sin(A(4,i));
   %localisation de X dans matrice B
   inf_x=fix(A(1,i)/pas_bath+1);
   %localisation de Y dans matrice B
   inf_y=fix(A(2,i)/pas_bath+1);
   %coeff du barycentre
   coeff_x=((A(1,i)/pas_bath+1)-inf_x);
   coeff_y=((A(2,i)/pas_bath+1)-inf_y);

   %calcul du barycentre
   r1=sqrt(coeff_x^2+coeff_y^2);
   r2=sqrt(coeff_x^2+(1-coeff_y)^2);
   r3=sqrt((1-coeff_x)^2+coeff_y^2);
   r4=sqrt((1-coeff_x)^2+(1-coeff_y)^2);

   somme=1/(r1^2)+1/(r2^2)+1/(r3^2)+1/(r4^2);
   profondeur=1/r1^2*B(inf_y,inf_x);
   profondeur=profondeur+1/r2^2*B(inf_y+1,inf_x);
   profondeur=profondeur+1/r3^2*B(inf_y,inf_x+1);
   profondeur=profondeur+1/r4^2*B(inf_y+1,inf_x+1);
   profondeur=profondeur/somme;

   %resolution du nouveau k
   k_old=A(3,i);
   kd=k_old*profondeur;

   if (kd)>3
      A(3,i)=omega^2/9.81;
      %'grands fonds'
   elseif (kd)<0.15
      A(3,i)=omega/sqrt(9.81*profondeur);
      %'faibles fonds'
   else
      fzero('snell',0);
      A(3,i)=abs(ans);
      %'fonds moyens'
   end

   %resolution du nouveau theta
   d_dy=(1-coeff_x)*(B(inf_y+1,inf_x)-B(inf_y,inf_x))/pas_bath;
   d_dy=d_dy+coeff_x*(B(inf_y+1,inf_x+1)-B(inf_y,inf_x+1))/pas_bath;

   d_dx=coeff_y*(B(inf_y+1,inf_x+1)-B(inf_y+1,inf_x))/pas_bath;
   d_dx=d_dx+(1-coeff_y)*(B(inf_y,inf_x+1)-B(inf_y,inf_x))/pas_bath;

   beta_old=A(5,i);

   if abs(d_dx)<= 0.001
      A(5,i)=pi/2;
      'problème'
   else
    A(5,i)=atan(d_dy/d_dx);
   end
   beta_moyen=(beta_old+A(5,i))/2;
   theta_si=k_old/A(3,i)*sin(A(4,i)-beta_moyen);

   A(4,i)=asin(theta_si)+beta_moyen;
   graph_x((l-1)*n+j,i)=A(1,i);
   graph_y((l-1)*n+j,i)=A(2,i);

end

%stockage des coordonnées des points pour le graphe

end
end

%affichage de la matrice
A

%partie graphique

for i=1:nb_point_vague
   plot(graph_x(1:n*nb_crete,i),graph_y(1:n*nb_crete,i))
   hold on
end
for j=1:n*nb_crete
   if mod(j,n)==0
   plot(graph_x(j,1:nb_point_vague),graph_y(j,1:nb_point_vague))
   hold on
   end
end

contour(By,Bx,B,50)
 
 
 
 

% programme de bathymetrie
clear all;
pas_bath=2.5;

%creation de la bathymetrie
for j=1:1400
   for i=1:500
      B(i,j)=-300*exp(-0.00016*sqrt(20*(i-250)^2+20*(j-1300)^2))+301;
   end
end

B;

for i=1:500
   Bx(i)=(i-1)*pas_bath;
end
for j=1:1400
   By(j)=(j-1)*pas_bath;
end

function y=snell(x)
global omega profondeur
y=sqrt(9.81*x*tanh(x*profondeur))-omega;