《MATLAB實踐教程》課件第5章

5.1MATLAB程序結(jié)構(gòu)及其運行機(jī)制

5.2MATLAB程序的流程控制

5.3編寫MATLAB函數(shù)

5.4編程案例

第5單元MATLAB程序設(shè)計

MATLAB程序由主程序和函數(shù)(亦稱子程序)兩部分組成，如圖5-1所示。編程采用模塊化程序設(shè)計方法，用戶任務(wù)分解為多個子任務(wù)，子任務(wù)編寫成用戶自行命名的MATLAB函數(shù)。主程序設(shè)計成若干個模塊，每個模塊實現(xiàn)特定的數(shù)據(jù)處理目標(biāo)，模塊調(diào)用MATLAB函數(shù)完成子任務(wù)，多個模塊順序銜接構(gòu)成主程序，在CommandWindow中鍵入主程序的存盤文件名則可執(zhí)行和完成全部任務(wù)。5.1MATLAB程序結(jié)構(gòu)及其運行機(jī)制圖5-1MATLAB程序的結(jié)構(gòu)與數(shù)據(jù)流5.2.1for循環(huán)

(1)利用循環(huán)向量化技術(shù)設(shè)計for循環(huán)，計算數(shù)列｛135…｝的前15項累加和與累乘積。程序如下：

5.2MATLAB程序的流程控制clc;closeall;clearall;

x=0;y=1; %賦初值

fork=1:2:2*15-1; %創(chuàng)建循環(huán)變量

x=x+k;y=y*k; %計算累加和與累乘積

end;

formatshorte;Results=[xy] %計算結(jié)果

程序執(zhí)行結(jié)果如下：

Results=

2.2500e+0026.1903e+015

(2)利用陣列預(yù)分配技術(shù)設(shè)計for循環(huán)，計算下面離散系統(tǒng)在正弦輸入下的響應(yīng)yk。

程序如下:

clc;closeall;clearall;tic %初始化和啟動秒表計時器

space=0.001;t=0:space:16;u=sin(2*pi*t); %計算正弦輸入

n=length(u); %計算u的數(shù)據(jù)個數(shù)

y=zeros(n,1); %按u的數(shù)據(jù)個數(shù)n、實施陣列y的預(yù)分配

y(1)=2*u(1);y(2)=0.75*y(1)+2*u(2);%計算前2個不能遞推計算的響應(yīng)

fork=3:n; %創(chuàng)建循環(huán)變量

y(k)=0.75*y(k-1)-0.125*y(k-2)+2*u(k); %用遞推公式計算響應(yīng)y

end;

Time=toc %測定程序運行時間

plot(t,y) %響應(yīng)y的計算結(jié)果展示

程序執(zhí)行結(jié)果如下，圖形如圖5-2所示。

Time=

5.1431e-003圖5-2離散系統(tǒng)的時間響應(yīng)y(t)5.2.2while循環(huán)

while循環(huán)程序設(shè)計應(yīng)注意兩點，一是循環(huán)變量累加，二是循環(huán)變量關(guān)系表達(dá)式，兩者構(gòu)成循環(huán)控制邏輯。

用關(guān)系表達(dá)式設(shè)計while循環(huán)，計算函數(shù)在區(qū)間上的面積。程序如下：

clc;closeall;clearall;tic %初始化和啟動秒表計時器

space=pi/30;theta=0:space:pi;y=sin(theta); %生成被積函數(shù)數(shù)據(jù)

x=0;k=1; %給面積變量x和循環(huán)變量k賦初值

whilek<=length(y); %用下標(biāo)k作循環(huán)變量，設(shè)計關(guān)系表達(dá)式

%whilek<length(y)+1; %可選的循環(huán)變量關(guān)系表達(dá)式x=x+y(k)*space; %累計計算面積和

k=k+1; %循環(huán)變量增1累加

end;

time＝toc

disp('k='),disp(k)

disp('aera='),disp(x)

程序執(zhí)行結(jié)果：

time=

0.0065

k=

32

aera=

1.99825.2.3if-else-end結(jié)構(gòu)

if-else-end結(jié)構(gòu)關(guān)系表達(dá)式根據(jù)檢驗結(jié)果有條件地執(zhí)行程序。

(1)編寫一個MATLAB函數(shù)用于判斷一個小于50的數(shù)是否為素數(shù)，并存盤為prime.m。程序如下：

