Traduction code MATLAB en Scilab

Bonjour, j'ai fais plusieurs codes sur les "CARACTÉRISATION DE SURFACE PAR LA MÉTHODE DU PONT CAPILLAIRE" en MATLAB,

j'aimerais savoir si quelqu'un pourrait traduire cela en SCILAB, svp !

La partie la plus dure est celle de la résolution d'équation différentielle de Laplace (du 2nd ordre). Savez vous résoudre celle ci sur SICLAB, sans passer par des solvers svp ?

Pour plus d'information vous pouvez m'envoyer un mail : janot92@hotmail.fr

Voici ci dessous les codes :

Fonction Dérivé de Y par rapport à z :

Code:

1
2
3
4
5
function [ v ] = DYdz( z,Y ) %v : paramètre de sortie ; z et y : paramètres d'entrée 
global lc; %lc = kappa 
v=[Y(2) 
  (1+Y(2)^2)/Y(1)+(1+Y(2)^2)^(3/2)*z/(lc^2)]; 
end

Fonction Point capillaire (en 2D) :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
function [z,Y] = ptcap2D (h,theta_deg) 
global lc; 
warning off; 
lc=2; 
R=100; 
z_lentille=linspace(0,10,101); 
r_lentille=sqrt(2*R.*(z_lentille-h)-((z_lentille-h).^2)); 
zpas=1; 
zini=h+0.01; 
for iii=1:1:4 
   zpas=zpas/(10); 
   for z0=zini:zpas:R/10 
       r0=sqrt(2*R.*(z0-h)-((z0-h)^2)); 
       alpha=asin(r0/R); 
       theta=theta_deg/180*pi; 
       rp0=cot(alpha+theta); 
       %solveur 
       [z,Y]=ode45('DYdz',[z0 0],[r0 rp0]); 
       taille=size(z); 
       zfinal=z(taille(1)); 
       taille2=size(Y); 
       y2final=Y(2,taille2(2)); 
       if (zfinal~=0 || abs(y2final)>100) 
           zini=z0-zpas; 
           [z,Y]=ode45('DYdz',[zini 0],[r0 rp0]); 
           break 
       end 
   end 
end 
%tracés (Attention : il faut mettre ces tracés en commentaire pour obtenir les autres courbes) 
hold on; 
 
plot(r_lentille,z_lentille); 
plot(-r_lentille,z_lentille); 
plot(Y(:,1),z); 
plot(-Y(:,1),z); 
axis equal; 
axis([-60 60 0 10]) 
end 
Fonction Pont capillaire (en 3D) : 
function [X,Y,Zlentille] = ptcap3D ( ) 
R=100; 
h=0; 
for zi=0:1:100 
      r_lentille=sqrt(2.*R.*(zi/10-h)-(zi/10-h).^2); 
   for the=0:100 
      X(zi+1,the+1)=r_lentille.*cos(the*3.6*pi/180); 
      Y(zi+1,the+1)=r_lentille.*sin(the*3.6*pi/180); 
   end 
end 
for xi=1:101 
   for yi=1:101 
      Zlentille(xi,yi)=((2*R+2*h)... 
         -sqrt(((-2*R-2*h)^2)-4.*(2*R*h+h^2+X(xi).^2+Y(yi).^2)))./2; 
   end 
end 
[Zpt,Ytmp]=ptcap2D(h,0);            % on se sert de ptcap2D 
tailleY=size(Ytmp); 
for zi=1:tailleY(1,1) 
   r_pt=Ytmp(zi,1); 
   for the=1:tailleY(1,1) 
      Xpt(zi,the)=r_pt.*cos((the-1)/tailleY(1,1)*360*pi/180); 
      Ypt(zi,the)=r_pt.*sin((the-1)/tailleY(1,1)*360*pi/180); 
   end 
end 
for xi=1:tailleY(1,1) 
   for yi=1:tailleY(1,1) 
      Zpt_tab(xi,yi)=Zpt(xi); 
   end 
end 
%tracés du graphe 3D 
figure( ) 
hold on 
mesh(X,Y,Zlent); % la function mess équivaut à un plot, mais en 3D 
mesh(Xpt,Ypt,Zpt_tab); 
colormap cool; %colormap est une fonction permettant d'avoir plusieurs coloris du graphe 
axis equal; 
 
hold off; 
end

Fonction variation de Thêta (pour obtenir la courbe A0=f(), pour h=0) :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
function variation_theta () 
h = 0; 
kkk = 0; 
for theta_deg=0:1:178 
   kkk = kkk +1; 
   [z,Y] =ptcap2D(h,theta_deg); 
   tab_theta(kkk,1)=theta_deg; % on stock les valeurs de theta_deg dans un tableau 
   Aire(kkk,1)=pi*Y(1,1)^2; 
   %pi*Y(1,1)^2 
end 
hold on 
plot(tab_theta,Aire,,'+r'); 
hold off 
end

Fonction variation de h (pour obtenir la courbe A=f(h) pour différents  fixé) :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
function variation_h() 
num_theta=0 ; 
for theta_deg=0:20:180 
   %faire varier h avec un thêta fixe qui se trouve dans argument grâce à un for 
% ici on varie thêta entre 0° et 180° 
num_theta= num_theta+1 ; 
hpas=0.5; 
   %Aire=zeros(100,1); 
   %tab_h=zeros(100,1); 
   jjj=0; 
   for jjj=0:1:1 
       kkk=0; 
       for h=0:hpas:100 
          kkk=kkk+1; 
          [z,Y] =ptcap2D(h,theta_deg); 
          %h 
          tab_h(kkk, num_theta)=h; 
          Aire(kkk, num_theta)=pi*Y(1,1)^2; 
          %pi*Y(1,1)^2 
          taille=size(z); 
          zfinal=z(taille(1)); 
          if(zfinal~=0) 
              if(kkk<=30) 
 
             hpas=hpas/10; 
             break; 
          else 
             break; 
          end 
       end 
        %utiliser ptcap2D pour avoir dernière valeur du tableau z qui correspond au pont capillaire 
       %si cette valeur est différente de 0 on arrête avec un break 
       %On doit prendre la première valeur de Y pour calculer l'aire de contact 
       %pi*r2=pi*Y2 => aire de contact 
       %tracer a(h) 
       if(kkk>30) 
          break; 
       end 
    end 
  end 
  hold on 
  plot(tab_h,Aire,'+r'); 
  hold off 
end 
end