双向解析光束法

双向解析光束法
光束法程序有问题，在Getelement这个函数里便出现索引超限，这个问题一直解决不了
光束法的流程：
1.根据同名像点对对相交理论求系数阵A，系数阵B和常数阵L
a11=(a1f+a3x)/Z; a12=(b1f+b3x)/Z; a13=(c1f+c3x)/Z; a14=ysin(omega)-[x/f(xcos(kappa)-ysin(kappa))+fcos(kappa)]cos(omega);
a25=-fcos(kappa)-y/f(xsin(kappa)+ycos(kappa));
a26=-x;
2.求方程的改化法方程求出外方位元素和物方坐标改正数
3.判断改正数的值，如果小于限差则输出结果
光束法是最严密的一种方法的原因：
在一张相片中，待定点与控制点的像点与摄影中心及相应地面点均构成一条光束，该方法是以每张相片所组成的一束光线作为平差的基本单元，已共线条件方程作为平差的基础方程，通过各个光束在空间中的旋转和平移，使模型之间公共点的光线实现最佳交汇，并使整个区域纳入到已知的地面控制点坐标系中，所以要建立全区域统一的误差方程，整体解求全区域内每张相片的六个外方位元素及所有待定点坐标，光束法区域网平差是基于摄影时像点，物点和摄站点三点共线提出来的。由单张相片构成区域，其平差的数学模型是共线条件方程，平差单元是单个光束，像点坐标是观测值，未知数是每张相片的外方位元素及所有待定点坐标。误差方程直接由像点坐标的观测值列出，能对像点坐标进行系统误差改正。
}
//求右片的Xs,Ys,Zs
Xsr = (strX[1] + strX[3]) / 2;
Ysr = (strY[1] + strY[3]) / 2;
L = Math.Pow(Math.Pow(strX[1] - strX[3], 2) + Math.Pow(strY[1] - strY[3], 2) + Math.Pow(strZ[1] - strZ[3], 2), 0.5);
//phi,omega,kappa为零，故旋转矩阵为单位阵
//像空间辅助坐标系中的坐标
Class1 uvw = XYZt.Multiply(H);
//求地面摄影测量坐标系中的坐标
for (int i = 0; i < 8; i++)
{
uvw.SetElement(0, i, (uvw.GetElement(0, i) + Xsr+UVW.GetElement (0,i))/2);
l = Math.Pow(Math.Pow(strXl[0] - strXl[2], 2) + Math.Pow(strYl[0] - strYl[2], 2), 0.5);
H = f * L / l;
Zsl = (strZ[0] + strZ[2]) / 2 + H;
//计算片左的加密点的物方坐标
strxy = new double[24];
Console.WriteLine("请输入右片加密点的相片坐标");
for (int i = 0; i < 8; i++)
{
//strxy[0 + 3 * i] = Convert.ToDouble(Console.ReadLine());
//strxy[1 + 3 * i] = Convert.ToDouble(Console.ReadLine());
for (int i = 0; i < 8; i++)
{
UVW.SetElement(0, i, UVW.GetElement(0, i) + Xsl);
UVW.SetElement(1, i, UVW.GetElement(1, i) + Ysl);
UVW.SetElement(2, i, UVW.GetElement(2, i) + Zsl);
Class1 xyz = new Class1(8, 3, strXY);
Class1 xyzt=xyz.Transpose ();
//phi,omega,kappa为零，故旋转矩阵为单位阵
//像空间辅助坐标系中的坐标
Class1 UVW = xyzt.Multiply(H);
//求地面摄影测量坐标系中的坐标
l = Math.Pow(Math.Pow(strXr[1] - strXr[3], 2) + Math.Pow(strYr[1] - strYr[3], 2), 0.5);
H = f * L / l;
Zsr = (strZ[1] + strZ[3]) / 2 + H;
//右片加密点
double[] strxy;
uvw.SetElement(1, i, (uvw.GetElement(1, i) + Ysr+UVW.GetElement (1,i))/2);
光束法的程序代码为：
//计算像片外方位元素，架设phi=0,omega=0,kappa=0求Xs,Ys,Zs
//求左片的Xs,Ys,Zs
Xsl = (strX[0] + strX[2]) / 2;
Ysl = (strY[0] + strY[2]) / 2;
L = Math.Pow(Math.Pow(strX[0] - strX[2], 2) + Math.Pow(strY[0] - strY[2], 2) + Math.Pow(strZ[0] - strZ[2], 2), 0.5);
a15=-fsin(kappa)-x/f(xsin(kappa)+ycos(kappa));
a16=y;
a21=(a2f+a3y)/Z; a22=(b2f+b3y)/Z; a23=(c2f+c3y); a24=-xsin(omega)-[x/f(xcos(kappa)-ysin(kappa))-fsin(kappa)]cos(omega)
strxy[2 + 3 * i] = -f;
}
for (int i = 0; i < 24; i++)
{
strxy[i] = strxy[i] * 0.000025;
}
//计算片左的加密点的物方坐标
Class1 XYห้องสมุดไป่ตู้ = new Class1(8, 3, strxy);
Class1 XYZt=XYZ .Transpose ();