functiony=prime(x)

ifrem(x,2)==0&x>2;

y='不是素數(shù),';

elseifrem(x,3)==0&x>3;

y='不是素數(shù),';

elseifrem(x,5)==0&x>5;

y='不是素數(shù),';

elseifrem(x,7)==0&x>7;

y='不是素數(shù),';

else;

y='是素數(shù)，';

end;編寫主程序指定要判斷的多個數(shù)字(向量x)，存盤名為myprime01.m。程序如下：

clc;closeall;clearall;

x=[123456789101112131415]; %指定要判別的數(shù)字

n=length(x);p=[]; %判別結(jié)果變量p賦初值

fori=1:n; %設(shè)計循環(huán)逐個判別

p=[p[num2str(x(i))prime(x(i))]]; %判別結(jié)果放到一個字符串變量p里

end;

disp(p);%顯示判別結(jié)果

向量素數(shù)判斷主程序的執(zhí)行結(jié)果如下：

1是素數(shù),2是素數(shù),3是素數(shù),4不是素數(shù),5是素數(shù),6不是素數(shù),7是素數(shù),8不是素數(shù),9不是素數(shù),10不是素數(shù),11是素數(shù),12不是素數(shù),13是素數(shù),14不是素數(shù),15不是素數(shù),

(2)將判斷一個矩陣各個元素是否為素數(shù)的任務(wù)編寫成一個MATLAB函數(shù)，并存盤為primematrix.m。程序如下：

functiony=primematrix(x)

x1=rem(x,2)==0&x>2;

x2=rem(x,3)==0&x>3;

x3=rem(x,5)==0&x>5;

x4=rem(x,7)==0&x>7;

y=find((x1+x2+x3+x4)==0);

編寫主程序指定要判斷的多個數(shù)字(矩陣x)，存盤名為myprime02.m。程序如下：

clc;closeall;clearall;

x=[12345;678910;1112131415]%指定要判別的矩陣

p=[primematrix(x)]'%判別結(jié)果為向量p，其元素值為矩陣x中素數(shù)元素的序號矩陣素數(shù)判斷主程序的執(zhí)行結(jié)果如下：

x=

12345

678910

1112131415

p=

134579135.2.4switch-case-end結(jié)構(gòu)

switch-case-end結(jié)構(gòu)根據(jù)狀態(tài)表達(dá)式的值執(zhí)行不同語句塊，相當(dāng)于多條if塊的嵌套使用。

將判斷一個數(shù)是否為素數(shù)的任務(wù)編寫成一個MATLAB函數(shù)并存盤為oneprime.m。程序如下：

functiony=oneprime(x)

u=(rem(x,2)==0&x>2)|(rem(x,3)==0&x>3)|(rem(x,5)==0&x>5)…

|(rem(x,7)==0&x>7);

switchu;

case0;

y='是素數(shù),';

case1;

y='不是素數(shù),';

otherwise;%這部分語句塊在本函數(shù)可省略

y='是素數(shù),';

end;編寫主程序指定要判斷的多個數(shù)字(向量x)，存盤名為myprime03.m。程序如下：

clc;closeall;clearall;

x=[123456789101112131415];%指定要判別的數(shù)字

n=length(x);p=[];%判別結(jié)果變量p賦初值

fori=1:n;%設(shè)計循環(huán)逐個判別

p=[p[num2str(x(i))oneprime(x(i))]];%判別結(jié)果放到一個字符串變量p里

end;

disp(p);%顯示判別結(jié)果

主程序執(zhí)行結(jié)果如下：

1是素數(shù),2是素數(shù),3是素數(shù),4不是素數(shù),5是素數(shù),6不是素數(shù),7是素數(shù),8不是素數(shù),9不是素數(shù),10不是素數(shù),11是素數(shù),12不是素數(shù),13是素數(shù),14不是素數(shù),15不是素數(shù),5.3.1MATLAB函數(shù)的結(jié)構(gòu)及其調(diào)用機(jī)制

MATLAB函數(shù)是一種特定結(jié)構(gòu)的程序，在MATLAB環(huán)境里的程序與一般程序一樣都是可視的文本，因其存盤名的擴(kuò)展名為“.m”，亦稱作M-file函數(shù)。

