在普通物理中的应用省公共课一等奖全国赛课获奖课件.pptx

第6章 在普通物理中应用第1页【例6-1-1】温度单位转换命题：写出一个程序，能把用户输入摄氏温度转为华氏，也可反求。解：建模两种温度之间转换公式为：摄氏变华氏 华氏变摄氏 程序中要先考虑由用户选择转换方向，再给数据。第2页程序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),度,%注意此语句编写方法第3页【例6-1-2】各种单位间换算写出一个程序，能把用户输入长度单位在厘米、米、千米、英寸、英尺、英里、市尺、市里之间任意转换。解：建模这里采取技巧是分成两步，先把输入量变换为米，第二步再把米变换为输出单位，另外，把变换常数直接表示为一个数组，选择单位序号也就成了数组下标；这么程序就比较简明易读。程序以下：第4页长度单位换算程序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;第5页程序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页【例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 第7页程序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画出曲线如右图。第8页6.2 力学基础力学基础【例6-2-1】目标相对于射点高度为 yf,给定初速和射角,计算物体在真空中飞行时间和距离。建模：这里目标和射点不在同一高度上，不好求封闭形式解，用MATLAB使整个计算和绘图过程自动化。其好处是可快速地计算其在不一样初速和射角下飞行时间和距离。关键在求落点时间tf时需要解一个二次线性代数方程由 解出t，她就是落点时间tf。它会有两个解，我们只取其中一个有效解。再求第9页程序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),%计算射程，画出轨迹第10页例6-2-1运行结果在检验曲线正确后,键入hold命令,把曲线保留下来，方便用一样初速，不一样射角，比较其曲线和最大射程。运行结果输入初始速度(m/s):50,输入初速方向(度):40输入目标高度（米）:8得xf=237.4738而初速方向为50度时，xf=241.0454所得曲线见图6-2-1.第11页例6-2-2 质点平面运动 给定质点沿x和y两方向运动规律 x(t)和 y(t),求其运动轨迹,并计算其对原点角动量。解解:建模：由用户输入解析表示式需要用到字符串输入语句,其第二变元为s,而运行这个字符串要用eval命令.当x(t)和 y(t)都是周期运动时,所得曲线就是李萨如图形.动量矩等于动量与向径叉乘(cross product).求速度需要用导数,可用MATLABdiff函数作近似导数计算。设角动量为 ，质点动量为 ，向径为 ，则在x-y平面上投影为第12页程序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;%画出轨迹及角动量随时间改变曲线第13页程序运行结果运行此程序，输入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)画出轨迹及角动量曲线第14页例6-2-3 质点系动力学物体A(质量为m1)在含有斜面物体B(质量为m2)上靠重力下滑,设斜面和地面均物摩擦力,求A沿斜面下滑相对加速度a1和B加速度a2,并求斜面和地面支撑力N1及N2。解：建模，对物体A，列出方程对物体B，列出方程第15页方程组矩阵建模四个方程包含四个未知数，将含未知数项移到等式左边，常数项移到等式右端，得到矩阵方程于是有 X=AB第16页程序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)第17页运行结果输入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矩阵除法就能够方便而准确求出其解.第18页例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动能守恒化为第19页碰撞问题方程由（3）（5）代入（4）（6）主动球能量损失为展开并整理多项式(6)，得可用roots命令求根,第20页程序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%绘图第21页程序运行结果数字结果为(省略了几行)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为正,说明主动球将继续沿原来方向运动.第22页例6-3-1 麦克斯韦速度分布律命题：求摄氏27度下氮气分子运动速度分布率，并求速度在300500m/s范围内分子所占百分比，讨论温度T及分子量mu对速度分布曲线影响。解:建模 麦克斯韦速度分布律为本例将说明怎样从复杂数学公式中绘制曲线并研究单个参数影响.先把麦克斯韦速度分布律列成一个子程序,方便经常调用,把一些惯用常数也放在其中,主程序就简单了.第23页麦克斯韦分布律子程序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;第24页主程序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)%求该范围概率积分第25页执行此程序所得曲线积分结果为：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，使速度分布向高端移动；这是与物理概念相一致。麦克斯韦分布曲线第26页6.4 静电场【例6-4-1】设电荷均匀分布于从z=-L到z=L经过原点线段上，其密度为q库仑/米，求出在xy平面上电位分布。解：建模 点电荷产生电位可表为V=Q/4r0其中r为电荷到测量点距离.线电荷所产生电位可用积分或叠加方法来求。为此把线电荷分为长为dLN段(在MATLAB中,dL应了解为L)。每段上电荷为q*dL.它产生电位为 然后对全部电荷求和即可。第27页程序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%对RN+1点循环计算Rk=sqrt(Lm.2+R(k)2);%测量点到电荷段向径长度Vk=C0*dL*q./Rk;%第k个电荷段产生电位V(k)=sum(Vk);%对各电荷段产生电位求和Endplot(R,V),grid 第28页程序运行结果(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，上图为半对数坐标，下列图为线性坐标。线电荷产生静电场分布第29页6-4-2 由电位表示式计算电场已知空间电位分布，画出等电位线和电场方向.解:建模假如已知空间电位分布 V=V(x,y,z),则空间电场等于电位场负梯度其中 分别为x,y,z三个方向单位向量。MATLAB中设有gradient函数，它是靠数值微分，所以空间观察点应取得密一些，以取得较高精度。第30页程序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);第31页程序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左电位分布曲面,右面是电场分布向量图.第32页Exn642运行结果图6-4-2 V(x,y)=log(x.2+y.2)电位三维立体图，等高线及电场分布图第33页6-5 恒稳磁场恒稳磁场例6-5-1 用毕奥萨伐定律计算电流环产生磁场解:建模载流导线产生磁场基本规律为，任一电流元 在空间任一点P处所产生磁感应强度 为以下向量叉乘积，即其中，为电流元到P点矢径，为导线元长度矢量，P点总磁场可沿载流导体全长积分各段产生磁场来求得。第34页程序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;%各段中点坐标分量第35页程序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第36页程序运行结果%用quiver 画磁场向量图clf;quiver(x,y,Bx,By);图形标注语句及在图上画出圆环位置语句略运行此程序所得图形见图6-5-1,读者可改变电流环直径来分析其影响，也可加上显示各点磁场强度语句来分析其强度分布。图6-5-1 电流环产生磁场分布图第37页例6-5-2 亥姆霍兹线圈一对相同共轴彼此平行载流圆线圈,当它们间距恰好等于其线圈半径时,称之为亥姆霍兹线圈.计算表明,亥姆霍兹线圈轴线附近磁场是十分均匀,而且都沿x方向.本题要求对这一论断进行验证。解:建模 本题计算模型与上例相同,只是把观察区域取在两线圈之间小范围内，如右图所表示。第38页程序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);第39页程序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);第40页程序exn652运行结果第41页6-6 振动与波振动与波例6-6-1 振动合成及拍频现象解:建模，将两个同方向振动相加，可得当 很靠近时，成为一个很低频率，称为拍频。用MATLAB程序得到图形和声音中能够很清楚地看出拍频现象。第42页程序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)第43页程序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 拍频现象第44页例6-6-2 多普勒效应例6-2-2 设声源从500m外以50m/s速度对听者直线开来,其轨迹与听者最小垂直距离为20m.声源角频率为1000弧度/s,试求出听者接收到信号波形方程并生成其对应声音。解:建模设声源发出信号为f(t)，传到听者处，被听者接收信号经历了声音传输拖延，拖延时间为 其中c为音速，为声源与听者之间距离。被接收信号形式为（不考虑声波传输衰减）。只要给出随t改变关系，即可求得并将它恢复为声音信号。第45页程序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);第46页程序exn662运行结果打开计算机声音系统,运行此程序将会听到类似于火车汽笛声音.第一声是火车静止时汽笛声,第二声是本题中听者听到运动火车汽笛声,它频率先高于,后低于原来汽笛声.程序中两个sound语句之间加pause语句是不可少,而且暂停时间要足够长，方便再打开声音系统，这个量与计算机硬件相关，在本书中用5。第47页6.7 光学【例6-7-1】单色光经过两个窄缝射向屏幕，相当于位置不一样两个同频同相光源向屏幕照射叠合，因为抵达屏幕各点距离（光程）不一样引发相位差，如图6-7-1所表示。叠合结果在有点加强，在有点抵消，造成干涉现象。考虑到纯粹单色光不易取得，通常都有一定光谱宽度，它对光 干涉会产生何站种效应，要求用MATLAB计算并仿真这一问题。第48页单色光双缝干涉模型解:建模 考虑两个离中心点距离各为d/2相干光源S1和S2到屏幕上任意点距离差引发相位差,先分析光程则光程差为L=L1-L2将L除以波长，并乘以2，得到相位差 。设两束相干光在屏幕上产生幅度相同,均为A0,则夹角为两个向量A0合成向量幅度为:A=2 A0 cos(/2)光强B正比于振幅平方,故有:B=4 B0 cos2(/2)依据这些关系式,能够编写出计算屏幕上各点光强程序,第49页程序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);%画图象第50页程序exn671运行结果运行exn671和程序所得屏幕光强图像见图6-7-2。光非单色性造成干涉现象减弱。光谱很宽光将不能形成干涉。图6-7-2单色光双缝干涉条纹及光强分布第51页考虑光非单色性再研究复杂一些问题,考虑到光非单色性对干涉条纹影响。此时波长将不是常数，必须对不一样波长光作分类处理再叠加起来。假定光源光谱宽度为中心波长正负10%，而且在该区域内均匀分布。在(0.91.1)之间，按均等间距近似取11根谱线，其波长分别为则上面求相位差计算式求出将是对不一样谱线11个不一样相位。计算光强时应把这11根谱线产生光强迭加取平均值,即第52页程序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;%叠加各波长影响计算光强其它不变，运行结果见右图。不纯光双缝干涉条纹及光强分布能够看出，光非单色性造成干涉现象减弱。光谱很宽光将不能形成干涉。第53页例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也是一个一样长度数组，其计算公式为第54页编程思绪以光源中心到屏幕中心之间距离z为基准，则其光程差等于L-z，对应光相位差为，它应该是ys和yPoint函数。在计算时，先取定屏幕上一点ys，让yPoint取全部值，把得到Phi相加起来，就算出了屏幕上这点光强。然后再换另一个屏幕点ys，再作循环。这么双重循环，计算量是很大，比如把光源数和屏幕点数都取51，则进行循环将为2500次以上。用手工计算是极难想象。下面就是按此思绪写出程序：第55页程序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是一个数组第56页程序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%屏幕上光强与位置关系曲线图形标注语句第57页程序exn672运行结果:依次输入波长Lambda=500nm,距离z=1M，缝宽为a=0.2,1,2mm三种情况，程序运行结果见图6-7-3.三种情况统称费涅耳衍射,只有最左边情况符合夫琅和费衍射条件:*a2/(4*Z)1这个现象也适合用于电磁波发射,天线设计和测量都要用这个概念.天线探测目标时通常符合夫琅和费衍射条件，形整天线远场波瓣。在天线测量时却希望近些才方便，这时往往处于近场，于是要建立二者之间转换关系，MATLAB程序能够发挥作用。第58页程序exn672运行结果(图)图6-7-3 缝宽a为0.2,1,2mm三种情况所得衍射光强曲线（左为夫琅和费衍射）第59页第60页