test.cpp

该程序是重力勘探中对剩余异常场的求取算法

#include"stdio.h"
#include"math.h"
#include"stdlib.h"


const float pi=3.141592654f;
void main()
{
	FILE *fp1,*fp2;

   int winXRadius=4;
   int winYRadius=7;
   int i=0;
   int j=0;
   int n=0;
   if((fp1=fopen("pole.grd","r"))==NULL)
   {
	   printf("Cannot open the file\n");
	   exit(0);
   }

   int row,col;
   char s[4];
   float rowStart,rowEnd,colStart,colEnd,min,max;
   fscanf(fp1,"%s\t %d\t %d\t %f\t %f\t %f\t %f\t %f\t %f\t",s,&col,&row,&rowStart,&rowEnd,&colStart,&colEnd,&min,&max);
   float* f=new float[row*col];


   while(!feof(fp1))
   {
      fscanf(fp1,"%f\t",&f[n++]);
   }
   float** data=new float*[row+winYRadius*2];
   for(i=0;i<row+winYRadius*2;i++)
		data[i]= new float[col+winXRadius*2];
   float** dataNew=new float*[row+winYRadius*2];
   for(i=0;i<row+winYRadius*2;i++)
		dataNew[i]= new float[col+winXRadius*2];
   n=0;
   for(i=0;i<row;i++)
   {
	   for(j=0;j<col;j++)
	   {
		   data[i][j]=f[n++];
	   }
   }

   
   /*for(i=0;i<row;i++)
   {
	   
		   for(j=0;j<winXRadius;j++)
		   {
			   dataNew[i][j]=data[i][0]*(1-cos(pi/winXRadius*j))/2;
		   }
		   for(j=winXRadius;j<winXRadius+col;j++)
		   {
			   dataNew[i][j]=data[i][j-winXRadius];
		   }
		   for(j=winXRadius+col;j<winXRadius+col+winXRadius;j++)
		   {
			   dataNew[i][j]=data[i][col-1]*(1-cos(pi/winXRadius*(2*winXRadius+col-j-1)))/2;
		   }
	   
   }

   for(j=0;j<col+2*winXRadius;j++)
   {
	   
		   for(i=0;i<winYRadius;i++)
		   {
			   data[i][j]=dataNew[0][j]*(1-cos(pi/winYRadius*i))/2;
		   }
		   for(i=winYRadius;i<winYRadius+row;i++)
		   {
			   data[i][j]=dataNew[i-winYRadius][j];
		   }
		   for(i=winYRadius+row;i<winYRadius+row+winYRadius;i++)
		   {
			   data[i][j]=dataNew[row-1][j]*(1-cos(pi/winYRadius*(2*winYRadius+row-i-1)))/2;
		   }
	   
   }*/
   for(i=0;i<row;i++)
   {
	       float g=0.0;
	       for(j=0;j<winXRadius;j++)
			   g=g+data[i][j]-data[i][j+1];
		   g=g/winXRadius;
	   
		   for(j=0;j<winXRadius;j++)
		   {
			   dataNew[i][j]=data[i][0]+g*(winXRadius-j);
		   }
		   for(j=winXRadius;j<winXRadius+col;j++)
		   {
			   dataNew[i][j]=data[i][j-winXRadius];
		   }
		   g=0.0;
		   for(j=col-winXRadius;j<col;j++)
				g=g+data[i][j]-data[i][j-1];
		   g=g/winXRadius;
		   for(j=winXRadius+col;j<winXRadius+col+winXRadius;j++)
		   {
			   dataNew[i][j]=data[i][col-1]+g*(j-col-winXRadius+1);
		   }
	   
   }

   for(j=0;j<col+2*winXRadius;j++)
   {
		   float g=0.0;
	       for(i=0;i<winYRadius;i++)
			   g=g+dataNew[i][j]-dataNew[i+1][j];
		   g=g/winYRadius;
	   
		   for(i=0;i<winYRadius;i++)
		   {
			   data[i][j]=dataNew[0][j]+g*(winYRadius-i);
		   }
		   for(i=winYRadius;i<winYRadius+row;i++)
		   {
			   data[i][j]=dataNew[i-winYRadius][j];
		   }
		   g=0.0;
		   for(i=row-winYRadius;i<row;i++)
				g=g+dataNew[i][j]-dataNew[i-1][j];
		   g=g/winYRadius;
		   for(i=winYRadius+row;i<winYRadius+row+winYRadius;i++)
		   {
			   data[i][j]=dataNew[row-1][j]+g*(i-row-winYRadius+1);
		   }
	   
   }

   FILE* fp3;
   fp3=fopen("temp.txt","w");
   for(i=0;i<row+2*winYRadius;i++)
   {
	   fprintf(fp3,"\n");
	   for(j=0;j<col+2*winXRadius;j++)
	   {
		   fprintf(fp3,"%f\t",data[i][j]);
	   }
   }
   fclose(fp3);

   float a,b,c,d,e,B,B1,deltB1x,deltB1y,deltBx,deltBy,q;
   float** tempData=new float*[row+winYRadius*2];
   for(i=0;i<row+winYRadius*2;i++)
		tempData[i]= new float[col+winXRadius*2];

   int times=0;

   while(times<5)
   {
	   for(i=0;i<row;i++)
	   {
		   for(j=0;j<col;j++)
		   {
			   B=data[i+winYRadius][j+winXRadius];
			   
			   B1=(data[i+winYRadius*2][j+winXRadius]+data[i][j+winXRadius]+data[i+winYRadius][j+winXRadius*2]+data[i+winYRadius][j])/4;
			   deltB1x=B-(data[i+winYRadius*2][j+winXRadius]+data[i][j+winXRadius])/2;
			   deltB1y=B-(data[i+winYRadius][j+winXRadius*2]+data[i+winYRadius][j])/2;
			   deltBx=data[i+winYRadius*2][j+winXRadius]-data[i][j+winXRadius];
			   deltBy=data[i+winYRadius][j+winXRadius*2]-data[i+winYRadius][j];
			   if(deltBx==0 && deltB1x==0)
			   {
				   b=1;
			   }
			   else
			   {
				   b=deltBx*deltBx/(deltB1x*deltB1x+deltBx*deltBx);
			   }
			   if(deltBy==0 && deltB1y==0)
			   {
				   c=1;
			   }
			   else
			   {
				   c=deltBy*deltBy/(deltB1y*deltB1y+deltBy*deltBy);
			   }			   
			   a=b+c;
			   d=1-a/2;
			   e=a/2;
			   q=d*B1+e*B;
			   data[i+winYRadius][j+winXRadius]=q;
			   tempData[i][j]=q;
			   
		   }
	   }
	   /*for(i=0;i<row;i++)
	   {
		   
			   for(j=0;j<winXRadius;j++)
			   {
				   dataNew[i][j]=tempData[i][0]*(1-cos(pi/winXRadius*j))/2;
			   }
			   for(j=winXRadius;j<winXRadius+col;j++)
			   {
				   dataNew[i][j]=tempData[i][j-winXRadius];
			   }
			   for(j=winXRadius+col;j<winXRadius+col+winXRadius;j++)
			   {
				   dataNew[i][j]=tempData[i][col-1]*(1-cos(pi/winXRadius*(2*winXRadius+col-j-1)))/2;
			   }
		   
	   }

	   for(j=0;j<col+2*winXRadius;j++)
	   {
		   
			   for(i=0;i<winYRadius;i++)
			   {
				   data[i][j]=dataNew[0][j]*(1-cos(pi/winYRadius*i))/2;
			   }
			   for(i=winYRadius;i<winYRadius+row;i++)
			   {
				   data[i][j]=dataNew[i-winYRadius][j];
			   }
			   for(i=winYRadius+row;i<winYRadius+row+winYRadius;i++)
			   {
				   data[i][j]=dataNew[row-1][j]*(1-cos(pi/winYRadius*(2*winYRadius+row-i-1)))/2;
			   }
		   
	   }*/
	   for(i=0;i<row;i++)
	   {
		   
			   float g=0.0;
			   for(j=0;j<winXRadius;j++)
				   g=g+tempData[i][j]-tempData[i][j+1];
			   g=g/winXRadius;
		   
			   for(j=0;j<winXRadius;j++)
			   {
				   dataNew[i][j]=tempData[i][0]+g*(winXRadius-j);
			   }
			   for(j=winXRadius;j<winXRadius+col;j++)
			   {
				   dataNew[i][j]=tempData[i][j-winXRadius];
			   }
			   g=0.0;
			   for(j=col-winXRadius;j<col;j++)
					g=g+tempData[i][j]-tempData[i][j-1];
			   g=g/winXRadius;
			   for(j=winXRadius+col;j<winXRadius+col+winXRadius;j++)
			   {
				   dataNew[i][j]=tempData[i][col-1]+g*(j-col-winXRadius+1);
			   }
		
		   
	   }

	   for(j=0;j<col+2*winXRadius;j++)
	   {
		       float g=0.0;
			   for(i=0;i<winYRadius;i++)
				   g=g+dataNew[i][j]-dataNew[i+1][j];
			   g=g/winYRadius;
		   
			   for(i=0;i<winYRadius;i++)
			   {
				   data[i][j]=dataNew[0][j]+g*(winYRadius-i);
			   }
		   
			   for(i=winYRadius;i<winYRadius+row;i++)
			   {
				   data[i][j]=dataNew[i-winYRadius][j];
			   }
			   g=0.0;
			   for(i=row-winYRadius;i<row;i++)
					g=g+dataNew[i][j]-dataNew[i-1][j];
			   g=g/winYRadius;
			   for(i=winYRadius+row;i<winYRadius+row+winYRadius;i++)
			   {
				   data[i][j]=dataNew[row-1][j]+g*(i-row-winYRadius+1);
			   }
			   
	   }
	   times++;
   }

   n=0;
   float** residualField=new float*[row];
   for(i=0;i<row;i++)
		residualField[i]= new float[col];
   float** regionField=new float*[row];
   for(i=0;i<row;i++)
		regionField[i]= new float[col];
   for(i=0;i<row;i++)
   {
	   for(j=0;j<col;j++)
	   {
		   float aa=f[n];
		   float bb=data[i+winYRadius][j+winXRadius];
		   residualField[i][j]=aa-bb;
		   regionField[i][j]=bb;
		   n++;
	   }
   }

   min=max=residualField[0][0];
   float min1=regionField[0][0];
   float max1=regionField[0][0];
   for(i=0;i<row;i++)
   {
	   for(j=0;j<col;j++)
	   {
		   if(residualField[i][j]<min)
			   min=residualField[i][j];
		   if(residualField[i][j]>max)
			   max=residualField[i][j];
		   if(regionField[i][j]<min1)
			   min1=regionField[i][j];
		   if(regionField[i][j]>max1)
			   max1=regionField[i][j];
	   }
   }

   if((fp2=fopen("residualField.grd","w"))==NULL)
   {
	  printf("Cannot open the file\n");
	  exit(0);
   }
   n=0;
   fprintf(fp2,"DSAA\n");
   fprintf(fp2,"%d\t %d\n",col,row);
   fprintf(fp2,"%f\t %f\n",rowStart,rowEnd);
   fprintf(fp2,"%f\t %f\n",colStart,colEnd);
   fprintf(fp2,"%f\t %f\n",min,max);
   for(i=0;i<row;i++)
   {
	   for(j=0;j<col;j++)
	   {
		   fprintf(fp2,"%f\n",residualField[i][j]);
	   }
   }
   fclose(fp2);
   if((fp2=fopen("regionField.grd","w"))==NULL)
   {
	  printf("Cannot open the file\n");
	  exit(0);
   }
   n=0;
   fprintf(fp2,"DSAA\n");
   fprintf(fp2,"%d\t %d\n",col,row);
   fprintf(fp2,"%f\t %f\n",rowStart,rowEnd);
   fprintf(fp2,"%f\t %f\n",colStart,colEnd);
   fprintf(fp2,"%f\t %f\n",min1,max1);
   for(i=0;i<row;i++)
   {
	   for(j=0;j<col;j++)
	   {
		   fprintf(fp2,"%f\n",regionField[i][j]);
	   }
   }
   fclose(fp2);
   fclose(fp1);

    for(i=0;i<row+winYRadius*2;i++)
	{
		delete[] data[i];
	    delete[] dataNew[i];
		delete[] tempData[i];
	}
	delete[] data;
	delete[] dataNew;
	delete[] tempData;

	for(i=0;i<row;i++)
		delete[] residualField[i];
	delete[] residualField;
	delete[] f;
}