一、试验要求
根据水平圆柱体和有限延伸厚板正演公式，计算给定参数二度体磁场X分量，Z分量，总场异常，并总结磁异常特征
二、实验课时：4课时
三、实验成果提交 ：各磁场分量及总场主剖面图
四、水平圆柱体及斜交磁化厚板正演公式

 
五、正演参数

1. 水平圆柱体磁异常图
   

（1）水平圆柱体水平X分量图（倾角0-90变化规律）

（2）水平圆柱体垂直Z分量图（倾角0-90变化规律）

（3）总磁异常图（倾角0-90变化规律）

磁异常特征：当垂直磁化时Za为两边有负值的对称曲线，而Ha为原点对称曲线；若为水平磁化，可以知道受磁化倾角的影响比Za大，若水平圆柱体斜磁化，则Ha，Za，均为非对称曲线一般正值小于Za，而负值大于Za。水平磁化的Z分量曲线和垂直磁化时的X分量曲线重合。

1. 厚板磁异常图
   

（1）、厚板水平X分量（角0-90变化规律）

（2）、厚板垂直Z分量（角0-90变化规律）

（3）厚板总磁异常图（角0-90变化规律）

磁异常特征分析：顺层磁化时，Za剖曲线对称于Z轴，且没有负值，Hax曲线对称于原点。相当于Za和Ha分别和不同角函数的合成。当角从零度变到九十度时，Za曲线有轴对称变成原点对称，Hax则由原点对称变为反轴向对称。Za为非对称曲线，其不对称性随角增大而增大。

七、源代码

1、水平圆柱体

#define
_CRT_SECURE_NO_WARNINGS

#include
"stdafx.h"

#include<stdio.h>

#include<math.h>

#include<stdlib.h>

#define
dcx1
"file1.dat"

#define
dcx2
"file2.dat"

#define
dcx3
"file3.dat"

int main()

{

    double Za, Ha, T, M = 1, is = 0.785,I=0.785,r=10, R = 300,A=0.785,is0,x; int i;

    FILE *fp1,*fp2,*fp3;

    is0 = atan(1.0 / cos(A) / tan(I));

    fopen_s(&fp1,dcx1, "a+");

    fopen_s(&fp2, dcx2, "a+");

    fopen_s(&fp3, dcx3, "a+");

    for (x = -500; x <= 500; x = x + 10)

    {

        fprintf(fp1, "%f\t", x);

        fprintf(fp2, "%f\t", x);

        fprintf(fp3, "%f\t", x);

        for (is = 0; is <= 90; is = is + 30)

        {

            Za = 1256 * r * r / 2.0 / pow((x*x + R*R), 2.0)*((R*R - x*x)*sin(is / 180.0 * 3.14) - 2 * R*x*cos(is / 180.0 * 3.14));

            Ha = -1256 * r * r / 2.0 / pow((x*x + R*R), 2.0)*((R*R - x*x)*cos(is / 180.0 * 3.14) + 2 * R*x*sin(is / 180.0 * 3.14));

            T = -1256 * r * r / 2.0 / pow((x*x + R*R), 2.0)*sin(I) / sin(is0)*((R*R - x*x)*cos(is / 180.0 * 3.14 + is0) + 2 * R*x*

                sin(is / 180.0 * 3.14));

            fprintf(fp1, "%f\t", Za);

            fprintf(fp2, "%f\t", Ha);

            fprintf(fp3, "%f\t", T);

            //printf("%f\t%f\t%f\t%f\n", x, Za, Ha, T);

            //fprintf(fp1, "%f\t%f\t%f\t%f\n", x, Za, Ha, T);

        }

        fprintf(fp1, "\n");

        fprintf(fp2, "\n");

        fprintf(fp3, "\n");

    }

    fclose(fp1);

    fclose(fp2);

    fclose(fp3);

    return 0;

}

2、厚板

#define
_CRT_SECURE_NO_WARNINGS

#include
"stdafx.h"

#include<stdio.h>

#include<stdlib.h>

#include<math.h>

#define
dcx1
"file1.dat"

#define
dcx2
"file2.dat"

#define
dcx3
"file3.dat"

int main( )

{

    float u0 = 1256, Ms = 1,a=45,x,r,b=10,h=300,Z,H,T,I=0.785,is;

    FILE *fp1, *fp2, *fp3;

    fopen_s(&fp1, dcx1, "a+");

    fopen_s(&fp2, dcx2, "a+");

    fopen_s(&fp3, dcx3, "a+");//非标准C库函数

        for (x = -500; x <= 500; x = x + 10)

        {

            fprintf(fp1, "%f\t", x);

            fprintf(fp2, "%f\t", x);

            fprintf(fp3, "%f\t", x);

            for (r = 0; r <= 90; r = r + 30)

            {

                is = a - r;

            Z = u0*Ms*sin(a / 180.0*3.14) / 2.0 / 3.14*(0.5*sin(r / 180.0*3.14)*log(((x - b)*(x - b) + h*h) / ((x + b)*(x + b) + h*h)) + cos(r / 180.0*3.14)*(atan((x + b) / h) - atan((x - b) / h)));

            H = u0*Ms*sin(a / 180.0*3.14) / 2.0 / 3.14*(0.5*cos(r / 180.0*3.14)*log(((x - b)*(x - b) + h*h) / ((x + b)*(x + b) + h*h)) - sin(r / 180.0*3.14)*(atan((x + b) / h) - atan((x - b) / h)));

            T = u0*Ms*sin(a / 180.0*3.14) / 2.0 / 3.14*sin(I) / sin(is / 180.0*3.14)*(0.5*cos(a / 180.0*3.14 - 2 * is / 180.0*3.14)*log(((x - b)*(x - b) + h*h) / ((x + b)*(x + b) + h*h)) - sin(1 - 2 * is / 180.0*3.14)*

                (atan((x + b) / h) - atan((x - b) / h)));//计算

            fprintf(fp1, "%f\t", Z);

            fprintf(fp2, "%f\t", H);

            fprintf(fp3, "%f\t", T);

             //printf("%f\t%f\t%f\t%f\n", x, Z, H, T);

            //fprintf(fp, "%f\t%f\t%f\t%f\n", x, Z, H, T);//输出到文件

        }

        fprintf(fp1, "\n");

        fprintf(fp2, "\n");

        fprintf(fp3, "\n");

    }

    fclose(fp1);

    fclose(fp2);

    fclose(fp3);

    return 0;

}