地球物理数据处理长方体重磁异常

#includestdio.h#includemath.h#definepi3.1416#defineG6.67E-3floatatg(floatz,floatx){floatc;if((fabs(x))(1E-15))c=atan(z/x);elsec=(pi/2);returnc;}voidmain(){FILE*fp,*fp1,*fp2,*fp3;fp=fopen(正长方体重力异常.txt,w);fp1=fopen(正长方体磁异常在X方向分量.txt,w);fp2=fopen(正长方体磁异常在Y方向分量.txt,w);fp3=fopen(正长方体磁异常在Z方向分量.txt,w);floatx0,y0,z0,a,b,c,q,R,dg,DX,DY,DZ,Mx,My,Mz,x[2],y[2],z[2],g[101],X[101],Y[101],Z[101];floatM=2000.0,p=2.67,dx=20.0,dy=20.0,xk,yk;printf(请输入正长方体中心坐标x0,y0,z0\n);scanf(%f,%f,%f,&x0,&y0,&z0);printf(请输入磁化率各分量Mx,My,Mz\n);scanf(%f,%f,%f,&Mx,&My,&Mz);printf(请输入正长方体长a宽b高c\n);scanf(%f,%f,%f,&a,&b,&c);x[0]=x0+a/2;x[1]=x0-a/2;y[0]=y0+b/2;y[1]=y0-b/2;z[0]=z0+c/2;z[1]=z0-c/2;inti,j,k,m,n;i=50;{for(j=0;j101;j++){g[j]=0.0;X[j]=0.0;Y[j]=0.0;Z[j]=0.0;for(k=0;k2;k++){for(m=0;m2;m++){for(n=0;n2;n++){q=pow(-1,k+m+n);xk=j*dx;yk=i*dy;R=sqrt(x[k]*x[k]+y[m]*y[m]+z[n]*z[n]);dg=(-G*p*q*((x[k]-xk)*log(y[m]-yk+R)+(y[m]-yk)*log(x[k]-xk+R)+(z[n]*atg(((x[k]-xk)*(y[m]-yk)),(z[n]*R)))));g[j]=g[j]+dg;DX=(q/(4*pi))*(-Mx*atg((y[m]-yk)*z[n],((x[k]-xk)*R))+My*log(R+z[n])+Mz*log(R+y[m]-yk));X[j]=X[j]+DX;DY=(q/(4*pi))*(-My*atg((x[k]-xk)*z[n],((y[m]-yk)*R))+Mx*log(R+z[n])+Mz*log(R+x[k]-yk));Y[j]=Y[j]+DY;DZ=(q/(4*pi))*(-Mz*atg((x[k]-xk)*(y[m]-yk),(z[n]*R))+Mx*log(R+y[m]-yk)-My*log(R+x[k]-yk));Z[j]=Z[j]+DZ;}}}fprintf(fp,%f\n,g[j]);fprintf(fp1,%f\n,X[j]);fprintf(fp2,%f\n,Y[j]);fprintf(fp3,%f\n,Z[j]);}}fclose(fp);fclose(fp1);fclose(fp2);fclose(fp3);}