求助：高斯模型出图
 

气体高斯扩散模型 怎么调出图来啊 等浓度曲线图 等浓度区域图
调的时候我就是自己试一些数据但是得不到较好的图 请高人高人指教


function gauss()
global Q H u;
Q=2.5;
H=20;
u=2.5;
%Q=str2double(get(handles.input_Q,'String'));%获得计算的毒气量值
%u=str2double(get(handles.input_u,'String'));%获得用户输入的风速值
%H=str2double(get(handles.input_H,'String'));%获得用户输入的毒源有效高度值
%xmin=str2double(get(handles.input_xmin,'String'));%设置要显示的x轴的范围
%xmax=str2double(get(handles.input_xmax,'String'));
%ymin=str2double(get(handles.input_ymin,'String'));%设置要显示的y轴的范围
%ymax=str2double(get(handles.input_ymax,'String'));
%z=str2double(get(handles.input_z,'String')); %获得要显示浓度分布距地面的垂直高度值
%n=str2double(get(handles.input_n,'String')); %获得用户确定显示的坐标轴精度值
xmin=100;xmax=500;ymin=-250;ymax=250;z=183;n=100;
[x,y]=meshgrid(xmin:((xmax-xmin)/(n-1)):xmax,ymin:((ymax-ymin)/(n-1)):ymax);%网格划分
%创建扩散参数系数矩阵
C=[0.924279 0.177154 0.917595 0.106803;...
0.885157 0.232123 0 0 ];
C1=[0.926849 0.143940 0.838628 0.126152;...
0.886490 0.189396 0.756410 0.235667;...
0 0 0.815575 0.136659];
%s=get(handles.select,'Value'); %选择大气稳定度等级
s=4;
switch (s)
case 4 %如果选择大气稳定度为B1级
if x<=500
a1=B1(1,1);b1=B1(1,2);a2=B1(1,3);b2=B1(1,4);
elseif x>500&&x<=1000
a1=B1(1,1);b1=B1(1,2);a2=B1(2,3);b2=B1(2,4);
elseif x>1000
a1=B1(2,1);b1=B1(2,2);a2=B1(2,3);b2=B1(2,4);
end
case 5 %如果选择大气稳定度为C级
if x>=1&&x<=1000
a1=C(1,1);b1=C(1,2);a2=C(1,3);b2=C(1,4);
elseif x>1000
a1=C(2,1);b1=C(2,2);a2=C(1,3);b2=C(1,4);
end
end
y1=b1*x.^a1; %根据上述选取的扩散参数系数计算y、z轴上的扩散参数
z1=b2*x.^a2;
c=Q./(2*pi*u*y1.*z1).*exp(-y.^2./(2*y1.^2)).*(exp(-(z-H)^2./(2*z1.^2))+...
exp(-(z+H)^2./(2*z1.^2))); %高斯方程计算扩散浓度

%绘制某浓度的危险区域图程序
% c0=str2double(get(handles.input_c0,'String'));
c0=0.5;
[I,J]=find(abs(c-c0)<0.01); %限定误差范围误差
for ii=1:length(I)
c(I(ii),J(ii))=NaN;
end
surf(x,y,c);
colorbar;
shading interp
xlabel('x'),ylabel('y'),zlabel('c');title('特定浓度值的范围及浓度分布');
view(0,90); %查看俯视图

%
% %绘制多条等值浓度曲线程序
% cm=str2double(get(handles.cmatrix,'String')); %获得用户输入的浓度值
% a=contour(x,y,c,[cm cm],'r'); %绘制等浓度曲线
% clabel(a)
% xlabel('x'),ylabel('y'),title('特定浓度的等值线');
% grid on
% hold on
%
% c=str2double(get(handles.input_c1,'String')); %获得用户输入的浓度值
% gx=str2double(get(handles.input_guess,'String')); %通过查看等浓度曲线获得地面轴线下风向距离的近似值
% ……%用 gx代替以上扩散参数系数选择代码中的 x，确定系数 a1,b1,a2,b2的值
% p=sym('Q/(2*pi*u*(b1*x^a1)*(b2*x^a2))*(exp(-(-H)^2/(2*(b2*x^a2)^2))+exp(-H^2/(2*(b2*x^a2)^2)))-c=0'); %方程字符化
% p=simple(p); %化简方程
% x=solve(p); %求解字符方程
% x=subs(x);%求得某浓度下风向地面轴线距离 x的精确值
% x=num2str(x);
% set(handles.show_x1,'String',x);
%
% x0=str2num(get(handles.input_x0,'String')); %获得用户输入的下风向距离值
% x=x0;
% ……% 省略部分参照前面程序
% c=num2str(c); % 将计算出的浓度值转化为字符型
% set(handles.show_cx,'String',c); % 输出该处浓度值
这是气体高斯扩散的模型描述，前面一部分注释可以直接赋值计算，后面注释主要能调通“等浓度曲线”的部分就好了，非注释部分已调通，并出图，只是离理论形状差的有点远，问题可能出在坐标轴上。
和这个帖有些类似的： [url]http://www.chinavib.com/forum/attachments/month_0903/20090324_6f387074816aefe7fbc9BYVzlFNPdKlE.jpg[/url][IMG]http://www.chinavib.com/forum/attachments/month_0903/20090324_6f387074816aefe7fbc9BYVzlFNPdKlE.jpg[/IMG] 得到这样的效果就好