與一般程序不同，M-file函數(shù)不能在CommandWindow中通過簡單鍵入它的存盤基本名執(zhí)行，而是采用按M-file函數(shù)定義行所指定的函數(shù)格式鍵入函數(shù)表達(dá)式的方法調(diào)用，主程序調(diào)用函數(shù)亦應(yīng)遵循指定的函數(shù)格式。M-file函數(shù)的標(biāo)準(zhǔn)結(jié)構(gòu)如圖5-3所示。5.3編寫MATLAB函數(shù)圖5-3M-file函數(shù)的結(jié)構(gòu)在Editor窗口編寫M-file函數(shù)，存盤名由“基本名+擴(kuò)展名”構(gòu)成，存盤時僅僅鍵入基本名，擴(kuò)展名默認(rèn)為“.m”，函數(shù)調(diào)用時使用存盤的基本名和指定的函數(shù)格式。注意，MATLAB7.1版本以后存盤函數(shù)基本名與函數(shù)定義行指定的函數(shù)名可以不一致，而早先的版本要求存盤基本名與函數(shù)定義行的指定函數(shù)名必須一致。下面的函數(shù)示例計算(y,outp)=f(x,inp)，其中輸入?yún)?shù)x為自變量，輸入?yún)?shù)inp為指定擬合的多項式階數(shù)，輸出參數(shù)y為因變量，且

，輸出參數(shù)outp為依據(jù)(x,y)擬合的inp階多項式。程序如下：

function[y,outp]=myfun(x,inp);

y=100*sin(x).*exp(-x.^2);

outp=polyfit(x,y,inp);命名所編寫的M-file函數(shù)并存盤。(如存盤為myfun.m)

在CommandWindow窗口鍵入函數(shù)表達(dá)式，且預(yù)先為輸入?yún)?shù)賦值；在主程序里編寫函數(shù)表達(dá)式調(diào)用函數(shù)，且預(yù)先為輸入?yún)?shù)賦值。命令如下：

x=1:7;order=3;

[y,p]=myfun(x,order)

執(zhí)行函數(shù)調(diào)用的結(jié)果如下：

y=

30.95601.66540.0017-0.0000-0.0000-0.00000.0000

p=

-0.813611.6055-51.532969.80295.3.2編寫匿名函數(shù)

匿名函數(shù)指在主程序中以表達(dá)式方式定義的函數(shù)，下面示例展示它的格式和調(diào)用方法，程序如下：

quade=@(x,a,b,c)a*x.^2+b*x+c; %定義匿名函數(shù)并賦給變量quade

x=[123];a=2;b=-1;c=5; %輸入?yún)?shù)賦值

y=quade(x,a,b,c)%以調(diào)用變量quade調(diào)用匿名函數(shù)，注意輸入?yún)?shù)的位置和意義

主程序執(zhí)行結(jié)果如下：

y=

611205.3.3編寫標(biāo)準(zhǔn)M-file函數(shù)

標(biāo)準(zhǔn)M-file函數(shù)，指函數(shù)體中的語句不調(diào)用其它任何函數(shù)，均采用直接編寫的算式。

下面函數(shù)在x的多個點上計算二次函數(shù)的值，存盤名為stdfun.m。程序如下：

functiony=stdfun(x,a,b,c)

y=a*x.^2+b*x+c;

編寫調(diào)用函數(shù)stdfun的主程序。程序如下：

clc;closeall;clearall;%初始化

x=[123];a=2;b=-1;c=5;

y=stdfun(x,a,b,c)

主程序調(diào)用函數(shù)的執(zhí)行結(jié)果如下：

y=

611205.3.4編寫私人函數(shù)

私人函數(shù)是一種特殊的M-file函數(shù)，特別放置在名為private的文件夾里，僅僅能被父文件夾(上一級文件夾)里的函數(shù)調(diào)用。父文件夾之外的函數(shù)無法調(diào)用私人函數(shù)，故私人函數(shù)可以與其他文件夾的函數(shù)名相同，且一個函數(shù)調(diào)用M-file函數(shù)時優(yōu)先搜尋私人函數(shù)。若用戶想限制某個M-file函數(shù)的使用權(quán)限或隱蔽調(diào)用，可將其放置在private文件夾里。

