摄影测量后方交会求解外方位元素程序.docVIP

摄影测量后方交会求解外方位元素程序 2010-04-11 15:14 ground.txt内容：36589.41 25273.32 2195.1737631.08 31324.51 728.6939100.97 24934.98 2386.5040426.54 30319.81 757.31image.txt内容：-86.15 -68.99-53.40 82.21-14.78 -76.6310.46 64.43程序如下：#include stdio.h#include math.h#define PI 3.1415926//求矩阵a的转置矩阵b，a为m行、n列void transpose(double *a, double *b, int m, int n);//矩阵a乘以矩阵b，结果存储在c中，a为m×n大小，b为n×l大小void multiply(double *a, double *b, double *c, int m, int n, int l);//求矩阵a的逆int inMerse1(double *a, int n);//输出m行、n列的矩阵avoid shuchu(double *a, int m, int n);//求矩阵a与b的差值矩阵，a、b、c大小均为m×nvoid subtrat(double *a, double *b, double *c, int m, int n);void main(){FILE *fp = NULL;FILE *fp1 = NULL;if((fp=fopen(image.txt,r)) == NULL){printf(Open file error!);return;}if((fp1=fopen(ground.txt,r)) == NULL){printf(Open file error!);return;}//像点坐标和地面点坐标double imagecontrol[4][2]={0.0};double groundcontrol[4][3]={0.0};//摄影比例尺分母double scale = 40000.0;double f=0.15324;long i,j,k;for(i=0; i4; i++){for(j=0; j2; j++){fscanf(fp, %lf, imagecontrol[i][j]);imagecontrol[i][j] /= 1000.0;}for(k=0; k3; k++){fscanf(fp1, %lf, groundcontrol[i][k]);}}fclose(fp);fclose(fp1);double Zs(0.0),Xs(0.0),Ys(0.0),p0(0.0),w0(0.0),k0(0.0);Zs=scale*f;for( i=0;i4;i++){Xs+=groundcontrol[i][0];Ys+=groundcontrol[i][1];}Xs/=4.0;Ys/=4.0;double R[3][3]={0.0};double L3=0.0,L1=0.0,L2=0.0;double L[8][1]={0.0},ATA1[6][6]={0.0},x=0.0,y=0.0;double A[8][6]={0.0},AT[6][8]={0.0},ATA[6][6]={0.0},ATL[6][1]={0.0},M[6][1]={0.0},B[6][8]={0.0};double AM[8][1]={0.0},V[6][1]={0.0};int n(0);do{R[0][0]=cos(p0)*cos(k0)-sin(p0)*sin(w0)*sin(k0);R[0][1]=(-1)*cos(p0)*sin(k0)-sin(p0)*sin(w0)*cos(k0);R[0][2]=(-1)*sin(p0)*cos(w0);R[1][0]=cos(w0)*sin(k0);R[1][1]=cos(w0)*cos(k0);R[1][2]=(-1)*sin(w0);R[2][0]=sin(p0)*cos(k0)+cos(p0)*sin(w0)*sin(k0);R[2][1]=(-1)*sin(p0)*sin(k0)+cos(p0)*sin(w0)*cos(k0);R[2][2]=cos(p0)*cos(w0);for(i=0,j=0;j4;i+=2,j++){L1=R[0][0]*(groundcontrol[j][0]-Xs)+