试验项目二 二度体磁异常正演

2016年1月28日 599 次阅读 0 条评论 0 人点赞

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

 



 
五、正演参数

水平圆

柱体

埋深

计算范围(X,Y)

有效磁矩ms

磁倾角(I)

磁偏角(A)

300

(-500,500)

314

45

45

无限延伸厚板

埋深

计算范围(X,Y)

板宽

(2b)

有效磁化强度

板体倾角

磁倾角(I)

磁偏角

300

(-500,

500)

20

1

45

45

45

 

  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;

}

 

 

菜鸟

文章评论(0)