假定用戶文件夾mydir在MATLAB的搜索路徑上，mydir的一個子文件夾命名為private，則private里的函數(shù)(稱作私人函數(shù))僅僅能被mydir里的函數(shù)調(diào)用。注意，是父文件夾里的函數(shù)調(diào)用私人函數(shù)，而不是父文件夾里的主程序調(diào)用。創(chuàng)建用戶自己的文件夾mydir，編寫下面函數(shù)myprog.m并存放在mydir文件夾下，該函數(shù)設(shè)計為調(diào)用私人函數(shù)sin，操作File菜單的setpath項將mydir文件夾添加到MATLAB搜索路徑。

functiony=myprog(x)

y=sin(x)

在mydir文件夾內(nèi)創(chuàng)建private文件夾，編寫私人函數(shù)sin.m并存放在private文件夾里。程序如下：

functiony=sin(x);

y=x.^2;%故意編寫一個錯誤表達(dá)式編寫下面主程序mytest.m并可存放在任何MATLAB搜索路徑的文件夾里。程序如下：

x=[123;456]

y_common=sin(x) %此表達(dá)式將調(diào)用MATLAB系統(tǒng)的函數(shù)sin.m

y_private=myprog(x) %此表達(dá)式將調(diào)用mydir文件夾下函數(shù)myprog.m

%myprog.m函數(shù)將調(diào)用private文件夾下私人函數(shù)sin.m

試比較如下主程序的執(zhí)行結(jié)果：

x=

123

456

y_common=

0.8415

0.90930.1411

-0.7568-0.9589-0.2794

y_private=

149

1625365.3.5編寫形如“主函數(shù)+子函數(shù)”的M-file函數(shù)

如果一個M-file函數(shù)的函數(shù)體內(nèi)調(diào)用一個或多個其它M-file函數(shù)，則稱其為主函數(shù)，被調(diào)用的M-file函數(shù)稱作子函數(shù)。在用戶任務(wù)里，主函數(shù)被主程序或命令行調(diào)用，而子函數(shù)僅僅被主函數(shù)或其它子函數(shù)調(diào)用。

