下面一段matlab程序,不同参数下面计算的x,y的结果都有两个,应该每一组参数下面都有一对x,y才对,求指导程序和运行结果如下所示：程序：function xnxp_1%%输入基本参数%两张相片内方位角元素q=[0

下面一段matlab程序,不同参数下面计算的x,y的结果都有两个,应该每一组参数下面都有一对x,y才对,求指导
程序和运行结果如下所示：
程序：
function xnxp_1
%%输入基本参数
%两张相片内方位角元素
q=[0.001;0.002];%q=[q1;q2]
w=[0.001;0.002];
k=[0.001;0.002];
%输入相片的内方位元素
x0=[0;0];%x0=[x01;x02]
y0=[0;0];
f=[100;200];
m=[1000;1000];
%两张照片的外方位线元素
Xs=[1000;2000];%含有Xs1,Xs2,两张照片的外方位线元素
Ys=[1000;2000];
Zs=[1000;2000];
%输入世界坐标系统下的三个点
XA=10;%里面有三个点ABC的X坐标
YA=20;
ZA=30;
%%利用共线条件方程生成虚拟相片
for j=1:2
%先输入参数信息
fprintf('Xs = %.2fm\n',Xs(j));
fprintf('Ys = %.2fm\n',Ys(j));
fprintf('Zs = %.2fm\n',Zs(j));
fprintf('q= %11.8f rad\n',q(j));
fprintf('w= %11.8f rad\n',w(j));
fprintf('k= %11.8f rad\n',k(j));
fprintf('f = %.2f\n',f(j));
fprintf('m = %.2f\n',m(j));
fprintf('x0 = %.2f\n',x0(j));
fprintf('y0 = %.2f\n',y0(j));
%求解转换矩阵R的系数
a1=cos(q(j))*cos(k(j))-sin(q(j))*sin(w(j))*sin(k(j));
a2=-cos(q(j))*sin(k(j))-sin(q(j))*sin(w(j))*cos(k(j));
a3=-sin(q(j))*cos(w(j));
b1=cos(w(j))*sin(k(j));
b2=cos(w(j))*cos(k(j));
b3=-sin(w(j));
c1=sin(q(j))*cos(k(j))+cos(q(j))*sin(w(j))*sin(k(j));
c2=-sin(q(j))*sin(k(j))+cos(q(j))*sin(w(j))*cos(k(j));
c3=cos(q(j))*cos(w(j));
%利用共线方程求解
x=(-f*(a1*(XA-Xs(j))+b1*(YA-Ys(j))+c1*(ZA-Zs(j)))/(a3*(XA-Xs(j))+b3*(YA-Ys(j))+c3*(ZA-Zs(j))));
y=(-f*(a2*(XA-Xs(j))+b2*(YA-Ys(j))+c2*(ZA-Zs(j)))/(a3*(XA-Xs(j))+b3*(YA-Ys(j))+c3*(ZA-Zs(j))));
fprintf('x = %.2fmm\n',x);
fprintf('y = %.2fmm\n',y);
%A=[x;y]
%fid=fopen('test.txt','w');
%fprint(fid,'%f\n %f\n',A);%%%%?
%fclose(fid);
%type test.txt
end
end
结果：
Xs = 1000.00m
Ys = 1000.00m
Zs = 1000.00m
q= 0.00100000 rad
w= 0.00100000 rad
k= 0.00100000 rad
f = 100.00
m = 1000.00
x0 = 0.00
y0 = 0.00
x = -102.47mm
x = -204.94mm
y = -101.23mm
y = -202.47mm
Xs = 2000.00m
Ys = 2000.00m
Zs = 2000.00m
q= 0.00200000 rad
w= 0.00200000 rad
k= 0.00200000 rad
f = 200.00
m = 1000.00
x0 = 0.00
y0 = 0.00
x = -101.83mm
x = -203.65mm
y = -100.91mm
y = -201.82mm

clc
alpha=5.2182;
n=1.55;
e=0;
f=0;
for l=1:50
syms x f1 f2
if l==1
f1=sin(x)./x-cos(x);
f1dot=diff(f1,x,l);
f2=(sin(x)+i*cos(x))./x-cos(x)+i*sin(x);
f2dot=diff(f2,x,l);
else
f1=(pi*x/2).*besselj(l+0.5,x);
f1dot=(pi*x/2).*besselj(l-0.5,x)-l*(pi*x/2).*besselj(l+0.5,x)./x;
f2=(pi*x/2).*(besselj(l+0.5,x)-i*bessely(l+0.5,x));
f2dot=(pi*x/2).*(besselj(l-0.5,x)-i*bessely(l-0.5,x))-l*(pi*x/2).*(besselj(l+0.5,x)-i*bessely(l+0.5,x))./x;
end
g1=inline(vectorize(f1),x);
g1dot=inline(vectorize(f1dot),x);
g2=inline(vectorize(f2),x);
g2dot=inline(vectorize(f2dot),x);
a1=(g1dot(n*alpha).*g1(alpha)-n*g1(n*alpha).*g1dot(alpha))./(g1dot(n*alpha).*g2(alpha)-n*g1(n*alpha).*g2dot(alpha));
b1=(n*g1dot(n*alpha).*g1(alpha)-g1(n*alpha).*g1dot(alpha))./(n*g1dot(n*alpha).*g2(alpha)-g1(n*alpha).*g2dot(alpha));
e=e+(2*l+1)*(abs(a1).^2+abs(b1).^2);
f=f+(2*l+1)*real(a1+b1);
end
kext=2/alpha.^2.*f
ksca=2/alpha.^2.*e