测量平差程序设计1.角度（度分秒）到弧度AngleToRadian#define PI 3.14159265double AngleToRadian（double angle）{int D，M；double S,radian,degree, angle,MS；D=int（angle+0.3）；MS=angle-D；M=int（（MS）*100+0.3）；S=（MS*100-M）*100;degree=D+M/60.0+S/3600.0;radian=degree*PI/180.0;return radian;}注意:防止数据溢出,要加个微小量,例如0.3.2.弧度换角度(度分秒) RadianToAngle#define PI 3.14159265double RadianToAngle(double radian){int D,M;double S,radian,degree,MS,angle;degree=radian*180/PI;D=int(degree);MS=degree-D;M=int(MS*60);S=(MS*60-M)*60;angle=D+M/100.0+S/10000.0;return angle;}3.已知两点求坐标方位角Azimuth#include <math.h>double Azimuth(double xi,double yi,double xj,double yj) {double Dx,Dy,S,T;Dx=xj-xi;Dy=yj-yi;S=sqrt(Dx*Dx+Dy*Dy);if(S<1e-10) return 0;T=asin(Dy/S);if(Dx<0) T=PI-T;if(Dx>0&&(Dy<0)||T<0) T=2*PI+T;return T;}4.开辟二维数组的动态空间的宏#include <malloc.h>#define NewArray2D(type,A,i,n,m){A=(type**)malloc(n*sizeof(type*));\for(i=0;i<m;i++)\A[i]=(type*)malloc(m*sizeof(type));\}5.释放开辟的二维数组的空间#define FreeSpace(A,i,m){for(i=0;i<m;i++)\free(A[ i]);\free(A);\}注意:释放空间与开辟空间相反,释放空间是先释放列,后释放行.6.矩阵求转置transformmatrixvoid transformmatrix（double **A，double **B，int i,int j）{int m,n;for(m=0;m<=i;m++)for(n=0;n<=j;n++){B[n][m]=A[m][n]:}}7.矩阵相乘(mulmatrix)void mulmatrix(double **A,double **B,double **C,int i,int j,int k) {int m,n,p;for(m=0;m<i;m++)for(n=0;n<j;n++){C[m][n]=0;for(p=0;p<k;p++){C[m][n]+=A[m][p]*B[p][n]:}}}8.矩阵求逆(countermatrix)#include <math.h>void countermatrix(double **T, double **s, double **r, double **Q,double **N, double **rt,int n){for(i=0;i<n;i++){s=N[i][i];for(k=0;k<i;k++){s-=T[k][i]*T[k][i];}T[i][i]=sqrt(s)for(j=i+1;j<n;j++){s=N[i][j];for(k=0;k<i;k++){s-=T[k][i]*T[k][j];}T[i][j]=s/T[i][i];}}for(i=0;i<n;i++)for(j=0;j<n;j++){T[i][j]=0;}for(i=n-1;i>=0;i++){r[i][i]=1/T[i][i];for(j=i+1;j<n;j++){s=0;for(k=i;k<j-1;k++){s-=r[i][k]*T[k][j];}r[i][j]=s/T[i][i];}}for(i=0;i<n;i++)for(j=0;j<n;j++){r[i][j]=0;}transformmatrix(r,rt,n,n)mulmatrix(r,rt,Q,n,n)}9.平差主程序之读入数据typedef struct POINT{char name[8];double x,y;int type;}POINT;typedef struct READV ALUE{POINT *begin;POINT *end;double value;}READV ALUE;POINT *GETPOINT(char *name,POINT *pPoint,int nPoint) {int i;for(i=0;i<nPoint;i++){if (strcmp(pPoint[i].name,name)==0)return (pPoint+i)}for(i=0;i<nPoint;i++){if(pPoint[i]=NULL)strcmp(pPoint[i].name,name);pPoint[i].type=0;return(pPoint+i);}}double AngleToRadian（double angle）{int D，M；double S,radian,degree, angle,MS；D=int（angle+0.3）；MS=angle-D；M=int（（MS）*100+0.3）；S=（MS*100-M）*100;degree=D+M/60.0+S/3600.0;radian=degree*PI/180.0;return radian;}main(){POINT *pPoint=NULL;READV ALUE *pDirect=NULL;READV ALUE *pDistance=NULL;int nPoint,nKnownPoint,nDirect,nDistance,i;double mo,mf,ms;char begin[8],end[8];FILE *fp=0;fp=fopen(“c:\\dat\\t1.txt”,”r”)fscanf(fp,”%d,%d,%d,%d\n”,&nPoint,&nKnowPoint,&nDirect,&nDistance) if(nPoint>0)pPoint=(POINT*)malloc(nDirect*sizeof(POINT));if(nDirect>0)pDirect=(READV ALUE*)malloc(nDirect*sizeof(READV ALUE));if(nDistance>0)pDistance=(READV ALUE*)malloc(nDistance*sizeof(RAADV ALUE));fscanf(fp,”%lf,%lf,%lf\n”,&mo,&mf,&ms);for(i=0;i<nKnownPoint;i++){fscanf(fp,”%s,%lf,%lf\n”,pPoint[i].name,&pPoint[i].x,&pPoint[i].y);type=1;}for( ;i<nPoint;i++){pPoint[i].name=NULL;pPoint[i].x=0;pPoint[i].y=0;pPoint[i].type=0;}for(i=0;i<nDirect;i++){fscanf(fp,”%s,%s,%lf\n”,begin,end,&pDirect[i].value);pDirect[i].begin=GetPoint(begin,pPoint,nPoint);pDirect[i].end=GetPoint(end,pPoint,nPoint);}for(i=0;i<nDistance;i++){fscanf(fp,”%s,%s,%lf\n”,begin,end,&pDistance[i].value);pDistance[i].begin=GetPoint(begin,pPoint,nPoint);pDistance[i].end=GetPoint(end,pPoint,nPoint);}fclose(fp);}10.角度检验（checkangle）#include <math.h>int checkangle(double angle){int M,S;double MS;if(angle>=0&&angle<360){MS=angle-(int)(angle);if(M<6){S=(int)(MS*1000);if(S%10<6){return 1;}}}return 0;}11.前方交会#define PI=3014159265/***此处调用程序角度换弧度AngleToRadian***/Qianfang(double XE, double YE, double XF, double YF, doubleDEG, double DEF, double DFG, double DFE, double *DFE, double *DFG){double C,A,B;C=DGE-DGF;A=DEF-DEG;B=DFG-DFE;if((C<-PI&&C>-2*PI)||(C>0&&C<PI){XG=(XE/tan(B)+XF/tan(A)-YE+YF)/(1/tan(A)+ 1/tan(B);YG=(YE/tan(B)+YF/tan(A)+XE-XF)/ (1/tan(A)+ 1/tan(B);}if((C>-PI&&C<0)||(C>PI&&C<2*PI)){XG=(XE/tan(B)+XF/tan(A)+YE-YF)/(1/tan(A)+ 1/tan(B);YG=(YE/tan(B)+YF/tan(A)-XE+XF)/ (1/tan(A)+ 1/tan(B);}}12.坐标概算全方向法子函数取出观测方向GetAllDirectint GetAllDirect(char *name,int nDirect,READV ALUE *pDirect, READV ALUE *pStation) {int i,nCount=0;for(i=0;i<nDirect;i++)if(strcmp(pDirect[i].begin->name,name)==0)){pStation[nCount].begin=p(pDirect[nCount].begin;pStation[nCount].end=p(pDirect[nCount].end;pStation[nCount].value=p(pDirect[nCount].value;nCount++;}return nCount;}坐标概算全方向法子程序实现流程(coordinate)coordinate (入口参数设置){READV ALUE pStation[50],pObject[50];int nCount,i,j,k,m,n,p,nobject;for(i=0;i<nPoint;i++){nCount=GetAllDirect(pPoint[i].name,nDirect,pStation)if((nCount>1)||( nCount=1)){for(j=0;j<nCount;j++){if(pStation[j].end->type==1){for(k=0;k<nCount;k++){if(pStation[k].end->type==0)nobject=GetAllDirect(pStation[j].end->name,nDirect,pDirect,pobject)m=-1;n=-1;for(p=0;p<nobject;p++){if(strcmp(pobject[p].end->name,pPoint[i].name)==0){m=p;}if(strcmp(pobject[p].end->name,pStation[k].end->name)==0){n=p;}if(m>=0&&n>=0){pPoint[i]=pStation[k].end-pStation[j].end;pStation[j].end=pObject[m].value-pObject[n].value;{Xe=pPoint[i].x;Ye=pPoint[i].y;Xf=pStation[j].end->x;Yf=pStation[j].end->y;Lef=pStation[j].value;Leg=pStation[k].value;Lfe=pObject[m].value;Lfg=pObject[n].value;Qianfang(Xe,Xf,Ye,Yf,Lef,Leg,Lfe,Lfg,*Xg,*Yg;)pStation[k].end->x=*xg;pStation[k].end->y=*yg;pStation[k].end.type=2;}}}}}}}}13.坐标增量法(calcoordinate)子函数由端点名称得边长值的函数GetDistancedouble GetDistance(char *begin,char *end,int nDistance,READV ALUE *pDistance){int i;for(i=0;i<nDistance;i++){if((strcmp(pDistance[i].begin->name,begin)==0&&strcmp(pDistance[i].end->name,end==0)| |(strcmp(pDistance[i].begin->name,end)==0&&strcmp(pDistance[i].end,begin)==0))return pDistance[i].value;}return -1;}/***函数取出观测方向GetAllDirect***/void calcoordinate(int nDirect,READV ALUE *pDirect,int nDistace,READV ALUE *pDistance,int nPoint,POINT *pPoint){int nPoint,nCount,nDirect,nDistance;int m=-1,i,j,k;double x1,y1,x2,y2,A0,A,S,dx,dy;READV ALUE*pDirect=NULL;READV ALUE pStation[50];for(i=0;i<nPoint;i++){if(pPoint[i].type>0){nCount=GetAllDirect(pPoint[i].name,nDirect,pDirect,pStation[50]);for(j=0;j<nCount;j++){if(pStation[j].end->type>0)m=j;if(m!=-1){for(k=0;k<nCount;k++){if(pStation[k].end->type==0){x1=pPoint[i].x;y1=pPoint[i].y;x2=pStation[j].end->x;y2=pStation[j].end->y;A0=Bearing(x1,y1,x2,y2);A=A0-(DMSToRAD(pStation[m].value)-DMSToRAD(pStation[k].value));if(A<0)A=A+2*PI;if(A>2*PI)A=A-2*PI;S=GetDistance(pPoint[i],pStation[k].end,nDistance,pDistance);if(S<0)continue;else{dx=S*cos(A);dy=S*sin(A);pStation[k].end->x=pPoint[i].x+dx;pStation[k].end->y=pPoint[i].y+dy;pStation[k].end->type=2;}}}}}}}}14.高斯正反算高斯正算：#include <math.h>#include <stdio.h>#define PI 3.14159265double DMSToRAD（double dDMS）{int L1,L2；double T,L3；L1=(int)（dDMS+0.3）；L2=(int)((dDMS-L1)*100+0.3)；L3=((dDMS-L1)*100-L2)*100；T=(L1+L2/60.0+L3/3600.0)*PI/180.0;return T;}void PreGausePositive(double B,double L,double L0, double a, double b, double *N, double *l, double *c, double *t, double *X,double *B1){double a0,a2,a4,a6,a8,m0,m2,m4,m6,m8;double e,e1;e=(sqrt(a*a-b*b))/a;e1=(sqrt(a*a-b*b))/b;B1=DMSToRAD(B);t=tanB1;c=sqrt(e1*e1*cosB1*cos*B1);l=L-L0;N=a/(sqrt(1-e*e*sinB1*sinB1));m0=a*(1-e*e);m2=3/2*e*e*m0;m4=5/4*e*e*m2;m6=7/6*e*e*m4;m8=9/8*e*e*m6;a0=m0+m2/2+3*m4/8+5*m6/16+35*m8/128;a2=m2/2+m4/2+15*m6/32+7/16*m8;a4=m4/8+3*m6/16+7*m8/32;a6=m6/32+m8/16;a8=m8/128;X=a0*B1-a2*(sin(2*B1))/2+a4*(sin(4*B1))/4-a6*(sin(6*B1))/6+a8*(sin(8*B1))/8;}Void BLToXY(double *x,double *y,double N,double l,double c,double t,double B1,double X)x=X+N*l*l*t*cosB1*cosB1*((3+l*l*cosB1*cosB1*(5-t*t+9*c*c+4*c*c*c*c)/4+l*l*\cosB1*cosB1*(61-58*t*t+t*t*t*t)/30))/6;y=N*l*cosB1(1+l*l*cosB1*((1+c*c-t*t)+l*l*cosB1*cosB1(5-18*t*t+t*t*t*t+14*c*c\-58*t*t*c*c)));}高斯反算void XYToBL(double x,double y,double L0,double a,double b,double q,double *B,\double *L){double Bf,c,t,y,N,e1,ee=(sqrt(a*a-b*b))/(a*a);e1=(sqrt(a*a-b*b))/(b*b);for(Bf=0;;){t=tanBf;c=e1*e1*cosBf;N=a/(sqrt(1-e*e*sinBf*sinBf));B=Bf-(1+c*c)*t*y*y/(2*N*N)*(1-y*y)/(12*N*N)*(15+3*t*t+c*c-9*t*t*c*c)-y*y/ (30*N*N)*(61+90*t*t+45*t*t*t*t);if(fabs(B-Bf)<q)break;Bf=B;}L=L0+y/(N*cosBf)*(1-y*y/(6*N*N))*((1+2*t*t+c*c)-y*y/(20*N*N)*(5+6*c*c)+28*t* t+24*t*t*t*t+8*t*t*c*c)-y*y/(42*N*N)*(61+662*t*t+1320*t*t*t*t+720*t*t*t*t*t*t));B=RADToDMS(B);L=RADToDMS(L);}附录资料：不需要的可以自行删除SHA算法的实现C语言程序：#include <iostream>#include <vector> //定义vector数组#include <string> //记录消息using namespace std;const int NUM = 8; //一个字由32比特(或者8个16进制数)const int BIT = 512; //消息认证码要以512比特一组//字常量string H0 = "67452301";string H1 = "EFCDAB89";string H2 = "98BADCFE";string H3 = "10325476";string H4 = "C3D2E1F0";//定义SHA1(安全哈希算法)类class SHA1{public://将一个字符串形式的字转化为vector数组vector<int> hex_into_dec(string word);//将vector转化为string字符串形式string num_into_message(vector<int> A);//两个字X和Y的逻辑"和"vector<int> word_AND(vector<int> A,vector<int> B);//两个字X和Y的逻辑"或"vector<int> word_OR(vector<int> A,vector<int> B);//两个字X和Y的逻辑"异或"vector<int> word_XOR(vector<int> A,vector<int> B);//两个字X和Y的逻辑"补"vector<int> word_COMPLEMENT(vector<int> A);//两个字X和Y的摸2^32整数加vector<int> word_ADD(vector<int> A,vector<int> B);//将字X循环左移s个位置vector<int> ROTL(vector<int> A,int s);//SHA-1的填充方案,我们设定msg由ASCII码组成vector<vector<int> > SHA_1_PAD(string msg);//将SHA-1压成以字为单位vector<vector<vector<int> > > compress(vector<vector<int> > result);//定义ft函数，每个ft函数都有B,C,D三个字作为输入，并产生一个字作为输出vector<int> Ft(int t,vector<int> B,vector<int> C,vector<int> D);//定义字常数Kvector<int> K(int t);//开始进行SHA-1(安全Hash算法)的加密vector<vector<int> > SHA_1(string msg);};//将vector转化为string字符串形式string SHA1::num_into_message(vector<int> A){int i;string msg = "";for(i = 0;i < A.size();i++){if(A[i] >= 0 && A[i] <= 9)msg += '0' + A[i];else if(A[i] >= 10 && A[i] <= 15)msg += 'A' + (A[i] - 10);}return msg;}//将一个字符串形式的字转化为vector数组vector<int> SHA1::hex_into_dec(string word){int i;vector<int> result(NUM,0);for(i = 0;i < NUM;i++){if(word[i] >= '0' && word[i] <= '9'){result[i] = word[i] - '0';}else if(word[i] >= 'A' && word[i] <= 'F'){result[i] = 10 + word[i] - 'A';}}return result;}//两个字X和Y的逻辑"和"vector<int> SHA1::word_AND(vector<int> A,vector<int> B) {vector<int> result(NUM,0);int i;for(i = 0;i < NUM;i++){result[i] = A[i] & B[i];}return result;}//两个字X和Y的逻辑"或"vector<int> SHA1::word_OR(vector<int> A,vector<int> B) {vector<int> result(NUM,0);int i;for(i = 0;i < NUM;i++){result[i] = A[i] | B[i];}return result;}//两个字X和Y的逻辑"异或"vector<int> SHA1::word_XOR(vector<int> A,vector<int> B) {vector<int> result(NUM,0);int i;for(i = 0;i < NUM;i++){result[i] = A[i] ^ B[i];}return result;}//两个字X和Y的逻辑"补"vector<int> SHA1::word_COMPLEMENT(vector<int> A) {vector<int> result(NUM,0);int i;for(i = 0;i < NUM;i++){result[i] = 15 - A[i];}return result;}//两个字X和Y的摸2^32整数加vector<int> SHA1::word_ADD(vector<int> A,vector<int> B) {vector<int> result(NUM,0);int i;for(i = NUM - 1;i >= 0;i--){result[i] = A[i] + B[i];if(i != 0){int temp = result[i] / 16;result[i-1] += temp;}result[i] %= 16;}return result;}//将字X循环左移s个位置vector<int> SHA1::ROTL(vector<int> A,int s){vector<int> result = A;vector<int> temp(NUM,0);int i,j;for(i = 0;i < s;i++){for(j = NUM - 1;j >= 0;j--){if(result[j] / 8 >= 1){temp[j] = 1;result[j] <<= 1;result[j] %= 16;if(j < NUM - 1)result[j] += temp[j + 1];}else if(result[j] / 8 == 0){temp[j] = 0;result[j] <<= 1;result[j] %= 16;}}result[NUM - 1] += temp[0];}return result;}//SHA-1的填充方案,我们设定msg由ASCII码组成vector<vector<int> > SHA1::SHA_1_PAD(string msg){int len = msg.length();int bit_num = len * 8;int i,j;int num,lest = bit_num % 512;if(lest != 0) //看消息长度是否超过512字节，我们需要将它补成512的倍数num = bit_num / 512 + 1;elsenum = bit_num / 512;//首先我们以8位字节为一组保存到vector里面，512比特为一组，即一组里面有64位元素vector<vector<int> > result;result.resize(num);for(i = 0;i < num;i++){result[i].resize(64);}for(i = 0;i < num;i++){for(j = 0;j < 64 && i * 64 + j < len;j++){result[i][j] = msg[i * 64 + j];}}//下面开始为未够512比特的消息分组进行补长度操作if(lest != 0){int x = num - 1,last_len = lest / 8;result[x][last_len] = 128; //先补一个"1"for(i = last_len + 1;i < 56;i++){result[x][i] = 0;}int last_l = lest;j = 63;while(j >= 56){result[x][j] = last_l % 128;last_l /= 128;j--;}}return result;}//将SHA-1压成以字为单位(三维数组有点复杂)vector<vector<vector<int> > > SHA1::compress(vector<vector<int> > result) {vector<vector<int> > rr;rr.resize(result.size());int i,j;for(i = 0;i < rr.size();i++){rr[i].resize(128);}for(i = 0;i < result.size();i++){for(j = 0;j < result[i].size();j++){rr[i][2 * j] = result[i][j] / 16;rr[i][2 * j + 1] = result[i][j] % 16;}}vector<vector<vector<int> > > rrr;rrr.resize(result.size());for(i = 0;i < rrr.size();i++){rrr[i].resize(16);}for(i = 0;i < rrr.size();i++){for(j = 0;j < 16;j++){rrr[i][j].resize(8);}}for(i = 0;i < rr.size();i++){for(j = 0;j < rr[i].size();j++){rrr[i][j / 8][j % 8] = rr[i][j];}}return rrr;}//定义ft函数，每个ft函数都有B,C,D三个字作为输入，并产生一个字作为输出vector<int> SHA1::Ft(int t,vector<int> B,vector<int> C,vector<int> D){vector<int> result;if(t >= 0 && t <= 19){vector<int> a1 = word_AND(B,C);vector<int> a2 = word_AND(word_COMPLEMENT(B),D);result = word_OR(a1,a2);}else if((t >= 20 && t <= 39) || (t >= 60 && t <= 79)){vector<int> a1 = word_XOR(B,C);result = word_XOR(a1,D);}else if(t >= 40 && t <= 59){vector<int> a1 = word_AND(B,C);vector<int> a2 = word_AND(B,D);vector<int> a3 = word_AND(C,D);vector<int> a4 = word_OR(a1,a2);result = word_OR(a4,a3);}return result;}//定义字常数Kvector<int> SHA1::K(int t){vector<int> result;if(t >= 0 && t <= 19){result = hex_into_dec("5A827999");}else if(t >= 20 && t <= 39){result = hex_into_dec("6ED9EBA1");}else if(t >= 40 && t <= 59){result = hex_into_dec("8F1BBCDC");}else if(t >= 60 && t <= 79){result = hex_into_dec("CA62C1D6");}return result;}//开始进行SHA-1(安全Hash算法)的加密vector<vector<int> > SHA1::SHA_1(string msg){vector<int> h0 = hex_into_dec(H0);vector<int> h1 = hex_into_dec(H1);vector<int> h2 = hex_into_dec(H2);vector<int> h3 = hex_into_dec(H3);vector<int> h4 = hex_into_dec(H4);vector<vector<int> > result1 = SHA_1_PAD(msg);vector<vector<vector<int> > > result2 = compress(result1);int n = result2.size();int i,j;for(i = 0;i < n;i++){vector<vector<int> > W;W.resize(80);for(j = 0;j < 16;j++){W[j] = result2[i][j];}for(j = 16;j < 80;j++){vector<int> a1 = word_XOR(W[j-3],W[j-8]);vector<int> a2 = word_XOR(a1,W[j-14]);vector<int> a3 = word_XOR(a2,W[j-16]);W[j] = ROTL(a3,1);}//将string转化为vector数组vector<int> A = hex_into_dec(H0);vector<int> B = hex_into_dec(H1);vector<int> C = hex_into_dec(H2);vector<int> D = hex_into_dec(H3);vector<int> E = hex_into_dec(H4);for(j = 0;j < 80;j++){vector<int> a1 = ROTL(A,5);vector<int> a2 = Ft(j,B,C,D);vector<int> a3 = word_ADD(a1,a2);vector<int> a4 = word_ADD(a3,E);vector<int> a5 = word_ADD(a4,W[j]);vector<int> temp = word_ADD(a5,K(j));E = D;D = C;C = ROTL(B,30);B = A;A = temp;}h0 = word_ADD(h0,A);h1 = word_ADD(h1,B);h2 = word_ADD(h2,C);h3 = word_ADD(h3,D);h4 = word_ADD(h4,E);}//返回结果(H0||H1||H2||H3||H4)vector<vector<int> > result;result.push_back(h0);result.push_back(h1);result.push_back(h2);result.push_back(h3);result.push_back(h4);return result;}int main(){SHA1 sha1; //定义SHA1算法类string message = "cryptographyisthepracticeandstudyoftechniquesforsecurecommunicationinthepresenceofthirdpartiesmoregenerallyitisaboutconstructingandanalyzingprotocolsthatovercome theinfluenceofadversariesandwhicharerelatedtovariousaspectsininformationsecuritysu chasdataconfidentialitydataintegrityauthenticationandnonrepudiationmoderncryptogra phyintersectsthedisciplinesofmathematicscomputerscienceandelectricalengineeringap plicationsofcryptographyincludeATMcardscomputerpasswordsandelectroniccommerc e";vector<vector<int> > result;result = sha1.SHA_1(message);cout << "消息为:" << endl << message << endl;cout << "利用填充方案SHA-1-PAD给出对消息的填充,得出SHA-1(x)得:" << endl;int i;for(i = 0;i < result.size();i++){cout << sha1.num_into_message(result[i]);}cout << endl;return 0;}程序运行结果：- 21 -。
