本来是去回答百度知道里的问题的，自己技术不很够，折腾几个小时才搞明白，提交时已经有大神解决了问题了。。。。

是这个问题：

然后又悲剧地发现自己的结果和大神的不一样。。。。

公式是：

我得承认大神的代码比我的还是简洁多了，写函数时用了arrayFun，好简洁的。我是一个个写的。。。。

但是觉得我做出来的这个还是有点特点的，所以贴出来，和大家交流一下：

我的思路是：

这个问题是在一定的xyz范围内，对每一个点(x,y,z)生成对应的dB并积分，然后作出图像。图中的那个积分式，对于每个(x,y,z)点都要计算一次，由它决定锥形的方向与大小(积分结果是一个向量，当然就决定了方向，然后其模长就是锥形的大小)。

我采用的数值如下：

R=0.1

I=0.1

x,y,z都在[-15,15]区间内，间隔为2.

图像是这样的，我觉得这个图好神奇，照了很多个角度：

代码如下：

clc

clear all

close all

R=0.1;

I=0.1;

k=10^(-7);%k=u0/(4*pi)

[x,y,z]=meshgrid(-15:2:15);

u_fun=@(f)k.*R.*I.*z.*cos(f)./(sqrt((x-R.*cos(f)).^2+(y-R.*sin(f)).^2+z.^2)).^3;

v_fun=@(f)k*R*I.*z*sin(f)./(sqrt((x-R*cos(f)).^2+(y-R*sin(f)).^2+z.^2)).^3;

w_fun=@(f)-k*R*I*(sin(f)*(y-R*sin(f))+cos(f)*(x-R*cos(f)))./(sqrt((x-R*cos(f)).^2+(y-R*sin(f)).^2+z.^2)).^3;

u=integral(u_fun,0,2*pi,'ArrayValued',true);

v=integral(@(f)v_fun(f),0,2*pi,'ArrayValued',true);

w=integral(@(f)w_fun(f),0,2*pi,'ArrayValued',true);

[cx,cy,cz]=meshgrid([-15:2:15]);

h=coneplot(x,y,z,u,v,w,cx,cy,cz,5);

set(h,'FaceColor','r','EdgeColor','none');

camlight;lighting gouraud;

grid on;box on;

大神的结果是：

代码是：

L = R;

I = 0.1; R = 0.1;

u0 = pi*4e-7;

dB = @(x,y,z,t) u0/4*pi * [ R*I*z*cos(t); R*I*z*sin(t); ...

R*I*(sin(t).*(y-R*sin(t))+cos(t).*(x-R*cos(t))) ] / ...

sqrt( (x-R*cos(t)).^2 + (y-R*sin(t)).^2 + z.^2 ).^3;

B = @(x,y,z) quadv(@(t)dB(x,y,z,t),0,2*pi);

[x,y,z] = meshgrid(linspace(-L,L,10));

B = arrayfun(B,x,y,z,'UniformOutput',false);

[m,n,l] = size(B);

[m,n,l] = meshgrid(1:m,1:n,1:l);

Bx = arrayfun(@(i,j,k)B{i,j,k}(1),m,n,l);

By = arrayfun(@(i,j,k)B{i,j,k}(2),m,n,l);

Bz = arrayfun(@(i,j,k)B{i,j,k}(3),m,n,l);

B = sqrt(Bx.^2+By.^2+Bz.^2);

clf reset

colordef(gcf,'black')

[cx cy cz] = meshgrid(linspace(-L,L,10));

h=coneplot(x,y,z,Bx,By,Bz,cx,cy,cz,B,0.8);

set(h,'EdgeColor', 'none');

view(35,25)

xlabel('x');

ylabel('y');

zlabel('z');

box on;

axis tight

camlight left;

camlight right;

lighting phong

关于用coneplot画图 - MATLAB 基础讨论 - MATLAB中文论坛 - Powered by Discuz!

http://www.ilovematlab.cn/thread-93141-1-1.html

MIT向量场绘图教程