一个matlab制作的水滴模型

基于malab，持续不断的有水滴滴下来，在平面上产生一圈一圈的波纹

function waterwave

n = 100;                 
g = 9.8;                
dt = 0.02;             
dx = 1.0;
dy = 1.0;
nplotstep = 8;           
dropstep = 600;         
D = drop（121）;    

[surfplot] = initial（n）;
   
   H = ones（n+2n+2）;   U = zeros（n+2n+2）;  V = zeros（n+2n+2）;
   Hx = zeros（n+1n+1）; Ux = zeros（n+1n+1）; Vx = zeros（n+1n+1）;
   Hy = zeros（n+1n+1）; Uy = zeros（n+1n+1）; Vy = zeros（n+1n+1）;
   nstep = 0;
   
   while true
      
       if mod（nstepdropstep） == 0
           w = size（D1）;
           i = ceil（50）+（1:w）;
           j = ceil（50）+（1:w）;
           H（ij） = H（ij） + D;
       end
       
       
       nstep = nstep + 1;
       i = 1:n+1;
       j = 1:n;

       Hx（ij） = （H（i+1j+1）+H（ij+1））/2 - dt/（2*dx）*（U（i+1j+1）-U（ij+1））;
   
       Ux（ij） = （U（i+1j+1）+U（ij+1））/2 -  ...
                 dt/（2*dx）*（（U（i+1j+1）.^2./H（i+1j+1） + g/2*H（i+1j+1）.^2） - ...
                            （U（ij+1）.^2./H（ij+1） + g/2*H（ij+1）.^2））;

       Vx（ij） = （V（i+1j+1）+V（ij+1））/2 - ...
                 dt/（2*dx）*（（U（i+1j+1）.*V（i+1j+1）./H（i+1j+1）） - ...
                            （U（ij+1）.*V（ij+1）./H（ij+1）））;
       

       i = 1:n;
       j = 1:n+1;
   

       Hy（ij） = （H（i+1j+1）+H（i+1j））/2 - dt/（2*dy）*（V（i+1j+1）-V（i+1j））;
   

       Uy（ij） = （U（i+1j+1）+U（i+1j））/2 - ...
                 dt/（2*dy）*（（V（i+1j+1）.*U（i+1j+1）./H（i+1j+1）） - ...
                            （V（i+1j）.*U（i+1j）./H（i+1j）））;

       Vy（ij） = （V（i+1j+1）+V（i+1j））/2 - ...
                 dt/（2*dy）*（（V（i+1j+1）.^2./H（i+1j+1） + g/2*H（i+1j+1）.^2） - ...
                            （V（i+1j）.^2./H（i+1j） + g/2*H（i+1j）.^2））;

       i = 2:n+1;
       j = 2:n+1;

       H（ij） = H（ij） - （dt/dx）*（Ux（ij-1）-Ux（i-1j-1）） - ...
                         （dt/dy）*（Vy（i-1j）-Vy（i-1j-1））;
       U（ij） = U（ij） - （dt/dx）*（（Ux（ij-1）.^2./Hx（ij-1） + g/2*Hx（ij-1）.^2） - ...
                         （Ux（i-1j-1）.^2./Hx（i-1j-1） + g/2*Hx（i-1j-1）.^2）） ...
                       - （dt/dy）*（（Vy（i-1j）.*Uy（i-1j）./Hy（i-1j）） - ...
                         （Vy（i-1j-1）.*Uy（i-1j-1）./Hy（i-1j-1）））;
       V（ij） = V（ij） - （dt/dx）*（（Ux（ij-1）.*Vx（ij-1）./Hx（ij-1）） - ...
                         （Ux（i-1j-1）.*Vx（i-1j-1）./Hx（i-1j-1））） ...
                       - （dt/dy）*（（Vy（i-1j）.^2./Hy（i-1j） + g/2*Hy（i-1j）.^2） - ...
                         （Vy（i-1j-1）.^2./Hy（i-1j-1） + g/2*Hy（i-1j-1）.^2））;
   

      if mod（nstepnplotstep） == 0
          C = abs（U（ij）） + abs（V（ij））;  
          set（surfplot‘zdata‘H（ij）‘cdata‘C）;
          drawnow
     end
       if all（all（isnan（H））） break end  
   end
close（gcf）;


% ------------------------------------

function D = drop（heightwidth）

   [xy] = ndgrid（-1:（2/（width-1））:1）;
   D = height*exp（-5*（x.^2+y.^2））;

% ------------------------------------

function [surfplot] = initial（n）

   clf
   shg
   s

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件       3249  2012-10-26 01:01  waterwave.m

----------- ---------  ---------- -----  ----

                 3249                    1