长安大学重磁报告三

0 重磁实验报告 三 复 杂 形 体 正 演 姓名 学号 专业 勘查技术与工程 指导教师 鲁宝亮 王万银 完成日期 2013 年 12 月 18 1 日 2 目录 一 基本原理 1 1 重力异常计算公式 1 2 磁力异常计算公式 2 3 化极磁力异常计算公式 3 二 输入 输出数据格式设计 3 1 输入数据文件名的格式设计 3 2 输出数据文件名的格式设计 4 3 重要变量的名称 4 三 总体设计 6 四 测试结果 7 五 结论及建议 9 附录 源程序代码 9 0 一 基本原理 1 O X A x y z Y dv Z 地质体重力异常的计算 重力异常计算公式 v z zyx dddz GV z V g 2 3222 式中 V 为地质体的剩余质量对 A 点的单位质量产生的引力位 v 为地质体 的体积 我们还可以推导出计算重力异常垂向梯度或重力垂向梯度异常的基本公式 为 ddd zyx yxz GV z g v zz 2 5222 222 2 计算重力异常水平梯度或重力水平梯度异常的基本公式为 ddd zyx xz GV x g v xz 2 5222 3 ddd zyx yz GV y g v yz 2 5222 3 一阶导数类边缘识别计算公式 垂向导数 z zyxg zyxVDR 总水平导数 22 y zyxg x zyxg zyxTHDR 解析信号振幅 22 VDRTHDRASM 1 2 O X x y z a b Y c 000 zyx Z 直立长方体示意图 均匀磁化直立长方体磁场各分量表达式 2 2 2 2 2 2 22 0 0 0 0 0 0 0 ln ln arctan 4 c z c z b y b y a x a x zy xax yrMzrM zzrx yx MzyxH 1 2 2 2 2 2 2 22 0 0 0 0 0 0 0 ln arctan ln 4 c z c z b y b y a x a x z yxay xrM zzry yx MzrMzyxH 2 2 2 2 2 2 2 0 0 0 0 0 0 0 arctan ln ln 4 c z c z b y b y a x a x Z yxa zr yx M xrMyrMzyxZ 3 把 1 2 3 式代入磁异常计算公式中得 T 2 2 2 2 2 2 6 22 5 22 4 321 0 0 0 0 0 0 0 arctan arctan arctan ln ln ln 4 c z c z b y b y a x a x zr yx k zzry yx k zzrx yx k zrkyrkxrkMzyxT 其中 为异常体的坐标 21222 zyxr 为地磁场的方向余弦 为总磁化强度 为其在 X Y Z 2 坐标方向的分量 设 M 为的模 I 为磁化强度倾角 D 为磁化强度磁偏角 且 满足关系式 DIMMxcoscos DIMMysincos IMMzsin 3 化极磁力异常磁极因子为 2 2 2 122 111 2 122 222 1 2 vuvui vuvui q q l l 式中 为原磁化方向的方向余弦 为新磁化方向的 1 1 1 222 方向余弦 磁化方向即为磁化强度方向 该公式适合于转换到任意磁化方向 二 输入 输出数据格式设计 1 输入数据文件名格式设计 1 输入数据及文件名对应如下 input source filename source dat input plane terrain filename AXYZL GRD input rugged terrain filename BXYZL GRD body num 3 input source filename 对应 source dat 文件 存放直立六面体相关参 数 input plane terrain filename 对应 AXYZL GRD 文件 存放平面地形数据 input rugged terrain filename 对应 BXYZL GRD 文件 存放曲面地形数据 body num 赋值为 3 表示地质体数量 2 场源参数保存在 source dat 中 第一列为剩余密度 g cm3 第二 列为磁化强度 10 6CGSM 第三列为磁化方向倾角 DEG 第四列为磁化方向与 x 轴的夹角 DEG 第五列为磁化方向与 y 轴的夹角 DEG 第六列 第七列为 x 坐标的起点和终点 km 第八列 第九列为 y 坐标的起点和终点 km 第十 列 第十一列为 z 坐标的起点和终点 km 向下为正 0 7 20000 50 85 5 9 3 4 8 2 7 0 8 34000 50 85 5 0 7 0 5 3 7 0 9 17000 50 85 5 4 8 8 3 2 7 3 起伏地形数据保存在 bxyzl grd 文件中 GRD 格式 平面坐标范 围与曲面坐标范围相同 但其 z 坐标值为 5 3km 形式如下 DSAA 27 27 26 26 26 26 3 5 3 1 7 0 15 0 0 15 0 3 0 4 0 5 0 6 0 7 0 85 可知 网格化数据按 GRD 格式保存 有 27 条线 每条线 27 个点 点 线 坐标从 26m 26m 点 线 距为 2m Zmin Zmax 分别表示 Z 坐标的最小值与 最大值 2 输出数据文件名格式设计 1 输出数据及文件名对应如下 plane gravity anomaly filename plane gravity anomaly grd plane magnetic anomaly filename plane magnetic anomaly grd rugged surface gravity anomaly filename rugged surface gravity anom aly grd rugged surface magnetic anomaly filename rugged surface magnetic an omaly grd plane pole magnetic anomaly filename plane pole magnetic anomaly grd rugged surface pole magnetic anomaly filename rugged surface pole m agnetic anomaly grd 2 对给定三维规则地质体 计算其在给定平面上的重力异常 plane gravi ty anomaly 曲面上的重力异常 rugged surface gravity anomaly 平面上的 磁力异常 plane magnetic anomaly 曲面上的磁力异常 rugged surface magne tic anomaly 平面化极磁力异常 plane pole magnetic anomaly 曲面化极磁 力异常 输出格式如下 DSAA 27 27 26 26 26 26 field min field max 可知 网格化数据按 GRD 格式保存 有 27 条线 每条线 27 个点 点 线 坐标从 26m 26m 点 线 距为 2m field min field max 分别表示网格化 数据的最小值与最大值 3 重要变量的名称 INTEGER m0 m1 m2 m3 m0 m3 为扩边后的起点和终点点位 m1 m2 为原始数据的起点和终点点位 4 INTEGER n0 n1 n2 n3 n0 n3 为扩边后的起点和终点线位 n1 n2 为原始数据的起点和终点线位 INTEGER mpoint nline mpoint 表示点数 nline 表示线数 REAL remaining density 剩余密度 g cm3 REAL magnetic intensity 磁化强度 10 6CGSM REAL magnetic inclination 磁化方向倾角 DEG REAL magnetic x angle 磁化方向与 x 轴的夹角 DEG REAL magnetic y angle 磁化方向与 y 轴的夹角 DEG REAL X source 二维数组 x 坐标的起点和终点 km REAL Y source 二维数组 y 坐标的起点和终点 km REAL Z source 二维数组 z 坐标的起点和终点 km REAL Z coordinate 二维数组 GRD 数据的 Z 坐标 REAL plane terrain data 平面地形数据 REAL rugged terrain data 曲面地形数据 REAL plane gravity anomaly 平面重力异常 REAL plane magnetic anomaly 平面磁力异常 REAL rugged surface gravity anomaly 曲面重力异常 REAL rugged surface magnetic anomaly 曲面磁力异常 REAL plane pole magnetic anomaly 平面化极磁力异常 REAL rugged surface pole magnetic anomaly 曲面化极磁力异常 REAL x original unit vector 原磁化方向方向余弦的 X 坐标 REAL y original unit vector 原磁化方向方向余弦的 Y 坐标 REAL z original unit vector 原磁化方向方向余弦的 Z 坐标 REAL x new unit vector 新磁化方向方向余弦的 X 坐标 REAL x new unit vector 新磁化方向方向余弦的 Y 坐标 REAL x new unit vector 新磁化方向方向余弦的 Z 坐标 REAL expand rugged surface magnetic anomaly 余弦扩边后的曲面磁力异常场 REAL field real FFT 变换后磁力异常场的实部 REAL field image FFT 变换后磁力异常场的虚部 REAL pole factor real 磁极因子实部 REAL pole factor image 磁极因子虚部 dx 点距 5 dy 线距 3 总体设计 A段流程图 A A 定义变量参数 定 义输入输出文件名 读取文件名与 文件参数 读取场源数据 计算平面重力异常 二维 FFT 正变换 计算磁极因子 输出平面 重力异常 输出曲面 重力异常 结束 开始 开始 定义变量参数 定 义输入输出文件名 乘积运算 二维 FFT 反变换 结束 读取地形数据 计算曲面重力异常 计算平面磁力异常 计算曲面磁力异常 输出平面 磁力异常 输出曲面 磁力异常 读入场值 二维余弦扩边 输出平面化 极磁力异常 输出平面化 极磁力异常 6 4 测试结果 25 20 15 10 50510152025 25 20 15 10 5 0 5 10 15 20 25 0 15 30 45 60 75 90 105 120 135 150 165 180 195 210 225 平面重力异常 25 20 15 10 50510152025 25 20 15 10 5 0 5 10 15 20 25 0 15 30 45 60 75 90 105 120 135 150 165 180 195 210 225 240 255 270 曲面重力异常 确定扩边点位 确定扩边点位 7 25 20 15 10 50510152025 25 20 15 10 5 0 5 10 15 20 25 350 300 250 200 150 100 50 0 50 100 150 200 250 300 350 400 450 500 550 600 650 700 750 800 850 900 平面磁力异常 25 20 15 10 50510152025 25 20 15 10 5 0 5 10 15 20 25 800 700 600 500 400 300 200 100 0 100 200 300 400 500 600 700 800 900 1000 1100 1200 曲面磁力异常 25 20 15 10 50510152025 25 20 15 10 5 0 5 10 15 20 25 150 90 30 30 90 150 210 270 330 390 450 510 570 630 690 750 810 8 平面化极磁力异常 25 20 15 10 50510152025 25 20 15 10 5 0 5 10 15 20 25 350 250 150 50 50 150 250 350 450 550 650 750 850 950 1050 1150 曲面化极磁力异常 五 结论及建议 结论 1 在计算重力和磁力异常时 我们使用了三个直立六面体 在用数组存放 六面体的 X Y Z 坐标起点与终点时 我们应该使用三个二维数组 这样在编 程时才更便于计算场点到源点的距离 R 2 平面重力异常与曲面重力异常相比较 等位线显得更加圆滑 曲面重力 异常更能突出局部异常 平面磁力异常与曲面磁力异常相比较 异常的范围更 大 曲面磁力异常更能突出局部异常 3 对磁力异常化极后 我们可以清楚地看到高值异常周围有四块低值区 域 而且分布比较对称 曲面化极磁力异常与平面化极磁力异常相比 有三个 分开的高值区域 更能清楚的突出三个异常体的存在 建议 我们每次做完实验 交完报告后 都不会得到任何来自老师的反馈 那我 们就很难知道自己在写实验报告 编程序过程中的不足 建议老师能指出我们 一些实验方面的不足和改正的办法 附录 源程序代码 复杂形体正演程序设计 PROGRAM COMPLEX BODY FORWARD INTEGER mpoint nline INTEGER body num REAL xmax xmin ymax ymin zmax zmin A B C eigval CHARACTER 120 input source filename input rugged terrain filename 9 CHARACTER 120 plane gravity anomaly filename plane magnetic anomaly filename rugged surface gravity anomaly filename CHARACTER 120 rugged surface magnetic anomaly filename CHARACTER 120 plane pole magnetic anomaly filename rugged surface pole magnetic anomaly filename REAL ALLOCATABLE remaining density REAL ALLOCATABLE magnetic intensity REAL ALLOCATABLE magnetic inclination REAL ALLOCATABLE magnetic x angle REAL ALLOCATABLE magnetic y angle REAL ALLOCATABLE X source REAL ALLOCATABLE Y source REAL ALLOCATABLE Z source REAL ALLOCATABLE Z coordinate REAL ALLOCATABLE plane terrain data REAL ALLOCATABLE rugged terrain data REAL ALLOCATABLE plane gravity anomaly REAL ALLOCATABLE plane magnetic anomaly REAL ALLOCATABLE rugged surface gravity anomaly REAL ALLOCATABLE rugged surface magnetic anomaly REAL ALLOCATABLE plane pole magnetic anomaly REAL ALLOCATABLE rugged surface pole magnetic anomaly REAL ALLOCATABLE field real REAL ALLOCATABLE field image REAL ALLOCATABLE pole factor real REAL ALLOCATABLE pole factor image REAL ALLOCATABLE random terrain call read parameter input source filename input rugged terrain filename plane gravity anomaly filename plane magnetic anomaly filename rugged surface gravity anomaly filename rugged surface magnetic anomal y filename plane pole magnetic anomaly filename rugged surface pole magnetic anomaly filename body num 读取文件名与文件参数 call read source input source filename body num remaining density magnetic intensity magnetic inclination magne tic x angle magnetic y angle X source Y source Z source 读取场源数据 call read parameter grd input rugged terrain filename mpoint nline Xmin Xmax Ymin Ymax Zmin Zmax call read terrain grd input rugged terrain filename mpoint nline Z coordinate 读取曲面地形数据 dx xmax xmin mpoint 1 确定间隔 dy ymax ymin nline 1 ALLOCATE remaining density 1 body num ALLOCATE magnetic intensity 1 body num ALLOCATE magnetic inclination 1 body num ALLOCATE magnetic x angle 1 body num ALLOCATE magnetic y angle 1 body num ALLOCATE X source 1 2 1 body num ALLOCATE Y source 1 2 1 body num ALLOCATE Z source 1 2 1 body num 10 ALLOCATE Z coordinate 1 mpoint 1 nline ALLOCATE plane terrain data 1 mpoint 1 nline ALLOCATE rugged terrain data 1 mpoint 1 nline ALLOCATE plane gravity anomaly 1 mpoint 1 nline ALLOCATE plane magnetic anomaly 1 mpoint 1 nline ALLOCATE rugged surface gravity anomaly 1 mpoint 1 nline ALLOCATE rugged surface magnetic anomaly 1 mpoint 1 nline plane terrain data Zmin rugged terrain data Z coordinate call gravity anomaly calculate remaining density X source Y source Z source Xmin Xmax Ymin Ymax plane terra in data mpoint nline body num plane gravity anomaly 计算平面重力异常 call gravity anomaly calculate remaining density X source Y source Z source Xmin Xmax Ymin Ymax rugged terr ain data mpoint nline body num rugged surface gravity anomaly 计算曲面重力异常 call magnetic anomaly calculate magnetic intensity magnetic inclination magnetic x angle magnetic y angle X so urce Y source Z source Xmin Xmax Ymin Ymax plane terrain data mpoint nline body num plane magnetic anomaly A B C 计算平面磁力异常 call magnetic anomaly calculate magnetic intensity magnetic inclination magnetic x angle magnetic y angle X so urce Y source Z source Xmin Xmax Ymin Ymax rugged terrain data mpoint nline body num rugged surface magnetic anomaly A B C 计算曲面磁力异常 eigval 1 701411e 38 call OUTPUT grd plane gravity anomaly plane gravity anomaly filename 1 mpoint 1 nline mpoint nline eigval xmi n xmax ymin ymax 输出平面重力异常 call output grd rugged surface gravity anomaly rugged surface gravity anomaly filename 1 mpoint 1 nline mpoint nline eigval xmin xmax ymin ymax 输出曲面重力异常 call output grd plane magnetic anomaly plane magnetic anomaly filename 1 mpoint 1 nline mpoint nline eigval xm in xmax ymin ymax 输出平面磁力异常 call output grd rugged surface magnetic anomaly rugged surface magnetic anomaly filename 1 mpoint 1 nline mp oint nline eigval xmin xmax ymin ymax 输出曲面磁力异常 11 call EX EDGE mpoint m0 m1 m2 m3 1 确定扩边点位 call EX EDGE nline n0 n1 n2 n3 1 ALLOCATE field real m0 m3 n0 n3 ALLOCATE field image m0 m3 n0 n3 ALLOCATE pole factor real m0 m3 n0 n3 ALLOCATE pole factor image m0 m3 n0 n3 ALLOCATE random terrain m0 m3 n0 n3 ALLOCATE plane pole magnetic anomaly m0 m3 n0 n3 ALLOCATE rugged surface pole magnetic anomaly m0 m3 n0 n3 call Read field plane magnetic anomaly filename m0 m1 m2 m3 n0 n1 n2 n3 plane magnetic anomaly 读入场值 call Read field rugged surface magnetic anomaly filename m0 m1 m2 m3 n0 n1 n2 n3 rugged surface magnetic an omaly call EXP 2D cos sub m1 m2 m3 m4 n1 n2 n3 n4 plane magnetic anomaly 1 二维余弦扩边 call EXP 2D cos sub m1 m2 m3 m4 n1 n2 n3 n4 rugged surface magnetic anomaly 1 平面化极磁力异常 field real plane magnetic anomaly 确定实部虚部 field image 0 0 call random pole calculate Field Real Field Image m0 m3 n0 n3 dx dy A B C 0 0 0 0 1 0 pole factor real pole facto r image random terrain plane pole magnetic anomaly random terrain call output grd sub plane pole magnetic anomaly filename mpoint nline Xmin Xmax Ymin Ymax m0 m1 m2 m3 n 0 n1 n2 n3 plane pole magnetic anomaly 曲面化极磁力异常 field real rugged surface magnetic anomaly 确定实部虚部 field image 0 0 call random pole calculate Field Real Field Image m0 m3 n0 n3 dx dy A B C 0 0 0 0 1 0 pole factor real pole facto r image random terrain rugged surface pole magnetic anomaly random terrain call output grd sub rugged surface pole magnetic anomaly filename mpoint nline Xmin Xmax Ymin Ymax m0 m1 m 2 m3 n0 n1 n2 n3 rugged surface pole magnetic anomaly DEALLOCATE remaining density DEALLOCATE magnetic intensity DEALLOCATE magnetic inclination DEALLOCATE magnetic x angle DEALLOCATE magnetic y angle 12 DEALLOCATE X source DEALLOCATE Y source DEALLOCATE Z source DEALLOCATE Z coordinate DEALLOCATE plane terrain data DEALLOCATE rugged terrain data DEALLOCATE plane gravity anomaly DEALLOCATE plane magnetic anomaly DEALLOCATE rugged surface gravity anomaly DEALLOCATE rugged surface magnetic anomaly DEALLOCATE plane pole magnetic anomaly DEALLOCATE rugged surface pole magnetic anomaly DEALLOCATE field real DEALLOCATE field image DEALLOCATE pole factor real DEALLOCATE pole factor image DEALLOCATE random terrain END PROGRAM COMPLEX BODY FORWARD 读取文件名与文件参数 SUBROUTINE read parameter input source filename input rugged terrain filename plane gravity anomaly filename plane magnetic anomaly filename rugged surface gravity anomaly filename rugged surface magnetic anomal y filename plane pole magnetic anomaly filename rugged surface pole magnetic anomaly filename body num CHARACTER input source filename input rugged terrain filename CHARACTER plane gravity anomaly filename plane magnetic anomaly filename rugged surface gravity anomaly filename rugged surface magnetic anomaly filename CHARACTER plane pole magnetic anomaly filename rugged surface pole magnetic anomaly filename INTEGER body num OPEN 11 file parameter txt status old READ 11 input source filename READ 11 input rugged terrain filename READ 11 plane gravity anomaly filename READ 11 plane magnetic anomaly filename READ 11 rugged surface gravity anomaly filename READ 11 rugged surface magnetic anomaly filename READ 11 plane pole magnetic anomaly filename READ 11 rugged surface pole magnetic anomaly filename READ 11 body num CLOSE 11 END SUBROUTINE read parameter 读取场源数据子程序开始 SUBROUTINE read source filename body num remaining density magnetic intensity magnetic inclination magnetic x angle m agnetic y angle X source Y source Z source INTEGER body num CHARACTER filename 13 REAL remaining density 1 body num magnetic intensity 1 body num magnetic inclination 1 body num magnetic x angle 1 body num magnetic y angle 1 body num X source 1 2 1 body num Y source 1 2 1 body num Z source 1 2 1 body nu m OPEN 20 file filename status old DO i 1 body num READ 20 remaining density i magnetic intensity i magnetic inclination i magnetic x angle i magnetic y angle i X source 1 i X source 2 i Y source 1 i Y source 2 i Z source 1 i Z source 2 i END DO CLOSE 20 END SUBROUTINE read source 读取场源数据子程序结束 读取地形数据子程序开始 SUBROUTINE read parameter grd filename mpoint nline Xmin Xmax Ymin Ymax Zmin Zmax CHARACTER filename INTEGER mpoint nline REAL Xmin Xmax Ymin Ymax Zmin Zmax OPEN 30 file filename status old READ 30 READ 30 mpoint nline READ 30 Xmin Xmax READ 30 Ymin Ymax READ 30 Zmin Zmax CLOSE 30 END SUBROUTINE read parameter grd SUBROUTINE read terrain grd filename mpoint nline Z coordinate CHARACTER filename INTEGER mpoint nline REAL Z coordinate 1 mpoint 1 nline OPEN 40 file filename status old DO j 1 nline DO i 1 mpoint READ 40 Z coordinate i j END DO END DO CLOSE 40 END SUBROUTINE read terrain grd 读取地形数据子程序结束 计算重力异常子程序开始 SUBROUTINE gravity anomaly calculate remaining density X source Y source Z source Xmin Xmax Ymin Ymax Z coordina te mpoint nline body num gravity anomaly IMPLICIT NONE REAL 8 G PARAMETER G 6 672e 11 INTEGER mpoint nline body num I J K L M N REAL remaining density 1 body num X source 1 2 1 body num Y source 1 2 1 body num Z source 1 2 1 body nu 14 m Z coordinate 1 mpoint 1 nline X coordinate 1 mpoint Y coordinate 1 nline gravity anomaly 1 mpoint 1 nline REAL Xmin Xmax Ymin Ymax sum r dx dy gravity anomaly 0 0 dx xmax xmin mpoint 1 dy ymax ymin nline 1 DO j 1 nline Y coordinate j Ymin j 1 dy 场点 Y 坐标 DO k 1 mpoint X coordinate k Xmin k 1 dx 场点 X 坐标 DO i 1 body num 地质体数量 sum 0 0 DO l 1 2 源点 X 坐标 DO m 1 2 源点 Y 坐标 DO n 1 2 源点 Z 坐标 R sqrt X source l i X coordinate k 2 Y source m i Y coordinate j 2 Z source n i Z coordinate k j 2 sum sum 1 l m n 1 X source l i X coordinate k log r Y source m i Y coordinate j Y source m i Y coordinate j log r X source l i X coordinate k Z source n i Z coordinate k j atan X source l i X coordinate k Y source m i Y coordinate j r Z source n i Z coordinate k j END DO END DO END DO gravity anomaly k j gravity anomaly k j G remaining density i sum END DO END DO END DO END SUBROUTINE gravity anomaly calculate 计算重力异常子程序结束 计算磁力异常子程序开始 SUBROUTINE magnetic anomaly calculate magnetic intensity magnetic inclination magnetic x angle magnetic y angle X so urce Y source Z source Xmin Xmax Ymin Ymax Z coordinate mpoint nline body num magnetic anomaly A B C IMPLICIT NONE INTEGER mpoint nline body num I J K L M N REAL magnetic intensity 1 body num magnetic inclination 1 body num magnetic x angle 1 body num magnetic y angle 1 body num REAL X source 1 2 1 body num Y source 1 2 1 body num Z source 1 2 1 body num REAL Xmin Xmax Ymin Ymax dx dy sum REAL R A B C k1 k2 k3 k4 k5 k6 REAL X coordinate 1 mpoint Y coordinate 1 nline Z coordinate 1 mpoint 1 nline magnetic anomaly 1 mpoint 1 nlin e magnetic anomaly 0 0 dx xmax xmin mpoint 1 dy ymax ymin nline 1 DO j 1 nline Y coordinate j Ymin j 1 dy 场点 Y 坐标 15 DO k 1 mpoint X coordinate i Xmin k 1 dx 场点 X 坐标 DO i 1 body num 地质体数量 sum 0 0 DO l 1 2 源点 X 坐标 DO m 1 2 源点 Y 坐标 DO n 1 2 源点 Z 坐标 R sqrt X source l i X coordinate k 2 Y source m i Y coordinate j 2 Z source n i Z coordinate k j 2 A cos magnetic inclination i cos magnetic x angle i B cos magnetic inclination i sin magnetic x angle i C sin magnetic inclination i k1 2 B C k2 2 A C k3 2 A B k4 A 2 k5 B 2 k6 C 2 sum sum 1 l m n k1 log R X source l i X coordinate k k2 log R Y source m i Y coordinate j k3 log R Z source n i Z coordinate k j k4 atan X source l i X coordinate k Y source m i Y coordinate j X source l i X coordinate k 2 r Z source n i Z coordinate k j Z source n i Z coordinate k j 2 k5 atan X source l i X coordinate k Y source m i Y coordinate j Y source m i Y coordinate j 2 r Z source n i Z coordinate k j Z source n i Z coordinate k j 2 k6 atan X source l i X coordinate k Y source m i Y coordinate j r Z source n i Z coordinate k j END DO END DO END DO magnetic anomaly k j magnetic anomaly k j sum magnetic intensity i 1e 7 END DO END DO END DO END SUBROUTINE magnetic anomaly calculate 计算磁力异常子程序结束 确定扩边点位子程序开始 SUBROUTINE EX EDGE mpoint m0 m1 m2 m3 flag INTEGER mpoint m0 m1 m2 m3 flag INTEGER factor m mu mtemp mtemp mpoint factor m 1 DO WHILE mod mtemp 2 eq 0 and mtemp ne 0 mtemp mtemp 2 End DO IF mtemp eq 1 THEN m3 mpoint 2 ELSE mu int log float mpoint 0 693147 factor m m3 2 mu if flag 1 then