#include"stdio.h"
#include"math.h"
#define PI 3.1415926535
#define G 6.67e-11
#define r 200.0
#define M 2000.0
#define thegama 2.67e+3
double atg(double n,double m)
{
double theta;
if (fabs(m)>pow(1.0,-15))
{
if (n/m>0.0)
{
if (m>0.0&&n>0.0) theta=atan(n/m);
else theta=-PI+atan(n/m);
}
if(n/m<0.0) theta=PI+atan(n/m);
}
else
{
if (n<0.0) theta=-PI/2;
else theta=PI/2;
}
return theta ;
}
void main()
{
double e,Mx,Mz,detag,detax,detaz,x1,z1;
double Q,P,H;
double x[1000],z[1000],xk[10000],zk[10000],x0=1000.0,z0=1000.0;
int i,j,l,k;
FILE *f1;
f1=fopen("重磁异常.txt","w");
//******输入拟合边数******//
printf("please put in the number of edges l:\n");
scanf("%d",&l);
//*****输入磁化角度******//
printf("please put in 磁化角度e:\n");
scanf("%lf",&e);
//***角度转换***//
e=e*PI/180;
//*****输入测点数********//
printf("please put in the measurement points k:\n");
scanf("%d",&k);
printf("end of job\n");
for(j=0;j<k;j++)
{
xk[j]=20*j;
zk[j]=0;
}
xk[j]='\0';
zk[j]='\0';
//***磁分量转换***//
Mx=M*cos(e);
Mz=M*sin(e);
fprintf(f1,"重力异常 磁异常X分量 磁异常Z分量\n");
for(i=0;i<k;i++)
{
Q=0.0;P=0.0;H=0.0;
printf("%lf\n",xk[i]);
for(j=0;j<l+1;j++)
{
x[j]=x0+r*sin(j*2*PI/l);
z[j]=z0-r*cos(j*2*PI/l);
x[j]=x[j]-xk[i];
z[j]=z[j]-zk[i];
// printf("%lf %lf\n",x[j],z[j]);
}
x[j]='\0';
z[j]='\0';
for(j=0;j<l;j++)
{
x1=x[j+1]-x[j];
z1=z[j+1]-z[j];
//printf("%lf %lf\n",xx,zz);
//********Q,P,H值计算********//
Q+=z1*(-x1)/(z1*z1+x1*x1)*(atg(z[j],x[j])-atg(z[j+1],x[j+1]))-
0.5*z1*z1/(z1*z1+x1*x1)*log((x[j+1]*x[j+1]+z[j+1]*z[j+1])/(x[j]*x[j]+z[j]*z[j]));
P+=z1*z1/(z1*z1+x1*x1)*(atg(z[j],x[j])-atg(z[j+1],x[j+1]))+
0.5*z1*(-x1)/(z1*z1+x1*x1)*log((x[j+1]*x[j+1]+z[j+1]*z[j+1])/(x[j]*x[j]+z[j]*z[j]));
H+=(0.5*z1*(x[j]*z[j+1]-x[j+1]*z[j])/(z1*z1+x1*x1)*log((x[j+1]*x[j+1]+z[j+1]*z[j+1])/(x[j]*x[j]+z[j]*z[j]))
+x1*(x[j]*z[j+1]-x[j+1]*z[j])/(z1*z1+x1*x1)*(atg(z[j],x[j])-atg(z[j+1],x[j+1])));
}
//printf("%lf\n",H);
//*******重磁异常计算*******//
detag=2*G*thegama*H;
detax=(Mx*P+Mz*Q)/(2*PI);
detaz=(Mx*Q-Mz*P)/(2*PI);
fprintf(f1,"%.8lf %lf %lf\n",detag,detax,detaz);
}
fclose(f1);
}