下面的主函數(shù)rank調(diào)用svd、nargin、max、size、sum等5個子函數(shù)，以允差tol計算矩陣A的秩，主函數(shù)rank的存盤名為rank.m。程序如下：functionr=rank(A,tol)%RANKMatrixrank.s=svd(A);ifnargin==1;tol=max(size(A)')*max(s)*eps;end;r=sum(s>tol);編寫主程序計算矩陣x的秩。程序如下：

x=[1524;3679;5293];

rx=rank(x,1e-005)

主程序調(diào)用主函數(shù)rank的執(zhí)行結(jié)果如下：

rx=

35.3.6編寫形如“主函數(shù)+嵌套函數(shù)”的M-file函數(shù)

“主函數(shù)+嵌套函數(shù)”是一種特殊結(jié)構(gòu)的M-file函數(shù)，主函數(shù)供主程序或其它函數(shù)調(diào)用，嵌套函數(shù)按包含關(guān)系逐層嵌套在主函數(shù)內(nèi)。一個嵌套函數(shù)可調(diào)用任何外部的M-file函數(shù)，但它自身只能在主函數(shù)內(nèi)部被調(diào)用，且遵循下述規(guī)則：(1)主函數(shù)或父嵌套函數(shù)可調(diào)用子嵌套函數(shù)；(2)主函數(shù)或嵌套函數(shù)不能越層調(diào)用深層嵌套函數(shù)；(3)底層嵌套函數(shù)可調(diào)用高層嵌套函數(shù)，包括越層調(diào)用；(4)一個嵌套函數(shù)可調(diào)用同一層其它嵌套函數(shù)。

編寫下面嵌套函數(shù)，存盤為mynested.m。

function[y,p]=mynested(x,a,b,c,order)

y=computing(x,a,b,c);

[p,x1,y1]=polyfiting(x,y,order);

ploting(x,y,x1,y1);

functiony_value=computing(x,a,b,c)

y_value=a*x.^2+b*x+c;

end

function[p_poly,x1,y1]=polyfiting(x,y,order)

p_poly=polyfit(x,y,order);

[x1,y1]=interposing(x,y);

function[x1,y1]=interposing(x,y)

x1=min(x):0.01:max(x);

y1=interp1(x,y,x1);

end

end

functionploting(x,y,x1,y1)

plot(x1,y1,'--',x,y,'o');

end

end

編寫下面主程序調(diào)用嵌套函數(shù)。

clc;closeall;clearall;

x=linspace(-1.5,3,21);a=2;b=-1;c=5;order=3;

[y,p]=mynested(x,a,b,c,order)主程序調(diào)用嵌套函數(shù)結(jié)果之一如下：

y=

Columns1through11

11.00009.52628.25507.18626.32005.65635.19504.93634.88005.02635.3750

Columns12through21

5.92636.68007.63628.795010.156311.720013.486315.455017.626320.0000

p=

0.00002.0000-1.00005.0000

主函數(shù)調(diào)用嵌套函數(shù)結(jié)果之二如圖5-4所示。圖5-4主程序調(diào)用嵌套函數(shù)的結(jié)果之二5.3.7編寫函數(shù)的函數(shù)

函數(shù)的函數(shù)(functionfunctions)是一種供求根(Zerofinding)、優(yōu)化(Optimization)、積分(Quadrature)和解常微分方程(Ordinarydifferentialequations)等工具函數(shù)調(diào)用的函數(shù)。詳見第4單元積分和解常微分方程的函數(shù)編程，亦可參見第8單元優(yōu)化設(shè)計的函數(shù)編程。5.3.8全局變量和局部變量

在M-file函數(shù)中定義一些變量x1,x2,…為全局變量，使用語句globalx1x2…。如果M-file函數(shù)中的變量沒有定義為全局變量(GlobalVariable)，則默認(rèn)全是局部變量(LocalVariable)。外部程序或外部M-file函數(shù)不能直接調(diào)用一個函數(shù)內(nèi)的局部變量，故變量名可以與外部變量同名，但卻可直接調(diào)用全局變量。編程采用全局變量，遵循先定義后使用的規(guī)則。盡量避免采用全局變量。

編寫含全局變量a、b、c的M-file函數(shù)quading.m。程序如下：

functiony_value=quading(x)

globalabc;

y_value=a*x.^2+b*x+c;

end

編寫含同名全局變量a、b、c的主程序，實現(xiàn)主程序與函數(shù)之間的參數(shù)傳遞，程序如下：

clc;closeall;clearall;

globalabc;

x=linspace(-1.5,3,21);a=2;b=-1;c=5;order=3;

y=quading(x)主程序調(diào)用函數(shù)的執(zhí)行結(jié)果如下：

y=

Columns1through11

11.00009.52628.25507.18626.32005.65635.19504.93634.88005.02635.3750

Columns12through21

5.92636.68007.63628.795010.156311.720013.486315.455017.626320.00005.3.9函數(shù)句柄

無論是系統(tǒng)的M-file函數(shù)或是自定義的M-file函數(shù)，都可以用函數(shù)句柄(FunctionHandle)來代表，也可用一個自定義的變量名替代，程序調(diào)用函數(shù)被調(diào)用函數(shù)句柄代替，這樣可使程序變得簡潔和更容易理解。

創(chuàng)建函數(shù)句柄并運用它。程序如下：

clc;clearall;closeall;

sqr=@(x)x.^2; %為匿名函數(shù)f(x)=x2創(chuàng)建一個函數(shù)句柄sqr

x=0:pi/3:2*pi;

y=sqr(5) %調(diào)用函數(shù)句柄sqr計算函數(shù)f(x)=x2的值

sqr=@sin; %為系統(tǒng)的M-file函數(shù)sin創(chuàng)建一個函數(shù)句柄sqr

a=sqr(x) %調(diào)用函數(shù)句柄sqr計算函數(shù)f(x)=sin(x)的值

b=quad(sqr,0,pi) %調(diào)用函數(shù)句柄sqr計算函數(shù)f(x)=sin(x)在區(qū)間[0,π]上的積分上面程序的執(zhí)行結(jié)果如下：

y=

25

a=

00.86600.86600.0000-0.8660-0.8660-0.0000

b=

2.0000利用函數(shù)句柄技術(shù)編寫一個“主函數(shù)+嵌套函數(shù)”結(jié)構(gòu)的M-file函數(shù)，主函數(shù)命名并存盤為get_plot_handle.m，主函數(shù)的輸入?yún)?shù)為plot圖所需的線型設(shè)置字符串變量LS、線寬數(shù)值變量LW、點標(biāo)記邊界顏色字符串變量MEC、點標(biāo)記表面顏色字符串變量MFC和點標(biāo)記尺寸數(shù)值變量MS。嵌套函數(shù)命名為draw_plot，它的輸入?yún)?shù)為x,y，函數(shù)體中的plot函數(shù)調(diào)用主函數(shù)變量LS、LW、MEC、MFC和MS。程序如下：

