Kringing.rar_ Kringing_kriging matlab co_kriging matlab code_kri

#include<stdio.h> #include<stdlib.h> #include<math.h> #include"time.h" #include"rinv.c" #define M 60 #define N 40 #define C0 0 #define C 0.1 #define a 5//sqrt(25/3) main() { int i,j,k,ii,jj; int x1[M],y1[M]; double x2[1681],y2[1681],x[1680],y[1680],z[1680]={0.0},z1[M]; double h[M][M],cov[M+1][M+1],h1[61][1680];//h2[61][1680];//(二维数组和一维的定义要分开，而且二维的行列表示形式要一样，要吗都是数，要吗都含有字母，而且只能是M、N。而且要刚好占满，否则定义新的数组时程序不能正常运行) double Y[M+1],R[M+1]={0.0}; FILE *fp1,*fp2,*fp3; fp1=fopen("tt.txt","w"); fp2=fopen("sj.txt","w"); fp3=fopen("MN.txt","w"); srand((unsigned)time(NULL)); //设随机值 for (i=0;i<M;i++) { x1[i]=rand() % (N); y1[i]=rand() % (N); z1[i]=sin(sqrt(x1[i]*x1[i]+y1[i]*y1[i]))/sqrt(x1[i]*x1[i]+y1[i]*y1[i]); //fprintf(fp1,"%d\t%d\t%9.4f\n",x1[i],y1[i],z1[i]);} } for(i=0;i<=N;i++)//给出所有坐标点 {for(j=40*i+i;j<1681;j++) {x2[j]=i;} } for(i=0;i<=N;i++) for(j=0;j<=N;j++) y2[i*(N+1)+j]=j; for(i=0;i<1680;i++) x[i]=x2[i+1],y[i]=y2[i+1]; /*for (j=0;j<1680;j++) { printf("%f\t%f\n",x[j],y[j]); fprintf(fp2,"%f\t%f\n",x[j],y[j]);} }*/ for (i=0;i<M;i++) for (j=0;j<M;j++) { h[i][j]=sqrt(pow((x1[j]-x1[i]),2)+pow((y1[j]-y1[i]),2)); } /*for (i=0;i<M;i++) {for (j=0;j<M;j++) {fprintf(fp3,"%f ",h[i][j]);} fprintf(fp3,"\n"); }*/ for (i=0;i<M;i++) for (j=0;j<M;j++) cov[i][j]=C-C0+C*(1-exp(-pow(h[i][j],2)/pow(a,2))); for (i=0;i<M;i++) cov[i][M]=1.0; for (j=0;j<M;j++) cov[M][j]=1.0; cov[M][M]=0; /*if(rinv(cov,(M+1))==1) for (i=0;i<=M;i++) {for (j=0;j<=M;j++) {fprintf(fp3,"%9.8f ",cov[i][j]);} fprintf(fp3,"\n",cov[i][j]);}}*/ //if (rinv(cov,M+1)==1) for (k=0;k<M;k++) for (ii=0;ii<1680;ii++) { h1[k][ii]=sqrt((x1[k]-x[ii])*(x1[k]-x[ii])+(y1[k]-y[ii])*(y1[k]-y[ii]));} /*for (k=0;k<M;k++) { for (ii=1;ii<1681;ii++) {fprintf(fp3,"%9.8f ",h1[k][ii]);} fprintf(fp3,"\n",h1[k][ii]); }}*/ for (k=0;k<M;k++) for (ii=0;ii<1680;ii++) {h1[k][ii]=C-C0+C*(1-exp(-pow(h1[k][ii],2)/pow(a,2)));} for (ii=0;ii<1680;ii++) h1[M][ii]=1.0; /*for (k=0;k<=M;k++) { for (ii=1;ii<1681;ii++) {fprintf(fp3,"%9.8f ",h1[k][ii]); } fprintf(fp3,"\n",h1[k][ii]);}}*/ for (k=0;k<M;k++) { Y[k]=z1[k];} Y[M]=0.0; /*for (k=0;k<=M;k++) {printf("%9.4f\n",Y[k]); fprintf(fp3,"%9.4f\n",Y[k]);}}*/ if (rinv(cov,M+1)==1) {{for (j=0;j<=M;j++) for (i=0;i<=M;i++) {R[j]+=Y[i]*cov[i][j];} /*for(j=0;j<=M;j++) {printf("%9.8f\n",R[j]); fprintf(fp3,"%9.8f\n",R[j]);} }*/ for(jj=0;jj<1680;jj++) for (j=0;j<=M;j++) {z[jj]+=R[j]*h1[j][jj];} for(i=0;i<1680;i++) {for(j=0;j<M;j++) if((x[i]==x1[j])&&(y[i]==y1[j])) {z[i]=z1[j];}} } for(j=0;j<120;j++) {printf("%f\t%f\t%9.8lf\n",x[j],y[j],z[j]);} for(j=0;j<1680;j++) {fprintf(fp3,"%f\t%f\t%9.8lf\n",x[j],y[j],z[j]);} } }