[ÇóÖú] ÇóÖú£¬°ïÎÒ·ÖÎöϳÌÐòÄÄÀï³öÁËÎÊÌ⣿ÇóÖ¸½Ì

function H=Ha(nx,nz,hm,epslong,elta,N,w2,Ps,d,L)    rj=d/2; u2=2*pi*N/60; deltax=2*pi*rj/(nx-1) deltaz=L/(nz-1) delta=2*pi/(nx-1) h=zeros(nx+1,nz+1); a=zeros(nx+1,nz+1); b=zeros(nx+1,nz+1); c=zeros(nx+1,nz+1); d=zeros(nx+1,nz+1); e=zeros(nx+1,nz+1); f=zeros(nx+1,nz+1); ff=zeros(nx+1,nz+1); k=1; Ax=12*elta*u2*L/(hm^2*Ps); Az=12*elta*w2*L/(hm^2*Ps); for i=1:nx          theta=delta*(i-1);           h(i,=hm*(1+epslong*cos(theta+delta/2)); end S=0; T=0; ERR=1e-3; GAP=1; while GAP>ERR     k=k+1         for j=2:1:nz-1         if i==1            a(i,j)=h(i,j)^3/deltax^2-3*h(i,j)^2*(h(i+1,j)-h(nx,j))/(4*deltax^2);         b(i,j)=h(i,j)^3/deltax^2+3*h(i,j)^2*(h(i+1,j)-h(nx,j))/(4*deltax^2);         c(i,j)=-2*h(i,j)^3/deltax^2-2*h(i,j)^3/deltaz^2;         d(i,j)=h(i,j)^3/deltaz^2;         e(i,j)=h(i,j)^3/deltaz^2;         ff(i,j)=(a(i,j)*f(nx,j)+b(i,j)*f(i+1,j)+d(i,j)*f(i,j-1)+e(i,j)*f(i,j+1))/(-c(i,j));         elseif i==nx         a(i,j)=h(i,j)^3/deltax^2-3*h(i,j)^2*(h(1,j)-h(i-1,j))/(4*deltax^2);         b(i,j)=h(i,j)^3/deltax^2+3*h(i,j)^2*(h(1,j)-h(i-1,j))/(4*deltax^2);         c(i,j)=-2*h(i,j)^3/deltax^2-2*h(i,j)^3/deltaz^2;         d(i,j)=h(i,j)^3/deltaz^2;         e(i,j)=h(i,j)^3/deltaz^2;         pp(i,j)=(a(i,j)*f(i-1,j)+b(i,j)*f(1,j)+d(i,j)*f(i,j-1)+e(i,j)*f(i,j+1))/(-c(i,j));         else           if pp(i,j)<0                     pp(i,j)=0;                 else                 end            end         end         for i=2:1:nx-1         if j==1         a(i,j)=h(i,j)^3/deltax^2-3*h(i,j)^2*(h(i+1,j)-h(i-1,j))/(4*deltax^2);         b(i,j)=h(i,j)^3/deltax^2+3*h(i,j)^2*(h(i+1,j)-h(i-1,j))/(4*deltax^2);         c(i,j)=-2*h(i,j)^3/deltax^2-2*h(i,j)^3/deltaz^2;         d(i,j)=h(i,j)^3/deltaz^2;         e(i,j)=h(i,j)^3/deltaz^2;         ff(i,j)=(a(i,j)*f(i-1,j)+b(i,j)*f(i+1,j)+d(i,j)*Ps^2+e(i,j)*f(i,j+1))/(-c(i,j));         elseif j==nz         a(i,j)=h(i,j)^3/deltax^2-3*h(i,j)^2*(h(i+1,j)-h(i-1,j))/(4*deltax^2);         b(i,j)=h(i,j)^3/deltax^2+3*h(i,j)^2*(h(i+1,j)-h(i-1,j))/(4*deltax^2);         c(i,j)=-2*h(i,j)^3/deltax^2-2*h(i,j)^3/deltaz^2;         d(i,j)=h(i,j)^3/deltaz^2;         e(i,j)=h(i,j)^3/deltaz^2;         ff(i,j)=(a(i,j)*f(i-1,j)+b(i,j)*f(i+1,j)+d(i,j)*f(i,j-1)+e(i,j)*f(i,j+1))/(-c(i,j));         else             if ff(i,j)<0                    ff(i,j)=0;                 else                 end         end         end         if i==1&&j==1            p(i,j)=Ps         elseif i==nx&&j==1             p(i,j)=Ps         elseif i==1&&j==nz         a(i,j)=h(i,j)^3/deltax^2-3*h(i,j)^2*(h(1,j)-h(i-1,j))/(4*deltax^2);         b(i,j)=h(i,j)^3/deltax^2+3*h(i,j)^2*(h(1,j)-h(i-1,j))/(4*deltax^2);         c(i,j)=-2*h(i,j)^3/deltax^2-2*h(i,j)^3/deltaz^2;         d(i,j)=h(i,j)^3/deltaz^2;         e(i,j)=h(i,j)^3/deltaz^2;         ff(i,j)=(a(i,j)*f(nx,j)+b(i,j)*f(i+1,j)+d(i,j)*f(i,j-1)+e(i,j)*f(i,nz))/(-c(i,j));           elseif i==nx&&j==nz         a(i,j)=h(i,j)^3/deltax^2-3*h(i,j)^2*(h(1,j)-h(i-1,j))/(4*deltax^2);         b(i,j)=h(i,j)^3/deltax^2+3*h(i,j)^2*(h(1,j)-h(i-1,j))/(4*deltax^2);         c(i,j)=-2*h(i,j)^3/deltax^2-2*h(i,j)^3/deltaz^2;         d(i,j)=h(i,j)^3/deltaz^2;         e(i,j)=h(i,j)^3/deltaz^2;         ff(i,j)=(a(i,j)*f(i-1,j)+b(i,j)*f(1,j)+d(i,j)*f(i,j-1)+e(i,j)*f(i,nz-1))/(-c(i,j));            else              if ff(i,j)<0                     ff(i,j)=0;                 else                 end              end   for i=2:1:nx     for j=2:1:nz             S=sum(sum(abs(ff-f)));             T=sum(sum(abs(f)));      end    end     GAP=S/T;     f=ff end

»Ø¸´´ËÂ¥