第6章 在普通物理中的应用【例6-1-1】温度单位转换命题：写出一个程序，能把用户输入的摄氏温度转为华氏，也可反求。解：建模两种温度之间的转换公式为：摄氏变华氏 华氏变摄氏 程序中要先考虑由用户选择转换的方向，再给数据。程序exn611k = input(选择1：摄氏变华氏；选择2：华氏变摄氏；键入数字1或2：);Tin = input(输入待变换的温度(允许输入数组):);if k=1 Tout = Tin*9/5 +32; % 摄氏转华氏elseif k=2 Tout = (Tin-32)*5/9; % 华氏转摄氏else disp(未给转换方向,转换无效),ends = 华氏;摄氏;s1 = 转换后的温度为,s(k,:),num2str(Tout),度, % 注意此语句的编写方法【例6-1-2】多种单位间的换算写出一个程序，能把用户输入的长度单位在厘米、米、千米、英寸、英尺、英里、市尺、市里之间任意转换。解：建模这里采取的技巧是分成两步，先把输入量变换为米，第二步再把米变换为输出单位，另外，把变换常数直接表示为一个数组，选择单位的序号也就成了数组的下标；这样程序就比较简明易读。程序如下： 长度单位换算程序ex612.mclear all; disp( 长度单位换算程序) fprintf(长度单位: n); % 选择输入输出的单位fprintf( 1) 厘米 2) 米 3) 千米 4) 英寸 n);fprintf( 5) 英尺 6) 英哩 7) 市尺 8) 市里 n);InUnits = input(选择输入单位编号: );OutUnits = input(选择输出单位编号: );% 令各种单位对米的变换常数数组为ToMeter ToMeter = 0.01, 1.00, 1000.0, 0.0254, 0.3048, 1609.3, 1/3, 500 ;程序ex612.m(续)FrmMeter= 1./ ToMeter; Value = input(输入待变换的值(0为退出): ); while( Value = 0 ) ValueinM = Value*ToMeter(InUnits); % 把输入值变为米 NewValue = ValueinM*FrmMeter(OutUnits); % 把米变为输出单位 fprintf(变换后的值是 %g n,NewValue); % 打印变换后的值 Value = input(输入待变换的值(0为退出): ); % 提问下个输入值end【例6-1-3】实验数据拟合命题：命题：设在某一实验中，给某元件加1,2,3,4,5v电压，测得的电流为 0.2339, 0.3812, 0.5759, 0.8153, 0.9742 ma。求此元件的电阻。解：解： 建模建模模型:设直线的方程为 y= a(1)x +a(2), 待定的系数是a(1),a(2)。将上述数据分别代入x,y,a(1)+a(2)=0.23392a(1)+a(2)=0.38123a(1)+a(2)=0.57594a(1)+a(2)=0.81535a(1)+a(2)=0.9742把这五个方程联立,用矩阵表述, 得 datax *a(1) + ones(N,1) * a(2) = datay 程序exn613这是一个超定方程组，写成A*a = B，其最小二乘解可以用左除运算符a = A B来求得。因此程序如下：lear, datax = 1:5;datay= 0.2339, 0.3812, 0.5759, 0.8153, 0.9742A = datax , ones(5,1); B = datay; a = AB, r=1/a(1) plot(datax,datay,o),hold on 运行结果为：a(1) = 0.1905a(2) = 0.0247画出曲线如右图。6.2 力学基础力学基础【例6-2-1】目标相对于射点的高度为 yf,给定初速和射角,计算物体在真空中飞行的时间和距离。 建模：这里目标和射点不在同一高度上，不好求封闭形式的解，用MATLAB使整个计算和绘图过程自动化。其好处是可快速地计算其在不同初速和射角下的飞行时间和距离。关键在求落点时间tf时需要解一个二次线性代数方程由 解出t，她就是落点时间tf。它会有两个解，我们只取其中一个有效解。再求程序exn621clear; y0 = 0; x0 = 0; % 初始位置vMag = input(输入初始速度 (m/s): ); % 输入初始速度vDir = input( 输入初速方向(度): );yf = input(输入目标高度（米）: ); % 输入目标高度yhvx0 = vMag*cos(vDir* (pi/180); % 计算x,y方向的初速vy0 = vMag*sin(vDir* (pi/180); % wy = -9.81; wx = 0; % 重力加速度 (m/s2)tf=roots(wy/2,vy0,y0-yf); % 解代数方程计算落点tftf=max(tf); % 去除tf两个解中的庸解t=0:0.1:tf;y = y0 + vy0*t + wy*t.2/2; % 计算轨迹x = x0 + vx0*t + wx*t.2/2;xf = max(x),plot(x,y), % 计算射程，画出轨迹例6-2-1运行结果在检查曲线正确后,键入hold命令,把曲线保留下来，以便用同样的初速，不同的射角，比较其曲线和最大射程。运行结果输入初始速度 (m/s): 50, 输入初速方向(度): 40输入目标高度（米）: 8得xf = 237.4738而初速方向为50度时，xf = 241.0454所得曲线见图6-2-1.例6-2-2 质点的平面运动 给定质点沿x和y两方向的运动规律 x(t) 和 y(t),求其运动轨迹,并计算其对原点的角动量。解解:建模：由用户输入解析表示式需要用到字符串的输入语句,其第二变元为s,而运行这个字符串要用eval命令.当x(t) 和 y(t)都是周期运动时,所得的曲线就是李萨如图形.动量矩等于动量与向径的叉乘(cross product).求速度需要用导数,可用MATLAB的diff函数作近似导数计算。设角动量为 ，质点的动量为 ，向径为 ，则在x-y平面上的投影为程序exn622x = input(: ,s); y = input(: ,s); % 读入字符串tf = input( tf= ); Ns=100; t=linspace(0,tf,Ns); dt=tf/(Ns-1); % 分Ns个点,求出时间增量dtxPlot=eval(x);yPlot=eval(y);% 计算各点x(t), y(t)的近似导数和角动量。 p_x = diff(xPlot)/dt; % p_x = M dx/dt p_y = diff(yPlot)/dt; % p_y = M dy/dt %求角动量 LPlot = xPlot(1:Ns-1).* p_y - yPlot(1:Ns-1).* p_x;% 画出轨迹及角动量随时间变化的曲线程序运行结果运行此程序，输入x=t.*cos(t)y=t.*sin(t)tf=20后，得出图6-2-2。如果输入x=cos(2*t)y=sin(3*t)图6-2 按方程x=tcos(t),y=tsin(t)画出轨迹及角动量曲线例6-2-3 质点系的动力学物体A(质量为m1)在具有斜面的物体B(质量为m2)上靠重力下滑, 设斜面和地面均物摩擦力,求A沿斜面下滑的相对加速度a1和B的加速度a2,并求斜面和地面的支撑力N1及N2。解：建模，对物体A，列出方程对物体B，列出方程方程组的矩阵建模四个方程包含四个未知数，将含未知数的项移到等式左边，常数项移到等式右端，得到矩阵方程于是有 X=AB程序exn623m1=input(m1=【千克】 );m2=input(m2=【千克】 );theta=input(theta【度】= );theta=theta*pi/180; g=9.81;A = m1*cos(theta),-m1, -sin(theta), 0;. m1*sin(theta), 0, cos(theta), 0;. 0 , m2, -sin(theta), 0;. 0 , 0, -cos(theta), 1 ;B = 0, m1*g, 0, m2*g; X=AB;a1=X(1),a2=X(2),N1=X(3),N2=X(4)运行结果输入m1=2【kg】,m2=4【kg】,及theta=30【deg】,得到a1 = 6.5400【m/s2】; a2 = 1.8879【m/s2】N1 = 15.1035【N】; N2 = 52.3200【N】静力学平衡和动力学中求力与加速度关系的问题,通常都可归结为线性方程组的求解,只要方程组列写正确,用MATLAB的矩阵除法就可以方便而准确的求出其解.例6-2-4 碰撞问题质量为m的小球以速度u0正面撞击质量为M的静止小球,假设碰撞是完全弹性的,即没有能量损失,求碰撞后两球的速度,及它们与两球质量比K=M/m的关系.解: 建模 设碰撞后两球速度都与u0同向,球m的速度为u, 球M的速度为v,列出动量守恒和能量守恒方程,则引入质量比K=M/m和相对速度ur= u/u0,vr=v/u0后, 有动量守恒mu0=mu+Mv动能守恒化为碰撞问题的方程由（3）（5）代入（4）（6）主动球的能量损失为展开并整理多项式(6)，得可用roots命令求根,程序exn624clearK=logspace(-1,1,11);%设自变量数组K,从K=0.1到10,按等比取for i=1:length(K)% 对各个K循环计算ur1=roots(1+1/K(i),-2/K(i),(1/K(i)-1);% 二次方程有两个解ur(i)=ur1(abs(ur1-1)0.001);% 去掉在1邻近的庸解endvr=(1-ur)./K;% 用（5）式求vr,用元素群运算em=1-ur.*ur;% 主动球损失的相对能量K,ur,vr,em% 显示输出数据semilogx(K,ur,vr,em),grid%绘图程序运行结果数字结果为(省略了几行) Kur vr em0.1000 0.8182 1.8182 0.33060.3981 0.4305 1.4305 0.81471.0000 0 1.0000 1.00002.5119 -0.4305 0.5695 0.814710.000 -0.8182 0.1818 0.3306绘出的曲线见图6-2-4.可以看出,当K1时,ur为负,即当静止球质量大于主动球质量时,主动球将产生回弹.K=1时ur=0, 即主动球将全部动能传给静止球. K1时,ur为正,说明主动球将继续沿原来方向运动.例6-3-1 麦克斯韦速度分布律命题：求摄氏27度下氮气的分子运动速度分布率，并求速度在300500m/s范围内的分子所占的比例，讨论温度T及分子量mu对速度分布曲线的影响。解:建模 麦克斯韦速度分布律为本例将说明如何从复杂数学公式中绘制曲线并研究单个参数的影响.先把麦克斯韦速度分布律列成一个子程序,以便经常调用,把一些常用的常数也放在其中,主程序就简单了.麦克斯韦分布律子程序mxwlfunction f=mxwl(T,mu,v)mu 分子量,公斤.摩尔-1（如氮为28103）v 分子速度（可以是一个数组）T 气体的绝对温度R=8.31;%气体常数k=1.381*10(-23);%玻尔茨曼常数NA=6.022*1023;%阿伏伽德罗数m=mu/NA;%分子质量%麦克斯韦分布率f=4*pi*(m/(2*pi*k*T).(3/2).*exp(-m*v.2./(2*k*T).*v.*v;主程序exn631T=300;mu=28e-3; % 给出T,muv=0:1500;% 给出自变量数组y=mxwl(T,mu,v); % 调用子程序plot(v,y),hold on % 画出分布曲线v1=300:500;% 给定速度范围y1=mxwl(T,mu,v1);% 该范围的分布fill(v1,500,300,y1,0,0,r)trapz(y1)% 求该范围概率积分执行此程序所得曲线积分结果为：ns = 0.3763可在程序中再加几句：% 改变T，画曲线T=200;mu=28e-3;y=mxwl(T,mu,v);plot(v,y)% 改变mu, 画曲线T=300;mu=2e-3;y=mxwl(T,mu,v);plot(v,y)可见减小T，使分子的速度分布向低端移动；减小分子量mu，使速度分布向高端移动；这是与物理概念相一致的。 麦克斯韦分布曲线6.4 静电场【例6-4-1】设电荷均匀分布于从z=-L到z=L通过原点的线段上，其密度为q库仑/米，求出在xy平面上的电位分布。解：建模 点电荷产生的电位可表为V = Q/4r0其中r为电荷到测量点的距离.线电荷所产生的电位可用积分或叠加的方法来求。为此把线电荷分为长为dL的N段(在MATLAB中, dL应理解为L)。每段上电荷为q*dL.它产生的的电位为 然后对全部电荷求和即可。程序exn641E0 = 8.85e-12; % 真空电介质常数C0 = 1/4/pi/E0;% 归并常数L0 = linspace(-L,L,N+1);% 将线电荷分N段L1 = L0(1:N); L2=L0(2:N+1);% 确定每段的起点和终点Lm = (L1+L2)/2;dL= 2*L/N;% 每段的中点和长度的数组R = linspace(0,10,Nr+1);% 将R分N+1点for k = 1:Nr+1% 对R的N+1点循环计算Rk = sqrt(Lm.2+R(k)2);% 测量点到电荷段的向径长度Vk = C0*dL*q./Rk;% 第k个电荷段产生的电位V(k) = sum(Vk);% 对各电荷段产生的电位求和Endplot(R,V),grid 程序运行结果(1)。q=1,L=5, N=50, Nr=50(2)。q=1,L=50, N=500,Nr=50所得结果为：电场的最大最小值：（1）1.0e+010* 9.3199 ,0.8654（2）1.0e+011 * 1.3461, 0.4159沿R的电场分布见图6-4-1，上图为半对数坐标，下图为线性坐标。 线电荷产生的静电场分布6-4-2 由电位的表示式计算电场已知空间的电位分布，画出等电位线和电场方向.解: 建模如果已知空间的电位分布 V=V(x,y,z),则空间的电场等于电位场的负梯度其中 分别为x,y,z三个方向的单位向量。MATLAB中设有gradient函数，它是靠数值微分的，因此空间的观测点应取得密一些，以获得较高的精度。程序exn642V = input(例如: log(x.2 + y.2): ,s); % 读入字符串,xMax = 5; NGrid = 20; % 绘图区从 x= -xMax 到 x= xMax,网格线数xPlot = linspace( -xMax, xMax, NGrid); % x,y取同样范围,生成二维网格x,y=meshgrid(xPlot);% 按给定的x,y执行输入的字符串VVPlot=eval(V); %电场等于电位的负梯度ExPlot, EyPlot = gradient(-VPlot);% 画出含等高线的三维曲面clf; subplot(1,2,1),meshc(VPlot); 程序exn642(续)% 规定等高线图的范围及比例% 建立第二子图subplot(1,2,2), axis(-xMax xMax -xMax xMax); % 画等高线,cs是等高线值，并加上编号cs = contour(x,y,VPlot); clabel(cs); hold on; % 在等高线图上加上电场方向quiver(x,y,ExPlot,EyPlot); % 画电场 E 的箭头图xlabel(x); ylabel(y);hold off; 运行 在输入电位方程V(x,y) = log(x.2 + y.2)时,得出图6-4-2左的电位分布曲面,右面是电场分布的向量图.Exn642的运行结果图6-4-2 V(x,y) = log(x.2 + y.2)的电位三维立体图，等高线及电场分布图6-5 恒稳磁场恒稳磁场例6-5-1 用毕奥萨伐定律计算电流环产生的磁场解:建模载流导线产生磁场的基本规律为，任一电流元 在空间任一点P处所产生的磁感应强度 为下列向量叉乘积，即其中， 为电流元到P点的矢径， 为导线元的长度矢量，P点的总磁场可沿载流导体全长积分各段产生的磁场来求得。程序xn651x=linspace(-3, 3, 20); y=x; % 确定观测点的x,y坐标数组Nh = 20;% 电流环分段数 % 计算每段的端点,环在x=0平面上,其坐标x1,x2均为零 theta0 = linspace(0,2*pi, Nh+1); % 环的圆周角分段 theta1 = theta0(1:Nh);% 注意theta1和theta2的差别 y1 = Rh*cos(theta1);% 环各段向量的起点坐标y1,z1 z1 = Rh*sin(theta1); theta2 = theta0(2:Nh+1); y2 = Rh*cos(theta2);% 环各段向量的终点坐标y2,z2 z2 = Rh*sin(theta2); dlx=0; dly = y2-y1; dlz = z2-z1; %计算dl的长度分量 xc=0; yc=(y2+y1)/2; zc=(z2+z1)/2; %各段中点的坐标分量程序xn651(续)% 循环计算各网格点上的B(x,y) 值for i=1:NGy for j=1:NGx % 对yz平面内的电流环分段作元素群运算,先算环上某段与观测点之间的向量r rx=x(j)-xc; ry=y(i)-yc; rz=0-zc; % 观测点在z=0平面上 r3 = sqrt(rx.2 + ry.2 + rz.2).3; % 计算r3 dlXr_x = dly.*rz - dlz.*ry; % 计算叉乘积 dlXr_y = dlz.*rx - dlx.*rz; % 把环各段产生的磁场分量累加Bx(i,j) = sum(C0*dlXr_x./r3); By(i,j) = sum(C0*dlXr_y./r3); end, end程序运行结果% 用quiver 画磁场向量图clf; quiver(x,y,Bx,By); 图形标注语句及在图上画出圆环位置的语句略运行此程序所得图形见图6-5-1,读者可改变电流环的直径来分析其影响，也可加上显示各点磁场强度的语句来分析其强度的分布。图6-5-1 电流环产生的磁场分布图例6-5-2 亥姆霍兹线圈一对相同的共轴的彼此平行的载流圆线圈,当它们的间距正好等于其线圈半径时,称之为亥姆霍兹线圈.计算表明, 亥姆霍兹线圈轴线附近的磁场是十分均匀的,而且都沿x方向.本题要求对这一论断进行验证。解: 建模 本题的计算模型与上例相同,只是把观测区域取在两线圈之间的小范围内，如右图所示。程序exn652clear all; % 初始化(给定环半径,电流,图形 )mu0 = 4*pi*1e-7; % 真空导磁率 (T*m/A)I0 = 5.0; Rh=1; %,在本题中不影响结果C0 = mu0/(4*pi) * I0; % 归并常数% 下面三行输入语句与上题不同, 观测范围x取-Rh,Rh,即线圈的左右都取,因为以后要% 把第一个线圈的右磁场与第二个线圈的左边磁场叠加, y也取-Rh,Rh.NGx =21 ;NGy = 21;% 设定观测点网格数x=linspace(-Rh,Rh, NGx); % 设定观测点范围及数组y=linspace(-Rh,Rh, NGy);程序exn652(续)主程序段同例主程序段同例6-5-1(从从Nh到最后一个到最后一个end)后两行是与上例不同的输出绘图语句% 把x0的区域Bax=Bx(:,11:21)+Bx(:,1:11); % 模仿右边线圈所增加的磁场Bay=By(:,11:21)+By(:,1:11);subplot(1,2,1),% 画出其Bx分布三维图mesh(x(11:21),y,Bax);xlabel(x);ylabel(y);subplot(1,2,2),plot(y,Bax),grid,xlabel(y);ylabel(Bx);程序exn652运行结果6-6 振动与波振动与波例6-6-1 振动的合成及拍频现象解:建模，将两个同方向的振动相加，可得当 很接近时， 成为一个很低的频率，称为拍频。用MATLAB程序得到的图形和声音中可以很清楚地看出拍频现象。程序exn661t=0:0.001:10;%10秒钟，分10000个点%输入两组信号的振幅和频率a1=input(振幅1=); w1=input(频率1=); a2=input(振幅2=); w2=input(频率2=); y1=a1*sin(w1*t); %生成两个正弦波y2=a2*sin(w2*t);y=y1+y2;%将两个波叠加subplot(3,1,1),plot(t,y1),ylabel(y1)%画出曲线subplot(3,1,2),plot(t,y2),ylabel(y2)subplot(3,1,3),plot(t,y),ylabel(y),xlabel(t)程序exn661运行结果pause %产生声音sound(y1);pause(2), sound(y2);pause(2),sound(y),pause程序运行结果:按a1=1.2; w1=300a2=1.8; w2=310运行的结果见图6-6-1，由于两个频率非常接近，产生了差拍频率.图6-6-1 拍频现象例6-6-2 多普勒效应例6-2-2 设声源从500m外以50m/s的速度对听者直线开来,其轨迹与听者的最小垂直距离为20m.声源的角频率为1000弧度/s,试求出听者接收到的信号波形方程并生成其对应的声音。解: 建模设声源发出的信号为f(t)，传到听者处，被听者接收的信号经历了声音传播迟延，迟延时间为 其中c为音速，为声源与听者之间的距离。被接收的信号形式为（不考虑声波的传播衰减） 。只要给出随t变化的关系，即可求得并将它恢复为声音信号。程序exn662x0=500;v=60;y0=30;% 设定声源运动参数c=330;w=1000;% 音速和频率t=0:0.001:30;%设定时间数组r=sqrt(x0-v*t).2+y0.2);% 计算声源与听者距离t1=t-r/c;% 经距离迟延后听者的等效时间u=sin(w*t)+sin(1.1*w*t); % 声源发出的信号u1=sin(w*t1)+sin(1.1*w*t1); % 听者接受到的信号% 先后将原信号和接受到的信号恢复为声音sound(u);pause(5);sound(u1);程序exn662运行结果打开计算机的声音系统,运行此程序将会听到类似于火车汽笛的声音.第一声是火车静止时的汽笛声, 第二声是本题中听者听到的运动火车的汽笛声,它的频率先高于,后低于原来的汽笛声. 程序中两个sound语句之间加的pause语句是不可少的,而且暂停的时间要足够长，以便再打开声音系统，这个量与计算机硬件有关，在本书中用5。6.7 光学【例6-7-1】单色光通过两个窄缝射向屏幕，相当于位置不同的两个同频同相光源向屏幕照射的叠合，由于到达屏幕各点的距离（光程）不同引起相位差，如图6-7-1所示。叠合的结果在有的点加强，在有的点抵消，造成干涉现象。考虑到纯粹的单色光不易获得，通常都有一定的光谱宽度，它对光的 干涉会产生何站种效应，要求用MATLAB计算并仿真这一问题。单色光双缝干涉模型解:建模 考虑两个离中心点距离各为d/2的相干光源S1和S2到屏幕上任意点的距离差引起的相位差,先分析光程则光程差为L=L1-L2将L除以波长，并乘以2，得到相位差 。设两束相干光在屏幕上产生的幅度相同,均为A0,则夹角为的两个向量A0的合成向量的幅度为:A = 2 A0 cos(/2)光强B正比于振幅的平方,故有:B = 4 B0 cos2(/2)根据这些关系式,可以编写出计算屏幕上各点光强的程序,程序exn671输入波长Lambda=500nm,光缝距离d=2mm,光栅到屏幕距离Z=1myMax = 5*Lambda*Z/d; xs = yMax; % 设定图案的y,x 向范围Ny=101;ys = linspace(-yMax,yMax,Ny); % y方向分成101点for i=1:Ny% 对屏上全部点进行循环计算 % 计算第一和第二个光源到屏上各点的距离 L1 = sqrt(ys(i)-d/2).2 + Z2 ); L2 = sqrt(ys(i)+d/2).2 + Z2 ); Phi = 2*pi*(L2-L1)/Lambda; % 从距离差计算相位差 B(i,:) = 4*cos(Phi/2).2; % 计算该点光强（设两束光强相同)endclf; figure(gcf); % 清图形窗，将它移到前面,准备绘图NCLevels = 255; % 确定用的灰度等级% 定标：使最大光强(4.0)对应于最大灰度级（白色）Br = (B/4.0) * NCLevels;subplot(1,2,1),image(xs,ys,Br);%画图象程序exn671运行结果运行exn671和程序所得的屏幕光强图像见图6-7-2。光的非单色性导致干涉现象的减弱。光谱很宽的光将不能形成干涉。 图6-7-2单色光双缝干涉条纹及光强分布考虑光的非单色性再研究复杂一些的问题,考虑到光的非单色性对干涉条纹的影响。此时波长将不是常数，必须对不同波长的光作分类处理再叠加起来。假定光源的光谱宽度为中心波长的正负10%，并且在该区域内均匀分布。在(0.91.1)之间，按均等间距近似取11根谱线，其波长分别为则上面求相位差的计算式求出的将是对不同谱线的11个不同相位。计算光强时应把这11根谱线产生的光强迭加取平均值,即程序exn671a运行结果 则在原程序exn671中Phi和B(I,:)两句程序要换成以下四句： Nl=11; dL=linspace(-0.1,0.1,Nl);%设光谱相对宽度正负10%, Lambda1=Lambda*(1+dL);%分11根谱线,波长为一个数组 Phi1 = 2*pi*(L2-L1)./Lambda1; % 从距离差计算各波长的相位差 B(i,:) = sum(4*cos(Phi1/2).2)/Nl; % 叠加各波长影响计算光强其他不变，运行结果见右图。不纯光双缝干涉条纹及光强分布可以看出，光的非单色性导致干涉现象的减弱。光谱很宽的光将不能形成干涉。例6-7-2 光的衍射模型把单色平行光通过的光缝当作N点干涉来计算。把单缝看作一排NPoints个等间隔光源，分布在-a/2+a/2区间内。若屏幕离光源的距离为z，取屏幕上一点ys，该处的光强应为这NPoints个光源照射结果的合成。设光源的坐标为yPoint，它是一个数组，其长度为NPoints，起点是-a/2，终点为a/2，用MATLAB语句可表为yPoint = linspace(-a/2,a/2,NPoints)。在则它到屏幕点ys的路程L也是一个同样长度的数组，其计算公式为编程的思路以光源的中心到屏幕的中心之间的距离z为基准，则其光程差等于L-z，对应的光相位差为，它应该是ys和yPoint的函数。在计算时，先取定屏幕上一点ys，让yPoint取所有的值，把得到的Phi相加起来，就算出了屏幕上这点的光强。然后再换另一个屏幕点ys，再作循环。这样的双重循环，计算量是很大的，比如把光源数和屏幕点数都取51，则进行的循环将为2500次以上。用手工计算是很难想象的。下面就是按此思路写出的程序： 程序exn672按提示从键盘输入波长、缝宽a和距离Z的输入语句从略ymax = 3*Lambda*Z/a;%屏幕范围（沿y向）Ny = 51; % 屏幕上的点数（沿y向）ys = linspace(-ymax,ymax,Ny);NPoints = 51; % 缝上的点数（沿y向）yPoint = linspace(-a/2,a/2,NPoints); %把点数设成数组for j=1:Ny% 对屏幕上y向各点作循环% 对光缝中各点作循环,计算缝点到屏幕位置的距离 L = sqrt(ys(j)-yPoint).2 + Z2 ); % L是一个数组程序exn672(续) Phi = 2*pi.*(L-Z)./Lambda; % 对于屏幕中心的相位差,也是一个数组 % 求每个分量的累加和 SumCos = sum(cos(Phi); % 数组求和 SumSin = sum(sin(Phi); % 求屏幕上的归一化光强; B(j) = (SumCos2 + SumSin2)/NPoints2; endclf,plot(ys,B,*,ys,B);grid% 屏幕上光强与位置的关系曲线图形标注语句程序exn672运行结果:依次输入波长Lambda=500nm, 距离z=1M，缝宽为a=0.2,1,2mm三种情况，程序运行结果见图6-7-3. 三种情况统称费涅耳衍射,只有最左边的情况符合夫琅和费衍射的条件:*a2/(4*Z)1这个现象也适用于电磁波的发射,天线的设计和测量都要用这个概念. 天线探测目标时通常符合夫琅和费衍射的条件，形成天线的远场波瓣。在天线测量时却希望近些才方便，这时往往处于近场，于是要建立两者之间的转换关系，MATLAB程序可以发挥作用。程序exn672运行结果(图) 图6-7-3 缝宽a为0.2,1,2mm三种情况所得衍射光强曲线（左为夫琅和费衍射）