functionh=get_plot_handle(LS,LW,MEC,MFC,MS)

h=@draw_plot;

functiondraw_plot(x,y)

plot(x,y,LS,'LineWidth',LW,...

'MarkerEdgeColor',MEC,...

'MarkerFaceColor',MFC,...

'MarkerSize',MS)

end

end采用函數(shù)句柄技術(shù)編寫主程序調(diào)用上面的函數(shù)。程序如下：

clc;clearall;closeall;

h=get_plot_handle('--rs',2,'k','g',10);%創(chuàng)建函數(shù)句柄h加載圖形設(shè)置數(shù)據(jù)

x=-pi:pi/10:pi; %計算繪圖數(shù)據(jù)x

y=tan(sin(x))-sin(tan(x)); %計算繪圖數(shù)據(jù)y

h(x,y) %調(diào)用函數(shù)句柄h，它按加載的圖形設(shè)置數(shù)據(jù)繪圖

主程序執(zhí)行的結(jié)果如圖5-5所示。圖5-5采用函數(shù)句柄技術(shù)編程的繪圖結(jié)果5.4.1分析圖像的灰度分布

編寫顯示原圖的M-file函數(shù)original_image.m。程序如下：

functionoriginal_image(image)

%ShowaOriginalImage

imshow(image)5.4編程案例編寫顯示圖像灰度直方圖的M-file函數(shù)gray_histogram.m。程序如下：

functiongray_histogram(image,N);

%ShowtheGrayscaleHistogramofaImage

ifnargin==1;

imhist(image);boxoff;

else;

imhist(image,N);boxoff;

end;

編寫主程序myimage_analysis.m，實現(xiàn)原圖和灰度頻數(shù)直方圖的顯示。程序如下：

clc;clearall;closeall;

char1='D:\users\pictures\';char2='sample01.jpg';

h1=@original_image;h2=@gray_histogram;

xx=imread([char1char2]);yy=rgb2gray(xx);

subplot(2,2,1),h1(xx);title('OriginalImage')

subplot(2,2,2),h1(yy);title('GrayImage')

subplot(2,2,3),h2(yy);title('GrayscaleHistogram')

subplot(2,2,4),h2(yy,99);title('GrayscaleHistogram')

主程序的執(zhí)行結(jié)果如圖5-6所示。圖5-6案例主程序的執(zhí)行結(jié)果5.4.2懸掛機(jī)組機(jī)具跟蹤性分析

輪式拖拉機(jī)與鏵式犁懸掛連接組成懸掛犁耕機(jī)組，當(dāng)拖拉機(jī)轉(zhuǎn)彎時，機(jī)具末端外側(cè)偏離理想軌跡的程度稱作機(jī)具跟蹤性。機(jī)組轉(zhuǎn)彎參數(shù)滿足下面的懸掛跟蹤方程：

試求機(jī)具偏角與拖拉機(jī)前輪轉(zhuǎn)角的函數(shù)關(guān)系=f?(θ)(相關(guān)參數(shù)如圖5-7所示)。圖5-7懸掛犁耕機(jī)組和轉(zhuǎn)彎參數(shù)

解:懸掛跟蹤方程是一個隱函數(shù)方程，不能直接獲得的解析表達(dá)式。因此，擬采用數(shù)值解法求出?

與??的n對觀察值，再用多項式擬合求出的近似函數(shù)關(guān)系。

算法流程設(shè)計如圖5-8所示。圖5-8懸掛犁耕機(jī)組機(jī)具跟蹤性分析程序的算法流程編寫懸掛跟蹤方程求根的M-file函數(shù)myequation.m，即給定θ解出。程序如下：

function[phi,fval]=myequation(theta)

L=16*0.3048;h=10*0.3048;b=4*0.3048;

n=length(theta);tol=1e-5;

m=floor(pi/(2*tol));

phi0=linspace(0,pi/2,m);

phi=zeros(n,1);index=zeros(n,1);fval=zeros(n,1);

fori=1:n;

f0=((L+h)*cos(phi0)-b*sin(phi0))*sin(theta(i))-L*sin(theta(i)+phi0);

[fval(i),index]=min(abs(f0));

ph